Quantum inverse iteration algorithm for programmable quantum simulators

We propose a quantum inverse iteration algorithm, which can be used to estimate ground state properties of a programmable quantum device. The method relies on the inverse power iteration technique, where the sequential application of the Hamiltonian inverse to an initial state prepares the approximate ground state. To apply the inverse Hamiltonian operation, we write it as a sum of unitary evolution operators using the Fourier approximation approach. This allows to reformulate the protocol as separate measurements for the overlap of initial and propagated wavefunction. The algorithm thus crucially depends on the ability to run Hamiltonian dynamics with an available quantum device, and can be used for analog quantum simulators. We benchmark the performance using paradigmatic examples of quantum chemistry, corresponding to molecular hydrogen and beryllium hydride. Finally, we show its use for studying the ground state properties of relevant material science models, which can be simulated with existing devices, considering an example of the Bose-Hubbard atomic simulator.


I. INTRODUCTION
Quantum computing offers the drastic speed up for certain computational problems, and has evolved as a unique direction in the theoretical information science [1]. However, the field of experimental quantum computing is yet at its infancy. The typical size of quantum chips for the reliable gate based quantum computation ranges from one to several tens of physical qubits, with the main limits posed by decoherence. Despite the imperfections, the algorithms of ever-increasing complexity were implemented on different platforms, with circuit depth exceeding a thousand gates [2,3], potentially allowing for quantum supremacy demonstration.
At the same time, the vox populi of quantum engineers says that while experimental setups are developed and mastered rapidly, the theorists in the field lag behind. Whereas by now textbook examples of quantum algorithms with exponential and quadratic speed up for factoring and search serve as a great motivation [1], the estimates of gate counts are daunting, making them distant goals for the future fault tolerant quantum computers [4]. The recent developments in this fast evolving field call for new short depth algorithms which can solve useful problems in the era of noisy intermediate scale quantum (NISQ) devices [5].
One of the most promising directions for quantum computation is the field of quantum chemistry and materials [4,6]. Targeting the access to the ground state properties of molecules and strongly correlated matter, it can offer huge gain for various technological applications, for instance helping to find a catalyst for the nitrogen fixation [7]. To date, different quantum theoretical protocols were developed, and several proof-of-principle experiments on various platforms were performed in the simplest cases. Examples include simulation of molecular hydrogen with linear optical setup [8], superconducting circuits [9][10][11][12], and trapped ions [13]. Finally, the variational simulation for larger molecules (LiH and BeH 2 ) were reported recently [10]. From the material science perspective, the use of cold atom quantum simulators has shown great promise, where simulations of Fermi-Hubbard lattice dynamics [14], large scale quantum Rydberg chain [15] and Ising model [16], and two-dimensional many-body localization [17] have been performed. However, in the latter cases the analog approach to simulation is taken, given an access to unitary dynamics, while precluding the study of ground state properties.
To access the ground state properties of a quantum chemical Hamiltonian, several routines can be used (see Refs. [18,19] for a review). First option corresponds to the quantum phase estimation algorithm (PEA) [20], which exploits the unitary dynamics of the system controlled by register qubits. Although this algorithm is efficient, giving logarithmic error and polynomial gate scaling, its implementation requires substantial circuit depth for currently available circuits [9]. Moreover, the controlled type of operations require the digitization of the circuit, thus complicating the use of analog quantum circuits for PEA. Another approach is the adiabatic quantum computing, which was already applied to quantum chemistry problems [21]. However, the required adiabaticity of dynamics typically results in the effectively long circuit depth. Finally, an alternative route to quantum chemistry and materials is offered by hybrid-classical variational approaches which were proposed recently [22]. They rely on term-by-term energy measurement for the prepared trial quantum state (ansatz) with consequent classical optimization, and are referred to as Variational Quantum Eigensolvers (VQE) [18,19,23]. It can use a chemically inspired ansatz [23], Hamiltonian variational ansatz [24], or rely on the variational imaginary time evolution [25]. In this case the depth of the quantum circuit is greatly reduced, though at the expense of increased number of measurements, being favourable strategy for NISQ devices. For VQE the number of variational parameters scales as O[(3N ) k ], where N is a number of qubits and k represents an approximation order [23].

arXiv:1901.09988v2 [quant-ph] 30 Sep 2019
While for quantum chemistry applications k = 2 suffices to give useful results, these approaches are yet to be tested for larger system sizes, where multi-variable optimization may raise the problems for genuine ground state estimation [26].
In the following, we propose the quantum inverse iteration algorithm for the estimation of the ground state energy (GSE) of a quantum system. It is inspired by the classical inverse power iteration algorithm for finding the dominant eigenstate of the matrix, where the computationally demanding part of matrix inversion and multiplication is performed with a quantum circuit. Previously, a direct iteration approach was considered as a general purpose quantum algorithm [27], aiming for a large scale fault tolerant implementation. Here, we present the protocol of the hybrid quantum-classical nature. It relies on performing quantum evolution for different propagation times and classical post-processing of the measured observables. The approach is applied to quantum chemistry examples (H 2 and BeH 2 molecules), showing favourable scaling with system parameters. Finally, when applied to the Bose-Hubbard quantum simulator, it allows to study its ground state properties, thus showing promise as a protocol for near-term quantum devices.

II. PROTOCOL
We start by considering a generic interacting system, which can be described by the second quantized Hamiltonian. It can be written as sum of two-body and four-body partsĤ whereâ † i (â i ) can correspond to the fermionic or bosonic creation (annihilation) operators, and cover broad range of models. In the fermionic case, the Hamiltonian (1) can describe the full configuration interaction problems in quantum chemistry, with operatorâ j corresponding to molecular orbital j. Using the existing mappings between fermionic and spin-1/2 systems, one can rewrite (1) in the form of a local HamiltonianĤ for interacting qubits which involves strings of Pauli operators. The task is then to find the lowest eigenvalue of the large matrixĤ, corresponding to GSE.

A. Inverse iteration
We propose a procedure which can be seen as a quantum version of the inverse power iteration algorithm for finding the dominant eigenvalue of the matrix, represented by the inverse of the Hamiltonian matrixĤ −1 , which is treated as a dimensionless matrix in this section. Given thatĤ is invertible, the order of eigenvalues is reversed and, with the appropriate shift of the diagonal to make eigenvalues positive, the power iteration allows to find GSE. Namely, starting with an initial state |ψ 0 which has non-zero overlap with the sought ground state |ψ gs , by repetitive application of the inverted matrix one can prepare (unnormalized) state, | ψ k = (Ĥ −1 ) k |ψ 0 , such that | ψ gs |ψ k | 2 < for sufficiently large number of iterations k ≥ K [28], where |ψ k = | ψ k / | ψ k (Fig. 1A). We note that generally this method has favourable logarithmic complexity in the iteration depth, being K = log[ sin −2 (θ 0 )/(λ 1 − λ n )]/[2 log(λ 2 /λ 1 )], where λ 1 , λ 2 , and λ n correspond to dominant, sub-dominant, and smallest eigenvalue of H −1 . Here sin 2 θ 0 parametrizes the overlap between |ψ 0 and |ψ gs , marking that convergence of the procedure depends on the initial guess, and generally can be made nonzero taking |ψ 0 as a random state. While classical power iteration methods generally have good convergence in the number of iterations, the main caveat comes from the complexity scaling with the system size. The requirement for K matrix multiplications leads to O[K2 2N ] operations (for dense matrices), yielding exponential scaling. Even worse situation is for the inverse matrix algorithm, where an overhead comes from theĤ −1 calculation, requiring extra O[2 N ] operations.

B. Fourier approximation
In the following we show that we can exploit the iterative procedure with logarithmic iteration depth in , while providing exponential speed up for the inverse Hamiltonian multiplication process. The latter comes from the approximation theory [29], observing that the inverse can be represented as an integral x −1 = +∞ 0 exp(−xy)dy, which by applying the trapezoidal rule can be written as a sparse sum of exponents. For quantum systems a similar idea was proposed in Ref. [30], where Fourier approximation of the Hamiltonian inverse was presented as a double integral of the unitary propagator. This was further used to design an efficient solver for the quantum linear equation system problem [30,31]. Here we extend the Fourier approximation to the k-th power of the inverse (k ≥ 1), which formally readŝ and N k is a normalization factor. The integral can be then discretized aŝ , where ∆ y,z correspond to the discretization steps for integration variables, and M y,z represent cutoffs for inte-quantum inverse iteration method where we have rewritten the double summation in Eq.
(3) using the superindex (j y , j z ), φ k, = (j y ∆ y )(j z ∆ z ) is a phase of evolution for parameters chosen to discretize k-th inverse, and L k = M y (2M z + 1). Here c k, represent purely imaginary coefficients for the series, and |c k, | define corresponding weights. The required size and number of discretization steps ∆ y,z and M y,z depends on and κ (see Methods, section A, for the details). Importantly, they set the maximal evolution phase φ max = (M y ∆ y )(M z ∆ z ), which serves as an equivalent of the total gate count for analog quantum simulation.
As we need to prepare the approximate ground state by applying generally non-unitary operatorĤ −K to the initial state |ψ 0 , we shall either introduce an ancillary register to perform it, or properly account for the normalization of the resulting wavefunction. The former option is an excellent strategy for the future fault-tolerant devices, and has beneficial scaling (Methods, sec. A), while the latter is more suitable for programmable quantum simulators. In the following we present the strategy which can be applied to estimate the ground state properties by sequential evaluation of terms in the series. Similarly to VQE approaches, this relies on performing large number of measurements, and thus adds an extra complexity as compared to the generic implementation of the inverse operator. At the same time, term-by-term readout offers better resilience to errors where even imperfect procedure can yield reasonable GSE estimate for quantum simulators.

C. Sequential energy estimation
Our final goal is to estimate observables of the system, provided that an approximate ground state is prepared. For any operatorÂ it can be retrieved from the measurement A = ψ gs |Â|ψ gs / ψ gs |ψ gs , where the normalization is accounted for explicitly. In particular, we are interested in the calculation of the ground state energy λ gs ≈ λ k , choosing the operatorÂ asĤ. This amounts to measurement of Hamiltonian expectation value for |ψ k =Ĥ −k |ψ 0 in the form We proceed by considering each propagated wavefunction separately, such that λ k can be related to wavefunction overlaps (see flowchart in Fig. 1). This is motivated by the Hamiltonian averaging procedure [32] used in VQE to reduce the circuit depth at the expense of larger number of sequential measurements. Using Fourier expansion of the inverse Hamiltonian (4), the estimated energy reads Note that expression (6) now includes overlaps between initial and evolved wavefunction for the fixed phase, which shall be calculated separately for the numerator ("energy") and denominator ("norm"). Finally, we note that the wavefunction overlap can be inferred using different approaches. One option corresponds to using the SWAP test [33,34], which represents a common digital strategy and requires system doubling. This can be conveniently realized in some near-term setups, being successfully demonstrated for cold atom lattices by manybody interferometry of two copies of a quantum state [35]. Additionally, in the Methods, section B, we describe an alternative approach which does not require extra qubits and relies on the measurements of observables, once the reference state for the system is chosen.

III. RESULTS: QUANTUM CHEMISTRY APPLICATIONS
In the previous sections we described the general algorithm and discussed its key properties, namely the scaling and sequential operation. To show its use for the ground state estimation and characterize the required resources for realistic problems, we apply it to quantum chemistry. There, new strategies are much required due to rapidly growing complexity with the number of qubits (orbitals) N .

A. Molecular hydrogen
We start with by now the standard example of molecular hydrogen, H 2 . As a test task we consider the spinful case. This allows to examine the protocol for a system of higher complexity (N = 4), comparable to lithium hydrate four-qubit simulation considered in Ref. [10]. The details of mapping of quantum chemical structure into qubits are presented in the Methods, section C. In the following we work with four-qubit molecular hydrogen HamiltonianĤ H2 , with all eigenenergies shifted to positive values. Starting from the Hartree-Fock (HF) energy λ 0 , the task is to estimate GSE λ gs , using the protocol described in the preceding section. This shall be done within the chemical precision , which is equal to = 0.0016 Hartree, and thus defines the relevant cut-off for the iteration procedure.
We start by benchmarking the inverse power procedure in its general form, and define how many iteration steps one needs to come close to the ground state. For this, we first perform the inverse Hamiltonian iteration in the ideal setting, assuming that an exact inverse is known. Then, we compare it to the quantum inverse iteration, which uses the Fourier approximation (4). GSE is estimated using the measurement of propagated and initial wavefunction overlaps. To quantify the performance two characteristics are employed. The first, and the most natural one, corresponds to the difference between estimated energy value λ k and true GSE λ gs , being ∆λ/J ≡ (λ k − λ gs )/J. It allows to observe the convergence and provides an indication of how well the procedure works for a given system. The second quantity corresponds to the trace distance between an idealized inverse iteration matrixĤ −k and its approximationĤ −k a , defined as a half of trace norm for the difference of two matrices. It reveals the actual success of mimicking the ideal inverse in full generality. At the same time, this is the quantity which cannot be straightforwardly observed in the experiment, and only serves for the analysis.
The results of the inverse power iteration for molecular hydrogen HamiltonianĤ H2 are shown in Fig. 2A as a function of iteration step k. The ideal version of inverse iteration is plotted in red and reveals exponential convergence to GSE. The chemical precision is achieved already at the second step of the iteration, as depicted by the blue shaded area starting at ∆λ = 1.6 × 10 −3 J. The idealized case is then compared to the quantum inverse iteration procedure with combined measurement of wavefunction overlaps as stated in Eq. (6). Here, we assumed that the genuine unitary evolution with Hamil-tonianĤ H2 is run in the analog simulation fashion. The case of digital evolution with associated Trotterization technique and its benchmarking is considered in the Supplemental Material, where we also present the circuit scheme for digital evolution. The approximation was performed using equal number of steps M z = M y = 30, and the discretization values ∆ z = ∆ y J were adjusted to match the maximal propagation phases of the propagation phase is taken to be dimensionless by absorbing energy unit prefactor J from the Hamiltonian). The corresponding curves show the improvement of the quantum power iteration estimation for increasing number of iteration steps. The convergence rate also depends on the maximal phase of the propagation. For small phases (top curves in Fig. 2A), the initial estimator does not give successful convergence, but comes closer to GSE for large k. As the propagation phase grows, the approximation λ (a) k , starts to resemble the idealized iteration procedure. However, this only happens up to a certain value of k past which the approximate energy grows, thus deviating from the ideal solution. From the point of view of process fidelity, the trace distance Tr[Ĥ −k ,Ĥ −k a ] between ideal and approximate inverse operators increases monotonically with k (Fig. 2B). The increase of φ max allows to reduce Tr[Ĥ −k ,Ĥ −k a ] at each k. The performance of the quantum inverse iteration procedure is further analyzed in Fig. 2C,D where energy distance to ground state and trace distance are shown as a function φ max for several fixed iteration steps (k = 2, 4, 7). The calculations were performed accounting for two different ways of arranging the phase. First, the approximation grid was fixed setting M y = M z = 30 while changing ∆ z = ∆ y J (solid curves in Fig. 2C,D). In the second case the fixed step size ∆ z = ∆ y J = 0.05 was combined with the increment of M z,y (dashed curves in Fig. 2C,D). For both energy distance (Fig. 2C) and trace distance ( Fig. 2D) we observe no difference between two approximation procedures, but clear indication of the importance of maximal propagation phase (time). For ∆λ one sees a non-monotonic dependence on φ max which starts with a decrease of the energy difference for increasing maximal phase (φ max /2π < 0.4). At larger phases the dependence experiences pronounced dips (note the log scale), which are more visible for many iterations. Overall the difference remains well-within chemical precision and experiences saturation. When the trace distance is considered, one sees that success of the approximation monotonically improves with φ max . At the same time, for fixed approximation parameters {M z,y , ∆ z,y } it is more difficult to represent inverse iteration operator faithfully, in-line with scaling analysis discussed in the Methods. Finally, the comparison of results in Fig. 2C and D allows to suggest that non-monotonicity in the spectroscopic sig-natures can come from the particular structure of the Hamiltonian and the initial state, where certain phases might be preferable (i.e. not all elements of the Hamiltonian matrix contribute equally to the inverse iteration procedure).
To decide on the optimal way to approximate the inverse, we consider different discretization steps for y and z auxiliary variables, characterized by the skewness parameter ∆ y J/∆ z . The calculation is done for M z = M y = 30 with the maximal phase fixed to φ max /2π = 0.92. The results are shown in Fig. 2E,F as a function of skew. The energy difference parameter shows that for approximating the inverse for small iteration numbers (k = 2 curve in Fig. 2E) larger skew factors are preferable, with z variable requiring finer approximation. However, for increased iteration number the optimum flows to ∆ y J/∆ z ∼ 1 values, suggesting close-to-equal spacing can work well for varied k. Examining the trace distance, we see that in unbiased setting the skew ratio of ∆ y J/∆ z ∼ 2 is preferable.
Finally, the very important issue to address is an influence of noise on the operation of quantum inverse iteration protocol. For this we have performed the analysis including relevant dephasing processes which influence the estimate for overlaps (see details in the Supplemen-  tal Material). Although noise makes the estimation of energies at large iteration step k less reliable, it is possible to estimate the energy within chemical precision using simple noise mitigation techniques.

B. Beryllium hydride
To test the scalability of the approach, we consider a molecule of bigger size which requires larger Hilbert space simulation. For this, we choose to simulate beryllium hydride (BeH 2 ) in the full spinful version using N = 8 qubits (see Methods, section D, for the details).
We proceed in the same manner as for H 2 molecule, and quantify the operation of quantum inverse iteration procedure for BeH 2 . The approximation parameters were chosen as ∆ y J = ∆ z = 0.05, with the number of discretization points M y,z adjusted accordingly to maintain maximal propagation phase. The results of the simulation are shown in Fig. 3. The first plot (Fig. 3A) shows that ideal iteration works well for the beryllium hydrate, with chemically precise GSE obtained already at k = 1 iteration step. The Fourier approximation for the inverse at small phases does not reach required accuracy, while for increased iteration step number and φ max /2π > 1 ground state estimate can be attained. Fig. 3B shows this behavior as a function of phase for several representative k's, and yields the same conclusion. The increase of the required propagation phase is attributed to the increased condition number for BeH 2 Hamiltonian matrix, being 39.2 as compared to 3.38 for H 2 .

C. Bose-Hubbard simulator
Finally, to provide an example where quantum inverse iteration algorithm can be largely beneficial, we consider the bosonic version of the Hubbard model (see sketch in Fig. 4A). The corresponding system Hamiltonian readŝ whereâ † i (â i ) corresponds to the bosonic creation (annihilation) operator at lattice site i. Here,n i =â † iâ i corresponds to the number operator. The first term in the Hamiltonian (7) describes the tunneling of bosons with rate J, which leads to their delocalization at the lattice. For simplicity we consider a one-dimensional lattice, and that tunneling happens only between neighbouring sites i, j . The second term in (7) denotes the contact interaction between bosons with strength U . The last term corresponds to the chemical potential µ. Importantly, the Bose-Hubbard model corresponds to the paradigmatic example of a hard material science problem [36], and received considerable attention from both theoretical [37,38] and experimental perspective [39]. In particular, the model was shown to be easily solvable in the so-called Mott insulating regime where U J and ground state corresponds to the product state of one atom per site, |ψ Mott = i |1 i . However, going into the superfluid regime, where J becomes comparable to U , the ground state is entangled, and it is generally difficult to study the low energy properties of this system.
We show that quantum inverse iteration algorithm can be applied to study the ground state properties of the Bose-Hubbard Hamiltonian in the wide range of parameters. For this, we consider an initial state corresponding to the Mott state, |ψ 0 = |ψ Mott , and perform the iteration as detailed in the "Protocol" section. The numerical results for the energy deviation are shown in Fig. 4B, where ∆λ/U ≡ |λ k − λ gs |/U is plotted as a function of iteration step k. We consider N = 5 lattice with µ/U = 0.5, different values of tunneling rate J, and similarly to H 2 and BeH 2 examples the energies were adjusted by trivial overall shift E 0 to ensure positive eigenvalues. The Fourier approximation is performed using M y = M z = 40 and ∆ z = ∆ y U = 0.075. While for relatively small tunneling J = 0.01U the product state is a good ground state approximation, thus giving small energy deviation (lowest curve in Fig. 4B), for larger J ∼ 0.1U the system enters a superfluid phase (critical value for µ = 0.5U approximately corresponds to J crit /U ≈ 0.13 [37]). Despite the qualitatively different initial state, the algorithm allows to distill correct energy properties even for large J.
Going beyond the ground state energy estimation, quantum inverse iteration can be also used to measure the correlations in the ground state of the model. We considered the correlation function of the form ψ k |â † c+râc |ψ k , where |ψ k corresponds to the propagated state, and overall procedure is similar to one described in Eq. (6). We fix c = 3 to be the central site in the lattice, and look for intra and intersite correlations with r = 0, 1, 2. The results are shown in Fig. 4C,D,E, where in the Mott insulator regime correlation fall-off rapidly (C), while in the superfluid case the quantum inverse iteration method allows to restore genuine correlations (blue dots and curve) at increased k.
Importantly, the described Bose-Hubbard simulation can be performed experimentally using the available controllable devices [35,39]. First, it allows to run unitary dynamics for sufficiently long times and large particle numbers. Second, the many-body interference technique enables the convenient measurement of the wavefunction overlap, where part of the system is evolved and later is interfered with the initial copy [35]. We note that our approach also bares similarities with recently introduced Quantum Virtual Cooling [40]. This study demonstrated for reducing the temperature in half, while the inverse iteration offers the access to close-to-zero temperature. Finally, the intriguing possibility is an application of quantum inverse iteration to 2D models [17] and Fermi-Hubbard model [14], which shall be possible with improved interference techniques [40].

IV. DISCUSSION AND CONCLUSIONS
We have presented the algorithm for the ground state energy (GSE) estimation of a quantum Hamiltonian. It is based on the iterative application of the Hamiltonian inverse to the initial state, and can be represented as a sum of unitary evolution operators. Targeting near-term quantum simulators, we described the protocol as a separate estimation of GSE contributions from the wavefunction overlap measurements. Then, the results from the quantum dynamical simulation are post-processed classically, and provide energy estimate for each iteration step.
The algorithm was applied to several quantum chemistry examples, being molecular hydrogen and beryllium hydrate. Using the four-qubit H 2 simulator, we benchmarked the performance of iteration and inverse approximation, showing that the most valuable resource for GSE estimation is a maximally available time for unitary evolution. Both digital and noisy operation was considered, and found to be sufficient for a GSE calculation with chemical accuracy.
As an outlook, we highlight that the approach can be beneficial for analog quantum simulators such as cold atoms lattices [17], Rydberg atom simulators [15], trapped ions [16], and superconducting devices [41]. For instance, the analog-type fermionic quantum chemistry simulator [42] would be much valued for the task. Future applications also include material science problems, with the main target being Fermi-Hubbard model [14]. For instance, we note that recently proposed approach of Quantum Virtual Cooling [40], which was experimentally applied to Bose-Hubbard model, has similar iterative structure and requires interferometric measurements. This poses the question of connection between the measurement-based cooling scheme and the dynamic protocol described in the current study.

A. Scaling and fault-tolerant implementation
In this section we consider the scaling for the quantum inverse iteration algorithm. For k = 1 it was shown in Ref. [M1] that the inverse Hamiltonian can be approximated up to an error setting the discretization steps to ∆y = Θ( / log(κ/ )), ∆z = Θ(1/κ log(κ/ )), and summing up to My = Θ(log(κ/ )κ/ ), Mz = Θ(κ log(κ/ )) (κ is a coordination number). The maximal evolution phase then scales as φmax = (My∆y)(Mz∆z) = O(κ log(κ/ )), and is an equivalent of the total gate count for analog quantum simulation. Generalizing the result to k-th inverse iteration, the upper limit on the maximal required phase shall be multiplied by K, where truncation of inverse iteration introduces an additional error. In the study we considered particular examples, and quantified the validity of Fourier approximation as a function of {My, ∆y, Mz, ∆z} parameters.
Considering the fault-tolerant implementation, the approximate ground state can be prepared by applying generally non-unitary operatorĤ −K to the initial state |ψ0 using an ancillary qubit register. One possible option here is the amplitude amplification approach [M2] which addresses the task of implementation of the sum of unitary operators, of the same type as the one in Eq. (4). Moreover, since we also require simulation of Hamiltonian dynamics for exp(−iφ k, Ĥ ), which may be not accessible in analog-type simulation, the subsequent use of Hamiltonian simulation [M3] or qubitization [M4] methods would lead to favourable resource scaling. The algorithm will require O(log(L) log(cφmax/ )/ log log(cφmax/ )) auxiliary qubits (c ≡ |c |) and same order of controlled unitaries. This scaling can be compared to the iterative modification of the quantum phase estimation procedure (IPEA), based on a small fixed register [M5] or a single auxiliary qubit [M6]. The latest represents conceptually the closest algorithm to the one described in the paper, and thus will serve as benchmark. The complexity of IPEA was discussed in Ref. [M6], showing the requirement of O[log( ) log(log( )/ )] phase iterations to approach an error of = 2 −m (energy is rescaled such that Ĥ < 2π, and m is the number of relevant bits of precision, typically limited to < 20 for quantum chemistry applications). Each k-th IPEA step then requires implementation of the c-U k operation, defined as implementation of (e −iĤ ) k , controlled on the register qubit. This leads to O[N 4 log( ) log(log( )/ )] gate count, comparable to the inverse iteration procedure described above.
Given the favorable scaling, the cost of the general purpose quantum inverse iteration algorithm can be small for future large scale quantum devices. However, generally there is no simple procedure to perform c-U operation, and it requires decomposition into a set of universal gates or multilayer swap technique [M7]. This enlarges the actual circuit depth (while being polynomial), and precludes the implementation of U = exp(−iφĤ) in analog fashion. Therefore, we target programmable devices with possible analog-type implementation and use sequential estimation strategy described in the main text.

B. Overlap measurement
To measure the overlap between the evolved and initial wavefunction, we propose to exploit a single eigenstate |ψR of the system as a reference, and measure the overlap with respect to its energy λR (usually set to zero). This nicely fits the task of GSE estimation for the fermionic Hamiltonian, as its Hilbert space includes a vacuum state with no fermions present (unless space reduction procedure was performed). Similar technique was used for extracting spectroscopic signatures of photon localization [M8], and the same approach was applied for the measurement of the density of states for the many-body system from the random state evolution [M9]. The main steps for the measurement are as follows. The task is formulated as finding ψ0|ψ0(t) , where |ψ0 is the initial state (typically corresponding to the Hartree-Fock solution). The state |ψj(t) =Û(t)|ψj is a time-propagated state with some unitaryÛ defined by the expansion. The HF state can be prepared from the reference |ψR (vacuum or other product state) using the product of local operators, and we note that these states are orthogonal. Then, the overlap probability is measured as an expectation value of the operatorM 0 = |ψ0 ψ0| for time-evolved wavefunction, which reads Next, the superposition of the vacuum and initial state shall be prepared as |ψ+ = (|ψR + |ψ0 )/ √ 2 and evolved to |ψ+(t) . Its overlap probability is measured as an expectation value ofM+ = |ψ+ ψ+| operator. This can be written as and provides an information about real and imaginary parts of O0,t. The same procedure is performed for measuring an expectation value of the operatorMi = |ψi ψi|, where |ψi = (|ψR + i|ψ0 )/ √ 2. An additional information is gained with and both real and imaginary part of O0,t can be found for the known reference λR from the system of Eqs. Note that we are mostly interested in the real part of the sought overlap Re{ ψ0|ψ(t) }, as both the "norm" and "energy" terms of λ k are real, and imaginary parts of the overlap cancel out [M10]. Thus it is also possible to deduce the real part indirectly as and its sign can be inferred from the measurement in Eq. (M4). We note that in principle these are two related ways to estimate the real part of the overlap, corresponding to Eqs. (M4) and (M6). While being equivalent in the noiseless case, the effects of decoherence make two approaches distinct (see Supplemental Material for the details). For the practical purposes we thus refer to the overlap estimate in (M4) as direct estimation approach and refer to (M6) as the indirect estimation. Physically, the measurement procedure resembles the Belltype measurement, where in the simple case of two-qubits a single CNOT operation and the Hadamard gate are required. For the larger system the measurement is generalized to GHZ-type, and consequently requires increasing number of two qubit operators, which depends on how much an initial HF state is different from the reference state.
Finally, we note that direct measurement is possible in atomic setups, where many-body interferometry is applied to two copies. This can be used in the analog circuits where system size doubling plays lesser role, and can compensate the absence of fully addressable individual gate operation.

C. Molecular hydrogen Hamiltonian
The fermionic HamiltonianĤH 2 is first written in the form (1), where coefficients vij and V ijkl are calculated by conventional quantum chemistry methods. Here, we exploited the OpenFermion package for Python [M11], which allows to extract the interfermionic interactions for four Gaussian orbitals fit via STO-3G method and perform the fermions-toqubits transformation. For the small N = 4 system we have chosen to use the Jordan-Wigner transformation, although other options may be used as the system size increases. Specifically, we consider the bond length for H2 to be 0.7414 (measured inÅngström) and consider full excitation space. For concreteness, we provide the full form for the Hamiltonian, beinĝ where Xj, Yj, Zj denote Pauli matrices for qubit j. Hartree-Fock (HF) solution for the problem is given by the approximate ground state |ψ0 = (↓, ↓, ↑, ↑) T , and associated HF energy is −1.116684J. As required by the Fourier approximation approach, the reference energy is then shifted towards positive values by adding constant term equal to E0/J = 2., and we refer to the shifted Hamiltonian asĤH 2 in the main text. Its HF energy is λ0 = 0.883316J after the shift. The task is then to estimate the ground state energy λgs of the HamiltonianĤH 2 , achieved by preparation of approximate ground state |ψ k . This shall be done within the chemical precision , which is equal to = 0.0016 Hartree, and thus defines the relevant cut-off for the iteration procedure.

D. Beryllium hydride Hamiltonian
The molecular data structure of beryllium hydride (BeH2) was generated using Psi4 quantum chemistry package [M12] considering equal Be−H distances equal to 1.33Ångström. While generically described by six spin orbitals, we set lowest and second excited orbital to be occupied, and set multiplicity of unity, such that the ground state energy lies close to the full configurational space solution (STO-3G basis). The fermionic Hamiltonian is then obtained using OpenFermion package, and as in the case of H2 the Jordan-Wigner transformation was used to rewrite it in the qubit form [M12]. The problem then can be solved using N = 8 qubits. The GSE from the exact diagonalization of original BeH2 Hamiltonian reads −1.806750 Hartree, and analogously to the molecular hydrogen the Hamiltonian matrix is shifted by constant energy term of 2 Hartree. The product state corresponding to Hartree-Fock solution reads |ψ0 = (↓, ↓, ↑, ↑, ↑, ↑, ↑, ↑) T , with associated energy for the shifted Hamiltonian being λ0/J = 0.203323.

ACKNOWLEDGMENTS
I would like to thank Anders S. Sørensen for useful discussions and suggestions, and Bart Olsthoorn for technical advise and reading the manuscript.

V. COMPETING INTERESTS
The author declares that there are no competing interests.

VI. AUTHOR CONTRIBUTIONS
O.K. proposed the original idea, preformed the calculations, analyzed the results, and written the manuscript.

VII. FUNDING
The work was supported by the Government of the Russian Federation through the ITMO Fellowship and Professorship Program.

VIII. DATA AVAILABILITY
The authors declare that the data supporting the findings of this study are available within the paper and its supplementary information file.

Supplemental Material: Quantum inverse iteration algorithm for programmable quantum simulators by Oleksandr Kyriienko
In this Supplemental Material we provide the detailed analysis of the performance of quantum inverse iteration algorithm for the case of digitized evolution, show the relevant digital operation scheme, and analyse the influence of noise.

IX. TROTTERIZATION
As shown by the analysis in the main text (Sec. III.A), the success of the ground state estimation protocol depends on the implementation of unitary operatorsÛ (φ) for various values of the phase. Thus it largely favors the analog-type simulation, where Hamiltonian evolution can be seen as a resource for ground state energy estimation. When this is not available, the corresponding unitary operator can be implemented approximately using the digital or Floquet strategies [1]. As the simplest approach for systems with limited resources we choose standard technique of Trotter expansion. In particular, we use the second-order expansion [2,3], where a single Trotter step unitary readŝ (1) where N Tr denotes the total number of Trotter steps.
Here, the Hamiltonian is considered as a sum of termŝ H = M mĤ m which can be implemented separately, and M is the total number of terms. For the molecular hydrogen Hamiltonian it counts fourteen Pauli terms (see Methods, sec. C, of the main text), M = 14. However, the first ten and the last four terms are mutually commuting, so effectively Trotterization procedure relies only on the application of two non-commuting partsĤ Z s andĤ strings . We note that an exact realization of corresponding unitaries is platform-dependent, and thus shall be implemented specifically to the hardware capabilities. One potential implementation is offered in the end of this section.
The performance of the quantum inverse iteration with approximate Trotterized unitary dynamics is shown in Fig. S1. First, we plot the difference between estimated value and ideal ground state ∆λ as a function of Trotter step number, considering fixed iteration step k = 4 and several values of maximal propagation phase (Fig. S1A). We see that for small φ max the estimate lies outside of chemical precision, as governed by Fourier approximation error for the inverse, while the dependence on Trotter step number is weak and shows quick convergence. For increasing phase (curves 0.43-0.92) the estimate improves, and even two Trotter steps may be sufficient, favoring NISQ device operation. At the same time, as expected for the Trotterization procedure, larger phase (evolution time) requires finer procedure, and shows stronger N Tr dependence. The discrepancy between analog and digital unitary dynamics can be seen in Fig. S1B, where the energy difference for both is plotted for various N Tr as a function of phase. While generally the increase of Trotter step number leads to convergence between the two (see curve 15 staying close to zero at all phases), for small N Tr the difference is a non-monotonous function of phase. In particular, at certain point the difference between λ (a) k and λ (a;Tr) k shrinks to zero, which can be We chosen the reference state to be |ψ R = (↓, ↓, ↓, ↓) T .
To account for the noisy operation we exploit the wavefunction Monte-Carlo (WFMC) approach well-known in quantum optics [7]. Contrary to a mixed state description [8,9] this allows to operate in the original vector space. Importantly, it relates directly to an experimental workflow, and provides extra intuition for running the algorithm in noisy setting. In a somewhat modified form, it was already applied in Refs. [10,11]. WFMC relies on calculating the evolution of the system subjected to noise in the form of quantum jump operatorĈ j for qubit j. Considering the uncorrelated dephasing processes for all qubits, we define jump operators asĈ j = √ γZ j , where γ is a dephasing rate. The expectation values are then measured over an ensemble of trajectories with different noise realization, and each probability overlap |O 0,q | 2 , |O +,q | 2 , |O i,q | 2 is measured separately in the numerical procedure, mimicking the experiment. We consider the concrete example of molecular hydrogen quantum inverse iteration with the Fourier approximation grid fixed to M y = M z = 5 and ∆ y J = ∆ z = 0.5, with φ max /2π ≈ 1. To perform energy estimation we consider explicitly all terms in Eq. (6) of the main text and observe that in many cases the propagation phases coincide. In the following the unique phase difference values are chosen, and we are interested only in their absolute values, as ±δφ q evolution yields equal real parts for the overlap for systems which respect time-reversal symmetry. Finally, δφ = 0 is accounted for trivially. This brings us to just 35 phase difference values to be used for quantum evolution. Next, the classical post-processing is performed, where measured noisy values for overlaps O q and O H q are multiplied by coefficients p k,q and summed together, taking different values of iteration step. An important information at this step is the prefactor in front of each overlap, as it gives a weight in the total estimate. Note that in the case of noisy operation this changes the success of the estimate, as it favors large weight for small δφ evolution and small weight for large δφ's.
In Fig. S3A we show the histogram of weights, taken as w k,q = ( i |p k,qi |)/( q,i |p k,qi |) for different iteration step number k. Here, i summation goes over all coefficients for the same δφ q , and each weight is normalized by the total sum of prefactors at given k. The histogram reveals that for small k the main contribution comes from small-to-intermediate size phase differences δφ/2π < 0.5, while higher lying weights are nearly negligible. As iteration step k grows, the estimates rely more on the overlaps for larger δφ (longer propagation times).
WFMC simulation was performed using 5000 trajectories and considering different noise rates γ ranging from 10 −4 J to 0.1J. Importantly, this is also compared to the average interaction constant for the Hamiltonian ξ = 0.1224J, taken as a norm ofĤ H2 . The case of γ ∼ξ then corresponds to the highly noisy limit, and genuine operation typically requires γ/ξ 1. The noise rate can effectively be changed (increased) in the simulation. As ξ is generally tunable, it sets the physical timescale for the operation. Simultaneously rescaling the physical couplings {ξ m } by changing J and physical evolution time, one can perform calculations for the same propagation phase δφ but with effectively different noise influence.
The examples of the measured overlap values are shown in Fig. S3B,C for phase differences δφ/2π = 0.68 (B) and δφ/2π = 1.4 (C). They were obtained using two approaches, which we refer as direct and indirect estimation. As explained in the Methods, sec. B, the direct estimation refers to measurement of Re{O 0,t } using Eq. (M4) and the indirect approach uses Eq. (M5). The corresponding result for the direct measurement is shown in blue dots, and the indirect corresponds to red dots. In both cases, the measured estimate coincides with the noiseless value at γ → 0, shown by the blue line (see Fig. (S3)B,C).
Next, in the spirit of the error mitigation technique [8,9], we choose an optimal point γ min where the dephasing rate is minimized (overall coupling strength J is maximized). This marks the best uncorrected results, and in the calculations we set γ min /J = 0.02, which is equal to γ min /ξ ≈ 0.16. Being only a one-sixth of the interaction strength, it certainly makes the measurement less feasible, and we choose it as an exaggerated case aiming to see if the quantum inverse iteration can deal with high-noise operation mode. Taking the measurements at effectively elevated noise rates with γ > γ min (bold circles), the overlap is extrapolated to γ/J = 0 value. This is done using linear (a = 1) and cubic (a = 3) spline fitting, which corresponds to thin dashed and thick solid lines in Fig. S3B,C. The gray shaded area corresponds to the inaccessible range of noise rates.
We observe that for smaller propagation phase δφ/2π = 0.68 (B) extrapolation of both direct and indirect measurement provides a decent estimates for the overlap O  (Fig. S3C). The procedure is complicated by nonlinear overlap dependence. Again, direct estimation yields better estimate, although we note that this is not a universal trend. Having collected the information, we also compose the weighted averaged overlap value as (0)] 2 . The reasoning behind the heuristic weighting procedure is as follows. First, the noise dependence for direct and indirect measurement of the overlap has opposite trends, with each deviating into higher or lower overlap values. This calls for a technique which can use the average of the two. Second, typically the cubic extrapolation with smaller curvature is more trustful; A Histogram of weights p k,q for different wavefunction overlaps generated by the evolution with phase differences δφq. Approximation grid is fixed for all panels, and is set by My = Mz = 5, ∆yJ = ∆z = 0.5, with φmax/2π ≈ 1. Blue, magenta, and red lines correspond to k = 1, 2, 3 iteration steps, respectively. Each case is normalized by the total weight at given k. B, C Real parts of the wavefunction overlaps Re{ ψq|ψ0 } = Re{Oq} as a function of the effective dephasing rate γ. The Monte-Carlo procedure was performed using 5000 trajectories for each point, and the phase difference was fixed to δφq/2π = 0.68 (B) and δφq/2π = 1.4 (C). Two types of overlap estimation are used, given by direct and indirect inference from Eqs. (M4) and (M6) of Methods sec. B, respectively. The noiseless case of γ = 0 is shown by the horizontal blue line. The error mitigation is performed by using results of extrapolation of linear (a = 1) and cubic (a = 3) order to γ = 0 value. The result for the weighted procedure [see Eq. (S10)] is shown by the red circle and highlighted as extrapolated value. Grey shaded region contains the results for small dephasing rate, and marks the inaccessible region excluded from the mitigation. D Combined results of quantum inverse iteration as a function of iteration step, shown for different data analysis strategies. Top black and grey curves correspond to noisy data with direct and indirect inference, respectively, taken at γmin/J = 0.02 (γmin/ξ ≈ 0.16). The red curve shows the error mitigated result. The lowest blue curve is for the noiseless operation for φmax/2π ≈ 1. In each case, the difference between true ground state and its estimate version is plotted (log scale), and blue shaded region denotes the chemical precision. but when it largely deviates from linear approximation we enter the regime of high noise, and its likely that actual solution lies in between direct and indirect extrapolation result. The final result of the weighted mitigation is depicted in Figs. S3B,C by circles labelled as "extrapolated value".
Collecting the overlap data and coefficients, in Fig. S3D we plot the results of ground state estimation procedure for different post-processing techniques. The noiseless quantum inverse iteration is shown by the lowest (blue) curve, and serves as a reference. For noisy operation, the two upper curves correspond to overlap estimates without error mitigation by taking their values at lowest achievable noise O (ind) (γ min ) and O (dir) (γ min ) for an indirect and a direct inference only. While the former stays outside of chemical precision for all iteration steps due to large overlap deviations (black curve), the direct inference technique allows to reach chemical precision even at high noise (γ min /ξ = 0.16). Applying the error mitigation technique as stated in Eq. (S10), the results can be improved for initial iteration steps, but approach unmitigated values for k > 10. This can be explained by the change of the weight distribution for large k, where overlaps at large times are important (Fig. S3A), and weighted extrapolation does not provide good estimate in this case.