Quantum computed moments correction to variational estimates

The variational principle of quantum mechanics is the backbone of hybrid quantum computing for a range of applications. However, as the problem size grows, quantum logic errors and the effect of barren plateaus overwhelm the quality of the results. There is now a clear focus on strategies that require fewer quantum circuit steps and are robust to device errors. Here we present an approach in which problem complexity is transferred to dynamic quantities computed on the quantum processor - Hamiltonian moments, $\langle H^n\rangle$. From these quantum computed moments, estimates of the ground-state energy are obtained using the ''infinum'' theorem from Lanczos cumulant expansions which manifestly correct the associated variational calculation. With system dynamics encoded in the moments the burden on the trial-state quantum circuit depth is eased. The method is introduced and demonstrated on 2D quantum magnetism models on lattices up to 5 $\times$ 5 (25 qubits) implemented on IBM Quantum superconducting qubit devices. Moments were quantum computed to fourth order with respect to a parameterised antiferromagnetic trial-state. A comprehensive comparison with benchmark variational calculations was performed, including over an ensemble of random coupling instances. The results showed that the infinum estimate consistently outperformed the benchmark variational approach for the same trial-state. These initial investigations suggest that the quantum computed moments approach has a high degree of stability against trial-state variation, quantum gate errors and shot noise, all of which bodes well for further investigation and applications of the approach.

The variational principle of quantum mechanics is the backbone of hybrid quantum computing for a range of applications. However, as the problem size grows, quantum logic errors and the effect of barren plateaus overwhelm the quality of the results. There is now a clear focus on strategies that require fewer quantum circuit steps and are robust to device errors. Here we present an approach in which problem complexity is transferred to dynamic quantities computed on the quantum processor -Hamiltonian moments, H n . From these quantum computed moments, estimates of the groundstate energy are obtained using the "infinum" theorem from Lanczos cumulant expansions which manifestly correct the associated variational calculation. With system dynamics encoded in the moments the burden on the trial-state quantum circuit depth is eased. The method is introduced and demonstrated on 2D quantum magnetism models on lattices up to 5×5 (25 qubits) implemented on IBM Quantum superconducting qubit devices. Moments were quantum computed to fourth order with respect to a parameterised antiferromagnetic trial-state. A comprehensive comparison with benchmark variational calculations was performed, including over an ensemble of random coupling instances. The results showed that the infinum estimate consistently outperformed the benchmark variational approach for the same trial-state. These initial investigations suggest that the quantum computed moments approach has a high degree of stability against trial-state variation, quantum gate errors and shot noise, all of which bodes well for further investigation and applications of the approach.

I. INTRODUCTION
Quantum computers represent a new paradigm for computing that is witnessing rapid advances in both hardware and software. Fully programmable devices are emerging, and evidence of information processing at a scale that competes with supercomputers for certain sampling problems has been reported in a quantum device comprising 53 qubits [1]. The major challenge of the field is to demonstrate "quantum advantage" for real-world problems. There are approaches to a range of potential application areas, including bioinformatics [2], chemistry [3,4], optimisation [5], finance [6,7] and machine learning [8], to name a few. However, in the short to medium term, quantum computer (QC) technology will be constrained to the so-called noisy intermediate scale quantum (NISQ) regime [9] -where the performance of QC devices will inevitably be dominated by the level of logic precision inherent in the hardware.
In NISQ devices, the quest for quantum advantage is challenged by errors in logic and read-out which place severe restrictions on the number of time-steps, or "depth", of any given quantum circuit before the results are scrambled. Hybrid quantum algorithms such as the Variational Quantum Eigensolver (VQE) [10] or the Quantum Approximate Optimisation Algorithm (QAOA) [5] adapt variational-style hybrid approaches to the problem cast in Hamiltonian form. However, the application to realworld problems generally requires relatively deep quantum circuits with the accumulation of quantum logic errors leading to barren plateau effects [11,12]. In NISQ * lloydch@unimelb.edu.au devices, quantum circuit depth is perhaps the most precious quantum resource.
In this work we introduce an alternative method for computing the lowest energy of a problem Hamiltonian system, H, based on minimising circuit depth by transferring complexity to the computation of moments of the Hamiltonian, H n , with respect to a given trial-state ( Figure 1). Such quantities are central to certain nonperturbative approximation schemes in many-body theory, but are generally difficult to compute classically as the problem scales. The key in our approach is to reserve the quantum resource for the direct computation of the Hamiltonian moments, which are then used to determine ground-state energy estimates that correct the variational result (first moment, H ). In our initial implementation of the quantum computed moments (QCM) method, the results obtained appear highly robust to device noise and provide a significant correction to the energy of the benchmark variational result. This paper is organised as follows: Section II describes the problem in Hamiltonian form and the variational approach in the context of hybrid quantum computing. In Section III we review the background of the QCM approach: the connection of Hamiltonian moments to the Lanczos expansion approach and the infinum theorem that provides energy estimates to a given moment order. In Section IV we define the problem Hamiltonian for the demonstration -quantum spin systems on 2D latticesand provide details on the tensor product basis set construction and scaling. Results obtained by computing moments on IBM Quantum devices and associated simulations are provided in Section V, and conclusions drawn in Section VI.

II. HAMILTONIAN FORM AND THE VARIATIONAL LIMIT
We consider problems which, when converted to a quantum context, can be reduced to finding the groundstate energy of an equivalent problem Hamiltonian H over strings of qubit operators {I, X, Y, Z} -this general class of problems encompasses a number of applications. We write the problem Hamiltonian over q qubits as: The sum is over Pauli strings [P i ] ≡ σ 1 (i)σ 2 (i)...σ q (i), where for the mth qubit we have σ m (i) = (I, X, Y, or Z) defined by the problem, and each term has a weight w i .
The solution is the lowest energy state of H, which we denote E 0 . This is, in general, a difficult task given the Hilbert space dimension grows as 2 q .
Approaches to solve such problems based on the wellknown variational principle in quantum mechanics have gained widespread appeal in the NISQ era of quantum computing. The variational quantum eigensolver (VQE) [10] is a hybrid approach to finding an approximation to E 0 on a quantum computer. VQE begins by setting up a quantum circuit to create the parameterised trial-state |φ trial ( θ) over a set of parameters θ. The expectation values [P i ] θ , in the state |φ trial ( θ) are estimated term by term via repeated initialisation and measurement of the QC, and summed to produce H θ . This procedure is incorporated into a classical loop minimising the re-sult with respect to the trial-state parameters to produce the lowest upper bound H θ on E 0 . In VQE, the state |φ trial ( θ) is implemented as a sequence of p sub-circuits involving mixing and entangling operations governed by parameters θ(k) in the kth sub-circuit block. As the number of circuit blocks increases, the total circuit depth grows, as does the set of parameters { θ(1), θ(2), ... θ(p)}, the trial-state becomes more and more complicated in order to better approximate the ground-state of the problem. Working against this, errors in a NISQ device accumulate in the output, thereby restricting the maximum depth of the trial-state circuit and limiting convergence of the variational procedure [11,12]. This is generally a critical barrier to overcome in practical applications.
While the variational procedure is an obvious and time-proven place to start, it is clear that the quantum resource must be used as efficiently as possible. We focus on valuing quantum circuit depth by exploiting dynamics in the Hamiltonian as the definitive generator of the ground-state estimate to ease the complexity of the trial-state. The quantum computed moments approach produces a correction to the variational result, effectively trading quantum circuit depth for an increased number of measurements and classical post-processing.

III. QUANTUM COMPUTED MOMENTS (QCM) APPROACH
The focus on using a quantum computer to directly compute moments of the Hamiltonian H n with respect to a trial-state is inspired by a cumulant expansion of the Lanczos method [13]. Uncovered for extensive systems in the context of lattice gauge field theory, the cumulant expansion of the Lanczos tri-diagonalised form allows for the ground-state of the diagonalised system to be obtained via an "infinum" theorem [14]. Typically, the computation of H n scales poorly with classical computing resources for non-homogenous H -the utility of the infinum approach has been limited to homogeneous systems. However, the emergence of quantum computers sheds new light on the overall approach and its possibilities.
We begin the description with the well-known Lanczos recursion method [15] for estimating the lowest energy eigenvalue(s) of a system. It has recently been considered in the quantum computing context with the introduction of "imaginary-time" evolution algorithms [16][17][18] and in error mitigation [19]. The transformation of the Hamiltonian into tri-diagonal form with respect to some initial trial-state |v 1 , proceeds according to the Lanczos recursion as [15]: The matrix elements of the Hamiltonian in the Lanczosbasis are given by In classical applications, the eigenvalues of the truncated tri-diagonal Hamiltonian matrix in the Lanczos-basis |v i are computed numerically. Although the approximates converge relatively rapidly to the low lying states of the original Hamiltonian, the classical computation of the α i and β i are generally limited by the matrix dimension of the problem. For the quantum computing context we consider the Lanczos recursion in the formalism of Hamiltonian moments H n ≡ v 1 |H n |v 1 , with respect to an appropriate initial state, |v 1 (e.g. in the ground state sector of the Hilbert space). Some time ago, it was found that there exists a general expansion of the Lanczos matrix elements with respect to the Hamiltonian moments [13,14]: where z is a continuous positive parameter related to the recursion index [14,20] (z ≡ i/V for extensive systems of volume V ), and the quantities c n are the cumulants derived from the moments H n : The explicit general expressions for the Lanczos matrix elements α(z) and β(z) allow powerful results from orthogonal polynomial theory to be brought into play [21,22]. As a result, at all orders in the z-expansion for the Lanczos matrix elements, the ground-state energy of the Hamiltonian system can be expressed directly via an infinum theorem [14]: This rather innocuous looking relationship actually diagonalises the tri-diagonal Lanczos system in its expanded form, allowing one to determine approximates, E (inf) 0 , to the solution E 0 by truncating the z-expansion to some maximum moment/cumulant order n max . At first order in z, a general "infinum" approximate for the ground state energy involving cumulants to order n max = 4 was derived [23]: Identifying the first order cumulant as c 1 = H , the second term in the infinum estimate, involving cumulants c 2 → c 4 , provides a direct means of improving on the variational calculation over a given trial-state |φ trial ≡ |v 1 .
In principle, one can determine higher order infinum estimates [24], however, here we will show the analytic E (inf) 0 expression for n max = 4 already provides a convenient and powerful correction to the variational calculation represented explicitly by c 1 . Dynamics of the system are neatly encapsulated in the cumulants, and since the correction to the variational estimate of E 0 obtained from the infinum theorem corresponds to the diagonalised Hamiltonian in the tridiagonal Lanczos basis, derived approximates sum the associated (truncated) dynamical effects to all orders. The Lanczos expansion and associated infinum theorem has been applied to a number of homogenous many-body systems, from quantum magnetism [25] to lattice gauge theory [26], however, the computation of moments for non-homogenous systems generally scale poorly on classical resources as the system size grows. For systems of interest we seek to directly compute these quantities on a quantum processor.
We thus arrive at the quantum computed moments (QCM) approach for estimating the lowest energy of the problem Hamiltonian:

Obtain infinum approximate E
Before moving to implementation, we make a few general remarks. The QCM approach places more emphasis on the number of circuit runs on the quantum computer and associated classical computation to achieve a better energy estimate rather than increasing the circuit depth in the variational approach. The actual complexity of the quantum information processing task involved in computing moments relative to classical resources is dependent on the details of the trial-state, as we note further on. At first sight, the procedure for exponentiating the Hamiltonian appears to scale badly, however, the number of terms can be tightly controlled by creating Tensor Product Basis (TPB) sets. Another issue that we focus on in the implementation is how shot noise and device errors flow through the arithmetical operations from TPB set measurements, to moments, to cumulants. Although we have articulated the approach in the context of obtaining infinum estimates for the ground-state energy, we note the quantum computed moments may be used beyond this context. The approach here is quite distinct to recently proposed Lanczos-based methods [16][17][18][19] and approximates in the connected moment expansion form of the t-expansion [27][28][29]. The infinum estimate corresponds to the diagonalisation of a truncated cumulant expansion of the Lanczos tri-diagonal basis, as opposed to the usual approach of computing and diagonalising the truncated Lanczos basis itself.

IV. HAMILTONIAN PROBLEM, OPERATOR REDUCTION AND SCALING
To show how the algorithm performs in practice, we will consider a non-trivial example system from quantum magnetism. The quadratic Hamiltonian (density) for q qubits is given by: where the sum is over a problem graph defined by the vertices (qubits) i = 1...q, edges connecting qubits { i j }, and couplings J (s) ij along each edge (s = x, y, z). Here we consider nearest-neighbour 2D lattices with freeboundary conditions. The uniform coupling case J ij is the well known 2D Heisenberg model, for which the exact ground state has been extensively studied numerically.
We first detail the Hamiltonian exponentiation and the scaling of the effective number of Pauli strings required for measurement. Initially, one concatenates and compresses products of Pauli strings at each level of H n : where the [P (n) j ] are q-length Pauli strings resulting from the product reductions, and A (n) j are the resulting weights. Naive counting suggests the number of Pauli strings in the expressions corresponding to powers of the Hamiltonian increases exponentially with n. However, by exploiting the properties of the Pauli matrices and their commutation relations, the number of strings required for measurement can be drastically reduced by finding tensor product basis (TPB) sets [Q k ] of Pauli strings that mutually qubit-wise commute (QWC) [4]. Thus, we rewrite the operator H n in terms of the TPB sets as: Finding the optimal TPB sets [Q k (n)] can be mapped to a minimum clique cover problem for the equivalent graph and solved heuristically [30]. Here we have instead determined the sets [Q k (n)] via an identity-operator sorting algorithm. The final TPB sets [Q (n) k ] of Pauli strings to be measured depends on the underlying problem graph { i j }. To show this, we perform the reduction process for three types of graphs -linear, heavy-honeycomb, and square lattice -for systems defined on up to 36 qubits. In Figure 2(c) we plot the growth with the number of qubits, q, for naive term counting. The dramatic effect of the TPB grouping process on the scaling is evident -for a given q, the number of Pauli strings to be measured drops by several orders of magnitude with sublinear scaling in q. The qubit-wise measurements of each of the concatenated strings [P j ] from which the moments H n are assembled.

V. RESULTS
As a first application of the method we consider the case of the 2D Heisenberg model defined on a 2 × 3 lattice with uniform coupling J (s) ij = 1. The first few TPB sets produced by the grouping algorithm at order H 4 are shown in Figure 2(b). Following a simplified version of the VQE construction Fig 3(a), we define a trial-state in single parameter form, |φ trial (θ) , as shown in Fig 3(b). This choice of trial-state includes the antiferromagnetic Néel state at θ = π. Away from θ = π the full set of 2 q states are engaged (e.g. θ = 0.7π shown). While the model itself is defined on a 2D lattice, we meet the challenge of restricted qubit array connectivity by defining the trial-state over a 1D array of qubits.
Results shown in Fig 4 correspond to the QCM algorithm run on the IBM Quantum processor ibmq montreal -the qubits used are indicated on the device map. Comparison simulations were run on the Qiskit QASM simulator. We have plotted, as a function of the trial-state parameter θ, the moments H n , and associated cumulants c n , assembled from the QC measurements of the TPB sets {[Q (n) k ]} up to n max = 4. Quantum calculations were carried out using 5 × 1024 shots per expectation value. Note: these are raw results with no attempt at error mitigation or improved sampling [32][33][34][35]. Compared to the exact/simulation results (solid lines), the moments computed on the quantum computer system are surprisingly free of shot noise, with deviations largely due to the device errors. The cumulants have higher statistical noise, as expected given their composition in terms of the moments. In Fig 4(d) we plot the infinum estimates E (inf) 0 obtained from the device runs together with variational results on H θ and simulations carried out for different noise levels (zero to 8× device default error model). We make the following observations: (i) The infinum estimate significantly improves on the variational result for the same trial-state; (ii) Through all the post-classical manipulation of measured quantities to assemble H n and c n , the overall statistical noise in the final QCM infinum results appears to be not too much greater than the variational results, and certainly much less than their difference; (iii) The quality of the infinum estimate derived from the trial-state on the 1D qubit array persists for a range of values of the trial-state parameter either side of the variational minimum (θ = π); (iv) The simulations indicate the infinum estimate is more robust to device noise than the variational calcula- tion on H θ .
To test these observations we move on to larger and more complex instances of the 2D model. The 1D trialstate form |φ trial (θ) is retained, but the model is generalised to the case of random couplings, {J (x,y,z) ij }. We note another important feature of the QCM approachonce the Pauli string reduction and measurements have been carried out for a particular problem graph, one need not repeat when the couplings in the problem Hamiltonian are changed. In effect, one only needs to run the moments computation once on the quantum computer -the infinum estimates for an arbitrary large ensemble of random instances can be computed efficiently using classical resources post-facto by recycling the quantum computed moments output. For problem instances up to 4 × 4 we compute and compare with exact results, however, at 5 × 5 the Hilbert space dimension of the problem is O(10 7 ) and begins to challenge convenient classical computation. As a reference, we compare with the 2D Heisenberg model case with uniform couplings for which the ground-state is known numerically [36] (Note the infinite lattice limit value is E 0 = −2.676 in our qubit notation).
In Fig 5(a) we show the results for square lattices 3×3, 4×4, and 5×5 (8192 shots per data point) with the trialstate implemented on qubit chains as shown on the IBM Quantum device ibmq manhattan Fig 5(b). Random coupling instances correspond to choosing J (x,y,z) ij ∈ [0, 1) (three decimal places). Across the board, the QCM infinum results improve on the variational benchmark and are remarkably close to the exact results given the 1D restriction of |φ trial (θ) . The data is accompanied by zero-noise simulations, which clearly show that while the variational data points obtained with device noise consistently move away from the true ground-state energy, the QCM infinum estimates around the Néel point are remarkably inert and maintain the robustness to both noise and change in the trial-state as shown in the q = 6 results. In Fig 5(c) we plot the fidelity ratios with respect to the exact results. The difference between the QCM results and the variational benchmark is cleardevice and shot noise are well under control for the QCM infinum estimates. For the largest 5 × 5 instance, the fidelity ratio with respect to the exact value (uniform 2D Heisenberg model [36]) for the QCM infinum estimates is 91%. The recycling of the one-time quantum output also works well, despite the assembly process for H n and c n , providing an ensemble of results over random coupling instances -the ensemble distributions (10 3 instances) are shown in Fig 5(d).

VI. CONCLUSION
The Quantum Computed Moments approach presented here shifts the focus of the representation of prob-lem complexity from the trial-state to the quantities being measured on the quantum computer -Hamiltonian moments. We demonstrated the method on models of 2D quantum magnetism using IBM Quantum processors for instances up to 25 qubits. At order H 4 , the data suggests the classical pre-processing into TPB scales sub-linearly for these models. For our investigations a single-parameter trial-state was chosen, implemented over relatively shallow depth circuits on the devices at hand, and which encompassed the antiferromagnetic Néel state. Hamiltonian moments to fourth order were quantum computed and the infinum estimate obtained was found to provide a significant correction to the variational estimate. The infinum results were stable over a significant range of the variational parameter either side of the minimum H (Néel point). The 5 × 5 instances begin to Solid lines correspond to zero-noise simulations (Qiskit QASM simulator) for variational (blue) and infinum (orange) estimates. The exact groundstate energy is plotted as a green dashed line. For each lattice size the trial state |φ trial (θ) (Fig 3(b)) was encoded using a linear chain of qubits, as shown in the device map (b). (c) Fidelity ratio with respect to the exact results for uniform coupling model and the random coupling ensemble average (10 3 instances). (d) Frequency distribution over the random coupling ensemble for QCM infinum estimates and variational benchmark, compared to the exact results calculated for 3 × 3 and 4 × 4.
surpass convenient classical verification, however, our results compared well to those reported in the literature for the uniform coupling case [36]. Our trial-state was chosen as a simple benchmark for comparison purposes, rather than the quest for precision in the final result. Given the relative robustness of the infinum results to device noise we expect there is scope to improve the precision further through more carefully designed trial-states and/or error mitigation on the measured quantities [32,33]. On the question of relative computation workload on the quantum processor: for the benchmark trial-state used here the classical computation of the moments can be cast as an efficient sampling problem given the structure of the trial-state circuit [38]. More complex trial-states can be constructed, however, our focus has been on demonstrating the robustness of the approach -in particular that away from the Néel point at θ = π where the trial-state is more populated, the estimates clearly survive the effect of device errors and shot noise on the arithmetical processes. Finally, we demonstrate an important feature of the hybrid calculation -once the TPB measurements are carried out on the quantum computer, the infinum estimate for any other Hamiltonian of the same form can be computed entirely in the classical post-processing. Our results for random coupling instances are obtained by recycling a one-time set of quantum measurements over the single-parameter trial state. Even though the trial-state is more suited to uniform couplings, where we were able to perform exact calculations, the relative precision of the QCM infinum results consistently holds as per the uniform Heisenberg model case.
In this work we have focused on the practicalities of the quantum computation of Hamiltonian moments to produce estimates of the ground-state energy of a given problem. In introducing the moments based approach we have demonstrated the relative improvement possible over variational calculations for the same relatively simple trial-state. Additionally, higher order approximates can be obtained from the cumulant expansion. Further studies would include a systematic analysis of different trial-states on the precision of the estimates obtained and the potential to provide stable results within the quantum volume [37] constraints of state-of-the-art devices. In terms of scaling, we have shown that the TPB set reduction offers significant savings in both classical computing requirements in the pre-processing step, and the number of measurements on the quantum system. For problems with a large number of terms in the base-level Hamiltonian, perturbative approaches could be usefully applied to maintain a workable scaling in H n . Going beyond ground state energy problems, there is scope for determining the solution configuration through the Lanczos approach, as well as the application of similar moment based results for excited states [39], and the application to ZZ optimisation problems by systematic inclusion of mixing terms in the corresponding Hamiltonian form of the problem -thereby expanding the class of problems the approach could be applied to.