Quantum algorithms for generator coordinate methods

This paper discusses quantum algorithms for the generator coordinate method (GCM) that can be used to benchmark molecular systems. The GCM formalism defined by exponential operators with exponents defined through generators of the Fermionic U(N) Lie algebra (Thouless theorem) offers a possibility of probing large sub-spaces using low-depth quantum circuits. In the present studies, we illustrate the performance of the quantum algorithm for constructing a discretized form of the Hill-Wheeler equation for ground and excited state energies. We also generalize the standard GCM formulation to multi-product extension that when collective paths are properly probed, can systematically introduce higher rank effects and provide elementary mechanisms for symmetry purification when generator states break the spatial or spin symmetries. The GCM quantum algorithms also can be viewed as an alternative to existing variational quantum eigensolvers, where multi-step classical optimization algorithms are replaced by a single-step procedure for solving the Hill-Wheeler eigenvalue problem.

This paper discusses quantum algorithms for the generator coordinate method (GCM) that can be used to benchmark molecular systems.The GCM formalism defined by exponential operators with exponents defined through generators of the Fermionic U(N) Lie algebra (Thouless theorem) offers a possibility of probing large sub-spaces using low-depth quantum circuits.In the present studies, we illustrate the performance of the quantum algorithm for constructing a discretized form of the Hill-Wheeler equation for ground and excited state energies.We also generalize the standard GCM formulation to multi-product extension that when collective paths are properly probed, can systematically introduce higher rank effects and provide elementary mechanisms for symmetry purification when generator states break the spatial or spin symmetries.The GCM quantum algorithms also can be viewed as an alternative to existing variational quantum eigensolvers, where multi-step classical optimization algorithms are replaced by a single-step procedure for solving the Hill-Wheeler eigenvalue problem.

I. INTRODUCTION
The rapid development of quantum technologies and quantum algorithms addresses long-standing computational challenges of many-body physics and quantum chemistry.While the primary target is to overcome exponential growth in complexity associated with approaching the exact limit in the simulations, the possibility of identifying physically meaningful solutions to problems of interest is equally important.Although Quantum Phase Estimation algorithms [1][2][3][4][5][6][7] have been designed in a way that adequately addresses both of these issues, their applicability is currently limited by the necessity of using complex quantum circuits with corresponding depths that preclude its practical applications on existing quantum computing platforms, dominated by Noisy Intermediate-Scale Quantum (NISQ) devices.Instead, hybrid algorithms such as various Variational Quantum Eigensolver (VQE) [8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25][26] are currently being intensively tested on NISQ quantum computers to characterize the properties of correlated quantum systems.In this effort, in many aspects, the VQE formulations, for example, based on the unitary coupled-cluster (UCC) [27][28][29][30][31] representation of the wave function, mirror the standard conventional formulations of CC theory, where a large number of excitations are included in the cluster operator and simultaneously optimized in the iterative process.This algorithm leads to another set of challenges associated with the potential problems with the convergence of iterative process (commonly referred to as the barren minimum problem) and representation of UCC Ansatz on quantum registers, which may result in long Trotter-like products of exponential operators defined by multi-qubit gates.Although several strategies mitigating these problems have recently been proposed, utilizing VQE-UCC formulations and extending these methods beyond specific system-size limits may be challenging and require reformulation of the quantum many-problem into recently introduced quantum flow equations, 32 where large subspaces of Hilbert space can be sampled through the constant-depth small-dimensionality coupled eigenproblems.
a) Electronic mail: karol.kowalski@pnnl.govInstead, in this paper, we explore the applicability of the Generator Coordinate Method (GCM) [33][34][35][36][37][38][39] as an alternative to popular VQE formulations.The main difference with the VQE method is that the GCM method avoids highlynonlinear parametrization of the wave function and provides an efficient mean for direct extension of the probed subspaces.Additionally, it offers an efficient utilization of Ansatzes represented by low-depth quantum circuits.
The GCM method was one of the first attempts to combine two distinct aspects of many-body theories: independent-particle models and theories describing collective phenomena, where the approximate eigenstates |Ψ GCM of the Hamiltonian H are expressed using a family of N-body wave functions |Φ(q) , where q is a set of collective variables that describe correlation effects in many-body systems, and usually, the corresponding |Φ(q) is represented a complicated linear combination of Slater determinants.The scalar f (q) is referred to as the weight function.The advantage of the GCM approach is obtaining ground states and classes of excited states described by the chosen set of generator coordinates.In general, in the analogy to the coherent state representation, 40,41 the family |Φ(q) forms an overcomplete basis.Upon substituting (1) into Schödinger equation one gets the so-called Hill-Wheeler integral equations for unknown f (q) coefficients dq [H(q, q )-ES(q, q )] f (q ) = 0 , where the integral kernels H and S are defined as In typical applications, the Hill-Wheeler equation is usually solved numerically by discretization, which transforms integral equation ( 3) into an algebraic eigenvalue problem.There are two categories of GCM formulations in applications to many-body quantum systems.The first category's purpose is to restore broken symmetries of the |Φ(q) states.For example, Bardeen-Cooper-Schrieffer states are not eigenstates of the particle number operator N. To project the Bardeen-Cooper-Schrieffer states onto the Hilbert space with the desired number of particles N 0 , one uses the projection operator expressed in terms of integral over gauge angle φ N 0 , In this case, the properties of f k (q) coefficients are determined by the properties of the symmetry (projection) operators (e.g., particle number projection operator).In the second category of GCM formulations, the unknown weight function is optimized to capture correlation effects encoded in generator coordinates.However, designing an adequate grid or path to probe the Hill-Wheeler equation is a rather empirical procedure that requires much intuition and prior knowledge of the sought-after many-body system.The appealing feature of the GCM method, especially from the point of possible quantum computing applications, is the possibility of combining low-depth representations of |Φ(q) functions with simple, one-step optimization conditions for weight function.In this approach, the role of quantum computing is to map a discrete number of states {|Φ(q p ) } M p=1 to quantum register and to evaluate matrix elements for Hamiltonian and overlap matrices (Eqs.( 3) and ( 4)).In contrast, solving a generalized eigenvalue problem in a discrete basis representation takes place only once on a classical computer, avoiding multiple instances of quantum-classical communication as in the VQE formalism.In the discussed formalism, we follow an "algebraic" GCM formulation discussed by Fukutome in Ref. 42, extend the standard GCM to multi-product exponential formulas, and provide an algorithm for sampling coordinate space in a way that provides selective approaching classes of excited Slater determinants.In this context, the multi-product GCM formulation alleviates some problems associated with the usage of high-fidelity of Trotter-type expansions.
We demonstrate the performance of the quantum GCM algorithm on the example of the H4 benchmark system in various configurations.We show that by a judicious choice of the |Φ(q) , one can recover high-level of accuracy both in weakly and strongly correlated regimes using |Φ(q) s that use the manifold of single excitations (or U(N) Lie algebra generators).The obtained level of accuracy is similar to the one obtained with the advanced VQE formulations.

II. THEORY
The GCM was initially introduced to describe collective effects in nuclei. 33,36,39,43Fukutome in his seminal paper 42 considered the GCM from the Lie-algebra theoretical standpoint.We use that language throughout this paper.Let us assume that the Fermion system is described by annihilation and creation operators a p and a † q that satisfy the following set of anticommutation relations: As discussed in Ref. 42 there are several Fermion algebras including U(N), SO(2N), and SO(2N + 1) Lie algebras, and Clifford algebras that can be used to characterized approximate many-body wave functions (here, N stands for the number of single particle states).These algebras can be defined by the following set of operators For example, • U(N): • SO(2N): {E p q , E pq , E pq } , • SO(2N + 1): {a p , a † q , E p q , E pq , E pq } .
The U(N) algebra is the only algebra where particle number operator commutes with all operators belonging to U(N) algebra.
Let us focus attention on the U(N) algebra.By Γ(Z) we designate an anti-Hermitian operator defined as where Z can be viewed as an anti-Hermitian matrix [z pq ] where z qp = −z pq or set of indexed parameter {z pq } satisfying z qp = −z pq .Canonical transformation U(Z) generated by Γ(Z) takes the form Standard canonical Thouless transformation of the Hartree-Fock (HF) determinant |Φ can be obtained by acting with U(Z) onto |Φ , i.e., The standard unitary CC model with singles (UCCS) is a special case of (10) where z i j = z ab = 0, where i, j, . . .and a, b, . . .correspond to spin-orbital indices occupied and unoccupied in the Slater determinant |Φ (where Z can be identified with the set of elements {z ia }).In comparison to |Φ , |Φ(Z) can contain certain classes of correlation effects (that are Z dependent; it can contain elements of symmetry breaking and effects responsible for the appearance of HF instabilities).The parametrized states (10) can be viewed as a nonorthogonal (in general overcomplete) basis in the Hilbert space.The generator coordinate method utilizes this fact by representing the wave function in the form where collective variables q of Eq.( 1) are now identified with Z matrix, i.e., q → Z .
A typical way for solving Hill-Wheeler equations ( 2) is through discretization of the Z domain.Let us introduce a set of the Z-points -Q = {Z i } M i=1 then equations (2) takes the form of non-orthogonal eigenvalue problem where H and S (M × M matrices) and f (M-dimensional vector) are defined as follows: Using this form of discretization, the optimal form of wave function (we assume the ground-state wave function in this paper) is given by the expansion which is reminiscent of recently discussed non-orthogonal variational approaches discussed in the context of quantum computing. 15If M is equal for a given basis set to a dimension of the full configuration interaction (FCI) problem (for a given spin and spatial symmetry) and e Γ(Z p ) |Φ are linearly independent, then the expansion (17) with optimized f i coefficients describes the exact electronic wave function.

III. MULTI-PRODUCT EXTENSION OF THE GENERATOR COORDINATE FORMALISM
A possible extension of the GCM expansion given by Eq.( 11) can be provided by the expansion involving multiple products (the product GCM formalism, abbreviated as PGCM (k) ) of k exponential operator (which is inspired by a recent progress in the development of dynamical GCM methods 39,[43][44][45][46].For example, one can introduce the following expansion where In the above representation, Γ(Z(i)) can belong to various Lie algebras.In the following we will focus on the case where all Γ(Z(i))'s (i = 1, . . ., k) belong to the same U(N) Lie algebra; that is, In fact, the PGCM formula may be viewed as a special case of the GCM where We will demonstrate that the U(N) case of PGCM lends itself for an efficient way of representing higher-rank excitations in quantum computing.In particular, this goal can be achieved by using a simple algorithm for discretization of Z(i) domains in the GCM method.In particular, we can show that PGCM (k) can be used to approximate 2k-tuple excitations in the configuration-interaction-type expansion.
Let us focus on the specific case when k = 2 (i.e., the PGCM (2) formalism).In this case we will represent the |Ψ PGCM wave function, given by the formulas: and The GCM algorithm can be adapted easily for the product representation of the trial wave functions.Now the discretization of the problem, in analogy to Eq.( 13), involves a Q set defined as Q = {Z I } I , where where we use composite index I = (p, q).Controllable sampling algorithms of the sub-spaces of the Hilbert space corresponding to higher-rank excitations with the PGCM formalism based on the U(N) algebras require careful selection of the sampling points, which will be discussed in the following section in the context of PGCM (2) method applications to the H4 system.

IV. H4 MODEL: THE CHOICE OF THE GCM SAMPLING POINTS
We use the H4 model of Ref.48 as a benchmark system.The system's geometry can be defined by a single parameter α.For α = 0.5, the H4 model corresponds to a linear chain of hydrogen atoms with distances between adjacent hydrogen atoms equal to 2.0 a.u.; for α = 0.005, the system is almost in a square configuration.The main difference between α = 0.5 and α = 0.005 H4 models is in the structure of the corresponding ground-state wave function.While for α = 0.5, the ground-state wave function is dominated by the restricted Hartree-Fock (RHF) determinant, for α = 0.005, the ground state is quasi-degenerate.The contribution of the RHF Slater determinant is almost the same as the doubly excited configuration where two electrons from the highest occupied orbital are promoted to the lowest lying virtual orbital.The spin-orbital numbering scheme is shown in Fig. 1.In the PGCM (2) formalism for H4, we will adopt the following sampling scheme, where the choice of the cluster operators is consistent with analyzing HF equation stability conditions based on the Thouless theorem [49][50][51][52] :

|R
(±) These sampling vectors can be naturally tied to the general form of the |Φ (2) (Z(1), Z(2)) basis given by Eq.( 23).For example, |R 0 corresponds to 1 (t 1 ) , Z(2) = 0, and z(1) 52 = z(1) 74 = −z(1) 25 = −z(1) 47 = t 1 while remaining matrix elements are equal to zero, etc.Another advantage of using this form of sampling vectors is that their combinations provide a rudimentary (yet not exact) mechanism for eliminating symmetry impurities when R i operators break the symmetry of the reference state |Φ .For example, if in general R (±) 1 (t 1 ) the operator is expressed in terms of excitations that produce a triplet state when acting on the reference function, then this singlyexcited impurity (or instability) is eliminated by taking a combination of where linear triplet "impurities" are eliminated.These combinations also can allow us to selectively approach doubly excited configurations, or using combinations of 1 (t 7 ) , quadruply excited ones using manifold of single excitations and very simple quantum circuits to represent e ±tR i operators.This analysis can be extended to higher-order PGCM (k) formulations to include higher-rank excitations.It should be stressed that t i parameters in Eqs.( 24)-( 31) can be chosen at random.The effects of random t i parameters on ground-state FCI energies are illustrated in Appendix A.

V. GENERAL OUTLINE OF QUANTUM ALGORITHM AND POST-PROCESSING ON CONVENTIONAL COMPUTERS
After selecting GCM sampling points, the remaining work is involves computing matrices H and S through Eq.( 14) and ( 15) using quantum computers.To accommodate gate-based quantum computers, the expectations are computed in the form of and where Hamiltonian matrix H is transformed to the linear combination of Pauli strings H = ∑ j h j P j under Jordan-Wigner (JW) transformation and |Φ is the HF state.The expectations can be easily evaluated using algorithms like the Hadamard test on fault-tolerant quantum computers.A pure quantum algorithm proposed for GCM and its complexity analysis is given in Appendix B, where we propose to exploit a block encoding for the operation S −1 and use the phase estimation to compute the eigenvalues for the nonorthogonal eigenproblem (13).Note that the pure quantum algorithm gives a favorable scaling but relies on the use of the controlled unitary circuits that makes the quantum simulations on NISQ devices challenging.
For near-term devices, we focus on hybrid quantumclassical approach.there is a trade-off between the depth of the quantum circuit and classical computation time.Note that, following the JW transformation, each of the operators R 1 to R 4 in Eq.( 25) to (28) can be transformed to four commuting Pauli strings.That is to say, the matrix exponential e ±t i R i for i ∈ {1, 2, 3, 4} are exactly the linear combinations of Pauli strings.So, for both H4 models, the calculations of expectation in Eq.( 33) and (34) are essentially equivalent to Φ|P|Φ for some Pauli string P, while qubit-wise commuting (QWC) terms are grouped to perform a simultaneous measurement.The only implemented form of the circuit is illustrated in Fig. 2 and the whole process is summarized in Alg.1.The entire simulated quantum computation is achieved using Qiskit. 53We put a brief discussion about trotterization under Qiskit in Appendix C, along with other implementation details related to duplicated Pauli strings and effects of random parameters in Appendix A.
Algorithm 1 Quantum GCM (QuGCM) for near-term devices Require: Hamiltonian matrix H = ∑ j h j P j , HF state |Φ , and a set {Γ(Z i )} M i=1 where the index i could be a composite up to k terms as in Eq. (23)   1: Transform all {Γ(Z i )} M i=1 using JW transformation 2: Generate unitaries {V i } M i=1 for V i := e Γ(Z i ) with Eq.( 9) and Eq.(20)   3: Trotterize each element in {V i } M i=1 to a linear combination of Pauli strings 4: for each V q in {V i } M i=1 do

5:
Compute ∑ j h j P j V q classically 6:

7:
Compute ∑ j h j V † p (P j V q ) classically 8: Compute V † p V q classically 9: Evaluate H pq := ∑ j h j Φ|V † p P j V q |Φ in a quantum device 10: Evaluate S pq := Φ|V † p V q |Φ in a quantum device

VI. COMPLEXITY ANALYSIS AND EFFECTS OF FINITE SAMPLINGS
Recall that M is the number of selected Z-points, which is equivalent to the number of sampling vectors.Let n be the number of spin orbitals.We know the Hamiltonian matrix, H, consist of at most O(n 4 ) Pauli strings, so operations for multiplying matrices and solving the generalized eigenvalue problem in Alg.As we have M 2 iterations in Alg.1 and the final generalized eigenvalue problem of M × M matrices takes O(M 3 ) operations, Hybrid PGCM (2) has the overall classical part of the worst-case scaling of O(c 2 n 5 M 2 + M 3 ) operations.
The specific value of c highly depends on the how users want to approximate the molecular models.It is affected by the number of creation and annihilation operators pairs in each single-excitation cluster operators, the value of k, and number of trotterization steps if it is necessary.As we discussed earlier, we can properly choose two pairs of creation and annihilation operators in each of single-excitation cluster operators without the requirements of trotterization.The results in Section VII are promising enough when only single and double excitations are considered, i.e., when k = 2.
Regarding the quantum part of the operation complexity, it includes the number of circuits and the number of measurements in every circuit.Naively, line 9 and line 10 in Alg.1 indicate there are O(n • c 2 • n 4 ) number of circuits and every circuit contains only O(n) number of 1-qubit gates.Many existing methods can be applied to reduce the number of terms associated with V † p (HV q ) and V † p V q .For example, to reduce the order O(n for the number of circuits in every iteration, the linear combination of unitaries technique, 54 amplitude amplification approach, 55 Hamiltonian simulation, 56 qubitization, 57 or the direct block-encoding 58 methods can be typically employed at the cost of introducing deeper circuits and implementing controlled unitary operations.Nevertheless, these approaches come with a probability of failure and require advanced circuit and error mitigation that might go beyond the capability of the current NISQ devices.Toward a more feasible NISQ approach, in light of the unitary partitioning scheme proposed by Izmaylov et al. 59 , Peng and Kowalski recently proposed a more efficient unitary partitioning approach guided by the single-reference trial state used in the simulation. 60In particular, through numerical tests over a wide range of molecules in different bases, they found that the non-unitary wave operators when acting on single-reference trial wave function (such as e Γ(Z i ) |Φ and He Γ(Z i ) |Φ in the present discussion) can be represented by a much more compact unitary basis, thus providing a more efficient route for performing the general non-unitary quantum simulations.
It is worth mentioning that the above discussion is focused on the number of terms/operations that can be efficiently reduced through groupings that feature the commutativity or anti-commutativity of Pauli strings, while the total number of measurements required from the number of groups also critically depends on the covariance between the contributing terms, cov P i , P j , and the desired precision ε.Given that we can always write a matrix operator as a linear combination of Pauli strings, the total number of measurements for evaluating the expectation value of an operator with respect to the trial wave function can then be expressed as [61][62][63] # of Measurements = ∑ G ∑ i, j,∈G h i h j cov P i , P j ε 2 (35)   with G indexing the groups.Therefore, it is likely that the number of groups decreases at the cost of introducing larger covariances that could essentially increase the total number of measurements required to achieve a desired precision.The variance of an individual Pauli string P i can be bounded by which provides bounds to the covariance of any contributing terms cov(P i , P j ) ≤ var(P i )var(P j ) ≤ 1.
Thus, according to Eq.( 35), the bounds for the standard deviations of the matrix entry H pq and S pq under a finite number of measurements are Let Hpq and Spq be the entries computed from finite number of measurements.We can empirically estimate the effects of the uncertainty by considering Hpq and Spq as normal random variables To illustrate the maximum possible fluctuations brought by finite-sampling errors, we set each ε H pq and ε S pq to their maximum as in Eq.( 38) and Eq.( 39), and summarize the results in Fig. 3.It is clear that to generally reach chemical accuracy, the α = 0.005 case requires about two orders of magnitude more measurements than the α = 0.500 case using Alg.1.
Figure 3. Differences between the ground-state energy estimations from PGCM (2) and FCI formalism in milli-Hartrees for α = 0.005 (up) and α = 0.500 (down) H4 models, respectively.The variance of each random sampling is set to the maximum according to Eq.( 38) and Eq.( 39).The red-shaded region is the range of chemical accuracy ±1.5936mHa.The blue-shaded region is the 95% confidence interval estimated from 100 simulations with the red bars marking the rough number of measurement per entry at which the energy difference would be within the chemical accuracy.
We also notice that there are also many other advanced measurement schemes proposed recently.One example is to simultaneously obtain expectation values of multiple observations by randomly measuring and projecting the quantum state into classical shadows [64][65][66][67][68][69][70][71][72][73] .In principle, the algorithm enables measurements of m low-weight observations using only O(log 2 m) samples.The practical performance of the algorithm for model and molecular Hamiltonians on NISQ devices, in terms of accuracy and efficiency, is still under intense study.

VII. RESULTS
The GCM results for the ground-state energies and excitation energies corresponding to low-lying states of the symmetry of the reference function (singlet A 1 states) are collated in Tables I and II, respectively.To evaluate the accuracy of ground-state simulations, we compared GCM results with the results obtained with the RHF, multiconfigurational self-consistent field (MCSCF) formalism for active space defined by four electrons and three active orbitals (MCSCF(4e,3o)), configuration interaction method with singles and doubles (CISD), CC method with singles and doubles (CCSD), VQE formalism, and FCI formalism.We used the equation-of-motion CC approach (EOMCCSD) and FCI formalism for excited states.Note that any quantum computation in VQE and GCM assume infinite number of measurements.
Inspection of the results in Table I indicates that for both geometries, the GCM results are in very good agreement with the FCI result despite the simplicity of the expansions given by Eqs.( 24)- (31).Interestingly, the GCM energies are significantly better in both cases than the VQE ones.The excellent performance of the GCM formalism is well illustrated by the strongly correlated α = 0.005 case, where CISD and CCSD methods struggle to capture needed correlation effects.For the weakly correlated variant of H4 (α = 0.500), the GCM formalism reproduces nearly FCIlevel accuracy with the energy error of 0.022 milli-Hartree.
In Table II we juxtaposed the excitation energies (ω 1 , ω 2 , and ω 3 ) obtained with the GCM approaches for three lowest-lying 1 A 1 symmetry states, of H4 model for α = 0.005 and α = 0.500.For the ω 1 excitation energies, the GCM approach provides consistently better estimates of their exact (FCI) values than the ubiquitous EOMCCSD approach.While for ω 2 GCM prediction is better than the EOMCCSD one only for the α = 0.005, for non-degenerate case (α = 0.500), the GCM prediction is by 0.6 eV off the FCI value.For the ω 3 excitation energies, the GCM is capable of providing estimates within 0.160 eV (α = 0.005) and 0.33 eV (α = 0.500) of error.Again, given the simplicity of the GCM formulations, one should view the GCM estimates of the excitation energies as quite satisfactory.

VIII. CONCLUSION
In this work, we explored use of the GCM in the context of quantum computing.For this purpose, we introduced the multi-product extension of the GCM formalism that enables one to construct state vectors in Hilbert space using various types of the Fermion Lie algebra and general quantum algorithms that allow one to perform GCM calculations in a way that can be viewed as a specific case of the quantum algorithms for configuration interaction formalisms involving non-orthogonal basis.In the present studies, we focused entirely on the U(N) algebra, where resulting state vectors can be interpreted in terms of the Thouless theorem.This analogy is essential in the sampling process of the parametrized unitary canonical transformations.It enables one to construct corresponding state vectors and corresponding linear space where higher-order excitations (e.g., double, triple, quadruple, etc.) excitations can be selectively approached using the language of single excitations to define generators (Γ(Z)) of the canonical transformations.The discussed procedures can be easily related to the searches of various type instabilities in independent particle formulations with the HF method as a specific example.
Using the H4 system as a benchmark, we showed that the quantum GCM algorithm could provide ground-state energies competitive to the VQE simulations involving the explicit form of the double excitations.In contrast to standard VQE algorithms, the GCM formalism also yields the values of excited-state energies.We showed that GCM excitation energies corresponding to the low-lying excited states could be competitive with the excitation energies obtained with the popular EOMCCSD approach.
Although in the present form, the GCM formalism falls into the category of hybrid formulations, unlike the VQE method, it avoids multiple quantum-classical machine data transfer and unstable and tedious iterative processes.Also, our near-term implementation of the quantum GCM algorithm can run in parallel easily using multiple quantum and classical nodes.
In future studies, we plan to use the full potential of the quantum algorithms based on the quantum CI formulations and explore the possibility of using other types of Fermionic algebras, i.e., SO(2N) and SO(2N + 1) Lie algebras.

X. DATA AVAILABILITY
The code and data in the paper are openly available in the GitHub repository. 74pendix A: Duplicated Pauli strings and selection of parameters For each of H4 models, while M = 15, the symmetries of H and S allows us to only compute 120 iterations instead of 225 iterations.Among those 120 iterations, as shown in Table III, we only need to measure 4,000 Pauli strings because more than 97% of them are duplicated.This brings the two orders of magnitude reduction on the total number of measurements in practice.
Meanwhile, we conducted the following experiments for α = 0.005 and α = 0.500 H4 models to demonstrate the influence of the random choices of t i parameters on the estimations of the ground-state FCI energies.For each molecular model, we generated 50 sets of {t i } 7 i=1 for t i in [0, 1), [0, 100), and [100, 100), respectively.Each of 300 sets of {t i } 7 i=1 produced an estimation of the ground-state FCI energy from Alg. 1 (assume infinite number of shots).The distributions of differences between the estimations from PGCM (2) and FCI formalism under various settings are illustrated in Fig. 4  Perhaps the most direct way to use quantum computing to solve the electronic structure problem using Generator Coordinate Methods is by simply solving the non-orthogonal eigenvalue problem by dilating it to a square matrix in a higher dimensional space.In this case, it is most convenient to express our eigenvalue problem as where we have assumed here that S is an invertible matrix.Now let us define an isometric extension of our original space.We do this to exploit a block encoding for the operation S −1 , which allows us to express it as a unitary operation in a higher dimensional space.Specifically, let be a unitary matrix (i.e., U † = U −1 ) for arbitrary matrices .Also let Inside this enlarged space, the eigenvalue equation reads for eigenvalue E for a set of unitary U (H) j .in various random {t i }'s and FCI formalism in milli-Hartrees for α = 0.005 (up) and α = 0.500 (down) H4 models, respectively.The horizontal red-shaded region shows the range of chemical accuracy ±1.5936 mHa.The fact that minimum value and the 25 th percentile are in the shaded region when t i ∈ [0, 1) in both H4 models indicates that random generations of t i in that range is relatively likely to provide a good approximation in our algorithm.
For us to use phase estimation to compute these eigenvalues, we also need a further modification.Note that e −iUJ is not necessarily unitary because UJ is not necessarily Hermitian.We address this issue by considering a further embedding: This expression also takes the form of a linear combination of unitary matrices; however, the coefficient sum now obeys ∑ j |h j | = 2α, which does not change the overall complexity.Finally note that K is anti-diagonal block matrix.We can therefore search for an eigenvector of the form h = 0 g † as the eigenvectors of K can be taken to have support only on one of the two blocks.The eigenvalues of K corresponding the the eigenvector g of PUJP must be ±E from this construction.To see this note that Thus h is an eigenvector of K and the eigenvalues of K in the support of h must be ±E.Thus, we can estimate the absolute value of E by using phase estimation on K.
To perform this simulation using qubitization ideas, we will need to propose prepare-and-select circuits for the coefficients.Specifically, B9) and further: This then forms a block encoding of our operator for the generalized eigenvalue problem: Using qubitization, we can convert this into a unitary with result 58 , we can perform phase estimation to learn E with a variance of at most ε 2 using a number of queries to U sel and U prep that is in O α S −1 /ε .However, each application of U sel requires O(1) queries to U, which is a block encoding of the inverse of S. A deeper analysis of the cost is possible; however, to do this, we first need to have consensus on the input model used for the simulation.
The easiest way to compute a matrix element for the overlap matrix is through the use of a controlled unitary.The overlap that, the estimate takes the form: where V p is the basis transform operation such that for the reference state |Φ , V p |Φ = |Φ(Z p ) .To use the Hadamard test to reconstruct this circuit we require a single application of V † p and a single application of V q .The probability that the control qubit that governs this circuit is 0 is:: If needed, the imaginary part can be similarly found by applying an S † gate to the control.If we apply amplitude amplification to the result, then we can construct a matrix AA(P) with eigenvalues inside the sector: and is equivalent to the following matrix in the two dimensional space spanned by the initial state and the marked state: Applying quantum signal processing we can apply a transformation u → 2u 2 − 1.Note that this is 1) an even degree polynomial and 2 it is between [−1, 1] for u in a similar range.This means that quantum signal processing can be used to apply this transformation on the top block of the matrix to prepare a block encoding of the form 58 : Where specifically we have that, the initial state used in the amplitude amplification routine block-encodes the overlap matrix S. In all, this process costs O(polylog(1/ε)) queries to the state preparation oracle to produce this block encoding.
Next we need to prepare a block encoding of S −1 .Using the results of Ref. 75 we can prepare such a block encoding using O( S S −1 polylog(1/ε)) queries where ε is our target accuracy.Thus, the overall query complexity is the number of queries made to K multiplied by the number of queries per U or U (H) j .The former result then shows that the final cost of the simulation using this approach (in terms of queries to the prepare and select oracles of H and the queries to Z p , Z q scale as: This shows that the query complexity of the simulation using this approach does not necessarily scale with the dimension of the space.It does, however, depend strongly on the one-norm of the coefficients of the Hamiltonian and the norm of the inverse of the overlap matrix.Thus, the worse conditioned the matrix is, the worse we expect the performance of the algorithm to be.In contrast, learning the matrix S in high dimensional spaces to perform the inverse can be expensive as noted in the main body.This approach gives a theoretical alternative in such cases that has favorable scaling asymptotically at the price of the algorithm requiring a large number of qubits.

Appendix C: Trotterization of a matrix exponential
This section aims to transform a matrix exponential into a linear combination of Pauli group matrices.Because the standard form of a matrix exponential in Qiskit is e −itA instead of e tA as we have in the main body for some parameter t and matrix A and we implemented the quantum part of our algorithm in Qiskit, we keep i = √ −1 in the power.Because we will use dot product between matrices, we do not omit the operator ⊗ when we do the Kronecker product between Pauli matrices for clarity (e.g., XY means the dot product of X and Y instead of X ⊗Y ).Let n be the number of qubits, I n the identity matrix in space C 2 n × C 2 n , and P ∈ {X,Y, Z, I} ⊗n an n-qubit Pauli group matrix.Then, for any n-qubit Pauli group matrix, we have: for some matrices A, B,C, D in the appropriate dimensions.Now, by doing Taylor expansion of e itP at t = 0, we obtain the following exact conversion:

Figure 2 .
Figure 2. Shallow circuit for Alg.1.State-preparation gates initialize the ground state to the HF state.Pauli gates are a single layer of 1-qubit Pauli gates.QWC rotation gates consist of 1-qubit rotation gates that allow us to measure each wire on Pauli-X,Y, or Z bases with the final Pauli-Z-basis measurement.So all gates in the circuit are 1-qubit gates and the total number of 1-qubit gates is O(n).
1 dominate the classical part of the operation complexity.If an exponent of a cluster operator has up to c number of Pauli strings after trotterization and the multiplication between two length-n Pauli strings takes at most O(n) operations, then the matrix multiplication in line 9 and 10 of Alg.1 takes O(n • c 2 • n 4 ) and O(n • c 2 ) operations, respectively, because all matrices are decomposed into the linear combinations of Pauli terms.
IX. ACKNOWLEDGMENTS Muqing Zheng, Bo Peng, Nathan Wiebe, Ang Li, and Karol Kowalski were supported by the "Embedding QC into Many-body Frameworks for Strongly Correlated Molecular and Materials Systems" project, which is funded by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, the Division of Chemical Sciences, Geosciences, and Biosciences.The Pacific Northwest National Laboratory is operated by Battelle for the U.S. Department of Energy under Contract DE-AC05-76RL01830.Xiu Yang was supported by National Science Foundation CAREER DMS-2143915.Muqing Zheng and Xiu Yang both also were supported by Defense Advanced Research Projects Agency as part of the project W911NF2010022: The Quantum Computing Revolution and Optimization: Challenges and Opportunities.
B5)Thus the entire product can be expressed asPUJP = ∑ j h j 2 P U I ⊗ U (H) j + U Z ⊗ U (H) j P,(B6) Thus the enlarged Hamiltonian J can be expressed as a linear combination of unitary matrices that further has the exact same value of α = ∑ j |h j |.

Figure 4 .
Figure 4. Differences between the ground-state energy estimations from PGCM(2) in various random {t i }'s and FCI formalism in milli-Hartrees for α = 0.005 (up) and α = 0.500 (down) H4 models, respectively.The horizontal red-shaded region shows the range of chemical accuracy ±1.5936 mHa.The fact that minimum value and the 25 th percentile are in the shaded region when t i ∈ [0, 1) in both H4 models indicates that random generations of t i in that range is relatively likely to provide a good approximation in our algorithm.

e+
−i ∑ m l=1 s l P l = m ∏ l=1 e −is l P l + O(m 2 s 2 ).(C4)where s := max l s l .So for large t, to control error, we need to separate the evolution into multiple steps:e −i ∑ m l=1 s l P l = m ∏ k=1 e −i(s l /r)P l r + O(m 2 s 2 /r).(C5)where r is the number of evolution steps.There is also a second-order formula for smaller errore −i ∑ m l=1 s l P l = O(m 3 s 3 /r 2 ).(C6)So, if we set the trotterization error level at ε, then we need to split the evolution into r ∈ O m 3/2 s 3/2 √ ε steps with the second-order formula.By choosing t = − s 2r , Eq.(C3) gives: , Eq.(C6) becomes: e −i ∑ m l=1 s l P l ≈ m this case, we approximate the exponential of the linear combination of the Pauli group matrices by another linear combination of Pauli group matrices, where the expectation of the latter one can be evaluated in a gate-based quantum computer easily.Note that in the r = 1 case, if we have m terms in the power, we end up with at most O(m 2 ) terms after the transformation, which is exactly the scenario we had in H4 examples in the main text.

Table III .
using box-plots.It worth noting that when we produce t i 's from [0, 1), about 25% random generations can provide an estimation of ground-state energy within the range of chemical accuracy.Number of unique Pauli terms and their ratios after grouping among all iterations.
is a scalar parameter.If we deal with more than a single Pauli group matrix, then we need to use Suzuki trotterization.For P l ∈ {X,Y, Z, I} ⊗n , the first-order Suzuki trotterization is: