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

Quantum computation in noisy, near-term implementations holds promise across multiple domains ranging from chemistry and many-body physics to machine learning, optimization, and finance. However, experimental and algorithmic shortcomings such as noise and decoherence lead to the preparation of imperfect states which differ from the ideal state and hence lead to faulty measurement outcomes, ultimately undermining the reliability of near-term quantum algorithms. It is thus crucial to accurately quantify and bound these errors via intervals which are guaranteed to contain the ideal measurement outcomes. To address this desire, based on the formulation as a semidefinite program, we derive such robustness intervals for the expectation values of quantum observables using only their first moment and the fidelity to the ideal state. In addition, to get tighter bounds, we take into account second moments and extend 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.

Quantum computation in noisy, near-term implementations holds promise across multiple domains ranging from chemistry and many-body physics to machine learning, optimization, and finance. However, experimental and algorithmic shortcomings such as noise and decoherence lead to the preparation of imperfect states which differ from the ideal state and hence lead to faulty measurement outcomes, ultimately undermining the reliability of near-term quantum algorithms. It is thus crucial to accurately quantify and bound these errors via intervals which are guaranteed to contain the ideal measurement outcomes. To address this desire, based on the formulation as a semidefinite program, we derive such robustness intervals for the expectation values of quantum observables using only their first moment and the fidelity to the ideal state. In addition, to get tighter bounds, we take into account second moments and extend 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 value 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, it is difficult to guarantee the degree of validity of the output. This unfortunate hindrance practically arises from various sources such as circuit noise and decoherence [1], lim- * aspuru@utoronto.ca † ce.zhang@inf.ethz.ch ‡ zhikuan.zhao@inf.ethz.ch ited 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 estimates of expectation values arising from imperfect approximations of a target state. To that end, based on analytical solutions to a semidefinite program (SDP), we present lower and upper bounds to quantum expectation values of observables which depend only on the fidelity with the target state and postprocessing of previously obtained measurement results. In addition, to obtain tighter bounds, we take into account higher moments of the observable by leveraging the Gramian method for pure states [12] and extend it to derive bounds which are valid for general density operators. This extension from pure to mixed states is particularly relevant in the NISQ era as it allows to derive robustness bounds which are valid in noisy scenarios. We then apply these bounds numerically and run noisy and noiseless simulations of VQE for finding ground state energies of electronic structure Hamiltonians of several molecules. Finally, we implement our results in 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
In this section, we present our main results, namely robustness intervals for quantum expectation values of observables. We first introduce the notation, summarize the technical results, and then proceed with presenting the robustness intervals.

A. Notation
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. 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 (∆A σ ) 2 = ⟨A 2 ⟩ σ − ⟨A⟩ 2 σ . We write for the trace norm of an operator A ∈ L(H). The fidelity between quantum states σ, ρ ∈ S(H) is defined as F(ρ, σ) = max ψρ, ψσ ⟨ψ ρ ψ σ ⟩ 2 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
In this work, we present lower and upper bounds on the true expectation value ⟨A⟩ σ of an observable A under a target state σ, in terms of an approximate state ρ. To that end, we employ two different techniques, each with its advantages and disadvantages in terms of efficiency and accuracy. The first technique is based on the formulation as an SDP 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], 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).
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
We present a robustness interval which is based on a formulation of the problem as a semidefinite program. This interval only requires the knowledge of the first moment of the operator A and a lower bound to the fidelity F(ρ, σ) between the target state σ and the approximation ρ. The idea is to define the problem as a semidefinite program (SDP) and then, using results from quantum hypothesis testing, derive a corresponding closed form solution. This yields upper and lower bounds on the true expectation ⟨A⟩ σ which can be expressed in terms of the fidelity F(ρ, σ) and the expectation ⟨A⟩ ρ under the approximation ρ. 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 ad-
versarial robustness of quantum classification. Formally, we present the following theorem: , the lower bound of ⟨A⟩ σ can be expressed as Similarly, for ≤ 1 2 (1 − ⟨A⟩ ρ ), the upper bound of ⟨A⟩ σ becomes In a nutshell, this result is obtained by first formulating the problem as an SDP in which the quantity ⟨Λ⟩ σ is maximized (or minimized in case of the lower bound) over all operators Λ with −1 ≤ Λ ≤ 1 and whose expectation value is equal to ⟨Λ⟩ ρ = ⟨A⟩ ρ . Formally, the upper bound can be expressed as where for the lower bound the maximization is replaced with a minimization. A closed form solution to these SDPs can then be derived by expressing it in terms of a binary quantum hypothesis testing problem. Thus, Theorem 1 provides an interval which is guaranteed to contain the true expectation value ⟨A⟩ σ of an observable A ∈ L(H) under a target state σ. In other words, it provides thresholds on the value of that expectation value when only an approximation of the state σ, the state ρ, is provided. 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⟩ ρ and whose expectation under σ saturates the bound. In Appendix A 2, we give a formal, constructive proof for this statement. Furthermore, in practice, it is typically not feasible to measure the exact expectation value ⟨A⟩ ρ due to finite sampling errors, measurement noise, and other experimental imperfections. For this reason, one needs to rely on bounds which contain the exact value with high probability. This can be accounted for in the bounds from Theorem 1 by noting that they are monotonic in ⟨A⟩ ρ , whose exact value can thus be replaced by a corresponding lower or upper bound. 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 of 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 ovelap ⟨ψ φ⟩ 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, det(G) ≥ 0. The function det(G) is a quadratic polynomial in R(⟨ψ A φ⟩) where (∆A φ ) 2 = ⟨φ A 2 φ⟩ − ⟨φ A φ⟩ 2 is the variance of A under φ⟩. This polynomial vanishes when R(⟨ψ A φ⟩) takes the values 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 and whithin this Gramian approach, a 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 this technique to mixed states.

Gramian method for mixed states
In the NISQ era, it is essential to have bounds which are valid for mixed states as these appear in noisy scenarios. 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: 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 a 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 take as starting the point Eq. (7) which is first extended to mixed state via purifications and 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 interval. 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. (3) 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 hence 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 With this, one can explicitly show that the error between each bound and the true expectation ⟨A⟩ σ scales with O( √ ) as → 0. For values of close to 1 on the other hand, we remark that both expectation value bounds tend towards the trivial bounds, namely 0 for the expectation value bound, and ±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, 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.
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 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]: In scenarios where knowledge of the lowest lying eigenvalues λ 0 and λ 1 is available, one can thus directly lowerbound the fidelity and use Eq. (11) 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 a lower degree of 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. (12) 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 be estimated via classical methods [33].
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 hold also in the degenerate case.
The fidelity bounds presented so far have in common that some degree of knowledge about the ground state in question is required. However, this degree of knowledge is not necessarily available in practice. One scenario in which this can potentially be circumvented is in the context of variational quantum imaginary time evolution (VarQITE). In general, VarQITE is a technique to find the ground state of a Hamiltonian H [34-36] by projecting the time evolution of the initial state to the evolution of the ansatz parameters. VarQITE 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-to-one 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. Thus, in this scenario, our bounds can be evaluated without requiring any knowledge on the true ground state.

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 a 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. For the expectation value bounds derived via SDP and via 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 thatH ∶= H + c1 ≥ 0 and calculate the bound forH, 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. (15) 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 →H k 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 A k = U k H k U † 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.

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], mad-ness [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 error probability 0.01, 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 scenario, 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 UpCCGSD 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 UpCCGSD. A further observation to be made is that in most cases, the lower bound is closer to the true ground state energy than the corresponding VQE, as can be seen from the error plots. The exception here   Theorem 3) for LiH (2,4)and BeH2 (4,8). The lower half in each figure shows the error between VQE and the true energy E0 computed as V QE − E0, and the error between E0 and the lower bound (LB) given by E0−LB. Here, an ideal scenario without noise is simulated and the approximation errors stem from the limited expressibility of the Ansatz state.
is BeH 2 with an UpCCGSD Ansatz where the noise level is too high for VQE to converge to a good enough approximation. Indeed, the ground state fidelities for different bond lengths vary between 0.01 and 0.1 for UpCCGSD and between 0.37 and 0.92 for the SPA Ansatz. 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 . Here we again see that in both cases the lower bound is closer to the true value than VQE, in some regions even reaching chemical accuracy where VQE does not. This is interesting, as it indicates that even though for Ansätze which lack the expressive power required to approximate the ground state, it is still possible to obtain an estimate of the ground state energy within chemical accuracy by calculating the lower bound. 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. Nevertheless, these results indicate that, in principle, it is possible to obtain very accurate lower bounds to ground state energies of electronic structure Hamiltonians.

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: 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 its formulation as an SDP, 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, to take into account further moments of A and get tighter bounds, we have revisited the Gramian method [12] 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 a given Hamiltonian. 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 . Within these experiments, we observed that the lower bound on the ground state energy is up to one order of magnitude closer to the true ground state energy compared to the respective VQE estimate, at times even reaching chemical accuracy where VQE does not. These results are promising as they indicate that our bounds not only provide sound guarantees for the performance of these algorithms, but can even be used in conjunction and deliver further accuracy im-provements.
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.

CODE AVAILABILITY
The code for the numerical simulations is publicly available at https://github.com/DS3Lab/ robustness-intervals-quantum-measurements.
Appendix A: Proof of Theorem 1

Auxiliary Lemma
Here we state a result from quantum hypothesis testing which serves as a central tool for deriving the bounds in Theorem 1. We remark that this has been shown in the appendix of Ref. [17], although in a slightly different form and for pure states. For completeness, we include the proof here and refer the reader to [17] for proofs of specific lemmas. Recall that the semidefinite program for the optimal type-II error probability in binary quantum hypothesis testing is given by The following lemma provides a lower bound to this SDP which, in the case of pure states is saturated.
Proof. It is well known that solutions to the SDP β * are attained by Helstrom operators, i.e. operators of the form Λ t = P t,+ +X t with 0 ≤ X t ≤ P t,0 where P t,0 = 1−P t,+ −P t,− and P t,± are projections onto the null, positive and negative eigenspaces of the operator σ − tρ for t ≥ 0. In particular, optimality of the Helstrom operators follows from the fact that for any 0 ≤ λ ≤ 1 we have For a formal proof of this, see Lemma 6 in the supplementary information of [17]. The idea is now to choose t and, correspondingly, X t such that Tr[(P t,+ +X t )ρ] = α 0 which, by the above implication yields β * (α 0 ; ρ, σ) = Tr[(1−Λ t )σ]. The next step is then to show that Tr[(1 − Λ t )σ] is equal to the RHS in Eq. (A2) for pure states. In the last step, this is extended to an inequality for mixed states using Uhlmann's Theorem and the monotonicity of β * under the partial trace.
In the sequel, define t 0 to be the nonnegative number This is well defined since the function t ↦ Tr[P t,+ σ] is non-increasing and continuous from the right (see Lemmas 2 & 3 in the supplementary information of [17] for a proof of this). From Lemma 5 in the supplementary information of [17] we have the following sandwich inequality which enables us to set Λ t0 ∶= P t0,+ + q 0 ⋅ P t0,0 with so that Tr[Λ t0 ρ] = α 0 as desired. It follows that β * (α 0 ; ρ, σ) = Tr[(1 − Λ t0 )σ] and it remains to show that We start with pure states and extend to mixed states in the subsequent step.
Theorem 1 (restated). 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 ≤ 1 2 (1 − ⟨A⟩ ρ ), the upper bound of ⟨A⟩ σ becomes Proof. The idea is to first formalize the problem as a semidefintie program which takes into account the first moment of A under ρ and the assumption that −1 ≤ A ≤ 1. This is then connected to the SDP Eq. (A1) where the second equality follows from the fact that replacing the equality with an inequality in the constraint of the SDP Eq. (A1) leads to the same solution. It then follows from Lemma 1 that and again use Lemma 1 and obtain the lower bound

Tightness of the SDP bound
Here we show that the SDP bounds presented in Theorem 1 are tight for pure states in the sense that, under the given assumptions, there exists an observable which saturates the bound in the case that F(ρ, σ) = 1 − . Let us first consider the upper bound and recall that in the proof of Theorem 1 we have shown where the last inequality follows from Lemma 1. Additionally, it follows from Lemma 1 that for pure states this inequality is indeed an equality, that is, we have shown that for pure states we have β * (α 0 ; ρ, σ) = α 0 (2 − 1) + (1 − ) − 2 α 0 (1 − )(1 − α 0 ) (A26) for 1− ≥ α 0 for arbitrary α 0 ∈ [0, 1]. This followed from constructing an operator Λ ⋆ with Tr[σ(1−Λ ⋆ )] = β * (α 0 ; ρ, σ). Let Λ ⋆ be such an operator for α 0 = Tr[Λ ⋆ ρ] and let A ⋆ ∶= 2Λ ⋆ − 1. We notice that which shows that the bound is indeed saturated for the observable A ⋆ . Tightness of the lower bound can be shown analogously.
and similarly Replacing A with A 2 in the above, we find Finally, we have where ⟨⋅, ⋅⟩ HS denotes the Hilbert-Schmidt inner product. Without loss of generality, we assume that ⟨Φ Ψ⟩ is real and non-negative since otherwise we can multiply each purificiation by a global phase. Applying the Gramian inequalities to the observable A ⊗ 1 and the purifications Ψ⟩, Φ⟩, we find Thus, we have shown that inequalities similar to Eq. (7) also hold for mixed states. To finish the proof, note that by assumption A ≥ 0 and hence A has a square root A = A Under the condition that F(ρ, σ) (1 − F(ρ, σ)) ≥ ∆A ρ ⟨A⟩ ρ , we can square both sides of the inequality 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.

a. Pure states
Suppose that ρ is a pure state ρ = φ⟩⟨φ . For Eckart's criterion, an analogous calculation leads to and hence so that Eckart's criterion holds in the degenerate case for this particular choice of eigenstate ψ⟩ and for pure approximation states φ⟩. Using again analogous calculations, we also obtain the extensions of the bounds Eq. (13) and Eq. (14) for pure approximation state φ⟩ in the degenerate case and for the same choice of eigenstate ψ⟩.

b. Counterexample for mixed states
If the approximation state ρ is allowed to be arbitrarily mixed, the above fidelity bounds do not hold in general. Indeed, consider the Hamiltonian for arbitrary λ, µ ∈ R with λ < µ and some arbitrary unitary U . Furthermore, let ρ be the maximally mixed state ρ = 1 3 1 3 and note that ⟨H⟩ ρ = 2λ+µ 3 . Thus, for any ψ⟩ we find that in violation of Eckart's criterion. To see that we can also construct a counterexample for the other two bounds, we calculate the variance and notice that and similarly ⟨ψ ρ ψ⟩ = 1 3