Boosting quantum amplitude exponentially in variational quantum algorithms

We introduce a family of variational quantum algorithms, which we coin as quantum iterative power algorithms (QIPAs), and demonstrate their capabilities as applied to global-optimization numerical experiments. Specifically, we demonstrate the QIPA based on a double exponential oracle as applied to ground state optimization of the H 2 molecule, search for the transmon qubit ground-state, and biprime factorization. Our results indicate that QIPA outperforms quantum imaginary time evolution (QITE) and requires a polynomial number of queries to reach convergence even with exponentially small overlap between an initial quantum state and the final desired quantum state, under some circumstances. We analytically show that there exists an exponential amplitude amplification at every step of the variational quantum algorithm, provided the initial wavefunction has non-vanishing probability with the desired state and that the unique maximum of the oracle is given by λ1>0 , while all other values are given by the same value 0<λ2<λ1 (here λ can be taken as eigenvalues of the problem Hamiltonian). The generality of the global-optimization method presented here invites further application to other problems that currently have not been explored with QITE-based near-term quantum computing algorithms. Such approaches could facilitate identification of reaction pathways and transition states in chemical physics, as well as optimization in a broad range of machine learning applications. The method also provides a general framework for adaptation of a class of classical optimization algorithms to quantum computers to further broaden the range of algorithms amenable to implementation on current noisy intermediate-scale quantum computers.


I. INTRODUCTION
Quantum computers promise exponential speedup over classical counterparts in solving certain tasks [1].When fault-tolerant general-purpose quantum computers become available, adiabatic state preparation and quantum phase estimation may become the standard quantum routines for determining the ground-state energy of sophisticated physical Hamiltonians [2][3][4][5].However, such schemes are very costly in terms of required overhead and hence are not suitable for the current era of noisy intermediate-scale quantum (NISQ) hardware [6][7][8][9][10].This limitation of quantum computers today shifts central attention towards low-depth hybrid quantumclassical algorithms, known as NISQ algorithms [11][12][13][14].The variational quantum eigensolver (VQE) [15,16] serves as a prototypical example, as an algorithm that computes the expectation value of a Hamiltonian, which is measured on a quantum machine, resulting in a cost function with a set of variational parameters, which are optimized using classical computers.The process is repeated until the cost function reaches its local minimum.
On the other hand, the variational quantum simulator [17] has been proposed for hybrid quantum-classical simulations of quantum dynamics based on the McLachlan's variational principle [18,19], including quantum imaginary time evolution (QITE) to prepare ground states [18][19][20][21][22][23][24][25][26].Here, we introduce the "quantum iterative power algorithm" (QIPA) inspired by the variational quantum simulator to provide an efficient method to the general problem of global optimization with near term quantum computers.
We recognize the strategy of tensor-train IPA can be implemented on quantum computers to enable global optimization of an even broader class of optimization problems.In tensor-train IPA [57], the optimization cost function of interest is taken to be a potential energy surface.A density is initialized in the potential energy surface, and an oracle is iteratively applied in a sifting approach akin to imaginary time propagation (with infinite mass) to localize the density as a delta function at the global minimum position.The expectation value of position then gives the location of the global minimum.Tensor-train IPA represents the density and potential energy surface as tensor trains to avoid calculation of the cost function everywhere in search space, which is efficient for representation of problems amenable to low-rank representations, such as prime factorization or molecular geometry optimization [57].However, the tensor-train strategy faces the roadblock that highly-coupled systems cannot be efficiently represented in low-rank tensor-train format.In contrast, quantum computers excel in the simulation of highly-coupled systems, as the coupling or entanglement between qubits is limited only by the choice of ansatz [58][59][60][61][62][63][64].
The quantum iterative power algorithm (QIPA) takes advantage of the high degree of entanglement possible on quantum computers with a hybrid variational scheme.In standard variational approaches such as the variational quantum eigensolver (VQE) [12, 14-16, 20, 65-68], classical optimizers are used to determine the parameters of a quantum circuit, which are used to prepare trial wavefunctions measured to obtain expectation values.Analogously, the variational quantum simulator [17] evolves the parameters that define the time-evolved wavefunction by using a classical computer that integrates the Euler-Lagrange equation obtained from the Schrödinger equation with the McLachlan's variational principle.Parameters required by the Euler-Lagrange equation are obtained with a quantum circuit with a small number of quantum operations.QIPA generalizes the variational quantum approach to evolve an arbitrary initial density distribution so that it would become localized at the global minimum of a given cost function (see Fig. 1).As in IPA, the propagator of QIPA is not limited to the imaginary time quantum propagator (with infinite mass) enabling the use of other propagators that are maximal at the minimum of the cost function.
Fig. 1 shows the overall work flow of the QIPA algorithm.First, we select parameters θ 1 , θ 2 , . . ., θ N θ corresponding to the initial state.Second, we use a quantum co-processor and the Hadamard test to calculate the parameters A k,m , and C k introduced by the Euler-Lagrange equation m A k,m θm = C k .Next, we update the parameters θ j by numerical integration of the Euler-Lagrange equation using a classical computer.Having updated the parameters, the process is iterated until convergence to obtain parameters θ j corresponding to a distribution function localized at the global minimum.The converged parameters are then used by a quantum circuit to prepare the final state and obtain the expectation value of the global minimum coordinates.
The paper is organized, as follows.First, we introduce the quantum iterative power algorithm from a general perspective.Next, we apply QIPA to a wide-range of applications and we compare the performance of QIPA with that of QITE, including a molecular Hamiltonian example, the H 2 complex, quantum computer-aided design, and biprime factorization.Finally, the Discussion and Conclusions section includes a discussion of results and ideas for further development of the QIPA algorithm.The Tequila package [69] has been used to implement all the quantum circuit operations, using Qulacs [70] as a quantum backend.

II. ORACLE FUNCTION
The main idea behind our proposed algorithm is the realization that an arbitrary cooling function [71,72] or strictly positive oracle function that is maximized at the location of the global minimum of a potential energy surface V can be used to amplify the global minimum amplitude of any initial state |ψ(0) with finite amplitude at the global minimum.Here, we show that oracles defined by concatenated exponential functions, with n > 1, α 0 (y) = y and real constants a 1 , . . ., a n = 0, provide effective algorithms based on a generalization of the McLachlan's variational principle (Appendix A).The calculations reported in this paper are based on the oracle defined by the double exponential α 2 (−τ V ) = e e −τ V , which we obtain for the choice n = 2 and a 2 = a 1 = 1.
A particular case of global optimization involves the search of the ground state of a Hamiltonian H, a problem that is typically solved by imaginary time propagation.QIPA can solve that problem analogously by simply replacing V by H in the definition of the oracle, α n .In that case, the normalized oracle function f (H; τ ) acts on the initial wavefunction |ψ(0) (onwards, = 1), as follows: where U n (τ ) = α n (−Hτ ).In particular, α 1 corresponds to the standard imaginary time evolution, which is widely used in quantum Monte Carlo algorithms [73,74].

Algorithm 1 Variational quantum iterative power algorithm
Refs. [18,19] show that one can perform imaginary time evolution [22] with unitary gates defined by Eq. ( 2) with n = 1 that evolve the initial state according to the Wickrotated Schrödinger equation, obtained by the McLachlan's variational principle: where Here, we introduce a family of near-term quantum algorithms defined by α n with n > 1 that evolve the initial state according to the generalized Wick-like-rotated Schrödinger equation (Appendix A): where With the choice n = 2, we arrive at a double exponential function and the following Wick-like-rotated Schrödinger equation: with According to the McLachlan's variational principle, when we constrain the equation of motion as such the result is equivalent in finding a solution of the linear equation where the entries of the symmetric and positive semidefinite matrix A and the right-hand side C can be computed on a quantum computer by deploying the Hadamard test.The value of θ is updated with θ for a short time step δτ > 0 according to the Euler method as θ(τ + δτ ) ≈ θ(τ ) + θ(τ )δτ .The underlying assumption is that we can approximate |ψ(τ ) are parameterized quantum circuits (PQCs), with θ = (θ 1 , . . ., θ N θ ) the corresponding real-valued parameter vector.
In general, for an N -qubit system with Hamiltonian H with N H ≥ 1 Pauli words and a parameterized wavefunction |φ(θ) (where τ dependency θ(τ ) is understood throughout) with N θ ≥ 1 parameters, the upper bound for the number of distinct measurements N A required to obtain the matrix A for QIPA via the Hadamard test and the number of gates required are N θ (N θ − 1)/2 and G N A ≥ N θ , respectively.Such an estimate can be understood as the number of times required to completely evaluate all the A matrix elements since A is symmetric.Moreover, to obtain the vector C, the number of measurements and gates required (assuming a second-order Taylor series expansion of the required function of the Hamiltonian) are N θ and and G N C holds when two-qubit gates are not parameterized, while '=' sign holds when they are parameterized.Assuming a polynomial scaling: ), with the same number of Hadamard test measurements required.In general, QIPA yields improved convergence in shorter times compared to QITE, requiring the same number of Hadamard test operations and a higher number of unitary gates.

A. Molecular ground-state search
It would be natural to use a fully fault-tolerant general purpose quantum computer to simulate large-scale  quantum chemical molecular systems.However, due to the high overhead requirements of such a quantum simulation, a full-scale quantum chemistry simulation has not yet been seen.Henceforth, it is important to continue to push the frontier of quantum simulation with limited quantum resources.
The Hamiltonian for a chemical system in the second quantization picture has the following general form where a † i is a creation operator that creates an electron on the ith orbital, a i is an annihilation operator which removes an electron from the ith orbital, and h ij and V ijkl are the one-electron and two-electron interaction coefficients, respectively, which are determined for specific systems.The antisymmetric property of electrons is fulfilled by the anti-commutation relation of the creation and annihilation operators {a i , a † j } = δ ij , {a † i , a † j } = 0.The above anti-commutation relation precludes direct encoding of a chemical Hamiltonian on a quantum computer, since the operating units of a quantum computer (i.e., qubits) obey the commutation relation of spins.The remedy to this discrepancy is to perform a fermion-spin mapping, such as the Jordan-Wigner transformation [75] presented here as an example.The transformation maps the fermionic creation and annihilation operators to qubit raising and lowering operators σ ± = X ±iY with a string of Z operators to enforce the fermionic anti-commutation properties: With the above transformation, the fermionic anticommutation relation is preserved.For other fermionspin mapping approaches, the reader is referred to the literature [76][77][78].An alternative method to decompose gates for molecular systems is provided in ref. [79].
Here, we show that one can efficiently search for the ground-state energy of hydrogen molecule across various bond stretching distances with QIPA with fewer time steps than QITE [18][19][20][21][22][23][24][25][26].Results are shown in Fig. 2(a) comparing QIPA and QITE.For consistency and fair comparison, we use the same time step for each bond distance for both QIPA and QITE runs.The error difference between the exact energy obtained from full configuration interaction (CI) calculations and the QIPA and QITE results can be seen in Fig. 2(b).QIPA features less error for all bond distances considered.

B. Quantum computer-aided designs
As the number of high-quality qubits inside a quantum processing unit (QPU) grows over time, it is expected that eventually no classical supercomputer will be able to simulate, verify, and cross-check the inner working mechanism and data obtained from the QPU.This is commonly known as "quantum advantage."Once such an event occurs, from a practical point-of-view, it is beneficial to make use of existing quantum hardware that is already well-calibrated to simulate subsets of new QPU designs.Quantum computer-aided designs of superconducting qubits [80] and photonic chips [81] have recently been proposed and experimentally realized in a superconducting qubit architecture [82], but not yet with imaginary-time-like evolution.Here, we show that, with the proposed QIPA, we are able to optimize for the ground-state-energy search of a flux-tunable superconducting transmon system.Given an arbitrary classical electrical circuit diagram composed of inductors, capacitors, and Josephson junctions, one can quantize such circuit into a quantum Hamiltonian [83] via the Legendre transformation.Once we obtain the quantum Hamiltonian, the task is to translate it into a language that a quantum computer can understand, such as Pauli words or strings.Let us consider the case of a flux-tunable transmon system shown in the main text Figure 1 prior to conversion as an example.The system has the following Hamiltonian: Here, e is the electron charge.The normalized external flux f = Φ ext /Φ 0 is derived from the external magnetic flux Φ ext that penetrates the loop formed by the two Josephson junctions of the transmon.The Josephson energy of the two junctions is equal and given by E J while C is the capacitance.The magnetic flux quantum Φ 0 is a fundamental constant that describes the smallest amount of flux that a superconducting loop can sustain.Here, ϕ and N are the phase and number operators, respectively and fulfill the commutation relation [ϕ, N ] = i.And, the following relations follow: where |n j are the eigenstates of N .We notice that the operators e ±iϕj are similar to the usual bosonic creation and annihilation operators, without the square root prefactor.They are, in fact, the Susskind-Glogower phase operators [84].When we write down the transmon Hamiltonian in the charge number basis, we use (11) and assign the operators as: In general, the number of Cooper pairs can take on infinitely many integer values.However, for practical purposes, we are only interested in low-lying energy states.
In that case, we can truncate the Hilbert space as described above by introducing a finite maximum number of excitations d = 2 k .The number of data qubits used in the quantum simulation of the transmon qubit is k ∈ N.
Next, we convert the eigenstates of the number operator N into the computational basis states of k data qubits by representing the integer charge number in a preferred encoding [85][86][87].This implies a truncation of the physical space to the subspace spanned by 2 k Cooper pair numbers.There are combinatorially many ways to map such a state space to a set of qubits.For all the numerical quantum simulation experiments presented in this work, we have employed the Gray code [80,81,87] due to its resource-efficient representation of tridiagonal quantum matrix operators (see Table I).After integer labeling, each level is encoded into a set of bits, which is then mapped to Pauli operators, as follows: Here, X, Y, Z are the usual Pauli matrices and I is the identity.Encodings for the operators seen in (10) truncated at d = 16 are given in Table I.As shown in Fig. 2(c), both QITE and QIPA successfully minimize the energy of a given transmon circuit.Moreover, QIPA requires significantly fewer iterations than QITE for global optimization of the energy.These results constitute the successful implementation of two forms of imaginarytime-like evolution to perform quantum computer-aided design.
To solve the prime factorization problem for a given biprime (product of two prime numbers) N = q * × p * , we consider the Hamiltonian H N (q, p) = d(N ; q, p) 2 , d(N ; q, p) = N − q × p, (15) defined on the space of prime numbers q, p ≤ √ N .This Hamiltonian is a non-negative function that attains its global minimum H N (q * , p * ) = 0 for the unique pair of solutions q = q * and p = p * .Binary representation [24,80,81,90,92,95] of the prime factors yields where (q L , . . ., q 0 ), (p L , . . ., p 0 ) ∈ {0, 1} L are the binary representations of q and p, respectively, and L = log 2 (N/2) + 1.Assuming that p * , q * > 2 (otherwise N is an even number), we further restrict the search space by setting q 0 = p 0 = 1 such that ( 16) can also be written in terms of the combined parameter x = (q 1 , . . ., q L , p 1 , . . ., p L ), as follows: where x 2 l = x l holds for all 1 ≤ l ≤ 2L since x l ∈ {0, 1}.Equivalently, in terms of the scaled spin parameters s l = 2x l − 1, we write where s 2 l = 1 holds since s l ∈ {−1, 1}.For example, for N = 15 we obtain L = 2 and d(15; s) = 15 with twin global minima.Alternatively, with previous information about the prime factors of 15, the integer can be factorized with fewer qubits via the test Hamiltonian [24] which features a unique ground state |011 that corresponds to the correct factorization of the number 15 to 3 and 5.The resulting circuit implementation is given in Fig. 7.
As shown in Fig. 3(a) and (b), we find QIPA identifies the prime factors of biprimes with fewer iterations than QITE for all integers factored.The number of required iterations varies depending on the biprime factored and the ansatz, with speedups shown here of up to 50%.In addition, QIPA successfully identifies both factors of biprimes.As expected, QIPA succeeds for stable integration with a small, converged time step, for which the twin solutions are readily identified as the two components of the final wavefunction with the largest and equal amplitude.Moreover, QIPA succeeds for unstable integration with a larger time step, for which the twin solutions are identified as the two largest components of the final wavefunction with unequal amplitude, as depicted in Fig. 3(c).To our knowledge this ability to identify multiple prime factors is unique among published quantum computing prime factorization results.

IV. DISCUSSION AND CONCLUSIONS
In summary, we have presented a family of generalized imaginary-time-like near-term quantum algorithms which we coin the "quantum iterative power algorithm," inspired by its classical counterpart.We have analyzed its convergence rate.
One caveat is that since the proposed algorithm relies heavily on the ansatz circuit used, its convergence rate is difficult to discern in the generic case.We have also determined QIPA's estimated resource count, and demonstrated it can outperform the quantum imaginary time evolution while it reduces the number of required iterations, at the cost of a moderate increase in the number of gates.Furthermore, we have used the three numerical case studies -quantum chemical molecular simulation of the hydrogen molecule for various bond dissociation distances, quantum computer-aided design of a superconducting transmon, and finding optimal solutions for double prime factorization -to highlight how QIPA outperforms QITE.
We would like to point out that an additional consideration for such an algorithm besides the choice of the ansatz circuit is setting the right parameter for the time step δt at each evolution step.A number of proposals [96,97] suggest the use of an adaptive time step to overcome such an issue.However, there exists an opportunity to develop a systematic way to adjust the time step δτ according to gradient descent, rather than with a heuristic argument on how to adjust δτ .A major drawback with quantum imaginary-time-like evolution algorithm is that it involves constructing the matrix elements of A with NISQ hardware.Since we are working with noisy quantum hardware, any large fluctuation in the matrix elements would result in suboptimal θ angles.Future work will explore preconditioning the matrix based on errors associated with performing the Hadamard tests.
The generality of the global-optimization method presented here invites further application to other problems that currently have not been explored with QITE-based quantum computing algorithms.Quantum approaches could facilitate identification of reaction pathways and transition states in chemical physics, as well as optimization in a broad range of machine learning applications.The method also provides a general framework for adaptation of a class of classical optimization algorithms to quantum computers to further broaden the range of algorithms amenable to implementation on current NISQ quantum computers.
The main code used to generate data presented here can be found at the following repository https://github.com/aspuru-guzik-group/QIPA/.obtained from the McLachlan's variational principle.Note that the denominator is a normalization factor.Since a key to our proposal is the realization that there is nothing except the requirement that there be a continuous, integrable, strictly positive oracle U (x) that is maximized at the global minima of H(x) = V (x) to prevent us from assuming a particular form of the oracle; we present quantum imaginary time evolution more generally in terms of non-unitary time evolution according to an oracle, as follows: where α could be any strictly increasing positive function.When α(y) = e y , we recover imaginary time evolution, which corresponds directly to the Wick rotation (τ = −it).With that choice of α(y), we obtain that the above quantum state satisfies the Wick-like rotated Schrödinger equation: where E 1 (τ ) = ψ(τ )| H |ψ(τ ) .Even though |ψ(τ ) is a valid wavefunction that can be represented in a quantum computer, the non-unitary time evolution cannot be straightforwardly mapped to a quantum circuit based solely on unitary gates.Here, the McLachlan's variational principle [18,19,101,102] comes to the rescue and demands that with • representing the L 2 -norm and δ the functional derivative.In the following, we intend to simulate the action of non-unitary dynamics, (A3), on a quantum computer via McLachlan's variational principle.
In variational quantum simulations, instead of directly encoding the quantum state |ψ(τ ) at time τ , we approximate it with a parameterized quantum circuit |ψ(τ ) ≈ |φ(θ(τ )) with a real-valued parameter vector θ(τ ) = (θ 1 (τ ), θ 2 (τ ), . . ., θ N θ (τ )).We assume that physically relevant quantum states span a restricted region of the full Hilbert space [103] for a given time interval, such that the trial state parameterized by θ is sufficient to prepare a desired quantum state by applying a sequence of parameterized unitary gates U to the initial state | 0 = |0 • • • 0 .Thus, we have |φ(θ) = U (θ) | 0 , where U (θ) is referred to as the ansatz, and U k (θ k ) is the kth unitary gate controlled by classical parameter θ k .Here, we are only concerned with single-or two-qubit gates, which is sufficient for universal quantum computing.We see that the resulting equation for double-exponential QIPA takes a similar form to the original Wick-rotated equation of motion, but with more rapid convergence to the ground state of H. (A27)

Figure 1 .
Figure 1.Sketch of the quantum iterative power algorithm.First, the physical problem is mapped into the language of the quantum computer.Second, initialization of a parameterized wavefunction |φ(θ(τ )) is achieved by using a certain quantum ansatz circuit.Third, based on the ansatz choice, Hadamard test measurements are performed to obtain the A matrix and C vector on the quantum computer.Fourth, the new θ parameters are obtained from A and C on the classical computer.If the desired convergence is obtained, the program is stopped and the global minima are identified.Otherwise, the step of evaluating A and C is repeated to obtain new angles θ.

Figure 2 .
Figure 2. H2 energy dissociation curve in the minimal basis set (sto-3g).(a) Exact diagonalization result/ full CI (black solid line) is seen with data points from QIPA (red squares) and QITE (black dots) runs for different bond distances.(b) Absolute energy difference between the exact energy and the QIPA and QITE results in Hatrees, corresponding to the data in (a).(c) Ground-state energy optimization plot for a flux tunable transmon at the external flux f = 0.25 as a function of the number of iteration steps for both QIPA and QITE.In all results, both QIPA and QITE are run with the same time step.Here QIPA runs require significantly fewer steps to reach the convergence criteria.

Figure 3 .
Figure 3. QIPA factorization of biprimes 55, 65, 77, and 91 for the (top) YZ and (center) Y Ansatz as compared to QITE for equal time steps, and (bottom) the amplitude of the wavefunction components corresponding to the prime factors of 15 for varying numbers of time steps in QIPA.The black dashed line represents to the final amplitude of the wavefunction component with the largest magnitude at convergence.Here, all the horizontal dashed lines represent pre-defined numerical threshold where our algorithm would stop running.