Multi-exponential error extrapolation and combining error mitigation techniques for NISQ applications

Noise in quantum hardware remains the biggest roadblock for the implementation of quantum computers. To fight the noise in the practical application of near-term quantum computers, instead of relying on quantum error correction which requires large qubit overhead, we turn to quantum error mitigation, in which we make use of extra measurements. Error extrapolation is an error mitigation technique that has been successfully implemented experimentally. Numerical simulation and heuristic arguments have indicated that exponential curves are effective for extrapolation in the large circuit limit with an expected circuit error count around unity. In this Article, we extend this to multi-exponential error extrapolation and provide more rigorous proof for its effectiveness under Pauli noise. This is further validated via our numerical simulations, showing orders of magnitude improvements in the estimation accuracy over single-exponential extrapolation. Moreover, we develop methods to combine error extrapolation with two other error mitigation techniques: quasi-probability and symmetry verification, through exploiting features of these individual techniques. As shown in our simulation, our combined method can achieve low estimation bias with a sampling cost multiple times smaller than quasi-probability while without needing to be able to adjust the hardware error rate as required in canonical error extrapolation.


I. INTRODUCTION
While fault-tolerant quantum computers promise huge speed-up over classical computers in areas like chemistry simulations, optimisation and decryption, their implementations remain a long term goal due to the large qubit overhead required for quantum error correction.With the recent rapid advance of quantum computer hardware in terms of both qubit quantity and quality, culminating with the "quantum supremacy" experiment [1], one must wonder is it possible for us to perform classically intractable computations on such Noisy Intermediate-Scale Quantum (NISQ) hardware without quantum error correction [2].Resource estimation has been performed for one of the most promising applications on NISQ hardware: the Fermi-Hubbard model simulation [3,4], realising that even with an optimistic local gate error rate of 10 −4 , the large number of gates needed for a classically intractable calculation will lead to an expected circuit error count of the order of unity.To obtain any meaningful results under such an expected circuit error count, it is essential to employ error mitigation techniques, which relies on making extra measurements, as opposed to employing extra qubits in the case of quantum error correction, to estimate the noise-free expectation values from the noisy measurement results.Three of the most well-studied error mitigation techniques are symmetry verification [5,6], quasi-probability and error extrapolation [7][8][9].
Previously all of these error mitigation techniques have been discussed separately.They make use of different in- * cai.zhenyu.physics@gmail.comformation about the hardware and the computation problems to perform different sets of extra circuit runs for error mitigation.Symmetry verification makes use of the symmetry in the simulated system and performs circuit runs with additional measurements.Quasi-probability makes use of the error models of the circuit components and performs circuit runs with different additional gates in the circuit.Error extrapolation makes use of the knowledge about tuning the noise strength via physical control of the hardware, and performs additional circuit runs at different noise levels.Consequently, the three error mitigation techniques are equipped to combat different types of noise with different additional sampling costs (number of additional circuit runs required).Hence, it is natural to wonder how these techniques might complement each other.For NISQ application, it may be essential to understand and develop ways for these error mitigation techniques to work in unison, to achieve better performance than the individual techniques in terms of lower bias in the noise-free expectation values estimates and/or lower sampling costs.Thus one key focus of our Article is on the methods for combining these error mitigation techniques, and trying to gauge their performance under different scenarios through analytical arguments and numerical simulations.
To achieve efficient combinations of these mitigation techniques, we will need to exploit certain features of these constituent techniques.As we will see later, we will show that quasi-probability can be used for error transformation instead of error removal and the circuit runs that fail the symmetry verification can actually be utilised instead of being discarded.It is also essential to understand the mechanism behind error extrapolation, especially in the NISQ limit in which the number of errors in the circuit will follow a Poisson distribution.Heuristic arXiv:2007.01265v2[quant-ph] 4 Mar 2021 arguments and numerical validations have been made by Endo et al. [9] on error extrapolation using exponential decay curves in this NISQ limit.However, it cannot be applied to certain situations arising in practice, for example when the data points have any increasing trend.Our Article will take this further and provide a more rigorous argument showing that single-exponential error extrapolation is just a special case of the more general multi-exponential error extrapolation framework, using which we can achieve a much lower estimation bias.

II. SYMMETRY VERIFICATION
Suppose we want to perform a state preparation and we know that the correct state must follow a certain symmetry S, i.e. we expect our end state to be the eigenstate of S with the correct eigenvalue s (or within a set of eigenvalues {s}).In such a case, we can perform S measurement on our output state and discard the circuit runs that produce states that violate our symmetry.This was first proposed and studied by McArdle et al. [5] and Bonet-Monroig et al. [6].Discarding erroneous circuit runs results in an effective density matrix that is the original density matrix ρ projected into the S = s subspace via the projection operator Π s : Here we have used Π s Π s = Π s .Now let us suppose we want to measure an observable O, which commutes with our symmetry S. Thus they can both be measured in the same run and we will discard the measurement results in the runs that failed the symmetry verification.The symmetry-verified expectation value of the observable O is then: in which we have used [S, O] = 0 ⇒ [Π s , O] = 0. Note that Π s measurement takes the value 1 if the symmetry verification is passed and 0 otherwise and hence Π s = Tr(Π s ρ) is just the fraction of circuit runs that fulfil the symmetry condition.We will use P d to denote the fraction of circuit runs that fail the symmetry verification, which gives Recall that ρ s is the effective density matrix of the nondiscarded runs, which as mentioned is a Tr(Π s ρ) fraction of the total number of runs.Therefore to obtain statistics from ρ s , we need a factor of more circuit runs than obtaining directly from ρ.
In the method discussed above, OΠ s is usually obtained through measuring O and S in the same run.However, sometimes this cannot be done due to, for example, inability to perform non-demolishing measurements.In such a case we need to break OΠ s into its Pauli basis [6] and reconstruct it via post-processing, this is discussed in Appendix A. In this Article, we will mainly be focusing on direct symmetry verification instead of postprocessing verification, but most of the arguments are valid for both methods besides discussions about costs.Now let us move on to see what errors are detectable by symmetry verification.We want to produce the state |ψ f which is known to fall within the S = s symmetry subspace: To produce the state, we usually start with a state |ψ 0 that follows the same symmetry and uses a circuit U that consists of components that conserve the symmetry: Suppose that some error E occurs during the circuit in between the symmetry-preserving components and it satisfies in which s, s are some eigenvalues of the symmetry operator S. We then have: Π s EΠ s = EΠ s means that E is a transformation within the S = s subspace, hence E is undetectable by the symmetry verification using S. Π s EΠ s = 0 means that E contains no components that map states in the S = s subspace back into the same subspace, hence E is (completely) detectable by the symmetry verification using S.
A general error will be a combination of detectable and undetectable error components.
In this Article, we will be focusing on Pauli errors and Pauli symmetries, for which Eq. ( 4) is reduced to:

III. QUASI-PROBABILITY
To describe the quasi-probability method, we will make use of the Pauli transfer matrix (PTM) formalism [10].Using G to denote the set of Pauli operators, any density operators can be written in the vector form by decomposing into the Pauli basis G ∈ G: where we have defined the inner product as: Note that we need to add a normalisation factor 1 / √ 2 N when we use the Pauli operators as basis, with N being the number of qubits.The quasi-probability method was first introduced by Temme et.al.[8], and the implementation details were later studied by Endo et.al.[9].Let us suppose we are trying to perform the operation U, but in practice we can only implement its noisy version In addition to U , we can also implement a set of basis operation {B n }.We can decompose the ideal operation U that we want to implement into a set of gates {B n U } that we can implement: In this way, we are essentially trying to simulate the behaviour of the inverse noise channel E −1 using the set of basis operations {B n }, which can undo the noise E.
If we have a state |ρ passing through the circuit U and we perform measurement O, then the observable we obtain during the experiment will be: in which Q = n |q n | and s n = sgn(q n ).This is implemented by sampling from the set of basis operations {B n } with the probability distribution { |qn| Q }.We will weight each measurement outcome by the sign factor s n and rescale the final expectation value by the factor Q. Now if we break down our computation into many components U = M m=1 U m , with noise associated with each component, then the observable that we want to measure is: but in reality we can only implement the noisy version: Each noise element can be removed by simulating the inverse channels using the set of basis operations {B n }: Hence, we can get back the noiseless observable using in which we have used n to denote the set of number {n 1 , n 2 , • • • , n M } and we have defined To implement this, we simply sample the set of basis operations {B n1 , B n2 , • • • , B n M } that we want to implement with the probability |q n | /Q and weight each measurement outcome by a sign factor s n = sgn(q n ), so that the outcome we get is an effective Pauli observable O Q .And the error-free observable expectation value can be obtained via: Hence, to estimate O by sampling O Q , we need C Q times more samples than sampling O directly, where the sampling cost factor C Q is: In this Article we will be mainly focusing on Pauli error channels, which can be inverted using quasi-probability by employing Pauli gates as the basis operations.
Using to denote a super-operator any Pauli channel can be written in the form: with G α G = 1.We use p to denote the total probability of all the non-identity components.This channel can be approximately inverted using the quasi-probability channel G −p since: Hence, to the first order approximation, the cost of inverting G p will be the cost of implementing G −p , which using Eq. ( 8) is Here we have only discussed approximately inverting a Pauli channel because the exact inverse channel can be hard to express in a compact analytical form.However, it can be obtained numerically by first obtaining the Pauli transfer matrix of the noise channel and then performing matrix inversion.
Instead of removing the error channel completely, quasi-probability can also be used to transform the form of an error channel.In the case of Pauli channels, suppose we want to transform a channel of the form in Eq. ( 9) to we can approximately achieve this transformation up to first order in q and p by applying the quasi-probability channel: This will incur the implementation cost: In the limit of small p and q , we have:

In the last step we have used
If we are suppressing all error components evenly, or if we are simply removing certain error components, we will have p α G ≥ q β G ∀G ∈ G − I.In this case, the cost of implementing the transformation using quasi-probability will simply be:

IV. GROUP ERRORS
Here we will introduce a special kind of error channel: group error channels, which enable us to make more analytical predictions about the error mitigation techniques that we have already discussed and also will help our understanding about error extrapolation later.
The group error J p,E of the group E is defined to be: By groups we mean the subgroups of the Pauli group with a composition rule that ignores all the irrelevant phase factors.For the case of p = 1, we will call J 1,E the pure group errors.Many physically interesting noise models like depolarising channels, dephasing channels, Pauli-twirled swap errors and dipole-dipole errors are all group errors.Now let us consider the effect of applying a set of Pauli symmetry checks S to the group error in Eq. (12).Using Eq. ( 5), S can remove and detect components in J p,E that anti-commute with any elements in S ∈ S. We look at the action of S on the subset of qubits affected by J p,E , and denote the set of these operators on the subset of qubits as S sub .The commutation relationship between S and E is equivalent to that of S sub and E. We denote their generators as S sub and E. Note that here S sub is not a group, by S sub we just mean the set of independent elements in S sub .For Pauli generators, we can choose E in such a way that for every S sub ∈ S sub , there will at most be only one element in E that anti-commutes with it.We will denote the elements in E that commute with all elements in S sub as Q: and it will generate the remaining error components in E that are not detectable, which we denote as Q.Hence the detectable error components are just E − Q Going back to our error channel in Eq. ( 12), the probability that the error gets detected is just the total probability of the detectable error components: Removing the detected errors in J p,E and renormalising the error channel by the factor 1 − p d gives the effective channel after verification, which is just another group error channel with An example of removing detectable errors for depolarising channels will be worked out later in Section VII.For a given general Pauli channel, we have only discussed its approximate inverse channel in Section III.This is because its exact inverse channel can be hard to express in a compact analytical form.However, for any group channel, we can easily write down the explicit form of its inverse channel.
As shown in Appendix B, it is easy to verify that the inverse of a group channel J p,E is just: with α = p 1−p .Using Eq. ( 15) and Eq. ( 8), the cost of using quasiprobability to invert J p,E is thus: which is the same as Eq. ( 10) with As shown in Eq. ( 14), for a given group error J p,E , the resultant error channel after symmetry verification is another group channel J r,Q where Q is a subgroup of E and r = |Q| |E| p.The remaining errors can then be completely removed by implementing J −1 r,Q using quasiprobability.
Similarly if we implement the quasi-probability inverse channel J −1 r,Q first and then perform symmetry verification, we can still completely remove the group error J p,E .The gates we need to implement in the inverse channel J −1 r,Q will not be detected and thus will not be affected by the symmetry verification.As shown in Appendix C, the resultant error channel after applying J −1 r,Q is: This is a channel that only contains the error components that are detectable by the symmetry verification with the error rate p d .As discussed in Section III and explicitly shown in Appendix C, we can implement additional quasi-probability operations to further reduce the error rate of the resultant channel to q ≤ p d .The resultant detectable error channel will be: = (1 − q) I +qV 1 q ≤ p d .

V. NISQ LIMIT
The number of possible error locations in the circuit, which is usually the number of gates in the circuit, will be denoted as M .These error locations might experience different noise with different error rates.From here on, instead of building our discussions around local gate error rates, we will see that the more natural quantity to consider in the context of NISQ error mitigation is the expected number of errors occurring in each circuit run, which is called the mean circuit error count and denoted as µ.In order to achieve quantum advantage using NISQ machine, we would generally expect the circuit size to be large enough to be classically intractable while the mean circuit error count should be not be too far beyond unity: This is also called the Poisson limit since using Le Cam's theorem, the number of errors occurring in each circuit run will follow the Poisson distribution -i.e. the probability that l errors occur will be: If we assume that every local error channel in the circuit can be approximated as the composition of an undetectable error channel and a detectable error channel, then symmetry verification will have no effects on the undetectable error channels and we can focus only on the detectable error channels.Alternatively, as we have seen in Section IV, we can use quasi-probability to remove all the local undetectable errors in the circuit, leaving us with only detectable error channels.The expected number of detectable errors occurring in each circuit run is denoted as µ d .Taking the NISQ limit and using Eq. ( 20), the probability that l detectable errors occur in the circuit is: Using Eq. ( 5), an odd number of detectable errors will anti-commute with the symmetry and get detected while an even number of detectable errors will commute with the symmetry and pass the verification.Therefore, the total probability that the errors in the circuit will be detected by the verification of one Pauli symmetry is thus: Note that this is upper-bounded by 1 2 , i.e. at most we can catch errors in half of the circuit runs.
Hence, using Eq. ( 2), the cost of implementing symmetry verification for one Pauli symmetry is: which is upper-bounded by 2.
After symmetry verification, the fraction of circuit runs that still have errors inside is: In Eq. (11) we have only been focusing on applying quasi-probability to one error channel.Assuming there are M such channels in the circuit, then using Eq.(11) and taking the NISQ limit (Eq.( 20)), the sampling cost factor of transforming all M error channels with error probability p into new error channels with error probability q ≤ p using quasi-probability is: At q = 0, we will obtain the sampling cost of removing all the noise using quasi-probability: Remember that we are focusing on Pauli errors and p is defined to be the total probability of all non-identity components.This is not equivalent to the error probability p because sometimes there are some identity components in our definition of error probability such as the group errors that we discussed in Section IV.Similar to the definition of p , we can denote the expected number of non-identity errors in each circuit run as µ .In the above cases, we have µ = M p and similarly we can define ν = M q .In the cases where different noise locations experience different noise, using Le Cam's theorem with negative probabilities and focusing on the zerooccurrence case, we can generalise Eq. ( 25) to: i.e. the sampling cost of quasi-probability transformation grows exponentially with the reduction in the error rate µ − ν .

VI. MULTI-EXPONENTIAL ERROR EXTRAPOLATION
The idea of amplifying the hardware error rate and performing extrapolation using the original result and the noise-amplified result was first introduced by Li et al. [7] and Temme et al. [8], and was later successfully realised experimentally using superconducting qubits [11].Endo et al. [9] provided heuristic arguments on why the exponential decay curve should be used for error extrapolation in the large circuit limit, whose improved performance over linear extrapolation was validated via numerical simulations in Ref. [9,12].
Using L to denote the set of locations that the errors have occurred, when l errors have occurred in the circuit, our observable O will be transformed to a noisy observable O |L|=l .Recall that in the NISQ limit, the probability that l errors happen (denoted as P l ) will follows the Poisson distribution in Eq. (21).Therefore, the expectation value of the observable O at the mean circuit error count µ is then: Hence, how O µ changes with µ is entirely determined by how O |L|=l changes with l.When we try to fit a nth degree polynomial of µ to O µ , for example performing a linear extrapolation [7,8], we are essentially assuming that O |L|=l is a nth degree polynomial of l using the expressions of the moments of the Poisson distribution.At l = 0, we have the error-free result O : At large error number l, in the case of stochastic errors, the circuit will move closer to a random circuit.Hence, for a Pauli observable O we will expect: A generic polynomial of l will not satisfy the above boundary conditions.Hence, to align with the above boundary conditions, we can instead assume an exponential decay of O |L|=l as l increase: in which γ is the observable decay rate that satisfies 0 ≤ γ ≤ 1.This will lead to an exponential function in µ: which is just the extrapolation curves employed in Ref. [9].Using Eq. ( 29), if we probe at the error rates µ and λµ, we can perform two-point exponential extrapolation and obtain the error-mitigated estimate of O , denoted as O 0 , using the following equation: As discussed in Appendix G 2, the sampling cost factor of performing such an extrapolation is Now let us try to gain a deeper insight about the reason behind the exponential decay of O |L|=l .If a Pauli error G occurs at the end of a circuit and we are trying to measure an Pauli observable O, then the expectation value is just: where η(G, O) is the commutator between G and O: If a pure group error J 1,E occurs at the end of the circuit, then the resultant expectation value is: Using the fact that the composition of the commutators of elements in a Pauli subgroup follows the same structure as the composition of the elements themselves: we can rewrite the above formula in terms of the generator of E: Hence, if a pure group error J 1,E occurs right before measuring a Pauli observable O, then the resultant expectation value Tr(J 1,E (ρ)O) will only remain unchanged if O commutes with all elements in E, otherwise the information about O is "erased" by the group error and the expectation value will be 0.
If we decompose the gates in a unitary circuit i.e. the circuit U can be viewed as the superposition of many Pauli circuits G j .
The expectation value of observing O after applying the circuit U on ρ is which is a linear sum of the measurement results for the set of effective Pauli observables G † i OG j for different i and j.Similar to Eq. ( 31), the information about G † i OG j will either be "erased" or perfectly preserved if a group error occurs in the circuit.Assuming the average fraction of group error locations in our circuit that can erase the information about G † i OG j is γ i, j , then as proven in Appendix D, in the limit of large M and non-vanishing 1 − γ i, j , the expectation value given l errors occurred can be approximated to be: which is just a multi-exponential decay curve.
If we consider the case in which many error locations in the circuit are affected by the same type of noise, and adding onto the fact that in practice there are usually many repetitions of the circuit structures along the circuit and across the qubits, we can expect many γ i, j of different i and j to be very similar.Hence, by grouping the terms with similar γ i, j together, Eq. (33) becomes where A k is the sum of 2 Re α * i α j Tr ρG † i OG j for some subset of i and j.Note that A k are independent of l and we have . So far we have only been considering group error channels.However, as shown in Appendix E, by approximating general Pauli channels as the composition of pure group channels, we can prove that decay of the expectation value under general Pauli noise can also be approximated by a sum of exponentials like in Eq. (34).
Eq. ( 34) can be translated into a multi-exponential decay of O µ over the mean circuit error count µ using Eq. ( 28): This can be rewritten as Comparing with Eq. (34), we see that the shape of O µ and O |L|=l are the same up to the leading order of γ k with the mean circuit error count µ in place of the circuit error count l.
The simplest shape that we can fit over O µ is a single exponential decay curve (K = 1), which is the one we used in exponential extrapolation.We see here that a natural extension of this will be probing O µ at more than two different error rates and trying to fit them using a sum of exponentials (K > 1).The estimate of the errorfree observable O can then be obtained by substituting µ = 0 into our fitted curve.
From Eq. (36), we see that the kth exponential component can only survive up to the mean circuit error count µ ∼ 1 γ k , thus we can only obtain information about this component by probing at the mean circuit error count µ we should be able to retrieve all exponential components if we can probe enough points within µ 1.In practice, there is a minimum mean circuit error count that we can achieve, which we denote as µ * .To have an accurate multi-exponential fitting, it is essential that for all components with non-negligible A k , we have µ * 1 γ k , i.e. none of the critical exponential components has died away at the minimal error rate that we can probe.

VII. NUMERICAL SIMULATION FOR MULTI-EXPONENTIAL EXTRAPOLATION
In this section, we will try to apply multi-exponential extrapolation to the Fermi-Hubbard model simulation circuit as outlined in Ref. [3], which consists of local twoqubit components that correspond to different interaction terms and the fermionic swaps.It can be used for both eigenstate preparation and time evolution simulation.We will assume there are M of these two-qubit components, and they all suffer from two-qubit depolarising noise of error probability p, which is just a group channel of the two-qubit Pauli group (without the phase factors).Using Eq. ( 17), we have Since we usually know the number of fermions in the system, we can try to verify the fermion number parity symmetry of the output state, which is simply: in the Jordan-Wigner qubit encoding.All the local twoqubit components in the circuit conserve the symmetry S σ .Hence when we start in a state with the right fermion number, the output state should also have the correct fermion number, enabling us to perform symmetry verification.
By checking S σ = i Z i , we can detect all error components with one X or Y in the local two-qubit depolarising channels since they anti-commute with S σ .We will remove the other error components in the local twoqubit depolarising channels using quasi-probability.The removed components are those can be generated from the set Q = {Z 1 , Z 2 , X 1 X 2 } following Section IV.Thus we have |Q| = 2 3 = 8 and using Eq. ( 13) we also have The resultant noise channel after the application of quasiprobability is given by Eq. ( 19), which is just a uniform distribution of the two-qubit Pauli errors that are detectable by S σ .We will call it the detectable noise.
In this section and later in Section IX, we will perform numerical simulations using the circuit for the 2 × 2 half-filled Fermi-Hubbard model, which consists of 8 qubits and 144 two-qubit gates.The two-qubit gates in the circuit that correspond to interaction terms are parametrised gates with the parameters indicating the strength of the interaction.In our simulation, we will obtain the results for a set of randomly chosen gate parameters (with additional results for another set of random parameters listed in Appendix I).We will also look at two different error scenarios: depolarising errors and detectable errors.One of them is a group channel while the other is a more general Pauli channel.The measurements that we perform will be the Pauli components of the Hamiltonian, from which we can reconstruct the energy of the output state.The simulations are performed using the Mathematica interface [13] of the high-performance quantum computation simulation package QuEST [14].
In Fig. 1, we have plotted the noisy expectation values O µ for each Pauli observable at the mean circuit error counts µ = 0.5, 1, 1.5, 2. When we perform multiexponential extrapolation on them, we find that all of the observables can be fitted using a sum of at most two exponentials, even though there should be very few symmetries in our circuits since they are generated from a set of random parameters.We now can proceed to compare the absolute bias in the estimate of the noiseless expectation values (µ = 0) using dual-exponential extrapolation against that using the conventional single-exponential extrapolation for all of the noisy observables in Fig. 1.The absolute biases of single-and dual-exponential extrapolation are denoted as 1 and 2 , respectively, and different colours in the plots correspond to different estimation bias ratios 1 / 2 .
For 32 out of the 34 observables we plotted, dualexponential extrapolation can achieve a smaller estimation bias than single-exponential extrapolation ( 1 / 2 > 1).Within the two cases that dual-exponential extrapolation is outperformed (the green curves in Fig. 1), one of them is the case in which single-and dual-exponential extrapolation both achieve very small estimation bias of similar order.The other remaining case with larger 2 relative to 1 is mainly due to the small magnitude of its true expectation value, which lead to large uncertainties in the fitting parameters.It is also because we are only using the bare minimum of 4 data points to fit a dualexponential curve with 4 free parameters and thus the problem may be alleviated by simply probing at more error rates to obtain more data points.On the other hand, there are also a few cases in which the 1 are exceptionally large (e.g.certain orange and red curves in Fig. 1 (a) and (b)).These are usually observables whose decay curves have extrema and/or crossing over the xaxis, thus it is impossible to get a good fit with a singleexponential curve.We have zoomed into one such observable in Fig. 2. For these observables, dual-exponential extrapolation can still perform extremely well and achieve  exponential extrapolation.This is shown in Table I, from which we see that by using dual-exponential extrapolation instead of single-exponential extrapolation, we can achieve tens or even a hundred times reduction in estimation bias across both noise models.Note that in Fig. 1 it appears to the eye that the true (noiseless) expectation values, marked by filled circles, never deviate from the dual-exponential (dashed) lines.In fact there are minute discrepancies as specified in Table I, but the extrapolation is remarkably successful.

VIII. COMBINATION OF ERROR MITIGATION TECHNIQUES
We have shown that extrapolating using multiexponential curve can be very effective assuming Pauli noise.Besides the shape of the extrapolation curve, the other key component to error extrapolation is the way to tune the noise strength.Previously, noise are boosted by increasing gate pulse duration [8], applying additional gates that cancel each other [12] or simulating the noise using random gate insertion [7].There are various practical challenges associated with these noiseboosting techniques, and furthermore as discussed at the end of Section VI, data points at boosted error rates may not contain enough information for effective extrapolation.Now as shown in Section III, we can actually use quasi-probability to reduce the error rate and obtain a set of data points with reduced noise strength, which can then be used for extrapolation.We should expect a smaller estimation bias by using data points with re-duced noise strength rather than boosted noise strength, but the sampling cost will also increase due to the use of quasi-probability.For the special case of two-point extrapolation using a single-exponential curve with the two data points at the unmitigated error rate µ and the quasi-probability-suppressed error rate ν, the total sampling cost as shown in Appendix G 3 is with µ = λν.We will call this special case quasiprobability with exponential extrapolation.
As discussed in Section VI, the number of exponential components in the multi-exponential extrapolation can be reduced with increased degree of symmetry in the circuit and/or if the error channels are group channels.Hence, besides using quasi-probability for error suppression, we can also use quasi-probability to transform the error channels into group errors and/or errors of similar form for easier curve fitting in the extrapolation process.
Moving on, we may wish to combine symmetry verification and quasi-probability.We can first apply quasiprobability to transform all the error channels in the circuit with the total mean error count µ into detectable error channels with a total mean error count µ d .After that, we can apply symmetry verification, but note that the additional quasi-probability operations may contain gates that take us from one symmetry space to another, for which we need to adjust our symmetry verification criterion accordingly.As discussed in Section V, an even number of occurrence of local detectable errors can still escape the symmetry test and lead to circuit errors.We can further suppress them by applying additional quasiprobability operations as shown in Appendix F. Alternatively, we can also try to remove these remaining errors by applying error extrapolation as we will see below.
After using quasi-probability to transform all local errors into local detectable errors with the mean circuit error count µ d , performing symmetry verification will split the circuit runs into two partitions, one has an even number of detectable errors occurring and will pass the symmetry test, the other has an odd number of detectable errors occurring and will fail the symmetry test.Consequently, the noisy observable expectation value (Eq.( 35)) can also be split into the weighted sum of these two partitions: (41) We will consider the case that the decay of our expectation value O µ d over increased mean circuit error count µ d follows a single exponential curve (K = 1) for simplicity, we then have: which gives: (43) Note that we have used the fact that O c,µ d and O have the same sign since 1−γ > 0 and µ d > 0.Here with the help of symmetry verification and quasi-probability, we can now obtain an estimate of the error-free expectation value O by combining the expectation value of the passed runs and failed runs at one error rate µ d instead of combining the expectation value of runs at different error rates in the conventional error extrapolation.Note that here we have assumed that we know the value of the mean detectable circuit error count µ d , which needs be known before we can apply the quasi-probability step anyway.The method we employed in Eq. (43) will be called hyperbolic extrapolation.
As derived in Appendix G 4, the sampling cost factor of hyperbolic extrapolation is To combine all three error mitigation techniques, we first use quasi-probability to remove the error components that are undetectable by symmetry verifications.Applying symmetry verification will then split the circuit runs into two sets: runs with even number of errors and runs with odd number of errors, obtaining two separate expectation values.Using our understanding about the decay of the expectation value from our study of error extrapolation, we can simply combine these two erroneous expectation values and obtain the error-free expectation value.The full process is called quasi-probability with hyperbolic extrapolation, and the corresponding total sampling cost factor can be obtained using Eq. ( 27) and Eq. ( 44): We note that this is always smaller than the cost of pure quasi-probability C Q,0 = e 4µ .Its cost saving over pure quasi-probability will increase with the increase of γ.

FIG. 3.
Comparison of the biases in the error-mitigated expectation values between quasiprobability with exponential extrapolation (QE) and quasi-probability with hyperbolic extrapolation (QH) in a 8-qubit simulation.Plots showing (a) observables following single-exponential decay and (b) observables following dual-exponential decay.Within each plot, different colours represent different observables.The solid lines denote QH, while the dashed lines denote QE.At the mean circuit error counts µ = 2, for each observable, we use markers to denote the method that has lower estimation bias out of the two.For a given observable, circular markers denote lower estimation bias when using QH, while triangle markers denote lower estimation bias when using QE.
In this section, we will compare the performance of quasi-probability with exponential extrapolation (QE) and quasi-probability with hyperbolic extrapolation (QH) discussed in Section VIII.Similar to Section VII, we will perform Fermi-Hubbard model simulation with local two-qubit depolarising noise with a mean circuit error count µ.The symmetry we will used in QH is the fermionic number parity symmetry, which means that the resultant mean detectable circuit error count after we apply the quasi-probability step in QH will be µ d = µ 2 following Eq.(38).For the quasi-probability in QE in this section, we will keep it at the same strength as that in QH, which means that they have the same resultant circuit error rate: ν = µ 2 .Note that even though resultant channels after the partial quasi-probability in both QE and QH give the same mean circuit error count ν = µ 2 = µ d , in one case the resultant noise is still depolarising while in the other case the resultant noise is locally detectable.In this section, we will assume the quasi-probability process is performed perfectly.Recall that for simplicity we have only explicitly derived QH under the assumption that the observable follows a single-exponential decay, so for a fair comparison the QE method in this section will also only employ singleexponential extrapolation.However as we will see later, even when we look at observables that follow a dualexponential decay, which breaks our assumptions above, QH can still achieve robust performance.
As shown in Fig. 1, for our example circuits, some observables can be fitted well enough using singleexponential decay curves while the other observables can only be fitted well using dual-exponential decay curves.We will call these two types of observables single-exponential observables and dual-exponential observables, respectively.In Fig. 3, we have plotted the absolute estimation bias est using the two different extrapolation techniques for the single-exponential and dualexponential observables.First, we can see that the estimation biases for the dual-exponential observables are almost one order of magnitude higher than the biases for the single-exponential observables.This should not come as a surprise since both single-exponential extrapolation and hyperbolic extrapolation are derived under the assumptions of single-exponential observables.At the mean circuit error counts µ = 2, for each observable, we have used markers to label the method that can achieve a lower estimation bias out of the two.We see that the number of single-exponential observables that can achieve a lower estimation bias using QE is comparable to that of QH.On the other hand, almost all dualexponential observables can achieve a lower estimation bias using QH.
In Table II, we further calculate the average estimation bias of single-exponential, dual-exponential and all observables separately at the circuit error rate µ = 1, 2, which re-confirm all of our observations above.We see that the estimation bias of QE is lower than that of QH for single-exponential observables, and on the other hand QH can achieve a lower estimation bias for dual-est/10   exponential observables.In other words, the performance of QH is more robust against whether the observable is single-exponential or not.When looking at the estimation bias averaged over all observables, we see the estimation bias of QH is always lower than QE and can be 4 times smaller than QE at µ = 2.The all-observable averages can be more indicative about the practical performance of the mitigation techniques since in experiments we do not know whether a given observable should be fitted with single-exponential or not beforehand.
There is another added layer of robustness when we try to apply QH instead of QE to multi-exponential observables when we look back at the hyperbolic extrapolation equation Eq. ( 43).We can see that if the shape of the observable is far off from a single-exponential decay, then this might lead to a negative number in the square root of Eq. ( 43), allowing us to realise that we need to probe at more error rates to perform multi-exponential extrapolation instead, and avoiding performing a bad extrapolation with a very large bias.In the simulation, we indeed identify a few observables that we cannot perform QH on.For these observables, we can still perform QE, but it will lead to huge biases in the estimates.These observables have been excluded in our comparison between QE and QH.
Using Eq. ( 45), Eq. (37) and Eq. ( 38), the sampling cost factor of performing QH in our example circuit is where γ is the decay rate of the observable expectation values under noise.Using Eq. (39), Eq. ( 37) and ν = µ d = µ 2 ⇒ λ = 2, the sampling cost factor of performing QE is: For comparison purpose, we also write down the sampling cost factor for removing all the errors using quasi-probability given by Eq. ( 26): The comparison between C Q,0 , C QE (γ) and C QH (γ) at different γ is plotted in Fig. 4. We can see that C QH (γ) is always lower than C Q,0 across all µ and γ, i.e. we can always get a sampling cost saving by applying QH instead of pure quasi-probability, which is also proven in Section VIII.On the other hand, at γ = 1, C QE is larger than C QH for all µ and larger than C Q,0 for µ < 2.4.As γ decreases, C QH (γ) will increase while C QE (γ) will decrease.Thus they naturally complement each other as QH will be more suitable for large-γ error mitigation while QE will be more suitable for small-γ error mitigation.At γ = 0, we see that C QE becomes lower than both C QH and C Q,0 at µ > 1.8.
The average fitted γ of all the single-exponential observables within each plot in Fig. 1 all lie within the range 0.5 − 0.6.Hence, we will now focus on the γ = 0.5 plot in Fig. 4 to get an indication of the practical sampling costs of implementing different mitigation techniques.
At µ = 1, the sampling cost factor of quasi-probability is 43.QE requires a higher sampling cost, thus there is no point performing QE since pure quasi-probability can remove all the noise perfectly in theory with a lower cost.Compared to quasi-probability, QH can reduce the cost by more than 70% while still achieving the small estimation bias QH ∼ 3 × 10 −3 shown in Table II.In order for quasi-probability to have any advantages over QH, we must sample enough times such that the shot noise of pure quasi-probability is smaller than the estimation bias of QH (more rigorous arguments in Appendix H), which will require N * ∼ C Q,0 / 2 QH ≈ 4.3 × 10 6 samples for each observable.Therefore in practice, QH could be the preferred method over pure quasi-probability as it is challenging to sample more than N * for each observable within reasonable runtime constraints [3].
At µ = 2, now the sampling cost factor of quasiprobability is 1800, which is hardly practical.QE can reduce this sampling cost by half while achieving estimation bias around 4 × 10 −2 (Table II), and QH can reduce this sampling cost by almost 90% while achieving estimation bias around 1 × 10 −2 (Table II), thus they both would be preferred over pure quasi-probability in practice following similar arguments in the µ = 1 case.We also see that QH outperforms QE in terms of both sampling cost and estimation bias at µ = 2, and thus would be preferred over QE.The cost of QE will only become lower than QH at µ = 3.9, however at this point, neither of their sampling costs are likely to be practical.

X. DISCUSSION
In this Article, we have recapped and studied the mechanism and performance of three of the most wellknown error mitigation techniques: symmetry verification, quasi-probability and error extrapolation under Pauli noise.By introducing the concepts of group errors and NISQ limits, we managed to prove that the change of the expectation value of a Pauli observable with increased Pauli noise strength can be approximated using multi-exponential decay, enabling us to extend exponential error extrapolation to multi-exponential extrapolation.We then performed 8-qubit numerical simulations for Fermi-Hubbard simulation under two different Pauli noise models, finding that the decay of their Pauli expectation values can all be fitted using single-or dualexponential curves, confirming our earlier proof of multiexponential decay.Using the same circuits, we performed dual-exponential extrapolation by probing at four differ-ent error rates, which is minimal number of data points required, and managed to obtain a low estimation bias of 10 −4 for almost all 34 observables except for one fringe case.In our simulations, the estimation bias of dualexponential extrapolation are on average ∼ 100 times lower than that of single-exponential extrapolation, with the maximum factor of bias reduction reaching ∼ 10 4 .
We then proceeded to combine different error mitigation techniques in the context of well-characterised local Pauli noise.Instead of using quasi-probability to completely remove all the noise, we can use it to suppress the noise strength and perform error extrapolation, which is named quasi-probability with exponential extrapolation (QE).Alternatively, we can use quasi-probability to remove the local undetectable noise and perform symmetry verification.Then instead of discarding all the circuit runs that fail the symmetry test, we have developed a way to recombine the expectation values of the 'failed' and 'passed' runs to obtain an estimate of the noiseless observable.The full combined method is called quasiprobability with hyperbolic extrapolation (QH).Note that both QE and QH are free of the requirement to adjust the hardware error rate despite the name 'extrapolation'.By performing 8-qubit Fermi-Hubbard model simulations under local depolarising noise and using the fermionic number parity symmetry, we found that QH outperforms QE in terms of both estimation bias and sampling costs for almost all cases.When compared to pure quasi-probability, QH can achieve factor-of-4 and factor-of-9 sampling cost savings at the mean circuit error count µ = 1 and µ = 2, respectively, while still maintaining low estimation bias of 10 −3 ∼ 10 −2 .Hence, QH would outperform pure quasi-probability in our examples unless we obtain an impractical number of samples (more than millions) per observable.
QH is derived under the assumption that the observables decay along single-exponential curves with increased noise.Our simulation shows that QH can be robust against violation of this assumption when applied to dual-exponential observables.However, such robustness may not persist with a further increase in the number of exponential components.A multi-exponential version of QH can be done through probing at more error rates and fitting Eq. (41) to the data.Alternatively, instead of probing at more error rates, we can also try to verify more symmetries.In such a way, we can obtain a set of expectation values corresponding to different verification syndromes for the multi-exponential version of the hyperbolic fitting.An example can be using the separate fermion number parity symmetries for each spin subspace, which will lead to expectation values corresponding to the four possible verification syndromes.However, how to recombine these expectation values in the case of multiple symmetries and how to use quasi-probability to transform the local error channels into the suitable forms for such a recombination is not a simple extension of the single-symmetry case we considered.
In our derivation, the number of exponential components in the expectation value decay curve in Eq. ( 34) is expected to scale exponentially with the number of gates.However, in our simulations, we fitted at most two exponential components for each of the observable decay curves.More analysis is needed to bridge the gap between the expected and the actual number of exponential components required, possibly based on the symmetry of the circuit.This will help us understand how the number of exponential components scales with the system size, enabling us to gauge the performance and the costs of scaling up the multi-exponential extrapolation method.It might be useful to draw ideas from non-Clifford randomized benchmarking [15][16][17], in which multi-exponential decay is also employed for the fitting of the fidelity curves.When applying multiexponential extrapolation in practice, we might want to develop Bayesian methods to determine whether we need to probe at more error rates, which error rates to probe, and whether to change the number of exponential components of our fitted curve based on the existing data.This has been done in the context of randomized benchmarking [18] and it would be interesting to see its performance in the context of multi-exponential extrapolation.
One combination of error mitigation techniques that we have not explored here is pairing symmetry verification with error extrapolation without using quasiprobability.The naive version of such a combination is discussed in Ref. [3].To make use of the results in this Article, one possible way is to approximate all the local error channels as the compositions of detectable and undetectable error channels, so that we can deal with them separately using hyperbolic extrapolation and exponential extrapolation.It would be very interesting to see the implementation details of such a method and how it compares to pure error extrapolation.
We have only considered Pauli noise in this Article, thus it will also be interesting to see whether our arguments can be extended to other error channels like amplitude damping or coherent errors.In practice, we can transform any error channels into Pauli channels using Pauli twirling [19,20] and then apply our methods.Note that we can even perform further twirling like Clifford twirling to transform the error channels into group channels, which can be better mitigated as we have observed.Ways to transform a given error channel into a group channel can be an interesting area of investigation.

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.

CODE AVAILABILITY
The code used in the current study are available from the corresponding author upon reasonable request.
we have: where D is the group generated by the union of of the generators of E and B: i.e.D is the group of elements that one can obtained by composing the elements in E and B.
In the case of E and B have no overlaps, i.e.E ∩ B = 0 (note not B here since there is a degree of freedom in choosing B), we have D = E + B. We will have the same result if we have E ∩ B = 0 instead.In such as case, for the subgroup E of D, B is the corresponding quotient group, and vice versa.
In the case of E is a subgroup of B, we have which leads to Using Eq. (B3) we then have: which is still a group error of the same group E with a modified error probability.
For J q,E to be the inverse of J p,E , we require: Hence, Appendix C: Removing Subgroup Components in a Group Channel

Form of the quasi-probability channel
For a group error channel of the group E, we would want to remove the error components Q in the channel where Q is a subgroup of E. i.e.Thus we want to transform the channel: Using Eq. (B2), we have along with Eq. (B4), the channel we need to perform to transform J p,E to V will be: Using Eq. (C1), we have: It turns Eq. (C2) into: |E| p as defined in Eq. ( 13), which is the probability of the errors in V ∈ E − Q occurring.
At q = p d , the channel we need to apply is just: i.e. we only need to apply components in Q in our quasiprobability channel.And using Eq.(C1), our resultant channel is simply: same as what we derived in Appendix F.

Cost of the quasi-probability channel
Decomposing the quasi-probability channel that we need to implement into its components: and compared to Eq. (C3), we have: The normalisation factor of implementing the quasiprobability channel VJ −1 p,E is: where we have used E∈E β E = 1.It means that the cost of implementing the quasiprobability channel is split into two case by the threshold error rate q = p d of the final channel V q .
Targeted error rate above or at the threshold This means: Now using Eq.(C5), we have The cost to implementing the transformation channel V q J −1 p,E is then i.e. the cost to implementing the quasi probability channel in this regime is independent of q.However, the exact quasi-probability channel we need to implement will still change with q.Uses Eq. ( 17) and Eq. ( 13), we have: Target error rate below the threshold This means: Now using Eq.(C5), we have The cost to implementing the transformation channel V q J −1 p,E is just where we have use the expression of p in Eq. ( 17).Here 1 + 4p is the cost to invert the entire channel J p,E .For each bit of remaining error probability q in our final channel V q , we can reduce the cost factor by 4q until we hit the threshold q = p d .This is the same as what we obtained in Eq. (11).Now if a pure group error J 1,E f happens between V f and V f +1 , and we denote then the expectation value of a Pauli observable O is: where in the last step we have use Eq.(31) and we have defined: Here X f ( i, j) takes the value 0 if the effective observable G † i f +1:M OG j f +1:M right after the noise does not commute with all elements in E f , and takes the value 1 otherwise.For the former case, we will say the information about G † i OG j is "erased" by the error at location f .Note that the effect of the noise does not relate to the gates implemented before it at all.For a given effective observable G † i OG j , the fraction of error locations that can "erase" its information is simply: Denoting O L as the expectation value we obtain when the set of error location is L, we have: Recall that O |L|=l is the expectation value we obtain when there are l errors in the circuit, regardless of the error location.By definition we have: As proven in Appendix D 2, in the limit of large M and non-vanishing 1 − γ i, j , we have: Thus Eq. (D1) can be approximated as: where we have use γ i, j = γ j, i from the definition of γ i, j .

Expansion of the sum of Bernoulli samples
We denote X f as the f th sample from a Bernoulli distribution.From the total M samples taken, we can estimate the success probability denoted r:

Now we define
Then Y 1 is just the sum of these samples: Using multinomial expansion, we have where the multinomial coefficient is f , which is the number of ways the distribute l distinct balls (the total power of l) into M distinct bins (M terms of X f for f ∈ {1, 2, 3, • • • , M }) such that the number of balls in bin f is n f (the power of X f is n f ).Now using X n f = X f for any n, we have where we have used l l = 1.
Hence, in the limit of large M and assuming the success sample fraction r is non-vanishing, we have: From the definition of Y l we have: The estimated failure probability of the Bernoulli distribution is: which gives: After using quasi-probability to remove all the local undetectable errors in the circuit, which reduces the mean circuit error count from µ to µ d as discussed in Section V, we can further suppress the remaining errors to achieve a mean circuit error count of ν ≤ µ d and follow the same arguments in Section V with ν in place of µ d .Using Eq. ( 27), the sampling cost factor of the required quasi-probability transformation is: where we have made use of ν = ν for the detectable error channels like V q in Eq. ( 19).
In the resultant circuit with only local detectable errors present, the sampling cost of applying symmetry verification is given by Eq. ( 23) with ν in place of µ d : .
Hence, the total sampling cost including quasiprobability is: which is smaller than the cost of using quasi-probability to remove all of our errors without the help of symmetry verification in Eq. ( 26 However, we need to note that with the pure quasiprobability method, we can remove all of the noise, while when we combine with symmetry verification, the fraction of erroneous circuit runs after applying our combined error mitigation is given by Eq. ( 24) with ν in place of µ d : Hence, we must choose a ν such that P circ is the circuit error rate that we can tolerate.The saving of sampling cost when combining symmetry verification with quasi-probability over pure quasiprobability can be written as: To achieve P circ 1, we have: Therefore by combining with symmetry verification, the factor of saving in the sampling cost over pure quasiprobability is C Q,0 /C QS,ν = 1.5 for the circuit error rate P circ = 10 −2 and C Q,0 /C QS,ν = 1.15 for the circuit error rate P circ = 10 −3 .To push the circuit error rate any lower, the saving in the sampling cost will be negligible.And in the limit of P circ = 0 ⇒ ν = 0, we just have pure quasi-probability which removes all the errors.We can of course apply hyperbolic extrapolation on top of this, which is the same as Section VIII, but with ν in place of µ d .The cost function is thus similar to Eq. ( 45) We need to note that the cost expression here is only valid for 0 ≤ ν ≤ µ d following from Eq. (F1), thus the minimum cost that we can achieve will be at ν = µ d .
The only reason to try to push ν below µ d is when such an action will result in easier (smaller K) and/or better fitting to our extrapolation curves.
Appendix G: Cost of Error Extrapolation

Cost for two-point extrapolation
Suppose for an observable O, we can estimate its expectation value by combining the expectation value of two other observables A and B. Then by denoting the estimate as O 0 , we have: for some estimation function f .Now suppose we take N samples in total and α is the fraction of A samples within, then using A, B to denote the sample averages, we can obtain the sample estimate of O 0 (and thus O ): Note that we have abuse the notation here since O 0 and O 0 are not the expectation value and the sample average of some observable O 0 .There is no such an observable.Rather, O 0 is exactly defined as the estimate of O after N total samples of A and B using the equation above, and similarly O 0 is defined as the case N → ∞.Note that O 0 ≡ O as well since O is an actual observable of the noiseless computation that lead us to O .We will try to compare the variance of the sample estimate O 0 against the variance of the noiseless sample average O.The variances of various sampling averages follow these equations:

Hence, assuming Var
, the variance of the sample estimate is: Hence, when sampling for O 0 instead of O, the additional cost factor is: We then have the extremal point being α ± = a a ± b .
Since α − = a a−b ≥ 1, we will only keep α + .We can also obtain C (α + ) ≥ 0, which means that it is a local minimum, whose value is ) i.e. the minimal cost factor is achieved when the fraction of A samples is α + .
For the naive version of evenly distributing all the samples: α = 0.5, we have When compared to the optimal distribution, we have: i.e. the saving in the number of samples by using optimal sample distribution will be less than half.The saving will be exactly half in the case of a b or b a.In practice, a and b is often unknown and thus it is hard to achieve the optimal sample distribution.Hence, in most of the cases in this Article, we just use the naive sample distribution, which should be of the same order of magnitude as the optimal case.

Cost of exponential extrapolation
For an observable O, we can estimate its expectation value by probing at the error rate µ and λµ and fitting with an exponential curve to get the zero-noise value: Note that O will only be exactly the same as O 0 if the observable follows strictly a single exponential decay with increased noise.Now following Appendix G 1 with A = O µ , B = O λµ and using O µ ≈ O 0 e −γµ , we have: Hence, the sampling cost factor of exponential extrapolation can be obtained using Eq.(G5): (G8)

Cost of quasi-probability with exponential extrapolation
As outlined in Section VIII, we use quasi-probability to suppress the error rate from µ to ν = µ λ and then use the point at µ = λν and ν to perform exponential extrapolation.Thus we have the same equation as Eq.(G6) with µ → ν: Here λν = µ is the original error rate and thus O λν will be obtained via direct sampling, while ν is the quasiprobability suppressed error rate, thus O ν will be obtain via quasi-probability.Thus using Eq. ( 7), we have O λν and using Eq.(G7) with ν in place of µ, we have: Using Eq. ( 8), we have Q = C Q,ν .From Eq. ( 27), we have C Q,ν = e 4(µ −ν ) .Since in extrapolation, we need to suppress all error components evenly, we have µ = λν just like the relation between µ and ν.Hence, we have: Hence, the sampling cost factor of quasi-probability with exponential extrapolation can be obtained using Eq.(G5): The fraction of O c,ν samples among all the samples can be obtained from Eq. ( 40): α = e −ν cosh(ν).

Appendix H: Comparison between Error Mitigation Techniques
For a given error mitigation technique whose average estimation bias for a Pauli observable is and whose sampling cost is C, we can define ∆ to be the effective observable that represents the error in each circuit run for obtaining the error-mitigated Pauli observable, which has: In another word, is the systematic bias of the error mitigation techniques while C is the strength of the random errors (shot noise) of the technique.
We will label the average error after N sample runs as ∆, for which we have: Hence, the expected square errors of applying the error mitigation technique with N samples is thus: Now suppose we have two different error mitigation techniques, and w.l.o.g.we will assume technique 1 will have lower systematic bias 1 ≤ 2 .If C 1 < C 2 , then technique 1 is obviously the better technique.If C 2 ≤ C 1 however, there is a break-even point which both methods have similar mean square errors: .
i.e.N * is roughly the number of samples needed using technique 1 to reach a shot noise level that is equal to the systematic bias of technique 2 ( 2 ).When the number of samples N N * , the shot noise will dominate and thus we will choose technique 2 over technique 1 due to the lower sampling cost.When the number of samples increase and reach N N * , then the systematic bias will dominate over shot noise and thus we will choose technique 1 over technique 2 due to the lower systematic bias.
Appendix I: Numerical result for another set of random parameters Here we repeat the numerical simulation in Section VII and Section IX, but with a different set of random circuit parameters.Fig. 5 and Table III show the comparison between single-exponential extrapolation and dualexponential extrapolation, which follow similar trends as we have discussed in Section VII.Fig. 6 and Table IV show the comparison between quasi-probability with exponential extrapolation and quasi-probability with hyperbolic extrapolation, which follow similar trends as we have discussed in Section IX.FIG. 6.Comparison of the biases in the error-mitigated expectation values between quasi-probability with exponential extrapolation (QE) and quasi-probability with hyperbolic extrapolation (QH) for (a) observables following single-exponential decay and (b) observables following dualexponential decay in a 8-qubit simulation.Within each plot, different colours represent different observables.The solid lines denote QH, while the dashed lines denote QE.At the mean circuit error counts µ = 2, for each observable, we use markers to denote the method that has lower estimation bias out of the two.For a given observable, circular markers denote lower estimation bias when using QH, while triangle markers denote lower estimation bias when using QE.

FIG. 1 .
FIG. 1.Comparison between single-exponential extrapolation and dual-exponential extrapolation in a 8qubit simulation.Plots showing the noisy expectation values of different Pauli observables under (a) depolarising noise and (b) detectable noise obtained at the four mean circuit error counts µ = 0.5, 1, 1.5, 2 (cross markers).The singleand dual-exponential extrapolation curves fitted to the data points are represented by the solid and dashed lines, respectively.The circular markers lie at µ = 0 and denote the true noiseless expectation values.Different colours represent different ratios between the estimation bias of using single-and dual-exponential extrapolation (e.g.orange means that for the given observable, the estimation bias of single-exponential extrapolation is between 10 3 to 10 4 times larger than that of dual-exponential extrapolation).

FIG. 2 .
FIG. 2.A noisy Pauli observable from Fig.1(b) that cannot be fitted well using single exponential decay since it crosses over the x-axis.Here we have plotted the noisy expectation values at the four mean circuit error counts µ = 0.5, 1, 1.5, 2 (cross markers).The single-and dualexponential extrapolation curves fitted to the data points are represented by the solid and dashed lines, respectively.The circular markers lie at µ = 0 and denote the true noiseless expectation values.
) in which e −µ d cosh(µ d ) and e −µ d sinh(µ d ) are the probability to have even and odd number of errors occurring in the circuit, respectively.And O c,µ d , O s,µ d are the corresponding expectation values in these cases with (a) µ = 1 and (b) µ = 2.The entries in (a) and (b) are in the units of 10 −4 and 10 −3 , respectively.

FIG. 4 .
FIG. 4.Comparison of the sampling cost factors between pure quasi-probability (Q), quasi-probability with exponential extrapolation (QE) and quasi-probability with hyperbolic extrapolation (QH).Plots showing different noisy observable decay rate (a) γ = 1, (b) γ = 0 and (c) γ = 0.5.We have labelled the values of the lines at the mean circuit error counts µ = 0.5, 1, 2, 4. We have also labelled the intersects between lines of different methods using red markers.

Appendix D : 1 m=M
Decay of Pauli Expectation Value under Group Channels 1. Overall derivationIn Section VI, we are looking at a circuit of the form U = V m with each gate decomposed into their Pauli components: V m = jm α mjm G mjm .
b where L is the subset of non-empty bins and l b is the Stirling number of the second kind, which is the number of ways to distribute l distinct balls into these b identical bins such that none are empty.And l b b! is just the number of ways of distribute l distinct balls into these b distinct bins such that none are empty.Hence,

FIG. 5 .
FIG. 5. Comparison between single-exponential extrapolation and dual-exponential extrapolation under (a) depolarising noise and (b) detectable noise in a 8-qubit simulation.Plots showing the noisy expectation values of different Pauli observables obtained at the four mean circuit error counts µ = 0.5, 1, 1.5, 2 (cross markers).The single-and dualexponential extrapolation curves fitted to the data points are represented by the solid and dashed lines, respectively.The circular markers lie at µ = 0 and denote the true noiseless expectation values.Different colours represent different ratios between the estimation bias of using single-and dualexponential extrapolation

TABLE II .
The biases in the error-mitigated estimates using QE and QH averaged over single-exponential, dualexponential and all observables at the mean circuit error counts