Toward Reliability in the NISQ Era: Robust Interval Guarantee for Quantum Measurements on Approximate States

Near-term quantum computation holds potential across multiple application domains. However, imperfect preparation and evolution of states due to algorithmic and experimental shortcomings, characteristic in the near-term implementation, would typically result in measurement outcomes deviating from the ideal setting. It is thus crucial for any near-term application to quantify and bound these output errors. We address this need by deriving robustness intervals which are guaranteed to contain the output in the ideal setting. The first type of interval is based on formulating robustness bounds as semi-definite programs, and uses only the first moment and the fidelity to the ideal state. Furthermore, we consider higher statistical moments of the observable and generalize bounds for pure states based on the non-negativity of Gram matrices to mixed states, thus enabling their applicability in the NISQ era where noisy scenarios are prevalent. Finally, we demonstrate our results in the context of the variational quantum eigensolver (VQE) on noisy and noiseless simulations.


I. INTRODUCTION
Today's quantum computers are characterized by a low count of noisy qubits performing imperfect operations in a limited coherence time.In this era of quantum computing, the noisy intermediate-scale quantum (NISQ) era [1], researchers and practitioners alike strive to heuristically harness limited quantum resources in order to solve classically difficult problems and also to benchmark and potentially develop new quantum subroutines.A typical pattern of these NISQ algorithms [2], exemplified by the seminal variational quantum eigensolver (VQE) [3] and quantum approximate optimization algorithm (QAOA) [4], consists of the preparation of ansatz states with a parameterized unitary circuit followed by useful classical output being extracted by means of quantum measurements, more generally as expectation values of quantum observables through repeated measurements.
The promising potential of these NISQ algorithms spans across a wide spectrum of applications, ranging from quantum chemistry, many-body physics, and machine learning to optimization and finance [2].However, as a consequence of their heuristic nature and the prevalent imperfections in near-term implementation, NISQ algorithms in practice typically produce outputs deviating from the exact and ideal setting.This unfortunate hindrance practically arises from various sources such as circuit noise and decoherence [1], limited expressibility of ansatze [5,6], barren plateaus during optimization in variational hybrid quantum-classical algorithms [7][8][9], measurement noise and other experimental imperfections [10,11].To determine the usefulness of a given NISQ application, it is thus crucial to quantify the error on the final output in the presence of a multitude of the aforementioned sources of imperfection.
In this work, we endeavour to systematically certify the reliability of quantum algorithms by deriving robustness bounds for expectation values of observable on approximations of a target state.To that end, based on analytical solutions to a semidefinite program (SDP), we present lower and upper bounds to expectation values of quantum observables which depend only on the fidelity with the target state and post-processing of previously obtained measurement results.Furthermore, we take into account higher statistical moments of the observable by generalizing the Gramian method for pure states [12] to generic density operators, thus extending its application to bounding output errors resulting from noisy circuits.Although the focus of our investigation is on errors arising from circuit imperfection, the underlying techniques are also valid for other sources of errors such as algorithmic shortcomings.Finally, we apply these bounds to numerically obtain robustness intervals on simulated noisy and noiseless VQE for ground state energy estimation of electronic structure Hamiltonians of several molecules.The robustness certification protocol resulting from this work is integrated with the open source Tequila [13] library.
The remainder of the paper is organised as follows.In Sec.II, we present our main results, namely, the bounds based on SDP and the Gramian method.In Sec.III, we present our numerical simulations and explain the applicability of our bounds in the context of VQE.Sec.III B highlights the implementation in Tequila and concluding remarks are given in Sec.IV.

II. ROBUSTNESS INTERVALS
The goal of this work is to provide techniques to compute intervals which are guaranteed to contain the expectation value of an observable A under an ideal, but unavailable, target state σ.Any such interval is referred to as a robustness interval.More formally, instead of having access to the state σ, we assume access to the approximate state ρ and which is further assumed to have at least fidelity 1 − with the target state σ.Given these assumptions, we define a robustness interval to be an interval I = [ χ, χ] ⊆ R for which it is guaranteed that and which is a function of the infidelity , the observable A, and the state ρ.
The Hilbert space corresponding to the quantum system of interest is denoted by H ≡ C d with dimension d = 2 n .We use the Dirac notation for quantum states, i.e. elements of H are written as kets ψ⟩ ∈ H with the dual written as a bra ⟨ψ .The space of linear operators acting on elements of H is denoted by L(H) and elements thereof are written in capital letters A ∈ L(H).The set of density operators on H is written as S(H) ⊂ L(H) and lower case greek letters are used to denote its elements σ ∈ S(H) which are positive semidefinite and have trace equal to 1.For an element A ∈ L(H) we write A ≥ 0 if it is positive semidefinite, A T stands for its transpose, and A † is the adjoint.We also use the Loewner partial order on the space of Hermitian operators, i.e. for two Hermitian operators A, B ∈ L(H), we write A ≥ B if and only if A − B ≥ 0. Expectation values of observables, i.e.Hermitian operators A ∈ L(H), are written as ⟨A⟩ σ = Tr[Aσ] for some σ ∈ S(H).The variance of an observable is given by where the maximum is taken over all purifications of ρ and σ.For pure states, the fidelity reduces to the squared overlap F( ψ⟩⟨ψ , φ⟩⟨φ ) = ⟨ψ φ⟩ 2 .Finally, the real part of a complex number z ∈ C is written as R(z) and the imaginary part as I(z).

B. Summary of technical results
We employ two different techniques to bound the true expectation value ⟨A⟩ σ , each with its advantages and disadvantages in terms of efficiency and accuracy.The first FIG. 1.Bond dissociation curves and robustness interval (RI) for Lithium Hydride in a basis-set-free approach [14,15].The exact, theoretical energies are shown in black, the energy estimates provided by a noisy VQE with an UpCCGSD Ansatz [16] is shown in blue.The robustness interval is guaranteed to contain the true ground state energy and is based on the Gramian eigenvalue bounds for mixed states (Theorem 3).
technique is based on the formulation of lower and upper bounds as SDPs and makes use of a closed form solution of optimal type-II error probabilities from quantum hypothesis testing [17].The second technique is based on the non-negativity of the determinant of Gramian matrices for a suitable collection of vectors.This second technique was initially proposed by Weinhold [12] in the context of pure states.Using Uhlmann's Theorem [18], which relates the fidelity between two mixed states to the trace norm, we extend these results to mixed states.This ultimately justifies their applicability in the current NISQ era, where the assumption of a closed quantum system is violated and one needs to make use of the density operator formalism to accurately model these states and their evolutions.
In Table I, we summarize all the results, together with the conditions under which they apply and the quantities that are covered.Fig. 1 shows the ground state energies of molecular Lithium Hydride in the basis-set-free approach of Refs.[14,15], with energy estimates provided by VQE with an UpCCGSD ansatz.The lower and upper bounds on the true energy are obtained via the Gramian method from Theorem 3.

C. Bounds via Semidefinite Programming
Here we derive a robustness interval which is based on expressing lower and upper bounds as a semidefinite program which we connect to optimal type-II error probabilities for binary quantum hypothesis testing (QHT).
a.Quantum Hypothesis Testing.Binary QHT can be formulated in terms of state discrimination where two states have to be discriminated through a measurement.On a high level, the goal is to decide whether a quantum system is either in the state ρ, referred to as the null hypothesis, or in the state σ, the alternative hypothesis.Any such hypothesis test is represented by an operator 0 ≤ Λ ≤ 1 which corresponds to rejecting the null ρ in favor of the alternative σ.The central quantities of binary QHT are the two different probabilities of making an error, namely the type-I and type-II error probabilities, defined as (type-II error) and which quantify the probability of rejecting the null hypothesis when it is true, and accepting the null when the alternative is true, respectively.One seeks a test Λ which minimizes the probability of making a type-II error under the constraint that the type-I error is below some predefined threshold α 0 .Formally, we have the SDP The following Lemma establishes a closed form solution for this SDP for pure states, and a lower bound for the general case of mixed states.
In the following, we use this result to get closed form solutions for robustness bounds formulated as SDPs.

Robustness Interval
Consider a bounded observable −1 ≤ A ≤ 1 and let ρ be the approximate state, corresponding to the alternative hypothesis, and let σ be the target state, corresponding to the null hypothesis.We can express lower and upper bounds to ⟨A⟩ σ as semidefinite programs which take into account measurements of ρ.Namely, we have the upper bound and the lower bound It is straight forward to see that these optimization problems are indeed valid lower and upper bounds to ⟨A⟩ σ by noting that the operator A is feasible.In addition, as shown in Appendix A 3, the tightness of the bounds is an immediate consequence of the formulation of the robustness interval as an SDP.We can rewrite these SDPs and express them in terms of optimal type-II error probabilities, so that the upper bound reads and, similarly, for the lower bound This establishes a close connection between state discrimination via hypothesis testing, and the robustness of expectation values under perturbations to states.Indeed, these robustness bounds formalize the intuition that states which are hard to discriminate, i.e., which admit higher error probabilities, will have expectation values which are closer together.Furthermore, this connection also has the interesting interpretation that, if the approximate expectation ⟨A⟩ ρ is close to the extreme −1, a statistical hypothesis test is restricted to have type-I error probability close to 0. This makes it harder for the corresponding optimal type-II error probability β * to be low and hence ⟨A⟩ σ will generally be closer to ⟨A⟩ ρ .Finally, Lemma 1 provides a closed form solution to the SDP β * , which only depends on the fidelity between ρ and σ and hence establishes a robustness interval of the form Eq. ( 1).This result is summarized in the following Theorem: Theorem 1.Let σ, ρ ∈ S(H d ) be density operators with F(ρ, σ) ≥ 1 − for some ≥ 0. Let A be an observable with −1 ≤ A ≤ 1 and with expectation value ⟨A⟩ ρ under ρ.For ≤ 1 2 (1 + ⟨A⟩ ρ ), the lower bound of ⟨A⟩ σ can be expressed as Similarly, for We remark that, although used in a different context, this technique is conceptually similar to the result presented in Theorem 1 of Ref. [17] for adversarial robustness of quantum classification.Stemming from the formulation as an SDP, these bounds are tight for pure states in the sense that, for each bound, there exists an observable A with expectation ⟨A⟩ ρ under ρ and whose expectation under σ saturates the bound.In Appendix A 3, we give a formal, constructive proof for this statement.Furthermore, in practice, it is typically not feasible to measure the exact value of ⟨A⟩ ρ due to finite sampling errors, measurement noise, and other experimental imperfections.For this reason, one needs to rely on confidence intervals which contain the exact value of ⟨A⟩ ρ with high probability.This can be accounted for in the bounds from Theorem 1 by noting that they are monotonic in ⟨A⟩ ρ , what allows us to replace the exact value by bounds which hold with high probability.Finally, it is worth noting that, if one has access to an estimate of the fidelity, i.e. some > 0 with F(ρ, σ) ≥ 1 − , this interval can be calculated by merely postprocessing previous measurement results, and hence does not cause any computational overhead.

D. Bounds via non-negativity of the Gramian
Here we employ a different technique to derive robustness bounds, taking into account the variance of the observable as an additional piece of information.The method is based on the non-negativity of the Gramian and was pioneered by Weinhold [12].We first give a brief overview of the Gramian method and then present the extension to mixed states.This extension is important as the restriction to pure states hinders the applicability of this method in practice and, in particular, in the current NISQ era, where one often has to deal with noisy states that are expressed as probabilistic ensembles of pure states in the density operator formalism.

Review of the Gramian method
Consider a Hermitian operator A ∈ L(H), a target state ψ⟩ and an approximation of this state φ⟩.The Gram matrix for the vectors ψ⟩, φ⟩, A φ⟩ is given by where, without loss of generality, it is assumed that the overlap ⟨ψ φ⟩ is real (otherwise multiply each state by a global phase).Since Gram matrices are positive semidefinite (e.g., Theorem 7.2.10 in [19]), their determinants are non-negative, and the non-negativity of G thus limits the permissible values of R(⟨ψ A φ⟩) to be within these boundaries.Since I(⟨ψ A φ⟩) 2 ≥ 0 we have the bounds Starting from these inequalities, bounds for expectation values of quantum observables have been derived for pure states [12,[20][21][22][23][24][25][26] and in the context of classical methods for quantum chemistry.While certainly useful, this leaves a gap for these bounds to be applied in the NISQ era where noise is prevalent and quantum states and their evolutions are described by density operators.In the following section, we fill this gap and extend the technique to mixed states.

Gramian method for mixed states
Here, we build on the Gramian method and derive bounds which are valid for mixed states and which, in contrast to Theorem 1, take into account the second moment of the observable of interest.In principle, as more information is included, one can expect that this results in a tighter bound at the cost of having to measure additionally the expectation value of A 2 .The following theorem provides a lower bound to expectation values of non-negative observables: Theorem 2 (Expectation values).Let σ, ρ ∈ S(H d ) be density operators with F(ρ, σ) ≥ 1 − for some ≥ 0. Let A ≥ 0 be an arbitrary observable with expectation value ⟨A⟩ ρ under ρ.For with (1 − ) ≥ ∆A ρ ⟨A⟩ ρ , the lower bound of ⟨A⟩ σ can be expressed as In the case where the target state σ is an eigenstate of an observable A, the Gramian method allows to derive a further bound.While the assumptions here are stronger, this bound is particularly useful in applications such as the variational quantum eigensolver and when the observable of interest commutes with a Hamiltonian H for which the target state is an eigenstate.Formally, we have the following result: Theorem 3 (Eigenvalues).Let ρ ∈ S(H d ) be a density operator and let A be an arbitrary observable with eigenstate ψ⟩ and eigenvalue λ, A ψ⟩ = λ ψ⟩.Suppose that ≥ 0 is such that F(ρ, ψ⟩) = ⟨ψ ρ ψ⟩ ≥ 1 − .Then, lower and upper bounds for λ can be expressed as The proofs for both Theorems presented in this section start from Eq. ( 13) which is first extended to mixed state via purifications and an application of Uhlmann's Theorem.The second step of the proof is a rewriting of the inequality in the case where the target state is an eigenstate of A (Theorem 3), and, in the case of Theorem 2, an application of the Cauchy-Schwarz inequality.The proofs are given in Appendix B and Appendix C.

E. Comparison of the bounds
We have seen three different methods to derive robustness intervals.Namely, the interval based on SDP given in Theorem 1, the expectation value lower bound based on the Gramian method from Theorem 2, and the robustness interval for eigenvalues from Theorem 3, which is also based on the Gramian method.As a first observation, we notice that the SDP bounds are dependent only on the first moment of the observable, while the bounds derived from the Gramian method take into account the second moment via the variance.In principle, this hints at a trade-off between accuracy and efficiency.That is, by taking into account higher moments, which comes at a higher computational cost, one can hope for an improvement in accuracy as more information is included.On the other hand, the SDP bounds can be calculated as a postprocessing step and thus do not require to measure additional statistics.However, as less information is included, this typically comes at the cost of lower accuracy.
On the practical side, one needs to consider that for the SDP bounds to be applicable, it is required that the observable lies between −1 and 1.In practice, however, this is not always the case and the observable needs to be appropriately rescaled, e.g. by using its eigenvalues.As the exact eigenvalues might not be available, one needs to use lower and upper bounds on these, which results in a loss in tightness.This is because the set of feasible points in the SDP from Eq. ( 5) and Eq. ( 4) becomes larger and hence loosens the bounds.A similar issue emerges in the lower bound for expectation values based on the Gramian method where the observable needs to be positive semidefinite.If this assumption is violated, one again needs to apply an appropriate transformation of the observable, leading to a potentially looser bound.Instead of scaling, one can decompose the observable into individual terms, each satisfying the constraints, and then bound each term separately and aggregate the bounds over the decomposition.In Sec.III we consider such a decomposition in the context of VQE.Specifically, we decompose the underlying Hamiltonian into groups of mutually commuting Pauli terms and bound the expectation of each group separately.In contrast, the eigenvalue bound based on the Gramian method does not suffer from these issues and it is applicable for general observables.It is worth remarking that this comes at the cost of less generality in the sense that the bound only applies to eigenvalues rather than general expectation values.
Assuming that the observable A = P is a projection, satisfying P 2 = P , allows for a direct comparison between the bounds.Note that, in this case, the variance is fully determined by the first moment via (∆P ρ ) 2 = ⟨P ⟩ ρ − ⟨P ⟩ 2 ρ and we expect that the Gramian Expectation bound should not be tighter than the SDP bound.First, we incorporate the knowledge that P is a projection in the SDP lower bound by applying it to the observable 2P − 1 so that we have the bound which is exactly the same as the lower bound derived via the Gramian method in Theorem 2 when applied to the projection P .We can also compare this bound to the Gramian eigenvalue bound from Theorem 3. Since the latter is less general, in the sense that it only holds for target states which are eigenstates, we expect this to be tighter than the expectation value bound.As can be seen from Fig. 2, this is indeed the case.Finally, we notice that all of the above bounds are faithful in the sense that, as the approximation error vanishes → 0, the bounds converge to the true expectation value ⟨A⟩ σ .To compare the rate of convergence, consider the case of pure states with the target state given by σ = ψ⟩⟨ψ and the approximation state ρ = φ⟩⟨φ with φ⟩ ±1 for the SDP bounds.This ultimately stems from the underlying assumptions required for the bounds to hold.In contrast, the Gramian eigenvalue bound has no assumptions on the observable A and the bounds diverge as approaches 1.

F. Fidelity estimation
All bounds presented so far have in common that they depend on the fidelity with the target state σ.However, in many practical settings, it is not possible to access the target state and thus difficult to obtain even a lower bound to the true fidelity.Here we seek to address this topic and present lower bounds on the true fidelity for the case where the target state is the ground state of a Hamiltonian H.
Let H be a Hamiltonian with eigenvalues λ 0 ≤ λ 1 ≤ . . .≤ λ d and assume that λ 0 has geometric multiplicity 1 so that the corresponding ground state ψ 0 ⟩ is unique.Let ρ be a possibly mixed state approximation of ψ 0 ⟩.If both λ 0 and λ 1 are known, one can make use of Eckart's criterion [27] to bound the fidelity via In scenarios where knowledge of the lowest lying eigenvalues λ 0 and λ 1 is available, one can thus directly lowerbound the fidelity and use Eq. ( 17) in the computation of the robustness intervals.In scenarios where one does not have full knowledge of these eigenvalues (or, in the least, corresponding bounds), Eckart's criterion cannot be directly applied.However, we can still use the inequality if less knowledge about the spectrum of H is available.If it is known that the energy estimate ⟨H⟩ ρ is closer to λ 0 than to λ 1 then, as an immediate consequence of Eckart's criterion, one finds that We remark that substituting Eq. ( 18) into the Gramian eigenvalue bounds from Theorem 3 yields the mixed state extension of the Weinstein bounds [28,29] in the nondegenerate case.If, in addition, a lower bound δ on the spectral gap is known such that λ 1 − λ 0 ≥ δ, then we have the bound derived in Ref. [30], which is a nontrivial lower bound whenever the variance is small enough such that ∆H ρ ≤ δ.With a similar technique, one obtains a further tightening of the bound: for variances with ∆H ρ ≤ δ 2. We note that this bound has also been reported in Refs.[31,32] and we provide an alternative proof in the Appendix.In practice, the bound δ on the spectral gap can also be estimated via classical methods, as for example truncated classical configuration interaction or density-matrix renormalization group techniques.In principle, also non-variational methods like truncated coupled-cluster (and the associated equationof-motion or linear-response variants for the excited state energies) could be applied.In either case, the idea is to use these classical methods to compute the ground state and the first excited state energies to get an estimate of the spectral gap which can then be used for the fidelity estimation.The classical method which is the best to choose will generally depend on the system of interest and the available computational time.We refer the reader to [33] for a detailed treatment over some of those methods.
The above bounds hold for Hamiltonians H whose lowest eigenvalue is non-degenerate.In Appendix D we consider the degenerate case and show that when the approximate state ρ is pure, then there always exists a state ψ⟩ which is an element of the eigenspace associated with the lowest (possibly degenerate) eigenvalue, and for which the above fidelity lower bounds hold.However, if ρ is a mixed state, this cannot in general be said, as is shown in the appendix with a counterexample.In summary, when the approximate state ρ is allowed to be mixed, then the fidelity bounds are applicable only when the underlying Hamiltonian has a non-degenerate ground state.If, on the other hand, ρ is pure, then the bounds also hold in the degenerate case.Finally, we remark that these fidelity bounds all require varying degrees of knowledge about the ground state and Hamiltonian in question.They thus can only partially address the topic of fidelity estimation in scenarios where such knowledge is not available.
At this point we would like to point out an interesting connection to variational quantum time evolution (Var-QTE).In general, VarQTE is a technique to find the ground state of a Hamiltonian H [34][35][36] by projecting the time evolution of the initial state to the evolution of the ansatz parameters.VarQTE typically comes with an approximation error, stemming from a limited expressibility of the Ansatz state or from noise.In Ref. [37], this approximation error is quantified in terms of an upper bound on the Bures distance between the evolved state and the true ground state.Since there is a one-toone correspondence between Bures distance and fidelity, these error bounds can be converted to a lower bound on the latter.This in turn can then be used to calculate the eigenvalue and expectation bounds presented in this work.

III. APPLICATIONS
In this section, we put into practice the theoretical results presented in the previous sections and calculate the robustness intervals for ground state energies of electronic structure Hamiltonians when the approximation of the ground state is provided by VQE.We remark that, while VQE serves as an example application, our results are not limited to ground state energies but can be used in a more general context where the goal is to calculate error bounds for expectation values.Consider a Hamiltonian H with Pauli decomposition and let ρ be an approximation to the true ground state ψ 0 ⟩.Given ≥ 0 such that ⟨ψ 0 ρ ψ 0 ⟩ ≥ 1 − , and an estimate of the variance ⟨H 2 ⟩ ρ − ⟨H⟩ 2 ρ , it is straightforward to evaluate the Gramian eigenvalue bounds from Theorem 3. In contrast, for the expectation value bounds derived via SDP and the Gramian method from Theorems 1 and 2), one needs to be more careful since the Hamiltonian H might violate the underlying assumptions.To evaluate the latter, we can account for this by adding a sufficiently large constant c such that H ∶= H +c1 ≥ 0 and calculate the bound for H, before reversing the translation in order to get the desired bound for H. Clearly, a valid choice for c is given by −λ 0 where λ 0 is the lowest eigenvalue of H.However, it is not always clear which constant c leads to the tightest lower bound.Similarly, to evaluate the SDP bounds, we need to apply Theorem 1 to operators which are bounded between ±1.If the full spectrum of H was known, one could normalize H using these eigenvalues.However, in the context of VQE, the spectrum is not a priori known as this is precisely the task that VQE is designed to solve, and we need a different approach for the expectation value bounds.The idea is to partition the terms in the Pauli decomposition from Eq. ( 21) into groups so that each term corresponding to a group can be normalized.To this end, we first partition H into groups of mutually qubit-wise commuting terms Given such a representation, the spectrum of each of the H k can be calculated classically in order to scale H k → Hk appropriately such that the assumptions for the bounds are satisfied.One can then compute the bounds for each of the terms in the summation and get the final bounds by aggregating the individual bounds.We further make use of the approach presented in Refs.[38,39] where one applies a unitary transform U k to each of the H k terms so that single-qubit measurement protocols can be used.Specifically, instead of measuring H k under the state ρ, one measures k under the unitarily transformed U k ρU † k .One can then scale each A k appropriately by classically computing its eigenvalues and apply the expectation value bounds (Theorems 1 and 2) to each term separately before aggregating.It is also worth noting the generality of Theorem 1: although in the preceding demonstration, the matrix A is generally taken to be a Pauli observable for measuring the output of a quantum circuit, the condition −1 ≤ A ≤ 1 is satisfied much more generally (e.g. by Fermionic operators).The application of this theorem in settings without explicit Pauli decomposition would be a fruitful ground for future research.

A. Numerical simulations
Here, we numerically investigate the different robustness bounds for the ground state energies for a set of electronic structure Hamiltonians, namely H 2 , LiH and BeH 2 molecules where the qubit Hamiltonians are obtained within the basis-set-free approach of Ref. [14] using directly determined pair-natural orbitals on MP2 level [15].All our experiments have been implemented in the Tequila [13] library using the qubit encodings from openfermion [40], optimizers from scipy [41], madness [42] as chemistry backend, qulacs [43] as the simulation backend for noiseless simulations and qiskit [44] for simulations which include noise.We model noise as a combination of bitflip channels acting on single qubit gates with 1% error probability, and depolarizing noise acting on two qubit gates, also with an error probability of 1%.
For a given Hamiltonian H, we first approximate its ground state ψ 0 ⟩ via VQE.That is, for an Ansatz state ρ θ with parameters θ one minimizes the objective ⟨H⟩ ρ θ and obtains optimal parameters θ * = arg min θ ⟨H⟩ ρ θ .It follows from the Rayleigh-Ritz variational principle [45,46] that the expectation ⟨H⟩ ρ * θ is an upper bound to the true ground state energy λ 0 .The such obtained ground state approximation ρ θ * is then used to evaluate the bounds by computing the relevant statistics, i.e. expectation values and variances of observables under this state.We notice that the quality of this state in terms of a distance to the true ground state is not easily obtainable without having some prior knowledge over the system of interest (see also Sec.II F in this regard).For this reason and in order to investigate and compare the bounds, here we assume knowledge of the true fidelity with the ground state.In practice, this is of course not realistic and, as discussed previously, one needs to approximate the true fidelity.Given the ground state approximation ρ θ * and the fidelity F(ρ θ * , ψ 0 ⟩⟨ψ 0 ), we then estimate the expectation values and variances under ρ θ * in order to evaluate the bounds.In the noiseless sce-nario, these statistics can be calculated exactly, whereas in the noisy scenario they need to be estimated due to finite sampling errors.Thus, in the noisy case, we repeat the calculation of the bounds 20 times and report one-sided 99%-confidence intervals.In Fig. 3, we consider the noisy scenario and compare the different types of bounds for H 2 (2, 4) and LiH(2, 4) with approximation states provided by VQE with an UpCCGSD Ansatz [16] and optimized fermionic gradients [47].For both molecules, we notice that the Gramian eigenvalue bound is the tightest, while the expectation value bounds are less tight.However, this is not surprising, as the eigenvalue bound is more suited for this task, compared to the other bounds which hold more generally for expectation values.In Fig. 4, we again consider the noisy scenario and compare the Gramian eigenvalue bounds for approximation states obtained via the SPA Ansatz [48] and via the UpCCGSD Ansatz for H 2 (2, 4), LiH(2, 4) and BeH 2 (4,8).We first notice that the SPA Ansatz is generally less vulnerable to noise, which stems from the associated shallow circuits, compared to the UpCCGSD Ansatz.In particular, SPA and UpC-CGSD have the same expressibility for H 2 but, since SPA uses more efficient compiling, its energy estimates and lower bounds are more accurate compared to Up-CCGSD.In Appendix E we show robustness intervals for the LiH(2, 4) molecule with the error rate increased to 10%.For the SPA ansatz, even with this error rate, the ground state fidelities vary between 0.51 and 0.65 while the UpCCGSD states have low ground state fidelities in the range 0.1.In other words, UpCCGSD fails to converge to states which are close to the true ground state.It is interesting to note that for lowest ground state fidelities, the expectation value bounds reduce to trivial bounds and the eigenvalue bound starts to diverge.In Fig. 5 we consider the noiseless scenario with an SPA ansatz for LiH and BeH 2 .In contrast to the noisy scenario, here the bounds based on the UpCCGSD Ansatz are tigther, compared to the ones based on the SPA Ansatz for large bond distances.This is due to the fact that SPA generally has more difficulties in approximating ground states for far stretched bond distances and hence results in lower ground state fidelities.Finally, it is worth remarking that these bounds are obtained under the assumption of having complete knowledge of the true ground state fidelity, an assumption which is idealistic and typically violated in practice.

B. Implementation
All our robustness intervals are implemented in the open source Tequila [13] library.In the following example, we run VQE for the H 2 Hamiltonian in a minimal representation (4 qubits), before computing the lower and upper bounds based on the optimized circuit, using the function robustness_interval:  for eigenvalues based on the Gramian method (Theorem 3) for LiH (2,4)and BeH2 (4,8).Here, an ideal scenario without noise is simulated and the approximation errors stem from the limited expressibility of the Ansatz state.
1 import tequila as tq 2 from tequila .apps .robustness import r o b u s t n e s s _ i n t e r v a l Used in this way, the function calculates the robustness interval for all three methods and returns the tightest bounds.Alternatively, one can specify the type of bound via the keywords kind and method where the former stands for which kind of interval is desired, that is expectation or eigenvalue, and the latter stands for the method used to obtain the bound (Gramian or SDP).For example, calculating a robustness interval for eigenvalues using the Gramian method, can be implemented as robustness_interval(..., kind="eigenvalue", method="gramian").In general, any type of expectation value can be used.Note that our implementation is agnostic with respect to the molecular representation, so that replacing n_pno=1 with basis_set="sto-3g" will lead to a 4 qubit Hamiltonian in a traditional basis set.
IV. DISCUSSION AND CONLUSIONS The current experimental stage of quantum computation offers the possibility to explore the physical and chemical properties of small systems and novel quantum algorithms are being developed to extract the most from this first generation of quantum devices.However, this potential for computational advantage, compared to classical methods, comes at a price of noisy and imperfect simulations stemming from low qubit counts and thus the lack of quantum error correcting qubits.The Variational Quantum Eigensolver (VQE) is the canonical example of these NISQ algorithms that allow us to obtain an approximation of Hamiltonian eigenstates by exploiting the variational principle of quantum mechanics.Besides the broad applications and promising results of this approach [3,30,49], its performance guarantees should be studied and understood.
In this work we have made first progress in this direction and have derived robustness intervals for quantum measurements of expectation values.For a target state σ, these intervals are guaranteed to contain the true expectation value ⟨A⟩ σ of an observable A when we only have access to an approximation ρ.Based on resource constraints, accuracy requirements, and depending on the task, we have seen three different types of robustness intervals.Firstly, based on the formulation of robustness bounds as SDPs, we have derived upper and lower bounds to ⟨A⟩ σ which take into account only the first moment of the observable A and can thus be obtained by post-processing measurements of ⟨A⟩ ρ together with the fidelity with the target state F(ρ, σ).Secondly, we have revisited the Gramian method [12] to take into account higher statistical moments of A and extended this technique to mixed states, thereby enabling their applicability in noisy scenarios which are prevalent in the NISQ era.This has led to a further lower bound to expectation values and, additionally, to lower and upper bounds on eigenvalues of observables.We have also implemented our bounds in the open source Tequila [13] library.Furthermore, we have validated our results with numerical simulations of noisy and noiseless scenarios with VQE as an example application to calculate robustness intervals for ground state energies of electronic structure Hamiltonians of H 2 , LiH and BeH 2 .For the molecules considered in these experiments, we have observed that the robustness intervals provide accurate estimates of the errors incurred by noise, in particular when the ground state approximation is close enough to the true ground state in terms of fidelity.
The main requirement of the bounds obtained is the knowledge of the fidelity between the target state and its approximation.Although such a quantity is not always experimentally accessible and hence poses a challenge in the practical applicability of these bounds, there exist algorithms, such as within the variational quantum imaginary time evolution [37] framework, which allow for a quantification of the required approximation error in terms of distances between the target and the approximate state.Nonetheless, our aim is to provide a formal framework to study the robustness of broadly used approaches as are the Variational Quantum Algorithms.There are still many questions around the applicability of these quantum algorithms and its robustness against noise.Within this work, we seek to unravel the uncertainties around these state-of-the-art quantum algorithms with the goal of improving its performance and applicability.

ACKNOWLEDGMENTS
The authors are grateful to Joseph Fitzsimons (Horizon Quantum Computing) and Nana Liu (Shanghai Jiao Tong University) for inspiring discussions on the topic of robustness of NISQ algorithms.for which we find eigenvalues The corresponding eigenvectors η 0 ⟩, η 1 ⟩ are determined by their coefficients z k,ψ , z k,φ for k = 0, 1, for which we find where the coefficient A k arises from requiring the eigenvectors η k ⟩ to be normalized (see Ref. [51], Section 8).Note that ∀t ≥ 0 we have η 0 ≤ 0 and η 1 > 0 so that P t,+ = η 1 ⟩⟨η 1 .Recall that we defined t 0 to be the positive number t 0 ∶= inf{t ≥ 0∶ Tr[P t,+ ρ] ≤ α 0 }.It follows that and notice that the right hand side is continuous in t over [0, ∞) whenever γ < 1.Since t ↦ Tr[P t,+ ρ] is nonincreasing, it's maximum is attained at t = 0 so that Tr[P t,+ ρ] ≤ γ 2 and hence t 0 = 0 if α 0 > γ 2 .In this case, we obtain β * (α 0 ; ρ, σ) = 0. If, on the other hand, α 0 ≤ γ 2 , then we solve the equation Tr[P t,+ ρ] = α 0 and obtain the expression for t 0 For t = t 0 we have η 0 < 0 and η 1 > 0 so that Λ t0 = η 1 ⟩⟨η 1 and η 1 ⟩ = −γA 1 ψ⟩ + (1 − η 1 )A 1 φ⟩.Hence Plugging t 0 into the expressions above yields Since the right hand side of Eq. (A22) is monotonically decreasing in γ 2 and γ 2 ≥ 1 − , the claim follows for pure states.To see that the above expression also constitutes a valid lower bound for mixed states, let Ψ and Φ be arbitrary purifications of σ and ρ respectively, both with purifying system H E .It is well known that β * is monotonically increasing under the action of any completely positive and trace preserving map E, i. where Tr E [⋅] denotes the partial trace over the purifying system.It follows from Uhlmann's Theorem that we can choose Ψ, Φ such that ⟨Ψ Φ⟩ 2 = F(ρ, σ).The claim now follows from the observation that the RHS of Eq. (A23) is monotonically decreasing in ⟨Ψ Φ⟩ 2 .This completes the proof.

Main Proof
Given Lemma 1, here we provide the full proof for Theorem 1.We first restate the result.

(B14)
Notice that the RHS is monotonically decreasing as F(ρ, σ) decreases.Hence, we can replace the true fidelity by a lower bound to it.In particular, for ≥ 0 with F(ρ, σ) ≥ 1 − and (1 − ) ≥ ∆A ρ ⟨A⟩ ρ we get which is the desired result.

0 FIG. 2 .
FIG.2.Difference between the SDP bounds from Theorem 1 and the Gramian eigenvalue bound from Theorem 3 as a function of the infidelity and the expectation value ⟨P ⟩ρ.The observable is assumed to be a projection P and the target state is an eigenstate of the observable.The difference is calculated by subtracting the SDP bound from the Gramian bound.a) shows the difference between lower bounds, b) shows the difference between the upper bounds.As can be seen from the figures, the Gramian eigenvalue bound is always more accurate than the expectation bound.Note that the Gramian expectation value lower bound (Theorem 2) equals the SDP lower bound under these assumptions.

FIG. 3 .
FIG.3.Comparison of the different lower bounds (LB) and robustness intervals (RI) presented in Sec.II for bond dissocation curves of H2(2, 4) and LiH(2,4).The approximation states are provided by VQE with an UpCCGSD ansatz.Both the VQE optimization and the evaluation of the bounds were simulated with bit flip and depolarization noise with 1% error probability.

FIG. 5 .
FIG. 5. Bond dissociation curves, robustness interval (RI)for eigenvalues based on the Gramian method (Theorem 3) for LiH(2,4)and BeH2(4, 8).Here, an ideal scenario without noise is simulated and the approximation errors stem from the limited expressibility of the Ansatz state.

TABLE I .
Overview of bounds for the true expectation values and eigenvalues of a Hermitian operator A, with σ the target state and ρ the approximation.For the eigenvalue bound, σ = ψ⟩⟨ψ is the density operator corresponding to the eigenstate ψ⟩ with eigenvalue λ = ⟨ψ A ψ⟩.
A.A.-G.acknowledges the generous support from Google, Inc. in the form of a Google Focused Award.This work is partly supported by the U.S. Department of Energy under Award No. DESC0019374 and the U.S. Office of Naval Research (ONS506661).A.A.-G.also acknowledges support from the Canada Industrial Research Chairs Program and the Canada 150 Research Chairs Program.A.A.-G.acknowledges generous support from Anders G. Fröseth and Sony Research.CZ and the DS3Lab gratefully acknowledge the support from the Swiss National Science Foundation (Project Number 200021_184628, and 197485), Innosuisse/SNF BRIDGE Discovery (Project Number 40B2-0_187132), European Union Horizon 2020 Research and Innovation Programme (DAPHNE, 957407), Botnar Research Centre for Child Health, Swiss Data Science Center, Alibaba, Cisco, eBay, Google Focused Research Awards, Kuaishou Inc., Oracle Labs, Zurich Insurance, and the Department of Computer Science at ETH Zurich.