The theory of variational hybrid quantumclassical algorithms
Abstract
Many quantum algorithms have daunting resource requirements when compared to what is available today. To address this discrepancy, a quantumclassical hybrid optimization scheme known as “the quantum variational eigensolver” was developed Peruzzo et al. (2014) with the philosophy that even minimal quantum resources could be made useful when used in conjunction with classical routines. In this work we extend the general theory of this algorithm and suggest algorithmic improvements for practical implementations. Specifically, we develop a variational adiabatic ansatz and explore unitary coupled cluster where we establish a connection from second order unitary coupled cluster to universal gate sets through relaxation of exponential splitting. We introduce the concept of quantum variational error suppression that allows some errors to be suppressed naturally in this algorithm on a prethreshold quantum device. Additionally, we analyze truncation and correlated sampling in Hamiltonian averaging as ways to reduce the cost of this procedure. Finally, we show how the use of modern derivative free optimization techniques can offer dramatic computational savings of up to three orders of magnitude over previously used optimization techniques.
I Introduction
Eigenvalue and more general optimization problems lie at the heart of applications and technologies ranging from Google’s Page Rank and aircraft design to quantum simulation and quantum chemistry Page et al. (1999); Cullum (1994); Golub and van der Vorst (2000). Quantum computers promise to provide ground breaking advances in our ability to solve these problems by offering solutions that may be exponentially faster than the classical equivalent in some cases. However, delivering on these promises may require overcoming considerable technological challenges.
Since the initial proposal by Richard Feynman Feynman (1982), a number of advances have been made in understanding how to use a quantum computer to help solve eigenvalue and optimization problems. The quantum simulation algorithms of Abrams and Lloyd Abrams and Lloyd (1997, 1999) showed how eigenvalues corresponding to some Hermitian operator could be extracted from eigenvectors exponentially faster with respect to dimension than the classical equivalent. Leveraging this idea, AspuruGuzik et. al. showed how one could perform exact quantum chemistry computations in polynomial time for some instances, pushing the boundaries of predictive quantum chemistry AspuruGuzik et al. (2005). These ideas have since been tested successfully in proofofprinciple quantum experiments using architectures such as quantum photonics, nitrogen vacancies in diamond, and ion traps Lanyon et al. (2010); Lu et al. (2011); AspuruGuzik and Walther (2012); Peruzzo et al. (2014); Wang et al. (2015).
In recent years, there has been a growing interest in the particular application of quantum chemistry on quantum computers. As a result, a number of efforts have been made to study the scaling and performance of various algorithms while simultaneously offering dramatic algorithmic improvements Wang et al. (2008); Whitfield et al. (2011); Kassal et al. (2011); Seeley et al. (2012); Wecker et al. (2014); Kais (2014); Hastings et al. (2015); Poulin et al. (2014); McClean et al. (2014); Babbush et al. (2015); Whitfield (2013, 2015); Babbush et al. (2015a, b); Veis and Pittner (2014); Toloui and Love (2013); Trout and Brown (2015); Huh et al. (2015). The original proposal of quantum chemistry on a quantum computer also introduced the idea of adiabatic state preparation, closely related to general adiabatic quantum computation. A number of advances in this field as well as extensions of adiabatic computation concepts to more general optimization problems have arisen as well Neven et al. (2008); Veis and Pittner (2014); Babbush et al. (2014).
Unfortunately, despite developments in quantum algorithms and optimization of resource requirements, many of the algorithms have hardware requirements far beyond the capability of nearterm quantum computers. Moreover, the overhead of some asymptotically optimal algorithms is such that even the first quantum computers competitive with classical supercomputers may not be able to run them. To this end, in 2014 Peruzzo and McClean et al. developed the variational quantum eigensolver (VQE), a hybrid quantumclassical algorithm designed to utilize both quantum and classical resources to find variational solutions to eigenvalue and optimization problems not accessible to traditional classical computers Peruzzo et al. (2014). This algorithm was originally implemented and tested on a photonic quantum chip and has since been extended both theoretically and experimentally to ion trap quantum computers Yung et al. (2014); Shen et al. (2015).
The VQE has the notable property that it can run on any quantum device, making it a candidate for exploring the performance of early quantum computers. Moreover, the algorithm is designed to take advantage of the strengths of a given architecture. That is, if some gates or quantum operations may be performed with higher fidelity, then the algorithm can leverage these strengths in the design of the quantum hardware ansatz. Perhaps one of the most interesting features of the algorithm is its ability to variationally suppress some forms of quantum errors, which is discussed later in this work. This intrinsic robustness to quantum errors in combination with low coherence time requirements has placed this algorithm as a potential candidate for the first to surpass a classical computer, using a prethreshold quantum device. Even in the event that some error correction is required to exceed current computational capabilities, this same robustness may translate to requiring minimal error correction resources when compared with other algorithms.
In this work we aim to present the hybrid quantumclassical variational approach in more detail, offering both theoretical and practical exposition on developments since the original hybrid quantumclassical proposal. Additionally, although a strength of the variational quantum eigensolver is its ability to adapt to the given hardware, this work will be the first to analyze VQE in the abstract, in a way that is completely general to any quantum device. We begin by reviewing background and notation as well as the outline of the variational quantum eigensolver algorithm. This is followed by a discussion of ansatz states that allow one to explore classically inaccessible regions of Hilbert space, including a variational formulation of adiabatic state preparation and unitary coupled cluster. We then explore how this approach may be used to variationally suppress certain types of quantum errors. Following this, we introduce several computational enhancements to the Hamiltonian averaging method for obtaining expectation values, including the truncation of unimportant terms and grouping terms by commutation and covariance. These enhancements are able to considerably reduce the cost of the procedure. Finally we cover aspects of the classical optimization procedure associated with the VQE and show how modern derivativefree optimization technique have the potential to greatly enhance the efficacy of the method.
Ii Background and Notation
ii.1 General Quantum Systems and the Variational Principle
Let us consider a quantum system composed of qubits which will act as our quantum computer, and a Hamiltonian of a different system that need have no relation to other than acting on a space of qubits. This Hamiltonian could be derived from a physical system such as a collection of interacting spins or the discretization of an interacting electronic system. Similarly it could come from the encoding of an optimization problem or the problem Hamiltonian in adiabatic quantum computation. In all of these instances, one is interested in the eigenvectors and eigenvalues, , of the Hamiltonian , and the goal will be to find and study these eigenvectors and eigenvalues using .
In the VQE approach, the eigenvectors are encoded by a set of parameters that can be used to prepare them on demand when other observables are desired. We order the eigenvectors by the eigenvalues such that . Indeed in many cases, the eigenvectors corresponding to the lowest few eigenvalues and their properties are of primary interest. In physical systems this is because lowenergy states play a dominant role in the properties of the system at modest temperatures, and in optimization problems they often encode the optimal solution.
Recall the expectation value of an operator with respect to a state
(1) 
We will assume normalization of the wavefunction, , for the remainder of the work, however attention should be paid in the case of leakage errors from the computational basis. Our attention is restricted to the class of operators whose expectation value can be measured efficiently on and mapped to . A sufficient condition for this property is that operators have a decomposition into a polynomial sum of simple operators as
(2) 
where is an operator than acts on , runs over a number of terms polynomial in the size of the system, is a constant coefficient, each has a simple measurement prescription on the system . This will allow for straightforward determination of expectation values of on by weighted summation of projective measurements on the quantum device . A simple example of this is the decomposition of a Hermitian operator into a sum of tensor products of Pauli operators weighted by constant coefficients.
Consider a set of real valued parameters , which we arrange into a vector , and the Hamiltonian of . If one prepares into a quantum state depending on these parameters, , then the variational theorem of quantum mechanics states that
(3) 
As a result, the optimal choice of to approximate the ground state (or eigenvector corresponding to the lowest eigenvalue) is the choice which minimizes . Note that the state is normalized for all choices of by the unitarity of quantum evolution or trace preservation under quantum operations in state preparation. The variational principle also extends to other eigenstates. If one has constructed an ordered orthonormal set of approximate eigenstates , such that , then
(4) 
where are the ordered true eigenvalues of the operator . Thus, repeated application of the variational principle under orthogonality constraints can yield an approximation to as much of the spectrum as desired, incurring additional cost for each eigenvalue. Alternatively, one can perform a spectral transform to the Hamiltonian and use the groundstate variational principle to find excited states, as in the folded spectrum method Wang and Zunger (1994). That is, minimize where and is some real parameter. In the transformed Hamiltonian, the ground state corresponds to the eigenvalue in the original Hamiltonian closest to .
More generally, the state preparation scheme may be influenced by an environment and would be better represented by an ensemble given by a density matrix . In an ideal scenario where the preparation is error free and a pure state is maintained, . In the density matrix formalism, the expectation value of an operator is given by
(5) 
and the ground state variational principle on the Hamiltonian still holds such that for any approximate density matrix , and for all choices of ,
(6) 
As a result, the optimal choice of to approximate the ground state is that which minimizes . The fact that this principle still holds for mixed states has important consequences for the robustness of the method to errors and environmental influence. By finding the set of parameters that minimizes the energy, one is in effect, finding a set of experimental parameters most likely to produce the ground state on the average, potentially affecting a blind purification of the state being produced. This ability to suppress errors without knowledge of the mechanism will be elaborated upon later in this work.
Another important quantity is the variance of an operator with respect to a state. For an operator and a general mixed state , this is given by
(7)  
(8) 
A variational principle on the variance exists as well, and has been used extensively for optimization in the context of quantum Monte Carlo Lester et al. (1994). Note that for any eigenstate of an operator , the variance is given by
(9) 
and for any approximate eigenstate , we have that
(10) 
ii.2 Fermionic Hamiltonians and Quantum Chemistry
While the VQE and its principles can be applied to general quantum problems, an application of particular recent interest is that of quantum chemistry and fermionic Hamiltonians. Given a set of nuclear charges and a number of electrons, the standard form of the electronic structure problem is to solve for the eigenvectors and eigenvalues of the electronic Hamiltonian , written as
(11) 
where atomic units have been used, are nuclear positions, electronic positions, and are nuclear masses. Due to large separations in the nuclear and electronic masses, an excellent approximation to this problem at the time and energy scales of chemical interest is to treat the nuclei as classical point charges under the BornOppenheimer approximation with fixed positions . The problem as written is referred to as the first quantized representation of the quantum chemistry problem. A number of algorithms have been developed for quantum computers to treat the problem directly within this framework Kassal et al. (2008); Toloui and Love (2013); Welch et al. (2014), however the focus in this work will be on the second quantized treatment.
To reach the practical form of the second quantized Hamiltonian, one must project the problem into a finite, orthogonal, spinorbital basis, of which we will denote members , and impose the requirements of fermion antisymmetry through the fermion creation and annihilation operators and . With these steps, the second quantized Hamiltonian takes the form
(12) 
with coefficients determined by the spin orbital basis as
(13)  
(14) 
where describes both the spatial position and spin of an electron as . The operators and obey the standard fermion commutation relations as
(15)  
(16) 
A crucial part of solving these problems on quantum computers is the mapping from fermions to qubits. The two most common mappings under current study are the JordanWigner transformation Jordan and Wigner (1928); Somma et al. (2002) and the BravyiKitaev transformation Bravyi and Kitaev (2000); Seeley et al. (2012); Tranter et al. (2015). In the case of the JordanWigner transformation, the mapping from fermion operators to qubits is
(17)  
(18)  
(19) 
ii.3 Reference States
Many traditional methods for electronic structure involve the concept of a reference state. A reference state is a product state that is used as a starting point to define a more general quantum state, and can allow for great formal simplification. Here we will briefly introduce why they are convenient and useful, and then how they are obtained.
An example spin reference and fermion reference state might be the general product states
(20)  
(21) 
where is the fermion vacuum state, is the number of sites a fermion can occupy, and is the number of qubits or fermions. Even though these are separable product states, their manipulation theoretically or preparation on a quantum computer can be cumbersome as written. However, because they are product states, there exist efficient, local unitary basis transformations and such that these states can be rotated into a simple form with weight on a single computational basis state. That is
(22)  
(23) 
and because the transformations are local, the transformation of the Hamiltonian to the new basis such that the physical problem remains unchanged is also efficient. In the case of quantum chemistry, this corresponds to a transformation of the integral terms and , which may be computed in a time O exactly.
These new simpler forms of the state have advantages both in theoretical manipulation, and in ease of preparation with quantum resources. For example, the preparation of the untransformed spin reference state could require at least local rotations, not including error correction on a quantum device to prepare from a computational basis state, whereas the new reference is simply the computational basis state from which most computations begin. Here we have traded modest classical effort in transforming the basis of the Hamiltonian for savings in quantum resources.
These reference states are typically obtained from mean field calculations, which are guaranteed to have product states, such as those given above, as solutions. In chemistry, this procedure is called HartreeFock, and the transformation of the state to the simplified form is known as the canonical condition in the solutions of the HartreeFock equations, resulting in the canonical molecular orbitals.
When the problem is well treated by meanfield theory, it can be shown through perturbation theory that the dominant corrections to the meanfield solution are given by quantum states “close” to the meanfield solution in the sense of fermion excitations Helgaker et al. (2014) or Hamming distance. This is the origin of the perturbative MP2 method, configuration interaction, and coupled cluster methods Bartlett and Musiał (2007); Helgaker et al. (2014), which all solve the problem close to a given reference and have been applied to both electronic and frustrated spinsystems Götze et al. (2011).
In some problems, particularly when correlation is strong, the meanfield description is a poor starting point for the problem. In this case, one may still use a referencelike formalism, but starting with an entangled state. These methods are called multireference methods in quantum chemistry Helgaker et al. (2014); Laidig and Bartlett (1984); Musiał et al. (2011), and carry considerably more theoretical and computational challenges with them. In this work, we will highlight how the generalization of methods on a quantum computer to the multireference case is often more natural than in the classical case.
ii.4 Algorithm Outline
To use a variational methodology to find approximations to the eigenvalues and eigenvectors of the Hamiltonian in a quantum computer, it is convenient to break the task into three distinct pieces and outline the algorithm very coarsely as

Prepare the state or on the quantum computer, where can be any adjustable experimental or gate parameter.

Measure the expectation value

Use a classical nonlinear optimizer such as the NelderMead simplex method to determine new values of that decrease

Iterate this procedure until convergence in the value of the energy. The parameters at convergence define the desired state.
In the coming sections we will elaborate on what is known about each of these steps and offer new algorithmic and conceptual improvements.
Iii State parameterization and preparation
The set of states a quantum computer can easily manipulate that a classical computer cannot is not yet fully understood Mora and Briegel (2005); Gross et al. (2009); Cai et al. (2015). Given the set of parameters , it’s clear that in order for a quantum computer to have an advantage, one would like the state to be good at describing the solution of interest, while also difficult to prepare and/or sample from classically using currently known methods. Here we will first discuss topics relevant to state preparation for all classes of states in the variational quantum eigensolver, independent of any notion of how difficult they are to prepare classically. We will then discuss some details concerning two classes of states currently believed to be both good at describing systems of interest and difficult to prepare and/or sample from classically, namely adiabatically parameterized states and (multireference) unitary coupled cluster states.
iii.1 Error bounds and distributions
Once a state has been prepared as a function of some set of parameters , one would like to know how close this state is to the solution of the problem being solved. In this work, we will say a measured value is known to precision based on a normal distribution approximation with standard deviation , which is reasonable given that most of our estimates will be derived from sums of random variates with finite variance, which by the central limit will rapidly converge to a normal distribution.
Suppose, for now, that the goal is to know an eigenvalue of to within a specified precision . Let be the eigenvalue of closest to . Under these assumptions on the eigenvalue the Weinstein inequalities Weinstein (1934); MacDonald (1934) hold
(24) 
As a result, a sufficient condition is to rigorously achieve the precision requirement on the eigenvalue is
(25) 
where as one approaches an eigenstate, the variance approaches . When considering only the ground state, one can derive a simple bound on the quality of the state. More specifically, in the zero variance limit, if has multiplicity , then the eigenstate corresponding to is reproduced as well. That is, if a bound on the gap to the first eigenstate is known in addition to the variance, such that , and , and we decompose the state into its eigenstate representation then we can quantify the quality of state preparation as a function of the measured variance
(26) 
For general excited states , one may find a similar bound exists based on a measurement of the variance of the operator and a known bound on the gap , such that
(27) 
where , and both bounds given here are derived in this appendix. If one has prior knowledge that a single eigenstate dominates the expansion, such that , and a lower bound , then Delos and Blinder Delos and Blinder (1967) showed through the method of moments that a tighter lowerbound on the eigenvalue is given by
(28) 
These bounds may be used to estimate the absolute accuracy the minimization procedure obtained within the given basis and decide if the eigenvalue has been determined to the desired accuracy and precision or if the state ansatz should be altered to adjust the cost or accuracy of the procedure.
iii.2 Adiabatically parameterized states
One type of quantum state that can be explored as a parametric ansatz is that produced by adiabatic state state preparation with a variable path. In adiabatic quantum computation Farhi et al. (2000, 2001); Boixo and Somma (2010) and adiabatic state preparation AspuruGuzik et al. (2005); Veis and Pittner (2014) one makes use of the adiabatic theorem Born and Fock (1928), which states loosely that if one prepares the lowest eigenstate of an initial Hamiltonian , by continuously changing the Hamiltonian from to a final problem Hamiltonian , one finishes in the lowest eigenstate of if the evolution was slow enough. In adiabatic computation, slow enough is quantified relative to the minimum eigenvalue gap between the ground and first excited states along the evolution. While many developments have occurred in the area of adiabatic quantum computation and modifications to the Hamiltonian, perhaps the most commonly considered form of evolution is defined by
(29) 
where , , and . The evolution is controlled by continuously changing the parameter as a function of time .
Consider the set of all paths of and from to as a function of time , and denote it , where is some finite time. Label one such path as . In a noiseless coherent situation at K, the unitarity of evolution dictates that the final state of the evolution is uniquely determined by the path . In this situation, we may write the final pure state as a higherorder function of the path , or . Thus any expectation values of the final state may be written as functionals of the path, , and by the variational principle
(30) 
such that the optimal path is the path in that minimizes the value of . This functional minimization may be changed into a standard minimization by parameterizing the path by a set of parameters , and performing an optimization on the parameters that determine the path. As such, adiabatic state preparation may be considered as an ansatz to be used in the variational hybrid quantumclassical approach, where the state parameters are the shape or nature of the path. The idea of refining the adiabatic path has been used before in the context of local adiabatic evolution Roland and Cerf (2002) with great success. The idea here is to achieve similar benefits in an entirely blackbox manner, guided only by a variational principle and measurements of the final point of the evolution.
As a simple example, consider a linear path in defined by a single parameter that controls how quickly the evolution is performed
(31)  
(32) 
and the parameter is restricted by membership in to . In the case of an ideal evolution with enough quantum resources such that the evolution is much longer than required by the problem gap, the adiabatic theorem implies that is optimal at the extremal point . Moreover, in the limit that , the adiabatic theorem implies that for any finitely gapped problem contains a path that prepares the exact ground state, and even the simplest linear paths, which are a subset of , are sufficient to do so.
Within this simple example, it’s not immediately clear why one would want the flexibility offered by the variational quantum eigensolver formulation, as one could choose the linear path with minimal without the need for any optimization of . However, a more realistic situation may be such that is smaller than the required time of evolution dictated by the problem gap, due to technological constraints or simply human time constraints in a hard problem. It might also be possible that no good estimate of the gap is known, and one must attempt several paths regardless to establish confidence that the evolution is not too fast to impair accuracy. One should exercise caution in such attempts however, as the probability of success does not necessarily increase monotonically with evolution time, especially when one is far short of the time required by the problem gap or when errors are present Bookatz et al. (2014). Moreover, it is known that for systems experiencing decoherence or dephasing on the timescale of evolution that the slowest possible evolution is not optimal in preparing the ground state of the final problem Hamiltonian Steffen et al. (2003); Åberg et al. (2005); Crosson et al. (2014). In all situations, the final density matrix is determined by the parameters of the path, such that determines a density matrix , and an optimal choice of parameters can be made without detailed knowledge of the gap or errors present in a system by minimizing as a function of .
The Hamiltonians may also be generalized to include intermediate operators Farhi et al. (2002); Hofmann and Schaller (2014); Crosson et al. (2014); Zeng et al. (2015) such as
(33) 
where one considers any number of intermediate Hamiltonians and with . The set of paths satisfying these boundary conditions with available intermediate Hamiltonians , , offers more flexibility, and again a guiding principle to select parameters defining the optimal paths is given by the variational principle.
From this discussion it is clear that adiabatic state preparation where the path of evolution is defined by some set of parameters is one choice of parametric ansatz for the variational quantum eigensolver. It can be inferred from the known capabilities of adiabatic quantum computation that this ansatz is capable of preparing states that cannot be efficiently prepared or sampled from classically using only a small number of parameters with currently known methods Aharonov et al. (2008). As seen in the simple linear example, the number of parameters to meet this condition may be as few as for a linear interpolation that is slow enough in ideal conditions.
Variational Adiabatic Path Example
To further illustrate the utility of a variational perspective on adiabatic quantum computational methods in a resource constrained setting, we consider here a simple 1qubit problem first studied in the adiabatic context in the original work of Farhi et al Farhi et al. (2000). In particular, we will consider this problem in a resource constrained context where the maximum evolution time is limited. In this problem, the Hamiltonian the initial and problem Hamiltonians are given by
(34)  
(35) 
If we take the following form of the schedule Hamiltonian
(36) 
then the eigenvalues of this problem undergo an avoided crossing with a gap determined by the size of the perturbation . For this example we choose and the resulting spectrum is plotted in Fig.1 as a function of . Suppose that we are attempting to prepare the ground state of our problem Hamiltonian in a situation where the total evolution time is limited.
We will consider two types of paths, the first of which is a fixed standard linear path as a function of time. That is with . The second type of path will be a parameterized path of two variables defined by the best cubic Bspline fit of the 4 points , where the the parameters are determined by a nonlinear minimization the expectation value of the final state in the (possibly non)adiabatic evolution with fixed maximum evolution time, . In this simple example we use the NelderMead simplex method to perform a derivative free optimization of , in analogy to how it might be performed on a quantum device. We use as an initial condition and in the optimization, which corresponds to the linear path.
The resulting variationally optimal adiabatic spline path is plotted alongside the standard linear path in Fig. 2, which shows that the method naturally finds a path which slows evolution near the closing gap, without any prior knowledge of the spectrum, and only measurements at the endpoint as opposed to the entire path. The effect of this on the success of preparing the ground state as a function of the total available evolution time is shown in Fig. 3. From this figure we observe that the variationally optimal adiabatic spline path is able to achieve similar results to a linear path with roughly times less evolution time. That is, at the cost of some classical minimization, we have reduced the quantum evolution time requirement by a factor of by slightly deforming the schedule in a blackbox manner relying only on measurements of the final state of the evolution and no prior knowledge of the problem. Moreover, even at this reduced evolution time, we achieve the desirable property that the success of the computation is a monotonically increasing function of , which is not true of the linear schedule in this case.
Pontryagin’s Principle and NonAdiabatic BangBang Quantum Computation
While adiabatic evolution or attempted adiabatic evolution is one way to prepare a desired state, it is certainly not the only option. Nonadiabatic evolution opens a different class of potential schedules for preparing a desired state guided by the variational principle. The form of the schedule Hamiltonian has a particularly interesting form, namely that it is a linear evolution problem with a control that effects a linear coupling. In the theory of optimal control, it is known through application of Pontryagin’s minimization principle that the optimal control setting for reaching a desired state of the controlled system when the system has a linear coupling to the control is to have the control at its extremal values Bellman et al. (1965). That is, becomes a sequence of step functions where it takes the values or and need not satisfy the previous boundary conditions and . This class of solutions to optimal control problems is known as a “bangbang” solution, and is obviously nonadiabatic by construction. This principle has been shown in quantum optimal control outside of the context of quantum computation, where a Monte Carlo minimization scheme was applied to determine the schedule of step functions, and a different variational principle was employed Rahmani et al. (2013). However this scheme could be straightforwardly adopted using the variational principle methods described here to engineer state preparation schedules for a state of interest, or to perform more general quantum computation.
iii.3 Unitary coupled cluster
Another method to parametrically explore the Hilbert space of possible quantum states is the unitary coupled cluster method developed in quantum chemistry Taube and Bartlett (2006); Bartlett and Musiał (2007). The projective nonunitary (and nonvariational) form of these equations form the basis for the goldstandard of classical quantum chemistry, coupled cluster with single and double excitations with perturbative triple excitations [CCSD(T)] Raghavachari et al. (1989); Bartlett and Musiał (2007) and has its origins in nuclear physics Coester and Kümmel (1960). The unitary form of these equations do not have a well defined truncation as the projective form does, and one must rely on perturbative arguments to handle the BCH expansion that break down when the parameters defining the states grow. This ansatz for electronic systems has been documented in classical quantum chemistry and in previous works on the variational quantum eigensolver Taube and Bartlett (2006); Bartlett and Musiał (2007); Peruzzo et al. (2014); Yung et al. (2014), and here we document its generalization to generic collections of interacting twolevel quantum systems, which include the antisymmetric electronic case as a specialization. We note that coupled cluster has been utilized before in the context of frustrated spin systems such as Kagome lattices Schmalfuß et al. (2006); Götze et al. (2011), but our treatment will extend beyond a fixed reference and also focus on the unitary variant of the method.
To conceptually introduce the approach, recall the introduction of reference states earlier in this work, and consider a single computational reference state of an qubit quantum system, . One way to parametrically explore Hilbert space is to consider the space of states “close” to in the sense of Hamming distance or bit flips. This method, sometimes called configuration interaction (CI) or state space restriction enumerates available states through the use of spin flip Helgaker et al. (2014); Kuprov et al. (2007). For example, all states flip away from may be written as
(37) 
where in this case are complex coefficients and is the qubit raising operator applied to qubit . This expansion can be extended systematically by including multiqubit spinflip operators to eventually parametrize all states in the Hilbert space, or full configuration interaction (FCI). While this parametric construction of states is straightforward, it has a number of deficiencies that render it nonoptimal. We will not attempt to explore all of those here, and note only that this ansatz is efficient to prepare and use classically for any truncation to a fixed number of spin flips , and it is not clear that there is an advantage to specifically preparing a linear truncated state on a quantum device.
An idea closely related to this is coupled cluster, which also uses the spinflip concept to explore states “close” to a reference, but as a generator used in exploration of the space. In the case of quantum computing, its unitary variant is of particular interest, as unitary state preparation is a natural operation on a quantum computer. Conventional implementations of coupled cluster often utilize a single, well defined reference state with all spins aligned, i.e. . With this assumption, one may explore all of quantum space through successive flips in the computational basis. As a simple example, if one is interested in only real wavefunctions, the space of single spin flips may be explored by
(38) 
and successively larger fractions of the space of real wavefunctions may be covered by introducing multiple spin flips. In the study of general quantum states however, it is sometimes necessary or more efficient to explore quantum state space from an arbitrary reference , which could be entangled or simply more complex than . These challenges have been studied in the context of multireference coupled cluster in quantum chemistry Laidig and Bartlett (1984); Musiał et al. (2011). Moreover in quantum computation one may not have perfect knowledge of the reference state, nor want to require it in their algorithm. For example the reference state could be prepared by some adiabatic state preparation procedure. In this situation one could accidentally have as a reference state with , from which no state exploration is possible with the above cluster operator. The space of nontrivial single qubit operators is spanned by . As such we want to generalize to a set of antiHermitian operators spanning the same space, given by
(39)  
(40)  
(41) 
For convenience we have introduced the standard Pauli operators in the numerical indexing scheme, that is , , , . As one is not typically interested in global phase factors, we implicitly ignore the identity operator in all equations going forward and with the remaining operators we may write the first order cluster operator as
(42) 
where are real, Roman indices indicate different qubits, and the Greek indices indicate different Pauli operator bases. More generally the ’th order cluster operator may be written as
(43) 
where , is a index tensor containing the variational parameters, and the full cluster operator up to order is written
(44) 
From this general cluster operator, we define the unitary coupled cluster state of order with reference as
(45) 
With this exposition it becomes clear that unitary coupled cluster generators for a totally general spin reference case at order are the antiHermitian algebra and the set of possible actions on the qubits are all possible unitary transformations on qubits that leave the global phase unchanged, or .
This represents a parametric state preparation with real parameters. While this has the potential to represent any known quantum operation at sufficient order and precision of implementation, practically speaking one often restricts to the case of , which has been found to be quite powerful in expressing states in quantum chemistry. This represents a powerful ansatz with a number of parameters that grows only quadratically in the size of the system. Additionally, the state preparation is manifestly unitary by construction, and has no known efficient classical preparation or method for sampling with arbitrary (possibly entangled) reference . As has been noted previously, this state can be prepared efficiently for any fixed order to a specified accuracy on a quantum device by using the SuzukiTrotter factorization of the unitary operator Peruzzo et al. (2014); Trotter (1959); Suzuki (1993). We note that as one is not trying to faithfully reproduce some dynamics as in many uses of the SuzukiTrotter factorization, that a coarse factorization may suffice, simply altering the definition of the ansatz, but still remaining difficult to simulate classically.
As an extension to the suggested implementation of spin unitary coupled cluster by SuzukiTrotter, one may use the connection to to take a more geometric approach and explore states through geodesic constructions as was done by Nielsen et al. Nielsen et al. (2006). Moreover if one allows values of different parameters at different Trotter steps, one may perform arbitrary 1 and 2 qubit gates at , which forms a universal gate set and the ansatz can be made equivalent to an arbitrary quantum circuit with a sufficient number of Trotter steps. To see this, consider the first order in a Trotter factorization with a second order cluster operator and a Trotter number of . One could prepare the desired state from a given reference as
(46) 
where we emphasize that it is more correct to consider the use of the exponential splitting as a redefinition of the ansatz than an approximation. Instead of following this precise splitting procedure, where the same parameters are used in each Trotter step, one can relax the parameters to have independent values at each time step, and to not split Pauli operators acting on the same two qubits within one time step. This results in an ansatz of the form
(47) 
The operator defined by
(48) 
can express an arbitrary element in and thus its exponential can be used to form an arbitrary two qubit gate on any two qubits, or said differently, an arbitrary element of on any two qubits. Arbitrary two qubit gates on any qubit are known to constitute a universal gate set Nielsen and Chuang (2010), and then clearly can be used to construct any desired universal gate set such as the CliffordT set. This establishes a clear connection between second order unitary coupled cluster and universal quantum computation through relaxation of parameters in an exponential operator splitting. This also opens the research direction of connecting states of this type to tensor networks where the network is defined by the action at each “timestep” of unitary coupled cluster Schwarz et al. (2012).
iii.4 Fermionic UCC
Due to particular interest in the quantum chemistry and other fermionic problems, it is worth discussing the specialization of this method to those cases. First taking again the case of a fixed computational reference, such as , in analogy to the spin case, the first and second order cluster operators conventionally take on a simple form, that is
(49)  
(50) 
with indexing the occupied orbitals, indexing the unoccupied orbitals, and higher orders defined in the obvious way of including more excitation operators. These generators are constructed to conserve particle number at all orders and parametrically depend on real parameters at order . In the case of a single reference, it should be noted that all the excitation operators commute as a direct consequence of the creation and annihilation operators being restricted to act on different subspaces. As a result, Trotter factorization of this ansatz may be performed to arbitrary times exactly that allows one to explore regimes where low order truncations of the BCH expansion are not accurate and thus may be difficult to sample from classically.
We can understand the equivalent action on qubits by mapping the fermion operators to spin operators via either the JordanWigner or BravyiKitaev transformations discussed earlier in this work. In the case of the JordanWigner mapping, as a result of the nonlocality of these mappings, at every fermion order , we find spin flips up to all spins and observe that the allowed operations on the qubits are a nontrivial subgroup of at every order . This demonstrates that it is key to develop the ansatz in the fermionic framework before mapping the problem to a spin representation. If one were to first map to spins, then use the spin coupled cluster formulation, the ansatz might explore many irrelevant or symmetry broken states, such as mixtures of different particle number states.
In analogy to our exposition on spins however, this type of cluster operator is reference state specific. That is, there are some reference states from which it will fail to parameterize the entirety of the fermion space and extensions to multireference states can require a different cluster operator for each reference. This can be seen from dimension counting in the vector space of the fermion excitation operators. For example at first order these operators only span a real vector space of dimension whereas the full space of all fermion linear operators has real dimension . In classical implementations of multireference coupled cluster there are many different approaches to solving this and related problems going by names such as “universal” or “state selective” multireference coupled cluster Laidig and Bartlett (1984); Musiał et al. (2011). In the case of unitary coupled cluster on a quantum computer, in analogy to how we generalized the distinguishable spin operators, we can generalize the fermion operators to treat arbitrary references without such concerns.
The operators and their tensor products where and run over all orbitals (instead of restricting them to occupied and unoccupied relative to a reference) form a basis for the real vector space of operators on fermion states. As a result, to allow arbitrary action on the space of fermions, the span of the generating operators used must match this. To span the same real vector space as these operators we use the following antiHermitian basis
(51)  
(52) 
and all possible fold tensor products of these operators. One can verify by dimension counting of the real vector space that these operators in fact span the entire space of possible fermion operators. With these operators, the first order fermion cluster operator can be written as
(53) 
where and run over all orbitals and indexes the antiHermitian fermion generators. Higher orders of the cluster operator can be built naturally from tensor products of these operators, such that at the ’th order we have
(54) 
where the same vector operator shorthand as the spin case has been used. With this construction the power of the cluster operator is state agnostic, and fermion number conserving. We term this the state universal quantum unitary coupled ansatz (SUQUCA). Again, in all cases the optimal choice of the parameters is determined through the application of the variational principle with respect to the Hamiltonian of interest.
iii.5 Quantum Error Suppression and Symmetries
A variational hybrid quantumclassical is designed to perform on prethreshold computers, where gates may be imperfect and random bit flip or phase errors may be introduced into the computation. Fortunately the variational formulation allows one to suppress certain types of errors naturally, which we will discuss here in the context of variational error suppression.
In the design of a parametric wavefunction ansatz, it is common to enforce known symmetry requirements for both theoretical and practical purposes. For example, in the fermionic unitary coupled cluster wavefunctions, the ansatz is designed to conserve the number of particles for all possible choices of the parameters . That is both the ansatz and the Hamiltonian commute with the number operator . While we haven’t explicitly done so here, it is also possible to adapt the cluster operators to conserve total spin Helgaker et al. (2014). In a fully error corrected quantum computer, this introduces no additional concerns and can simplify the problem under consideration. However in a prethreshold device or any with only partial error correction this must be taken into consideration.
Consider the preparation of an ansatz from some initial state, which we denote as . In a prethreshold, nonerror corrected quantum device, there can be a distinction between the formal specification of the ansatz preparation as a gate or operation sequence and the operation sequence actually performed on the system with inputs , which we will denote . We call an error in such an implementation suppressible if there exists a correction input vector such that for a specified , and further denote it variationally suppressible if the corrected vector also corresponds to an optimum on the parameter surface. In such a case, the variational quantum eigensolver can suppress these errors naturally without detailed knowledge of the error mechanism.
A troublesome nonsuppressible case is when an error violates a symmetry of the ansatz. More explicitly, if we denote the symmetries of the ansatz as the set of operators such that for all , then for any symmetry violating error such that , there does not exist any correction vector such that the desired preparation can be performed.
To be more concrete, consider the two examples given in this section, parameterized adiabatic state preparation and coupled cluster. In these cases, some symmetries of the ansatz can be trivially determined by the generating operators. In adiabatic state preparation, the symmetries will be given by the set of operators such that for all Hamiltonians , including the initial, problem, and intermediate Hamiltonians. In the case of coupled cluster, this will be the set of operators such that for all excitation type operators , such as the number operator. These represent sufficient conditions for for every possible choice of . In the case of fermionic coupled cluster, the generating operators are specifically designed to conserve particle number, such that one symmetry of the system is the number operator . In a JordanWigner qubit representation, this simply counts the number of qubits in state . As such, if a random error of the form is acted on any qubit, this error is not suppressible.
This particular error can be made suppressible by extending the set of generating operators to include spin flips (e.g. and ) or fermionic nonnumber conserving operators, e.g. and as well as all tensor products of these operators with the rest of the generating set. With the addition of these operators, this error become suppressible, however the error will only be variationally suppressible if the desired symmetry state of the ansatz corresponds to an energetic minimum. In the event that it does not, one can construct an auxiliary Lagrangian of the form
(55) 
where are penalty multipliers and are constants corresponding to the desired expectation values of the operators . In order to be efficient, measurements corresponding and must be also be efficient. Using this construction, one may minimize with respect to expectation values instead of , and in the limit that the symmetries will be exactly preserved while allowing variational error suppression under action by the extended operator set.
This methodology also allows for access to excited states that correspond to an energetic minima of a given symmetry. An example of this could be the lowest triplet energy state of a molecule with a natural singlet ground state, or the ionic state of a molecule after photodissociation. Use of this construction may allow easier access to these particularly important excited states, as compared to a more general excited state approach.
Iv Operator Averaging
Once a trial state has been prepared, the next crucial step in the VQE is the evaluation of the objective function corresponding to the problem operator , . One possibility is to use the quantum phase estimation algorithm Abrams and Lloyd (1997, 1999); AspuruGuzik et al. (2005). If is an eigenstate, then the value is obtained after a single state preparation with a cost in the desired precision of . Unfortunately, to achieve this precision, all of the operations must be coherent which is a prohibitive technological requirement for current and nearterm quantum computers. Moreover, if the state is instead a mixture of many eigenstates, it will still require repetitions of the entire procedure to converge the value to a precision . The use of quantum phase estimation done to a precision surpassing opens the possibility to instead minimize the minimal value found in a projective measurement of the energy in a sequence of phase estimation runs. However we do not explore that option further here.
In 2014, Peruzzo and McClean et. al Peruzzo et al. (2014) suggested a way to retain the advantage of preparing classically inaccessible states while removing the overwhelming coherence time requirements to measure the energy. This method is called Hamiltonian averaging and has been discussed recently in more detail McClean et al. (2014).
The original formulation used the fact that tensor products of Pauli operators form a basis for the space of Hermitian operators. As such any Hermitian operator may be written as
(56) 
and by linearity the expectation value as
(57) 
As a result, all that is required is the weighted sum of the results from simple Pauli measurements. This is an operation requiring coherence time assuming parallel qubit rotation and readout are possible, otherwise the coherence time required is where is the locality of the term to be measured. Previously, some scaling analysis of this procedure was done in the context of locality McClean et al. (2014), but here we detail more specifically how to perform the averaging and verify the error on the fly in a simulation of a general state.
Consider the Hamiltonian decomposed as
(58) 
where each is a Hermitian operator with associated measurement outcomes and , of which Pauli operators are a special case. In order to get the desired precision in a normal distribution approximation, we require a variance of in the estimator of , which we denote with a large hat as . The estimator we have described is constructed as a sum of independent estimators ,
(59) 
each of which is a built a sequence of independent measurements . As the measurements are taken from independent state preparations, we have that the covariance between the individual estimators on the measurements is or and thus the variance of the total estimator is the sum of the variances of the individual estimators
(60) 
The individual estimators are constructed as the mean of a sequence of independent measurements corresponding to the operator on independent preparations of the state . Each measurement of the total operator requires a state preparation and measurement for each individual term, and thus the total number of expected state preparations and measurements to achieve a precision of in is
(61) 
While this offers insight into how many measurements one expects to take, it does not yet constitute a practical algorithm, as the true value of the variances in general will be unknown except in toy examples. Instead, one has access to the sample mean and unbiased sample variance as the measurements are taken. That is, after measurements of the operator have been taken on , one computes
(62) 
and continues taking measurements until , and moves on to the next term. While straightforward, this methodology suffers from some ambiguities when using a small number of measurements or when the state represents an eigenstate of the operator . In particular, how many measurements are required to confirm that the variance is to the desired precision. This is related to how unobserved events are addressed in a frequentist perspective of probability. In practical implementations these issues are often left unaddressed rigorously in stochastic sampling methods and a reasonable minimum number of measurements is chosen such as or before the estimates of are taken to be reliable, trusting that after a number of samples that it is well represented by a normal distribution and the higher moments associated with errors in estimates of the variance vanish rapidly. An alternative perspective that addresses such concerns from the outset is a Bayesian perspective, which has been investigated in the context of quantum phase estimation Wiebe and Granade (2015), and we now explore in the context of Hamiltonian averaging.
iv.1 Bayesian Perspective
In a Bayesian perspective, we start from an uninformative prior for the distribution . In the case of two measurement outcomes, the likelihood function is the binomial likelihood, and the posterior distributions after measurement can be worked out analytically when used with a conjugate Beta prior. These distributions are welldefined even for small numbers of measurements or when is close to an eigenstate of , resulting in potentially unobserved events in a sequence of measurements.
Consider a sequence of independent measurements with two possible outcomes , such as the quantum measurement of a Pauli operator. The likelihood of observing the sequence of measurements is completely defined by a single variable , and is written
(63) 
with being the total number of measurements and being the number of measurements equal to . The value defines the probability of observing and will be directly related to . Our current knowledge of is defined by the prior distribution . Many choices for the form of the prior distribution can be made, but an analytical result can be obtained by choosing the conjugate prior to the Binomial distribution, which is the Beta distribution
(64) 
The Beta distribution is a function of two parameters and , and these are the parameters we will seek to update with a Bayes inference scheme. Simply put, given the measurements with instances of , the posterior distribution is given by
(65) 
From and , one can determine both the mean value and variance in our desired quantity as
(66)  
(67) 
and the expected value and variance of may be used in the estimators associated with . In particular
(68)  
(69) 
A reasonable choice of initial prior in this situation before any measurements are taken is the uniform prior (sometimes called the Bayes’ prior probability in this case) . Thus a practical strategy in the Bayes setting is to let , then take measurements. One then updates and to and according to eq. 65, and continues taking measurements until , which is simply computed as a function of the new and through the above formulae. We note that if one has a good reference state, a prior distribution can be constructed from it to yield an informative prior. This has the potential to reduce the cost and will converge to the same result under most reasonable conditions. However one must be careful as this may introduce a bias for poor reference states with a small number of measurements
After using either the frequentist or Bayesian approach to check convergence of for all , under a normal distribution approximation the final estimation of is precise to the desired precision .
An alternative to the normal approximation confidence intervals may be used in the Bayesian approach if desired. As the measurements are taken for each of the operators in the Bayesian approach, the associated probability distribution is known. The probability distribution of a sum of independent random variables is known to be the convolution of the individual probability distributions, such that