Simulating Many-Body Systems with a Projective Quantum Eigensolver

We present a new hybrid quantum-classical algorithm for optimizing unitary coupled-cluster (UCC) wave functions deemed the projective quantum eigensolver (PQE), amenable to near-term noisy quantum hardware. Contrary to variational quantum algorithms, PQE optimizes a trial state using residuals (projections of the Schr\"{o}dinger equation) rather than energy gradients. We show that the residuals may be evaluated by simply measuring two energy expectation values per element. We also introduce a selected variant of PQE (SPQE) that uses an adaptive ansatz built from arbitrary-order particle-hole operators, offering an alternative to gradient-based selection procedures. PQE and SPQE are tested on a set of molecular systems covering both the weak and strong correlation regimes, including hydrogen clusters with 4-10 atoms and the BeH$_2$ molecule. When employing a fixed ansatz, we find that PQE can converge disentangled (factorized) UCC wave functions to essentially identical energies as variational optimization while requiring fewer computational resources. A comparison of SPQE and adaptive variational quantum algorithms shows that - for ans\"{a}tze containing the same number of parameters - the two methods yield results of comparable accuracy. Finally, we show that SPQE performs similar to, and in some cases, better than selected configuration interaction and the density matrix renormalization group on 1-3 dimensional strongly correlated H$_{10}$ systems in terms of energy accuracy for a given number of variational parameters.


I. INTRODUCTION
Efficient quantum algorithms for determining the ground and excited states of many-body systems are of fundamental interest to chemistry, condensed matter physics, and materials science [1][2][3][4]. The ability of quantum devices to represent N -body states with qubits scaling linearly in N make them particularly appealing for representing highly entangled states, as is common in systems with strongly correlated electrons. Therefore, quantum (and hybrid quantum-classical) algorithms offer an alternative to methods such as the density matrix renormalization group [5] (DMRG), selected configuration interaction [6,7] (SCI), determinant-based Monte-Carlo [8], variants of coupled-cluster (CC) theory [9,10] amenable to treating strong correlation [11][12][13][14], and multireference CC (MRCC) methods [15][16][17][18][19][20]. Although these classical algorithms can accurately predicted energies and properties of certain classes of strongly correlated systems, they still have high-order polynomial or exponential cost in the general case.
Since Feynman's proposal to use a controlled quantum system to carry out simulations [21], significant algorithmic and experimental advances have been made. The earliest demonstrations of quantum simulation for small molecules [2] utilized the quantum phase estimation algorithm [22][23][24] with Suzuki-Trotter decomposed time evolution [25,26] of an adiabatically prepared trial state. It is believed that some combination of these techniques will permit the efficient simulation [27][28][29][30] of certain classes of Hamiltonians [31], but that they may require deep circuits with high fidelity, a requirement incompatible with * nstair@emory.edu † francesco.evangelista@emory.edu current noisy intermediate-scale quantum (NISQ) devices [32]. Several low-depth quantum-classical hybrid algorithms have been developed for NISQ hardware. These algorithms prepare and measure properties of many-body states on a quantum device, but they store and optimize the parameters that define such states on a classical computer. The variational quantum eigensolver (VQE) approach [33][34][35][36] has been used in several landmark experiments, demonstrating quantum calculations on non-trivial molecular systems [3,[37][38][39][40][41]. In VQE, the ground state is approximated by a normalized trial state |Ψ ⟩ = Û (t) |Φ 0 ⟩, in which the unitary operator Û (t) depends on the parameter vector t and Φ 0 is (usually) an unentangled reference state. The VQE energy (E VQE ) is then obtained by minimization of the trial state energy expectation value as The VQE scheme employs an optimization algorithm running on a classical computer to minimize the energy expectation value, with all inputs (energy/gradients) being evaluated with the help of a quantum computer. An important advantage of VQE over classical manybody methods is the ability to use trial states that cannot be represented efficiently on a classical computer. VQE was initially implemented with an exponential operator ansatz inspired by unitary coupled-cluster (UCC) theory [42][43][44][45][46][47][48], but has more recently been extended to hardware-efficient [3] and qubit-space [49] UCC variants as well. We exclusively use the abbreviation UCC to refer to unitary coupled-cluster theory, and not unrestricted formulations of conventional coupled-cluster methods [50], which historically share this abbreviation. Other promising hybrid approaches include quantum imaginary time evolution [51,52], and quantum subspace diagonalization techniques [51,[53][54][55][56].
Despite the indisputable importance of VQE in the field of quantum simulation, there are a few drawbacks that make its practical application challenging to largescale problems. One such issue is the slow convergence of VQE due to noise in the measured energy and gradients, and the large-scale nonlinear nature of the optimization problem. These issues are compounded by the sizable number of total measurements needed for operator averaging [57]. Another challenge is the potentially large number of classical parameters and resulting circuit depth necessary to predict sufficiently accurate energies. These two problems are likely exacerbated in systems with strongly correlated electrons.
Progress addressing these deficiencies of VQE has been made on several fronts. For example, grouping commuting Pauli operators [3,35,[58][59][60], utilizing integral factorization strategies [61], and employing alternative bases [29,62] have been shown to reduce the number of measurements needed for operator averaging. Concurrently, computationally feasible approaches for measuring analytical gradients with quantum devices using the parameter-shift rule [63], or its recent lower-cost variant [64], have allowed gradient-based VQE to become potentially realizable on NISQ hardware. Other advances, of particular importance to this work, include VQE ansätze constructed iteratively, as done in the adaptive derivative-assembled pseudo-Trotterized VQE [36] (ADAPT-VQE) and the iterative qubit coupledcluster [65] (iQCC) methods. The primary advantage of ADAPT-VQE and iQCC is their ability to produce compact ansätze that result in fewer classical parameters, and shallower quantum circuits than those from UCC truncated to a given particle-hole rank. However, these advantages come at the cost of a greater number of energy and gradient evaluations for optimizing and selecting new unitary operators. Investigating more efficient ways to select important operators is an ongoing area of research [66,67].
In this work, we present an alternative to VQE for optimizing the amplitudes of a factorized form of the UCC ansatz (often referred to as Trotterized [35] or quantum [68] UCC), given by a product of exponential operators rather than the exponential of a sum of operators. We refer to this ansatz as disentangled UCC (dUCC)-a terminology borrowed from the field of Lie theory-to reflect the fact that it is not an approximation of UCC [69]. Inspired by the projective approach used in classical coupled-cluster theory [10,15], we propose an alternative trial state optimization algorithm that we deem the projective quantum eigensolver (PQE). PQE does not rely on variational minimization and therefore does not require evaluation of the energy gradients. Instead, PQE requires only the evaluation of residuals, that is, projections of the Schrödinger equations onto a linearly independent basis. As shown in this paper, residuals may be easily measured on NISQ devices with similar or fewer measurements than analytical gradients, and require quantum circuits that contain only one additional exponential term. We also propose a new selection scheme for identifying important operators based on the residual vector. This selected variant of PQE (SPQE) requires no pre-defined operator pool and employs only a small number of measurements to identify important operators. To demonstrate the practical advantages of PQE, we perform a comparison of VQE and PQE using a fixed dUCC ansatz for several molecular systems in the regime of weak and strong correlation, also considering the effect of stochastic noise. Finally, we compare SPQE against the ADAPT-VQE approach, selected configuration interaction, and the density matrix renormalization group.

A. The Projective Quantum Eigensolver approach
In this work, we propose to obtain the ground state of a general many-body system using a projective approach. Like in VQE, we approximate the ground state using a trial state |Ψ (t)⟩ = Û (t) |Φ 0 ⟩. After inserting the definition of the trial state in the Schrödinger equation and left-multiplying by Û † (t), we obtain the condition Projection onto the reference state Φ 0 yields the PQE energy (E PQE ) a quantity that is still an upper bound to the exact ground state energy. Projections onto the complete set of orthonormal many-body basis functions complementary to Φ 0 , here denoted as Q = {Φ µ }, yields a set of residual conditions where r µ is an element of the residual vector and µ runs over all elements of the many-body basis. Eqs. (3) and (4) form a system of nonlinear equations in the parameter vector t, that may be solved via a classical iterative solver. For an approximate ansatz with number of parameters less than the dimension of the Q space, Eq. (4) can be enforced only for a subset of the residuals. Then, the complete projection space Q can be partitioned into two sets: i) R, the space of basis functions Φ µ for which r µ = 0 is enforced, and ii) S = Q \ R the complementary space for which r µ may not be null. Figure 1 illustrates the connection between the PQE residual condition and the uncertainty in the groundstate energy estimated via Eq. (3). By the Gershgorin circle theorem, the difference between the exact (E) and the PQE (E PQE ) ground-state energy, |E PQE − E|, is bound by the radius ρ = ∑︁ µ̸ =0 |r µ |, where µ runs over the entire many-body basis, excluding the reference determinant. Therefore, when the residual is null (ρ = 0), the PQE energy is exact. When the PQE equation is satisfied only by a subset of the many-body basis-as in the case of an approximate trial state-the error |E PQE − E| is bound by the sum of the absolute value of the residual elements |r ν | with Φ ν ∈ S, for which the PQE equation is not satisfied.
Note that the residual condition [Eq. (4)] is satisfied by any eigenstate, and the Gershgorin circle theorem error bound applies also to excited states. A potential disadvantage is that PQE could converge on an excited state (an issue we did not experience in this study). However, this feature could be used to formulate excited state algorithms based on PQE, which use the residual condition as a criterion for convergence and do not require costly measurement of the variance, as is commonly done in VQE [35]. The residual vector is the first column of the transformed Hamiltonian matrix (first element excluded). (c) The difference between the approximate ground-state PQE energy (EPQE) and the exact eigenvalue (E) is bound by the radius ρ, which is equal to the 1-norm of the residual vector. The part of rµ that corresponds to states in the projection manifold R is null for the PQE solution, while elements involving projections on the space S = Q \ R is generally nonzero and determines the value of ρ.
The PQE is a general approach, however, in the following, we will focus on its applications to interacting fermions using a disentangled (factorized) form of the unitary coupled-cluster ansatz. We assume that our system is described by the general two-body Hamiltonian where â p (â † q ) is a fermionic annihilation (creation) operator, while h pq and v pqrs are one-electron and antisym-metrized two-electron integrals, respectively [70].

B. Traditional and disentangled unitary coupled-cluster ansätze
In UCC, the reference state is an easily-prepared single determinant |Φ 0 ⟩ = |ψ 1 ψ 2 · · ·⟩ specified by the occupied spin orbitals {ψ i }. A UCC unitary is parameterized using a pool of anti-Hermitian operators P = {κ µ : µ = 1, . . . , N pool op }. A generic anti-Hermitian operator κ µ = τ µ − τ † µ is defined in terms of the particle-hole excitation operators τ µ ≡ τ ab··· ij··· = â † a â † b · · · â j â i . Note that we have re-interpreted µ as the multi-index µ ≡ ((i, j, ..), (a, b, ..)) of unique excitations from hole/occupied (ψ i ψ j · · · ) to particle/unoccupied (ψ a ψ b · · · ) spin orbitals. Using this parameterization, when a cluster operator κ µ acts on the reference, it generates elements of the many-body basis (excited determinants) of the form and since in the case of a UCC (or dUCC) ansatz there is a 1-to-1 correspondence between operators and determinants, we may label them with the same index. Note that this operator basis satisfies the orthonormality condition In traditional UCC [42][43][44][45][46], the wave function is generated by an exponential operator assuming the cluster amplitudes t µ are real. In principle it is possible to construct a circuit that exactly implements the action of the UCC operator defined in Eq. (7), but in practice it is common to use a unitary with a simpler, and shallower, circuit. This is frequently accomplished using a factorized (disentangled) form of the UCC ansatz Because the operators κ µ do not commute, an ansatz of the disentangled form is uniquely defined by an ordered set (or subset) of operators A = (κ µi : i = 1, . . . , N op ) built from the pool P. The operators in A are then used to form an ordered product of exponential unitaries where t µi is the amplitude corresponding to the operator κ µi . The disentangled UCC ansatz may be viewed as a firstorder Suzuki-Trotter approximation to UCC; however, it was recently shown [69] that dUCC is substantially different from conventional UCC [42][43][44][45][46]. An arbitrary quantum state can always be represented in the form of Eq. (9) using particle-hole excitation/de-excitation operators [69]. However, only certain operator orderings A of a complete operator pool (containing up to N -body particle/hole excitations) can represent any quantum state. Nevertheless, orderings that do not satisfy this condition can exactly represent states that are close to the reference. In this work we assume that all operators appear at most once in A, but the more general case in which A contains repetitions has also been considered in other contexts [36,71].
For PQE formulated using a UCC or dUCC trial state, it is possible to show (see Appendix A) that for an exact ansatz the residual condition [Eq. (4)] and the VQE energy stationarity condition (∂E VQE /∂t µ = 0) are equivalent. However, for an approximate ansatz we find that the gradient of the PQE energy contains a contribution due to the nonzero residual elements corresponding to the subspace S: where t is a solution of the PQE equations in the projection space R. Suppose that R is chosen to be the space of single and double excitations and S its complement. Then this result shows that even if r µ = 0 for all singles and doubles, the gradient of E PQE with respect to singles and doubles may not be zero because residuals r ν corresponding to triple and higher excitations and the term ⟨Φ 0 | Û † (t)∂Û (t)/∂t µ |Φ ν ⟩ are generally not null. Therefore, PQE and VQE energies obtained from approximate ansätze will be different.
We note that the combination of PQE and dUCC produces energies that are additive for non-interacting fragments (size consistent) when using a localized basis. This property follows from the fact that dUCC excitation operators for non-interacting fragments act only on orbitals localized on each fragment, and therefore, commute. Consequently, the dUCC wave function is multiplicatively separable (as long as the order of the operators within a fragment is preserved). We have verified numerically that dUCC with singles and doubles trial states optimized with PQE energy are size consistent by performing calculations on H 4 + H 2 separated at a 1000 A and verified to within numerical convergence (10 −10 E h ) that the energy is additive in the fragments.

C. Numerical solution of the dUCC-PQE amplitude equation
To realize the PQE scheme on a quantum computer, the reference state, the Hamiltonian, and the unitary must be represented in a qubit basis via a fermionic mapping. After such a transformation, the Hamiltonian is a sum of the form Ĥ = ∑︁ ℓ h ℓ Ô ℓ , where h ℓ is an electron integral multiplied by a coefficient and is a unique product of σ x , σ y , or σ z Pauli operators acting on qubits q ℓi . Similarly, each term in the unitary exp(t µ κ µ ) may be implemented using a combination of one-and two-qubit operators following standard approaches [33][34][35].
To solve the PQE equations we measure the residuals corresponding to the operators contained in A on a quantum computer and update the parameter vector using a simple quasi-Newton iteration approach where the superscript "(n)" indicates the amplitude at iteration n. The quantities ∆ µ are standard Møller-Plesset denominators ∆ µ ≡ ∆ ab··· ij··· = ϵ i +ϵ j +. . .−ϵ a −ϵ b . . . where ϵ i are Hartree-Fock orbital energies. This update equation is derived in Appendix B using Newton's method and taking the leading contributions to the Jacobian to be the diagonal elements of the Fock operator [72]. It is further assumed that the amplitudes are small, so that the Jacobian can be approximated by terms linear in the operators κ µ and issues with non-commuting operators are avoided. Therefore, convergence of this quasi-Newton scheme is not mathematically guaranteed if one or more amplitudes are large. We found it useful to improve numerical stability and speed up convergence of the PQE equations, to combine amplitude updates via Eq. (11) with the direct inversion of the iterative subspace (DIIS) convergence accelerator algorithm [73,74].
It is important to note that the current formulation of PQE is compatible with any ansatz such that the metric matrix is the identity. In more general cases (e.g., when S µiµj is non-diagonal or singular), the PQE formalism requires a generalization of the amplitude update equations or the use of residual norm minimization instead of Eq. (4). These variants of PQE would allow one to consider ansätze that contain repeated operators in A, employ general many-body operators, or a basis of qubit operators. These extensions go beyond the scope of this work and will be considered in future studies. There are two advantages to the combination of the PQE and dUCC described above. Firstly, as we will show in the following subsection, one element of the residual vector (r µ ) can be evaluated with essentially the same resources required to measure the energy in VQE. Secondly, the magnitude of the residuals provides an indication of the importance of an excitation operator κ µ , which in turn may be used to define a selection procedure to form the sequence of unitiaries that enter in Û (t). The next two subsections describe these two points in detail.

D. Efficient measurement of the residual elements
For a trial state built from the ordered pool A, the number of the residual elements that must be evaluated to solve the PQE equations is equal to the size of the pool |A|. The PQE residuals can be expressed as the off-diagonal matrix elements of the operator H = Û † (t)Ĥ Û (t) as r µ = ⟨Φ µ | H |Φ 0 ⟩ (we use this notation throughout the paper, but we note that we never explicitly form the operator H on a classical computer). Then, in principle, the residuals can be measured on a quantum computer using a variant of the Hadamard test [75], but we have found an ancilla-free procedure in which these matrix elements are computed by measuring diagonal quantities. Acting on the reference with the operator e θκµ yields the state (13) noting that the above expression is valid because [47,48]). Taking the expectation value of the similarity transformed Hamiltonian with respect to Ω µ (θ) using θ = π/4, and noting that the wave function is real, leads to the following equation for the residual elements where All of these quantities are expectation values of H with respect to reference states that are easily generated with short quantum circuits. The evaluation of the exact residual via Eq. (14) has a cost similar to the evaluation of exact gradients in VQE via the shift rule [63,64].

E. Efficient operator selection
In this section, we generalize the PQE method to utilize a flexible dUCC ansatz built iteratively using a full operator pool. As shown in the case of VQE [36,65], significantly more compact and flexible approximations may be achieved if the operators that define the unitary are selected according to an importance criterion. To formulate a selected version of the PQE approach, we propose to combine information about the residual with a cumulative importance criterion. Since the residuals r µ are zero for an eigenstate, we propose to estimate the importance of the operators κ µ using the magnitude of the residual (|r µ |). However, instead of evaluating the importance of all the operators in the pool via operator averaging (like in gradient-based selection schemes [36,65]), we propose to sample a quantum state whose probability amplitudes encode the importance of all operators up to rank N .
Suppose that we have determined a unitary Û that satisfies the residual condition r µ = 0 for all κ µ in the current ordered set A. In our approach we prepare a (normalized) quantum state of the form |r⟩ = r 0 |Φ 0 ⟩+ ∑︁ µ r µ |Φ µ ⟩, where the quantities r µ are approximately proportional to the residuals r µ . When |r⟩ is represented in a qubit basis, there is a one-to-one mapping between elements of the computational basis and the states Φ µ . Therefore, a measurement of the state |r⟩ in the qubit basis will yield one of the states Φ µ with probability P µ = |r µ | 2 . Repeated measurement of the state |r⟩ provides a way to approximately determine the elements of the residual r µ with the largest magnitude, and the corresponding operators κ µ that should be included in the unitary. When this strategy is combined with an efficient way to prepare the state |r⟩, it is much more cost effective than evaluating all the elements of r µ not included in A via operator averaging [Eq. (14)].
Construction of the state |r⟩ would require one to apply the Hamiltonian, which is not a unitary operator. Therefore, we evaluate r using the unitary operator e −i∆tĤ = 1 − i∆tĤ + O(∆t 2 ) instead of Ĥ . By choosing a small time step, we can ensure that the nonlinear terms and errors due to the approximate implementation of e −i∆tĤ (e.g., via Trotterization) do not affect the residual to leading order in ∆t. The residual state can then be defined as The time-evolution operator may be approximated via Trotterization [25,26] in combination with low-rank representations of the Hamiltonian [76]. With a sufficiently large number of measurements M of the state |r⟩, we may approximate the values of the (normalized) squared residuals as where N µ is the number of times the state |Φ µ ⟩ is measured. Encoding the residual in a single quantum state allows us to efficiently sample the entire operator pool without the need to generate and store individual elements of the residual vector in memory. This distinctive feature makes it possible to employ this selection procedure with an operator pool that includes particle-hole excitation/de-excitation operators of arbitrary order.
To select important operators, we adopt a cumulative threshold approach that allows us to add a batch of operators at a time. Our goal is to iteratively construct a unitary that contains the fewest operators. This is realized with a selection procedure that adds the operators with the largest value of r µ to A (as motivated by the Gershgorin circle theorem bounds discussed in Sec. II A), and excludes all other operators in such a way that sum of their residuals squared is less than a threshold Ω 2 . Specifically, we enforce that where we have used the fact that |r µ | ≈ ∆t|r µ |. In practice, we sort the operators in ascending order according to |r µ | 2 , and starting from the first element, we discard operators until Eq. (17) is satisfied. The remaining operators are appended to the end of A in order of decreasing |r µ | 2 . The resulting operator ordering is consistent with the following renormalization transformation of the Hamiltonian which begins with the largest many-body rotation and gradually continues with smaller ones. Note that this ordering is the reverse of the one used in ADAPT-VQE and iQCC, where operators with the largest gradient are applied first to the reference. In applications of this selection scheme, we found that our operator ordering is more numerically robust and accurate than its reverse. This selection procedure is easily integrated in the PQE algorithm by performing a series of computations with increasingly larger ordered sets. When no new operators are added to A, the computation is considered converged, and the final operator set satisfies Eq. (17). The details of the selected PQE algorithm are discussed in Sec. II F. Throughout this work, we ignore errors that arise from finite measurement of the approximate residual |r⟩. However, since operator selection only requires an approximate determination of the |r µ | 2 values, it is not necessary to perform a large number of measurements. Indeed, in Appendix E, we discuss a simple strategy based on sampling |r⟩ a fixed number of times that performs as well as the exact scheme.

F. Outline of the selected PQE algorithm
The combination of the PQE approach with the selection procedure described in Sec. II E leads to a very efficient flexible ansatz quantum algorithm, which we refer to as selected PQE (SPQE). The selected PQE algorithm requires the simultaneous solution of the residual conditions and the selection of important excitation operators. To realize this scheme, we alternate micro-iterations to converge the residual equations (for the current ordered operator set) with macro-iterations that perform importance selection of new operators. The SPQE procedure is illustrated in Fig. 2 and consists of the following steps:  |Φ 0 ⟩. The number of times the state Φ µ is measured is accumulated in the variable N µ . These numbers are used to estimate the square residuals via Eq. (16), which are in turn used to select important excitations not included in the current pool A k . The current operator set A k and all the new selected operators are included in the new set A k+1 . When forming this new ordered set, we append the new operators-sorted in decreasing value of the approximate squared residual-to A k . If the sum in Eq. (17)  ], as well as the updated energy [E (k+1) ]. Increase the macro-iteration number by one (k ← k + 1) and go to Step 2.
All VQE and PQE methodologies were implemented in a development branch of the open source package QForte , and utilize its state-vector simulator. All VQE calculations use a micro-iteration convergence threshold ω g = 10 −5 E h for the gradient norm ∥g∥ and all PQE calculations use a micro-iteration threshold ω r = 10 −5 E h for the residual norm ∥r∥. Note that when the κ µ operators are mapped to a qubit basis they can be expressed as a sum of Pauli operator strings that commute [77], and therefore, a single operator e tµκµ can be implemented exactly as a product of exponentials without invoking the Trotter approximation.

A. Comparison of PQE and VQE with a disentangled UCC ansatz
Our initial goal is to compare the performance of PQE and VQE using a unitary coupled-cluster trial state truncated at a given particle-hole excitation level. We test these two methods on a family of linear hydrogen chains (with identical nearest-neighbor distance) ranging from four to ten atoms, both near their equilibrium geometries (r H-H = 0.75Å) and stretched geometries (r H-H = 1.5Å). Hydrogen models such as these have been studied experimentally with VQE [78] and have recently been used as a benchmark for both quantum [36,55] and classical [79][80][81] algorithms. Figure 3 shows the energy convergence of PQE and VQE using a disentangled UCC ansatz with singles and doubles (dUCCSD). All calculations employed restricted Hartree-Fock (RHF) orbitals from the quantum chem-   istry package Psi4 [82], and Pauli-operator Hamiltonians obtained via the Jordan-Wigner transformation implemented in QForte [83]. To achieve optimal performance for both VQE and PQE, we employ the Broyden-Fletcher-Goldfarb-Shannon (BFGS) algorithm [84][85][86][87] (as implemented in the SciPy [88] scientific computing library) with analytical gradients for VQE, and DIIS [73,74] to accelerate amplitude convergence of PQE. These computations use the same operator ordering for both approaches, with all amplitudes initialized to zero. The ordering of the operators e tµ i κµ i entering Eq. (9) is defined by the binary representation of the corresponding determinants |Φ µi ⟩ = τ µi |Φ 0 ⟩ in the occupation number representation. Because the disentangled UCCSD state cannot exactly parameterize an eigenstate of the Hamiltonian, the numerically converged PQE and VQE energies are not identical. Nevertheless, for all the cases we examined the converged PQE and VQE energies differ by less than 10 −6 E h .
Near the equilibrium geometry, we find that the PQE energy converges significantly faster than the VQE energy with number of residual vs. gradient evaluations, respectively. For example, to converge the near-equilibrium H 10 energy to 10 −6 E h , PQE requires only seven residual evaluations, while VQE requires approximately 23 gradient evaluations. In the case of VQE, we also observe that the number of required gradient evaluations grows with system size, with H 10 taking twice as many gradient vector evaluations than H 4 to converge. On the contrary, PQE computations converge with similar speed for all equilibrium hydrogen systems. Plots of the energy change vs. the norm of the residual/gradient vector show similar trends and are reported in Appendix C.
At the stretched geometry, strong correlation effects cause the disentangled UCCSD trial state to perform more poorly, with both PQE and VQE, yielding energy errors that range from 1.39 mE h (for H 4 ) to 13.59 mE h (for H 10 ). We see that PQE converges slightly more slowly in the stronger correlation regime, with stretched H 10 requiring 13 residual evaluations (instead of seven) to converge the energy to 10 −6 E h . However, PQE always converges faster than VQE, with the latter requiring 19 gradient-vector evaluations to converge stretched H 10 to the same accuracy level. We also find similar trends in PQE and VQE convergence for BeH 2 , for which convergence data can be found in Appendix C.
In summary, this initial set of results suggests that for a given trial state, optimization via PQE is faster than VQE and less dependent on the number of parameters to optimize. We expect this to be the case also for VQE based on numerical gradients or gradient-free optimization methods, since these two variants are known to be slower compared to the the BFGS approach adopted here [77].

B. Effect of stochastic errors on the convergence of PQE and VQE
The results presented in the previous section assumed error-free quantum gates and arbitrarily precise measurements. In practice, calculations performed on NISQ hardware are affected by decoherence errors, poor gate fidelity, readout errors, and loss of precision due to insufficient measurements. These sources of error will lead to incorrect gradients and residuals that are then passed to a classical optimizer. Therefore, it is interesting to compare the resilience of PQE and VQE procedures when the residuals and gradients are affected by stochastic errors.
To model the presence of errors, we modify the PQE procedure by adding to the residual vector a stochastic error sampled from a Gaussian distribution with standard deviation σ [N (0, σ 2 )] For VQE, we similarly add stochastic noise to the exact energy gradients. Using Eq. (19) as a noise model mainly emulates errors that arise from finite measurement. Because the inexact residuals r measured µ are used to update the cluster amplitudes via Eq. (11), this noise model also gives rise to control errors, or errors that refer to the difference between the unitiaries for noiseless updated amplitudes Û (t) and noisy updated amplitudes Û (t + ∆t). Control errors due to finite measurement have been modeled this way in previous studies [77] and will always propagate through optimization on physical hardware. We note, however, that using Eq. (19) exclusively as a noise model is insufficient to capture more nuanced or device-specific errors such as decoherence. We tested the performance of PQE and VQE under noise by performing a batch of 50 computations on the linear H 4 molecule with nearest-neighbor distance set to 1.0Å. We use the disentangled UCC ansatz with up to quadruple excitations, which spans the full operator set for this system. Figure 4 shows a comparison of PQE and VQE optimized with various levels of noise, controlled by the magnitude of σ in Eq. (19). We find that for all values of σ > 0, the energy convergence of PQE is essentially identical to that of noiseless PQE until some point, after which the energy error hovers around a finite value. Similar behavior can be seen for convergence of the residual vector. We find that VQE has similar characteristics to PQE in the presence of noise, but that it is able to achieve slightly more accurate energies for a given σ value. With the same noise level, however, VQE generally requires two to three times the number of gradient evaluations as the number of residual evaluations required by PQE. An important aspect of this comparison is that, for a given σ, both the residual and gradient vectors yield comparable asymptotic errors in PQE and VQE, respectively. Conversely, PQE and VQE computations of the comparable energy accuracy require similar precision in the measurement of the residual and gradient vectors. In Appendix D we use this result to estimate the relative cost of PQE and VQE via a formal analysis.

C. Selected PQE based on a full dUCC operator pool
Here we compare the results of selected PQE (SPQE) with an arbitrary order particle-hole operator pool, and ADAPT-VQE, which unless otherwise noted, uses a generalized singles and doubles pool (operators of the form â † p â q and â † p â † q â s â r , where the indices p, q, r, s run over all spin orbitals). To test these methods, we compute the energy as a function of bond distances for: 1) the symmetric dissociation of the linear BeH 2 molecule and 2) the symmetric dissociation of a chain of six hydrogen atoms. In both cases, there is a build up of strong correlation effects as the bond length increases. For each system, we report two sets of results for SPQE using the cumulative threshold Ω = 10 −1 and 10 −2 E h , and two sets of ADAPT-VQE results with gradient threshold 10 −1 and 10 −3 E h (ϵ 1 and ϵ 3 in the original notation used by the authors). Since we later compare these methods to classical approaches formulated in a determinant basis, we do not employ a pool of spin-adapted operators.
The dissociation curve of BeH 2 shown in Fig. 5 (a) demonstrates that both SPQE and ADAPT-VQE are able to achieve significantly smaller energy errors than dUCCSD, while using only 10-20 more parameters. Although the two approaches employ different selection schemes, the ADAPT-VQE(ϵ 3 ) and SPQE(Ω = 10 −2 ) approaches produce compact trial states with a similar number of classical parameters and comparable errors. The same trends are seen in symmetric dissociation curve of H 6 , which is shown in Fig. 5 (b). However, for this system we find that achieving sub-mE h accuracy-particularly with the onset of strongly correlation at r H−H values greater than 1.5Å, requires a number of parameters that approaches the size of the full Hilbert space (200) for both ADAPT-VQE and SPQE. The need to saturate Hilbert space to accurately describe H 6 is likely due to how small this example is, and it speaks FIG. 4. Energy and residual/gradient norm (∥r∥)/(∥g∥) convergence of dUCCSDTQ wave functions optimized with PQE/VQE with various amounts of stochastic noise added to the residuals/gradients. Energy error is relative to FCI. σ controls the degree of noise and is the standard deviation of the normal distribution, centered at the exact residual/gradient value, from which all residuals/gradients used in the calculations are randomly sampled [see Eq. (19)]. Values at each PQE/VQE iteration are averages over 50 runs on H4 at rH−H = 1.0Å. Error bars denote one standard deviation. more to the ability of the trial states to produced compact representations rather than the performance of these algorithms in optimizing such ansätze.
The most noticeable difference between SPQE and ADAPT-VQE can be seen in the bottom panel of Fig. 5 (a)-(b): for trial states with comparable number of parameters, SPQE requires significantly fewer residual element evaluations than gradient element evaluations in ADAPT-VQE. For example, at a Be-H bond distance of 1.65Å, both methods produce very similar energy errors using almost the same number of parameters, but ADAPT-VQE(ϵ 3 ) requires the evaluation of 35155 elements of the gradient, whereas SPQE(Ω = 10 −2 ) requires only 1220 elements of the residual. Importantly, we note that the bottom panels of Fig. 5 exclusively count the number of elements of the gradient or residual required by the optimization, and do not include the additional measurements required for operator selection.
Since ADAPT-VQE and SPQE select new operators from their pools using different importance criteria, it is not possible to perform a direct comparison of their performance using fixed thresholds. To facilitate this comparison, in Table I we report SPQE results using two values of Ω (10 −1 and 10 −2 E h ) together with ADAPT-VQE results obtained using an ansatz with the same number of parameters as SPQE. We include ADAPT-VQE results obtained using both a generalized singles and doubles operator pool (GSD), and a particle-hole singles and doubles pool (SD). The results in Tab The evaluation of fewer residual elements in SQPE will correspond to approximately the same savings in total number of measurements (see Appendix D). A final important aspect to compare between SPQE and ADAPT-VQE is the number of native CNOT gates, which we consider as a proxy for circuit depth. Tab. I reports the number of CNOT gates for the converged trial states. These numbers overestimate the actual gate count since they ignore optimizations such as the cancellation of Jordan-Wigner strings [89], especially for three-and higher-body operators. We see that at the larger threshold values (Ω = 10 −1 for SPQE and ϵ 1 = 10 −1 for ADAPT) the number of CNOTs for all three approaches are relatively similar, and are generally within a factor of 2 of one another. However, at tighter thresholds the SPQE circuit contains significantly more CNOT gates than the one for ADAPT-VQE. For example, r H−H = 2.0Å), 114 of the 169 operators used in SPQE are three-body or higher, while ADAPT-SD and ADAPT-GSD of course only contain only up to two-body operators. Consequently, the SPQE(Ω = 10 −2 ) circuit contains more CNOT gates (226768) than ADAPT-VQE-GSD (14776). As discussed in Appendix D, this large difference in CNOT count is due to the growth in cost to implement the exponential of n-body second-quantized operators as a function of n (and the lack of circuit optimization). Nevertheless, it is important to note that the systems studied here are not large enough to draw definitive conclusions about the relative performance of SPQE and ADAPT-VQE. For example, the pool of generalized singles and doubles (GSD) for H 6 contains 870 operators, a number significantly larger than the size of the full Hilbert space (200). This scenario is unlike most systems of interest, where the number of generalized singles and doubles is much less than the number of particle-hole operators. Unfortunately, our attempts to compare numbers for systems larger than BeH 2 and H 6 were unsuccessful due to the high computational cost of simulating ADAPT-VQE.
A second aspect we investigate, is the ability of the selected PQE approach to compactly represent wave functions for systems displaying strong correlation effects, were many-body methods commonly fail due to the breakdown of the mean-field approximation. We compare the performance of SPQE with the adaptive configuration interaction [90] (ACI) and DMRG using data generated in our recent benchmark study of hydrogen systems [81]. In this work, we characterize the resource requirements of a computational method X with the accuracy volume (V X ), defined as the smallest number of parameters necessary to achieve a given energy error per electron (here taken to be 10 −4 E h /electron). We consider three H 10 model systems: a one dimensional (1D) linear chain, a two dimensional (2D) triangular lattice and a three dimensional (3D) close-packed pyramid. Using a minimal STO-6G basis and a 1 qubit to 1 spinorbital mapping, the 1D, 2D, and 3D H 10 models are represented with 2 20 computational basis states, but are restricted to 31752 and 15912, and 15912 determinants, respectively, after accounting for spin and Abelian point-group symmetries. Figure 6 (a) displays the SPQE, ACI, and DMRG energy error as a function of classical variational parameters for the 1D H 10 system at a stretched bond length of r H−H = 1.50Å. DMRG, which is ideally suited to simulate gapped 1D systems, affords the most compact wave function for the 1D chain, with an accuracy volume of only 176 variational parameters (V 1D DMRG = 176). The SPQE exponential ansatz is less compact than the DMRG matrix product state, with an accuracy volume V 1D SPQE = 3510 parameters. We observe that the ACI wave function-a linear ansatz built by selecting the determinants with the largest energy contributions-gives the least compact representation, such that V 1D ACI = 22989 at the same level of accuracy. When the ACI energy is augmented with a second-order perturbative correctionwhich accounts for determinants excluded from the wavefunction expansion-a more compact ansatz is sufficient to achieve an energy error of 10 −3 E h . This suggests that it might be valuable to formulate a classical perturbative correction to the SPQE energy to help capture even more correlation energy.
Results for the 2D H 10 lattice [ Fig. 6 (b)] show that SPQE yields a more compact representation than DMRG or ACI until an energy error of approximately 1 mE h . For the 2D system the accuracy volume follows the order SPQE (4193) ≈ DMRG (4218) < ACI (11122). At energy errors less than 1 mE h DMRG still affords the most compact representation for the 2D system. Finally, for the 3D H 10 pyramid lattice shown in Fig. 6 (c), we find that DMRG yields the most compact representation (V 3D DMRG = 3549), but not nearly by the same margin as for the 1D system. Again, the SPQE has a significantly more compact wave function than ACI such that V 3D SPQE = 5544 and V 3D ACI = 10519. While further investigation of larger strongly-correlated systems will be necessary, it is encouraging to see that SPQE performs similarly to or better than two state-of-the-art classical methodologies when applied to the 2D and 3D H 10 sys-

tems.
We have also included in Fig. 6 energy errors for several classical coupled-cluster variants. Specifically we have included results from Ref. [81] for CC with singles and doubles excitations [91] (CCSD), CCSD with perturbative triples [92] [CCSD(T)], and the completely renormalized CC with perturbative triples [12] [CRCC (2,3)]. These CC methods are (generally) more computationally affordable than the other methods. In the 2D and 3D systems, CCSD performs comparably to dUCC-SQPE (with the same number of parameters), while CCSD(T) and CRCC(2,3) produce more accurate energies. However, all coupled-cluster energies are more than 10 mE h off from the FCI energy. Fig. 6 also reports results computed using Mukherjee MRCC with singles and doubles [93,94] (Mk-MRCCSD) and Mk-MRCCSD augmented with perturbative triples [95] [Mk-MRCCSD(T)] using an active space containing the highest occupied and lowest unoccupied Hartree-Fock orbitals. The Mk-MRCC results improve upon single-reference coupled-cluster methods at the cost of doubling the number of cluster amplitudes. The Mk-MRCC methods are particularly accurate for the 1D system, where they produce errors of the order of 2-3 mE h . Despite the improvement shown by the multireference CC methods, it is important to note that these methods have a computational cost that scales exponentially with the number of active space determinants.

IV. CONCLUSIONS
In this work, we present a new NISQ-friendly algorithm-the projective quantum eigensolver (PQE)to compute the ground state of a many-body problem using disentangled (factorized) unitary coupled-cluster trial states. The PQE approach consists of a nonlinear optimization problem whose solution requires the evaluation of projections of the Schrödinger equation onto a manybody basis (residual vector) but still gives energies that are a variational upper bound to the ground state energy. We show how to efficiently evaluate the residual vector via measurement of simple expectation values, with a cost that is twice that of an energy evaluation (per element). For small molecular systems, we find that PQE and VQE with a fixed dUCCSD trial state converge to nearly identical energies; however, the number of residual evaluations required by PQE is smaller than that of the gradient evaluations needed by VQE. PQE shows similar resiliency as VQE and still converges more rapidly in the presence of stochastic noise.
To treat strongly correlated electrons, we introduce a selected variant of PQE in which the trial state is constructed iteratively by adding batches of important operators. The resulting SPQE algorithm can construct efficient unitary circuits like ADAPT-VQE, but it requires orders of magnitude fewer residual element evaluations than gradient element evaluations for the latter. In SPQE, the selection of new operators is done according to the magnitude of the elements of the residual vector, and is performed by sampling a quantum state that directly encodes in its probability amplitudes the importance of the entire operator pool. Because the selection cost in SQPE is not affected by the size of the operator pool, the unitary can include operators of rank up to the total number of particles.
Finally, we compare the energy convergence with the number of parameters for 1D, 2D, and 3D H 10 lattices using SPQE and two classical methods well suited to treat strong electronic correlation: the adaptive configuration interaction and the density matrix renormalization group. Given a target accuracy of up to approximately 1 mE h , we find that SPQE produces significantly more compact trial states for the 2D system than ACI and is comparable to DMRG. However, DMRG affords the most compact wave function parameterization in 1D, for accuracies below 1 mE h in 2D, and by a much smaller margin, in 3D. Taken together, PQE and SPQE are very promising tools for studying many-body systems both in the strong and weak correlation regimes using NISQ hardware.
In summary, the PQE approach is a viable and more economical alternative to variational quantum algorithms. In its current formulation, PQE can be applied to any trial state generated by exponentiating a set of linearly independent operators with identity metric matrix, as it is the case for disentangled unitary coupledcluster ansätze. For these trial states, methods to reduce the number of measurements and exploit symmetries [59,60,96] developed for VQE can likewise be used to improve PQE. Interesting extensions of PQE include a generalization to unitaries that may contain repeated operators, that use generalized excitations/de-excitations pools, and hardware-efficient ansätze. In particular, a promising research direction is the formulation of a selected PQE using a basis of general one-and two-body operators, which could yield trial states with lower circuit depth than the current formulation. Within the greater ecosystem of quantum algorithms, PQE could be used to determine initial guess amplitudes for subsequent optimization via VQE. Additionally, similarly to VQE, PQE can be used as an alternative to adiabatic approaches to prepare initial states for quantum phase estimation. Although we only explore applications of PQE for quantum many-body simulation, the framework outlined by Eqs. (3) and (4) can be used to solve a variety of eigenvalue problems. With appropriate modifications, for example, PQE could be employed to diagonalize covariance matrices (after quantum encoding [97,98]) for use in machine learning or data-analysis. Moreover, because there is no requirement that the PQE (or SPQE) trial states have low-entanglement, PQE could be used as an alternative to methods such as quantum principle analysis [98], or variational quantum state diagonalization [99]. Most importantly, the most immediate impact of PQE could be speeding up quantum computations on current or near-term devices.
The second term in Eq. (A3) is null due to the residual condition. Applying these two simplifications we arrive at the gradient expression [Eq. (10)] reported in Sec. II B, containing only the last term of Eq. (A3). The term ∂tµ |Φ 0 ⟩ that multiplies the residual is generally non-null for both the UCC and dUCC ansätze. In the case of dUCC, the term corresponding to the derivative with respect to i-th amplitude t µi , is given by: The states to the left and right of the operator κ µi are general many-body states that potentially span the entire Hilbert space, and therefore, this quantity is gener-ally non-null. For completeness, we also report the same expression for the case of traditional UCC, which may be obtained using the derivative of the exponential map: ), the number of residual/gradient evaluations (N · X ), and the norm of the residual/gradient (∥ · ∥). The FCI energies at 1.0 and 2.0Å are −15.65068726 and −15.60861964 E h , respectively. All calculations use a STO-6G basis.  The preferred approach to minimize the energy in VQE is via gradient-based algorithms. The gradient of the dUCC energy with respect to a cluster amplitude t µ is given by Equation (D1) has the form of an off-diagonal matrix element of the Hamiltonian, and Romero et al. [77] showed that it may be measured on a quantum device using one ancilla qubit and four time the cost of measuring the en-ergy. With this algorithm, computing the gradient vector so that the 2-norm of the error is within precision ε grad requires a total number of measurements (m grad ) that is bound by the inequality [77] m grad ≤ 4N par The total number of measurements, and proportionally the runtime, required for a VQE calculation is dominated by the computation of the gradient vector and it is proportional to the number of gradient vector evaluations (N VQE grad ) required to converge the energy. This implies that the total number of VQE measurements (m VQE ) is bounded approximately by This estimate ignores the evaluation of the VQE energy, which has a cost inferior to that of computing one element of the gradient vector. More recently, Kottmann et al. [64] showed that the analytic gradient may be computed with a cost essentially equal to that of two energy evaluations using the so-called parameter-shift-rule [63]. This procedure avoids the use of an ancilla qubit and the number of measurements required still satisfies the bound expressed in Eq. (D2).
In the case of PQE, the number of measurements m res needed to compute the residual vector with precision ε res has an upper bound given by This estimate takes into account the fact that the residual can be evaluated as the sum of three terms (with different prefactors), and that E 0 in Eq. (14) only needs to be measured once. The total number of PQE measurements (m PQE ) is bounded by where N PQE res is the number of PQE residual vector evaluations. Assuming that the energy gradients in VQE and the residual in PQE are measured with the same precision (ε grad ≈ ε res ), we estimate that m res ≈ 3 4 m grad . This result suggest that PQE should have a similar or perhaps slightly smaller cost per iteration than VQE. However, a more important factor in determining the relative performance of VQE and PQE is the number of residual/gradient evaluations required, which as shown in Sec. III, favors PQE over VQE.
A detailed comparison of the adaptive variants of VQE (e.g., ADAPT-VQE) and PQE is more complex due to the significant differences in the form of the ansatz and the selection procedure used in these two methods. At iteration k, the ADAPT-VQE approach selects the operator κ µ with the largest absolute energy gradient. This selection scheme requires evaluating the gradient g (k) µ for all the operators in the pool Because most sub-terms of the Hamiltonian will commute with κ µ , a relatively small number of ⟨κ µ Ô ℓ ⟩ and ⟨Ô ℓ κ µ ⟩ terms need to be measured. Both Ĥ and the operator pool {κ µ } contain of the order N 4 elements (assuming a pool of general one and two-body operators). However, it has recently been pointed out [67] that if the quantity [Ĥ , κ µ ] is decomposed in terms of the reduced density matrices, then one can determine the pool gradients [Eq. (D6)] with the evaluation of a number of Pauli terms that scales as N 6 (again assuming a pool of general one and two-body operators). In ADAPT-VQE this cost must then be multiplied by the number of iterations performed. It is important to note that the ADAPT-VQE macro-iteration convergence threshold ϵ α = 10 −α is based on the norm of the vector of pool gradients [Eq. (D6)], such that ADAPT-VQE is considered converged when ∥g pool ∥ ≤ ϵ α . The number of measurements of the approximate residual vector (|r⟩) for the purpose of selection in SPQE is a parameter of a computation. In Appendix E, we show that a probabilistic estimate for the number of measurements required to converge SPQE with a threshold of Ω is of the order of (∆tΩ) −2 .
A trade-off of using three-and higher-body operators in SPQE is a greater circuit depth since an operator κ ab··· ij··· of many-body rank n becomes a sum of 2 2n−1 Pauli strings after Jordan-Wigner mapping to the qubit basis. Since all Pauli strings that are generated in this mapping commute, the unitary exp(θ κ ab··· ij··· ) can be written as the product of 2 2n−1 exponentials of operators containing Pauli strings of length 2n.
To analyze the compromise between the rank of the operator pool and the compactness of the ansatz in ADAPT-VQE and SPQE we consider the example of a three-body operator exp(t κ abc ijk ). In ADAPT-VQE this term may be approximated using general two-body operators as e t κ abc ijk = e t [κ ae ij ,κ bc ke ] ≈ e t κ ae ij e t κ bc ke e −t κ ae ij e −t κ bc ke , (D7) where the last term is a lowest-order Trotter approximation of the exponential of the commutator [κ ae ij , κ bc ke ] (with e ̸ = a, b, c). The last term in Eq. (D7) is implemented as a circuit that contains four different parameters and a product of 32 exponentials of Pauli strings of length four. The same three-body excitation is represented in SPQE using a single parameter and a longer circuit as a product of 32 exponentials of Pauli strings of length six. This comparison suggests, in accordance with the results of our study, that higher-body operators are represented less efficiently with an arbitrary particle-hole operator pool than a general singles and doubles operator pool.
Appendix E: Reduced-cost estimation of the approximate residual in selected PQE This appendix explores various methods to reduce the number of measurements M required to compute the approximate the residuals r µ used in the selected PQE method (see Sec. II E). Selection requires the identification of the elements of the approximate residual that corresponds to projections onto excited determinants. However, for small values of ∆t, the state |r⟩ is dominated by the reference determinant Φ 0 , and consequently, the measurement of important missing excitations may become inefficient. In practice one can consider an alternative convergence criterion for SPQE based on performing a fixed number of measurements M Ω on the state |r⟩. In such an approach, at each iteration k, the operators κ µ whose corresponding determinants |Φ µ ⟩ are measured at least once (over all M Ω measurements) are added to A. Because the residual magnitudes go to zero as the eigenstate is better approximated in successive k iterations, it becomes increasingly unlikely that any determinants besides the reference |Φ 0 ⟩ will be measured. The SPQE algorithm can then be considered converged when all M Ω measurements yield the reference state. Starting from Eqs. (16) and (17), and making the assumption that at convergence only a single determinant Φ µ is measured that is not the reference (i.e., ∑︁ µ N µ = 1), one can use a number of measurements to probabilistically test convergence of the residual vector to within the threshold Ω. In practice we find that using Eq. (E1) works well compared to the exact threshold given in Eq. (17). Figure 8 compares the energy convergence with number of selected operators/parameters using both convergence criterion corresponding to the same value of Ω.