Can Error Mitigation Improve Trainability of Noisy Variational Quantum Algorithms?

Samson Wang
Theoretical Division, Los Alamos National Laboratory, Los Alamos, NM 87545, USA
Imperial College London, London, UK

Piotr Czarnik
Theoretical Division, Los Alamos National Laboratory, Los Alamos, NM 87545, USA

Andrew Arrasmith
Theoretical Division, Los Alamos National Laboratory, Los Alamos, NM 87545, USA
Quantum Science Center, Oak Ridge, TN 37931, USA

M. Cerezo
Theoretical Division, Los Alamos National Laboratory, Los Alamos, NM 87545, USA
Center for Nonlinear Studies, Los Alamos National Laboratory, Los Alamos, NM, USA

Quantum Science Center, Oak Ridge, TN 37931, USA

Lukasz Cincio
Theoretical Division, Los Alamos National Laboratory, Los Alamos, NM 87545, USA
Quantum Science Center, Oak Ridge, TN 37931, USA

Patrick J. Coles
Theoretical Division, Los Alamos National Laboratory, Los Alamos, NM 87545, USA
Quantum Science Center, Oak Ridge, TN 37931, USA

Variational Quantum Algorithms (VQAs) are widely viewed as the best hope for near-term quantum advantage. However, recent studies have shown that noise can severely limit the trainability of VQAs, e.g., by exponentially flattening the cost landscape and suppressing the magnitudes of cost gradients. Error Mitigation (EM) shows promise in reducing the impact of noise on near-term devices. Thus, it is natural to ask whether EM can improve the trainability of VQAs. In this work, we first show that, for a broad class of EM strategies, exponential cost concentration cannot be resolved without committing exponential resources elsewhere. This class of strategies includes as special cases Zero Noise Extrapolation, Virtual Distillation, Probabilistic Error Cancellation, and Clifford Data Regression. Second, we perform analytical and numerical analysis of these EM protocols, and we find that some of them (e.g., Virtual Distillation) can make it harder to resolve cost function values compared to running no EM at all. As a positive result, we do find numerical evidence that Clifford Data Regression (CDR) can aid the training process in certain settings where cost concentration is not too severe. Our results show that care should be taken in applying EM protocols as they can either worsen or not improve trainability. On the other hand, our positive results for CDR highlight the possibility of engineering error mitigation methods to improve trainability.

§ INTRODUCTION

The prospect of obtaining quantum computational advantage for practical problems, such as simulating systems in chemistry and materials science, has generated much excitement. The past few years have witnessed tremendous progress towards this end, with significant focus on algorithm development for Noisy Intermediate-Scale Quantum (NISQ) computers. In particular, Variational Quantum Algorithms (VQAs) are a leading algorithmic approach because they adapt to the constraints of NISQ devices. Specifically, VQAs minimize a cost function by training a parameterized quantum circuit via a classical-quantum feedback loop <cit.>. The cost is computed efficiently on a quantum computer whilst the parameter optimization is carried out classically. Different implementations of this versatile framework have been proposed for a broad spectrum of problems from dynamical quantum simulation <cit.> to machine learning <cit.> and beyond <cit.>.

A central challenge in the NISQ regime is to combat the effects of noise as full error correction is not possible <cit.>. Decoherence, gate errors, and measurement noise all conspire to limit the complexity of quantum circuits that can be implemented on NISQ devices. While VQAs themselves offer some strategy to mitigate the impact of noise <cit.>, it is widely viewed that VQAs alone will not be enough, and additional strategies will be needed to obtain quantum advantage in the face of noise. This has spawned the field of error mitigation (EM), and many researchers believe that VQAs combined with EM techniques will be the path forward. Indeed, EM methods like Zero-Noise Extrapolation <cit.>, Clifford Data Regression <cit.>, Virtual Distillation <cit.>, Probabilistic Error Cancellation <cit.> and others <cit.> have been demonstrated to reduce errors of observable expectation values, sometimes by orders of magnitude. Hence, there has been hope that one can simply train the VQA in the presence of noise, and then after training, one can apply an EM method to extract the correct cost value (e.g., the ground state energy in the case of the variational quantum eigensolver <cit.>).

However, new challenges have recently been discovered for this approach <cit.>. It is now recognized that noise impacts the trainability of VQAs, that is, the ability of the classical optimizer to find the global cost minimum. For ansatzes (i.e., parameterized quantum circuits) with depth linear or superlinear in the number of qubits and local Pauli noise, the cost function landscape exponentially flattens, leading to an exponentially vanishing cost gradient, a phenomenon known as Noise-Induced Barren Plateaus (NIBPs) <cit.>. Thus, noise impedes the training process of VQAs, as in such a setting one requires an exponential number of shots per optimization step to resolve the cost landscape against finite sampling noise. As with other barren plateau effects <cit.>,  this exponential scaling does not only arise for gradient-based optimizers but also impacts gradient-free methods <cit.> and optimizers that use higher-order derivatives <cit.>. NIBPs represent a serious issue for VQA scalability, and could ultimately be a roadblock for near-term quantum advantage. It is therefore crucial to investigate potential methods to mitigate them.

Given the great success of EM methods in suppressing error in observable expectation values, it is natural to ask whether EM methods could address NIBPs. More generally, one could simply ask: does it help to use error mitigation during the training process for VQAs? This question is precisely the topic of our article. We remark that error mitigation has been successfully implemented during the VQA training process for a small-scale problem <cit.>. However, it is an open question as to whether or not EM can resolve large-scale trainability issues associated with cost concentration. This is due to the fact that even though EM can reverse the concentration of cost values, it also increases the stochastic uncertainty in the mitigated quantities. As summarized in Figure <ref>, this is a trade-off that should be carefully considered. Thus, it is a non-trivial question as to whether or not EM improves the resolvability of cost function values which is a key factor in determining the trainability of the landscape.

In this work, we investigate the effects of error mitigation on the resolvability of the cost function landscape. First, we consider a broad class of error mitigation protocols and show that, under the class of local depolarizing noise that is known to cause NIBPs, in order to reverse exponential cost concentration any such protocol needs to spend resources (e.g., shot resources or number of state copies) scaling at least exponentially in the number of qubits. This suggests that NIBPs are a serious scaling issue that cannot be simply resolved with error mitigation.

Error mitigation can impair the resolvability of cost function landscapes. (a): A central primitive in training VQAs is the task of comparing two cost function values ($C(\thv_1 )$ and $C(\thv_2 )$) on the cost landscape in parameter space. Ideally (with infinite sampling), these cost values correspond to the mean values of some probability distributions (left panel). However, in an experimental setup, one only has a finite shot budget and by collecting measurement statistics one obtains an estimate of the mean values by sampling from these distributions (right panel). (b): The effect of certain types of noise models is to concentrate cost function values. This impedes trainability as any two cost function values ($\widetilde{C}(\thv_1 )$ and $\widetilde{C}(\thv_2 )$)  have small separation and require many shots to accurately distinguish. (c): Error mitigation can mitigate many effects of noise and potentially recover key features of the noise-free cost function. In an ideal scenario, the separation of the mitigated cost values ($C_m(\thv_1 )$ and $C_m(\thv_2 )$) closely resembles that of the noise-free landscape. However, the caveat is that the variance of statistical outcomes can increase greatly. The effect of this is that the two cost function points can often require even more shots to resolve accurately, compared to the unmitigated case.

Error mitigation can impair the resolvability of cost function landscapes. (a): A central primitive in training variational quantum algorithms is the task to compare two cost function values on the cost landscape in parameter space. These cost values correspond to mean values of probability distributions (left panel). In a experimental setup, one has a finite shot budget and by collecting measurement statistics one obtains an estimate of the mean values by sampling from these distributions (right panel). (b): The effect of certain types of noise models is to concentrate cost function values. This impedes trainability as any two cost function values have small separation and require many shots to accurately distinguish. (c): Error mitigation can mitigate many effects of noise and potentially recover many features of the noise-free cost function. In an ideal scenario, the separation of the mitigated cost values closely resembles that of the noise-free landscape. However, the caveat is that the variance of statistical outcomes can increase greatly. The effect of this is that the two cost function points can often require even more shots to resolve accurately, compared to the unmitigated case.

Second, we study four specific error mitigation protocols in further detail: Zero Noise Extrapolation, Virtual Distillation, Probabilistic Error Cancellation, and strategies that implement a linear ansatz which includes Clifford Data Regression. We find that Virtual Distillation can actually decrease the resolvability of the noisy cost landscape, and impede trainability. Under more restrictive assumptions on the cost landscape, we find a similar result for Zero Noise Extrapolation. We also show that any improvement in the resolvability after applying Probabilistic Error Cancellation under local depolarizing noise exponentially degrades with increasing number of qubits. Finally, for strategies that use a linear ansatz such as Clifford Data Regression, we show that there is no change to the resolvability of any pair of cost values if the same ansatz is used. However, we do observe numerically that Clifford Data Regression increases trainability in some settings. This last observation provides some hope that a careful choice of error mitigation method can be useful. It also suggests that researchers could design and engineer error mitigation methods to enhance VQA trainability.

The rest of the manuscript is structured as follows. Section <ref> introduces the framework and notation for our work. We present our theoretical results in Section <ref> and our numerical results in Section  <ref>. Finally, our concluding discussions are presented in Section <ref>. The proofs for our main results are presented in the Appendix.

§ FRAMEWORK

§.§ Variational Quantum Algorithms

The main goal of Variational Quantum Algorithms (VQAs) is to solve an optimization problem by minimizing a cost function that can be efficiently estimated on a quantum computer. In this work we consider settings where the cost function takes the form

\begin{equation}\label{eq:costfn}
C(\boldsymbol{\theta}) = \Tr\left[U(\thv)\rho_{in} U\ad(\thv)O\right]\,.
\end{equation}

In the above, given some Hilbert space $\HC$, we define the set of density operators $\SC(\HC)$ and set of bounded linear operators $\BC(\HC)$. We then denote $\rho_{in} \in \SC(\HC)$ as the input state, $U(\thv)\in \BC(\HC)$ as a unitary that corresponds to a parametrized quantum circuit with trainable parameters $\thv$, and $O\in \BC(\HC)$ is a Hermitian operator. The Variational Quantum Eigensolver <cit.>, variational quantum compiling <cit.>, quantum autoencoders <cit.>, and several other VQAs fit under the framework of Eq. (<ref>).

A quantum computer is employed to evaluate the cost function, or gradients thereof, and part of the computational complexity of the algorithm is designated to a classical computer that leverages the power of classical optimizers to solve the problem

\begin{equation}\label{eq:optimization}
\argmin_{\thv}    C(\boldsymbol{\theta})\,.
\end{equation}

The optimization task defined in Eq. (<ref>) has been shown to be NP-hard <cit.>. Moreover, on top of the typical difficulties associated with solving classical non-convex optimization problems, there are challenges that arise when training the parameters of a VQA due to the quantum nature of the problem itself.

As quantum mechanics is intrinsically  a probabilistic theory, one has to deal with shot noise arising from finite sampling when estimating the cost function (or its gradient). This has led to the development of several quantum-aware optimizers that are frugal in the number of shots <cit.>.  Additionally,  it has been recently shown that certain properties of the cost function can induce so-called barren plateaus, originating due to highly expressive ansatzes <cit.>, global cost functions <cit.>, high levels of entanglement <cit.>, or the controllability of $U(\thv)$ <cit.>. When a cost function exhibits a barren plateau, with high probability the cost function partial derivatives are exponentially suppressed across the landscape. This means that an exponentially large number of shots are needed to navigate the flat landscape and determine a cost-minimizing direction <cit.>.

In this work we investigate the effect of noise and error mitigation techniques in solving the optimization task of Eq. (<ref>). For this purpose we investigate the task of resolving two points on the cost function landscape, as presented in Fig. <ref>. This is a central primitive in the training process that is utilized at each optimization step, regardless of whether one is using gradient-based or gradient-free methods. In gradient-based methods, a common strategy is to use the parameter shift rule, which constructs partial derivatives from two cost function values <cit.>. Gradient-free methods such as simplex-based methods also compare two or more cost function values at each optimization step <cit.>. Thus, this task is a key step for both gradient-based and gradient-free optimizers, and it reflects the ability of the optimizer to find a cost-minimizing direction at each step of the optimization.  As discussed below, under a finite shot budget this task becomes harder under cost concentration, leading to trainability issues.

§.§ Effect of noise on the training landscape

Hardware noise can impact the cost function landscape in a variety of ways such as changing the optimal cost function value, shifting the position of minima, and demoting a global minimum to a local minimum. All of the above present further challenges in the training of VQAs. In this section we briefly review some of the literature on the effect of noise on VQAs cost function landscapes. We summarize some of these effects in Fig. <ref>.

§.§.§ Noise resilience

Certain cost functions have been demonstrated to show optimal parameter resilience under particular noise models <cit.>. This is a phenomenon where the position of the global cost minimum of the cost landscape is invariant under the action of noise. This has important consequences for trainability. There are many VQAs where the goal is to obtain optimal parameters, rather than the optimal cost value, such as when solving combinatorial optimization problems with the Quantum Approximate Optimization Algorithm <cit.>. If such cost landscapes display optimal parameter resilience, this leaves open the possibility of noisy training even if the cost value of the global minimum is altered by the noise. However, noise can also severely affect the trainability of the landscape in a number of ways, which we summarize below.

§.§.§ Noise-induced cost concentration and noise-induced barren plateaus

Here we summarize the phenomenon of noise-induced cost concentration and noise-induced barren plateaus (NIBPs), as well as introduce some notation that we will use throughout the rest of this manuscript.
This was formulated in Ref. <cit.> for a general class of VQAs and a class of Pauli noise that includes as a special case local depolarizing noise.

(See also Refs. <cit.> for other discussions of the impact of noise.)
In this section we follow the treatment for local depolarizing noise.
Consider a model of noise acting through a depth $L$ circuit with $n$-qubit input state $\rho_{in}$ as

\begin{equation} \label{eq:noisystate}
\widetilde{\rho} = (\mathcal{N}\circ \mathcal{U}_L\circ \cdots \circ \mathcal{N}\circ \mathcal{U}_1\circ \mathcal{N}\big) (\rho_{in})
\end{equation}

where $\{\UC_k\}^L_{k=1}$ denote unitary channels that describe collections of gates that act together in a layer, and $\NC=\bigotimes_{i=1}^n \DC_i$ is an instance of local depolarizing channels, with depolarizing probabilities $\{p_i \}_{i=1}^n$. We denote a noisy cost function as

\begin{equation}
\widetilde{C} = \Tr\left[O\widetilde{\rho}\right]\,,
\end{equation}

where $O$ is some Hermitian measurement operator (throughout the article we will use a tilde to denote noisy quantities). In Ref. <cit.> it was shown that

\begin{equation}\label{eq:costconcentration}
\Big\vert \widetilde{C} - \frac{1}{2^n}\Tr[O]\Big\vert\, \leq\, D(q,n) \,,
\end{equation}

where $q=\max\{(1-p_i)\}_{i=1}^n <1$ and $D(q,n)\in \mathcal{O}(q^{\alpha n})$ for some positive constant $\alpha$ if $L \in \Omega(n)$. Thus, in the presence of the class of noise models considered, the noisy cost function exponentially concentrates on a fixed value if the depth scales linearly or superlinearly in the number of qubits.

The gradients across the cost function landscape show similar scaling <cit.>, demonstrating a phenomenon known as NIBPs. This implies that the task of accurately determining gradients or cost function differences during the training process requires an exponential number of shots due to the need to resolve quantities to an exponentially small precision.

Schematic of different effects due to noise on cost landscapes. We present a 1-dimensional slice of a simplified cost landscape corresponding to a single parameter $\theta$. a) Depending on the parameterization strategy, some ansatzes can have degenerate minima. b) Certain types of local Pauli noise can cause the cost landscape to exponentially concentrate on a fixed value. Some can problems display optimal parameter resilience (OPR), where the location of the optimal parameters are invariant under action of the certain noise models. c) Aside from cost concentration, noise can also corrupt the cost landscape by breaking the degeneracy of optimal parameters, and shifting the location of minima.

§.§.§ Cost corruption

In general, a noise model that exhibits cost concentration and NIBPs would not simply uniformly flatten the cost landscape. Instead, we expect noise to additionally alter the cost landscape in many non-trivial ways. We refer to any additional adverse effects on the landscape as cost corruption.
For example, it was shown in Ref. <cit.> that non-unital noise can break the degeneracy of exponentially-occurring global minima, thus proliferating local minima and impacting trainability. In addition, cost functions that do not exhibit optimal parameter resilience <cit.> limit the quality of noisy optimization, as the optimal parameters of $\widetilde{C}(\thv)$ do not correspond to the optimal parameters of $C(\thv)$.

§.§ Error Mitigation Techniques

We finish the discussion of our framework with a summary of the key features of the error mitigation techniques that we study in this article.

For a more detailed review, readers can refer to Refs. <cit.>.

Consider the effects of noise on the cost function in Eq. (<ref>). We suppose the noise can be characterized by a single parameter $\varepsilon$ and we denote the corresponding noisy state and cost function as $\widetilde{\rho}(\boldsymbol{\theta},\varepsilon)$ and $\widetilde{C}(\boldsymbol{\theta},\varepsilon) = \Tr[\widetilde{\rho}(\boldsymbol{\theta},\varepsilon)O]$ respectively. The goal of error mitigation is to construct an experimental protocol which obtains a mitigated cost function estimator $C_m(\boldsymbol{\theta})$ that approximates the noise-free value $C(\boldsymbol{\theta})$. The protocol to obtain $C_m(\boldsymbol{\theta})$ generally consists of running circuits that modify the original circuit of interest by inserting additional gates, preparing multiple copies of a state, changing the measurement operator, and classical post-processing of the expectation values of these circuits. These different utilizations of resources are summarized in a schematic in Fig. <ref>.

Error mitigation protocols often lead to a larger variance in the statistical outcomes of each experiment, and thus more shots are required to estimate the error-mitigated cost value $C_m(\boldsymbol{\theta},\varepsilon)$ to a desired precision compared to the unmitigated noisy value $\widetilde{C}(\boldsymbol{\theta},\varepsilon)$. This is often quantified by the error mitigation cost, which is defined below.

We define the error mitigation cost as

\begin{equation}
\gamma (\boldsymbol{\theta},\varepsilon) = \frac{\Var[{C}_m(\boldsymbol{\theta},\varepsilon)]}{\Var[\widetilde{C} (\boldsymbol{\theta},\varepsilon)]}\,,
\end{equation}

where $\widetilde{C} (\boldsymbol{\theta},\varepsilon)$ denotes the noisy cost function value corresponding to vector of parameters $\boldsymbol{\theta}$ at noise level $\varepsilon$, and ${C}_m(\boldsymbol{\theta},\varepsilon)$ denotes the corresponding error-mitigated quantity.

In certain settings we encounter, $\gamma(\boldsymbol{\theta},\varepsilon)$ is independent of $\boldsymbol{\theta}$. In other cases, where we have to compare $\gamma(\boldsymbol{\theta},\varepsilon)$ for two different parameters, we will seek parameter-independent bounds.

We now summarize the error mitigation techniques that we study in this article. We note that recently, unified error mitigation techniques have also been proposed that combine two or more of the protocols that we discuss in this section <cit.>. Our results are also applicable to such strategies, however, we will only review the root strategies here.

Schematic of resource use in error mitigation. Noise is indicated by the shaded orange region. (a) Cost function values are obtained by taking input state $\rho_{in}$, applying parameterized gates which we denote as a unitary channel $\mathcal{U}(\thv)$, and measuring the resulting state $\mathcal{U}(\thv)(\rho_{in})$ with observable $O$. (b) Noise can corrupt the gates in the circuit and the measurement process. (c) Error mitigation aims to obtain a good approximation to the noise-free cost $C(\thv)$ by employing a number of strategies such as: modifying the gates implemented $\mathcal{U}(\thv) \rightarrow \mathcal{V}(\thv)$ or the input state $\rho_{in} \rightarrow \sigma_{in}$, utilizing multiple copies of the quantum circuit, modifying the measurement operator $O \rightarrow X$, and utilizing clean ancillary qubits at the end of the circuit. Many such circuits with different hyperparameters can be run, with their expectation values combined in a post-processing step, in order to construct the final error mitigated cost value $C_m(\thv)$.
One feature that distinguishes the approaches to error mitigation studied here from error correction is that error correction allows global access to the larger Hilbert space from the start of the circuit, whereas the error mitigation techniques studied here only allow the possibility for global operations at the end of the circuit.

§.§.§ Zero Noise Extrapolation

The goal of Zero Noise Extrapolation is to run a given circuit of interest at $m+1$ increasing noise levels $\varepsilon < a_1\varepsilon < ... < a_m\varepsilon$, and to use information from the resulting expectation values to obtain an estimate of the zero-noise result. Here we summarize the key features of a protocol using Richardson extrapolation <cit.>, and exponential extrapolation <cit.>.

Richardson Extrapolation. Suppose that $\widetilde{C}(\boldsymbol{\theta}_i,\varepsilon)$ admits a Taylor expansion in small noise parameter $\varepsilon$ as

\begin{equation} \label{eq:ZNE_taylor}
\widetilde{C}(\boldsymbol{\theta}_i,\varepsilon)=\widetilde{C}(\boldsymbol{\theta}_i,0)+\sum_{k=1}^{m} p_{k}(\boldsymbol{\theta}_i) \varepsilon^{k}+\mathcal{O}(\varepsilon^{m+1})\,,
\end{equation}

where $p_k$ are unknown parameters and $\widetilde{C}(\boldsymbol{\theta}_i,0)=C(\thv)$ is the zero-noise cost function. By considering the equivalent expansion of $\widetilde{C}(\boldsymbol{\theta}_i,a_1 \varepsilon)$ and combining the two equations one obtains

\begin{equation}\label{eq:ZNE_TaylorCm}
C_m(\boldsymbol{\theta}_i) = \frac{a_1 \widetilde{C}(\boldsymbol{\theta}_i,\varepsilon) - \widetilde{C}(\boldsymbol{\theta}_i,a_1 \varepsilon) }{ a_1-1} = \widetilde{C}(\boldsymbol{\theta}_i,0) +\mathcal{O}(\varepsilon^2)\,,
\end{equation}

which is a higher-order approximation of $\widetilde{C}(\boldsymbol{\theta}_i,0)$ compared to simply using $\widetilde{C}(\boldsymbol{\theta}_i,\varepsilon)$. This process can be repeated iteratively $m$ times to obtain an estimator which is accurate up to $\mathcal{O}(\varepsilon^{m+1})$ error.

Exponential extrapolation. In some cases the noisy behavior may not be well-depicted by a Taylor expansion. As an alternative one can consider an exponential model

\begin{equation}\label{eq:ZNE_exp}
\widetilde{C}(\boldsymbol{\theta}_i,\varepsilon)= r(\boldsymbol{\theta}_i,\varepsilon)^{-t(\boldsymbol{\theta}_i,\varepsilon)} \left(
\sum_{k=0}^{m} p_{k}(\boldsymbol{\theta}_i) \varepsilon^{k}+\mathcal{O}(\varepsilon^{m+1}) \right)\,,
\end{equation}

for some $r$ and $t$ which in general can be functions of $\varepsilon$. We can also construct an extrapolation strategy that is tailored towards noisy cost function values that are dominated by NIBP scaling as in Eq. (<ref>), where we model the effects of noise as

\begin{align}\label{eq:ZNE_NIBP}
\widetilde{C}(\boldsymbol{\theta}_i,q)&= \widetilde{C}(\boldsymbol{\theta}_i,0) + q^{L} \bigg( B(\boldsymbol{\theta}_i) +\sum_{k=1}^{m} p_{k} (1-q)^{k} \\
& \qquad\qquad\qquad\qquad\quad +\mathcal{O}\left((1-q)^{m+1}\right) \bigg)\,.
\end{align}

Here, $\widetilde{C}(\boldsymbol{\theta}_i,0)$ is the fixed point of the noise and $B(\boldsymbol{\theta}_i)+\widetilde{C}(\boldsymbol{\theta}_i,0)$ is the noise free cost value.
For these two strategies we can similarly construct $C_m(\boldsymbol{\theta}_i)$ as linear combinations of $\{\widetilde{C}(\thv_i,a_i\varepsilon) \}_{i=0}^{m}$ to achieve $\mathcal{O}(\varepsilon^{m+1})$ approximations of $\widetilde{C}(\boldsymbol{\theta}_i,0)$. We detail these constructions in Section <ref> of the Appendix.

§.§.§ Virtual Distillation

Virtual Distillation, also known as Error
Suppression by Derangement, was proposed concurrently in Refs. <cit.> and <cit.>. In this article we consider the two error mitigation protocols in Ref. <cit.> (denoted "A" and "B") to respectively prepare

\begin{equation}\label{eq:VD_A}
C_m^{(A)}(\boldsymbol{\theta}_i)=\Tr[\widetilde{\rho}_i^M O]/\Tr[\widetilde{\rho}_i^M]\,,
\end{equation}

\begin{equation}\label{eq:VD_B}
C_m^{(B)}(\boldsymbol{\theta}_i) = \Tr[\widetilde{\rho}_i^M O]/\lambda_i^M\,,
\end{equation}

where $\lambda_i$ is the dominant eigenvalue of $\widetilde{\rho}_i \equiv \widetilde{\rho}(\boldsymbol{\theta}_i)$. The operator $\widetilde{\rho}_i^M$ can be obtained by preparing $M$ copies of $\widetilde{\rho}_i$ in a tensor product state $\widetilde{\rho}_i^{\otimes M}$ and applying a cyclic shift operator. We note that protocol B presumes access to the dominant eigenvalue beforehand, which could potentially be computed via the techniques of Ref. <cit.>.

§.§.§ Probabilistic Error Cancellation

Probabilisitic Error Cancellation utilizes many modified circuit runs in order to construct a quasiprobability representation of the noise-free cost function <cit.>. We assume that the effect of the noise can be described by a quantum channel $\mathcal{N}$ that occurs after a gate that we denote with unitary channel $\mathcal{U}$. Here we make the simplifying assumption that this is the only gate in the circuit, and we treat the general case in Section <ref> of the Appendix, as well as provide a more detailed exposition. The goal of this protocol is to simulate the inverse map $\mathcal{N}^{-1}$. Note that, in general, this will not always correspond to a CPTP map. Despite this fact, if one has a basis of (noisy) quantum channels $\{\mathcal{B}_\alpha \}_{\alpha}$, corresponding to experimentally available channels, one can expand the inverse map in this basis as $\mathcal{N}^{-1} = \sum_{\alpha} q_{\alpha}  \mathcal{B}_{\alpha}$,  for some set of $q_{\alpha} \in \mathbb{R}$. By defining a probability distribution $p_{\alpha} = |q_{\alpha}|/G_\mathcal{N}$ where $G_\mathcal{N} = \sum_{\alpha} |q_{\alpha}|$, the noise free expectation value can then be written as a quasiprobability distribution

\begin{align}
{C}_{\mathcal{U}(\rho)} &= G_\mathcal{N} \sum_{\alpha} \mathrm{sgn}(q_{\alpha})\, p_{\alpha}\, \Tr \big[ \BC_{\alpha}\NC\UC(\rho_{in})O\big] \,,
\end{align}

where $\rho_{in}$ is the input state, $O$ is the measurement operator, and $\mathrm{sgn}(q_{\alpha})$ denotes the sign of $q_{\alpha}$. The idea is that if one has access to the set of CPTP maps $\{ \mathcal{B}_{\alpha}\}_\alpha$ in the noisy native hardware gate set, then one can obtain an estimate of the noise free cost ${C}_{\mathcal{U}(\rho)}$ as follows: (1) With probability $p_\alpha$, prepare the circuit of interest with additional gate $\mathcal{B}_\alpha$ in order to obtain the expectation value $\Tr \left[ \BC_{\alpha}\NC\UC(\rho_{in})O\right]$. (2) Multiply the result by $\mathrm{sgn}(q_{\alpha})G_\mathcal{N}$. (3) Repeat process many times and sum results.

§.§.§ Clifford Data Regression (CDR) and linear ansatz methods

The main idea of linear ansatz methods is to assume that we can approximately reverse the effects of noise with an affine map, and thus we construct a linear ansatz of the form

\begin{equation}\label{eq:CDR_ansatz}
C_m(\thv,\boldsymbol{a}) = a_1(\thv)\widetilde{C}(\thv) + a_2(\thv)\,,
\end{equation}

where $\boldsymbol{a}(\thv) = (a_1(\thv),a_2(\thv))$ is a vector of parameters to be determined.  In general we expect $\boldsymbol{a}$ to be highly dependent on $\boldsymbol{\theta}$. In Ref. <cit.>, the authors use data regression to learn the optimal parameters $\vec{a}^*(\thv)$ with training data comprising of pairs of noise-free and corresponding noisy cost function values $\mathcal{T}_{\boldsymbol{\theta}} = \{(C_j, \widetilde{C}_j )\}_j$, where the circuits are predominantly constructed from Clifford gates. The noise-free cost values can be simulated efficiently on a classical computer whilst the noisy cost values can be evaluated directly on the quantum computer. This strategy is known as Clifford Data Regression.

Other methods have been proposed to learn the optimal parameters $\vec{a}^*(\thv)$. In Ref. <cit.> the authors further develop the idea of training-based error mitigation by considering alternative training data comprising of fermionic linear optics circuits. One can also model the noise as global depolarizing noise.

Under this assumption, $\vec{a}^*(\thv)$ has an exact solution in terms of a single noise parameter.
Subsequently, various techniques can be used to estimate the noise parameter <cit.>.

§ THEORETICAL RESULTS

We present two sets of theoretical results.
First, in Section <ref> we show that a broad class of error mitigation techniques cannot undo the exponential resource requirement that exponential cost concentration presents. This has implications for both the trainability of noisy VQAs, as well as the accurate estimation of noise-free cost function values in general. Second, in Section <ref>, we work predominantly in the non-asymptotic  regime (in terms of scaling in $n$) and investigate to what extent different error mitigation strategies can improve the resolvability of the noisy cost landscape, assuming that some cost concentration has occurred. For these purposes we introduce a class of quantities which quantify the improvement of the resolvability of the cost function landscape after error mitigation, which we call the relative resolvability (see Defs. <ref>-<ref>). Using these quantities we study Zero Noise Extrapolation (Sec. <ref>), Virtual Distillation (Sec. <ref>), Probabilistic Error Cancellation (Sec. <ref>) and linear ansatz methods which include Clifford Data Regression (Sec. <ref>).  In the settings that we consider, we find that in many cases error mitigation impedes the optimizer's ability to find good optimization steps, and is worse than performing no error mitigation.

§.§ Asymptotic scaling results (exponential estimator concentration)

In this section we show that full mitigation of exponential cost concentration is not possible for a general class of error mitigation strategies. Specifically, we show that one cannot remove the exponential scaling that local depolarizing noise incurs without investing exponential resources elsewhere in the mitigation protocol.

We start by remarking that, as summarized in Fig. <ref>, all of the strategies presented in Sec. <ref> consist of preparing linear combinations of expectation values of the form

\begin{align}\label{eq:thm1_estimator}
E_{\sigma,X,M,k} = \Tr\left[X \left(\sigma^{\otimes M}\otimes \ket{0}\bra{0}^{\otimes k}\right)\right]\,,
\end{align}

for some  $n$-qubit quantum state $\sigma \in S(\mathcal{H})$ that in general can be prepared by a different circuit to that of the state of interest, for $\ket{0}\bra{0} \in S(\mathcal{H}')$ and for some $X\in B(\mathcal{H}^{\otimes M}\otimes \mathcal{H}'^{\otimes k})$. That is, one can prepare multiple copies of a state, prepare different quantum circuits, and apply general measurement operators. In order to generalize the setting further, we also allow the possibility to ultilize multiple clean ancillary qubits at the end of the circuit. By considering linear combinations of such quantities, one also accounts for the ability to post-processing measurement results classically with a linear map, such as is the case with Probabilistic Error Cancellation. In the following theorem we show how quantities of the form (<ref>) concentrate under local depolarizing noise.

Consider an error mitigation strategy that, as a step in its protocol, estimates $E_{\sigma,X,M,k}$ as defined in Eq. (<ref>).
Suppose that $\sigma$ is prepared with a depth $L_{\sigma}$ circuit and experiences local depolarizing noise according to Eq. (<ref>). Under these conditions, $E_{\sigma,X,M,k}$ exponentially concentrates with increasing circuit depth on a state-independent fixed point as

\begin{align}
\left\vert E_{\sigma,X,M,k} - \Tr\left[X \left( \frac{\emph{\id}^{\otimes M}}{2^{Mn}} \otimes \ket{0}\bra{0}^{\otimes k}\right)\right]\right\vert\, &\leq\, G_{\sigma,X,M}(n) \,, \label{eq:prop1_conc}
\end{align}

where $\emph{\id} \in S(\mathcal{H})$ is the $n$-qubit identity operator and

\begin{equation}
G_{\sigma,X,M}(n) =  \sqrt{\ln 4}\,\big\|X\big\|_\infty M n^{1/2} q^{L_{\sigma}+1}\,,
\end{equation}

with noise parameter $q\in[0,1)$.

We remark that a similar result can be obtained for local Pauli noise, provided that each Pauli error occurs with a non-zero probability. Theorem <ref> shows that quantities of the form (<ref>) exponentially concentrate in the depth of the circuit. As we summarize in the schematic in Fig. <ref>, such quantities generalize expectation values that are prepared by many different error mitigation protocols. We now explicitly demonstrate how Theorem <ref> affects the mitigated cost values that these protocols output.

Consider an error mitigation protocol that approximates the noise-free cost value $C(\thv)$ by estimating the quantity

\begin{equation}\label{eq:cor1_estimator}
C_m(\thv) = \sum_{(\sigma(\boldsymbol{\theta}),X,M,k)\in T} a_{X,M,k}E_{\sigma(\boldsymbol{\theta}),X,M,k}\,,
\end{equation}

where each $E_{\sigma,X,M,k}$ takes the form (<ref>). We denote $M_{max}$ and $a_{max}$ as the maximum values of $M$ and $a_{X,M,k}$ respectively accessible from a set $T$ defined by the given protocol. Assuming $\|X\|_\infty \in \OC(\poly(n))$, there exists a fixed point $F$ independent of $\boldsymbol{\theta}$ such that

\begin{equation}
\big|C_m(\boldsymbol{\theta})-F\big|\in\OC(2^{-\beta n}a_{max}|T|M_{max})\,,
\end{equation}

for some constant $\beta\geq 1$ if the circuit depths satisfy

\begin{equation}
L_{\sigma(\boldsymbol{\theta})} \in \Omega(n)\,,
\end{equation}

for all $\sigma(\boldsymbol{\theta})$ in the construction (<ref>). That is, if the depths of the circuits scale linearly or superlinearly in $n$ then one requires at least exponential resources to distinguish $C_m$ from its fixed point, for instance by requiring an exponential number of shots, or by requiring an exponential number of state copies $M_{max}$.

We note that the assumption $\|X\|_\infty \in \OC(\poly(n))$ is satisfied in most settings, and in particular is satisfied for all error mitigation protocols discussed in Sec. <ref>. For instance, in the case of Virtual Distillation, $X$ corresponds to a cyclic shift operator followed by a Pauli observable, and thus $\|X\|_\infty \in \OC(1)$. Corollary <ref> implies that under conditions that generate a NIBP, in order to obtain an estimate of a noise-free cost value up to some arbitrary additive error, one requires resource consumption that scales exponentially in the number of qubits.  In Appendix <ref> we present a more detailed statement that explains how such resources may be consumed.

Whilst the use of clean ancillary qubits as part of an error mitigation protocol, utilized as in Equation (<ref>) and Fig. <ref>, has not been widely studied, Corollary <ref> rules out the possibility that such resources used at the end of the circuit would offer advantage in countering the exponential scaling effects due to cost concentration. Indeed, upon inspecting (<ref>), the ancilla appear explicitly in the form of the fixed point. This highlights a key difference between many error mitigation strategies and error correction, as error correction utilizes resources (such as a larger Hilbert space) in the middle of the computation, whilst the error mitigation protocols considered here are based on processing states obtained at the end of a noisy computation. Our result leaves open the possibility that novel error mitigation protocols that move beyond the framework of (<ref>) and Fig. <ref> can have hope of countering the exponential scaling of exponential cost concentration and NIBPs.

§.§ Non-asymptotic protocol-specific results

In this section we present predominantly non-asymptotic results for Zero-Noise Extrapolation (Sec. <ref>), Virtual Distillation (Sec. <ref>), Probabilistic Error Cancellation (Sec. <ref>), and methods which use a linear ansatz such as Clifford Data Regression (Sec. <ref>). For each protocol, we
investigate the effect of error mitigation on the resolvability of the cost landscape, for different classes of noisy states. To the end, we first define a class of resolvability measures which quantify how many shots it takes to resolve the cost landscape at some fixed precision after applying error mitigation, compared to no mitigation at all.

§.§.§ Definitions

Consider two locations in parameter space $\boldsymbol{\theta}_{1}$, $\boldsymbol{\theta}_{2}$ and their corresponding points on the cost landscape. Denote the number of shots to resolve these two points up to some fixed precision with and without error mitigation as $N_\text{EM}$ and $N_\text{noisy}$ respectively. We define the relative resolvability for $\boldsymbol{\theta}_{1}$ and $\boldsymbol{\theta}_{2}$ at error level $\varepsilon$ as

\begin{align}\label{eq:partial_EM}
\rchi({\boldsymbol{\theta}_{1,2}},\varepsilon) &= \frac{N_\text{noisy}({\boldsymbol{\theta}_{1,2}},\varepsilon)}{N_\text{EM}({\boldsymbol{\theta}_{1,2}},\varepsilon)} \\
&= \frac{1}{\gamma(\varepsilon)}\, \left(\frac{\Delta C_m({\boldsymbol{\theta}_{1,2}},\varepsilon)}{\Delta \widetilde{C}({\boldsymbol{\theta}_{1,2}},\varepsilon)}\right)^2 \,,
\end{align}

where we have used the shorthand notation $\rchi({\boldsymbol{\theta}_{1,2}}) = \rchi({\boldsymbol{\theta}_{1}},{\boldsymbol{\theta}_{2}})$, $\gamma$ is the error mitigation cost as defined in Definition $\ref{def:em_cost}$, and where we denote

\begin{align}
\Delta \widetilde{C}({\boldsymbol{\theta}_{1,2}},\varepsilon) &= \widetilde{C}({\boldsymbol{\theta}_{1}},\varepsilon) - \widetilde{C}({\boldsymbol{\theta}_{2}},\varepsilon)\,, \\
\Delta {C}_m({\boldsymbol{\theta}_{1,2}},\varepsilon) &= {C}_m({\boldsymbol{\theta}_{1}},\varepsilon) - {C}_m({\boldsymbol{\theta}_{2}},\varepsilon)\,.
\end{align}

We can see that if we have $\rchi({\boldsymbol{\theta}_{1,2}},\varepsilon) > 1$, then error mitigation has succesfully increased the resolvability of the cost values corresponding to the cost values at $\thv_1$ and $\thv_2$. Note that this criterion is a necessary but not sufficient condition for error mitigation to reverse the effects of cost concentration on the cost landscape. Namely, it does not require the mitigated landscape to accurately reflect the noise-free landscape, and it does not account for other trainability issues such as proliferation of minima. If $\rchi({\boldsymbol{\theta}_{1,2}},\varepsilon) < 1$ then ${N_\text{noisy}({\boldsymbol{\theta}_{1,2}},\varepsilon)}<{N_\text{EM}({\boldsymbol{\theta}_{1,2}},\varepsilon)} $. This implies that error mitigation has exacerbated the resolvability issues associated with cost concentration and NIBPs, and it has been counterproductive in fixing these trainability issues.

For a general cost functions, the relative resolvability of cost function points after mitigation may vary significantly across the landscape, or be different for different choices of ansatzes and noise models. This motivates us to seek averaged measures of resolvability. We consider two types of averaging: first, an average over cost function points generated by a given ansatz, noise and cost; second, an average over a set of noisy states.

Denote the vector of parameters that corresponds to the global cost minimum at noise parameter $\varepsilon$ as $\boldsymbol{\theta}_{*}$. We then define the averaged relative resolvability as

\begin{equation} \label{eq:av_resolvability}
\overline{\rchi}(\varepsilon) = \frac{1}{\gamma(\varepsilon)}\, \left( \frac{\langle\Delta C_m({\boldsymbol{\theta}_{i,*}},\varepsilon)\rangle_i}{\langle\Delta \widetilde{C}({\boldsymbol{\theta}_{i,*}},\varepsilon)\rangle_i} \right)^2\,,
\end{equation}

where $\langle\cdot\rangle_i$ denotes the mean over all parameter vectors $\boldsymbol{\theta}_{i}$ accessible with the given ansatz of consideration, and where we denote

\begin{align}
\Delta {C}_m({\boldsymbol{\theta}_{i,*}},\varepsilon) &= {C}_m({\boldsymbol{\theta}_{i}},\varepsilon) - {C}_m({\boldsymbol{\theta}_{*}},\varepsilon)\,,\\
\Delta \widetilde{C}({\boldsymbol{\theta}_{i,*}},\varepsilon) &= \widetilde{C}({\boldsymbol{\theta}_{i}},\varepsilon) - \widetilde{C}({\boldsymbol{\theta}_{*}},\varepsilon)\,.
\end{align}

Averaging across a given cost landscape gives a result that is particular to the choice of ansatz, measurement operator and noise model. In order to evaluate the performance of error mitigation in a more general setting, we consider a broader average over noisy states that have the same spectrum. This choice of class of noisy states is motivated by the fact that a central mechanism of cost concentration and NIBPs under unital Pauli noise is the loss of purity. When we consider the second averaged relative resolvability for Virtual Distillation in Section <ref>, it will turn out to be bounded by a function of the purity for such states.

Consider a  normalized spectrum $\boldsymbol{\lambda} \in \mathbb{R}_+^{2^n}$ which corresponds to the eigenspectrum of some noisy state. We define the 2-design-averaged relative resolvability as

\begin{equation}
\dbloverline{\rchi}_{\boldsymbol{\lambda}} = \frac{1}{\gamma(\boldsymbol{\lambda})} \frac{\big\langle\big(C_m(\rho_{\boldsymbol{\lambda}},U_i) - \Tr[O]/2^n \big)^2\big\rangle_{U_i}}{\big\langle\big(\widetilde{C}(\rho_{\boldsymbol{\lambda}},U_i) - \Tr[O]/2^n \big)^2\big\rangle_{U_i} }\,,
\end{equation}

where $\langle\cdot\rangle_{U_i}$ denotes an average over unitaries $U_i$ drawn from a unitary 2-design, $\rho_{\boldsymbol{\lambda}}$ is an arbitrarily chosen reference state with spectrum $\boldsymbol{\lambda}$, and where we denote

\begin{align}
\widetilde{C}(\rho_{\boldsymbol{\lambda}},U_i) &= \Tr[U_i\rho_{\boldsymbol{\lambda}}U^\dag_i O]\,, \\
{C}_m(\rho_{\boldsymbol{\lambda}},U_i) &= \Tr[\mathcal{M}(U_i\rho_{\boldsymbol{\lambda}}U^\dag_i) O]\,,
\end{align}

where $\mathcal{M}: S(\mathcal{H})\mapsto B(\mathcal{H})$ is a map that describes the action of the error mitigation protocol.

In Fig. <ref> we present a a schematic of our second averaged relative resolvability.

Schematic for average relative resolvability II. In Definition <ref> we consider a broader averaged resolvability measure where the average is taken over a class of noisy states, rather than a cost landscape generated by a particular ansatz. This is constructed from the following game: (1) Prepare a reference noisy state $\rho_{\boldsymbol{\lambda}}$ with spectrum $\boldsymbol{\lambda}$ and conjugate by unitary $U_i$ drawn from a 2-design. (2) Pass the resulting state through the considered error mitigation protocol, and evaluate the resolvability from the fixed point of the noise. (3) Do the same, without error mitigation. (4) Average over the 2-design and compare the averaged resolvabilities.

In the results that follow we do not constrain ourselves to investigating specific models of Pauli noise. Instead, we simply suppose that there exists a noise model that causes some concentration of the cost function onto some state-independent fixed point. Any further assumptions on the effects of the noise model are explicitly specified in the statement of the particular result in question. As a precursor to looking at more generic models of noisy states, we will also find it useful to investigate the performance of error mitigation strategies in the presence of global depolarizing noise of the form

\begin{equation}\label{eq:gdepol_noise}
\rho \xrightarrow{\mathcal{D}} \widetilde{\rho} = (1-p) \rho + p \frac{\id}{2^n}\,,
\end{equation}

where $\mathcal{D}$ is the global depolarizing channel and $p$ is the depolarizing probability. Our justification for studying this noise model is twofold. First, global depolarizing noise provides a clean model of cost concentration with no other cost corrupting effects of the noise. Therefore, if a given error mitigation strategy is to mitigate the effects of cost concentration and NIBPs, we expect it to be able to perform well on this noise model. Second, the structure of many error mitigation strategies is directly motivated by the model of global depolarizing noise <cit.>. Indeed, many such strategies have been shown to achieve good or perfect performance with this noise model in mitigating noisy cost function values <cit.>. However, we stress that trainability may simultaneously get worse, which is what we will now investigate.

§.§.§ Zero Noise Extrapolation

In this section we present our results on Zero Noise Extrapolation. First, we consider the simple model of global depolarizing noise.

Consider a circuit with $L$ instances of global depolarizing noise of the form Eq. (<ref>).
Consider a Richardson extrapolation strategy based on Eq. (<ref>), an exponential extrapolation strategy based on Eq. (<ref>) and a NIBP extrapolation strategy based on Eq. (<ref>). We presume access to an augmented noisy circuit where the error probability is exactly increased by factor $a_1>1$ as $p\rightarrow a_1p$. Then, we have

\begin{align}
&{\rchi_{depol}}\leq \frac{\left(c-\frac{(1-a_1p)^L}{(1-p)^L}\right)^2}{c^2+1}\,,
\end{align}

where $\rchi_{depol}$ is the relative resolvability (see Definition <ref>) for global depolarizing noise, and where

\begin{equation}
c = \begin{cases}
\quad a_1\quad &\text{for Richardson extrapolation,} \\
\frac{a_1r(\varepsilon)^{t(\varepsilon)}}{r(a_1\varepsilon)^{t(a_1\varepsilon)}} \quad &\text{for exponential extrapolation,} \\
\; a_1^{-(L+1)} \quad &\text{for NIBP extrapolation}\,.
\end{cases}
\end{equation}

Thus, $\rchi_{depol}\leq 1$ for all of the above extrapolation strategies with access to 2 noise levels.

We see that for all the above techniques, Zero Noise Extrapolation with access to 2 noise levels decreases the resolvability of the cost function under global depolarizing noise. Further, if one attempts to directly reverse the exponential scaling of NIBPs that global depolarizing noise incurs, one obtains an exponentially worse relative resolvability. We now consider how resolvability behaves under Zero Noise Extrapolation on average across the cost landscape, given a generic noise model.

Consider a Richardson extrapolation strategy based on Eq. (<ref>), an exponential extrapolation strategy based on Eq. (<ref>) and a NIBP extrapolation strategy based on Eq. (<ref>). We presume perfect access to an augmented noisy circuit where the noise rate is increased by factor $a_1>1$. We denote $\boldsymbol{\theta}_{\varepsilon*}$ as the parameter corresponding to the global cost minimum at base noise parameter $\varepsilon$.
Further denote $ \frac{\langle\Delta\widetilde{C}(\thv_{i,\varepsilon*}, a_1\varepsilon)\rangle_{i}}{\langle\Delta\widetilde{C}(\thv_{i,\varepsilon*},\varepsilon)\rangle_{i}} = z$.
Any such noise model has an average relative resolvability

\begin{equation}
\overline{\rchi} \leq  \frac{(z-c)^2}{c^2+1}\,,
\end{equation}

\begin{equation}
c = \begin{cases}
\quad a_1\quad &\text{for Richardson extrapolation,} \\
\frac{a_1r(\varepsilon)^{t(\varepsilon)}}{r(a_1\varepsilon)^{t(a_1\varepsilon)}} \quad &\text{for exponential extrapolation,} \\
\; a_1^{-(L+1)} \quad &\text{for NIBP extrapolation}\,.
\end{cases}
\end{equation}

Thus, under the assumption that $z\leq1$ and $\langle\Delta \widetilde{C}(\thv_{i,\varepsilon*},a_1\varepsilon)\rangle_i\geq0$, $\overline{\rchi}\leq 1$ for all of the above extrapolation strategies with access to 2 noise levels.

Proposition <ref> shows that under mild assumptions of the effect of the noise on the cost landscape, Zero Noise Extrapolation with access to 2 noise levels impairs the resolvability of the cost landscape. These assumptions have physical meaning: $z\leq1$ implies that on average the cost concentrates when the noise parameter is boosted, whilst $\langle\Delta \widetilde{C}(\thv_{i,\varepsilon*},a_1\varepsilon)\rangle\geq0$ implies that the landscape is not heavily corrupted after boosting the noise parameter so that the global minimum at the base noise level remains below the average cost value. We also see that in the presence of exponential cost concentration and NIBPs, the relative resolvability is exponentially small if one attempts to directly reverse the exponential scaling of NIBPs.

In the Appendix, we study a modification of the averaged resolvability in Definition <ref> and find that this is bounded by a function of the purity of the noisy states, such that the resolvability decreases if purity decreases with increasing noise level.
This result, along with the proofs of the above propositions, can be found in Appendix <ref>.
Finally, we remark that in the above results we consider a scenario where the Richardson, exponential or NIBP extrapolation strategies utilize expectation values from only two noise levels. In Appendix <ref> we show that similar results may be obtained for Richardson extrapolation with access to 3 distinct noise levels.

§.§.§ Virtual Distillation

Here we present our results on Virtual Distillation. In the following proposition we start again with the simple model of global depolarizing noise.

Consider global depolarizing noise of the form in Eq. (<ref>)
acting on some $n$-qubit pure state $\rho$ with error probability $p$. We consider the two error mitigation protocols of Ref. <cit.> (denoted "A" and "B") to respectively prepare (<ref>) and (<ref>), using $M$ copies of a quantum state. The relative resolvabilities to resolve any two arbitrary cost function points satisfy

\begin{align}
\rchi_{depol}^{(A)} \leq \rchi_{depol}^{(B)} = \Gamma(n,M,p)\,,
\end{align}

for all $n\geq1$, $M\geq2$, $p \in [0,1]$, and where

\begin{equation}
\Gamma(n,M,p)\leq 1 \,,
\end{equation}

is a monotonically decreasing function in $M$ (with asymptotically exponential decay) for $n\geq 1$, $M\geq 2$. Within this region the bound is saturated as $\Gamma(1,2,p)= 1$ for all $p$.

Proposition <ref> shows that Virtual Distillation decreases the resolvability of cost landscapes suffering from global depolarizing noise. Moreover, as the number of state copies $M$ increases, the effect worsens. We find similar results in the following proposition when considering averaged resolvabilities over a class of noisy states.

Consider an error mitigation protocol that prepares estimator $C_m(\boldsymbol{\theta}_i) = \Tr[\widetilde{\rho}_i^M O]/\Tr[\widetilde{\rho}_i^M]$ from some noisy parameterized quantum state $\widetilde{\rho}_i \equiv \widetilde{\rho}(\boldsymbol{\theta}_i)$. Consider the average relative resolvability $\dbloverline{\rchi}_{\boldsymbol{\lambda}}$ for noisy states of some spectrum $\boldsymbol{\lambda}$ with purity $P$ as defined in Definition <ref>. We have

\begin{equation}
\dbloverline{\rchi}_{\boldsymbol{\lambda}}\leq G(n,M,P)\,\leq\,1\,,
\end{equation}

where $G(n,M,P)$ is a monotonically decreasing function in $M$ (with asymptotically exponential decay) for all $n\geq 1$, $M\geq 2$.
Within this region the bound is saturated as $G(1,2,P)=1$ for all $P$, and as $G(n,M,1)=1$ for all $n\geq 1, M\geq 2$.

We present the explicit forms of $\Gamma(n,M,p)$ and $G(n,M,P)$ as well as a proof of the above propositions in Sec. <ref> of the Appendix. We note that in most cases the bound given by Proposition <ref> decreases with decreasing noisy state purity. This indicates that within such settings, the greater the loss of purity due to noise, the worse the impact on resolvability is after error mitigation with Virtual Distillation.

§.§.§ Probabilistic Error Cancellation

Here we present our results for Probabilistic Error Cancellation. We utilize the optimal quasiprobability decompositions studied in Ref. <cit.>, and the proofs can be found in Section <ref> of the Appendix.

Consider a quasi-probability method that corrects global depolarizing noise of the form (<ref>). For any pair of states corresponding to points on the cost function landscape, the optimal quasiprobability scheme gives

\begin{equation}
\rchi_{depol} = \frac{2^{2n}}{2^{2n}-p(2-p)} \geq 1\,,
\end{equation}

for all $n\,\geq1,\,  p\in[0,1]$, which is achieved with access to noisy Pauli gates.

Proposition <ref> shows that for the special case of global depolarizing noise, Probabilistic Error Cancellation actually improves the resolvability of the noisy cost landscape. However, this improvement is generally small and is decreasing quickly with the number of qubits $n$. For instance, for $n=1$, $\rchi_{depol}$ has maximum value $4/3$ (achieved in the limit of maximum depolarization probability). For for $n=2$, $\rchi_{depol}$ has maximum value $\approx 1.07$. In the limit of large $n$, $\rchi_{depol}$ tends to 1.

We extend this study to local depolarizing noise in Appendix <ref>.

We find that for a single instance of local depolarizing noise, Probabilistic Error Correction can either improve resolvability or worsen it, depending on the strength of concentration of the cost.
In addition, we show in the following proposition that if one wishes to mitigate all the noisy gates in the circuit and one has NIBP scaling, the improvement due to Probabilistic Error Cancellation degrades exponentially, and ultimately for large problem sizes this impairs resolvability.

Consider local depolarizing noise with depolarizing probability $p$ to act in $L$ layers through a depth $L$ circuit as in Eq. (<ref>). Suppose that the effect of this noise is to cause cost concentration

\begin{equation}\label{eq:QP_localscaling-assumption}
\langle\Delta \widetilde{C}(\boldsymbol{\theta}_{i,*})\rangle_i = Aq^L \langle\Delta C(\boldsymbol{\theta}_{i,*})\rangle_i\,,
\end{equation}

for some constant $A$ and noise parameter $q\in [0,1)$. The optimal quasiprobability method to mitigate the depolarizing noise in the circuit yields

\begin{equation}
\overline{\rchi} = \frac{1}{A^2q^{2L}}\left(Q(p) \right)^{nL} \,,
\end{equation}

where $Q(p)=1 - \frac{3p(2-p)}{4-p(2-p)} \in [0,1)$ for $p\in (0,1]$. Thus, the average relative resolvability has unfavourable scaling with system size.

We note that (<ref>) gives the best possible scaling of noisy cost differences allowed by (<ref>) under local depolarizing noise. Thus, Proposition <ref> shows that if has NIBP scaling under local depolarizing noise, one still has unfavourable scaling for resolvability.

§.§.§ Linear ansatz methods

In Proposition <ref> we consider a scenario where the same linear ansatz (<ref>) is applied to two points on the noisy cost landscape. For Clifford Data Regression this is a reasonable assumption in scenarios where one is comparing two points that are close in parameter space, for instance, when a simplex-based optimizer is exploring a small local region. However, we remark this is not always true in general settings.

Consider any error mitigation strategy that mitigates noisy cost function value $\widetilde{C}(\boldsymbol{\theta})$ by constructing an estimator $C_m(\boldsymbol{\theta})$ of the form
(<ref>). For any two noisy cost function points to which the same ansatz is applied, we have

\begin{equation}
\rchi = 1\,,
\end{equation}

for any noise process.

Under global depolarizing noise, the optimal linear ansatz gives $\rchi = 1$ for any pair of cost function points.

Corollary <ref> comes simply by noting that the optimal choice of linear ansatz under global depolarizing noise corrects the noise exactly and is state independent <cit.>. The above results imply that in some settings CDR has a neutral effect on the resolvability of the cost function landscape. This opens up the possibility that in practical settings CDR can improve the trainability of cost landscapes, if it can remedy other cost corrupting effects due to noise outside of cost concentration. This motivates our numerical studies of CDR, which we present in the following section.

§ NUMERICAL RESULTS

As discussed in Sec. <ref>, in many settings, current state-of-the-art error mitigation methods do not mitigate the effects of cost concentration. Nevertheless, as discussed in Sec. <ref>, trainability of VQAs is also affected by other cost-corrupting effects. We expect that error mitigation can reverse some of the effects due to cost corruption that affect the trainability of VQAs when the effects of cost concentration are not too severe.
In this section, we numerically investigate the effects of error mitigation on trainability in such a setting to provide possible evidence towards beneficial effects of error mitigation. To this end, we focus on CDR as in some settings it does not worsen the effects of cost concentration, as shown in Sec. <ref>.

We perform our numerical experiments by simulating the Quantum Approximate Optimization Algorithm (QAOA) <cit.> for 5-qubit MaxCut problems.
We use a realistic noise model of an IBM quantum computer <cit.>, which has been obtained by gate set tomography of IBM's Ourense quantum device. Further, we assume linear connectivity of the simulated quantum computer.

A MaxCut problem is defined for a graph $G=(V,E)$ of nodes $V$ and edges $E$. The problem is to find a bipartition of the nodes into two sets which maximizes the number of edges connecting the sets. This problem can be reformulated as finding the ground state of a  Hamiltonian

\begin{equation}
H_{\rm MaxCut}= -\frac{1}{2} \sum_{ij\in E}   (\mathds{1} - Z_i Z_j),
\end{equation}

where $Z_i,Z_j$ are Pauli $Z$ matrices.
Here we consider graphs with $n=5$ vertices, with $36$ randomly generated instances according to the Erdös-Rényi model <cit.>, where for each pair of vertices in the graph there is a connecting edge with probability $0.5$.

To  approximate the ground state of $H_{\rm MaxCut}$ we simulate the QAOA for number of rounds $p=1$ to $8$. The QAOA ansatz applied to the input state is given as

\begin{equation}
\prod_{j=p,p-1\dots,1} e^{i \beta_j H_M} e^{i \gamma_j H_{\rm MaxCut}}  (|+\rangle)^{\otimes n},
\end{equation}

where $H_M =   \sum_{j} X_j$, $X_j$ are Pauli $X$ matrices, we denote  $\ket{+}=(\ket{0}+\ket{1})/\sqrt{2}  $ , and $\beta_j$, $\gamma_j$ are variational parameters. We minimize the cost function $\langle H_{\rm MaxCut} \rangle$ using the Nelder-Mead algorithm <cit.>. We perform the optimization with shot budgets ranging from $N_{\rm tot}=10^7$ to $1.5\times 10^8$. We define $N_{\rm tot}$ as total number of shots spent on the optimization. We detail the optimization procedure in Appendix <ref>.  In our numerics, the values of $N_{\rm tot}$ are chosen to enable implementation of the optimization with current quantum computers.   To quantify the quality of the noisy or CDR-mitigated optimization we compute approximation ratios of the solutions using the exact expectation value of $\langle H_{\rm MaxCut} \rangle$. The approximation ratio is defined here as the ratio of a given solution's energy to the true ground state energy.

We gather our numerical results for CDR in Fig. <ref>.  In the figure we plot the approximation ratio averaged over $36$ randomly chosen graphs versus $N_{\rm tot}$.  We compare the quality of the solutions of noisy (unmitigated) and CDR-mitigated optimization and find that  CDR-mitigated optimization outperforms noisy optimization for all considered $p$ and  $N_{\rm tot}$ values. We observe that the solutions for  $p=2$  outperform those for $p=1$ for both CDR-mitigated and noisy optimization. The quality of $p>2$ solutions decline with increasing $p$ for noisy optimization, while it remains approximately the same for $p=2$ to $6$ for CDR-mitigated optimization. With CDR-mitigated optimization we see a decrease in quality of solution for the largest considered $p=8$.

The numerical results presented here are obtained for circuits shallow enough to be trainable while using the CDR-mitigated cost function. Therefore they are outside of the NIBP scaling regime.  As discussed in Section <ref>, even outside the NIBP regime noise may adversely impact trainability by corrupting the cost function landscape, which error mitigation has a chance to remedy.

Our results give hope that  CDR-mitigated optimization  may overall offer a trainability advantage for problems with such cost function landscape corruption.

As  discussed in Section <ref>  optimizing  an error mitigated cost function is not guaranteed to  outperform its noisy optimization even outside the NIBP regime. Indeed, we find numerically that for $p=2,4$ optimization with Virtual Distillation does not outperform noisy optimization for the considered implementations (see Appendix <ref>).

Comparing CDR-mitigated and noisy optimization for  $\boldsymbol{5}$ qubit Max-Cut QAOA. We plot the approximation ratio of solutions for noisy (red circles)  and CDR-mitigated  (blue diamonds) optimization of Max-Cut QAOA for $5$ qubits. Different panels show results for different numbers of QAOA rounds $p$ plotted versus total number of shots $N_{\rm tot}$ spent on the optimization of a MaxCut problem. Here, we compute the approximation ratios using exact $H_{\rm MaxCut}$ energies  to benchmark quality of the noisy and CDR mitigated optimization.  For each $p$ we average  the approximation ratio over $36$ MaxCut graphs chosen randomly   from   Erdös-Rényi ensemble. For all  $p$ and $N_{\rm tot}$ values we see an advantage of the CDR-mitigated optimization over noisy optimization.

§ DISCUSSION

Noise can exponentially degrade the trainability of linear (or superlinear) depth Variational Quantum Algorithms (VQAs) by flattening the cost landscape, thus requiring an exponential precision in system size to resolve its features <cit.>. This limits the scope for achieving possible quantum advantage with VQAs. At present there are no known strategies to avoid this exponential scaling completely aside from pursuing algorithms with sublinear circuit depth, and current strategies to mitigate this effect consist only of reducing hardware noise rates. Thus, is it a pressing challenge to search for possible solutions to this problem. Error mitigation strategies emerge as a natural candidate to tackle this problem under near-term constraints.

In this work we investigate the effects of error mitigation on the trainability of noisy cost function landscapes in two regimes. First, we work in the asymptotic regime (in terms of scaling with system size) and find that if a VQA is suffering from exponential cost concentration, requiring an exponential number of shots to accurately resolve cost values, then a broad class of error mitigation strategies (including as special cases Zero Noise Extrapolation, Virtual Distillation, Probabilistic Error Cancellation, Clifford Data Regression) cannot remove this exponential scaling. Within the considered paradigm, this exponential scaling implies that at least an exponential number of resources needs to be spent in order to extract accurate information from the cost landscape in order to find a cost-minimizing optimization direction. In Corollary <ref> we identify circuit samples (or shots) as well as number of copies of a quantum state as two such resources.

Second, we move out of the asymptotic regime and investigate whether or not particular error mitigation protocols can improve the resolvability of noisy cost landscapes. Should such a landscape be burdened with exponential cost concentration, this would correspond to an improvement in the coefficient in the exponential scaling. Our results indicate that some error mitigation protocols can worsen the resolvability, and ultimately the trainability, of cost landscapes in certain settings. In particular, in Propositions <ref> and <ref> we show analytically that Virtual Distillation impairs resolvability with worsening resolvability as the number of state copies increases. We obtain similar results for Zero Noise Extrapolation in Propositions <ref> and <ref> under some assumptions of the cost landscape. Numerical analysis of a particular MaxCut problem indicates that trainability is overall similarly impaired for Virtual Distillation.

Clifford Data Regression (CDR) distinguishes itself from the other error mitigation techniques considered in this article, as in contrast to the other protocols it does not necessarily increase the statistical uncertainty of cost values more than it reverses their concentration. This is reflected in the fact that under a global depolarizing noise model, CDR has neutral impact on resolvability (Corollary <ref>). However, it is also known that CDR can remedy the effects of more complex noise models. This indicates that CDR could resolve trainability issues arising due to corruptions of the cost function outside of cost concentration, whilst having a neutral effect on cost concentration itself, and thus overall improve trainability. In the numerical example studied, presented in Fig. <ref>, we observe this to be the case. This points to deeper future work studying the mechanisms that allow error mitigation to improve the trainability of noisy cost landscapes.

Finally, we identify that the broad class of error mitigation protocols we study in our asymptotic analysis all only consist of post-processing expectation values of noisy circuits, as summarized in Fig. <ref>. This gives intuition as to why they cannot escape the exponential scaling of noise-induced barren plateaus (NIBPs). However, the theory of error correction indicates that with sufficient resources NIBPs can indeed be avoided. This gives hope that there can exist novel error mitigation strategies that move beyond the framework of the protocols considered in this article and thereby avoid the exponential impairment to trainability that NIBPs present.

§ ACKNOWLEDGEMENTS

This work was supported by the Quantum Science Center (QSC), a National Quantum Information Science Research Center of the U.S. Department of Energy (DOE). SW was supported by the U.S. DOE, Office of Science, Office of Advanced Scientific Computing Research, under the Quantum Computing Application Teams program. SW was also partially supported by the Samsung GRP grant. Piotr C. and AA were supported by the Laboratory Directed Research and Development (LDRD) program of LANL under project numbers 20190659PRD4 (Piotr C.) and 20210116DR (AA and Piotr. C). MC acknowledge support from the Center for Nonlinear Studies at Los Alamos National Laboratory. MC and LC were also initially supported by the LDRD program of LANL under project number 20190065DR. PJC also acknowledges initial support from the LANL ASC Beyond Moore's Law project. This research used resources provided by the Los Alamos National Laboratory Institutional Computing Program, which is supported by the U.S. Department of Energy National Nuclear Security Administration under Contract No. 89233218CNA000001.

0pt Appendix :

Road map of appendices. In Appendix
<ref> we present some notation and definitions that we need in order to prove our main results, as well as provide further details on the error mitigation protocols studied in this article. In Appendix <ref> we derive some useful lemmas that are required for our proofs. In Appendix <ref> we present the proof for our asymptotic results on the exponential concentration of estimators. In Appendix <ref> we present our protocol-specific results on the change in resolvability of the cost landscape under error mitigation. Finally, in Appendix  <ref> we discuss details of our numerical implementations for Clifford Data Regression and Virtual Distillation.

§ PRELIMINARIES

§.§ Further details on error mitigation techniques

In this section we expand on our discussion in Section <ref> and provide further details on the Zero Noise Extrapolation and Probabilistic Error Cancellation protocols.

§.§.§ Zero Noise Extrapolation (ZNE)

For convenience in this section we recall the key points of Zero Noise Extrapolation as summarized in Section (<ref>). We also detail the explicit forms of the estimators that can be constructed for exponential extrapolation and an extrapolation strategy tailored towards NIBP effects, which will be required in order to prove our results.

Richardson Extrapolation. We suppose that $\widetilde{C}(\boldsymbol{\theta}_i,\varepsilon)$ admits a Taylor expansion in small noise parameter $\varepsilon$ as

\begin{equation} \label{eq:appdx_ZNE_taylor}
\widetilde{C}(\boldsymbol{\theta}_i,\varepsilon)=\widetilde{C}(\boldsymbol{\theta}_i,0)+\sum_{k=1}^{m} p_{k}(\boldsymbol{\theta}_i) \varepsilon^{k}+\mathcal{O}(\varepsilon^{m+1})\,,
\end{equation}

where $p_k$ are unknown parameters and $\widetilde{C}(\boldsymbol{\theta}_i,0)$ is the zero-noise cost function. By considering the equivalent expansion of $\widetilde{C}(\boldsymbol{\theta}_i,a_1 \varepsilon)$ and combining the two equations we can obtain

\begin{equation}\label{eq:appdx_ZNE_TaylorCm}
C^{(2)}_m(\boldsymbol{\theta}_i) = \frac{a_1 \widetilde{C}(\boldsymbol{\theta}_i,\varepsilon) - \widetilde{C}(\boldsymbol{\theta}_i,a_1 \varepsilon) }{ a_1-1} = \widetilde{C}(\boldsymbol{\theta}_i,0) +\mathcal{O}(\varepsilon^2)\,,
\end{equation}

which is a higher-order approximation of $\widetilde{C}(\boldsymbol{\theta}_i,0)$ compared to simply using $\widetilde{C}(\boldsymbol{\theta}_i,\varepsilon)$. This process can be repeated iteratively $m$ times to obtain an estimator which is accurate up to $\mathcal{O}(\varepsilon^{m+1})$ error. It can be shown that the general form for the estimator that uses $k$ noise levels can be written as

\begin{equation}\label{eq:appdx_ZNE_general_estimator}
C^{(k)}_m(\boldsymbol{\theta}_i) = \sum^k_{j=0} \beta_j \widetilde{C}(\thv_i,a_j\varepsilon)\,,
\end{equation}

where the coefficients $\beta_j$ satisfy the linear system of equations $\sum_{j=0}^k\beta_j=1$ and $\sum_{l=0}^k\beta_la_l^t=0$ for all $t\in \{1,...,k \}$ <cit.>.
For 3 noise levels, (<ref>) explicitly gives

\begin{equation}\label{eq:appdx_ZNE_TaylorCm3}
C^{(3)}_m(\boldsymbol{\theta}_i) = \frac{a_1a_2(a_2-a_1) \widetilde{C}(\boldsymbol{\theta}_i,\varepsilon) - a_2(a_2-1)\widetilde{C}(\boldsymbol{\theta}_i, a_1\varepsilon) + a_1(a_1-1)\widetilde{C}(\boldsymbol{\theta}_i, a_2\varepsilon) }{(a_1-1)(a_2-1)(a_2-a_1)} = \widetilde{C}(\boldsymbol{\theta}_i,0) +\mathcal{O}(\varepsilon^3)\,.
\end{equation}

Exponential extrapolation. We can also consider an exponential model

\begin{equation}\label{eq:appdx_ZNE_exp}
\widetilde{C}(\boldsymbol{\theta}_i,\varepsilon)= r(\boldsymbol{\theta}_i,\varepsilon)^{-t(\boldsymbol{\theta}_i,\varepsilon)} \left(
\sum_{k=0}^{m} p_{k}(\boldsymbol{\theta}_i) \varepsilon^{k}+\mathcal{O}(\varepsilon^{m+1}) \right)\,,
\end{equation}

for some $r$ and $t$, which in general can be functions of $\varepsilon$. Following in similar steps to Richardson extrapolation we can consider the same expansion at an augmented noise level $a_1 \varepsilon$. This enables us to construct the estimator

\begin{align}\label{eq:appdx_ZNE_expCm}
C_m(\boldsymbol{\theta}_i) = \frac{1}{a_1-1} \Big( &a_1 r(\boldsymbol{\theta}_i,\varepsilon)^{t(\boldsymbol{\theta}_i,\varepsilon)} \widetilde{C}(\boldsymbol{\theta}_i,\varepsilon) \nonumber \\
&-r(\boldsymbol{\theta}_i,a_1\varepsilon)^{t(\boldsymbol{\theta}_i,a_1\varepsilon)} \widetilde{C}(\boldsymbol{\theta}_i,a_1 \varepsilon) \Big)\,,
\end{align}

which approximates $\widetilde{C}(\boldsymbol{\theta}_i,0)$ to a higher order in $\varepsilon$ compared to $\widetilde{C}(\boldsymbol{\theta}_i,\varepsilon)$.

NIBP extrapolation. We can construct an alternative Zero Noise Extrapolation strategy that is tailored towards noisy cost function values that are dominated by NIBP scaling as in Eq. (<ref>). We model the effects of noise as

\begin{equation}\label{eq:appdx_ZNE_NIBP}
\widetilde{C}(\boldsymbol{\theta}_i,q)= \widetilde{C}(\boldsymbol{\theta}_i,0) + q^{L} \left( B(\boldsymbol{\theta}_i) +\sum_{k=1}^{m} p_{k} (1-q)^{k}+\mathcal{O}\left((1-q)^{m+1}\right) \right)\,,
\end{equation}

for all noisy cost function points $\widetilde{C}(\boldsymbol{\theta}_i,q)$, where $\widetilde{C}(\boldsymbol{\theta}_i,0)$ is the fixed point of the noise and $B(\boldsymbol{\theta}_i)+\widetilde{C}(\boldsymbol{\theta}_i,0)$ is the noise free cost value. (Note that for NIBPs we cannot consider lower-order polynomials of $q$, as else the NIBP condition would be broken for small $q$.)
We construct estimators for any given parameter $\boldsymbol{\theta}_i$, as

\begin{equation}\label{eq:appdx_ZNE_Cm_NIBP}
{C}_m(\boldsymbol{\theta}_i) = \frac{ a^{L+1}q^{-L} (\widetilde{C}(\boldsymbol{\theta}_i,q/a)-\widetilde{C}(\boldsymbol{\theta}_i,0)) -  q^{-L}(\widetilde{C}(\boldsymbol{\theta}_i,q)-\widetilde{C}(\boldsymbol{\theta}_i,0))}{a-1} + K\,,
\end{equation}

where $\widetilde{C}(\boldsymbol{\theta}_i,0)$ is the parameter-independent fixed value obtained at maximum noise level, and $K=\widetilde{C}(\boldsymbol{\theta}_i,0)-\sum_k p_k$ is an additive constant. As we are only interested in cost function differences for our results, this will cancel out.

§.§.§ Probabilistic Error Cancellation

General idea. Probabilisitic Error Cancellation utilizes many modified circuit runs in order to construct a quasiprobability representation of the noise-free cost function <cit.>. We assume that the effect of the noise can be described by a quantum channel $\mathcal{N}$ that occurs after a gate that we denote with unitary channel $\mathcal{U}$. For now we assume this is the only instance of noise in the circuit, however, we will later generalize to many instances of noise. The goal of Probabilistic Error Cancellation is to simulate the inverse map $\mathcal{N}^{-1}$. Note in general this will not always correspond to a CPTP map. Despite this, if we have a basis of (noisy) quantum channels $\{\mathcal{B}_\alpha \}_{\alpha}$, corresponding to experimentally available channels, we can expand the inverse map in this basis as

\begin{equation}
\mathcal{N}^{-1} = \sum_{\alpha} q_{\alpha}  \mathcal{B}_{\alpha} \,,
\end{equation}

for some set of $q_{\alpha} \in \mathbb{R}$. Then, the channel that describes the noiseless gate can be written as

\begin{align}
\mathcal{U} &= \NC^{-1}\circ \NC\circ \UC \\
&=\sum_{\alpha} q_{\alpha} \mathcal{K}_{\alpha} \,,\label{eq:appdx_QuasiNoiselessgate}
\end{align}

where we have defined $\mathcal{K}_{\alpha} = \mathcal{B}_{\alpha}\circ\mathcal{N}\circ\mathcal{U}$. Denote the input state to the gate as $\rho_{in}$ and a measurement operator as $O$. The expectation value can be written

\begin{equation}\label{eq:appdx_QuasiEV}
{C}_{\mathcal{U}(\rho)} = \Tr \left[ \mathcal{U}(\rho_{in})O\right] = \sum_{\alpha} q_{\alpha} \widetilde{C}_{\mathcal{K}_{\alpha}(\rho)} \,,
\end{equation}

where for simplicity we first assume that $\mathcal{U}$ is the only gate in the circuit, and $\widetilde{C}_{\mathcal{K}_{\alpha}(\rho)} \equiv \Tr \left[ \mathcal{K}_{\alpha}(\rho)O\right]$. Finally, we can explicitly define a probability distribution $p_{\alpha} = |q_{\alpha}|/G_\mathcal{N}$ where $G_\mathcal{N} = \sum_{\alpha} |q_{\alpha}|$. This gives us an alternative way to write (<ref>) and (<ref>) as

\begin{align}
\mathcal{U} &= G_\mathcal{N} \sum_{\alpha} \mathrm{sgn}(q_{\alpha})\, p_{\alpha}\, \mathcal{K}_{\alpha} \,,\label{eq:appdx_QuasiNoiselessgate2} \\
{C}_{\mathcal{U}(\rho)} &= G_\mathcal{N} \sum_{\alpha} \mathrm{sgn}(q_{\alpha})\, p_{\alpha}\, \widetilde{C}_{\mathcal{K}_{\alpha}(\rho)} \,,
\end{align}

where $\mathrm{sgn}(q_{\alpha})$ denotes the sign of $q_{\alpha}$. We call this the quasiprobability representation of the gate $\mathcal{U}$. The idea is that if we have access to the set of CPTP maps $\{ \mathcal{B}_{\alpha}\}_\alpha$ in our noisy native hardware gate set, then we can obtain an estimate of the noiseless expectation value ${C}_{\mathcal{U}(\rho)}$ as follows: (1) With probability $p_\alpha$, prepare the circuit corresponding to $\mathcal{K}_{\alpha}(\rho)$ and obtain $\widetilde{C}_{\mathcal{K}_{\alpha}(\rho)}$. (2) Multiply the measurement result by $\mathrm{sgn}(q_{\alpha})G_\mathcal{N}$. (3) Repeat process many times and sum results.

Correcting many gates. So far we have only considered a circuit with a single gate $\mathcal{U}$. We can generalize (<ref>) to a general circuit $\prod_k^{N_g} \mathcal{U}_k$ with ${N_g}$ gates with the quasiprobability representation

\begin{equation}\label{eq:appdx_QuasiCircuit}
\prod_k \mathcal{U}_k = G^{tot}_\mathcal{N} \sum_{\vec{i}} \mathrm{sgn}(q_{\vec{i}}) p_{\vec{i}} \mathcal{K}_{\vec{i}} \,,
\end{equation}

where $G^{tot}_\mathcal{N} = \prod_k G_k$,  $\vec{i}=(i_1, ..., i_{N_g})$,  $q_{\vec{i}} = \prod_k q_{i_k}$,  $p_{\vec{i}} = \prod_k p_{i_k}$,  $\mathcal{K}_{\vec{i}} = \prod_k \mathcal{K}_{i_k}$. Thus, a similar procedure can be carried out in order to mitigate the noise on each individual gate in the circuit.

§ USEFUL LEMMAS

§.§ Noise Induced Cost Concentration

Consider a parameterized noisy cost function $\widetilde{C}(\boldsymbol{\theta})=\Tr\left[\widetilde{\rho}(\boldsymbol{\theta})O\right]$, where $\widetilde{\rho}(\boldsymbol{\theta})$ is an $n$-qubit noisy state given by (<ref>) and $\thv \in \Theta$ is drawn from some set of accessible parameters $\Theta$.  Suppose the cost is suffering from exponential cost concentration according to Ref. <cit.>, that is

\begin{equation}\label{eq:appdx_costconcentration}
\left| \widetilde{C}(\boldsymbol{\theta}) - \frac{\Tr[O]}{2^n} \right| \leq q^L B(n)\left\|\rho_{in} - \frac{\emph \id}{2^n} \right\| \,,
\end{equation}

for all $\thv \in \Theta$, where $\rho_{in}$ is the input state in (<ref>), $0\leq q<1$ is some noise parameter, and $B(n) \in \poly(n)$.  Then, $\exists A(n) \in \poly(n)$ such that

\begin{equation}
\Big|\widetilde{C}(\boldsymbol{\theta}_1,q) - \widetilde{C}(\boldsymbol{\theta}_2,q)\Big| \leq A(n) q^L\,,
\end{equation}

for any two sets of parameters $\thv_1, \thv_2 \in \Theta$.

Starting from Eq. (<ref>) can simply write

\begin{equation}
\Big|\widetilde{C}(\boldsymbol{\theta}_1,q) - \widetilde{C}(\boldsymbol{\theta}_2,q)\Big|\, \leq\, 4B(n) q^{L} \,,
\end{equation}

for any two sets of parameters $\boldsymbol{\theta}_1, \boldsymbol{\theta}_2$, where we have used the triangle inequality in $1D$ and the fact that the trace distance has a maximum value of 1.

For the next lemma we will consider the $n$-qubit channel

\begin{equation}
\mathcal{W} = \mathcal{U}_k\circ \mathcal{N}\circ \cdots \circ \mathcal{N} \circ \mathcal{U}_2 \circ \mathcal{N} \circ \mathcal{U}_1\circ\mathcal{N}\,,
\end{equation}

where $\{\UC_k\}^L_{k=1}$ denote unitary channels that describe collections of gates that act together in a layer, and $\NC=\bigotimes_{i=1}^n \DC_i$ is an instance of local depolarizing channels, with depolarizing probabilities $\{p_i \}_{i=1}^n$. We denote $q=\max\{(1-p_i)\}_{i=1}^n <1$.

Consider $\mathcal{W}$ acting on some input state $\rho_{in}$. Then we have

\begin{equation}
\left\| \mathcal{W}(\rho_{in}) - \frac{\emph{\id}}{2^n} \right\|_1 \leq q^{k} n^{1/2} \sqrt{2\ln 2} \,.
\end{equation}

We have

\begin{align}
\left\| \mathcal{W}(\rho_{in}) - \frac{{\id}}{2^n} \right\|_1 &\leq \sqrt{2 \ln 2\, D\Big(\mathcal{W}(\rho_{in})\Big\|\frac{\id}{2^n}\Big)} \\
&\leq \sqrt{2 \ln 2\,q^{2k} D\Big(\rho_{in}\Big\|\frac{\id}{2^n}\Big)} \\
&\leq q^{k} n^{1/2}\sqrt{2\ln 2}\,,
\end{align}

where $D(\cdot\|\cdot)$ denotes the relative entropy. The first inequality is Pinsker's <cit.>, the second inequality is due to $k$ repeated applications of the strong data processing inequality of Ref. <cit.> and the unitary invariance of the relative entropy, and the third inequality is due to an upper bound on the relative entropy of purity.

§.§ Averages over unitary 2-designs

Consider the cost function value $C_\sigma(U)=\Tr[U\sigma U^\dag O]$ where $U$ is a $d\times d$ unitary matrix and $\sigma\in S(\HC)$ is some quantum state. Consider expectation values over $U_i\in Y$ where $Y\subset\mathcal{U}(d)$ is a unitary 2-design and $\mathcal{U}(d)$ is the unitary group of degree $d$. Denote such expectation values as $\langle \cdot \rangle_U$.
Then, we have

\begin{align}
&\langle C_\sigma \rangle_U = \frac{1}{d} \Tr[O]\,, \\
&\langle C_\rho C_\sigma \rangle_U = \frac{\Tr[O^2]\left(d\Tr[\rho\sigma] -  1 \right) - \Tr[O]^2\left(\Tr[\rho\sigma]-d \right)}{d(d^2-1)} \,,
\end{align}

for any two operators $\sigma,\rho \in B(\mathcal{H})$, $\mathrm{dim}(\mathcal{H}) = d$. This implies

\begin{equation}
\Var[C_\sigma] = \frac{\left(\Tr[O^2]- \frac{1}{d}\Tr[O]^2\right)  \left(\Tr[\sigma^2]-\frac{1}{d} \right)}{d^2-1} \,.
\end{equation}

We use the following standard expressions for integration with respect to the Haar measure over the unitary group of degree $d$:

\begin{align}
\int_{\mathcal{U}(d)} d \mu(W) w_{i, j} w_{p, k}^{*} &= \frac{\delta_{i, p} \delta_{j, k}}{d}\,, \\
\int_{\mathcal{U}(d)} d \mu(W) w_{i_{1}, j_{1}} w_{i_{2}, j_{2}} w_{i_{1}', j_{1}'}^{*} w_{i_{2}', j_{2}'}^{*} &=\frac{1}{d^{2}-1}\left(\delta_{i_{1}, i_{1}'} \delta_{i_{2}, i_{2}'} \delta_{j_{1}, j_{1}'} \delta_{j_{2}, j_{2}'}+\delta_{i_{1}, i_{2}'} \delta_{i_{2}, i_{1}'} \delta_{j_{1}, j_{2}'} \delta_{j_{2}, j_{1}'}\right) \\
&\quad\, -\frac{1}{d\left(d^{2}-1\right)}\left(\delta_{i_{1}, i_{1}'} \delta_{i_{2}, i_{2}'} \delta_{j_{1}, j_{2}'} \delta_{j_{2}, j_{1}'}+\delta_{i_{1}, i_{2}'} \delta_{i_{2}, i_{1}'} \delta_{j_{1}, j_{1}'} \delta_{j_{2}, j_{2}'}\right)\,, \nonumber
\end{align}

where $w_{i,j}$ are the matrix elements of the unitary $W \in \mathcal{U}(d)$. Then, the expectation values over 2-designs can be evaluated as

\begin{equation}
\langle C_\sigma \rangle_U = \frac{1}{d}  \Tr[O]\,,
\end{equation}

\begin{align}
\langle C_\rho C_\sigma \rangle_U =& \frac{1}{d^2-1}\Tr[O]^2 + \frac{1}{d^2-1}\Tr[\rho\sigma]\Tr[O^2]\\ &\ - \frac{1}{d(d^2-1)}\Tr[\rho\sigma]\Tr[O]^2 - \frac{1}{d(d^2-1)}\Tr[O^2]\,, \nonumber
\end{align}

where we have used $\Tr[\rho]=\Tr[\sigma]=1$. The final statement comes by noting that

\begin{align}
\Var[C_\sigma] &= \langle C_\sigma^2 \rangle_U - \langle C_\sigma \rangle_U^2 \\
&= \frac{\Tr[O^2]\left(d\Tr[\sigma^2] -  1 \right) - \Tr[O]^2\left(\Tr[\sigma^2]-d \right)}{d(d^2-1)} - \frac{1}{d^2} \Tr[O]^2 \\
&= \frac{\Tr[O^2]\left(d\Tr[\sigma^2] - 1 \right)}{d(d^2-1)} - \frac{\Tr[O]^2\left(d\Tr[\sigma^2]-1 \right)}{d^2(d^2-1)}\,,
\end{align}

which can be factorized to give the desired result.

§ EXPONENTIAL ESTIMATOR CONCENTRATION

We present a proof of Theorem <ref>, and restate the result here for convenience.

Consider an error mitigation protocol prepares the quantity

\begin{align}\label{eq:appdx_thm1_estimator}
E_{\sigma,X,M,k} = \Tr\left[X \left(\sigma^{\otimes M}\otimes \ket{0}\bra{0}^{\otimes k}\right)\right]\,,
\end{align}

for some quantum state $\sigma \in S(\mathcal{H})$, for $\ket{0}\bra{0} \in S(\mathcal{H}')$ and for some $X\in B(\mathcal{H}^{\otimes M}\otimes \mathcal{H}'^{\otimes k})$. We suppose $\sigma$ is prepared with a depth $L_{\sigma}$ circuit and experiences noise according to Eq. (<ref>). Under these conditions, $E_{\sigma,X,M,k}$ exponentially concentrates on a state-independent fixed point in the depth of the circuit as

\begin{align}
\left\vert E_{\sigma,X,M,k} - \Tr\left[X \left( \frac{\emph{\id}^{\otimes M}}{2^{Mn}} \otimes \ket{0}\bra{0}^{\otimes k}\right)\right]\right\vert\, &\leq\, G_{\sigma,X,M}(n) \,, \label{eq:prop1_conc_copy}
\end{align}

where $\emph{\id} \in S(\mathcal{H})$ is the $n$-qubit identity operator and

\begin{equation}\label{eq:appdx_thm1_G}
G_{\sigma,X,M}(n) =  \sqrt{\ln 4}\,\big\|X\big\|_\infty M n^{1/2} q^{L_{\sigma}+1}\,,
\end{equation}

with noise parameter $q\in[0,1)$.

\begin{align}
\left| \Tr\Big[ \left(\sigma^{\otimes M}\otimes \ket{0}\bra{0}^{\otimes k}\right) X\Big]- \Tr \left[\left(\frac{\id}{2^n}\otimes\sigma^{\otimes {M-1}} \otimes \ket{0}\bra{0}^{\otimes k}\right)\, X \right]\right|  &= \left|\Tr \left[\left(\big(\sigma-\frac{\id}{2^n}\big)\otimes\sigma^{\otimes M-1}\otimes \ket{0}\bra{0}^{\otimes k}\right) X \right]\right|  \\
&\leq \left\| \big(\sigma-\frac{\id}{2^n}\big)\otimes\sigma^{\otimes M-1} \otimes \ket{0}\bra{0}^{\otimes k} \right\|_1 \big\|X\big\|_\infty \\
&=  \Big\| \sigma-\frac{\id}{2^n} \Big\|_1 \Tr[\sigma]^{M-1} \Tr[\ket{0}\bra{0}]^k \big\|X\big\|_\infty\\ \label{eq:ancillagoeshere}
&\leq q^{L_\sigma+1} n^{1/2} \sqrt{2\ln 2}\, \big\|X\big\|_\infty\,.
\end{align}

The first inequality is due to the matrix Hölder's inequality, and the second inequality follows from Lemma <ref> Similarly, we have $M-1$ further such equations, which we display with the original equation:

\begin{align}
\Tr \Big[\sigma^{\otimes M} X \Big] - \Tr \left[\left(\frac{\id}{2^n}\otimes\sigma^{\otimes {M-1}}\right)\, X \right]
&\leq q^{L_{\sigma}+1} n^{1/2} \sqrt{2\ln 2}\, \Big\|X\Big\|_\infty\,, \\
\Tr \left[\left(\frac{\id}{2^n}\otimes\sigma^{\otimes {M-1}}\right)\, X \right] - \Tr \left[\left(\frac{\id}{2^n}\otimes \frac{\id}{2^n}\otimes \sigma^{\otimes{M-2}}\right)\, X \right]
&\leq q^{L_{\sigma}+1} n^{1/2} \sqrt{2\ln 2}\, \Big\|X\Big\|_\infty\,,\\
\Tr \left[\left(\frac{\id}{2^n}\otimes \frac{\id}{2^n}\otimes \sigma^{\otimes{M-2}}\right)\, X \right] - \Tr \left[\left(\frac{\id}{2^n}\otimes \frac{\id}{2^n}\otimes\frac{\id}{2^n}\otimes \sigma^{\otimes{M-3}}\right)\, X \right]
&\leq q^{L_{\sigma}+1} n^{1/2} \sqrt{2\ln 2}\, \Big\|X\Big\|_\infty\,,\\
\Tr \left[\left(\Big(\frac{\id}{2^n}\Big)^{\otimes M-1}\otimes\sigma\right)\, X \right] - \Tr \left[\Big(\frac{\id}{2^n}\Big)^{\otimes M}\, X \right]
&\leq q^{L_{\sigma}+1} n^{1/2} \sqrt{2\ln 2}\, \Big\|X\Big\|_\infty\,.
\end{align}

The summation of these equations gives

\begin{equation}
\Tr\Big[ X \left(\sigma^{\otimes M}\otimes \ket{0}\bra{0}^{\otimes k}\right) \Big] - \Tr\left[X \left( \frac{{\id}^{\otimes M}}{2^{Mn}} \otimes \ket{0}\bra{0}^{\otimes k}\right)\right]
\leq M q^{L_{\sigma}+1} n^{1/2} \sqrt{2\ln 2}\, \Big\|X\Big\|_\infty\,,
\end{equation}

which gives the desired bound.

We now present a more detailed statement of Corollary <ref>, which explains how one can spend an exponential number of resources in different ways in order to resolve concentrated cost values.

Consider an error mitigation protocol that approximates the noise-free cost value $C(\thv)$ by estimating the quantity

\begin{equation}\label{eq:appdx_cor1_estimator}
C_m(\thv) = \sum_{(\sigma(\boldsymbol{\theta}),X,M,k)\in T} a_{X,M,k}E_{\sigma(\boldsymbol{\theta}),X,M,k}\,,
\end{equation}

where each $E_{\sigma,X,M,k}$ takes the form (<ref>). We denote $M_{max}$ and $a_{max}$ as the maximum values of $M$ and $a_{X,M,k}$ respectively accessible from a set $T$ defined by the given protocol. Assuming $\|X\|_\infty \in \OC(\poly(n))$, there exists a fixed point $F$ independent of $\boldsymbol{\theta}$ such that

\begin{equation}\label{eq:appdx_cor1_estimator}
\big|C_m(\boldsymbol{\theta})-F\big|\in\OC(2^{-\beta n}a_{max}|T|M_{max})\,,
\end{equation}

for some constant $\beta\geq 1$ if the circuit depths satisfy

\begin{equation}
L_{\sigma(\boldsymbol{\theta})} \in \Omega(n)\,,
\end{equation}

for all $\sigma(\boldsymbol{\theta})$ in the construction (<ref>). That is, if the depth of the circuits scale linearly or greater then one requires at least exponential resources to distinguish $C_m$ from its fixed point, for instance in one of the following ways:

* $a_{max}|T|M_{max} \in \OC(\poly(n))$ and an exponentially large number of shots are used to distinguish two quantities with exponentially small separation

* $a_{max}|T| \in \Omega(2^{\beta' n})$ for some constant $\beta' \geq 1$ and an exponentially large number of shots are required to distinguish two quantities with exponentially large statistical uncertainty, as measurement outcomes are multiplied by $a_{max}|T|$.

* $M_{max} \in \Omega(2^{\beta'' n})$ for some constant $\beta'' \geq 1$ and an exponentially large number of copies of some quantum state $\sigma$ are required.

Explicitly applying the results of Theorem (<ref>) to the construction of $C_m(\boldsymbol{\theta})$ in (<ref>) we have

\begin{align}
\left|C_m(\boldsymbol{\theta})-\sum_{(\sigma(\boldsymbol{\theta}),X,M,k)\in T} a_{X,M,k}\,\Tr\left[X \left( \frac{\id^{\otimes M}}{2^{Mn}} \otimes \ket{0}\bra{0}^{\otimes k}\right)\right] \right| &\leq \sum_{(\sigma(\boldsymbol{\theta}),X,M,k)\in T} a_{X,M,k}\, G_{\sigma(\thv),X,M}(n)\\
& \in \OC\Bigg(\sum_{(\sigma(\boldsymbol{\theta}),X,M,k)\in T} a_{X,M,k} \big\|X\big\|_\infty M n^{1/2} q^{L_{\sigma(\thv)}+1} \Bigg)
\end{align}

where in the second line we have used (<ref>). If $L_{\sigma(\boldsymbol{\theta})} \in \Omega(n)$ then $q^{L_{\sigma(\thv)}+1}\in\OC(2^{-\beta(\thv)n})$ for some $\beta(\thv)\geq1$. Thus, we can write

\begin{equation}
\left|C_m(\boldsymbol{\theta})-\sum_{(\sigma(\boldsymbol{\theta}),X,M,k)\in T} a_{X,M,k}\,\Tr\left[X \left( \frac{\id^{\otimes M}}{2^{Mn}} \otimes \ket{0}\bra{0}^{\otimes k}\right)\right] \right| \in \OC(2^{-\beta n}a_{max}|T|M_{max})
\end{equation}

as required, where we can denote $\beta = \min_{\thv}\beta(\thv)$ and the fixed point as $F$, noting that $F$ is indeed parameter independent. From here, we can inspect the three presented cases:

* If $a_{max}|T|M_{max} \in \OC(\poly(n))$ then $C_m$ has exponentially small separation from $F$.

* There exists choice $\beta' \geq 1$ such that if $a_{max}|T| \in \Omega(2^{\beta' n})$ such that $C_m$ is not exponentially concentrated on $F$, however, $C_m$ now has an exponentially large statistical uncertainty, as measurement outcomes are multiplied by coefficients of order $a_{max}|T|$.

* There exists choice of $\beta'' \geq 1$ such that $M_{max} \in \Omega(2^{\beta'' n})$ and $C_m$ is not exponentially concentrated on $F$.

§ PROTOCOL-SPECIFIC RESULTS

§.§ Zero Noise Extrapolation

In this section we present our results for Zero Noise Extrapolation. As discussed in Section <ref> of the main text, we will consider a Richardson extrapolation strategy based on Eq. (<ref>), an exponential extrapolation strategy based on Eq. (<ref>) and a NIBP extrapolation strategy based on Eq. (<ref>). As we deal with two types of noise parameters, throughout this section we will adopt the unifying notation

\begin{equation}\label{eq:appdx_ZNE_notation}
\widetilde{C}(\thv,a) = \begin{cases}
\widetilde{C}(\thv,a\varepsilon) &\text{for Richardson/exponential extrapolation} \\
\widetilde{C}(\thv,q/a) \quad &\text{for NIBP extrapolation}\,,
\end{cases}
\end{equation}

for all $a\geq 1$. Thus, $\widetilde{C}(\thv,a)$ denotes the noisy cost value at the boosted noise level, and $\widetilde{C}(\thv,1)$ denotes the noisy cost value at the base noise level.

Throughout this section, in order to estimate the sample cost of error mitigation we will make the key assumption that

\begin{equation}\label{eq:appdx_ZNE_assumption}
\Var[\widetilde{C}(\thv,a)] \geq \Var[\widetilde{C}(\thv,1)]\,,
\end{equation}

for all $\thv$ and $a\geq 1$ that is, the statistical fluctuations in measurement outcomes at the boosted noise level are no smaller than that at the base noise level. We note that similar assumptions are made in the literature for Zero Noise Extrapolation and Quasi-Probability Methods <cit.>. Indeed, for noise models with a maximally mixed fixed point, we expect that high noise rates will tend to lead to larger variances. For example, in the simple scenario of a local Pauli measurement, the variance of measurement outcomes takes the form $\frac{(1-p_0)(p_0)}{N}$, where $p_0$ is the probability of obtaining a $"0"$ outcome and $N$ is the number of shots. This variance is maximized for $p_0=\frac{1}{2}$.

§.§.§ Relative resolvability under global depolarizing noise

Consider a circuit with $L$ instances of global depolarizing noise of the form (<ref>).
Consider a Richardson extrapolation strategy based on Eq. (<ref>), an exponential extrapolation strategy based on Eq. (<ref>) and a NIBP extrapolation strategy based on Eq. (<ref>). We presume access to an augmented noisy circuit where the error probability is exactly increased by factor $a_1>1$ as $p\rightarrow a_1p$. Then we have

\begin{align}
&{\rchi_{depol}}\leq \frac{\left(c-\frac{(1-a_1p)^L}{(1-p)^L}\right)^2}{c^2+1}\,,
\end{align}

where $\rchi_{depol}$ is the relative resolvability (see Definition <ref>) for global depolarizing noise, and where

\begin{equation}
c = \begin{cases}
\quad a_1\quad &\text{for Richardson extrapolation,} \\
\frac{a_1r(\varepsilon)^{t(\varepsilon)}}{r(a_1\varepsilon)^{t(a_1\varepsilon)}} \quad &\text{for exponential extrapolation,} \\
\; a_1^{-(L+1)} \quad &\text{for NIBP extrapolation}\,.
\end{cases}
\end{equation}

Thus, $\rchi_{depol}\leq 1$ for all of the above extrapolation strategies with access to 2 noise levels.

Upon inspecting Eqs. (<ref>), (<ref>) and (<ref>), one can verify that the Richardson, exponential and NIBP extrapolation strategies all take the form

\begin{equation}\label{eq:ZNE_generalCm}
C_m(\thv) = \frac{A\cdot\widetilde{C}(\thv,1)-B\cdot\widetilde{C}(\thv,a)}{D} + E
\end{equation}

where $A,B\geq 0$ (note that for NIBP extrapolation $E$ contains the state-independent cost value that represents the fixed point of the noise) and where we have adopted the notation defined in (<ref>). We note that under $L$ instances of global depolarizing noise (of the form (<ref>)) with error probability $p$, noisy cost differences are given by

\begin{equation}
\Delta \widetilde{C}(a) = {(1-ap)^L} \Delta C\,,
\end{equation}

for any pair of cost function points, where $\Delta C$ is the corresponding noise-free cost difference.

The error-mitigated cost function difference $\Delta C_m = {C}_m(\boldsymbol{\theta}_1) - {C}_m(\boldsymbol{\theta}_2)$ between two arbitrary points is given by

\begin{align}
\Delta{C}_m &= \frac{A\cdot\Delta\widetilde{C}(1)-B\cdot\Delta\widetilde{C}(a)}{D}\label{eq:appdx_ZNE_dCm} \\
&= \frac{A\cdot(1-p)^L\Delta{C}-B\cdot(1-ap)^L\Delta{C}}{D}\,.
\end{align}

Inpsecting (<ref>), we see that the error mitigation cost can be bounded simply as

\begin{align}
\gamma &= \frac{A^2+B^2\frac{\Var[\widetilde{C}(\thv,a)]}{\Var[\widetilde{C}(\thv,1)]}}{D^2}\\
&\geq \frac{A^2+B^2}{D^2}
\end{align}

where the inequality comes from our core assumption (<ref>).
Inserting $\gamma$, $\Delta{C}_m$ and $\Delta\widetilde{C}(1)$ into Definition <ref>, we have

\begin{align}
\rchi_{depol} = \frac{1}{\gamma}\left(\frac{\Delta{C}_m}{\Delta\widetilde{C}(1)}\right)^2 &\leq \frac{\left(A(1-p)^L -B(1-ap)^L\right)^2}{(A^2+B^2)(1-p)^{2L}}\\
&= \frac{\left(c-\frac{(1-ap)^L}{(1-p)^L}\right)^2}{c^2+1}\,,
\end{align}

where we have denoted $c=A/B$. By inspecting the specific values of $A$ and $B$ for the Richardson, exponential and NIBP extrapolation strategies respectively, we obtain the results for each strategy.

§.§.§ Average relative resolvability

Consider a Richardson extrapolation strategy based on Eq. (<ref>), an exponential extrapolation strategy based on Eq. (<ref>) and a NIBP extrapolation strategy based on Eq. (<ref>). We presume perfect access to an augmented noisy circuit where the noise rate is increased by factor $a_1>1$. We denote $\boldsymbol{\theta}_{\varepsilon*}$ as the parameter corresponding to the global cost minimum at base noise parameter $\varepsilon$.
Further denote $ \frac{\langle\Delta\widetilde{C}(\thv_{i,\varepsilon*}, a_1\varepsilon)\rangle_{i}}{\langle\Delta\widetilde{C}(\thv_{i,\varepsilon*},\varepsilon)\rangle_{i}} = z$.
Any such noise model has an average relative resolvability

\begin{equation}
\overline{\rchi} \leq  \frac{(z-c)^2}{c^2+1}\,,
\end{equation}

\begin{equation}
c = \begin{cases}
\quad a_1\quad &\text{for Richardson extrapolation,} \\
\frac{a_1r(\varepsilon)^{t(\varepsilon)}}{r(a_1\varepsilon)^{t(a_1\varepsilon)}} \quad &\text{for exponential extrapolation,} \\
\; a_1^{-(L+1)} \quad &\text{for NIBP extrapolation}\,.
\end{cases}
\end{equation}

Thus, under the assumption that $z\leq1$ and $\langle\Delta \widetilde{C}(\thv_{i,\varepsilon*},a_1\varepsilon)\rangle_i\geq0$, $\overline{\rchi}\leq 1$ for all of the above extrapolation strategies with access to 2 noise levels.

As in the previous proof, we can inspect Eqs. (<ref>), (<ref>) and (<ref>), and see that the Richardson, exponential and NIBP extrapolation strategies all take the form

\begin{equation}\label{eq:ZNE_generalCm}
C_m(\thv) = \frac{A\cdot\widetilde{C}(\thv,1)-B\cdot\widetilde{C}(\thv,a)}{D} + E
\end{equation}

where $A,B, D\geq 0$ (note that $E$ contains the state-independent cost value that represents the fixed point of the noise) and we have adopted the notation of (<ref>). The average mitigated cost differences (averaged over accessible parameters $\{\thv_i\}_i$) can be written

\begin{equation}
\langle\Delta C_m(\thv_{i,{\varepsilon*}})\rangle_i = \frac{A\cdot\langle\Delta \widetilde{C}(\thv_{i,{\varepsilon*}},1)\rangle_i- B\cdot\langle\Delta \widetilde{C}(\thv_{i,{\varepsilon*}},a)\rangle_i}{D} \,.
\end{equation}

Thus, we have

\begin{align}
\frac{\langle\Delta C_m(\thv_{i,{\varepsilon*}})\rangle_i}{\langle\Delta \widetilde{C}(\thv_{i,{\varepsilon*}})\rangle_i}  &=   \frac{A  - B\frac{\langle\Delta \widetilde{C}(\thv_{i,{\varepsilon*}},a)\rangle_i}{\langle\Delta \widetilde{C}(\thv_{i,{\varepsilon*}},1)\rangle_i}  }{D} \\
&=\frac{A  - Bz }{D} \,.
\end{align}

finally, by noting once again that the error mitigation cost is simply bounded as $\gamma \geq \frac{A^2+B^2}{D^2}$ due to (<ref>), we have

\begin{equation}
\rchi \leq \frac{(A-Bz)^2}{A^2 + B^2} \,,
\end{equation}

where we can obtain the desired form by defining $c=A/B$. Finally, the specific values of $c$ for each extrapolation strategy can be read off by inspecting Eqs. (<ref>), (<ref>) and (<ref>).

We now introduce a modification of the 2-design-averaged relative resolvability in Definition <ref> that we will use to prove an additional result for Zero Noise Extrapolation.

Consider a spectrum $\boldsymbol{\lambda}\in \mathbb{R}_{\geq0}^{2^n}$ with unit $\ell_1$-norm, which corresponds to a noisy reference state. Then define the unitarily-averaged relative resolvability as

\begin{equation}
\widehat{\overline{\rchi}}_{\boldsymbol{\lambda}} = \frac{1}{\gamma} \frac{\big\langle\big(\widehat{C}_m(\rho,U_i,O_{\boldsymbol{\lambda}}) - \Tr[O_{\boldsymbol{\lambda}}]/2^n \big)^2\big\rangle_{U_i}}{\big\langle\big(\widetilde{C}(\rho,U_i,O_{\boldsymbol{\lambda}}) - \Tr[O_{\boldsymbol{\lambda}}]/2^n \big)^2\big\rangle_{U_i} }\,,
\end{equation}

where $\langle\cdot\rangle_{U_i}$ denotes an average over $U_i$ drawn from a unitary 2-design, and where we denote

\begin{align}
\widetilde{C}(\rho_{\boldsymbol{\lambda}},U_i,O_{\boldsymbol{\lambda}}) &= \Tr[U_i\rho_{\boldsymbol{\lambda}}U^\dag_i O_{\boldsymbol{\lambda}}] \\
\widehat{C}_m(\rho_{\boldsymbol{\lambda}},U_i,O_{\boldsymbol{\lambda}}) &= \Tr[U_i\mathcal{M}(\rho_{\boldsymbol{\lambda}})U^\dag_i O_{\boldsymbol{\lambda}}]
\end{align}

where $\mathcal{M}: S(\mathcal{H})\mapsto B(\mathcal{H})$ is the map that describes the action of the error mitigation protocol.

Note that for Zero Noise Extrapolation we would generally expect to apply the map that describes the error mitigation after conjugation by the unitary $U_i$, counter to the definition presented above. Definition <ref> presumes that augmenting the noise level and conjugating the noisy state by a unitary are two processes that approximately commute, and as such its use is not justified for generic cases. However, below we still present a result using this definition as a mathematical curiosity.

[Average relative resolvability III with Zero Noise Extrapolation]
Consider a Richardson extrapolation strategy based on Eq. (<ref>), an exponential extrapolation strategy based on Eq. (<ref>) and a NIBP extrapolation strategy based on Eq. (<ref>). We presume perfect access to an augmented noisy circuit where the noise rate is increased by factor $a>1$. Denote the output state at the base and augmented noise levels as $\rho(1)$ and $\rho(a)$ respectively. Then we have

\begin{align}
\widehat{\overline{\rchi}}_{\boldsymbol{\lambda}} \leq \frac{c^2+\frac{ P(a)-1/2^n}{P(1)-1/2^n}}{c^2+1}\,,
\end{align}

where $\widehat{\overline{\rchi}}_{\boldsymbol{\lambda}}$ is the averaged relative resolvability defined in Definition <ref>, $P(1)$ is the purity of $\rho(1)$, $P(a)$ is the purity of $\rho(a)$, and

\begin{equation}
c = \begin{cases}
\quad a\quad &\text{for Richardson extrapolation} \\
\frac{ar(\varepsilon)^{t(\varepsilon)}}{r(a\varepsilon)^{t(a\varepsilon)}} \quad &\text{for exponential extrapolation} \\
\; a^{-(L+1)} \quad &\text{for NIBP extrapolation}\,.
\end{cases}
\end{equation}

Thus, $\widehat{\overline{\rchi}}_{\boldsymbol{\lambda}}\leq1$ when $P(a)\leq P(1)$.

We denote reference states $\widetilde{\rho}(\varepsilon)$ and $\widetilde{\rho}(a\varepsilon)$ as states with purity $P(\varepsilon)$ and $P(a\varepsilon)$ respectively. Moreover, denote the noisy cost function values $\widetilde{C}(U_i,\varepsilon) = \Tr[U_i\widetilde{\rho}(\varepsilon)U_i^\dag O]$ and $\widetilde{C}(U_i,a\varepsilon) = \Tr[U_i\widetilde{\rho}(a\varepsilon)U_i^\dag O]$ and further denote $C_m(U_i)$ as the corresponding error mitigated estimator. We start again by noting that in all three Zero Noise Extrapolation strategies the estimator takes the form

\begin{equation}\label{eq:ZNE_generalCmIII}
C_m(U_i) = \frac{A\cdot\widetilde{C}(U_i,1)-B\cdot\widetilde{C}(U_i,a)}{D} + E
\end{equation}

where $A,B\geq 0$ (see Eqs. (<ref>), (<ref>) and (<ref>)) and we have adopted the notation of (<ref>).
We first evaluate the relevant expectation values which correspond to integrals over the Haar distribution over the unitary group of degree $2^n$. We now proceed to derive the result for Richardson/exponential extrapolation, however, we note that the proof follows in a similar way for NIBP extrapolation with the simple substitution $a\varepsilon \mapsto q/a$. Utilizing Lemma <ref>, we have

\begin{align}
&\langle \widetilde{C}(U_i,\varepsilon) \rangle_{U_i} = \frac{1}{2^n}\Tr[\widetilde{\rho}(\varepsilon)]\Tr[O] = \frac{1}{2^n}\Tr[O]\,, \\[.2em]
&\big\langle (\Delta\widetilde{C}(U_i,\varepsilon))^2\big\rangle_{U_i} = \big\langle\big(\widetilde{C}(U_i,\varepsilon) - \langle \widetilde{C}(U_j,\varepsilon) \rangle_{U_j} \big)^2\big\rangle_{U_i} = \frac{\left(\Tr[O^2]- \frac{1}{2^n}\Tr[O]^2\right)  \left(\Tr[\widetilde{\rho}^2(\varepsilon)]-\frac{1}{2^n}\Tr[\widetilde{\rho}(\varepsilon)]^2 \right)}{2^{2n}-1} \\
&\hspace{22.13em} = \frac{\Tr[O^2]- \frac{1}{2^n}\Tr[O]^2 }{2^{2n}-1} \left(P(\varepsilon)-\frac{1}{2^n} \right)\,, \label{eq:ZNE_noisyvar2} \\[1em]
&\big\langle (\Delta\widetilde{C}(U_i,\varepsilon))(\Delta\widetilde{C}(U_i,a\varepsilon))\big\rangle_{U_i} = \frac{\left(\Tr[O^2]- \frac{1}{2^n}\Tr[O]^2\right)  \left(\Tr[\widetilde{\rho}(\varepsilon)\widetilde{\rho}(a\varepsilon)]-\frac{1}{2^n}\Tr[\widetilde{\rho}(\varepsilon)]\Tr[\widetilde{\rho}(a\varepsilon)] \right)}{2^{2n}-1}\\
&\hspace{12.6em}\geq 0 \,, \label{eq:ZNE_noisycov2}
\end{align}

where the inequality comes by observing that $\Tr[\widetilde{\rho}(\varepsilon)]=\Tr[\widetilde{\rho}(a\varepsilon)]=1$ and further applying Cauchy-Schwarz to $\Tr[\widetilde{\rho}(\varepsilon)\widetilde{\rho}(a\varepsilon)]$ and noting that the purity of an $n$ qubit state is lower bounded by $1/2^n$. Inspecting Eq. (<ref>)  we have

\begin{align}
&\langle C_{m}(U_i) \rangle_{U_i} = \frac{1}{2^n}\frac{A-B}{D}\Tr[O] \,,\\[0.2em]
&\big\langle\big(C_{m}(U_i) - \langle C_{m}(U_j) \rangle_{U_j} \big)^2\big\rangle_{U_i} = \left\langle \left( \frac{A\cdot\widetilde{C}(U_i,\varepsilon)-B\cdot\widetilde{C}(U_i,a\varepsilon)}{D}+E - \left( \frac{1}{2^n}\frac{A-B}{D}\Tr[O]+E\right) \right)^2 \right\rangle_{U_i} \\
&\hspace{8.5em} =\left\langle \left( \frac{ A \cdot\Delta\widetilde{C}(U_i,\varepsilon) -B\cdot\Delta\widetilde{C}(U_i,a\varepsilon)}{D} \right)^2 \right\rangle_{U_i} \\
&\hspace{8.5em} = \frac{A^2\big\langle \big(\Delta\widetilde{C}(U_i,\varepsilon)\big)^2\big\rangle_{U_i} + B^2\big\langle \big(\Delta\widetilde{C}(U_i,a\varepsilon)\big)^2\big\rangle_{U_i} - 2AB \big\langle \Delta\widetilde{C}(U_i,\varepsilon)\Delta\widetilde{C}(U_i,a\varepsilon)\big\rangle_{U_i}}{D^2} \\
&\hspace{8.5em} \leq \frac{\Tr[O^2]- \frac{1}{2^n}\Tr[O]^2 }{D^2(2^{2n}-1)}\left( A^2\Big(P(\varepsilon)-\frac{1}{2^n}\Big) + B^2\Big(P(a\varepsilon)-\frac{1}{2^n} \Big) \right)\,. \label{ZNE_estvar2}
\end{align}

The inequality comes by substituting in the expressions for $\big\langle \big(\Delta\widetilde{C}(U_i,\varepsilon)\big)^2\big\rangle_{U_i}$ and $\big\langle \big(\Delta\widetilde{C}(U_i,a\varepsilon)\big)^2\big\rangle_{U_i}$ obtained in (<ref>), and dropping the third term in the numerator, where we have used Eq. (<ref>).
Finally, we note that (<ref>) gives $\gamma^{-1} = \frac{D^2}{A^2 + B^2}$. Substituting the obtained expressions for $\gamma^{-1}$, (<ref>) and (<ref>) into Definition <ref> we obtain

\begin{equation}
\dbloverline{\chi}_{\boldsymbol{\lambda}} \leq \frac{A^2 + B^2\frac{P(a\varepsilon)-1/2^n}{P(\varepsilon)-1/2^n}}{A^2 + B^2}\,,
\end{equation}

where we can define $c=A/B$ to obtain the desired result. Further, the explicit form of $c$ for Richardson, exponential and NIBP extrapolation can be respectively found by inspecting Eqs. (<ref>), (<ref>) and (<ref>).

§.§.§ Richardson extrapolation with 3 noise levels

In this section we focus on Richardson extrapolation (see Appendix <ref> for review) and investigate the change in resolvability under an extrapolation strategy that utilizes 3 distinct noise levels.

[Relative resolvability of Richardson extrapolation with global depolarizing noise, 3 noise levels]
Consider $L$ instances of global depolarizing noise of the form (<ref>) acting through a circuit.
Consider a Richardson extrapolation strategy based on Eq. (<ref>), an exponential extrapolation strategy based on Eq. (<ref>) and a NIBP extrapolation strategy based on Eq. (<ref>) in the appendix. We presume access to two augmented noisy circuits where the error probability is perfectly increased by factors $a_2>a_1>1$ as $p\rightarrow a_1p$ and $p\rightarrow a_2p$ respectively. Then for all three extrapolation strategies and any such choices of $a_2$ and $a_1$, we have

\begin{align}
&{\rchi_{depol}}\leq 1\,,
\end{align}

where $\rchi_{depol}$ is the relative resolvability (see Definition <ref>) for global depolarizing noise.

We start by noting that under $L$ instances of global depolarizing noise with error probability $p$ (of the form (<ref>)), noisy cost differences are given by

\begin{equation}\label{eq:appdx_gdepol2}
\Delta \widetilde{C}(a) = {(1-ap)^L}\Delta C\,,
\end{equation}

for any noise augmentation factor $a$ and any pair of cost function points, where $\Delta C$ is the corresponding noise-free cost difference.

The error-mitigated cost function difference $\Delta C_m (\thv_{1,2}) = {C}_m(\boldsymbol{\theta}_1) - {C}_m(\boldsymbol{\theta}_2)$ between two arbitrary points constructed under Richardson extrapolation with 3 noise levels is given by

\begin{align}
\Delta{C}_m &= \frac{a_1a_2(a_2-a_1) \Delta\widetilde{C}(p) - a_2(a_2-1)\Delta\widetilde{C}( a_1p) + a_1(a_1-1)\Delta\widetilde{C}( a_2p) }{(a_1-1)(a_2-1)(a_2-a_1)} \\
&= \frac{a_1a_2(a_2-a_1) {(1-p)^L}\Delta C - a_2(a_2-1){(1-a_1p)^L}\Delta C + a_1(a_1-1){(1-a_2p)^L}\Delta C }{(a_1-1)(a_2-1)(a_2-a_1)}\,,
\end{align}

where in order to obtain the first equality we have used (<ref>). The second equality comes by substituting in (<ref>). Inspecting (<ref>), we see that the error mitigation cost can be bounded simply as

\begin{align}
\gamma &= \frac{a_1^2a_2^2(a_2-a_1)^2+a^2_2(a_2-1)^2 \frac{\Var[\widetilde{C}(\thv,a_1p)]}{\Var[\widetilde{C}(\thv,p)]} + a^2_1(a_1-1)^2\frac{\Var[\widetilde{C}(\thv,a_2p)]}{\Var[\widetilde{C}(\thv,p)]}}{(a_1-1)^2(a_2-1)^2(a_2-a_1)^2}\\
&\geq \frac{a_1^2a_2^2(a_2-a_1)^2+a^2_2(a_2-1)^2 + a^2_1(a_1-1)^2 }{(a_1-1)^2(a_2-1)^2(a_2-a_1)^2}
\end{align}

for any $\thv$, where the inequality comes from our core assumption (<ref>).
Inserting our expressions for $\gamma$, $\Delta{C}_m$ and $\Delta\widetilde{C}(1)$ into Definition <ref>, we have

\begin{align}
\rchi_{depol} = \frac{1}{\gamma}\left(\frac{\Delta{C}_m}{\Delta\widetilde{C}(1)}\right)^2 &\leq \frac{\left(a_1a_2(a_2-a_1) {(1-p)^L} - a_2(a_2-1){(1-a_2p)^L} + a_1(a_1-1){(1-a_2p)^L}\right)^2}{\left(a_1^2a_2^2(a_2-a_1)^2+a^2_2(a_2-1)^2 + a^2_1(a_1-1)^2\right)(1-p)^{2L}}\\
&= \frac{\left(a_1a_2(a_2-a_1)  - a_2(a_2-1)\frac{(1-a_1p)^L}{(1-p)^L} + a_1(a_1-1)\frac{(1-a_2p)^L}{(1-p)^L}\right)^2}{a_1^2a_2^2(a_2-a_1)^2+a^2_2(a_2-1)^2 + a^2_1(a_1-1)^2}\,.
\end{align}

The desired result can be observed by noting that $a_2(a_2-1)>a_1(a_1-1)$ and that $\frac{(1-a_1p)^L}{(1-p)^L}>\frac{(1-a_2p)^L}{(1-p)^L}$.

[Average resolvability of Richardson extrapolation, 3 noise levels]
Consider a Richardson extrapolation strategy based on Eq. (<ref>), an exponential extrapolation strategy based on Eq. (<ref>) and a NIBP extrapolation strategy based on Eq. (<ref>) in the appendix. We presume perfect access to two augmented noisy circuits where the noise rate is increased by factors $a_2>a_1>1$. We denote $\boldsymbol{\theta}_{\varepsilon*}$ as the parameter corresponding to the global cost minimum at base noise parameter $\varepsilon$. Further denote $ \frac{\langle\Delta\widetilde{C}(\thv_{i,{\varepsilon*}},a_1\varepsilon)\rangle_{i}}{\langle\Delta\widetilde{C}(\thv_{i,{\varepsilon*}},\varepsilon)\rangle_{i}} = z_1$ and $\frac{\langle\Delta\widetilde{C}(\thv_{i,{\varepsilon*}},a_2\varepsilon)\rangle_{i}}{\langle\Delta\widetilde{C}(\thv_{i,{\varepsilon*}},\varepsilon)\rangle_{i}} = z_2$. Any such noise model has an average relative resolvability

\begin{equation}
\overline{\rchi} \leq  \frac{\left(a_1a_2(a_2-a_1)  - a_2(a_2-1)z_1 + a_1(a_1-1)z_2\right)^2}{a_1^2a_2^2(a_2-a_1)^2+a^2_2(a_2-1)^2 + a^2_1(a_1-1)^2}\,,
\end{equation}

where $\overline{\rchi}$ is the averaged relative resolvability (see Definition <ref>).
Thus, under the assumption that $z_2\leq z_1\leq1$ (on average the cost concentrates with increasing noise level) and $\langle\Delta \widetilde{C}_{i,{\varepsilon*}}(a_1\varepsilon)\rangle_i, \langle\Delta \widetilde{C}_{i,{\varepsilon*}}(a_2\varepsilon)\rangle_i \geq0$ (boosting the noise level does not shift the cost value of the global minimum above the average cost value), then $\overline{\rchi}\leq 1$.

The averaged error-mitigated cost function difference $\langle\Delta C_m(\thv_{i,{\varepsilon*}})\rangle_i = \langle{C}_m(\boldsymbol{\theta}_i) - {C}_m(\boldsymbol{\theta}_{\varepsilon*})\rangle_i$ between two arbitrary points constructed under Richardson extrapolation with 3 noise levels is given by

\begin{align}
\langle\Delta C_m(\thv_{i,{\varepsilon*}})\rangle_i &= \left\langle\left(\frac{a_1a_2(a_2-a_1) \Delta\widetilde{C}(\thv_{i,{\varepsilon*}},p) - a_2(a_2-1)\Delta\widetilde{C}(\thv_{i,{\varepsilon*}}, a_1p) + a_1(a_1-1)\Delta\widetilde{C}(\thv_{i,{\varepsilon*}}, a_2p) }{(a_1-1)(a_2-1)(a_2-a_1)}\right)_{i,{\varepsilon*}}\right\rangle_i \\
&= \frac{a_1a_2(a_2-a_1) \langle\Delta\widetilde{C}(\thv_{i,{\varepsilon*}},p)\rangle_i - a_2(a_2-1)\langle\Delta\widetilde{C}((\thv_{i,{\varepsilon*}},a_1p)\rangle_i + a_1(a_1-1)\langle\Delta\widetilde{C}(\thv_{i,{\varepsilon*}},a_2p)\rangle_i }{(a_1-1)(a_2-1)(a_2-a_1)}\,.
\end{align}

As in the previous proof, we can inspect (<ref>) and we see that the error mitigation cost can be bounded simply as

\begin{align}
\gamma &= \frac{a_1^2a_2^2(a_2-a_1)^2+a^2_2(a_2-1)^2 \frac{\Var[\widetilde{C}(\thv,a_1p)]}{\Var[\widetilde{C}(\thv,p)]} + a^2_1(a_1-1)^2\frac{\Var[\widetilde{C}(\thv,a_2p)]}{\Var[\widetilde{C}(\thv,p)]}}{(a_1-1)^2(a_2-1)^2(a_2-a_1)^2}\\
&\geq \frac{a_1^2a_2^2(a_2-a_1)^2+a^2_2(a_2-1)^2 + a^2_1(a_1-1)^2 }{(a_1-1)^2(a_2-1)^2(a_2-a_1)^2}
\end{align}

for any $\thv$, where the inequality comes from our core assumption (<ref>).
Inserting our expressions for $\gamma$ and $\Delta{C}_m$ into Definition <ref>, we have

\begin{align}
\overline{\rchi} = \frac{1}{\gamma}\left(\frac{\Delta{C}_m}{\Delta\widetilde{C}(1)}\right)^2 &\leq \frac{\left(a_1a_2(a_2-a_1)  - a_2(a_2-1)\frac{\langle\Delta\widetilde{C}(\thv_{i,{\varepsilon*}},a_1\varepsilon)\rangle_{i}}{\langle\Delta\widetilde{C}(\thv_{i,{\varepsilon*}},\varepsilon)\rangle_{i}} + a_1(a_1-1)\frac{\langle\Delta\widetilde{C}(\thv_{i,{\varepsilon*}},a_2\varepsilon)\rangle_{i}}{\langle\Delta\widetilde{C}(\thv_{i,{\varepsilon*}},\varepsilon)\rangle_{i}}\right)^2}{a_1^2a_2^2(a_2-a_1)^2+a^2_2(a_2-1)^2 + a^2_1(a_1-1)^2}\,,
\end{align}

and the desired result comes by denoting $ \frac{\langle\Delta\widetilde{C}(\thv_{i,{\varepsilon*}},a_1\varepsilon)\rangle_{i}}{\langle\Delta\widetilde{C}(\thv_{i,{\varepsilon*}},\varepsilon)\rangle_{i}} = z_1$ and $\frac{\langle\Delta\widetilde{C}(\thv_{i,{\varepsilon*}},a_2\varepsilon)\rangle_{i}}{\langle\Delta\widetilde{C}(\thv_{i,{\varepsilon*}},\varepsilon)\rangle_{i}} = z_2$.

As with the results of Proposition <ref> we see that $\overline{\rchi}$ decreases with increasing cost concentration.

§.§ Virtual Distillation

§.§.§ Bounds on error mitigation cost

We recall the two error mitigation protocols of Ref. <cit.>, denoted "A" and "B" respectively, to prepare

\begin{equation}\label{eq:VD_A_copy}
C_m^{(A)}(\boldsymbol{\theta}_i)=\Tr[\widetilde{\rho}_i^M O]/\Tr[\widetilde{\rho}_i^M]\,,
\end{equation}

\begin{equation}\label{eq:VD_B_copy}
C_m^{(B)}(\boldsymbol{\theta}_i) = \Tr[\widetilde{\rho}_i^M O]/\lambda_i^M\,,
\end{equation}

where $\lambda_i$ is the dominant eigenvalue of $\widetilde{\rho}_i \equiv \widetilde{\rho}(\boldsymbol{\theta}_i)$. The protocols considered explicitly construct these quantities as

\begin{align}
\Tr[\widetilde{\rho}_iO] &= 2\mathrm{prob}_{1,i} -1 \,,\\
\Tr[\widetilde{\rho}_i^MO] &= 2\mathrm{prob}_{M,i} -1 \,,\\
\Tr[\widetilde{\rho}_i^M] &= 2\mathrm{prob}'_{M,i} -1 \,,
\end{align}

where $\mathrm{prob}_{1,i},\ \mathrm{prob}_{M,i}$ and $\mathrm{prob}'_{M,i}$ are expectation values of a Pauli-$Z$ measurement on a qubit ancillary subsystem. In order to obtain our results we will hereon make the core assumptions

\begin{align}
\Var[\mathrm{prob}_{1,i}] &= \Var[\mathrm{prob}_{1,j}]\quad\forall i\neq j \,,\\
\Var[\mathrm{prob}_{M,i}] &= \Var[\mathrm{prob}_{M,j}]\quad\forall i\neq j \,,\\
\Var[\mathrm{prob}_{M,i}]&\geq \Var[\mathrm{prob}_{1,i}] \quad \forall i,\, M\geq 2\,,
\end{align}

that is, the statistical uncertainty of the measurement outcomes of the circuit that prepares $\widetilde{\rho}_i$ and $\widetilde{\rho}^M_i$ can be approximated to be state independent, and the statistical uncertainty of the measurement outcomes of the circuit that prepares $\widetilde{\rho}^M_i$ are at best equal to that of $\widetilde{\rho}_i$. In the case of large $M$ we expect $\Var[\mathrm{prob}_{M,i}]$ to be large, as for any (non-pure) $\widetilde{\rho}$, the quantity $\Tr[\widetilde{\rho}^MO]$ is close to zero for large $M$. This corresponds to $\mathrm{prob}_{M,i}=\frac{1}{2}$, which maximizes the variance for a binomial distribution.

Denote the error mitigation cost (see Definition <ref>) corresponding to (<ref>) and (<ref>) as $\gamma^{(A)}$ and $\gamma^{(B)}$ respectively. We have

\begin{align} \label{VD_gammas}
\gamma^{(A)}\geq\frac{1}{(\Tr[\widetilde{\rho}^M])^2}\,,\quad \gamma^{(B)}=\frac{1}{\lambda^{2M}}\,.
\end{align}

For $\gamma^{(A)}$ and $\gamma^{(B)}$ we need to compute the variances of the estimators of $C_m^{(A)}$, $C_m^{(B)}$ respectively and likewise $\widetilde{C}=\Tr[\widetilde{\rho} O]$. We have

\begin{align}
\Var[\widetilde{C}] =  \Var[\Tr[\widetilde{\rho} O]] &= \Var[2\textrm{prob}_1 -1]\,,\\
\end{align}

\begin{align}
\Var[C_m^{(B)}] = \Var\Big[\frac{\Tr[\widetilde{\rho}^M O]}{\lambda^M}\Big] &= \frac{1}{\lambda^{2M}}\Var[2\textrm{prob}_M -1]\,,\\
\end{align}

\begin{align}
\Var[C_m^{(A)}] = \Var\Big[\frac{\Tr[\widetilde{\rho}^M O]}{\Tr[\widetilde{\rho}^M]}\Big] &= 4\Var[\textrm{prob}_{M}]\, \left(\mathbb{E}\left[\frac{1}{2\textrm{prob}'_{M}-1}\right]\right)^2+ 4\Tr[\widetilde{\rho}^MO]^{2}\,\Var\left[\frac{1}{2\textrm{prob}'_{M}-1}\right] \label{eq:VD_varEA} \\
&\quad\ + 4\Var[\textrm{prob}_{M}]\,\Var\left[\frac{1}{2\textrm{prob}'_{M}-1}\right]  \nonumber\\
&\geq 4\Var[\textrm{prob}_{M}]\, \left(\mathbb{E}\left[\frac{1}{2\textrm{prob}'_{M}-1}\right]\right)^2 \\
&\geq 4\Var[\textrm{prob}_{M}]\, \frac{1}{\left(\mathbb{E}\left[2\textrm{prob}'_{M}-1\right]\right)^2} \\
&= 4\Var[\textrm{prob}_{M}]\, \frac{1}{\left(\Tr[\widetilde{\rho}^M]\right)^2}\,.
\end{align}

Equation (<ref>) comes from the standard formula for the variance of the product of two independent random variables. To obtain the first inequality we simply drop the second and third terms, which are positive. The second inequality is an application of Jensen's inequality. Recalling the definition of error mitigation cost (Definition <ref>), the above three equations enable us to write

\begin{align} \label{VD_gammas}
\gamma^{(A)}\geq\frac{1}{(\Tr[\widetilde{\rho}^M])^2}\,,\quad \gamma^{(B)}=\frac{1}{\lambda^{2M}}\,,
\end{align}

where we have used our core assumption that $\Var[\textrm{prob}_1]\leq\Var[\textrm{prob}_M]$.

§.§.§ Relative resolvability for global depolarizing noise

Here we present a proof of Proposition <ref>, in which we upper bound the relative resolvability for Virtual Distillation, for any two cost function points under global depolarizing noise.

Consider $l$ instances of global depolarizing noise $\mathcal{D}$ of the form

\begin{equation}
\rho \xrightarrow{\mathcal{D}} \widetilde{\rho} = q^l \rho + (1-q^l) \frac{\emph{\id}}{2^n}
\end{equation}

acting on some pure state $\rho$ with some noise parameter $q\in [0,1)$. We consider the two error mitigation protocols of Ref. <cit.> (denoted "A" and "B") to respectively prepare (<ref>) and (<ref>). The relative resolvability of any pair of arbitrary cost function points satisfies

\begin{align}
\rchi^{(A)} \leq \rchi^{(B)} = \Gamma(n,M,q^l)
\end{align}

for all $n\geq1$, $M\geq2$, $q^l \in [0,1]$, and where

\begin{equation}
\Gamma(n,M,q^l)\leq 1 \,,
\end{equation}

is a monotonically decreasing function in $M$ (with asymptotically exponential decay) in the quadrant $n\geq 1$, $M\geq 2$. The bound is saturated as $\Gamma(1,2,p)= 1$ for all $p$.

In this proof we consider arbitrary cost function differences, that is, given two arbitrarily chosen points in parameter space $\thv_1$ and $\thv_2$, we consider

\begin{equation}
\Delta C = C(\thv_1) - C(\thv_2)\,,
\end{equation}

and the respective differences for the noisy cost $\widetilde{C}(\thv)$ and the mitigated costs $C^{(A)}_m(\thv)$, $C^{(B)}_m(\thv)$. In order to evaluate $\rchi_A$ and $\rchi_B$ we need to first evaluate the following quantities:

\begin{equation}
\Delta \widetilde{C}\,,\ \Delta C_m^{(A)} \,,\ \Delta C_m^{(B)} \,,\  \gamma^{(A),(B)}
\end{equation}

that is, the noisy cost function difference between two points, the difference between the virtual distillation estimators for the same points for both protocols, and the error mitigation rate for both protocols. The noisy cost function difference under global depolarizing noise is simply related the noiseless difference as

\begin{equation}
\Delta \widetilde{C} = q^{l} \Delta C\,. \label{eq:VD_deltaC}
\end{equation}

To evaluate the other quantities we note that

\begin{align}
\widetilde{\rho} &= \Big[q^l + \frac{1}{2^n}(1-q^l)\Big] \rho + \Big[\frac{1}{2^n}(1-q^l) \Big] (\id - \rho)\,, \\
\widetilde{\rho}^M &= \left[ \Big[q^l + \frac{1}{2^n}(1-q^l)\Big]^M - \Big[\frac{1}{2^{n}}(1-q^l)\Big]^M \right]\rho + \frac{2^n}{2^{nM}}(1-q^l)^M \frac{\id}{2^n}\,, \\
\Tr[\widetilde{\rho}^M] &= \Big[q^l + \frac{1}{2^n}(1-q^l)\Big]^M + \frac{2^n-1}{2^{nM}}(1-q^l)^M \,,\quad \label{eq:stateindependent} \\
\Tr[\widetilde{\rho}^M O] &= \left[ \Big[q^l + \frac{1}{2^n}(1-q^l)\Big]^M - \Big[\frac{1}{2^{n}}(1-q^l)\Big]^M \right]\Tr[\rho O] + \frac{1}{2^{nM}}(1-q^l)^M \Tr[O]\,.
\end{align}

In particular, we highlight that the expression for $\Tr[\widetilde{\rho}^M]$. As the dominant noisy eigenvalue $\lambda$ and $\Tr[\widetilde{\rho}^M]$ are state independent we have

\begin{align}
\Delta C_m^{(A)} &= \frac{1}{\Tr[\widetilde{\rho}^M]} \left[ \Big[q^l + \frac{1}{2^n}(1-q^l)\Big]^M - \Big[\frac{1}{2^{n}}(1-q^l)\Big]^M \right] \Delta C\,, \label{eq:VD_deltaEA} \\
\Delta C_m^{(B)} &= \frac{1}{\lambda^M} \left[ \Big[q^l + \frac{1}{2^n}(1-q^l)\Big]^M - \Big[\frac{1}{2^{n}}(1-q^l)\Big]^M \right] \Delta C \,, \label{eq:VD_deltaEB}
\end{align}

where the choice of $\widetilde{\rho}$ is arbitrary.

Now, using Definition <ref> and combining (<ref>), (<ref>), (<ref>) along with the result of Lemma <ref>, we have

\begin{equation}
\chi^{(A)} \leq \chi^{(B)} = \Gamma(n,M,q^l) \,, \label{eq:VD_SM}
\end{equation}

where we define the function

\begin{equation}\label{eq:appdx_Gamma}
\Gamma(n,M,q^l) = \frac{1}{q^{2l}} \left[ \Big[q^l + \frac{1}{2^n}(1-q^l)\Big]^M - \Big[\frac{1}{2^{n}}(1-q^l)\Big]^M \right]^2\,.
\end{equation}

First, note that for $M=2$

\begin{align}
\Gamma(q^l,n,2) &= \frac{1}{q^{2l}} \left[ \Big[q^l + \frac{1}{2^n}(1-q^l)\Big]^2 - \Big[\frac{1}{2^{n}}(1-q^l)\Big]^2 \right]^2 \\
&= \left( q^l + \frac{2}{2^n} (1-q^l) \right)^2\,.
\end{align}

For $n=1$, $\Gamma(n,2,q^l)=1$. For all $n>1$, $\Gamma(n,2,q^l)<1$ as $(1-q^l)>0$. Thus,

\begin{equation}\label{eq:gamma(M=2)}
\Gamma(n,2,q^l)\leq1\quad \forall\ n\geq1 \,.
\end{equation}

We complete the proof by showing that $\Gamma(n,M,q^l)$ monotonically decreases with $M$ for all $n\geq1$, $M\geq2$. This can be seen by inspecting the partial derivative (making the decomposition $\Gamma = (\Gamma^{1/2})^2$ due to the square in (<ref>))

\begin{equation}
\frac{\partial \Gamma}{\partial M} = 2 \Gamma^{1/2} \frac{\partial \Gamma^{1/2}}{\partial M}\,.
\end{equation}

We investigate when this quantity is negative. As $\Gamma^{1/2}$ is always positive, negativity is determined by the sign of $\frac{\partial \Gamma^{1/2}}{\partial M}$. Denoting $\delta = \frac{1}{2^{n}}(1-q^l)$, we have

\begin{align}
\frac{\partial \Gamma^{1/2}}{\partial M} &= \frac{1}{q^l} \left[ (q^l + \delta)^M \ln (q^l + \delta) - \delta^M \ln \delta \right] \\
&= \frac{1}{q^l} \left[  \delta^M\left(\ln (q^l+\delta) - \ln \delta \right) + \ln(q^l + \delta)\left( (q^l + \delta)^M - \delta^M \right) \right] \\
&\leq \frac{1}{q^l}\left[ q^l\delta^{M-1} + (q^l + \delta -1)\left( (q^l+\delta)^M - \delta^M \right) \right] \\
&=\frac{1}{q^l}\left[ q^l\delta^{M-1} - (2^n-1)\delta\left( (q^l+\delta)^M - \delta^M \right) \right] \\
&\leq \frac{1}{q^l}\left[ q^l\delta^{M-1} - (2^n-1)\delta\left( Mq^l\delta^{M-1} + \frac{1}{2}M(M-1)q^{2l}\delta^{M-2} \right) \right] \\
&= \frac{1}{q^l}\left[ q^l\delta^{M-1} - (2^n-1)\delta\left( \frac{1}{2^n}M(1-q^l)q^l\delta^{M-2} + \frac{1}{2}M(M-1)q^{2l}\delta^{M-2} \right) \right] \\
&=  \delta^{M-1}\left[ 1 - \left(1-\frac{1}{2^n}\right)M  - \frac{1}{2}(2^n-1)M\left(M-1-\frac{2}{2^n}\right)q^{l} \right] \\
&\leq \delta^{M-1}\left[ 1 - \frac{1}{2}M  - \frac{1}{2}M\left(M-2\right)q^{l} \right]\quad \forall\ n\geq1 \,,
\end{align}

where in order to obtain the first inequality we use the inequalities $\ln(q^l+\delta) - \ln \delta \leq q^l/\delta$ and $\ln(q^l+\delta) \leq q^l+\delta-1$. The second inequality comes from observing that the expansion of $\left( (q^l+\delta)^M - \delta^M \right)$ is a sum of positive terms, and considering only two such terms. The above implies that

\begin{equation}
\frac{\partial \Gamma}{\partial M} \leq 0 \quad \forall\ n\geq1, M\geq2\,,
\end{equation}

that is, $\Gamma$ is monotonically decreasing with $M$ in the quadrant $n\geq1$ $M\geq2$. Combined with (<ref>), we have the proof as required.

§.§.§ Average relative resolvability

Here we present a proof of Proposition <ref>, in which we upper bound the 2-design-averaged relative resolvability for Virtual Distillation.

Consider an error mitigation protocol that prepares estimator $C_m(\boldsymbol{\theta}_i) = \Tr[\widetilde{\rho}_i^M O]/\Tr[\widetilde{\rho}_i^M]$ from some noisy parameterized quantum state $\widetilde{\rho}_i \equiv \widetilde{\rho}(\boldsymbol{\theta}_i)$. Consider the average relative resolvability for noisy states of some spectrum $\boldsymbol{\lambda}$ with purity $P_{\boldsymbol{\lambda}}$ as defined in Definition <ref>. We have

\begin{equation}\label{eq:appdx_vd_avbound}
\dbloverline{\rchi}_{\boldsymbol{\lambda}}\leq G(n,M,P)\,\leq\,1\,,
\end{equation}

where $G(n,M,P)$ is a monotonically decreasing function in $M$ (with asymptotically exponential decay) for all $n\geq 1$, $M\geq 2$.