Adaptive Variational Quantum Imaginary Time Evolution Approach for Quantum Chemistry Calculations
Niladri Gomes, Anirban Mukherjee, Feng Zhang, Thomas Iadecola, Cai-Zhuang Wang, Kai-Ming Ho, Peter P. Orth, Yong-Xin Yao
AAdaptive Variational Quantum Imaginary Time Evolution Ap-proach for Quantum Chemistry Calculations
Niladri Gomes Anirban Mukherjee Feng Zhang Thomas Iadecola Cai-Zhuang Wang Kai-MingHo Peter P. Orth Yong-Xin Yao ∗ N. Gomes, A. Mukherjee, F. ZhangAmes Laboratory, Ames, Iowa 50011, USAT. Iadecola, C.-Z. Wang, K.-M. Ho, P.P. Orth, Y.-X. YaoAmes Laboratory, Ames, Iowa 50011, USADepartment of Physics and Astronomy, Iowa State University, Ames, Iowa 50011, USA ∗ [email protected]: quantum computation, quantum algorithms, quantum chemistry An adaptive variational quantum imaginary time evolution (AVQITE) approach is introduced that yields efficient representations ofground states for interacting Hamiltonians on near-term quantum computers. It is based on McLachlan’s variational principle ap-plied to imaginary time evolution of variational wave functions, and avoids difficult nonconvex high-dimensional optimization thatplagues other variational approaches such as the variational quantum eigensolver. The variational parameters evolve according toequations of motions that minimize the difference to the exact imaginary time evolution, which is quantified by the McLachlan dis-tance. Rather than working with a fixed variational ansatz, where the minimal McLachlan distance is constrained by the quality ofthe ansatz, the adaptive method iteratively expands the ansatz along the dynamical path to ensure the McLachlan distance remainsbelow a chosen threshold. AVQITE is used to prepare ground states of H , H O and BeH molecules and yields compact variationalans¨atze and ground state energies beyond chemical accuracy. Finally, quantum Lanczos calculations can also be naturally performedalongside AVQITE without additional quantum resource costs. Quantum computers promise to solve a subset of classically NP-hard problems in polynomial time ata bounded error probability, with quantum simulation as an important example [1]. In the long term,given access to fault-tolerant quantum computers, adiabatic state preparation followed by quantum phaseestimation may become the standard algorithm to determine the ground state energy of a quantum chem-istry Hamiltonian [2, 3, 4]. The required circuit depth, however, is beyond capabilities of near-term noisyintermediate-scale quantum (NISQ) devices, making low-depth hybrid quantum-classical algorithm suchas the variational quantum eigensolver (VQE) much more promising to achieve quantum advantage [5, 6,3, 4]. VQE takes the expectation value of a Hamiltonian, which is measured on a quantum device, as acost function with a set of variational parameters that are optimized using classical algorithms. Excited-state calculations using VQE have also been proposed by modifying the energy cost function or low en-ergy subspace expansion through linear response [5, 7, 8, 9].Meanwhile, quantum imaginary time evolution (QITE) has been developed as an alternative approachto prepare ground states on quantum computers. QITE inherits the advantage of classical imaginarytime evolution algorithms, which allow correlations to build faster than would be allowed by the Lieb-Robinson bound that governs real time evolution [10]. In QITE, the difficult high-dimensional noncon-vex optimization of the energy cost function that arises in VQE is replaced by evolving the variationalparameters along the imaginary time axis. The development of the QITE method includes two direc-tions: the unitarized QITE approach represents the imaginary time propagator by a unitary one via least-square fitting [11], while the variational QITE (VQITE) approach is based on variational quantum sim-ulation with fixed ans¨atze [12, 13]. The accuracy of the unitarized QITE method is tied to the size ofcorrelation domains associated with the Hamiltonian. Specifically, the circuit depth grows exponentiallywith the domain sizes (being roughly the system’s correlation length) and linearly with the number ofimaginary time steps. For practical implementations, strategies to reduce circuit complexity, e.g. by uti-lizing symmetries and effectively combining unitaries [14, 15, 16, 17], have been proposed. The VQITE a r X i v : . [ phy s i c s . c h e m - ph ] F e b ethod has the advantage of fixed circuit depth along the imaginary-time path, but its accuracy is lim-ited by the fidelity of the variational ansatz in representing the ground state. While various strategies toconstruct variational ans¨atze have been reported during the development of VQE [18, 6, 19, 20, 21, 22],their accuracy can be system-dependent and the variational circuits can often be suboptimal [23, 24].One promising strategy is to perform VQE with an adaptively generated ansatz, where operators aredrawn from a predefined operator pool and iteratively added to the ansatz during the calculation. It wasshown that adaptive VQE can yield highly accurate and compact variational ans¨atze for specific prob-lems [25, 26, 27, 28, 29].In this work, we develop an adaptive VQITE (AVQITE) method for variational quantum imaginary timeevolution to efficiently prepare ground states of interacting fermion systems with high accuracy. Themethod generalizes the recently proposed adaptive variational quantum dynamics simulation (AVQDS)method [30] from real to imaginary time. Like AVQDS, the variational ansatz in AVQITE is automati-cally generated by choosing optimal multi-qubit Pauli rotation gates along the dynamical path to keepa measure of ansatz quality, the McLachlan distance, within a desired accuracy. We demonstrate thecapabilities of AVQITE calculations by preparing the ground states of an H chain and H O and BeH molecules at representative bond lengths with increasing electron correlation effects. We find total ener-gies to chemical accuracy with compact ans¨atze similar to qubit-ADAPT-VQE results, yet without re-sorting to the complex optimization of parameters in a high-dimensional nonconvex energy landscape.Furthermore, Quantum Lanczos calculations can be carried out together with AVQITE, leading to fasterconvergence to the ground state. We envision AVQITE, with its compact variational circuits and avoid-ance of explicit high-dimensional optimization, as a viable way to efficiently prepare ground states of in-teracting fermion systems (e.g., molecules) on NISQ devices. The AVQITE algorithm generalizes the recently introduced AVQDS approach [30] from real to imagi-nary time evolution. By time evolving a quantum state in imaginary time, it yields the ground state ofa system as the final state. While the derivation of AVQITE resembles that of AVQDS, there are a fewkey differences that we point out in the following.
The theory of VQITE has been developed in reference [12] and has been used to find ground-state en-ergies of H and LiH molecules on classical simulators [12, 13]. Here, we review the VQITE formalismwithin the density matrix approach, which is insensitive to the global phase of the quantum state. Con-sider a system with Hamiltonian ˆ H in a pure state | Ψ (cid:105) . The evolution of the density matrix ˆ ρ ≡ | Ψ (cid:105)(cid:104) Ψ | under the imaginary-time propagator e − τ ˆ H is governed by the Liouville–von Neumann-type equation [12,31] d ˆ ρdτ = L [ ˆ ρ ] , (1)where the superoperator L [ ˆ ρ ] = −{ ˆ H , ˆ ρ } + 2 (cid:104) ˆ H(cid:105) ˆ ρ with the anticommutator { ˆ H , ˆ ρ } = ˆ H ˆ ρ + ˆ ρ ˆ H , andthe Hamiltonian expectation value (cid:104) ˆ H(cid:105) = Tr (cid:104) ˆ ρ ˆ H (cid:105) . The so-called imaginary time τ ∈ R is a real posi-tive parameter. For a generic variational ansatz | Ψ[ θ ] (cid:105) with a real parameter vector θ of dimension N θ ,the squared McLachlan distance L is defined by the Frobenius norm of the difference between the varia- .1 Variational Quantum Imaginary Time Evolution method tional and exact state propagations along the imaginary time axis, i.e. [12]: L ≡ (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) N θ (cid:88) µ =1 ∂ ˆ ρ [ θ ] ∂θ µ ˙ θ µ − L [ ˆ ρ ] (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) = (cid:88) µν M µν ˙ θ µ ˙ θ ν − (cid:88) µ V µ ˙ θ µ + Tr (cid:2) L [ ˆ ρ ] (cid:3) , (2)which is a quadratic function of the time derivatives { ˙ θ µ ≡ ∂θ µ /∂τ } . The real symmetric N θ × N θ matrix M is specified as M µν = 2 Re (cid:20) ∂ (cid:104) Ψ[ θ ] | ∂θ µ ∂ | Ψ[ θ ] (cid:105) ∂θ ν + ∂ (cid:104) Ψ[ θ ] | ∂θ µ | Ψ[ θ ] (cid:105) ∂ (cid:104) Ψ[ θ ] | ∂θ ν | Ψ[ θ ] (cid:105) (cid:21) . (3)The real vector V of dimension N θ is defined as V µ = 2 Re (cid:20) − ∂ (cid:104) Ψ[ θ ] | ∂θ µ ˆ H | Ψ[ θ ] (cid:105) + (cid:104) Ψ[ θ ] | ∂ | Ψ[ θ ] (cid:105) ∂θ µ (cid:104) ˆ H(cid:105) θ (cid:21) = 2 Re (cid:20) − ∂ (cid:104) Ψ[ θ ] | ∂θ µ ˆ H | Ψ[ θ ] (cid:105) (cid:21) (4)with the abbreviation (cid:104) ˆ H(cid:105) θ ≡ (cid:104) Ψ[ θ ] | ˆ H | Ψ[ θ ] (cid:105) . The second term in the first line of equation (4) vanishesdue to the normalization (cid:104) Ψ[ θ ] | Ψ[ θ ] (cid:105) = 1. The final term in Eq. (2) can be expressed in terms of theenergy variance as Tr (cid:2) L [ ˆ ρ ] (cid:3) = 2 (cid:16) (cid:104) ˆ H (cid:105) θ − (cid:104) ˆ H(cid:105) θ (cid:17) = 2 var θ [ ˆ H ] . (5)The McLachlan variational principle amounts to minimizing the quadratic cost function L with respectto { ˙ θ µ } . This leads to the following linear equations of motion: (cid:88) ν M µν ˙ θ ν = V µ . (6)The second term in equation (3) originates from the global phase of the wavefunction [13, 12]. Below weadopt a purely real wave function ansatz in a pseudo-Trotter form (with odd number of Pauli-Y terms)for which the global phase contribution vanishes. Consequently, the density-matrix and wavefunction-based derivations reach exactly the same results.The optimal McLachlan distance L of the variational ansatz Ψ[ θ ] given by L = 2 var θ [ ˆ H ] − (cid:88) µν V µ M − µν V ν . (7)VQITE is formulated in the exactly same form as real time variational quantum dynamics simulations(VQDS) [12, 30], with exception of the definition of vector V in Equation (4) due to the different super-operator L [ ˆ ρ ] in the Liouville–von Neumann-type equation. A typical VQITE calculation, which integrates the equation of motion (6) with a constant time step ∆ τ according to the Euler method, is illustrated in the green charts of Figure 1. Given a fixed variationalansatz initialized to a reference state | Ψ (cid:105) that is easily prepared on a quantum computer, the energy ex-pectation value is first measured. The convergence criterion, such as the energy difference between twoconsecutive steps, is checked. If the convergence condition is not satisfied, the real symmetric matrix M in Equation (3) and vector V in Equation (4) are determined, with which the step size ∆ θ of the vari-ational parameter vector θ at time step ∆ τ is calculated. The ansatz state | Ψ[ θ ] (cid:105) with the updated pa-rameters triggers another iteration until convergence is reached. Note that the energy variance var θ [ ˆ H ]provides a quality measure of the variational ansatz in approximating the ground state. .2 Adaptive Variational Quantum Imaginary Time Evolution Method ൿ ۧ |Ψ[𝜽] = |Ψ 𝐿 = 𝐿 𝜈2 = min 𝐿 , 𝐿 , ⋯ , 𝐿 𝑁 𝑝 ۧ|Ψ[𝜽] → 𝑒 −𝑖𝜃 ′ መ𝒜 𝜈 ۧ|Ψ[𝜽] | 𝜃 ′ =0 𝜽 → 𝜽 + ∆𝜽 ∆𝜽 = 𝑀 −1 𝑉 ∙ ∆𝜏 evaluate
𝑀, 𝑉
𝐸 = Ψ[𝜽] ℋ Ψ[𝜽] converged?
Yes
No Exitadaptively expand ansatz
Ψ[𝜽] evaluate 𝐿 𝐿 < 𝐿 cut2 Yes መ𝒜 መ𝒜 𝑒 −𝑖𝜃 ′ መ𝒜 ۧ|Ψ 𝜽 | 𝜃 ′ =0 መ𝒜 𝑁 𝑝 ⋯ operator pool NoEvaluate 𝐿 Evaluate 𝐿 Evaluate 𝐿 𝑁 𝑝 𝑒 −𝑖𝜃 ′ መ𝒜 ۧ|Ψ 𝜽 | 𝜃 ′ =0 𝑒 −𝑖𝜃 ′ መ𝒜 𝑁𝑝 ۧ|Ψ 𝜽 | 𝜃 ′ =0 Figure 1:
Schematic illustration of variational quantum imaginary time evolution algorithm, with an addi-tional module to adaptively expand the ansatz.
The green flowchart on the left shows a typical VQITE calculation.In AVQITE, a module (blue) is introduced to adaptively expand the variational ansatz by selectively appending parametricrotation gates to keep the McLachlan distance L under a threshold L along the imaginary-time evolution path. The VQITE formalism is presented for a generic wave function ansatz. In the development of VQE, vari-ational ans¨atze | Ψ[ θ ] (cid:105) = ˆ U [ θ ] | Ψ (cid:105) with two different forms of the unitary operator U [ θ ] acting on a ref-erence state | Ψ (cid:105) have been proposed. In hardware efficient ans¨atze, the unitary operator ˆ U [ θ ] is a prod-uct of parametrized native gates of the real device, e.g., single-qubit rotation gates plus two-qubit entan-gling gates [18]. Alternatively, at a higher algorithmic level, ˆ U [ θ ] can be expressed as a product of N θ multi-qubit rotation gates in a pseudo-Trotter form: | Ψ[ θ ] (cid:105) = N θ − (cid:89) µ =0 e − iθ µ ˆ A µ | Ψ (cid:105) , (8)where ˆ A µ are Hermitian operators. The unitary coupled cluster ansatz and its variants [6, 19, 20], theHamiltonian variational ansatz [21], and the ansatz in the quantum approximate optimization algorithm(QAOA) all belong to this category [22]. In AVQITE, we adopt a variational ansatz in the above form (8),and allow the number of ˆ A -operators to be dynamically expanded along the imaginary-time evolutionpath to maintain high accuracy in representing the evolving quantum state, as shown in the blue moduleof Figure 1.The initial steps of an AVQITE calculation are the same as those of a VQITE calculation, with the ex-ception of great flexibility in the choice of initial ansatz | Ψ[ θ ] (cid:105) . While the simulation accuracy is tied tothe initial ansatz in VQITE calculations, an AVQITE calculation monitors the quality of | Ψ[ θ ] (cid:105) in rep-resenting the evolving quantum state and adaptively expands the form of | Ψ[ θ ] (cid:105) to maintain a fixed ac-curacy. In fact, AVQITE calculations can simply take any | Ψ (cid:105) as the initial ansatz. For convenience,we use a product state as our initial state in all our calculations below, which is easily prepared on aquantum processor unit (QPU). In the first iteration where no variational parameters are present, theMcLachlan distance L in Equation (7) is determined by the energy variance var[ ˆ H ] in state | Ψ (cid:105) , which .2 Adaptive Variational Quantum Imaginary Time Evolution Method is generally larger than the threshold L , since | Ψ (cid:105) is not an eigenstate of ˆ H . As a result, a predefinedoperator pool is scanned and an operator ˆ A ν is chosen to construct a unitary e − iθ (cid:48) ˆ A ν to be appended tothe variational ansatz | Ψ[ θ ] (cid:105) , which produces minimal McLachlan distance, as illustrated in Figure 1.The additional parameter θ (cid:48) associated with the new operator is always initialized to zero to keep theimaginary time evolution of the variational state continuous. Nevertheless, the McLachlan distance canstill be reduced by addition of e − iθ (cid:48) ˆ A ν , because it involves derivatives of the ansatz. The adaptive pro-cedure of selectively appending a new unitary to the ansatz continues until the updated McLachlan dis-tance satisfies L < L . The expanded variational parameter vector θ is subsequently updated at theimaginary time step as in VQITE, and new iterations proceed until energy convergence is reached. The accuracy of AVQITE calculations is controlled by the McLachlan distance threshold L , while theability to reach a compact final ansatz | Ψ[ θ ] (cid:105) is tied to the operator pool. For quantum chemistry cal-culations, different ways to efficiently construct operator pools and pool completeness conditions havebeen extensively discussed in the context of adaptive approaches to VQE [25, 26, 27]. The fermionic op-erator pool proposed in reference [25] is composed of the single excitation operators and double excita-tion operators with respect to a Hartree-Fock (HF) reference state. When translated to a qubit repre-sentation using encoding methods such as Jordan-Wigner (JW) mapping [32], a single excitation op-erator can result in a weighted sum of two Pauli strings, and six Pauli strings for a double excitation.Here a Pauli string is defined as a product of single qubit Pauli operators, namely X, Y and Z . Alter-natively, a qubit operator pool can also be constructed directly with rudimentary Pauli strings [26, 27].Compared with fermionic operator pools, adaptive-VQE calculations with qubit operator pools generatevariational ans¨atze with considerably shallower circuits at the price of more variational parameters. Asthe complex nonconvex optimization problem in the high-dimensional parameter space of VQE is com-pletely avoided in AVQITE, the qubit operator pool is therefore very appealing and adopted in the fol-lowing AVQITE calculations. In the AVQITE calculations, we construct a qubit operator pool by choos-ing all the Pauli strings present in the fermionic single and double excitation operators of the unitarycoupled cluster ansatz [6, 19]. The parity mapping is used to transform the fermionic excitation opera-tors to qubit operators [33, 34], as fermion to qubit mappings other than JW have not been studied be-fore in the context of operator pool construction [26]. Since the fermionic excitation operators are real,every Pauli string in the operator pool contains an odd number of Pauli Y operators and is thereforeantisymmetric. The unitaries in equation (8) are thus real and an initially real wave function (i.e. withreal coefficients) remains real when evolving along the imaginary time path. Consequently, the expres-sion ∂ (cid:104) Ψ[ θ ] | ∂θ µ | Ψ[ θ ] (cid:105) = (cid:104) Ψ µ − | i ˆ A µ | Ψ µ − (cid:105) in Eq. (3), where | Ψ µ [ θ ] (cid:105) = (cid:81) µµ (cid:48) =0 e − iθ µ (cid:48) ˆ A µ (cid:48) | Ψ (cid:105) , vanishes for any µ .Therefore, the second term of M in Eq. (3) which originates from the global phase of the wavefunctionvanishes, and the density-matrix and wavefunction-based approaches lead to exactly the same results.As discussed previously, the VQITE and real time VQDS calculations amount to integrate an equationof motion of the same form (6). While sufficiently small step size is necessary to maintain the high accu-racy of VQDS results along the dynamical path, relatively bigger step size can be adopted in VQITE be-cause only the final state is of interest [14]. In the following AVQITE calculations of molecules in Hartreeatomic units, we find ∆ τ = 0 . The implementation of AVQITE amounts to measuring the symmetric matrix M in Eq. (3), the vec-tor V in Eq. (4) and the scalars (cid:104) ˆ H(cid:105) θ , (cid:104) ˆ H (cid:105) θ on a quantum device, all of which has been discussed inthe context of VQITE in reference [6, 35]. The energy (cid:104) ˆ H(cid:105) θ can be obtained as a weighted sum of ex-pectation values of the Pauli strings in the Hamiltonian, a procedure often termed “Hamiltonian aver-aging” [6]. The expectation value (cid:104) ˆ H (cid:105) θ can be obtained similarly by replacing ˆ H with ˆ H . The circuit mplementation to measure V , which is essentially the energy gradient, has also been discussed in thecontext of VQE optimization [36, 37]. The proposed indirect measurement circuit introduces an ancil-lary qubit with a Hadamard-type test. Two controlled-unitary gates, specifically, controlled multi-qubitPauli gates, need to be implemented. As the implementation of controlled-unitary gates can be chal-lenging for NISQ devices, general strategies to replace indirect measurements by direct measurementshave also been proposed [38]. It is especially appealing to use the parameter-shift rule to evaluate theenergy gradient V [39, 40, 41, 42], as it only requires some additional measurements of Hamiltonian ex-pectation values at shifted parameters, without necessity of introducing new circuits. For the matrix el-ement M µν = 2 Re (cid:104) ∂ (cid:104) Ψ[ θ ] | ∂θ µ ∂ | Ψ[ θ ] (cid:105) ∂θ ν (cid:105) in Eq. (3), where we have the vanishing global phase term with thepseudo-Trotter ansatz (8), the diagonal element can be simplified to M µµ = 2 (cid:104) Ψ µ [ θ ] | ˆ A µ | Ψ µ [ θ ] (cid:105) , where | Ψ µ [ θ ] (cid:105) = (cid:81) µµ (cid:48) =0 e − iθ µ (cid:48) ˆ A µ (cid:48) | Ψ (cid:105) . This can be measured with the Hamiltonian averaging method for theHermitian operator ˆ A µ . Because ˆ A µ is a single Pauli string for the qubit operator pools adopted here,ˆ A µ is an identity operator and M µµ = 2. The off-diagonal elements of M can be simplified to M µν =2 Re (cid:104) − ∂ (cid:104) Ψ ν [ θ ] | ∂θ µ i ˆ A ν | Ψ ν [ θ ] (cid:105) (cid:105) with µ < ν , which bears the same form as V in Eq. (4) with ˆ H → i ˆ A ν andcan also be measured using the parameter-shift rule [39]. The focus of the QITE approach lies on the final quantum state, which converges to the ground state forlarge enough time. Within QITE, all previous quantum states along the path are discarded. In contrast,the quantum Lanczos (QL) method was developed to achieve a more efficient calculations of ground andexcited state energies by exploiting information contained in all quantum states along the path [11]. Theessence of QL is to diagonalize the Hamiltonian within the Krylov subspace spanned by a subset of imaginary-time states {| Ψ n (cid:105) = C n e − n ∆ τ ˆ H | Ψ (cid:105)} with even time step indices n = 0 , , . . . N and normalizationconstants C n . In the classical Lanczos algorithm, the reduced Hamiltonian in the Krylov subspace isbrought to a tridiagonal form by sequentially applying the Hamiltonian on orthonormalized Krylov ba-sis vectors [43]. While for convenience, QL directly represents the reduced Hamiltonian ˆ H in the normal-ized basis {| Ψ n (cid:105)} characterized by an overlap matrix S nn (cid:48) = C n C n (cid:48) C n + n (cid:48) ) / . As the normalization coefficient C ( n + n (cid:48) ) / is used, only the even- n imaginary-time states (for which ( n + n (cid:48) ) / H nn (cid:48) = S nn (cid:48) E ( n + n (cid:48) ) / , where the expectation value E n ≡ H nn = (cid:104) Ψ n | ˆ H | Ψ n (cid:105) . The nor-malization coefficient C n can be calculated recursively as C − n = C − n − (1 − τ E n − + O (∆ τ )) with C = 1. Therefore, only the expectation values of the Hamiltonian { E n } are needed to set up the gener-alized eigenvalue equation in the Krylov subspace for QL. As { E n } are readily available in an AVQITEcalculation, a QL calculation can be efficiently performed alongside AVQITE to get the energy eigenval-ues with no additional quantum resource costs. Nevertheless, as the Lanczos eigenvector is expressed asa linear combination of imaginary time states, it can be much more involved to measure expectation val-ues of other observables beyond energy on a quantum computer [44]. The ab initio nonrelativistic molecular electron Hamiltonian is given byˆ H = (cid:88) pq (cid:88) σ h pq ˆ c † pσ ˆ c qσ + 12 (cid:88) pqrs (cid:88) σσ (cid:48) h pqrs ˆ c † pσ ˆ c † rσ (cid:48) ˆ c sσ (cid:48) ˆ c qσ , (9)where the one-electron core part h pq = (cid:82) d r φ ∗ p ( r )( T + V ion ) φ q ( r ), and the two-electron Coulomb integralis obtained as h pqrs = (cid:90) (cid:90) d r d r (cid:48) φ ∗ p ( r ) φ ∗ r ( r (cid:48) ) V ee ( | r − r (cid:48) | ) φ s ( r (cid:48) ) φ q ( r ) . (10) E ( H a ) (a) H AVQITEQLanczosExact 0 2 40510152025 N (b) 0 2 410 (c) L cut L var [ ]0 10 10 e rr o r ( k c a l / m o l ) Figure 2:
AVQITE calculation for H chain. Along the imaginary time path of τ = N ∆ τ , the total energies of H from AVQITE and quantum Lanczos calculations are shown in panel (a), the number of variational parameters N θ in (b),and the McLachlan distance L and energy variance var θ [ ˆ H ] in (c). The exact full configuration interaction result is shownas the dashed line in (a) for reference. The total energy errors E − E Exact are plotted in the inset of (b), with the shadedarea denoting chemical accuracy. The chosen threshold L = 5 × − is indicated by a dashed line in (c). Here p, q, r, s are composite indices for atom and orbital, and σ is the spin index. T , V ion and V ee are thekinetic energy, ionic potential operator and Coulomb interaction operator, respectively. The standardSTO-3G minimal basis set is adopted for the basis orbital functions { φ ( r ) } . In the following AVQITEcalculations, the PySCF quantum chemistry package is first used to generate the molecular Hamilto-nian (9) and produce the restricted Hartree-Fock(HF) solution [45]. A basis transformation from atomicorbitals to molecular orbitals is performed to the Hamiltonian (9) for the convenience of preparation ofthe initial HF state as a tensor product state on the quantum computer. The parity transformation isapplied to get the qubit representation of the molecular Hamiltonian, where two qubits are tapered dueto the Z symmetry from total electron number and total spin conservation [33, 46, 47].In figure 2 we illustrate a detailed numerical AVQITE calculation of an H chain molecule at a uniformbond length R = 1 . | Ψ (cid:105) , the total energy monotonically decreases and converges toward the exactresult with increasing imaginary time τ = N ∆ τ at ∆ τ = 0 .
1, as shown by the blue circles in figure 2(a).The associated error, which is defined as the difference between the AVQITE energy and the exact resultfrom a full configuration interaction (FCI) calculation E − E Exact , is shown in the inset of figure 2(b) on alog scale. The AVQITE calculation reaches chemical accuracy of 1 kcal/mol after τ = 4 .
7. The quantumstate fidelity [1], which is defined as the squared overlap between the ansatz state and the exact groundstate |(cid:104) Ψ[ θ ] | Ψ Exact (cid:105)| , goes beyond 99 . N θ , or equivalently thenumber of multi-qubit Pauli rotation gates, increases at several imaginary time steps and flattens at 24.For comparison, the qubit-ADAPT VQE calculation of H using the same operator pool generates a fi-nal ansatz of 30 variational parameters upon reaching chemical accuracy. Specific information about theMcLachlan distance L is plotted in figure 2(c), where L is reduced below the threshold L = 5 × − by adaptively expanding the variational ansatz whenever the initial value L > L at any time step.Together with L , the energy variance var θ [ ˆ H ] along the imaginary path is shown to be almost linear ona semi-log scale, implying an exponential convergence. In fact, the AVQITE total energies E ( τ ) can alsobe fitted with an exponential function, E ( τ ) = E ∞ + e − aτ , with E ∞ = − . E (var θ [ ˆ H ]) = E + b var θ [ ˆ H ] due to the approximately linear relation between E and var θ [ ˆ H ]in the numerical results, where E = − . R HH ( Å )2.01.51.0 E ( H a ) H (a) ExactAVQITE 1.0 1.5 2.0 R HO ( Å )75.0074.7574.5074.25 H O(b) 1 2 3 R BeH ( Å )15.515.014.5 BeH (c)0 100 20010 E rr o r ( k c a l / m o l ) (d)
25 50ADAPT-VQE2550 AV Q I T E N (e)
25 50ADAPT-VQE2550 AV Q I T E (f)
25 50ADAPT-VQE2550 AV Q I T E nu m b e r o f C N O T s (g) 0 2 40200400 (h) 0 2 40200400 (i) Figure 3:
AVQITE calculations of molecules with varying electron correlations.
Upper panels: AVQITE totalenergies of H , H O and BeH molecules at three bond lengths colored by red, green and blue. The exact potential energycurves from FCI calculations are also shown for reference. Middle panels: AVQITE energy error, E − E Exact , as a functionof imaginary time step τ = N ∆ τ for the molecules at bond lengths with the same color coding. The shaded area indicateserrors smaller than chemical accuracy. Insets: Number of variational parameters N θ of the AVQITE ansatz vs. that forqubit-ADAPT-VQE when the energy accuracy reaches 0.1 mHa. The symbol “ × ” indicates qubit-ADAPT-VQE does notconverge. Lower panels: Dependence of the number of CNOT gates in the AVQITE state preparation circuit as a functionof imaginary time for each set of molecules. error 0.1 kcal/mol.The quantum Lanczos results, which are conveniently calculated along with AVQITE at no additionalquantum resource cost, are plotted as orange squares in figure 2(a) for the total energy, with the errorshown in the inset of panel (b). The QL calculation converges relatively faster than AVQITE, since thereduced Lanczos Hamiltonian encodes information beyond the latest quantum state. Because the QLmethod is generally susceptible to numerical inaccuracies and the way to evaluate the Hamiltonian andoverlap matrix in the Krylov subspace used in the QL method accumulates sizable errors, excited statescannot be directly accessed in the current QL calculations.To demonstrate the general applicability of the AVQITE approach in generating compact ground stateans¨atze, we perform AVQITE calculations for H chains at bond length R HH = 0 . . . O at R OH = 0 . . . at R BeH = 0 . . . at three bond lengths colored in red, green and blue, which agree with the FCI results be-yond chemical accuracy. The detailed error convergence, E − E Exact , as a function of imaginary time step = N ∆ τ is shown in panel (d). With increasing bond length or correlation energy, the critical imagi-nary time τ c , which is defined as the time where chemical accuracy is reached, generally increases. Herethe correlation energy is defined as the difference between the HF and FCI energies, which correspondsto the initial AVQITE energy error as the AVQITE ansatz starts with the HF state. At R HH = 2 . τ c ≈ τ formolecules at larger bond length to reduce the number of imaginary time steps N . In this case, ∆ τ is in-creased from 0.1 to 0.5 to speed up the convergence, although sizable energy fluctuations are present inthe initial time steps. The number of controlled-NOT gates (CNOTs) generally increases in the initialsteps and levels off for τ above 3, as shown in panel (g). The maximal number of CNOTs reaches 250 forH at R HH = 2 . e − iθ ˆ σ with Pauli string ˆ σ of length p needs 2( p −
1) CNOTs (assuming all-to-all connectivity) [1].Moving on to the AVQITE calculations of H O and BeH at representative bond lengths, the error con-vergence behavior generally remains similar to the H calculations, with final energies beyond chemi-cal accuracy. Likewise, the number of CNOTs in the AVQITE state preparation circuit generally in-creases initially and levels off for τ >
3. The total number of CNOTs remains close to or under 400.Remarkably, the positive correlation between the number of CNOTs in the AVQITE ansatz and cor-relation energy, which seems to exist in the calculations of H , does not apply to H O and BeH . TheAVQITE state preparation circuits in the test-set are generally over one order of magnitude shallowerthan the original UCCSD ansatz. In the insets of the middle panels, we show that the number of varia-tional parameters N θ of the AVQITE ansatz is generally quite close to that of the qubit-ADAPT-VQEapproach at the energy accuracy of 0.1 mHa, albeit with appreciable difference in the case of BeH . TheH chain at R HH = 2 . ab initio Hamiltonian in Eq. (9) scales as N for molecules with N orb basis orbitals due to the presence of two-body Coulomb interactions. Theresource cost for straightforward evaluation of (cid:104) ˆ H (cid:105) θ then scales as N . The operator pool size growsas N if the Pauli operators are obtained from the single and double excitations of the UCCSD ansatz.Therefore, the resource cost for an AVQITE calculation, limited by measuring V in Eq. (4) and (cid:104) ˆ H (cid:105) θ ,scales as N , because all the measurements can be achieved using the Hamiltonian averaging methodwith the help of the parameter-shift rule as discussed previously in Sec. 2.2.3. Multiple techniques fromthe literature can be utilized to reduce this high-order polynomial scaling. The low-rank tensor factor-ization of the Hamiltonian coefficients h pq and h pqrs reduces the number of distinct measurement cir-cuits for (cid:104) ˆ H(cid:105) θ from N to linear scaling, with small overhead of basis transformation through Givensrotation [51, 52, 53]. Likewise, the resource cost for measuring (cid:104) ˆ H(cid:105) θ can be reduced to quadratic scal-ing. Scaling reduction based on reduced density matrices is another interesting approach [54]. Variousapproaches to generate efficient operator pools, including minimal pools whose size scale linearly withthe qubit-basis dimension, have been discussed [23, 26, 27]. Finally, AVQITE can be applied to reduced low energy” models in quantum chemistry such as those that arise within the complete active spaceself-consistent field method (CASSCF) [55]. More generally, it can be used as impurity solver for quan-tum embedding approaches like the rotationally invariant Gutzwiller embedding and density-matrix em-bedding methods [56, 57, 58, 59, 60, 61, 62], extending its application to larger molecules and solid statematerials. The adaptive variational quantum imaginary time evolution (AVQITE) approach is developed by gener-alizing the recently proposed adaptive quantum dynamics simulation (AVQDS) method from real timeto imaginary time. We present benchmark quantum chemistry calculations on H , H O and BeH withdifferent chemical bonding character and atomic spin multiplicity upon dissociation limit. They demon-strate the general applicability of AVQITE to finding accurate and compact variational ans¨atze for inter-acting many-electron models. The key advantage of AVQITE over adaptive VQE approaches is that itbypasses the complicated nonconvex optimization problem in high-dimensional parameter space by per-forming energy cost function minimization along a single imaginary time axis. With the favorable poly-nomial scaling of the AVQITE calculations which automatically generates compact variational ans¨atzeof high accuracy, we envision feasible follow-up applications of this method on NISQ devices to calculatemolecular ground-state properties, with excited states also accessible if combined with techniques likethe folded spectrum method [5, 6, 63, 64]. Acknowledgements
This work was supported by the U.S. Department of Energy (DOE), Office of Science, Basic Energy Sci-ences, Materials Science and Engineering Division. The research was performed at the Ames Laboratory,which is operated for the U.S. DOE by Iowa State University under Contract No. DE-AC02-07CH11358.
References [1] M. A. Nielsen, I. Chuang,
Quantum computation and quantum information , American Associationof Physics Teachers, .[2] A. Aspuru-Guzik, A. D. Dutoi, P. J. Love, M. Head-Gordon,
Science , , 5741 1704.[3] Y. Cao, J. Romero, J. P. Olson, M. Degroote, P. D. Johnson, M. Kieferov´a, I. D. Kivlichan,T. Menke, B. Peropadre, N. P. Sawaya, et al., Chem. Rev. , , 19 10856.[4] S. McArdle, S. Endo, A. Aspuru-Guzik, S. C. Benjamin, X. Yuan, Rev. Mod. Phys. , Nat. Commun. , New J. Phys. , , 2 023023.[7] J. R. McClean, M. E. Kimchi-Schwartz, J. Carter, W. A. De Jong, Phys. Rev. A , , 4042308.[8] R. Santagati, J. Wang, A. A. Gentile, S. Paesani, N. Wiebe, J. R. McClean, S. Morley-Short, P. J.Shadbolt, D. Bonneau, J. W. Silverstone, et al., Sci. Adv. , , 1 eaap9646.[9] O. Higgott, D. Wang, S. Brierley, Quantum , Phys. Rev. B , , 9 094434.[11] M. Motta, C. Sun, A. T. Tan, M. J. O’Rourke, E. Ye, A. J. Minnich, F. G. Brand˜ao, G. K.-L.Chan, Nat. Phys. , , 2 205. EFERENCES [12] X. Yuan, S. Endo, Q. Zhao, Y. Li, S. C. Benjamin,
Quantum , npj Quantum Inf. , , 1 1.[14] N. Gomes, F. Zhang, N. F. Berthusen, C.-Z. Wang, K.-M. Ho, P. P. Orth, Y.-X. Yao, J. Chem.Theory Comput. , , 10 6256.[15] K. Yeter-Aydeniz, R. C. Pooser, G. Siopsis, npj Quantum Inf. , , 1 1.[16] H. Nishi, T. Kosugi, Y.-i. Matsushita, arXiv:2005.12715 .[17] S.-N. Sun, M. Motta, R. N. Tazhigulov, A. T. Tan, G. K. Chan, A. J. Minnich, arXiv:2009.03542 .[18] A. Kandala, A. Mezzacapo, K. Temme, M. Takita, M. Brink, J. M. Chow, J. M. Gambetta, Nature , , 7671 242.[19] P. K. Barkoutsos, J. F. Gonthier, I. Sokolov, N. Moll, G. Salis, A. Fuhrer, M. Ganzhorn, D. J. Eg-ger, M. Troyer, A. Mezzacapo, et al., Phys. Rev. A , , 2 022322.[20] J. Lee, W. J. Huggins, M. Head-Gordon, K. B. Whaley, J. Chem. Theory Comput. , , 1 311.[21] D. Wecker, M. B. Hastings, M. Troyer, Phys. Rev. A , , 4 042303.[22] E. Farhi, J. Goldstone, S. Gutmann, arXiv:1411.4028 .[23] H. R. Grimsley, S. E. Economou, E. Barnes, N. J. Mayhall, In Nat. Commun. [25], 1–9.[24] I. G. Ryabinkin, T.-C. Yen, S. N. Genin, A. F. Izmaylov,
J. Chem. Theory Comput. , , 126317.[25] H. R. Grimsley, S. E. Economou, E. Barnes, N. J. Mayhall, Nat. Commun. , , 1 1.[26] H. L. Tang, E. Barnes, H. R. Grimsley, N. J. Mayhall, S. E. Economou, arXiv:1911.10205 .[27] Y. S. Yordanov, V. Armaos, C. H. Barnes, D. R. Arvidsson-Shukur, arXiv:2011.10540 .[28] A. G. Rattew, S. Hu, M. Pistoia, R. Chen, S. Wood, arXiv , arXiv–1910.[29] M. Ostaszewski, E. Grant, M. Benedetti, arXiv preprint arXiv:1905.09692 .[30] Y.-X. Yao, N. Gomes, F. Zhang, T. Iadecola, C.-Z. Wang, K.-M. Ho, P. P. Orth, arXiv:2011.00622 .[31] M. Berman, R. Kosloff, Comput. Phys. Commun. , , 1-3 1.[32] P. Jordan, E. P. Wigner, In The Collected Works of Eugene Paul Wigner , 109–129. Springer, .[33] S. B. Bravyi, A. Y. Kitaev,
Ann. Phys. , , 1 210.[34] J. T. Seeley, M. J. Richard, P. J. Love, J. Chem. Phys. , , 22 224109.[35] Y. Li, S. C. Benjamin, Phys. Rev. X , , 2 021050.[36] J. Romero, R. Babbush, J. R. McClean, C. Hempel, P. J. Love, A. Aspuru-Guzik, Quantum Sci.Technol. , , 1 014008.[37] G. Ortiz, J. E. Gubernatis, E. Knill, R. Laflamme, Phys. Rev. A , Phys. Rev. Research , , 1 013006.[39] A. Mari, T. R. Bromley, N. Killoran, arXiv:2008.06517 .[40] J. Li, X. Yang, X. Peng, C.-P. Sun, Phys. Rev. Lett. , , 15 150503. EFERENCES [41] K. Mitarai, M. Negoro, M. Kitagawa, K. Fujii,
Phys. Rev. A , , 3 032309.[42] M. Schuld, V. Bergholm, C. Gogolin, J. Izaac, N. Killoran, Phys. Rev. A , , 3 032331.[43] J. W. Demmel, Applied numerical linear algebra , SIAM, .[44] A. M. Childs, N. Wiebe, arXiv:1202.5822 .[45] Q. Sun, T. C. Berkelbach, N. S. Blunt, G. H. Booth, S. Guo, Z. Li, J. Liu, J. D. McClain, E. R.Sayfutyarova, S. Sharma, et al.,
Wiley Interdiscip. Rev.: Comput. Mol. Sci. , , 1 e1340.[46] A. Tranter, S. Sofia, J. Seeley, M. Kaicher, J. McClean, R. Babbush, P. V. Coveney, F. Mintert,F. Wilhelm, P. J. Love, Int. J. Quantum Chem. , , 19 1431.[47] S. Bravyi, J. M. Gambetta, A. Mezzacapo, K. Temme, arXiv:1701.08213 .[48] J. R. McClean, S. Boixo, V. N. Smelyanskiy, R. Babbush, H. Neven, Nat. Commun. , , 1 1.[49] E. Grant, L. Wossnig, M. Ostaszewski, M. Benedetti, Quantum , arXiv:2012.12658 .[51] B. Peng, K. Kowalski, J. Chem. Theory Comput. , , 9 4179.[52] M. Motta, E. Ye, J. R. McClean, Z. Li, A. J. Minnich, R. Babbush, G. K. Chan, arXiv:1808.02625 .[53] W. J. Huggins, J. McClean, N. Rubin, Z. Jiang, N. Wiebe, K. B. Whaley, R. Babbush, arXiv:1907.13117 .[54] J. Liu, Z. Li, J. Yang, arXiv:2012.07047 .[55] C. J. Cramer, Essentials of computational chemistry: theories and models , John Wiley & Sons, .[56] Y.-X. Yao, F. Zhang, C.-Z. Wang, K.-M. Ho, P. P. Orth, arXiv:2003.04211 .[57] N. Lanat`a, Y.-X. Yao, C.-Z. Wang, K.-M. Ho, G. Kotliar,
Phys. Rev. X , , 1 011008.[58] N. Lanat`a, Y. Yao, X. Deng, V. Dobrosavljevi´c, G. Kotliar, Phys. Rev. Lett. , Physical review letters , , 18 186404.[60] G. Knizia, G. K.-L. Chan, Journal of chemical theory and computation , , 3 1428.[61] Q. Sun, G. K.-L. Chan, Acc. Chem. Res. , , 12 2705.[62] T.-H. Lee, T. Ayral, Y.-X. Yao, N. Lanata, G. Kotliar, Phys. Rev. B , , 11 115129.[63] J. MacDonald, Phys. Rev. , , 9 828.[64] L.-W. Wang, A. Zunger, J. Chem. Phys. , , 3 2394., 3 2394.