A quantum algorithm for obtaining the lowest eigenstate of a Hamiltonian assisted with an ancillary qubit system

# A quantum algorithm for obtaining the lowest eigenstate of a Hamiltonian assisted with an ancillary qubit system

Jeongho Bang Center for Macroscopic Quantum Control, Department of Physics and Astronomy, Seoul National University, Seoul, 151-747, Korea Department of Physics, Hanyang University, Seoul 133-791, Korea    Seung-Woo Lee Center for Macroscopic Quantum Control, Department of Physics and Astronomy, Seoul National University, Seoul, 151-747, Korea    Chang-Woo Lee Department of Physics, Texas A&M University at Qatar, PO Box 23874, Doha, Qatar Center for Macroscopic Quantum Control, Department of Physics and Astronomy, Seoul National University, Seoul, 151-747, Korea    Hyunseok Jeong Center for Macroscopic Quantum Control, Department of Physics and Astronomy, Seoul National University, Seoul, 151-747, Korea Centre for Quantum Computation and Communication Technology, School of Mathematics and Physics, University of Queensland, Brisbane, Queensland 4072, Australia
July 29, 2019
###### Abstract

We propose a quantum algorithm to obtain the lowest eigenstate of any Hamiltonian simulated by a quantum computer. The proposed algorithm begins with an arbitrary initial state of the simulated system. A finite series of transforms is iteratively applied to the initial state assisted with an ancillary qubit. The fraction of the lowest eigenstate in the initial state is then amplified up to . We prove that our algorithm can faithfully work for any arbitrary Hamiltonian in the theoretical analysis. Numerical analyses are also carried out. We firstly provide a numerical proof-of-principle demonstration with a simple Hamiltonian in order to compare our scheme with the so-called “Demon-like algorithmic cooling (DLAC)”, recently proposed in [Nature Photonics 8, 113 (2014)]. The result shows a good agreement with our theoretical analysis, exhibiting the comparable behavior to the best ‘cooling’ with the DLAC method. We then consider a random Hamiltonian model for further analysis of our algorithm. By numerical simulations, we show that the total number of iterations is proportional to , where is the difference between the two lowest eigenvalues, and is an error defined as the probability that the finally obtained system state is in an unexpected (i.e. not the lowest) eigenstate.

03.67.Ac

## I Introduction

It is essential to deal with complex systems whose dynamics are typically described by many-body Hamiltonians in various fields of quantum information and computation. One important issue that has been raised in quantum computation is the eigen-problem of a Hamiltonian whose Hilbert-space is very large (possibly, even several hundreds or thousands). In particular, one of the frequently faced but quite formidable problem 111Such a task belongs to the class of “Non-deterministic Polynomial” (NP), or its quantum generalization, called “Quantum-Merlin-Arthur” (QMA) Liu07 (). is how to find the lowest energy eigenstate of a complex system Fischer91 (); Rojas96 (); Albert02 ().

A possible way to attack this problem is to use an approach of so-called algorithmic quantum cooling (AQC), which is often referred to as a systematic technique of descreasing (increasing) the fractions of the higher (lower) energy eigenstates of the system Boykin02 (); Fernandez04 (); Verstraete09 (); Xu14 (). The AQC is particularly useful when it is desired to initialize the part of higher energy states, say ‘hot’ qubits (e.g., a macroscopic number of spins Gershenfeld97 ()), to lower energy states in ensemble quantum computation King98 (); Roos99 (); James00 (); Popp06 (); Christ07 () or quantum simulation Lewenstein07 (); Bloch12 (). Most recently, a simple but powerful ‘pseudo’ AQC method, “Demon-like algorithmic cooling (DLAC)”, has been proposed and experimentally demonstrated Xu14 (). The DLAC method can drive a given initial state to the lowest eigenstate using a quantum-circuit module with an ancillary qubit. The core process in this method is the measurement of the ancillary qubit to discard the ‘heated’ state and to leave only the ‘cooled’ state (like Maxwell’s famous ‘demon’ Leff90 (); Lloyd97 ()).

In this work, we propose a quantum algorithm to obtain the lowest energy eigenstate of an arbitrary Hamiltonian simulated by a quantum computer. Similarly to some AQC methods Boykin02 (); Xu14 () (or other variational methods in classical computation Kosloff86 (); Lehtovaara07 ()), we start from an arbitrary initial state of the simulated system that (usually) contains very small fractions of the lower energy eigenstates. Applying a finite series of transformations, each of which consists of the quantum Householder reflection and the unitary of the system dynamics, the initial state is allowed to evolve amplifying the fraction of the lowest energy eigenstate up to . We note that our algorithm also employs an ancillary qubit system, similarly to DLAC method. However, no measurements are performed on the ancillar-qubit during the algorithm process; namely, Demon is not needed in our algorithm. The measurement is performed only once at the end of the algorithm to get the final state of the system removing the ancillary qubit 222We can also observe the cooling-like and the heating-like behaviors when we consider the whole system Hamiltonian involving the ancillary qubit system (See Appendix B for details).. In our theoretical analysis, we prove that our algorithm faithfully works for any given Hamiltonian. Numerical analyses are also carried out. Firstly, we provide a numerical proof-of-principle demonstration for a simple Hamiltonian whose eigenvalues are equally spaced. The result is quite consistent with our theoretical analysis, and in particular it exhibits the behavior comparable to the best ‘cooling’ available with the DLAC method. We then consider a random Hamiltonian model for further analysis of our algorithm. We presume that the required iterations to achieve an accuracy () is proportional to with , where is difference of the two lowest eigenvalues, and is tolerable error defined as the probability that the finally obtained state of the system is to be unexpected (i.e. not the lowest) eigenstate. By numerical simulations, it is found that and .

## Ii Problem & Method

To begin, we state the problem as follows: Consider a Hamiltonian of -dimensional Hilbert-space. The energy eigenvalues of are scaled as

 0<λ0<λ1≤…≤λN−1≤1, (1)

where we assume no perfect degeneracy between and (i.e. ). The energy eigenstate associated with the energy eigenvalue is given as . Note that the energy eigenvalues and the eigenstates are completely unknown. In this circumstance, the problem that we focus on here is: How to obtain the lowest energy eigenstate .

In solving this problem, we start with an initial system state in -dimensional Hilbert-space of the given Hamiltonian,

 |φ0⟩=N−1∑j=0aj∣∣vj⟩, (2)

where the computational bases and the coefficients are known to us. We assume that the fraction is not equal to zero but vanishingly small. We then consider an ancillary system of a clean qubit. By adopting an ancilla-qubit state , we prepare a composite initial state such that

 |ψ0⟩=|ϕ0⟩⊗|φ0⟩=|0⟩+|1⟩√2⊗N−1∑j=0aj∣∣vj⟩. (3)

Here, if we expand in terms of the eigenstates , the composite initial state is rewritten as

 |ψ0⟩=N−1∑k=0γ0,k√2(|0,λk⟩+|1,λk⟩), (4)

where , , and . Here, . Note that the coefficients cannot be evaluated, because the energy eigenstates are unknown.

Then, we construct a finite number of transformations, , according to the following recursive relation: For ,

 ^Ri=^Ti^Ri−1^T†i, and ^Ti=^Ri−1^U^Ri−1^U†, (5)

where is “quantum Householder reflection”, defined as . Such an operation has widely been used in quantum search Grover97 (), or other tasks Ivanov06 (); Ivanov07 (); Ivanov08 (). We also use a -dimensional unitary , defined as

 ^U=(|0⟩⟨0|⊗^A(τ))+i(|1⟩⟨1|⊗^A(τ)†), (6)

where is a unitary of the system’s dynamics 333Here, we omit the conventional minus (‘’) sign in the exponent., and is a scaling constant concerning the time of Hamiltonian action.

Then, the process of our algorithm can simply be thought of as a series of transformations applied on the initial composite state such that (See Fig. 1):

 ^Tnc…^T2^T1|ψ0⟩=|ψnc⟩. (7)

At the final step of , measurement is performed only once on the ancillary qubit. Tracing out the ancilar-qubit state by the measurement, the remaining state of the system is supposed to be close to the lowest energy eigenstate .

## Iii Theoretical analysis

We now analyze the process of our algorithm. By using Eqs. (4)-(7), we can describe the evolution of the composite state (See Appendix A):

 |ψ0⟩ ^T1⟶ |ψ1⟩=N−1∑k=0γ0,k√2(γ1,k|0,λk⟩+γ∗1,k|1,λk⟩) (8) ^T2⟶ |ψ2⟩=N−1∑k=0γ0,k√2(γ1,kγ2,k|0,λk⟩+γ∗1,kγ∗2,k|1,λk⟩) ⋮ ^Ti⟶ |ψi⟩=N−1∑k=0γ0,k√2(Γi,k|0,λk⟩+Γ∗i,k|1,λk⟩) ⋮

where (thus, ), and is given as a function of th-ordered energy eigenvalue (Here, ). The explicit form of is given as

 γi,k=(4|Wi−1|2−1)−2|Wi−1|e−iπ(1−τλk)4. (9)

where the factor is defined as

 |Wi|=∣∣⟨ψi|^U|ψi⟩∣∣=N−1∑k=0|γ0,k|2|Γi,k|2cosπ(1−τλk)4. (10)

Then, after sufficiently large step, we perform a von-Neumann measurement with on the ancillary qubit. The ancillary qubit is then removed, and we get the final state involving the fractions of the energy eigenstates. More specifically, the final state is given as follows (Note that ).

 |φnc⟩→⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩N−1∑i=0γ0,kΓnc,k|λk⟩ if |0⟩ is measured (with 12 % probability),N−1∑i=0γ0,kΓ∗nc,k|λk⟩ if |1⟩ is measured (with 12 % probability), (11)

We now show that the final state in Eq. (11), either or , becomes close to in the limit of . In particular, we will show that

 f0(λ0)≤f1(λ0)≤⋯

To this end, we state and sketch the proof of the following two propositions: [P.1] The coefficient factor of the lowest eigenstate is always larger than or equal to , and thus, for all . [P.2] The other probabilities of the higher eigenstates (excepting the lowest one) goes to zero with when . The proof of the latter is particularly important to guarantee that can go up to .

First, we provide the proof of [P.1]. To start, let us calculate from Eq. (9). With the simple algebra, we have

 (13)

where

 Δi−1,k=|Wi−1|−cosπ(1−τλk)4. (14)

By observing Eq. (13), we can know that the proof of [P.1], i.e. , is straightforward if

 {(a)|Wi−1|2≥12 (∀i),(b) Δi−1,0≥0 (∀i). (15)

Note that is larger than . To verify Eq. (15), we give the theoretical lower and upper bound on the value of , by using Eq. (10), as

 cosπ(1−τλ0)4≤|Wi−1|=⟨cosπ(1−τλk)4⟩i−1≤cosπ(1−τλN−1)4, (16)

where ‘’, which is an expectation value of at ()th step. Here, the (mathematical) condition to meet the lower bound is that and , whereas the upper bound is given when and . Note that Eq. (16) is always satisfied for the given eigenvalues () scaled as in Eq. (1). Thus, () () is true (because ). Then, by applying the lower bound value ‘’ to Eq. (14) with , we can directly verify that () () also holds. Therefore, [P.1] is always the case.

Next, let us consider [P.2]. To proceed, we assume that the eigenvalues () are divided into the two groups and at any ()th step, each of which is characterized by

 {g1:fi−1(λk)≤fi(λk) for λk≤ξi−1,g2:fi−1(λk)>fi(λk) for λk>ξi−1, (17)

where is a boundary factor of the ()th step. Then, by using Eqs. (13)-(15), we can show that

 |γi,k|2≥1⟺Δi−1,k≥0⟺λk≤ξi−1, (18)

where the explicit form of the boundary factor can be found as

 ξi−1=1τ(1−4πarccos|Wi−1|). (19)

Here, noting [P.1] and the property , we can verify that (), since the increment of necessarily results in the decrements of any other probabilities . This allows us to prove that [by using Eq. (19)]

 ξi≤ξi−1 for all i=1,2,…,∞. (20)

By using Eqs. (16) and (19), we can give the theoretical lower and upper bound of the boundary factor , for any , as

 λ0≤ξi≤λN−1. (21)

We note that the theoretical lower bound in Eq. (21) is an alternative expression of [P.1]; namely, it is true that the lowest eigenvalue always belongs to the group . By using Eq. (19), we can also obtain that 444Note that , and ().

 ξi−1−ξi = (22) = 4τπ∞∑l=0(2l)!4l(2l+1)(l!)2(|Wi−1|2l+1−|Wi|2l+1).

Here, (mathematically) if , then , i.e. does not converge but becomes smaller as increasing [as in Eq. (20)]. However, physically, must converge to a finite value larger than or equal to with when . Thus, we assume that is converged to a value in between and , where is a specific (integer) number of the eigenvalue index. From this, we assume further that, for ,

 fi(λk≤k′)→uk and fi(λk>k′)→0, (23)

where . Then, by using Eq. (10), we obtain that

 |Wi−1|−|Wi|→∑k≤k′(1−|γi,k|2)ukcosπ(1−τλk)4, for i→∞. (24)

In the circumstance, we can infer that the only possible solution for the physically reasonable condition (i.e. ) is given by

 |γi,0|2→1 & ∣∣γi,k≠0≤k′∣∣2<1, and, uk≠0≤k′=0, (25)

when . This solution yields that (for ),

 ⎧⎨⎩fi(λ0)→u0=1 (% because ∑k≤k′uk=1),|Wi|→cosπ(1−τλ0)4, % and thus, ξi→λ0 [from Eqs.~{}(???), % and (???)], (26)

which are consistent with Eq. (25) again 555i.e. from , it is verified that [from Eq. (14)], and thus, [as in Eq. (25)].. Here, if we assume any nonzero value(s) of with in Eq. (23) (and hence, ), it is in contradiction to the assumption (For detailed proof, see Appendix B). On the basis of the above description, it is certain that all the probabilities of the higher eigenstates, excepting the lowest one, go to zero [i.e. from Eq. (23) &  from Eq. (25)] when .

Therefore, we can show that the fraction of the lowest energy eigenstate only reaches close to unity after a sufficiently large number of the iterations [as in Eq. (12)], whereas all other probabilities go to zero.

## Iv Numerical analysis

We firstly provide a numerical proof-of-principle demonstration for a simple physical system where the energy levels are equally spaced (sometimes, called “Wannier-Stark ladder” Wannier60 (); Mendez88 (); Voisin88 ()). We thus consider a -dimensional Hamiltonian with the equidistant eigenvalues (). Here, the ground-state energy has a finite value less than . In the simulation, we set the dimension of the Hilbert-space as , and, for simplicity, the initial system state is chosen with for all . Note that the chosen state contains a very small fraction of the ground state, i.e. . In Fig. 2(a), we plot and with increasing the iteration , where reaches from to close to , and decays down to . For a comparison, we also include the data of and in Fig. 2(a), assuming that we use the DLAC method. Here we consider, particularly, the best ‘cooling’ available with the DLAC method (i.e. the case where the cooled results are only appeared in the Demon’s measurements). We note that, in DLAC method, the ‘cooling’ factor (similar to the factor in our algorithm) is given, for the system dynamics , as Xu14 ()

 1−sinξk≈e−π4λkτ, (27)

where is associated with the energy eigenvalues , and the approximation is done for small time evolution 666This ‘cooling’ behavior is also similar to that of a classical method, called “imaginary time propagation (ITP)” Kosloff86 (); Lehtovaara07 (). Actually, if we consider that the given Hamiltonian is diagonalized, we can associate these two methods (See Supplementary Information of Ref. Xu14 ()).. Thus, we adopt the same initial state to make the comparison as convincing as possible. It is observed that, in this case, also grows up to , and the behavior is quite similar to that of our algorithm. We also plot the coefficients and in Fig. 2(b), where is always larger than , whereas becomes smaller than leading to the decrease of . These are consistent with the above theoretical analyses in Eq. (25).

From the analyzed behaviors in the previous section and already known classical methods Abrams99 (); Yung12 (), we presume here that the number of iterations for an accuracy () is dominated by the difference of the two lowest eigenvalues, , and the tolerable error . If we assume that and are so small that all of the higher-order terms (e.g. , and ) can be negligible, we can explicitly calculate that is upper bounded as . Thus, more generally, we conjecture that with . Here, is a constant factor.

With the above prediction in mind, we consider a model of randomly generated Hamiltonian for more general analysis. We perform numerical simulations, and find for a given accuracy level. In the simulations, number of energy eigenvalues are randomly generated 777Here, we construct the random Hamiltonian in such a way: Firstly, we make a diagonal matrix with a randomly generated energy eigenvalue (but, is always be a certain predetermined value). We, then, construct an Hamiltonian by rotating such that , where the unitary is given by the randomly chosen real parameter vector and SU() group generators Hioe81 (). in (], but and are chosen such that the difference becomes a function . Here we consider three cases: . The iterations are continued until (or equivalently, ). In Fig. 3(a), we present versus graphs on a log-log scale. Each data is averaged over simulations. The data are fitted to , and we find that (, ) when , (, ) when , and (, ) when . Note here that the fitting parameters are very well matched to the parameter . These results allow us to estimate the value of : i.e. .

We also perform numerical simulations for the random Hamiltonian model to investigate for different accuracy (i.e. ) condition. In the simulations, we set as , and consider three cases by taking , , and . In Fig. 3(b), we give the graphs of versus on a log-log scale. Each data of point is also averaged over simulations. The data are well fitted to . The fitting parameters and that we found are: (, ) when , (, ) when , and (, ) when . The parameters of are similar to in all the cases; thus, we estimate .

## V Summary & Discussion

We have proposed a quantum algorithm to obtain the lowest energy eigenstate of a Hamiltonian. The process of the proposed algorithm could simply be regarded as a finite series of the transformations that consist of quantum Householder reflections and unitary of the system dynamics. Our algorithm is also assisted with an ancillary qubit system, similarly to the recently proposed method called “Demon-like algorithmic (pseudo) cooling (DLAC)”, but, in our algorithm, the measurement is performed only once at the end of the algorithm to remove the ancillary qubit.

In the theoretical analysis of Sec. III, we proved that our algorithm can faithfully work for any given Hamiltonian. We also carried out numerical analyses in Sec. IV. First, we provided a numerical proof-of-principle demonstration for a simple Hamiltonian whose eigenvalues were equally spaced. The results showed a good agreement with our theoretical analysis, exhibiting the comparable behavior to the best ‘cooling’ available with the DLAC method. We then considered a random Hamiltonian model for further analysis. We presumed that the total iterations required for an accuracy would be proportional to with , where was difference between the two lowest eigenvalues, and is a constant factor. In the simulations, we estimated that and .

In addition, we may compare our algorithm with the “pseudo-cooling” method Xu14 (), which is also aimed to increase the fractions of the lower energy eigenstates, even though the temperature could not be well-defined in the process. In contrast to the method described in Ref. Xu14 () using a demon-like selective process by measurement, our algorithm does not contain any post-selection by measurement, but it requires a measurement at the end of the algorithm to remove the ancillary qubit system (See Appendix C for detailed argument).

Our algorithm may be useful to deal with many problems arising in the studies of complex systems, which may require to reach the lowest eigenstate of a given Hamiltonian.

###### Acknowledgements.
The authors thank Sunwhan Jo, Chanhyoup Lee, and Junghee Ryu for helpful discussions. We acknowledge the support of the Basic Science Research Program through National Research Foundation of Korea (NRF) funded by the Ministry of Science, ICT & Future Planning (No. 2010-0018295 and No. 2010-0015059).

## Appendix A Proof of Eqs. (8)-(10)

In order to prove Eqs. (8)-(10), let us consider the transformation at the first step. From Eq. (5), we can write as

 ^T1 = ^R0^U^R0^U† (28) = (^11−2|ψ0⟩⟨ψ0|)(^11−2^U|ψ0⟩⟨ψ0|^U†) = ^11−2^U|ψ0⟩⟨ψ0|^U†−2|ψ0⟩⟨ψ0|+4⟨ψ0|^U|ψ0⟩|ψ0⟩⟨ψ0|^U†.

Applying the above to the composite initial state [as in Eq. (4)], the output state is computed such that

 |ψ1⟩ = |ψ0⟩−2W∗0^U|ψ0⟩−2|ψ0⟩+4|W0|2|ψ0⟩ (29) = (4|W0|2−1)|ψ0⟩−2W0^U|ψ0⟩ = N−1∑k=0γ0,k√2{[(4|W0|2−1)−2W∗0eiπ4τλk]|0,λk⟩ +[(4|W0|2−1)−2W∗0ei(π2−π4τλk)]|1,λk⟩},

where we let . Here, we represent the complex number as the polar form,

 W0 = N−1∑k=0|γ0,k|2⎛⎜⎝eiπ4τλk+ei(π2−π4τλk)2⎞⎟⎠ (30) = N−1∑k=0|γ0,k|2(1√2cosπτλk4+1√2sinπτλk4)(1+i√2) = N−1∑k=0|γ0,k|2[cosπ(1−τλk)4]eiπ4.

Thus, using the above Eq. (30), we rewrite the result of Eq. (29) as

 |ψ1⟩=N−1∑k=0γ0,k√2(γ1,k|0,λk⟩+γ∗1,k|1,λk⟩), (31)

where .

Based on the results, we can generalize the above-described computations to the higher th step (), using the Householder reflection and unitary . In such generalization, we can easily find the expression of the th output state [as in Eq. (8)] and the coefficients [as in Eq. (9)] with the factor [as in Eq. (10)].

## Appendix B The uniqueness of the solution Eq. (25)

As mentioned in the main text, one may consider a more general situation, where

 ∣∣γi,k≤k′∣∣2→1, and, uk≤k′≠0 for i→∞. (32)

Here, it is obvious that . This can yield that when , i.e. the probability of the lowest eigenstate cannot reach to .

In such a general assumption, we prove that Eq. (25) is the unique solution. To this end, we first give [using Eqs. (10) and (23)]

 |Wi|→∑k≤k′ukcosπ(1−τλk)4 (for i→∞). (33)

Thus, for any specific eigenvalue () and , we obtain [using Eq. (14)]

 Δi,l → ∑k≤k′ukcosπ(1−τλk)4−cosπ(1−τλl)4 (34) = ∑k≠l≤k′ukcosπ(1−τλk)4+(ul−1)cosπ(1−τλl)4 = ∑k≠l≤k′uk(cosπ(1−τλk)4−cosπ(1−τλl)4) = ∑k≤k′uk(cosπ(1−τλk)4−cosπ(1−τλl)4),

where we used , or . Here, we note that

 |γi,l|2→1, only when Δi,l→0 (for i→∞), (35)

which is verified by Eqs. (13) and (16). Then, from Eq. (34), we find that the solution in Eq. (32) is possible when all the following conditions are satisfied:

 u0(c0−c0)+u1(c1−c0)+ ⋯ +uk′(ck′−c0)=0 (for l=0), u0(c0−c1)+u1(c1−c1)+ ⋯ +uk′(ck′−c1)=0 (for l=1), ⋮ u0(c0−ck′)+u1(c1−ck′)+ ⋯ +uk′(ck′−ck′)=0 (for l=k′), (36)

where we let () just for convenience. We rewrite the above conditions by adding all () equations in Eq. (36) as

 u0⎛⎝(k′+1)c0−∑k≤k′ck⎞⎠+u1⎛⎝(k′+1)c1−∑k≤k′ck⎞⎠+⋯+uk′⎛⎝(k′+1)ck′−∑k≤k′ck⎞⎠=0. (37)

Therefore, we can represent the condition for the existence of the solution Eq. (32) as

 1k′+1∑k≤k′ck=c0=c1=⋯=ck′ for uk≤k′≠0. (38)

However, we can directly see that this Eq. (38) can never be satisfied, excepting the case and as in Eq. (25) (as long as [P.1] is true).

## Appendix C Comparison between our algorithm and “pseudo-cooling” method in Ref. Xu14 ()

Here we discuss that our algorithm is quite distinct from the “pseudo-cooling” method proposed in Ref.