Variational Quantum Eigensolver with Fewer Qubits

TL;DR

This paper introduces a qubit-efficient variational quantum eigensolver that uses fewer qubits to study ground states of quantum many-body systems, enabling high-fidelity results on near-term quantum hardware.

Contribution

It presents a novel tensor network approach with symmetry-preserving circuits that reduces qubit requirements and enhances training efficiency for quantum simulations.

Findings

Achieved over 97% fidelity in ground state estimation with only 6 qubits.

Demonstrated the ability to measure long-range correlations after optimization.

Showed that symmetry-respecting circuits improve training speed and mitigate gradient issues.

Abstract

We propose a qubit efficient scheme to study ground state properties of quantum many-body systems on near-term noisy intermediate scale quantum computers. One can obtain a tensor network representation of the ground state using a number of qubits smaller than the physical degrees of freedom. By increasing the qubits number, one can exponentially increase the bond dimension of the tensor network variational ansatz on a quantum computer. Moreover, we construct circuits blocks which respect U(1) and SU(2) symmetries of the physical system and show that they can significantly speed up the training process and alleviate the gradient vanishing problem. To demonstrate the feasibility of the qubit efficient variational quantum eigensolver in a practical setting, we perform first principle classical simulation of differentiable programming of the circuits. Using only $6$ qubits one can obtain the ground state of a $4\times 4$ square lattice frustrated Heisenberg model with fidelity over 97%. Arbitrarily long ranged correlations can also be measured on the same circuit after variational optimization.

Figures

Tables

Table 1. Table 1: Parameters used in the numerical experiment in Fig.  8 . The number of shots used in sampling is 1024 1024 1024 , the optimizer is Adam with learning rate 0.1 0.1 0.1 . E 𝐸 E and E G subscript 𝐸 𝐺 E_{G} are the lowest ground state energy obtained in training and exact ground state energy respectively.

Lattice	$4\times 4$	$4\times 4$	$6\times 6$	$6\times 6$
Model	Q-MPS	Q-PEPS	Q-MPS	Q-PEPS
Gates	180	180	420	450
Depth $d$	3	5	2	5
$R+V$	6	8	8	12
$E/N$	-0.4596	-0.4681	-0.4147	-0.4707
$E_G/N$	-0.4691	-0.4691	-0.4791	-0.4791

Table 2. Table 3: Variational energy expectation and fidelity of frustrated Heisenberg model with J 2 = 0.5 subscript 𝐽 2 0.5 J_{2}=0.5 for various block depth d 𝑑 d shown in Fig.  3 (c). We only count the number of parametrized gates in the table.

Depth	1	2	3	4	5	Exact
Gates	60	120	180	240	300	–
Energy	-0.454	-0.458	-0.463	-0.464	-0.463	-0.469
Fidelity	0.917	0.923	0.968	0.971	0.968	1.0

Equations

$$H=\sum _{i\alpha }{h}_{\alpha }^i{\sigma }_i^{\alpha }+\sum _{ij\alpha \beta }{h}_{\alpha \beta }^{ij}{\sigma }_i^{\alpha }{\sigma }_j^{\beta }+\dots ,$$

$$⟨H{⟩}_{\mathbf{\theta }}=\sum _{i\alpha }{h}_{\alpha }^i⟨\psi (\mathbf{\theta })|{\sigma }_i^{\alpha }|\psi (\mathbf{\theta })⟩+\sum _{ij\alpha \beta }{h}_{\alpha \beta }^{ij}⟨\psi (\mathbf{\theta })|{\sigma }_i^{\alpha }{\sigma }_j^{\beta }|\psi (\mathbf{\theta })⟩+\dots$$

$$\frac{\partial }{\partial {\theta }_i}⟨H{⟩}_{\mathbf{\theta }}=\frac{1}{2}\left(⟨H{⟩}_{\mathbf{\theta }+\frac{\pi }{2}{\mathbf{e}}_i}-⟨H{⟩}_{\mathbf{\theta }-\frac{\pi }{2}{\mathbf{e}}_i}\right).$$

$$\begin{array}{rl}H=& \frac{1}{4}[\sum _{⟨i,j⟩}{\sigma }_i^x{\sigma }_j^x+{\sigma }_i^y{\sigma }_j^y+{\sigma }_i^z{\sigma }_j^z\\ & +J_2\sum _{\text{set@color}⟨\text{leavevmode}\text{set@color}⟨i,j\text{set@color}⟩\text{leavevmode}\text{set@color}⟩}{\sigma }_i^x{\sigma }_j^x+{\sigma }_i^y{\sigma }_j^y+{\sigma }_i^z{\sigma }_j^z],\end{array}$$

$${\sigma }_g^2=\underset{\begin{array}{c}{\theta }_i\end{array}}{\mathrm{Var}}\left(\frac{\partial ⟨H⟩}{\partial {\theta }_i}\right),$$

$$\mbox{$\left|\psi_{cl}\right\rangle$}=\prod_{i=1}^{N-1}CZ_{i,i+1}\mbox{$\left|+\right\rangle$}^{\otimes N}.$$

$$\mbox{$\left|\psi_{cl}\right\rangle$}=\sum_{i_{1}i_{2}\cdots i_{N}}\mathrm{tr}\mbox{$\left[A^{(1)}_{i_{1}}A^{(2)}_{i_{2}}\cdots A^{(N)}_{i_{N}}\right]$}\mbox{$\left|i_{1}i_{2}\cdots i_{N}\right\rangle$},$$

Full text

I Introduction

Studying ground state properties of quantum many-body systems is a promising native application of quantum computers. Given limited qubit resources and noisy realizations of near-term quantum devices Preskill ; Boixo et al. (2018), a practical approach is to employ the variational quantum eigensolver (VQE) Peruzzo et al. (2014); McClean et al. (2016); Wecker et al. (2015a, b); McArdle et al. ; Cao et al. , which runs in a classical-quantum hybrid mode. In this scheme, a parametrized quantum circuit provides a variational ansatz for the ground state. A classical optimizer tunes the circuit parameters to reduce the expected energy of the target Hamiltonian of the output quantum state. There were already several small scale experimental demonstrations of VQE for molecules and quantum magnets Shen et al. (2017); O’Malley et al. (2016); Kandala et al. (2017); Colless et al. (2018); Hempel et al. (2018). These early experiments mostly employed gradient free or Bayesian approaches for classical optimization. Recent progress on unbiased gradient estimation on quantum circuits Li et al. (2017); Mitarai et al. (2018); Liu and Wang (2018); Verdon et al. ; Schuld et al. ; Vidal and Theis ; Bergholm et al. ; Guerreschi and Smelyanskiy ; Farhi and Neven ; Romero et al. (2018); Harrow and Napp ; Dallaire-Demers et al. breaks the information bottleneck between classical and quantum processors, thus providing a route towards scalable optimization of circuits with a large number of parameters.

There are nevertheless more challenges in the training of variational quantum circuits. The gradients of an unstructured, randomly parametrized circuit vanish exponentially as a function of the number of parameters McClean et al. (2018) due to the concentration of measure in high dimensional spaces  Gross et al. (2009); Bremner et al. (2009). Intuitively, this could be understood by the fact that the overlap between a random initial quantum state and a target state is exponentially small in the many-body Hilbert space. This difficulty motivates one to design the circuit architecture and initialize the circuit parameters with insights from classical tensor networks Huggins et al. (2019); Kim and Swingle ; Tianyi Peng and Wu and quantum chemistry ansatz Yung et al. (2014); Lee et al. (2018); Barkoutsos et al. (2018). Furthermore, since the number of required qubits is the same as the problem size in the standard VQE applications, one has to push up the number of controllable qubits way beyond the current technology to convincingly surpass the classical simulation approach in finding the ground states of quantum many-body systems. Related approaches such as the quantum approximate optimization algorithm Farhi et al. and related field such as quantum machine learning Ciliberto et al. (2018); Mitarai et al. (2018); Liu and Wang (2018) suffer from the same problem.

We address these problems by adopting the qubit efficient circuit architecture Monras et al. ; Huggins et al. (2019) for the variational quantum eigensolver. By measuring qubits sequentially and reusing the measured qubit, one can produce quantum states for an arbitrarily large system with a fixed number of qubits. This approach amounts to generating matrix product states (MPS) Östlund and Rommer (1995); Vidal (2003); Schollwöck (2011) on a quantum computer Schön et al. (2005); Cramer et al. (2010). Note that, although it is well known that an MPS with small bond dimension can be efficiently simulated classically Vidal (2003), having quantum resource allows one to reach exponentially large bond dimension that is inaccessible to classical computers. Despite its one-dimensional geometry, MPS is a versatile variational ansatz that has been successfully applied to systems with diverse lattice geometry and topology Stoudenmire and White (2012). The success lies in the fact that many physics and chemistry systems of interest exhibit relatively small entanglement entropy in their ground states Verstraete et al. (2004); Eisert et al. (2010).

A natural extension of MPS to higher dimensions is projected entangled pair states (PEPS) Verstraete et al. (2006), which is able to represent many area law entangled two dimensional quantum states with only polynomial number of parameters. However, classical algorithm for exact contraction of a finite size PEPS shows exponentially complexity with the problem size. We propose a PEPS-inspired quantum circuit ansatz for the ground states of two-dimensional quantum systems. The ansatz can capture area law entanglement entropy. Although classical simulation of this ansatz is exponential difficult, one can estimate the unbiased energy expectation of the ansatz on a quantum computer efficiently.

Variational optimization of MPS generated on an actual quantum device has previously been demonstrated in an experiment Eichler et al. (2015). The experiment exploits the fact that the radiation field of a cavity QED naturally realizes Barrett et al. (2013) the continuous MPS Verstraete and Cirac (2010) with a few tuning parameters. Here, we focus on variational MPS and PEPS calculation on programmable gate model quantum computers. This setup provides us with more systematic and precise control of the bond dimension and the number of variational parameters. Moreover, crucial technical advances such as gradient-based learning Li et al. (2017); Mitarai et al. (2018); Liu and Wang (2018); Verdon et al. ; Schuld et al. ; Vidal and Theis ; Bergholm et al. ; Guerreschi and Smelyanskiy ; Farhi and Neven ; Romero et al. (2018); Harrow and Napp ; Dallaire-Demers et al. and quantum number preserving circuit design greatly speed up the training process, and make it practically useful for solving challenging quantum many-body problems.

This paper is organized as follows. In Sec. II we introduce qubit efficient scheme of preparing MPS/PEPS on quantum circuits and a gradient-based variational training approach for obtaining the ground state of generic quantum many-body Hamiltonians. In Sec. III, we demonstrate the utility of this scheme by numerically simulating the VQE of a frustrated Heisenberg model using fewer qubits than the system size. In Sec. IV, we carry out a detailed gate counting to estimate the timings of actual experiments, followed by a discussion in Sec. V which point to future research directions. Codes and pretrained circuits parameters can be found at the Github repository git .

II Circuit architecture and training approach

Considering a generic quantum many-body Hamiltonian written in terms of the Pauli operators

[TABLE]

where $i,j=1,2\ldots N$ are site indices with $N$ being the system size, $\alpha,\beta=x,y,z$ are indices of the Pauli axis. We only show the first few terms for the sake of concreteness although higher order polynomials of Pauli operators are allowed. In the variational quantum eigensolver approach Peruzzo et al. (2014), the ground state $|\psi({\boldsymbol{\theta}})\rangle$ is represented by the output of a parametrized quantum circuit, where ${\boldsymbol{\theta}}=\left\{\theta_{i}\right\}$ are the circuit parameters. The variational energy is a summation of expectation values of Hamiltonian terms in the variational ground state,

[TABLE]

To estimate the expected energy, one can identify maximally commuting sets of Hamiltonian operators and measure all the commuting terms together on the corresponding bases.

One can prepare an MPS/PEPS as a variational state of $N$ qubits using a smaller number $R+V\ll N$ of qubits Schön et al. (2005); Cramer et al. (2010); Huggins et al. (2019), which will be dubbed Q-MPS/Q-PEPS hence force. The idea is to treat $R$ (for Q-MPS, $R=1$) of them as the physical qubit and use the remaining $V$ qubits as virtual degrees of freedom to mediate quantum entanglement. Sequentially measuring and reusing the physical qubits allow for producing an arbitrarily long MPS/PEPS. The circuit structure is illustrated schematically in the upper left corner of Fig. 1. First, one initializes the $R+V$ qubits to the product state \mbox{\left|\mathbf{0}\right\rangle}_{R}\otimes\mbox{\left|\mathbf{0}\right\rangle}_{V} and applies a circuit block parametrized by ${\boldsymbol{\theta}}_{1}$ to all qubits. Then, one measures the physical qubits on Pauli basis ${\sigma^{\alpha_{1}}_{1}}$ and stores the output ${q^{\alpha_{1}}_{1}}$ to a classical memory. Next, one recycles the measured qubit and reset it to state \mbox{\left|\mathbf{0}\right\rangle}_{R}. One then entangles them with the remaining $V$ qubits again by applying a second circuit block with parameters ${\boldsymbol{\theta}}_{2}$. After repeating these procedures until one has collected $N-(R+V)$ bits of classical information, one measures all qubits to collect the last $R+V$ bits. This sequential measure-and-reuse scheme is equivalent to sampling this $N$-qubit MPS/PEPS on the same basis. Comparing to the qubit saving scheme for generative modeling in LABEL:Huggins2019, here one to perform measurement on multiple different bases for the VQE calculation.

The proposed scheme directly applies to Hamiltonians with arbitrarily long-range interaction. In particular, fermionic systems can be easily studied using the Jordan-Wigner transformation Ortiz et al. (2001). A general quantum chemistry problem is more challenging than the quantum spin problem considered here since they contain $\mathcal{O}(N^{4})$ terms in the Hamiltonian. Nevertheless, the total number of measurements can be reduced using the techniques of Cao et al. ; Babbush et al. (2018); Motta et al. . In Appendix A, we provide a concrete example of preparing and sampling a cluster state in the qubit efficient scheme.

Given the variational circuit, we aim at solving the optimization problem ${\boldsymbol{\theta}}_{\rm opt}={\rm argmin}_{{\boldsymbol{\theta}}}\langle{H}\rangle_{\boldsymbol{\theta}}$. Gradient-based optimization algorithms are crucial to scaling to a large number of variational parameters Liu and Wang (2018). Suppose that all the parameters of the quantum circuit appear in the form $e^{-i\theta_{i}\Sigma/2}$ with $\Sigma^{2}=1$. The analytical expression of the gradient with respect to the parameter $\theta_{i}$ reads Mitarai et al. (2018),

[TABLE]

One can thus estimate the energy gradient by tuning the parameters to $\theta_{i}\pm\pi/2$ and use it for gradient descend optimization of the energy. Unlike numerical differentiation, Eq. (3) is an exact gradient estimator, which is crucial for unbiased stochastic optimization with a noisy estimate of the gradients Tucker et al. (2017).

To recapitulate, the key point of the proposed variational algorithm is to estimate the energy gradient of an $N$-qubit Hamiltonian with respect to an Q-MPS/Q-PEPS with fewer $(R+V\ll N)$ qubits. The steps are shown in Fig. 1. \raisebox{-.9pt} {1}⃝ Tune a selected circuit parameter $\theta_{i}$ to $\theta_{i}+\frac{\pi}{2}$, and collect bit strings by repeated measurements on various bases according to the Hamiltonian terms. Then we repeat for $\theta_{i}-\frac{\pi}{2}$. \raisebox{-.9pt} {2}⃝ Estimate the energy expectation value by assembling the statistics of all Hamiltonian terms \raisebox{-.9pt} {3}⃝ Estimate gradient of all parameters via Eq. (3). Feed the gradient information into a classical optimizer. \raisebox{-.9pt} {4}⃝ Update the circuit parameters according to suggestions of the classical optimizer. This completes one training epoch. The training stops when a prescribed convergence criterion is met. After reaching convergence, one may measure physical observables of interests on the optimized circuits.

III Application to Heisenberg model

As a concrete example, we apply the approach detailed in the previous section to the frustrated Heisenberg model on a square lattice:

[TABLE]

where $\langle i,j\rangle$ and \mathopen{\hbox{\set@color{\langle}}\kern-1.94444pt\leavevmode\hbox{\set@color{\langle}}}i,j\mathclose{\hbox{\set@color{\rangle}}\kern-1.94444pt\leavevmode\hbox{\set@color{\rangle}}} denote nearest and next-nearest neighbors pairs, respectively. $J_{2}>0$ is the strength of the frustration term that suppresses the Neel order. The energy expectation value and its gradient can be efficiently evaluated by sampling the circuit output on three bases $\sigma^{x}$, $\sigma^{y}$ and $\sigma^{z}$. In the following discussion, we consider the model on an open square lattice of the size $N=4\times 4$ with $J_{2}=0.5$. These sites are zigzag ordered in our ansatz as shown in the lower right of Fig. 1.

Frustrated quantum spin models are crucial to the study of quantum magnets with many open problems Zhou et al. (2017); Savary and Balents (2017). Classical computational approaches to these problems are either limited by the sign problem Troyer and Wiese (2005) or high computational cost at larger bond dimensions Orus . Variational optimization of Q-MPS/Q-PEPS on near-term quantum computers is a promising approach which may deliver valuable insights into open problems in this field.

III.1 MPS inspired ansatz with conserved quantum numbers

Figure 3(a) shows a general internal structure of the variational circuit which is efficient to be implemented on quantum hardware Kandala et al. (2017). Each layer contains $3(V+1)$ parameters in the rotational gates $R^{x}_{\theta}$ and $R^{z}_{\theta}$. We use ${\rm CNOT}$ gates with no variational parameters as the entanglers to generate entanglement between qubits. We repeat this construction for $d$ times within each circuit block. Thus there are $M=3d(V+1)$ parameters in each block. As we show below, taking into account the physical symmetries in designing of the VQE ansatz can reduce the number of parameters and increase the training performance.

The Heisenberg model (4) has a ${\rm U}(1)$ symmetry with good quantum number $S_{z}$. To preserve this symmetry Jiang et al. (2018); Barkoutsos et al. (2018), we construct a circuit block consists of $e^{-i\theta\sigma^{z}_{i}/2}$ and $e^{-i\theta{\rm SWAP}(i,j)/2}$ gates 111We note that these two gates are not sufficient to represent all ${\rm U}(1)$ symmetric operations, e.g. the CZ gate. A general ${\rm U}(1)$ ansatz which covers the whole symmetry sector will need all generators such as $\prod\limits_{i}\sigma^{z}_{i}$ and $\prod\limits_{ij}{\rm SWAP}(i,j)$, where $i$ runs over a subset of qubit indices.. The latter gate is equivalent up to a phase factor to the ${\rm SWAP}^{\alpha}$ gate with $\alpha=\theta/\pi$ Williams (2011). Viewing the setup as a wide circuit with $N$ qubits, it is clear that the quantum number of the initial state is conserved during the evolution. To obtain the ground state in $S_{z}=0$ sector, we prepare a spin-balanced initial state for the variational calculation. Fig. 3 (b) shows that, by applying an additional X gate before the variational gates in the odd steps, one has an anti-ferromagnetic product state $|1010\ldots 10\rangle$ as the initial state.

One can further exploit the full ${\rm SU}(2)$ symmetry of the Heisenberg model (4). While there are sophisticated approaches Weichselbaum (2012); Schmoll et al. to implement this non-Abelian symmetry in classical simulations, the implementation is straightforward on quantum circuits. As shown in Fig. 3 (c), we first prepare the input state in the total spin $\mathbf{S}^{2}=0$ sector, where the simplest choice is the singlet product state $|\psi_{0}\rangle=\bigotimes\limits_{i=1}^{N/2}|\uparrow\downarrow-\downarrow\uparrow\rangle$. To prepare $|\psi_{0}\rangle$, we use an additional ancilla qubit to carry the entanglement of the physical qubit in the odd and even steps 222The ancilla qubit increases the bond dimension of the MPS to $2^{V+1}$.. In the odd step, we prepare a spin singlet between the physical qubit and the ancilla qubit $|{\uparrow\downarrow-\downarrow\uparrow}\rangle=\overline{{\rm CNOT}(1,{\rm a})}{\rm H}(1){\rm X}(1)|0\rangle_{1}\otimes|0\rangle_{\rm a}$. Here, $\overline{{\rm CNOT}(1,{\rm a})}$ is the inverse controlled-NOT gate, which flips the ancilla qubit when the physical qubit is in state $|0\rangle$. In the even step, we swap the ancilla and physical qubits. The physical qubits in the odd and even steps thus form a spin singlet. Then, we repeatedly apply parametrized ${\rm SU}(2)$ symmetric operations to the initial state to generate the variational output. We choose the generators to be the ${\rm SWAP}^{\alpha}$ gate Williams (2011) between a collection of qubits pairs 333We note that the variational ansatz may have a connection to the Bethe Ansatz Murg et al. (2012) construction of the Q-MPS, where the scattering matrix has the same form as the parametrized SWAP gates.. Figure 2 shows the expanded view of the ${\rm SU}(2)$ symmetric quantum circuit. Parallel measurement of the final quantum state yields the same outcomes as the sequential measurement in the qubit efficient scheme shown in Fig. 3(c).

The ${\rm SU}(2)$ symmetric variational state reads $|\psi({\boldsymbol{\theta}})\rangle=\prod\limits_{\{i,j\}}e^{-i\theta_{i,j}{\rm SWAP}(i,j)/2}|\psi_{0}\rangle$, where the product is ordered by the circuit architecture. This variational ansatz resembles the classical variational ansatz for quantum spins in the valence bond basis Liang et al. (1988). However, in general, the state could not be sampled efficiently using the classical Monte Carlo method due to the appearance of complex weights Lou and Sandvik (2007). Moreover, since the swap operations are not commuting within each other, there is an additional difficulty in devising an efficient classical Monte Carlo scheme to sample from the variational ansatz. Therefore, variational optimization of this ansatz on a quantum device highlights the possible quantum advantage of the proposed qubit efficient VQE scheme.

To assess the feasibility of the qubit efficient VQE scheme on near-term quantum devices, we perform a faithful classical simulation of the training process. We simulate a circuit of $V+1$ qubits instead of an equivalent $N$ qubits circuit without qubit reusing. Therefore, even in the classical simulation, we do not have direct access to the final wavefunction but only to the measured bit strings. We sample the energy and its gradient on samples of batch size $4096$. Note that we purposely do not exploit the classical backpropagation algorithm (which reduces the complexity of gradient estimation from $\mathcal{O}(M^{2})$ to $\mathcal{O}(M)$, $M$ is the number of variational parameters) to be in line with the repetitive experimental measurement Liu and Wang (2018).

We use $V=4$ qubits for the virtual degrees of the MPS. The maximum entanglement entropy of the ansatz is thus $4\ln 2$ given the full capacity of the variational blocks. We employ the Adam optimizer with a learning rate $0.1$ Kingma and Ba for the stochastic gradient descent training. We compared three different circuit blocks shown in Fig. 3. All of them have fixed depth $d=5$. The variational parameters are random initialized with uniform distribution in $[0,\pi]$. As shown in Fig. 1, there are in total $N-V=12$ circuit blocks.

The general circuit structure shown in Fig. 3(a) contains $M=3(V+1)(N-V)d=900$ variational parameters. In $500$ steps of training, the energy per site decreases to $-0.416$. For comparison, the exact ground state energy per site is $E_{\mathrm{exact}}=-0.46909731$, the density matrix renormalization group result for bond dimension $2^{V}=16$ is $E_{\mathrm{DMRG}}=-0.46532670$ 444For bond dimension $32$, $E_{\mathrm{DMRG}}=-0.46868064$, which represents the entanglement upper bond for the ${\rm SU}(2)$ symmetry preserving ansatz. As the energy decreases, its fidelity with respect to the exact ground state increases from $5.7\times 10^{-3}$ to $0.69$.

Next, the ${\rm U}(1)$ symmetric circuit structure in Fig. 3(b) contains $10$ single qubit gates and $5$ two-qubit gates in each layer. Hence the number of circuit parameters is also $900$. However, the training efficiency increases significantly as shown in Figure 4. The ground state energy per site reaches $-0.454$, with a ground state fidelity of $0.92$.

Finally, the ${\rm SU}(2)$ symmetric circuit structure in Fig. 3(c) gives the best variational energy despite that it only has $M=(V+1)(N-V)d=300$ variational parameters. Using the same hyperparameters for training, the variational energy decreases to $-0.463$, and the fidelity reaches $0.97$. After obtaining the variational state, we can measure physical observables on the circuit. For example, the spin-spin correlation in $z$ direction $\langle{\sigma^{z}_{i}}{\sigma^{z}_{j}}\rangle$ measured on the ${\rm SU}(2)$ symmetry preserving circuit shown in Fig. 5 (a). As a comparison, using the same training hyperparameters we obtain a ground state with fidelity $0.98$ using ${\rm SU}(2)$ symmetric variational circuit ansatz Fig. 3(c) for unfrustrated Heisenberg lattice with $J_{2}=0$. As shown in Fig. 5 (b), the checkerboard pattern for the antiferromagnetic correlation is more visible in the unfrustrated case. These results suggest that, even with a moderate number of qubits, the qubit efficient VQE scheme is able to offer useful physical insights.

This result benefits from the fact that Q-MPS structure alleviates the gradient vanishing problem for studying low entangle physical systems. The gradients of an unstructured, randomly parametrized circuit vanish exponentially as a function of the number of parameters McClean et al. (2018) due to the concentration of measure in high dimensional spaces  Gross et al. (2009); Bremner et al. (2009). Intuitively, this could be understood by the fact that the overlap between a random initial quantum state and a target state is exponentially small in the many-body Hilbert space. Related approaches such as the quantum approximate optimization algorithm Farhi et al. and related field such as quantum machine learning Ciliberto et al. (2018); Mitarai et al. (2018); Liu and Wang (2018) suffer from the same problem.

We inspect the variance of the gradient signal for various system sizes to investigate the gradient vanishing problem in the training variational quantum circuits McClean et al. (2018). To compute the gradient variance, we sample $1000$ gradients of random circuit parameters.

First, we consider an unfrustrated Heisenberg model on an open chain of length $N$. Here $N$ should be regarded as the effective circuit width since the output bit strings lie in the Hilbert space of size $2^{N}$. Fig. 6(a) shows the variance of the gradient for circuit blocks with ${\rm U}(1)$ and ${\rm SU}(2)$ symmetries with $V=4$. Interestingly, the variance of the gradient shows a power-law decay in contrast to the exponential decay found in a circuit with generic structure McClean et al. (2018). Therefore, it appears that the Q-MPS structure alleviates the gradient vanishing problem at least for the problem under consideration. We attribute this to the fact that the low entropy variational ansatz captures the right inductive bias for the ground state of the target problem.

Next, for an $N=20$ Heisenberg chain, we examine the scaling of the gradient variance with the number of virtual qubits $V$. Again, we see that using symmetry greatly enhance the gradient in Fig. 6(b). We observe an exponential decrease of the gradient in the regime $V\ll N$. While the gradient increases with $V$ for $V\lesssim N$. Their values are still much smaller than the values in the small $V$ limit, which shows that the Q-MPS are easier to train compared to an unstructured quantum circuit of a generic structure.

Exponential decreasing gradient with respect to $V$ also warns us that there is no free lunch. Although a quantum circuit is able to represent an MPS with exponentially large bond dimension, the number of parameters in the ansatz should also scale exponentially in order to compensate the vanishing gradient. What’s worse, the sampling error can smear out too small gradient signals. To really make the ansatz scalable for a highly entangled two dimensional system, we need a parameter efficient design.

III.2 PEPS inspired ansatz

PEPS is much more parameter efficient for representing the ground states of a two-dimensional quantum lattice Hamiltonians. It is able to represent many area law entangled high dimensional states with a polynomial number of parameters with respect to the system size. However, there is no polynomial time algorithm to contract a PEPS exactly for energy expectation values. Moreover, there is no efficient scheme to prepare a Q-PEPS on a quantum circuit. This task is probably impossible due to the computational complexity argument Haferkamp et al. , otherwise, a quantum computer would be able to solve the $\#P$-hard problems Biamonte et al. (2015); Johnson et al. (2013); Biamonte and Bergholm . Nevertheless, it is possible to design a variational ansatz that shares the appealing properties of PEPS, such as the 2D area law entanglement entropy.

The SU(2) symmetric Q-PEPS ansatz for a $4\times 4$ square lattice is shown in Fig. 7.

The entanglement in this ansatz satisfies two-dimensional area law, which can be seen from the red dashed box enclosing a $2\times 2$ region, where the number of bonds crossing the boundary is proportional to the circumference of the box. In the entangle layer enclosed in the black dashed box, we first entangle each physical qubit with its own virtual qubit(s), then neighboring (periodic boundary condition) physical qubits, and finally neighboring virtual qubits. This completes a single entangle layer. We repeat this layer for $d$ times to increase the number of trainable parameters. Unlike the case in Q-MPS, where we do not distinguish between virtual qubits, here we assign a constant number of virtual qubits for each physical qubit so that $V/R$ stay as a constant as the system size grows. These physical qubits can be measured in parallel instead of one by one.

The simulation results reported in Fig. 8 demonstrates the scalability of this structure on a 2D lattice. In the simulation, we fix the ratio $V/R$ to $1$ and depth $d$ to $5$. So the number of parameters contained in a $4\times 4$ Q-PEPS is $180$, which is the same as that in a depth $d=3$ Q-MPS. The comparative study of two models in Fig. 8 (a) shows similar performance (or slightly better). However, when we scale up the lattice size to $6\times 6$, the ground state energy obtained by Q-PEPS is much lower, indicating better scalability in solving two-dimensional lattice Hamiltonians. One should notice that by increasing the lattice size from $4\times 4$ to a $6\times 6$, the size of full Hilbert space is increased by a factor of $2^{20}$.

The batch size in the above simulation is $B=1024$, as the it grows, the standard deviation of measured gradients $\sigma_{s}$ will decrease as $\sim 1/\sqrt{B}$. We can also define the variance of gradients as

[TABLE]

where circuit parameters ${\boldsymbol{\theta}}_{i}$ is randomly sampled. If this values is too small, then the gradient vanishes. In a momentum based classical optimizer, the learning rate can be atomatically adjusted to utilize small gradients. However, what matters here is not the absolute value of gradients, but the ratio between the standard deviation of gradients $r_{\rm gs}=\sigma_{g}/\sigma_{s}$. When $r_{\rm gs}\ll 1$, the gradient information will be smeared out and the training will probably fail, then we have to either increase the number of shots $B$ or change the ansatz. Scalable training relies on how does $r_{\rm gs}$ vary as the system size grows. If we fix the batch size to $4096$, at a system size $4\times 4$, $r_{gs}=0.142/0.021$, while at a system size $6\times 6$, it is $0.088/0.032$, i.e. the gradient is decreasing and sampling error is increasing. Due to the limit of classical computational power, it requires a quantum device to verify the scalability and fully exploits its power.

IV Complexity analysis by Gate Counting

The wall clock time of both classical simulation and actual experiment can be estimated by counting the gate operations, as summarized in Table 2. In the classical simulation, the operation time $T_{\mathrm{gate}}$ scales exponentially with the number of virtual qubits $V$. We employ the differentiable quantum simulator Yao.jl Yao and its GPU backend CuYao.jl CuY in order to boost the classical simulation efficiency. With native CUDA programming support Bezanson et al. ; Besard et al. (2018), one can parallelize the circuit simulation both for the Hilbert space (of size $2^{V+1}$) and for the batch dimension (of size 4096 in our case). It takes $11-48$ hours (varies for various blocks shown in Fig. 3) for a typical parameter reported in this letter, $V=4,d=5$ and $N=16$, on a single Nvidia Titan V GPU card.

The gate time of an actual quantum device shows constant scaling with respect to the qubit number. For experiment on SQUID qubits, $T_{\mathrm{gate}}\approx 25ns$ Chow (2010) for a ${\rm SWAP}^{\alpha}$ gate. Thus, solving the $4\times 4$ Heisenberg model with circuit depth $d=1$ and the same hyperparameters as our numerical experiments will take approximately $18$ minutes on a single quantum processing unit (QPU) without considering readout time. On the other hand, the typical gate time for ion traps is $T_{\mathrm{gate}}\approx 0.5\mu s-250\mu s$ Linke et al. (2017); Schäfer et al. (2018), which means longer time is required for solving the same model.

Furthermore, we note that gate operations in the lines $2-5$ of Table 2 are also trivially parallelizable on QPUs. Therefore, we envision that building a cluster of QPUs Wecker et al. (2015a, b) may provide further advantage for high-throughput gradient estimation of the VQE calculation. In this case, one only needs classical communications to collect the gradients measured on all QPUs. Technically, having intermediate scale quantum circuits running in parallel is also easier than building a fully entangled large scale quantum computer. In this way, we expect running the variational algorithm on parallel QPUs will soon win over classical processors with exponential gate time.

A practical issue is that the total running time $MT_{\mathrm{gate}}$ should be less than the coherence time of the qubits, which limits the block depth of the variational ansatz to be shallow circuits on near-term devices. Table 3 shows the variational energy and fidelity for various depths $d$ obtained for the $N=4\times 4$ frustrated Heisenberg model at $J_{2}=0.5$. One sees that it is possible to reach fidelity $0.917$ with only $60$ parametrized gates, which is within the reach of the current day quantum technology. For larger problem size, one will need to increase the circuit depth linearly with $N$. Assuming area law entanglement entropy scaling of the system, one also needs to scale the circuit width $V$ linearly with the boundary size for an accurate variational description of the ground state.

A crucial step for the qubit efficient Q-MPS preparation scheme is the measure and reset operation. The cost of this step is device dependent. It is straightforward for trapped ions. However, for SQUIDs, a single qubit measurement and reset can take several microseconds, which is even slower than applying a gate. Fabricating low-latency quantum circuits that support fast measure and control is a rewarding direction in light of the proposed qubit efficient VQE scheme. Alternatively, one can employ the same circuit architecture without reusing the measured qubits. In this case, one still has the benefit of enhanced gradient signal.

V Discussions

Classical quantum many-body computation approaches provide valuable insights to quantum algorithms. Using the proposed qubit efficient VQE scheme, one can access the ground state properties of quantum systems using fewer qubits than the system size. Moreover, by exploiting the physical symmetries in the quantum circuit architecture design, one can alleviate the gradient vanishing problem and speed up the convergence to the ground state.

The use of conserved quantum number in circuit construction also allows one to have access to excited states in various quantum number sectors. With this regard, it is interesting to consider what are the universal gate sets with respect to various physical symmetry constraints. In addition to internal symmetries, the spatial translational symmetry may be taken into account via parameter sharing in the circuit blocks. Then, it naturally raises the question of whether one can study infinite large periodic systems with a finite number of qubits.

Another interesting direction is to perform time evolution or measure time-dependent quantities in the qubit efficient scheme. Since one does not have access to the full wave function directly in the qubit efficient scheme, the Trotter decomposition based time evolution Trotter (1959); Lloyd (1996) may not be directly applicable. Variational quantum algorithms for time evolution Li and Benjamin (2017); Mitarai and Fujii appears to be a good candidate for this purpose.

Acknowledgements.

We thank Miles Stoudenmire for sharing the idea of Ref. Huggins et al. (2019) prior to its publication. We thank Yun-Fei Pu and Ding-Shun Lv for providing valuable information on experimental feasibility. We thank Norbert Schuch for comment on the ${\rm U}(1)$ preserving circuit construction and Yan-Xia Liu for helpful discussions on Bethe Ansatz. We thank Pan Zhang for generous allocation of GPU hours and Xiu-Zhe Luo for contribution to Yao.jl Yao . The authors are supported by the National Natural Science Foundation of China under the Grant No. 11774398, the Strategic Priority Research Program of Chinese Academy of Sciences Grant No. XDB28000000 and the research funding from Huawei Technologies under the Grant No. YBN2018095185.

Appendix A Qubit efficient scheme for cluster state

We provide a concrete example of preparing and sampling a cluster state in the qubit efficient scheme. The one-dimensional (1D) cluster state Schön et al. (2007) can be prepared by circuit shown in Fig. 9. One first prepares the $N$ qubits to $\left|+\right\rangle$ by applying Hadamard gate on $\left|\mathbf{0}\right\rangle$ and then sequentially applies CZ gate between the nearest neighbor qubits. Thus, $N$-qubit 1D cluster state can be written as

[TABLE]

This 1D cluster state has an equivalent MPS representation with bond dimension two D’Souza (2014)

[TABLE]

where the left and right boundaries are A^{(1)}_{0}=\frac{1}{\sqrt{2}}\mbox{\left\langle+\right|},A^{(1)}_{1}=\frac{1}{\sqrt{2}}\mbox{\left\langle-\right|} and A^{(N)}_{0}=\mbox{\left|0\right\rangle},A^{(N)}_{1}=\mbox{\left|1\right\rangle}; while for the remaining ones $A^{(i)}_{0}=\frac{1}{\sqrt{2}}H,A^{(i)}_{1}=\frac{1}{\sqrt{2}}HZ\ (2\leq i\leq N-1)$.

One can prepare the same $N$ qubit 1D cluster state using only two qubits. The key is to mediate the entanglement between the physical qubits in different steps using one virtual qubit. As shown in Fig. 9 (b) we initialize both qubits to \mbox{\left|0\right\rangle}\otimes\mbox{\left|0\right\rangle}, then apply Hadamard gates to create a \mbox{\left|+\right\rangle}\otimes\mbox{\left|+\right\rangle} state. Applying CZ gates on them we can measure the physical qubit. We then reset the physical qubit to $\left|0\right\rangle$ and reuse it. The SWAP gate then exchanges its state with the virtual qubit. We then apply the Hardmard and CZ gates and measure the physical qubits. Repeating these steps to the end, we will obtain the same statistics of the measured bits as in the case of $N$ physical qubits Fig. 9(a).

All two-point Pauli correlation functions vanish on the 1D cluster state. The non-zero correlation functions will emerge if one performs the stochastic local operations and classical communication (SLOCC) Dür et al. (2000), e.g., measuring the $i$-th qubit in $Z$ basis and only considering the states with measurement outcome 0. In the general case, for an $N$-qubit 1D cluster state with an SLOCC operator $S$ (S=DHR_{z}(\gamma),D=\cos\theta\mbox{\left|0\right\rangle}\mbox{\left\langle 0\right|}+\sin\theta\mbox{\left|1\right\rangle}\mbox{\left\langle 1\right|}) acting on $i$-th qubit, second-nearest-neighbour correlation function $\langle\sigma_{i-1}^{z}\sigma_{i+1}^{z}\rangle=\cos 2\theta\sin\gamma$ D’Souza (2014), and other two-point correlations vanish. In order to measure the correlation function, e.g., $\langle\sigma_{2}^{z}\sigma_{4}^{z}\rangle$, the estimation protocol is: \raisebox{-.9pt} {1}⃝ measuring the $2nd$ and $4th$ qubits in $Z$ basis; \raisebox{-.9pt} {2}⃝ performing $HR_{z}(\pi/2)$ operation on the $3rd$ qubit, measuring it in $Z$ basis and only keeping the states with measurement outcome 0; \raisebox{-.9pt} {3}⃝ the measuring bases of the remaining qubits are randomly chosen for the final estimation. Repeat the above procedure until a good estimation of $\langle\sigma_{2}^{z}\sigma_{4}^{z}\rangle\approx{\mathop{\mathbb{E}}\limits_{\begin{subarray}{c}\end{subarray}}\left[\sigma_{2}^{z}\sigma_{4}^{z}\right]}$ is obtained.