Variational ansatz-based quantum simulation of imaginary time evolution

Imaginary time evolution is a powerful tool for studying quantum systems. While it is possible to simulate with a classical computer, the time and memory requirements generally scale exponentially with the system size. Conversely, quantum computers can efficiently simulate quantum systems, but not non-unitary imaginary time evolution. We propose a variational algorithm for simulating imaginary time evolution on a hybrid quantum computer. We use this algorithm to find the ground-state energy of many-particle systems; specifically molecular hydrogen and lithium hydride, finding the ground state with high probability. Our method can also be applied to general optimisation problems and quantum machine learning. As our algorithm is hybrid, suitable for error mitigation and can exploit shallow quantum circuits, it can be implemented with current quantum computers.

Imaginary time is an unphysical, yet powerful, mathematical concept.It has been utilised in numerous physical domains including: quantum mechanics, statistical mechanics, and cosmology.Often referred to as performing a 'Wick rotation' [1], replacing real time with imaginary time connects Euclidean and Minkowski space [2], quantum and statistical mechanics [3], and static problems to problems of dynamics [4].In quantum mechanics, propagating a wavefunction in imaginary time enables: the study of finite temperature properties [5][6][7], finding the ground state wavefunction and energy [8][9][10][11], and simulating real time dynamics [12,13].For a system with Hamiltonian, H, evolving in real time, t, the propagator is given by e −iHt .The corresponding propagator in imaginary time, τ = it, is given by e −Hτ ; a non-unitary operator.
Using a classical computer, we can simulate imaginary time evolution by evaluating the propagator and applying it to the system wavefunction.However, because the dimension of the wavefunction grows exponentially with the number of particles, classical simulation of manybody quantum systems is limited to small or specific cases [14].While efficient variational trial states have been developed for a number of applications [15], powerful trial wavefunctions typically require classical computational resources which scale exponentially with the system size [11].
Quantum computing can naturally and efficiently store many-body quantum states, and hence is suitable for simulating quantum systems [16].We can map the system Hamiltonian to a qubit Hamiltonian, and simulate real time evolution (as described by the Schrödinger equation) by realising the corresponding unitary evolution with a quantum circuit [17].Using Trotterization [18], the real time propagator can be decomposed into a sequence of single and two qubit gates [19].The ability to represent the real time propagator with a sequence of gates stems from its unitarity.In contrast, because the imaginary time operator is non-unitary, it is impossible to decompose it into a sequence of unitary gates using Trotterization, and thus directly realise it with a quantum circuit.As a result, alternative methods are required to implement imaginary time evolution using a quantum computer.
Classically, we can simulate real (imaginary) time evolution of parametrised trial states by repeatedly solving the (Wick-rotated) Schrödinger equation over a small timestep, and updating the parameters for the next timestep [8,9,11,[20][21][22][23].This method has recently been extended to quantum computing, where it was used to simulate real time dynamics [24].Closely related are the variational quantum eigensolver (VQE) [25][26][27][28][29][30][31] and the quantum approximate optimisation algorithm (QAOA) [32], which update the parameters using a classical optimisation routine, to find the minimum energy eigenvalue of a given Hamiltonian.As 'hybrid quantumclassical methods', these algorithms use a small quantum computer to carry out a classically intractable subroutine, and a classical computer to solve the higher level problem.The quantum subroutine may only require a small number of qubits and a low depth circuit, presenting a potential use for noisy intermediate-scale quantum hardware [33].
In this letter, we propose a method to simulate imaginary time evolution on a quantum computer, using a hybrid quantum-classical variational algorithm.The proposed method thus combines the power of quantum computers to efficiently represent many-body quantum states, with classical computers' ability to simulate arbitrary (including unphysical) processes.We discuss using this method to find the ground state energy of manybody quantum systems, and to solve optimisation problems.We then numerically test the performance of our algorithm at finding the ground state energy of both the Hydrogen molecule (H 2 ) and Lithium Hydride (LiH).We compare our results for LiH to those obtained using the VQE with gradient descent optimisation.As our algorithm only requires a low depth circuit, it can be realised with current and near-term quantum processors.
Variational simulation of imaginary time evolution.-Wefocus on many-body systems that are described by Hamiltonians H = i λ i h i , with real coefficients, λ i , and observables, h i , that are tensor products of Pauli matrices.We assume that the number of terms in this Hamiltonian scales polynomially with the system size, which is true for many physical systems, such as molecules [34].Given an initial state |ψ , the normalised imaginary time evolution is defined by |ψ(τ ) = A(τ )e −Hτ |ψ(0) , where In the instance that the initial state is a maximally mixed state, the state at time τ is a thermal or Gibbs state ρ T =1/τ = e −Hτ /Tr[e −Hτ ], with temperature T = 1/τ .When the initial state has a non-zero overlap with the ground state, the state at τ → ∞ is the ground state of H. Equivalently, the Wick rotated Schrödinger equation is, where the term E τ = ψ(τ )|H|ψ(τ ) results from enforcing normalisation.Even if |ψ(τ ) can be represented by a quantum computer, the non-unitary imaginary time evolution cannot be naively mapped to a quantum circuit.
In our variational method, instead of directly encoding the quantum state |ψ(τ ) at time τ , we approximate it using a parametrised trial state |φ( θ(τ )) , with θ(τ ) = (θ 1 (τ ), θ 2 (τ ), . . ., θ N (τ )).This stems from the intuition that the physically relevant states are contained in a small subspace of the full Hilbert space [35].The trial state is referred to as the ansatz.In condensed matter physics and computational chemistry, a wide variety of ansatze have been proposed for both classical and quantum variational methods [11,16,36,37].
Using a quantum circuit, we prepare the trial state, |φ( θ) , by applying a sequence of parametrised unitary gates, V ( θ) = U N (θ N ) . . .U k (θ k ) . . .U 1 (θ 1 ) to our initial state, | 0 .We express this as |φ( θ) = V ( θ) | 0 and remark that V ( θ) is also referred to as the ansatz.We refer to all possible states that could be created by the circuit V as the 'ansatz space'.Here, U k (θ k ) is the k th unitary gate, controlled by parameter θ k , and the gate can be regarded as a single or two qubit gate.
To simulate the imaginary time evolution of the trial state, we use McLachlan's variational principle [38,39], δ (∂/∂τ + H − E τ ) |ψ(τ ) = 0, where ρ = Tr[ ρρ † ] denotes the trace norm of a state.By replacing |ψ(τ ) with |φ(τ ) = |φ( θ(τ )) , we effectively project the desired imaginary time evolution onto the manifold of the ansatz space.The evolution of the parameters is obtained from the resulting differential equation where ∂θi h α |φ(τ ) , and h α and λ α are the Pauli terms and coefficients of the Hamiltonian, as described above.The derivation of Eq. ( 2) can be found in the Supplementary Materials.As both A ij and C i are real, the derivative θj is also real, as required for parametrising a quantum circuit.Interestingly, although the average energy term E τ appears in Eq. ( 1), it does not appear in Eq. ( 2).This is because the ansatz applied maintains normalisation, as it is composed of unitary operators.
Imaginary time evolution with quantum circuits.-Byfollowing a similar method to that introduced in Refs.[24,40,41], we can efficiently measure A ij and C i using a quantum computer.
We assume that the derivative of a unitary gate U i (θ i ) can be expressed as , with unitary operator σ k,i .The derivative of the trial state is given by ∂ |φ(τ There are typically only one or two terms resulting from each derivative.As an example, when All of these terms are of the form aℜ(e iθ 0| U | 0 ) and can be evaluated using the circuit in Supplementary Materials.
With A(τ ) and C(τ ) at time τ , the imaginary time evolution over a small interval δτ can be simulated by evaluating ˙ θ(τ ) = A −1 (τ ) • C(τ ), and using a suitable update rule, such as the Euler method, (3) By repeating this process N T = τ total /δτ times, we can simulate imaginary time evolution over a duration τ total .When there are redundant parameters, the rank of A can be less than the number of parameters, which suggests a non-unique solution for Eq. ( 2).This can be circumvented by fixing the values of some of the parameters until A becomes invertible [11].This can be implemented algorithmically using truncated singular value decomposition.We refer to the Supplementary Materials for further discussion of the invertibility of A.
Ground state energy via imaginary time evolution.-Weapply our method to the problem of finding the ground state energy of a many-body Hamiltonian, H.As with the VQE, our goal is to find the values of the parameters, θ, which minimise the expectation value of the Hamiltonian where | 0 is our variational trial state.The VQE solves this problem by using a quantum computer to construct a good ansatz and measure the expectation value of the Hamiltonian, and a classical optimisation routine to obtain new values of the parameters.To date, the VQE has been used to experimentally find the ground state energies of several small molecules [25][26][27][28][29][30][31].In order to preserve the exponential speedup of the VQE over classical methods, the trial state is constructed using a number of parameters that scales polynomially with the system size.However, because we may need to consider many possible values for each parameter, the total size of the parameter space still scales exponentially with the system size.Moreover, many optimisation algorithms, such as gradient descent, are liable to becoming trapped in local minima.This combination can make the classical optimisation step of the VQE very difficult [42].
As described above, if the initial state has a non-zero overlap with the ground state, the propagation in imaginary time will evolve the system into the ground state, in the limit that τ → ∞.Classically, this has been leveraged as a powerful tool to find the ground state energy of quantum systems [8,9,11].Using our method, we can efficiently simulate imaginary time evolution to find the ground state, using a quantum computer.In the numerical simulations described below, we use the Euler method to solve differential equations, which corresponds to the update rule for the parameters shown in Eq. ( 3).We prove in the Supplementary Materials that when δτ is sufficiently small, the average energy of the trial state, E(τ ) = φ(τ )| H |φ(τ ) , always decreases when following the Euler update rule; E(τ + δτ ) ≤ E(τ ).
The update rule for the conventional gradient descent method is given by where a is a constant scaling factor, G(τ ) = −∇E(τ ) is the gradient of E(τ ), and C(τ ) ≡ −∇E(τ ) is the same vector in Eq. ( 2).The distinction between these two update rules is the presence of the A matrix in Eq. (3).Thus, the gradient descent method only considers information about the average energy, and not about the ansatz itself, which is encoded in A. We see from the toy problem depicted in Fig. 1 that variational evolution in imaginary time is qualitatively different to gradient descent.While gradient descent becomes trapped in local minimum as a result of moving downhill, imaginary time moves in a direction which takes it to the global minimum.We now explore the consequences for problems of practical interest.
Y is a controlled Y rotation with control qubit 0 and target qubit 1, R 0 X (θ1) is a rotation of qubit 0 along the X axis, and the rotation along the j axis is Rσ j (θ) = e −iθσ j /2 with Pauli matrices σj. ).Given a random initial parameter, gradient descent moves downhill and becomes trapped in a local minimum, unlike imaginary time evolution.
Simulation of H 2 and LiH.-Consider the problem to find the ground state energy of the H 2 and LiH molecules in their minimal spin orbital basis sets.We compare the LiH results to those obtained using the VQE, with gradient descent used for parameter optimisation.Gradient based methods have previously been found to perform well in the VQE, outperforming direct search optimisation algorithms (when hardware and algorithmic errors are neglected) [40].We calculate the gradient directly using a quantum circuit [40,41], which eliminates errors arising from a finite difference approximation of the gradient.We first map the fermionic Hamiltonian into a qubit Hamiltonian using the Bravyi-Kitaev encoding [43] and utilise a procedure for Hamiltonian reduction [37,44,45].We first focus on the H 2 results.After encoding, the two qubit Hamiltonian of H 2 [27] is where g i are real coefficients determined by the distance R between the two nuclei, and X, Y , Z denote the Pauli matrices.Here, we consider an internuclear distance of R = 0.75 Å and hence g 0 = 0.2252, g 1 = 0.3435, g 2 = −0.4347g 3 = 0.5716, g 4 = 0.0910, g 5 = 0.0910.There are numerous possible choices for the ansatz circuit, and we consider two of these for H 2 : the unitary coupled cluster (UCC) [37] and hardware efficient [31] ansatze.The UCC ansatz constructs a chemically motivated trial state, |φ(θ) = e −iθX0Y1 |0 1 |1 0 , which can be prepared using the quantum circuit in Fig. 2(a).In The quantum circuits for the H2 molecule in the minimal spin orbital basis set.Measurements are in the Z basis, and rotation gates are applied before measurement to measure in the X or Y basis.(a) The quantum circuit [27] for preparing the UCC ansatz.Here, R .In practice, the gates in the dashed box may be omitted.The other terms of A and C can be measured using similar circuits.
contrast, the hardware efficient ansatz attempts to generate a flexible wavefunction using as few gates as possible.There are different ways of designing a hardware efficient ansatz, and the ansatz we use is shown in Fig. 2(c).
As the UCC ansatz has only one parameter θ, both A and C in Eq. ( 2) are scalars.The value of C can be measured using the circuit in Fig. 2(b), and A = 0.25.As A is a constant, the imaginary time evolution is equivalent to the gradient descent method, with the gradient calculated by a quantum circuit [40,41].The optimal value of the parameter is θ opt = 5.2146 and the minimum energy is E min = −1.1457Hartree.When consid- ering the hardware efficient ansatz, we have eight parameters, and A is no longer a scalar.The simulation results of the variational imaginary time evolution are shown in Fig. 3(a).We see that the variational imaginary time evolution closely approximates the exact evolution and converge to the ground state energy for both ansatze.
We now discuss the LiH results.In the minimal basis, the LiH molecule has 12 spin-orbitals, but the Hamilto-nian can be reduced to act on six qubits, as described in Ref. [46] and the Supplementary Materials.We use a hardware efficient ansatz shown in the Supplementary Materials for our simulation, with 42 parameters.We use random initial values for these parameters.This may cause the trial state to have an exponentially small overlap with the ground state.Moreover, hardware efficient ansatze have recently been shown to become exponentially less powerful as the number of qubits increases [47].However, we do not expect either of these issues to be significant in our six qubit simulations.Furthermore, we have chosen a hardware efficient ansatz because of its greater experimental feasibility, compared with a UCC ansatz for six qubit LiH (although we note that an experimental simulation of simplified, three qubit LiH with a UCC ansatz has recently be performed [46]).In addition, we consider starting from a random initial state to be a more thorough test of both methods than starting from a good initial state, such as the Hartree-Fock state.We consider an internuclear distance of R = 1.45 Å, for which the exact ground state energy is E min = −7.8807Hartree.We compare the performance of variational imaginary time evolution with time step δτ = 0.1 with that of the VQE with gradient descent with step-size aδτ = 0.3.We applied both methods to 112 different trials, each initialised with different, random starting values of the parameters.As exemplified by one specific instance, shown in Fig. 3(b), we find that our method can almost always (with more than 90 % probability) find the ground state energy, while gradient descent often becomes trapped in local minima (as shown in the Supplementary Materials).
Discussion.-In this letter, we have proposed a method to efficiently simulate imaginary time evolution using a quantum computer.Our hybrid algorithm combines the ability of quantum computers to construct parametrised trial wavefunctions, with the ability of classical computers to update the parameter values.We have applied our method to finding the ground state energy of quantum systems, and have tested its performance on H 2 and LiH.As imaginary time evolution found the ground state with a significantly higher probability than gradient descent, we believe our method provides a competitive alternative to the VQE.We expect that our method would also be suitable for solving general optimisation problems, in place of the QAOA.
Our method can also be used to prepare a thermal (Gibbs) state, ρ T = e −H/T /Tr[e −H/T ] of Hamiltonian H at temperature T .Sampling from a Gibbs distribution is an important aspect of many machine learning algorithms, and so we believe that our method is applicable to problems in quantum machine learning.Moreover, while previous methods to prepare the Gibbs state [48,49] require long gate sequences (and hence, fault tolerance), our method can be implemented using a shallow circuit.Our algorithm can also be combined with recently proposed error mitigation techniques [24,50,51], which pro-tect shallow circuits from error accumulation, without the additional qubit resources required for full error correction.As a result, it is suitable for implementation on current and near-term quantum hardware.
Although exact imaginary time evolution deterministically propagates a good initial state to the ground state in the limit that τ → ∞, our variational method may still become trapped in local minima, if the chosen ansatz is not sufficiently powerful.In future works, we will investigate how our method may be optimally applied to a variety of tasks in chemistry, optimisation and machine learning.This will include developing suitable ansatze for a range of problems.Another interesting avenue for further study is a more general imaginary time evolution, e −f (H,τ ) , where f is a function, such as f (H, τ ) = (H − λ) 2 τ .This may enable the investigation of excited states, or expedite convergence to the ground state.

Variational simulation of imaginary time evolution
McLachlan's variational principle [38], applied to imaginary time evolution, is given by where and E τ = ψ(τ )|H|ψ(τ ) .For a general quantum state, McLachlan's variational principle recovers the imaginary time evolution If we consider a subspace of the whole Hilbert space, which can be reached using the ansatz |φ(τ ) = |φ(θ 1 , θ 2 , . . ., θ N ) , we can project the imaginary time evolution onto the subspace using McLachlan's variational principle.Replacing |ψ(τ ) with |φ(τ ) , yields Focusing on θi , we obtain Considering the normalisation condition for the trial state |φ(τ ) , we have and the derivative is simplified to where McLachlan's variational principle requires which is equivalent to the differential equation of the parameters Denoting E(τ ) = φ(τ )| H |φ(τ ) , we can show that the average energy always decreases by following our imaginary time evolution algorithm, for a sufficiently small stepsize; The third line follows from the definition of C i ; the fourth line follows from the differential equation of θ; the last line is true when A −1 is positive.First, we show matrix A is positive.We consider an arbitrary vector x = (x 1 , x 2 , . . ., x N ) T , and calculate ∂θi , then the first term equals i,j Similarly, we can show that the second term is also nonnegative.Therefore, x † • A • x ≥ 0, ∀x and A is nonnegative.In practice, when A has eigenvalues with value zero, A is not invertible.However, in our simulation, we define the inverse of A to be only the inverse of the nonnegative eigenvalues.Suppose U is the transformation that diagonalises A, i.e., G i,j = (U AU † ) i,j = 0, ∀i = j.Then, we define G −1 by The inverse of A is thus defined by Because A has nonnegative eigenvalues, G, G −1 , and hence A −1 all have nonnegative eigenvalues.
Evaluating A and C with quantum circuits In this section, we review the quantum circuit that can efficiently evaluate the coefficients A and C introduced in Ref. [24,40,41].
Without loss of generality, we can assume that each unitary gate U i (θ i ) in our circuit depends only on parameter θ i (since multiple parameter gates can be decomposed into this form).Suppose each U i is a rotation or a controlled rotation gate, its derivative can be expressed by with unitary operator σ k,i and scalar parameters f k,i .The derivative of the trial state is with In practice, there are only one or two terms, f k,i σ k,i , for each derivative.For example, when , and the derivative of the trial state ∂ |φ(τ ) /∂θ i can be prepared by adding an extra Z gate with a constant factor −i/2.
Therefore, the coefficients A ij and C i are given by All the terms of the summation follow the general form aℜ(e iθ 0| U | 0 ) and can be evaluated by the circuit in Fig. 1 in the main text.
In practice, we do not need to realize the whole controlled-U gate and can instead use a much easier circuit.For example, for the term and ℜ(e iθ 0| Ṽ † k,i Ṽl,j | 0 ) can be measured by the circuit in Fig. 4. The terms for C can be measured similarly.

Computational chemistry background
One of the central problems in computational chemistry is finding the ground state energy of molecules.This calculation is classically intractable, due to the exponential growth of Hilbert space with the number of electrons in the molecule.However, it has been shown that quantum computers are able to solve this problem efficiently [34].The Hamiltonian of a molecule consisting of M nuclei (of mass M I , position R I , and charge Z I ) and N electrons (with position r i ) is Because the nuclei are orders of magnitude more massive than the electrons, we apply the Born-Oppenheimer approximation, and treat the nuclei as classical, fixed point charges.After this approximation, the eigenvalue equation is Hermitian, the X gates acting on the ancilla qubit can be also omitted.
we seek to solve (in atomic units) is given by where |ψ is an energy eigenstate of the Hamiltonian, with energy eigenvalue E.
To solve this problem using a quantum computer, we first transform it into the second quantised form.We project the Hamiltonian onto a finite number of basis wave functions, {φ p }, which approximate spin orbitals.Electrons are excited into, or de-excited out of, these orbitals by fermionic creation (a † p ) or annihilation (a p ) operators, respectively.These operators obey fermionic anti-commutation relations, which enforces the antisymmetry of the wavefunction, a consequence of the Pauli exclusion principle.In the second quantised representation, the electronic Hamiltonian is written as h pq a † p a q + 1 2 p,q,r,s h pqrs a † p a † q a r a s , with where x is a spatial and spin coordinate.This Hamiltonian in general contains N 4 SO terms, where N SO is the number of orbitals considered.This fermionic Hamiltonian must then be transformed into a Hamiltonian acting on qubits.This is achieved using the Jordan-Wigner, or Bravyi-Kitaev transformations, which are described in Ref. [52].

Hydrogen
In our simulations, we consider the Hydrogen molecule in the minimal STO-6G basis.This means that only the minimum number of orbitals to describe the electrons are considered.'STO-nG' means that a linear combination of n Gaussian functions are used to approximate a Slater-type-orbital, which describes the electron wavefunction.Each Hydrogen atom contributes a single 1S orbital.As a result of spin multiplicity, there are four spin orbitals in total.
We are able to construct the molecular orbitals for H 2 by manually (anti)symmetrising the orbitals.These are where the subscripts on the 1S orbitals denote the spin of the electron in that orbital, and which of the two hydrogen atoms the orbital is centred on.By following the procedure in Ref. [52], the qubit Hamiltonian for H 2 in the BK representation can be obtained.This 4 qubit Hamiltonian is given by As this Hamiltonian only acts off diagonally on qubits 0 and 2 [27,31], it can be reduced to which only acts on two qubits.

Lithium Hydride
In our simulations, we consider Lithium Hydride in the minimal STO-3G basis.The Lithium atom has 3 electrons, and so contributes a 1S, 2S, 2P x , 2P y and 2P z orbital to the basis, while the Hydrogen atom contributes a single 1S orbital.With spin multiplicity, this makes 12 orbitals in total.However, we are able to reduce the number of spin orbitals required by considering their expected occupation.This reduces the qubit resources required for our calculation.In computational chemistry, the subset of spin orbitals included in a calculation is called the active space.
We first obtain the one electron reduced density matrix (1-RDM) for LiH, using a classically tractable CISD (configuration interaction, single and double excitations) calculation.
There are only six rows and columns in the 1-RDM because the spin-up and spin-down orbitals have been combined.The diagonal elements of the 1-RDM are the occupation numbers of the corresponding canonical orbitals (the orbitals contributed by the individual atoms, eg.1S Li ).In order to reduce our active space, we first perform a unitary rotation of the 1-RDM, such that it becomes a diagonal matrix, This gives the 1-RDM in terms of natural molecular orbitals (NMOs).The diagonal entries are called the natural orbital occupation numbers (NOONs).The Hamiltonian of LiH must also be rotated, using the same unitary matrix used to diagonalise the 1-RDM.This is equivalent to performing a change of basis, from the canonical orbital basis to the natural molecular orbital basis.
As can be seen, the first orbital has a NOON close to two, and so is very likely to be doubly occupied.As a result, we 'freeze' this core orbital, and consider it to always be doubly filled.We can then remove any terms containing a † 0 , a 0 , a † 1 , a 1 from the LiH fermionic Hamiltonian.We also notice that the fourth orbital has a NOON close to zero.As a result, we assume that this orbital is never occupied by either a spin-up or spin-down electron, and so remove another two fermion operators from the Hamiltonian.This leaves a fermionic Hamiltonian acting on 8 spin orbitals.We then map this fermionic Hamiltonian to a qubit Hamiltonian, using the BK transformation.The BK transformation enables us to further reduce the Hamiltonian by two qubits, using the spin and electron number symmetries of the molecule [31,53].This results in a Hamiltonian for LiH which acts on 6 qubits.All of these steps were carried out using OpenFermion [54], an electronic structure package to transform computational chemistry problems into a form that is suitable for investigation using a quantum computer.A similar encoding and reduction procedure for LiH is also described in Ref. [46].

Hardware-efficient ansatz
A hardware-efficient ansatz was recently realised using superconducting qubits [31].Following a similar, but slightly different structure, the general hardware-efficient ansatz we considered in this work is shown in Fig. 5 and the ansatz for simulating LiH is shown in Fig. 6.

Depth 1
Depth 2 5. General hardware efficient ansatz.We repeat the circuit structure of the first block to depth M .In the first block, we set the parameters of the first Z rotation gates to be 0.

Additional LiH simulation results
The exact ground state energy of LiH, at an internuclear distance of R = 1.45 Å, is E min = −7.8807Hartree (when working in the STO-3G basis).We compare the performance of our variational imaginary time evolution method with that of gradient descent, with time step δτ = 0.01.We applied both methods to 112 different trials, each initialised with different, random starting values of the parameters.As shown in Fig. 7, our method can almost always (with more than 90% probability) find the ground state energy, while gradient descent often becomes trapped in local minima.

FIG. 1 .
FIG. 1. (Color online) Comparing imaginary time and gradient descent for finding the ground state energy of a two qubit Hamiltonian H = diag(1, 2, 3, 0).The ansatz used is |ψ(θ1, θ2) = CR 0,1 Y (θ2)R 0 X (θ1) |00 (where CR 0,1 σ j α = e −iασ j /2 .The parameter θ is updated in order to minimise the total energy.(b) The quantum circuit that estimates C in Eq. (2).Here, hi corresponds to one of the six terms in the Hydrogen Hamiltonian.As hi is a tensor product of Pauli matrices, the control-hi gate can be realised by applying the individual control-Pauli gates in turn.(c) The quantum circuit for preparing the hardware efficient ansatz with eight parameters.In this ansatz, one parameter is redundant.(d) The circuit to measure A2,7 = ℜ ∂ φ(τ )| ∂θ 2 ∂|φ(τ ) ∂θ 7

FIG. 3 .
FIG. 3. (Color online) Simulation results for H2 and LiH with random initial parameters.The blue line is the exact ground state energy.(a) Simulation result for H2.Time step is δτ = 0.05.The dashed lines are exact simulations of imaginary time evolution under Eq.(1).The dot points are simulation results of the variational imaginary time evolution.UCC: unitary coupled cluster; HE: hardware efficient.(b) Simulation results for LiH.Time step is δτ = 0.1 for variational imaginary time evolution (red dots) and aδτ = 0.3 for gradient descent (black dashed), which converges to a local minimum.
The 1-RDM for a distance of 1.45 Å is shown below,

FIG. 6 .
FIG.6.The hardware efficient ansatz for the LiH simulation results presented in this paper.In our simulation, we used this two block circuit for the LiH molecule.In total, there are 42 parameters.

FIG. 7 .
FIG. 7. (Color online) Simulation results of the 112 trials with different, random starting values of the parameters.(a) Histogram of the 112 trials.The variational imaginary time evolution reaches within Chemical Accuracy (10 −3 Hartree) of the true ground state E0 = −7.8807Hartree in 103 out of 112 trials.In contrast, gradient descent generally gets trapped in local minima.(b) Scatter graph of the 112 trials.The dashed line in the inset plot marks results which fall within Chemical Accuracy of the ground state energy.