Practical Quantum Error Mitigation for Near-Future Applications
Suguru Endo Simon C. Benjamin Ying Li
quant-ph








Department of Materials, University of Oxford, Oxford OX1 3PH



Graduate School of China Academy of Engineering Physics, Beijing 100193, China
Department of Materials, University of Oxford, Oxford OX1 3PH


It is vital to minimise the impact of errors for near-future quantum devices that will lack the resources for full fault tolerance. Two quantum error mitigation (QEM) techniques have been introduced recently, namely error extrapolation <cit.> and quasi-probability decomposition <cit.>. To enable practical implementation of these ideas, here we account for  the inevitable imperfections in the experimentalist's knowledge of the error model itself. We describe a protocol for systematically measuring the effect of errors so as to design efficient QEM circuits. We find that the effect of localised Markovian errors can be fully eliminated by inserting or replacing some gates with certain single-qubit Clifford gates and measurements. Finally, having introduced an exponential variant of the extrapolation method we contrast the QEM techniques using exact numerical simulation of up to 19 qubits in the context of a `SWAP test' circuit. Our optimised methods dramatically reduce the circuit's output error without increasing the qubit count or time requirements.


Practical Quantum Error Mitigation for Near-Future Applications
    Ying Li
    December 30, 2023
===============================================================





§ INTRODUCTION


Controlling noise in quantum systems is crucial for the development of practical technologies. Such noise can occur due to unwanted interactions of a passive qubit with the environment, or due to imperfections in the use of circuit elements that compose the algorithm (qubit initialisation, gates, and measurement). In all cases the result is errors occurring at the level of physical qubits. The theory of quantum fault tolerance (QFT) reveals that the introduction of logical qubits, composed of numerous physical qubits, can allow one to detect and correct errors at the physical level; however this capacity comes at an enormous multiplicative cost in resources. A recent estimate suggests that a Shor algorithm operating on a few thousand logical qubits would require several million physical qubits <cit.>. While it is encouraging to know that such techniques exist, hardware on this scale is probably at least a decade away. The timely (indeed, urgent) question is, to what extent can we control the impact of errors in computing devices that are too small to support full QFT?

It may prove to be the case that deep quantum algorithms, such as Shor's factoring algorithm and Grover's search algorithm, cannot be successfully executed on classically-intractable problems without the support of QFT. However, fortunately there are other algorithms of potential practical significance that focus on shallow circuits, with the output typically being fed into a classical supervising algorithm so as to form a hybrid system. Such approaches have been proposed for the simulation to aid discovery in chemistry and materials science; see Refs. <cit.> for examples. Hybrid systems may be capable of yielding significant results, surpassing conventional computers, even when finite error rates are present because of their error resilience <cit.>. In order to achieve this it is desirable to suppress or mitigate errors to the greatest extent possible while keeping the qubit count ideally unchanged, or increasing only modestly compared to the high cost of full QFT. 

Recently two techniques were introduced for quantum error mitigation (QEM) in generic hybrid quantum algorithms where the expected value of an observable – say, a z-basis measurement of a given qubit – is the quantity of interest. The goal is to estimate the value that this observable would take given an error-free circuit, despite the reality that the real experimental system cannot perform operations with less than a certain error rate. Ref. <cit.> introduced a hybrid algorithm simulating quantum dynamics, which featured an active error minimisation technique involving extrapolation. The experimentalist would execute the circuit with all errors at their minimum achievable severity, obtain the expected value of the observable, and then repeat the exercise having deliberately increased the physical error rate (or having applied additional quantum gates to achieve the same effect). By noting the effect of the increased errors on the observable, the experimentalist would be able to make an extrapolated estimate of the zero-error value, presuming that the error sources had scaled proportionately. The technique was found to be very advantageous in the numerical simulations of few-qubit experiments presented in that paper (see e.g. Fig. 5 in Ref. <cit.>). 

A paper that appeared online at almost the same time was Ref. <cit.> by the IBM-based team of Temme, Bravyi and Gambetta. This paper presented a comprehensive analysis of the extrapolation technique, which the authors had independently conceived, and moreover it introduced a second technique with using (what we will call) a `quasi-probability' formalism. The authors explained that by replacing operations in the quantum circuit and assigning parity ±1 to each operation following a certain probability distribution dependent on the noise, an experimentalist can obtain the unbiased estimator, at the cost of an increase in the variance. Their method was shown to be applicable to specific noise types including homogeneous depolarizing errors and damping errors. The authors found both methods to be promising in few-qubit numerical simulations (see e.g. Fig. 2 in Ref. <cit.>).
 
As exciting as these studies were, open questions remained to be answered before these two techniques could be considered to be fully practical. First, both techniques rely on the full knowledge of the error model, whereas an experimentalist will have imperfect knowledge and the real noise will generally differ from the canonical types considered in these first papers. Second, we need an explicit method to derive the QEM circuits, i.e. a specification of how to algorithmically increase the error rate in the error extrapolation or how to sample circuits in the quasi-probability decomposition. In this paper, we solve these two problems. We find that gate set tomography (GST) <cit.> provides sufficient information to enable full elimination of the impact of localised Markovian errors. As with other process tomography protocols, GST cannot determine the exact physical error model due to noise associated with state preparation and measurement. However, we determine that preparation and measurement noise in GST is not harmful to the overall QEM approach. We also find that single-qubit Clifford gates and measurements are universal in computing expected values. Each quantum operation is a linear map, and single-qubit Clifford gates and measurements yield a complete set of linearly-independent maps (quantum operations). Therefore, any error can be simulated or subtracted by decomposing the error using the complete operation set, which is the standard linear decomposition. We prove that, by combining GST and the complete set decomposition, any localised and Markovian errors in the quantum computer can be systemically mitigated, so that the error in the final computational output is only due to unbiased statistical fluctuation.

For the quasi-probability method, we provide an upper bound of the cost in QEM, and we describe the utility of `twirling' operations <cit.> in minimising this cost. For the extrapolation method, which is a relatively straightforward technique, our optimisation is to observe that typically for the classes of noise most common in experiments it is appropriate to assume that the expected value of the observable will decay exponentially with the severity of the circuit noise. Adopting this underlying assumption, rather than a polynomial (e.g. linear) fit, proves to be quite advantageous. 

Having thus optimised both the quasi-probability and the extrapolation techniques, we make a series of numerical simulations to study their efficacy. We opt for a specific circuit, a realisation of the `SWAP test' that is often employed in quantum algorithms as a means for estimating the similarity of quantum states <cit.>. Our SWAP test operates on 2n+1 qubits, and we simulate a total of 15 qubits over a comprehensive set of cases as well as 19 qubits for specific cases. We numerically simulate the actions of the experimentalist, who must perform many circuit trials in order to make a single estimate of the observable (we choose 10^4 trials). But in order to evaluate our QEM techniques we must then repeat this entire process to determine the distribution of values that the experimentalist might obtain. We perform at least 10^3 repetitions so that the distribution becomes clear, thus at least 10^7 individual numerical experiments are performed for each of the curves that we presently report.



§ ERROR MITIGATION




In this paper, we focus on computing the expected value of an observable in a state (the final state of a quantum circuit) using a quantum computer. It is typical of a number of quantum algorithms and subroutines that the desired output is the expected value of a qubit or qubits – the SWAP-test <cit.> itself, which is a component of algorithms including the recently-introduced auto encoder <cit.>, and several proposed hybrid algorithms for simulating chemical or materials systems <cit.>.

Without using QEM as shown in Fig. <ref>(a), the quantum circuit is repeated for many times, and the measurement outcome μ of each time is collected. Then, we can calculate the average μ as our best estimate of the expected value. Given that the number of repetitions is finite, the value of μ is a random variable with an associated distribution. Because the implementation of the quantum circuit is imperfect, it is likely that the distribution of μ is not even centered at the ideal value, i.e. the exact expected value when the quantum circuit is perfectly implemented without error.

When we use QEM as shown in Fig. <ref>(b), instead of the original quantum circuit, we implement a set of modified circuits. The scheme depicted in the figure is relevant to the quasi-probability method for QEM, but can also apply to the extrapolation method as a means to deliberately boost errors. Each modified circuit is determined by a set of random numbers l. The distribution of random numbers, i.e. modified circuits, depends on the error model, which is measured using GST before the quantum computing. In each run of the quantum experiment, firstly the random number set l is generated, then depending on l a specific circuit is implemented, and finally the measurement outcome μ is collected. Rather than calculating the average μ, we use both l and μ to calculate the average of an effective outcome μ_ eff(l,μ), which will be given explicitly later. If QEM is successful, the distribution of μ_ eff(l,μ) is centered at the ideal value, but the distribution is wider than μ. Thus only error due to the statistical fluctuation remains, although it is amplified. By repeating the quantum experiment enough times, we can obtain an accurate computing result of the expected value.

In Sec. <ref>, we explicitly give the effective outcome μ_ eff(l,μ). Modified circuits and their distribution are given in Sec. <ref>.



§ NOTATION FOR STATES, OPERATORS AND OPERATIONS


We use the notation commonly used in quantum tomography (e.g. in Refs. <cit.>).

In quantum theory, a quantum state is usually represented by a density matrix ρ, and an observable is represented by a Hermitian operator Q. The expected value of the observable quantity in the state is Q = (Qρ). An operation is a map on the space of states, (ρ) = ∑_k E_k ρ E_k^†, expressed in the Kraus form.

Because an operation is a linear map, we can always express the operation  as a square matrix, e.g. using the Pauli transfer matrix representation, acting on the state expressed as a column vector ρ. Similarly, an observable can be expressed as a row vector Q, and the expected value is Q = Qρ. Throughout this paper, we use the Pauli transfer matrix representation, and see Appendix <ref> for details. In quantum tomography, usually we focus on observables that are POVM operators, which is not necessary here.

In the Pauli transfer matrix representation, vectors representing states or observables and matrices representing operations are all real. For n qubits, vectors and matrices are 4^n-dimensional. The expected value of the observable Q in the state ρ going through a sequence of operations _1, …, _N reads as follows:

    [Q_N ∘⋯∘_1(ρ)] = Q_N ⋯_1 ρ.




§ QUANTUM COMPUTING BY SAMPLING CIRCUITS



We suppose that the initial state is ρ^(0), which goes through a sequence of operations _1^(0), …, _N^(0), and in the final state the observable Q^(0) is measured. Each time the experimentalist implements this circuit, the measurement returns an eigenvalue of Q^(0), and the probability distribution of eigenstates is determined by the final state. By repeating such a circuit for many times, she can estimate the expected value Q^(0) = Q^(0)_ tot.^(0)ρ^(0) =  E[μ^(0)], where _ tot.^(0) = _N^(0)⋯_1^(0), and μ^(0) is the measurement outcome. Generally in this paper we will use the superscript 0 to denote the ideal noise-free realisation of a state, operation or observable quantity.

In the case that the quantum computation has errors, the actual initial state is ρ, actual operations are _1, …, _N, and the actually measured observable is Q. As a result, the estimation of the expected value converges to Q = Q_ tot.ρ rather than Q^(0). Here, _ tot. = _N ⋯_1, and we have assumed that errors are Markovian.

The central idea introduced by the IBM team in Ref. <cit.> is that one can exactly compensate for the effect of errors by sampling from a set of (real, error-burdened) circuits, each labelled _ tot.^(l) for l=1,2..., provided that their outputs satisfy 

    Q^(0) = ∑_l q_l Q^(l)_ tot.^(l)ρ^(l).
 
Ref. <cit.> describes how the real numbers { q_l } which represent quasi-probabilities can be efficiently derived given specific error models, assuming that the experimentalist has full knowledge of the model. Note that each _ tot.^(l) denotes the total operation composed by a sequence of operations in the l^th circuit.

We can use the Monte Carlo method to compute Q^(0). We note that Q^(l)_ tot.^(l)ρ^(l) =  E[μ^(l)], where μ^(l) is the measurement outcome in the l^th circuit. Then Q^(0) = ∑_l q_l E[ sgn(q_l)μ^(l)]. To compute Q^(0), we randomly choose a circuit to implement, and the l^th circuit is chosen with the probability p_l = q_l/C, where C = ∑_l q_l. Then, the computing result is given by the expected value of effective measurement outcomes, i.e. Q^(0) = C  E[μ_ eff], where the effective outcome is μ_ eff =  sgn(q_l)μ^(l) if the l^th circuit is chosen to be implemented, and μ^(l) is the outcome directly obtained in the l^th circuit.



§ PER-OPERATION ERROR CORRECTION



We can correct errors in each operation using the quasi-probability method, which will be the primary focus for the following several sections. We suppose that we have a set of initial states satisfying ρ^(0) = ∑_l_ in q_l_ in^[ in]ρ^(l_ in), and a set of operations satisfying _i^(0) = ∑_l_i q_l_i^[i]_i^(l_i) for each error-free operation _i^(0), and a set of observables satisfying Q^(0) = ∑_l_ out q_l_ out^[ out]Q^(l_ out). Then, computing with error mitigation can be expressed as

    Q^(0)   =   ∑_l_ in∑_l_1⋯∑_l_N∑_l_ out
    q_l_ in^[ in] q_l_1^[1]⋯ q_l_N^[N] q_l_ out^[ out]
          ×Q^(l_ out)_N^(l_N)⋯_1^(l_1)ρ^(l_ in).

When we sample circuits to compute Q^(0) = C  E[μ_ eff], the initial state is ρ^(l_ in) with  probability p_l_ in^[ in] = q_l_ in^[ in]/C_ in, the i^th operation is _i^(l_i) with  probability p_l_i^[i] = q_l_i^[i] / C_i, and the observable is Q^(l_ out) with  probability p_l_ out^[ out] = q_l_ out^[ out]/C_ out. Here, C_α = ∑_l_αq_l_α^[α], and C = C_ in C_1 ⋯ C_N C_ out accordingly. To calculate μ_ eff, we use sgn(q_l_ in^[ in]⋯ q_l_ out^[ out]) =  sgn(q_l_ in^[ in]) ⋯ sgn(q_l_ out^[ out]).



§ VARIANCE AMPLIFICATION IN QUASI-PROBABILITY DECOMPOSITION


The presence of quasi-probabilities taking negative values amplifies the variance of the expected value of the observable. We consider the case that Q^(l) is a Pauli operator (maybe with error) and the measurement reports two kinds of outcomes denoted by ± 1, respectively. In this case, the distribution is binomial. The standard deviation of the average of outcomes in the Monte Carlo calculation is σ = C √((1- E[μ_ eff]^2)/N_ r)≤ C/√(N_ r). Here, N_ r is the total number of samples, i.e. the total number of circuits of all kinds which the experimentalist performs is N_ r. We compare this to the error-free computing, i.e. the ideal original circuit Q^(0)_ tot.^(0)ρ^(0) is repeated for N_ r^(0) times to estimate Q^(0). For the error-free computing, the standard deviation is given by σ^(0) = √((1- E[μ^(0)]^2)/N_ r^(0)). Therefore, to achieve the same accuracy, i.e. σ = σ^(0), the error-mitigated computation needs N_ r/N_ r^(0) = (C^2- Q^(0)^2)/(1-Q^(0)^2) times more samples than the error-free computation. Here, we have used the fact that the error-mitigated computation and the error-free computation should converge to the same value of Q^(0), i.e. E[μ^(0)] = C  E[μ_ eff].

In order to limit the standard deviation to be σ∼ϵ, we can choose N_ r∼ (C/ϵ)^2. Therefore, if the factor C is larger, the computing takes longer.

Because C = C_ in C_1 ⋯ C_N C_ out if errors are corrected for each operation, we call C_α - 1 the cost for mitigating error in the corresponding operation. The overall cost therefore increases with the number of operations, thus it is important to reduce the operation number, e.g. in quantum computer with qubits fully connected <cit.>, operations for communication are not required, which may significantly reduce the cost.



§ UNIVERSAL OPERATION SET




The set of operations including measurement and single-qubit Clifford gates is universal in computing expected values of observables. The relevant measurement operation reads [π] = [1/2( + σ^ z)], which projects a qubit to the state |0⟩ . Here, [U](ρ) = Uρ U^† denotes a superoperator. Such a non-destructive measurement can be realised using a destructive measurement followed by initialising the qubit in the state |0⟩. Single-qubit Clifford gates include the Hadamard gate [H] = [1/√(2)(σ^ x + σ^ z)], the phase gate [S] = [1/√(2)( -i σ^ z)] and all other single-qubit Clifford gates can be derived from these two.

The measurement superoperator [π] also means post-selection, i.e. if the outcome of the measurement corresponding to [π] (which is not the final measurement on the observable Q^(l)) is |1⟩ in a trial, the value of the observable Q^(l) is noted as μ^(l) = 0, but the trial is counted in the total number of samples in the Monte Carlo calculation. If Q^(l) has two values ± 1, we can estimate the value of Q^(l)_ tot.^(l)ρ^(l) by calculating (N_+1^(l) - N_-1^(l))/(N_0^(l) + N_+1^(l) + N_-1^(l)). Here, we have supposed that the circuit is implemented for total N_0^(l) + N_+1^(l) + N_-1^(l) times; for N_0^(l) times the circuit does not pass post-selections (i.e. μ^(l) = 0), and for N_± 1^(l) times the circuit passes all post-selections and reports Q^(l) = ± 1 (i.e. μ^(l) = ± 1). It is the same when we compute Q^(0) using the Monte Carlo method. If the effect outcome is μ_ eff = 0 with the probability P_0, then the standard deviation of the Monte Carlo calculation becomes σ = C √([(1-P_0)^2- E[μ]^2]/(1-P_0)N_ r)≤ C√((1-P_0)/N_ r).

In Table <ref>, we list sixteen linearly independent operations that can be derived from the minimum universal operation set { [π], [H], [S] }. In the following, we use {_i^(0)| i=1,…,16} to denote these sixteen operations. Because they are linearly independent, any single-qubit operation , which is a 4× 4 real matrix, can be expressed as a linear combination of sixteen basis operations, i.e. = ∑_i=1^16 q_i _i^(0). Similarly, multi-qubit operations can be expressed as a linear combination of tensor products of basis operations. Using the quasi-probability method, any computation of expected values of observables can be realised using this operation set.



Note that these basis operations are universal, as one can verify by constructing a non-Clifford gate or an entangling gate: We can decompose T gate using our basis operations as [T] = 1/2[] - √(2)-1/2[σ^ z] + √(2)/2[R_ z^3] (see Appendix <ref> for controlled-NOT as a second example). However this construction would not be used in practice – it is not an efficient means to actually implement a desired T in the basic circuit, since the corresponding cost C = √(2) would imply and unacceptably steep exponential in the time overhead, as one would expect from e.g. Refs. <cit.>. Instead we rely on the assumption that the experimental system can directly implement a universal set of gates (including entangling and non-Clifford gates) with a reasonably high fidelity. Then rather than fully synthesising any of the basic gates using our basis, we need only compensate for slight imperfections. The cost for doing so, for each imperfect gate, is then ∼ C = 1+δ as we presently discuss.

Having obtained the complete operation set in Table <ref> we can use it in deriving the protocol that will compensate for errors. In this paper, we focus on the case that errors are localised: An (error-free) operation that is applied on a set of qubits S is a 4^S-dimensional real matrix, then the corresponding operation in real (i.e. error-burdened)  is also a 4^S-dimensional real matrix acting on the same set of qubits. The overall operation on the entire system can be expressed as _S̅⊗, where _S̅ is the identity acting on all other qubits. It is similar for the initialisation and measurement. If each qubit is initialised individually, the overall initial state is ⊗_m ρ_m, where ρ_m is a 2-dimensional real vector representing the m^th qubit's initial state. Similarly, individual measurement of qubits implies that the overall measured observable is ⊗_m Q_m, where Q_m is a 2-dimensional real vector representing the measured observable for the m^th qubit. In this case, a single-qubit operation with error can still be expressed using a 4× 4 real matrix. We suppose that for a qubit, sixteen basis operations with errors are {_i | i=1,…,16}, which are all 4× 4 real matrices. When errors are not significant, these sixteen bases should still be linearly independent, i.e. the set of basis operations with errors is still universal.

To make this statement more precise, we consider the 16× 16 real matrix

    A = [
    [  (_1)_∙,1         ⋯ (_16)_∙,1;  (_1)_∙,2         ⋯ (_16)_∙,2;  (_1)_∙,3         ⋯ (_16)_∙,3;  (_1)_∙,4         ⋯ (_16)_∙,4 ]].

Here, (_i)_∙,j denotes the j^th column of the matrix of the basis operation _i. Sixteen basis operations are linearly independent if the matrix A is invertible. We use ϵ_ max = max{_i - _i^(0)_ max| i=1,…,16 } as the measure of the error severity in basis operations. When ϵ_ max < 1/32(13-3√(17)) ≃ 0.0351, A is always invertible (see Appendix <ref>). We remark that even if ϵ_ max exceeds the threshold, basis operations are still likely to be linearly independent.



§ ERROR MITIGATION USING BASIS OPERATIONS





Given an operation with error , we can use sixteen basis operations to correct the error, i.e. realise the operation without error ^(0). There are two ways for correcting the error.

Compensation method. The operation  is close to ^(0). Therefore, we can keep the correct component of  and only decompose the error component using basis operations. We decompose the operation without error as ^(0) = λ + ∑_i q_i _i, where λ is an arbitrary real number. If basis operations are linearly independent, the decomposition always exists, and there is only one solution of coefficients {q_i} when λ is determined.

Inverse method. If the matrix ^(0) is invertible, we can express  as ^(0) followed by a noise operation, i.e. = ^(0), where the noise operation = ^(0)-1. In order to correct the error, we can decompose the inverse of the noise as ^-1 = ^(0)^-1 = ∑_i q_i _i. By applying the inverse of the noise after the operation , we can realise the operation without error, i.e. ^(0) = ^-1 = ∑_i q_i _i. Similar to the compensation method, if basis operations are linearly independent, the decomposition always exists, and there is only one solution of coefficients {q_i}. However, the inverse method can only be applied if the matrix  is invertible.

For multi-qubit operations, the decomposition is performed using tensor products of basis operations, as described explicitly in Appendix <ref>. Although basis operations are not entangling, we can use basis operations to efficiently mitigate multi-qubit errors and errors that can entangle qubits. As an example, we show how to decompose the controlled-NOT gate only using basis operations in Appendix <ref>, which suffices to imply that any error in the form of the controlled-NOT gate can be mitigated using basis operations.

Initialisation and measurement errors can also be corrected using basis operations. Taking first the case of initialisation errors: If ρ is the error-burdened initial state, and it is a non-zero vector, we can always find a transformation  that satisfies ρ^(0) = ρ where ρ^(0) is the error-free initial state. Thus by decomposing  using basis operations and applying it after the initialisation, we can prepare the initial state without error. Actually, given an initial state that is close to |0⟩, we can generate a complete set of linearly independent vectors {ρ_k} using basis operations. With these vectors, we can decompose the initial state  without error as ρ^(0) = ∑_k q_k ρ_k.

A similar approach yields the corresponding result for measurement: For an observable Q there will be some Q^(0) = Q where Q^(0) is the error-free quantity. If an observable is close to σ^ z then a linearly independent set {Q_j} can be generated; then the error-free observable Q^(0) = ∑_j q_j Q_j.

Circuits for QEM are shown in Fig. <ref>(a). Given quasi-probabilities, we can compute the corresponding probability in sampling circuits as shown in Sec. <ref>. More details of QEM using basis operations are given in Appendix <ref>.

Using the same technique, we can also increase the error in an operation, as required by the alternative error extrapolation method for QEM. Instead of decomposing the error-free operation ^(0) using  and basis operations, we can also decompose the error-boosted operation _ b(r) = (1-r)^(0) +r (r>1) using  and basis operations. It is similar for initial states and observables. We have noted that in the decomposition of an error-free operation, there are always some negative quasi-probabilities, i.e. the C factor is greater than 1, which leads to greater time costs. But fortunately when we merely wish to decompose an error-boosted operation we can do so without introducing negative quasi-probability, e.g. by boosting Pauli errors using Pauli gates <cit.>.



§ QUANTUM GATE SET TOMOGRAPHY


We can measure a set of initial states {ρ̅_k}, observables {Q̅_j} and operations {_i } (including basis operations) using GST <cit.>. These vectors and matrices with the bar notation describe the actual physical system. Because there are errors in both initial states and observables, and initialisation and measurement errors cannot be distinguished, we may not obtain exactly these vectors and matrices describing the actual physical system. Instead, the vectors and matrices obtained using GST are {ρ̂_k}, {Q̂_j} and {_i }, which are estimations of {ρ̅_k}, {Q̅_j} and {_i }, respectively.

If we know {ρ̅_k}, {Q̅_j} and {_i } because the physical system is well understood, we can directly use them in QEM. If our knowledge about the physical system is not enough, we can use GST to obtain {ρ̂_k}, {Q̂_j} and {_i }. We will show that, although the estimations may not be exact, we can exactly correct errors by using these estimations in QEM.

Using the protocol in Refs. <cit.> (also see Appendix <ref>), the estimation of an operation and the actual physical operation are similar matrices, i.e. _i = TM̅^ in -1_i M̅^ in T^-1, where M̅^ in is a matrix determined by initial states (i.e. M̅^ in_σ, k = σρ̅_k), and T is an arbitrary invertible matrix. We note that T and M̅^ in are independent of the operation _i, and M̅^ in cannot be determined by GST. By choosing T, we can obtain different estimations of the operation set. Similarly, ρ̂_k = TM̅^ in -1ρ̅_k and Q̂_j = Q̅_jM̅^ in T^-1.

All operations are transformed by the same similarity transformation, and initial states and observables are also transformed accordingly. As a result, these estimations obtained by GST can exactly predict the expected value of an observable, i.e. Q̅_j_N ⋯_1 ρ̅_k = Q̂_j_N ⋯_1 ρ̂_k. Therefore, we can directly use these estimations in QEM, and the similarity transformation does not lead to any computing error.

Using GST estimations in QEM, the actual operations realised in this way differ from operations without error, but the computing result is correct. To correctly obtain Q^(0)^(0)ρ^(0), we decompose the initial state, the observable and the operation using {ρ̂_k}, {Q̂_j} and {_i } respectively. Here, ρ^(0), Q^(0), ^(0) and GST estimations are all known to us. The decompositions are ρ^(0) = ∑_k q_k ρ̂_k, Q^(0) = ∑_j q_j Q̂_j and ^(0) = ∑_i q_i _i. Accordingly, we actually realise ρ̅^(0) = ∑_k q_k ρ̅_k, Q̅^(0) = ∑_j q_j Q̅_j and ^(0) = ∑_i q_i _i in the physical system. We have ρ̅^(0) = M̅^ inT^-1ρ^(0), Q̅^(0) = Q^(0) TM̅^ in -1 and ^(0) = M̅^ in T^-1^(0) TM̅^ in -1. Therefore, the physical system gives the computing result Q̅^(0)^(0)ρ̅^(0) = Q^(0)^(0)ρ^(0), i.e. the desired error-free output. The cost of this adaption lies in the potential increase to the number of samples required, as shown in Fig. <ref> and discussed in the caption.

We would like to remark that, when errors in actual operations are small, errors in estimations of operations are also small. If we take a proper strategy for choosing T, and errors in initial states and observables are small, the estimation of an operation  is close to the operation without error ^(0) when the actual operation  is close to ^(0). It is similar for estimations of initial states and observables. See Appendix <ref> for details.



§ ESTIMATION OF THE COST




In general, when the error in an operation is more significant, there is a higher cost for mitigating the error (to a given level of suppression). We take ϵ_ =  - ^(0)_ max as the measure of the error severity in the operation, where  (^(0)) is the n-qubit operation with (without) error. An upper bound of the cost for correcting error in  is

    C_ - 1 ≤16^2nϵ_/[s_ min(A^(0)) - 16ϵ_ max]^n.

Here, ϵ_ max is the maximum error in all basis operations for all n qubits, and s_ min(A^(0)) = 1/2(13-3√(17)) ≈ 0.315. Similar upper bounds can be obtained for correcting errors in initial states and observables. See Appendix <ref> for details.

There are several ways for reducing the cost. The upper bound of the cost is obtained using the compensation method and taking λ = 1. In general, we can optimise the value of λ or use the inverse method to minimise the cost. For example, for the depolarising error model (see Appendix <ref>), the cost of using the inverse method is lower than using the compensation method [see Fig. <ref>(a)]. We remark that, to obtain data for the compensation method in Fig. <ref>, we have optimised the value of λ. If we use estimations obtained from GST to correct errors, we can optimise the T matrix to minimise the cost. In Fig. <ref>(a), we can find that, without optimising T matrices, the cost using estimations obtained from GST is higher than using actual operations. If we choose the matrix in the form T = ⊗_m=1^n T_m, where T_m is a 4-dimensional real matrix corresponding to the m^th qubit, there are total 16n parameters to be optimised for a n-qubit quantum computer, which is a non-trivial task. Under some reasonable conditions, we can also use the Pauli twirling <cit.> to reduce the cost.



 §.§ Pauli twirling


In many quantum computing systems, e.g. superconducting qubits <cit.> and ion traps <cit.>, the fidelity of single-qubit gates is much better than the fidelity of two-qubit gates, and usually a state can be initialised with a high fidelity while the fidelity of measurement is worse. In this section, we focus on the case that error rates of initialisation and single-qubit gates are much lower than error rates of two-qubit gates and measurement.

If the error rate of initialisation is low (much lower than the error rate of measurement), we know how to choose T so that the estimation of an operation obtained from GST is close to the actual operation. We cannot exactly estimate operations using GST, because we cannot distinguish initialisation and measurement errors. If we treat all errors in the initialisation and measurement as measurement error (which corresponds to T = M^ in(0)⊗ n in Appendix <ref>), the difference between the estimation and the actual operation is only determined by the initialisation error. Therefore, if the initialisation is high-fidelity, the estimation obtained in this way and the actual operation are close.

Because the set of basis operations includes Pauli gates, it is easy to use basis operations to correct Pauli errors. By using the Pauli twirling, we can convert the error in a two-qubit entangling Clifford gate to Pauli error <cit.>, which is achieved by applying Pauli gates before and after the two-qubit gate. This treatment of the error is feasible only if the fidelity of Pauli gates is much better than the two-qubit gate, otherwise Pauli gates cause significant new errors, which may not be Pauli error, on the two-qubit gate. In Fig. <ref>(b), we can find that the cost can be significantly reduce by using the Pauli twirling for the over-rotation error model (see Appendix <ref>).

In Fig. <ref>(c), costs of different error models are compared, including the depolarising model, pure-dephasing model, amplitude-damping model and the over-rotation model. We also randomly generated many other error models, please see Appendix <ref> for details of these error models. For a random-operation model, we randomly generate an operation close to the ideal error-free operation, and we find that the cost is approximately the cost of the depolarising model. For a random-field model, we randomly generate a Hamiltonian that drives the erroneous evolution, and the cost is between the depolarising model and over-rotation model.

From Fig. <ref>(c) we see that the cost of quantum error mitigation varies according to the error model but is generally upper-bounded by the case of depolarising noise, over the range of noise levels shown here. (Note that other models can exceed the cost of the the depolarising model if we use even lower fidelity gates). For the depolarising model, the cost for mitigating error in a two-qubit entangling gate is C-1 ≃ aϵ, where ϵ is the error rate and the factor a is between 2 and 3 [see Fig. <ref>(a)]. If errors in initialisation and single-qubit gates are negligible, or if the matrix T is optimised to minimise the cost, the factor a can approach 2. Accepting the depolarising model as an approximate upper bound, we can estimate the overall cost in a quantum algorithm. Suppose the total number of gates in a quantum algorithm is N, the overall amplification of the standard deviation (uncertainty of the computing result) is (1+2ϵ)^N. Therefore, (1+2ϵ)^2N times more repetitions of the experiment are required in order to reduce the standard deviation. We are interested in the case that N is large but ϵ is small, therefore, (1+2ϵ)^2N∼ e^4Nϵ. As a rule of thumb we might take Nϵ = 2 as a limit for acceptable scenarios, since then e^4Nϵ≈ 3,000. However, larger overhead factors may be acceptable depending on the speed of the quantum computer.



§ NUMERICAL SIMULATION




In our numerical simulation, we apply QEM to the SWAP-test circuit <cit.> shown in Fig. <ref>, in which we realise each controlled-SWAP gate using Toffoli gates and realise each Toffoli gate using T gates, T^† gates, Hadamard gates and controlled-NOT gates <cit.>. We note with interest that very recently, the implementation of a SWAP test using shallow circuit has been proposed <cit.>. However, for present purposes it is not essential to use an optimised realisation of the SWAP circuit; its role is simply to act as a real test case for our technique and indeed the considerable depth of our non-optimal circuit is helpful here. The number of gates scales as 23 N_ q - 21, where N_ q is the number of qubits (e.g. N_ q= 7 in Fig. <ref>). Without error, the expected value of the observable Z (σ^ z of the probe qubit) in the SWAP-test circuit in Fig. <ref> is 0.5.



We consider error models according to which the same noise ℰ is applied after the initialisation to the state |0⟩, before the measurement, and before and after each gate. For the controlled-NOT gate, the noise applied is ℰ⊗ℰ on two qubits. We remark that basis operations are also affected by noise likewise. We consider two types of noise: inhomogeneous Pauli error and leakage error, which can be respectively described as

    ℰ_ inh = (1-p_ x-p_ y-p_ z)[] + p_ x[σ^ x] + p_ y[σ^ y] + p_ z[σ^ z]

and

    ℰ_ leak = [|0⟩⟨0| + √(1-p)|1⟩⟨1|],

where p_α is the probability of the error [σ^α], and p is the probability of the leakage error from the state |1⟩. It is worth mentioning that the leakage error is a non-trace-preserving error. In our simulations, we set p_ x = p_ y = 0.0001, p_ z = 0.0006, and p = 0.0008. Thus in both models the total error rate is 0.08% for initialisation and measurement, 0.16% for single-qubit gates and 0.32% for two-qubit gates, which is achievable with two-qubit gates in ion traps <cit.> and can be far surpassed for one-qubit gates <cit.>. Moreover, with these numbers the expected total number of error events in circuits of the depth and breadth that we consider here is approximately unity; this is a challenging domain for error mitigation.

In addition to quasi-probability decomposition (see Appendix <ref> for an instruction of the implementation), we also study the extrapolation technique introduced in Ref. <cit.>. The expected value of Z obtained by running the SWAP-test circuit in a quantum computer with noise depends on the error rate, i.e. it is a function that can be denoted as Z(ϵ), where ϵ is the overall error rate. For our first set of numerical experiments we consider linear extrapolation to the error-free value Z(0) as follows: We obtain the expected value Z(ϵ_0) with the lowest attainable error rate ϵ_0, and by increasing error rate to rϵ_0 with r>1, we obtain another expected value Z(rϵ_0). Using these two values, we can infer Z(0) = (rZ(ϵ_0) - Z(rϵ_0))/(r-1) as shown in Fig. <ref>(b), which is the final estimation of Z. Here, we set r = 2.

The first set of numerical results are shown in Fig. <ref>. We assume that the experimentalist makes her overall estimate of the Z after she performs 10^4 individual experiments. We take this number of runs as a fixed constraint (effectively, we are constraining her overall time resource), and she may choose to employ those runs using one of three alternative approaches: no error correction, linear extrapolation, and quasi-probability decomposition (using basis operations and incorporating GST). In each experiment the SWAP-test circuit or its variant for the purpose of QEM is implemented. Because of the finite number of samples, the estimation is stochastic. Therefore, in our numerical simulation we perform the appropriate series of 10^4 experiments, mirroring the actions of the experimentalist, and then we repeat ≥ 1,000 times in order to determine the distribution of final estimations that may be obtained. The distribution for each case is plotted in Fig. <ref>.



We can observe that both QEM approaches can improve the result, i.e. the corresponding distributions are shifted closer to the ideal value 0.5 compared to the approach without QEM. For the inhomogeneous Pauli error model, the means of distributions are at 0.1961, 0.3415, and 0.5011 for the three approaches, respectively. The distribution of the quasi-probability approach is centered at the ideal value, which clearly shows its desirable property of completely removing any systematic bias. However, the distribution is wider (as we expected) compared to other two approaches. A more fair metric would be the expected absolute error versus ideal value (i.e. Z-0.5). Given an ideal error-free computer and 10^4 trials, this metric would evaluate to 0.006910. Using the error-prone computer with our three protocols the three corresponding values are 0.3039, 0.1853 and 0.0491. Similarly, for the leakage error model, the means for three approaches now lie at 0.3819, 0.4710, 0.5007, while the expected absolute error evaluates to  0.1181, 0.0294, and 0.0434.

From these results it may appear that (given a large but reasonable number of samples) the quasi-probability technique outperforms the extrapolation method, with the latter unable to approach the mean of the error-free circuit. However, here the extrapolation method was limited to linear interpolation whereas the physical error rates are high enough that the linear assumption is poor. One could fit a higher order polynomial using more data points (here, we have only used two: one derived from the actual lowest possible error rate and one boosted to twice the error rate); however since we are limiting the total number of experimental runs to 10^4 this would lead to greater noise in each data point. Moreover, as we now argue, the underlying tend is likely to be well-approximated by an exponential decay rather than a polynomial one (i.e. the expected value of the observable falls exponentially with the physical error rate) and two data points will suffice to estimate the zero-error observable under that assumption.

In Fig. <ref> we show the results  when the experimentalist indeed assumes that the expected value Z(ϵ) changes exponentially with respect to the error rate ϵ and converges to 0 in the limit of ϵ→∞. Then she will infer the error-free value as

    Z(0) = Z(ϵ_0)^r/r-1Z(rϵ_0)^1/1-r

Here we take r=2.

As shown in Fig. <ref>, the distribution of the final result using the exponential extrapolation approaches the ideal value of Z (which is 0.5 for the SWAP-test circuit) much better than the linear extrapolation. Given the same 10^4 experimental runs, the mean of the experimentalist's estimate is now 0.5111 for the inhomogeneous Pauli error model and 0.4986 for the leakage error model. These numbers almost rival those of quasi-probability technique but do so with a smaller variance. The expected absolute error for inhomogeneous Pauli error and leakage error are 0.06501 and 0.01882, respectively. For the latter, the expected absolute error comes within a factor of three of the shot-noise limit that would be achieved by error-free ideal hardware (0.00691). This is despite the fact that our error-burdened circuits have error rates corresponding to at least one error event per circuit. We emphasise that this suppression results purely from the QEM protocol i.e. it is achieved at no cost in terms of the qubit count or the total number of runs (constrained to 10^4).

Due to the limited power of classical computer we utilised, our exact numerical simulations did not involve more than 19 qubits. However, it is of course very interesting to assess the relevance of our techniques to quantum computing using over 50 qubits, which is in the so called `quantum supremacy' regime. Therefore, we estimate the cost of quantum error mitigation in the SWAP-test circuit, using the same error models in our numerical simulation and error rates achievable in ion trap experiments <cit.>, i.e. the error rate of two-qubit gate is 0.1% and error rates of single-qubit operations are 0.01%. Take for example the SWAP test with N_ q = 51 qubits (the number of gates is 1,152). For the inhomogeneous Pauli error model, the overall cost is C = 2.956, which implies that we can attain the same computing precision as the ideal case if we have C^2 = 8.738 times more repetitions of the experiment, which is experimentally feasible. For the leakage error model, the cost for the 51-qubit SWAP test is C = 4.338, which means C^2 = 18.818 times more repetitions. A plot showing how the cost scales versus qubit count is shown in Fig. <ref>.

We also evaluate C^2 for a fully paralleled circuit, whose circuit depth is N_ q, and each layer has N_ q/2 single qubit gates and N_ q/4 controlled-NOT gates, which means the quantum circuit has N_ q^2/2 single qubit gates and  N_ q^2/4 controlled-NOT gates. As single qubit gates, we use T gate, S gate and Hadamard gate, because these gates plus controlled-NOT gate constitute a universal gate set, and we equally assign the number of qubits to these three types of single qubit gates. We plot C^2 versus the number of qubits for the SWAP-test circuit and the fully paralleled circuit in Fig. <ref>. We observe that for the SWAP test circuit, it is feasible to venture into the `supremacy' regime with today's best fidelities; for the more demanding case of full parallelism (so that the gate count scales as N_ q^2) we see that today's error rates would not suffice much beyond 50 qubits, but that error rates ten times lower would easily suffice for 80 qubits and beyond.





§ INTUITION FOR EXPONENTIAL EXTRAPOLATION



Intuitively, the explanation for the success of the exponential extrapolation is as follows. We express the i^th noise event occurring in the quantum circuit as

    ℰ_i(ϵ) = (1-ϵ) [I] + ϵℰ'(ϵ)_i,

where ℰ'(ϵ) is the error component. The only assumption is that the error component only weakly depends the error rate ϵ (see Appendix <ref>). Now, for simplification we ignore the computing operations, which do not affect our general argument. The total noise that the entire quantum circuit experiences is

    ∏_i=1^N ℰ_i = ∏_i=1^N [(1-ϵ) [I] + ϵℰ'_i].

Here, N is the total number of the noise-burdened operations. Expanding the overall noise, we get

    ∏_i=1^N ℰ_i = ∑_n=1^N Nn (1-ϵ)^N-nϵ ^n 𝒳_n,

where 

    𝒳_n = Nn^-1×(
    [    the sum of terms where; ℰ'  appears for  n  times ]).

Note that the coefficient of 𝒳_n in the overall noise corresponds to a binomial distribution, which can be approximated by the Poisson distribution. We have 

    ∏_i=1^N ℰ_i = e^-Nϵ∑_n=0^N (Nϵ)^n/n!𝒳_n,

We can find that the impact of the overall noise on the expected value of some observable is proportional to e^-Nϵ, which implies that exponential extrapolation works better than linear extrapolation.



§ CONCLUSIONS


We have demonstrated that, following our protocol step by step, an experimentalist can derive an algorithm to run on a noisy quantum computer so as to estimate an output observable with zero bias versus the ideal observable. The experimentalist does not require any prior knowledge of the physical property of the noise, and the only condition is that the noise is localised and Markovian. For this purpose, we have shown that quantum gate set tomography is a perfect tool for measuring the noise in a quantum computer, if the aim is only to compensate the effect of the noise in quantum computing; and we also have shown that single-qubit Clifford gates and measurement can derive a complete set of operations that can compensate any noise in quantum computing. 

The price of using such a systematic method to negate computing errors is that the quantum computation needs to run for a longer time than an error-free system. We verify the protocol with numerical simulations of up to 19 qubits, in which an alternative method, i.e. exponential error extrapolation, is introduced and studied. We find that the estimation using exponential error extrapolation is also very accurate, while the computing time could be shorter. An approach combining two methods may optimise both accuracy and efficiency. 

In Appendix <ref> we describe in detail the steps that an experimentalist would take in order to realise the quasi-probability method. We hope that this compact summary, presented in a single section, will indeed be useful to researchers who are interested in demonstrating the QEM technique with their hardware. 


Our general conclusion is that these quantum error mitigation techniques can dramatically enhance the performance of quantum computers, especially at the small-to-medium scale where full code-based quantum error correction is impossible. Our simulations have considered circuits up to 19 qubits, but with error rates considerably worse than the state of the art. Extrapolating from the trends that we observe in these smaller systems, we anticipate that hybrid algorithms involving 50+ qubits, i.e. beyond the reach of classical emulation, will benefit from QEM techniques if the hardware fidelity matches today's state-of-the-art error or modestly improves upon it.


This work was supported by the EPSRC National Quantum Technology Hub in Networked Quantum Information Technologies. SE is supported by Japan Student Services Organization (JASSO) Student Exchange Support Program (Graduate Scholarship for Degree Seeking Students). YL is also supported by NSAF (Grant No. U1730449). The authors would like to acknowledge the use of the University of Oxford Advanced Research Computing (ARC) facility in carrying out this work. http://dx.doi.org/10.5281/zenodo.22558.






§ PAULI TRANSFER MATRIX



A state ρ can be expressed as a real column vector

    ρ = [
    [   ⋯ ρ_σ   ⋯ ]]^ T,

where the vector element is

    ρ_σ = (σρ),

σ∈{, σ^ x, σ^ y, σ^ z}^⊗ n is a Pauli operator, and d = 2^n is the dimension of the Hilbert space. Similarly, an observable (i.e. Hermitian operator) Q can be expressed as a real row vector

    Q = [
    [   ⋯ Q_σ   ⋯ ]],

where the vector element is

    Q_σ = d^-1(σ Q).

Here, we use notations · and · to denote real row and column vectors, respectively. A physical operation  (i.e. (ρ) = ∑_k E_k ρ E_k^†) can be expressed as a real square matrix

    _σ, τ = d^-1[σ(τ)],

where σ, τ∈{, σ^ x, σ^ y, σ^ z}^⊗ n are Pauli operators. If ρ' =  (ρ), we have ρ' = ρ.



§ DECOMPOSITION OF CONTROLLED-NOT GATE USING BASIS OPERATIONS



The controlled-NOT gate reads

    Λ_ X =  + σ^ z/2⊗ +  - σ^ z/2⊗σ^ x.

The controlled-NOT gate can be decomposed as

    [Λ_ X]    =   1/2( [⊗σ^ x] + [σ^ z⊗] - [⊗ R_ x] 
           - [R_ z⊗] - [σ^ z⊗ R_ x] - [R_ z⊗σ^ x]) 
           + [σ^ z⊗σ^ x] + [R_ z⊗ R_ x] + [⊗π_ x] 
           + [π_ z⊗] - [σ^ z⊗π_ x] - [π_ z⊗σ^ x].

Then, the corresponding cost is given by C = 9.



§ ERROR THRESHOLD OF BASIS OPERATIONS



For two real matrices A^(0) and A and a non-zero real vector x, we have

    A^(0)x_2 = √(x^ T A^(0)  T A^(0) x)≥ s_ min(A^(0)) x_2,

where s_ min(A^(0)) is the minimum singular value of A^(0). We also have

    (A-A^(0))x_2 ≤A-A^(0)_2 x_2.

Therefore,

    Ax_2    ≥   A^(0)x_2 - (A-A^(0))x_2 
       ≥    (s_ min(A^(0)) - A-A^(0)_2) x_2.

If A-A^(0)_2 < s_ min(A^(0)), Ax_2 is always positive (non-zero), i.e. A is invertible.

Now, A is the matrix formed by basis operations with error as defined in Eq. (<ref>), and A^(0) is the matrix formed by basis operations without error. Because (A^(0)) = 16, A^(0) is invertible, i.e. basis operations without error are linearly independent. The minimum singular value is s_ min(A^(0)) = 1/2(13-3√(17)). Because A-A^(0)_2 ≤ 16 Ã - A_ max = 16ϵ_ max, the matrix A is invertible if ϵ_ max < 1/16s_ min(A).



§ DECOMPOSITION USING BASIS OPERATIONS



We consider the n-qubit operation . For each qubit, there is a set of basis operations {_m,i| i=1,…,16}, where m=1,…,n is the label of the qubit. For each set of basis operations, there is a matrix A as defined in Eq. (<ref>). We use A_m to denote the matrix of the m^th qubit.

The operation  is decomposed as

    = ∑_i_1 = 1^16⋯∑_i_n = 1^16 q_i_1,…,i_n_1,i_1⊗⋯⊗_n,i_n.

Coefficients form a 16^n-dimensional vector

    q = (
    [     q_1,1,⋯ ,1,1;                ⋮;    q_1,1,⋯ ,1,16;                ⋮;  q_16,16,⋯ ,16,1;                ⋮; q_16,16,⋯ ,16,16 ]),

Therefore, the decomposition is given by q = (A_1 ⊗⋯⊗ A_n)^-1 E, where E is a 16^n-dimensional vector corresponding to .

We choose the order of Pauli operators, i.e. the order of bases of Pauli transfer matrices {_m,i| i=1,…,16}, as , σ^ x, σ^ y and σ^ z (which are also denoted as I, X, Y and Z, respectively). Then, to be consistent with A_1 ⊗⋯⊗ A_n, we have

    E = (
    [ _I_1I_2⋯ I_n-1I_n, I_1I_2⋯ I_n-1I_n; _I_1I_2⋯ I_n-1I_n, I_1I_2⋯ I_n-1X_n;                                   ⋮; _I_1I_2⋯ I_n-1Z_n, I_1I_2⋯ I_n-1Y_n; _I_1I_2⋯ I_n-1Z_n, I_1I_2⋯ I_n-1Z_n;                                   ⋮; _Z_1Z_2⋯ Z_n-1I_n, Z_1Z_2⋯ Z_n-1I_n; _Z_1Z_2⋯ Z_n-1I_n, Z_1Z_2⋯ Z_n-1X_n;                                   ⋮; _Z_1Z_2⋯ Z_n-1Z_n, Z_1Z_2⋯ Z_n-1Y_n; _Z_1Z_2⋯ Z_n-1Z_n, Z_1Z_2⋯ Z_n-1Z_n ]).

Here, α_m (α = I,X,Y,Z) is a Pauli operator of the m^th qubit.

The state of a qubit is represented by a 4-dimensional real vector. To decompose the initial state of a qubit without error ρ^(0), we need four linearly independent initial states. If the qubit can be initialised in the state |0⟩, we can choose the set of four states as {ρ_k^(0)} = {|0⟩, |1⟩, 1/√(2)(|0⟩+|1⟩), 1/√(2)(|0⟩+i|1⟩) }. These four states can be obtained by applying basis-adjusting operations (Clifford gates) { [], [R_ x], [R_ x]^2, [R_ z][R_ x] } on the initial state |0⟩. Because of the error in the state |0⟩ and errors in basis-adjusting operations, the prepared four states {ρ_k } are not exactly states {ρ_k^(0)}. When the overall error is small, states {ρ_k } are still linearly independent. We introduce the matrix M^ in_σ, k = σρ_k, and M^ in (0) is the matrix corresponding to {ρ_k^(0)}. States {ρ_k } are linearly independent if M^ in is invertible. Similar to the analyse of the linear independence of basis operations (i.e. the invertibility of the matrix A, see Appendix <ref>), we have that M^ in is always invertible if M^ in - M^ in (0)_ max < 1/4s_ min(M^ in (0)) = 1/8√(5-√(17)/2)≃ 0.0828. The initial state without error is decomposed as ρ^(0) = ∑_k=1^4 q_k ρ_k. Coefficients form a 4-dimensional column vector q = [q_1   q_2   q_3   q_4]^ T. The decomposition is given by q = M^ in -1ρ^(0).

Similarly, an observable of a qubit is also represented by a 4-dimensional real vector. To decompose the observable of a qubit without error Q^(0), we need four linearly independent observables. If σ^ z can be measured, we can choose the set of four observables as Pauli operators { Q_j^(0)} = {, σ^ x, σ^ y, σ^ z}. The operator  denotes a trivial measurement, i.e. the outcome is always +1. Measurements of other three Pauli operators can be obtained by applying basis-adjusting operations (Clifford gates) { [], [R_ x], [R_ z]^3[R_ x][R_ z] } before the measurement of σ^ z. Because of the error in the measurement of σ^ z and errors in basis-adjusting operations, the measured observables { Q_j } are not exactly { Q_j^(0)}. When the overall error is small, observables { Q_j } are still linearly independent. We introduce the matrix M^ out_j, σ = Q_jσ, and M^ out (0) is the matrix corresponding to { Q_j^(0)}. observables { Q_j } are linearly independent if M^ out is invertible. We have that M^ out is always invertible if M^ out - M^ out (0)_ max < 1/4s_ min(M^ out (0)) = 1/4. The initial state without error is decomposed as Q^(0) = ∑_j=1^4 q_j Q_j. Coefficients form a 4-dimensional row vector q = [q_1   q_2   q_3   q_4]. The decomposition is given by q = Q^(0) M^ out -1.



§ QUANTUM GATE SET TOMOGRAPHY



To measure a set of operations {_1, …, _N } on n qubits using GST, we need to choose a set of 4^n linearly independent initial states {ρ̅_k } and a set of 4^n linearly independent observables {Q̅_j }. Given these initial states and observables, we measure expected values

    _j, k = Q̅_jρ̅_k.

Here,  is one of operations {_1, …, _N }.

The matrix  is equivalent to  up to a transformation. Because _σ, τ = στ and ∑_σσσ = (the sum is taken over all Pauli operators), we have

    = M̅^ outM̅^ in,

where M̅^ in and M̅^ out are matrices defined as M̅^ in_σ, k = σρ̅_k and M̅^ out_j, σ = Q̅_jσ. We remark that initialisation error and measurement error are included in M̅^ in and M̅^ out, respectively. We cannot measure matrices M̅^ in and M̅^ out independently, therefore we cannot determine  using GST. By taking  as the identity operation (i.e. =) in Eq. (<ref>), we can measure

    g = M̅^ outM̅^ in.


The estimation of  is given by

    = T g^-1 T^-1 = T M̅^ in -1M̅^ in T^-1.

Here, g and  are obtained by measuring the expected values of observables, and T is an arbitrary invertible matrix. If M̅^ in and T are different,  is different from , but they are always similar matrices. The estimations of states ρ̅_k and observables Q̅_j are given by

    ρ̂_k   =    T_∙,k = TM̅^ in -1ρ̅_k, 
    Q̂_j   =    (g T^-1)_j,∙ = Q̅_jM̅^ in T^-1.

Here, M_∙,k (M_j,∙) denotes the k^th column (j^th row) of the matrix M.

We introduce matrices M̂^ in and M̂^ out defined as M̂^ in_σ, k = σρ̂_k and M̂^ out_j, σ = Q̂_jσ, respectively. Then M̂^ in = T and M̂^ out = g T^-1.

For a sequence of operations _1, …, _N, because _i and _i are similar matrices up to the same transformation independent of the operation (i.e. the index i), we have

    Q̅_j_N ⋯_1 ρ̅_k = Q̂_j_N ⋯_1 ρ̂_k.

Therefore, although estimations {ρ̂_k, Q̂_j, _i } may be different from their correspondences {ρ̅_k, Q̅_j, _i }, they can always provide the correct prediction for the expected value of an observable in an initial state going through a sequence of operations.



§ STABILITY OF THE QUANTUM GATE SET TOMOGRAPHY



We define

    ε̅_ in   =   max{M̅^ in_m - M^ in(0)_2 | m=1,…,n}, 
    ε̅_ out   =   max{M̅^ out_m - M^ out(0)_2 | m=1,…,n}, 
    ε̅_   =    - ^(0)_2,

which describe severities of the initialisation error, measurement error and operation error, respectively. Here, M̅^ in_m and M̅^ out_m are matrices corresponding to the m^th qubit. The overall matrices of n qubits are M̅^ in = ⊗_m=1^n M̅^ in_m and M̅^ out = ⊗_m=1^n M̅^ out_m.

Similar to the analyse of the linear independence of basis operations (i.e. the invertibility of the matrix A, see Appendix <ref>), we have that M̅^ in and M̅^ out are always invertible, i.e. g = M̅^ outM̅^ in is always invertible, if ε̅_ in < s_ min(M^ in(0)) and ε̅_ out < s_ min(M^ out(0)). Choosing {ρ_k^(0)} and Q^(0) as in Appendix <ref>, we have s_ min(M^ in(0)) = 1/2√(5-√(17)/2)≃ 0.3311 and s_ min(M^ out(0)) = 1.

We choose T = M^ in(0) ⊗ n, then M̂^ in = ⊗_m=1^n M̂^ in_m and M̂^ out = ⊗_m=1^n M̂^ out_m, where M̂^ in_m and M̂^ out_m are matrices corresponding to the m^th qubit. The severity of errors in estimations of initial states is

    ε̂_ in   =   max{M̂^ in_m - M^ in(0)_2  | m=1,…,n} = 0,

and the severity of errors in estimations of observables is

    ε̂_ out   =   max{M̂^ out_m - M^ out(0)_2  | m=1,…,n}
       ≤    (ε̅_ outε̅_ in + M^ in(0)_2ε̅_ out + M^ out(0)_2ε̅_ in) 
          ×M^ in(0)-1_2.

Here, we have used that

    M̂^ out_m - M^ out(0)
       =    (M̅^ out_mM̅^ in_m - M^ out(0)M^ in(0)) M^ in(0)-1
       =    [ (M̅^ out_m - M^ out(0))(M̅^ in_m - M^ in(0)) 
          + (M̅^ out_m - M^ out(0))M^ in(0)
          + M^ out(0)(M̅^ in_m - M^ in(0)) ] M^ in(0)-1.

Choosing {ρ_k^(0)} and Q^(0) as in Appendix <ref>, we have M^ in(0)_2 = 1/2√(5+√(17)/2)≃ 1.0679, M^ out(0)_2 = 1 and M^ in(0)-1_2 = s_ min^-1(M^ in(0)) = 2√(2/5-√(17))≃ 3.0204.

The severity of the error in the estimation of a n-qubit operation is

    ε̂_   =    - ^(0)_2 
       ≤    - _2 +  - ^(0)_2 
       ≤   2ε̅_ in^(n)/[s_ min(M^ in(0))]^n - ε̅_ in^(n) (^(0)_2 + ε̅_) 
          + ε̅_,

as we will show next. Here,

    ε̅_ in^(n) = (M^ in(0)_2 + ε̅_ in)^n - M^ in(0)_2^n.


For an invertible matrix A, we have

    A^-1_2    =   sup_x≠ 0A^-1x_2/x_2 = sup_y≠ 0y_2/Ay_2.

Then, using the inequality (<ref>), we have

    A^-1_2 ≤1/s_ min(A^(0)) - A-A^(0)_2.


We have the expression

    - 
       =    M^ in(0)⊗ nM̅^ in-1M̅^ in (M^ in(0)⊗ n)^-1 - 
       =    (M^ in(0)⊗ n - M̅^ in) M̅^ in-1
          × (M̅^ in - M^ in(0)⊗ n) (M^ in(0)⊗ n)^-1
          + (M^ in(0)⊗ n - M̅^ in) M̅^ in-1
          +  (M̅^ in - M^ in(0)⊗ n) (M^ in(0)⊗ n)^-1.

First, we have _2 ≤^(0)_2 + ε̅_. Second, using

    A⊗ B - C⊗ D_2 
       =   A⊗ B - A⊗ D + A⊗ D - C⊗ D_2 
       ≤   A_2 B - D_2 + A - C_2 D_2,

we have

    M̅^ in - M^ in(0)⊗ n_2 
       ≤   ε̅_ in∑_h = 1^n M^ in(0)_2^n-h∏_m=1^h-1M̅^ in_m_2 
       ≤   ε̅_ in∑_h = 1^n M^ in(0)_2^n-h (M^ in(0)_2 + ε̅_ in)^h-1
       =   ε̅_ in^(n).

Third, using the inequality (<ref>), we have

    M̅^ in-1_2 ≤1/[s_ min(M^ in(0))]^n - ε̅_ in^(n).


We remark that for a d-dimensional matrix M, M_ max≤M_2 ≤ dM_ max.



§ UPPER BOUND OF THE COST



We consider compensation method and take λ = 1, i.e. the n-qubit operation without error is realised as ^(0) =  +, where  is decomposed using basis operations as shown in Eq. <ref>. Then the cost for correcting the error in  is determined by

    C_ = 1 + ∑_i_1,…,i_nq_i_1,…,i_n.

Decomposition coefficients are determined by q = (A_1 ⊗⋯⊗ A_n)^-1 E, where q and E are defined in Eqs. (<ref>,<ref>). Here, q and E are 16^n-dimensional vectors, and A_1 ⊗⋯⊗ A_n is a 16^n-dimensional matrix. Therefore, for each element of q,

    q_i_1,…,i_n≤ 16^n (A_1 ⊗⋯⊗ A_n)^-1_ max_ max.

Here, we have used that the maximum absolute value of an element of E is _ max. Because = ^(0) -, we have _ max = ϵ_. Using the inequality (<ref>), we have

    (A_1 ⊗⋯⊗ A_n)^-1_ max
       =   ∏_l=1^n A_l^-1_ max≤∏_l=1^n A_l^-1_2 
       ≤   ∏_l=1^n 1/s_ min(A^(0))-A_l-A^(0)_2
       ≤   ∏_l=1^n 1/s_ min(A^(0))-16A_l-A^(0)_ max
       ≤   1/[s_ min(A^(0))-16ϵ_ max]^n.

Here, ϵ_ max = max{A_l-A^(0)_ max| l = 1,…,n }. There are total 16^n decomposition coefficients, therefore

    C_≤ 1 + 16^2nϵ_/(s_ min(A^(0))-16ϵ_ max)^n.

Here, we have assumed that ϵ_ max < s_ min(A^(0))/16.

We consider using the set of initial states with errors {ρ_k} to realise the initial state ρ_k_0^(0), which is in the set of initial states without error {ρ_k^(0)}. The initial state can be decomposed as ρ_k_0^(0) = ρ_k_0 + ∑_k (q_k-δ_k,k_0)ρ_k (see Appendix <ref>). We use ϵ_ in = M^ in - M^ in(0)_ max as the measure of the error severity in initial states. Using ρ_k_0 - ρ_k_0^(0)_ max≤ϵ_ in and M^ in-1_ max≤ 4[s_ min(M^ in(0)) - ϵ_ in]^-1, we have the cost for correcting errors in initial states

    C_ in   =   ∑_k q_k≤ 1 + ∑_k q_k - δ_k,k_0
       ≤    1 + 4^2ϵ_ in/s_ min(M^ in(0)) - ϵ_ in.


It is similar for observables. We consider using the set of observables with errors {Q_j} to realise the observable Q_j_0^(0), which is in the set of observables without error {Q_j^(0)}. We use ϵ_ out = M^ out - M^ out(0)_ max as the measure of the error severity in observables. Then, the cost for correcting errors in measured observables is

    C_ out≤ 1 + 4^2ϵ_ out/s_ min(M^ out(0)) - ϵ_ out.




§ ERROR MODELS



We consider a quantum computer with the following operations. The initialisation ^(0) = [π] + [πσ^ x], which prepares the state |0⟩. The projective measurement [π]. Single-qubit Clifford gates [R_ x] and [R_ z]. The single-qubit non-Clifford gate [T], where T = cosπ/8 + iσ^ xsinπ/8. Two-qubit maximally entangling gate [Λ], where Λ = 1/√(2)( + iσ^ z⊗σ^ z), which is equivalent to the controlled-NOT gate and controlled-phase gate up to single-qubit gates. The sixteen basis operations can be realised as shown in Table <ref>. In order to perform GST, we choose initial states and observables as in Appendix <ref>, and we choose the invertible matrix T = M^ in(0) ⊗ n for n qubits.

For the initialisation, the state prepared is ρ_0 rather than |0⟩⟨0|. We can always express the initialisation operation with error as = _ i, where _ i(|0⟩⟨0|) = ρ_0.

A POVM is defined by a set of operators { E_k } satisfying ∑_k E_k^† E_k =. In a POVM, the state is mapped to E_k ρ E_k^† when the outcome is k. When the measurement has error, we may not be able to obtain all the information k. Usually there are only two outcomes corresponding to |0⟩ and |1⟩, respectively. In this case, maybe several k values correspond to the same outcome ν = 0,1. Therefore, we model the projective measurement [π] with error as ρ = ∑_k∈ K_0 E_k ρ E_k^†, where K_ν is the set of k corresponding to the measurement outcome ν.

For a gate without error ^(0), the gate with error can be expressed as = _ a^(0)_ b. Any noisy gate can be expressed in this form: Because ^(0) is invertible, we can always take _ b = [] and _ a = ^(0)-1.

We suppose that time costs of the measurement [π] and the two-qubit gate [Λ] are the same, and time costs of single qubit gates are negligible.

We distinguish the identity operation and the memory operation. Without error, both of them are the same operation []. In any case, the identity operation is [], which means that the next operation is performed immediately, so it takes no time and there is not any memory error. When the memory operation is performed, the qubit waits for the next operation, so memory errors may occur on it. We apply the identity operation for measuring the matrix g (see Appendix <ref>). In the basis operation set, the operation [] is replaced by the memory operation.

We set the cycle time of the computing as the time cost of the measurement and the two-qubit gate. In one cycle, only one operation is performed on a qubit. If the operation is a single-qubit gate, the gate is performed at the middle of the cycle, i.e. the overall operation is _ m_ m, where _ m denotes memory noise. If no gate or measurement is performed in the cycle, the overall operation is _ m^2, which is the error version of the operation [] in the basis operation set.

We suppose that the single-qubit noise is described by the superoperator ^(1)(ϵ), and the two-qubit noise is described by the superoperator ^(2)(ϵ). Here, ϵ is a parameter describing the intensity of the noise. Then, the initialisation noise is _ i = ^(1)(ϵ/10), and the measurement with noise is M̃ = ℰ^(1)(ϵ/2)[π]^(1)(ϵ/2). For single-qubit gates, _ a = _ b = ^(1)(ϵ/20). For the two-qubit gate, _ a = _ b = ^(2)(ϵ/2). For the memory operation, _ m = ^(1)(ϵ/200).



 §.§ Depolarising Error


The single-qubit depolarising noise is

    ^(1)(ϵ) = (1-4ϵ/3) []
    + ϵ/3∑_α = 0^3 [σ^α],

where σ^0, σ^1, σ^2 and σ^3 correspond to , σ^ x, σ^ y and σ^ z, respectively. The two-qubit depolarising noise is

    ^(2)(ϵ) = (1-16ϵ/15) [⊗]
    + ϵ/15∑_α,β = 0^3 [σ^α⊗σ^β].

The x-axis (error rate) in Fig. <ref>(a) is ϵ of the two-qubit gate.



 §.§ Dephasing Error


The single-qubit dephasing noise is

    ^(1)(ϵ) = (1-ϵ) []
    + ϵ [σ^ z].

The two-qubit dephasing noise is

    ^(2)(ϵ)    =    (1-ϵ) [⊗]
    + ϵ/3( [⊗σ^ z] 
          +  [σ^ z⊗] +  [σ^ z⊗σ^ z]).




 §.§ Damping Error


The single-qubit damping noise is

    ^(1)(ϵ)    =    [+σ^ z/2 + √(1-ϵ)-σ^ z/2] 
          + [√(ϵ)σ^ x+iσ^ y/2].

The two-qubit damping noise is

    ^(2)(ϵ) = ℰ^(1)(ϵ/2)
    ⊗ℰ^(1)(ϵ/2).




 §.§ Over-rotation error


Noise is gate dependent. Initialisation, measurement and memory operation are perfect, i.e. ^(1) = [] for these operations. Only gates have noise. For gate R_ x, ^(1)(ϵ) = [cosϵπ/4 + iσ^ xsinϵπ/4]. For gate R_ x, ^(1)(ϵ) = [cosϵπ/4 + iσ^ zsinϵπ/4]. For gate T, ^(1)(ϵ) = [cosϵπ/8 + iσ^ zsinϵπ/8]. For gate Λ, ^(2)(ϵ) = [⊗cosϵπ/4 + iσ^ z⊗σ^ zsinϵπ/4]. The x-axis (over rotation) in Fig. <ref>(a) is ϵ of the two-qubit gate.



 §.§ Random-field error


Noise is gate dependent. For each operation, the noise ^(1,2)(ϵ) = [e^-iϵπ H] is determined by a Hamiltonian. Here, H = (h+h^†)/2, and each element of h is randomly generated with a uniform distribution in the unit circle. We remark that the noise is time independent, i.e. the noise is the same for the same gate implemented at different times.



 §.§ Random-operation error


The operation without noise is ^(0). The operation with noise is (ϵ), which depends on the error parameter. As the same as other models, the error parameter ϵ is operation dependent. Each operation can be expressed using a χ-matrix <cit.>. We suppose the χ-matrix corresponding to (ϵ) is χ, and the χ-matrix corresponding to ^(0) is χ^(0). To generate χ, firstly we generate a Hermitian matrix around χ^(0), which is χ' = χ^(0) + ϵ H, where H is generate as the same as the random Hamiltonian. Second, if ^(0) is not measurement, (ϵ) should be trace preserving. However, χ' may correspond to a non-trace preserving operation. In this case, we project χ' to the subspace in the matrix space that corresponds to trace preserving operations, i.e. χ” is the matrix closest to χ' and corresponds to a trace preserving operation. If ^(0) is measurement, χ” = χ'. Third, χ” may not be positive semi-definite. Therefore, we take χ”' = χ”+λ_ min if the minimum eigenvalue λ_ min of χ” is negative; otherwise χ”' = χ”. Finally, χ = fχ”', where the factor f makes sure that the operation is still trace preserving and the maximum eigenvalue of χ is smaller than 1.



§ INSTRUCTION OF THE IMPLEMENTATION OF THE QUASI-PROBABILITY METHOD



This section is a self-contained description of how to implement QEM using the quasi-probability decomposition. There are three steps: first, implement GST; second, compute the quasi-probability decomposition; third, implement the quasi-probability decomposition using the Monte Carlo approach.
 


 §.§ Implementation of gate set tomography


General discussions of GST are given in the main text and Appendix <ref>, therefore, here we describe GST in a more concrete way. GST is implemented to measure all gates used in the quantum computation. We discuss how to measure single qubit gates first and two-qubit gates afterwards.

To measure a single-qubit gate using GST, we prepare initial states |0⟩, |1⟩, |+⟩, and |y+⟩, where |+⟩, and |y+⟩ are the eigenstates of Pauli operators σ^ x and σ^ y with the eigenvalue +1, respectively. We denote these states ρ̅_1,  ρ̅_2, ρ̅_3 and ρ̅_4, respectively. These initial states can be noisy states, which is the essential advantage of GST, i.e. GST can tolerate state preparation and measurement errors. Then, we apply the gate that we want to measure, for instance, Hadamard gate, T gate and T^† gate in the SWAP-test circuit. Here we will use  (superoperator acting on a reduced density matrix) to denote the gate to be measured, which has noise. Subsequently, we measure expectation values for four observables, , σ^ x, σ^ y and σ^ z, respectively. Here  is a trivial observable whose measurement outcome is always +1. We denote these observables as Q̅_1, Q̅_2, Q̅_3, Q̅_4, and measurements of these observables can also be noisy. Then, by repeating the experiment to compute the mean value of observables, we can construct the 4× 4 matrix , and matrix elements are

    _j,k = ( Q̅_j ρ̅_k ).

Similarly, we can obtain the 4× 4 matrix g by choosing not to apply any gate to initial state, so that the matrix elements are

    g_j,k = ( Q̅_j ρ̅_k ).

This process is implemented for each qubit and each type of single-qubit gate (including all basis operations).

One will find that there is a freedom in the specification of the gate . Legitimate variants can be obtained as

    = T g^-1 T^-1,

where T is an invertible 4× 4 matrix. The matrix T can be different for different qubits but must be the same for all gates on the same qubit. We can choose T to minimise the cost in QEM. In the case that the error rate of preparing initial states is low, we can take

    T = [  1  1  1  1;  0  0  1  0;  0  0  0  1;  1 -1  0  0;    ],

which approximately minimise the cost according to our experience.

Estimations of the initial state ρ̅_k and the observable Q̅_j are respectively

    ρ̂_k   =    T_∙,k, 
    Q̂_j   =    (g T^-1)_j,∙.

Here, M_∙,k (M_j,∙) denotes the k^th column (j^th row) of the matrix M.

To measure a two-qubit gate using GST, the procedure is basically the same. The only difference is that there are 16 initial states and 16 observables to be measured. Initial states are the tensor products of single-qubit initial states, i.e. ρ̅^(1)_k_1⊗ρ̅^(2)_k_2, and observables are tensor products of single-qubit observables, i.e. Q̅^(1)_j_1⊗Q̅^(2)_j_2. Here, the superscript is the label of the qubit. Accordingly, the matrix g = g^(1)⊗ g^(2), which is the tensor product of g matrices of two qubits, and similarly the matrix T = T^(1)⊗ T^(2), which is the tensor product of T matrices of two qubits. We need to implement two-qubit gate GST for each pair of qubits that the two-qubit gate may be performed on.



 §.§ Quasi-probability decomposition


Using results obtained from GST, we can compute the quasi-probability decomposition. From GST we obtain estimations of initial states, observables to be measured, gates (including basis operations), and they are

    ρ̂_k    :    Initial state
    Q̂_j    :    Observable
        :    Gate
    _i     :    Basis operation

These estimations are utilised to compute the quasi-probability decomposition.

Now, we focus on the inverse method. We use ^(0) to denote the Pauli transfer matrix of the ideal gate (without error). The estimation of the Pauli transfer matrix of the actual gate with noise (i.e. ) is , which is obtained in GST. To compute the decomposition, first, we compute the ideal matrix ^(0); second, we compute the inverse of the noise

    ^-1 = ^(0)^-1;

and finally, we solve the equation (for single-qubit gate)

    ^-1 = ∑_i q_,i_i

to determine quasi-probabilities q_,i of the gate . We need to compute quasi-probabilities for each qubit and each gate. For example, for the SWAP-test circuit, we need to compute the decomposition for Hadamard gate, T gate, T^† gate of each qubit and controlled-NOT gate of each pair of qubits that the controlled-NOT gate may be performed on.

For two-qubit gates, the procedure is the same but tensor products of single-qubit basis operations, i.e. ^(1)_i_1⊗^(2)_i_2 (where the superscribe is the label of the qubit), are used to decompose the inverse of the noise ^-1.

In order to mitigate errors in initial states and measurements of observables, we should solve the following equations for the quantities q^(m)_ρ,k and q^(m)_Q,j:

    ρ^(0)   =   ∑_k q^(m)_ρ,kρ̂_k, 
    Q^(0)   =   ∑_k q^(m)_Q,jQ̂_j

for each qubit. Here, m is the label of the qubit, ρ^(0) is the column vector representing the ideal initial state |0⟩⟨0|, and Q^(0) is the row vector representing the ideal observable σ^ z. 

Before implementing the quasi-probability decomposition on a quantum computer, we compute

    C^(m)_ρ   =   ∑_k q^(m)_ρ,k, 
    
    C^(m)_Q    =   ∑_k q^(m)_Q,j

for each qubit and

    C_   =   ∑_i q_,i

for each gate.



 §.§ Monte Carlo implementation of the quasi-probability decomposition


It is vital to note that we use estimations _i to decompose the inverse of the noise, but usually there is a difference between _i and the actual basis operation _i. This difference does not cause any error in the final computing result, because the computing result is invariant under a similarity transformation as we have explained in the main text.

Now, we describe how to implement the quasi-probability decomposition on a quantum computer. We suppose the circuit is sequentially performing gates _1, _2, …, _N on the initial state |00… 0⟩, and the first qubit is measured in the σ^ z basis to read the computing result. The procedure can be generalised to the case of measuring multiple qubits.

First, we generate a set of random integers: for each qubit m, we randomly select an integer k_m such that each integer would be selected with corresponding probability q^(m)_ρ,k_m/C^(m)_ρ; similarly for each gate l, we generate random integer i_l with corresponding probability q__l,i_l/C__l; and finally we generate random integer j_1 with the probability q^(1)_Q,j_1/C^(1)_Q.

Second, on the quantum computer, we implement the following quantum computing for once: we initialise the qubit m in the state ρ̅^(m)_k_m; then we sequentially perform gates _1, B̅_i_1, _2, B̅_i_2, …, _N, B̅_i_N; finally we measure the observable Q̅_j_1. The measurement outcome is μ.

Third, we compute the effective measurement outcome

    μ_ eff =  sgn( ∏_m q^(m)_ρ,k_m∏_l q__l,i_l q^(1)_Q,j_1) μ.


By repeating these three steps, we can obtain the mean of effective outcomes E[μ_ eff]. The final computing result is C  E[μ_ eff], where

    C = ∏_m C^(m)_ρ∏_l C__l C^(1)_Q.





§ ERROR COMPONENT OF PAULI ERROR AND LEAKAGE ERROR



Taking ϵ = p_ x + p_ y + p_ z, it is obvious that in the inhomogeneous Pauli error model, the error component does not depend on the error rate, i.e. ℰ' = ϵ^-1(p_ x[σ^ x] + p_ y[σ^ y] + p_ z[σ^ z]). We remark that ratios p_α/ϵ do not change with ϵ.

For the leakage error model, we take ϵ = p. Then the error component is

    ℰ'(ρ)    =   π_0 ρπ_0 + √(1-p) - (1-p)/p (π_0 ρπ_1 + π_1 ρπ_0) 
       =   π_0 ρπ_0 + [1/2+𝒪(p)] (π_0 ρπ_1 + π_1 ρπ_0),

where π_0 = |0⟩⟨0| and π_1 = |1⟩⟨1|. Therefore, ℰ' varies slowly with p.

9

Li2017 Y. Li and S. C. Benjamin,
Efficient variational quantum simulator incorporating active error minimization,
Phys. Rev. X 7, 021050 (2017).

Temme2017 K. Temme, S. Bravyi, and J. M. Gambetta,
Error mitigation for short-depth quantum circuits,
Phys. Rev. Lett. 119, 180509 (2017).



Joe2016 J. O'Gorman and E. T. Campbell,
Quantum computation with realistic magic state factories,
Phys. Rev. A 95, 032338 (2017).



Peruzzo2014 A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik, and J. L. O'Brien,
A variational eigenvalue solver on a photonic quantum processor,
Nat. Commun. 5, 4213 (2014).

Wecker2015 D. Wecker, M. B. Hastings, and M. Troyer,
Progress towards practical quantum variational algorithms,
Phys. Rev. A 92, 042303 (2015).

McClean2016 J. R. McClean, J. Romero, R. Babbush, and A. Aspuru-Guzik,
The theory of variational hybrid quantum-classical algorithms,
New J. Phys. 18, 023023 (2016).

Bauer2016 B. Bauer, D. Wecker, A. J. Millis, M. B. Hastings, and M. Troyer,
Hybrid quantum-classical approach to correlated materials,
Phys. Rev. X 6, 031045 (2016).

Kreula2016 J. M. Kreula, S. R. Clark, and D. Jaksch,
Non-linear quantum-classical scheme to simulate non-equilibrium strongly correlated fermionic many-body dynamics,
Sci. Rep. 6, 32940 (2016).

Kreula2016EPJ J. M. Kreula, L. García-Álvarez, L. Lamata, S. R. Clark, E. Solano, and D. Jaksch,
Few-qubit quantum-classical simulation of strongly correlated lattice fermions,
EPJ Quantum Technology 3, 11 (2016).



Shen2017  Y. Shen, X. Zhang, S. Zhang, J.-N. Zhang, M.-H. Yung, and K. Kim,
Quantum implementation of unitary coupled cluster for simulating molecular electronic structure,
Phys. Rev. A 95, 020501 (2017).

OMalley2016 P. J. J. O'Malley et al.,
Scalable quantum simulation of molecular energies,
Phys. Rev. X 6, 031007 (2016).

Kandala2017 A. Kandala, A. Mezzacapo, K. Temme, M. Takita, M. Brink, J. M. Chow, and J. M. Gambetta,
Hardware-efficient variational quantum eigensolver for small molecules and quantum magnets,
Nature 549, 242 (2017).



McClean2017  J. R. McClean, M. E. Kimchi-Schwartz, J. Carter, and W. A. de Jong,
Hybrid quantum-classical hierarchy for mitigation of decoherence and determination of excited states,
Phys. Rev. A 95, 042308 (2017).

Colless2018  J. I. Colless, V. V. Ramasesh, D. Dahlen, M. S. Blok, J. R. McClean, J. Carter, W. A. de Jong, I. Siddiqi,
Computation of molecular spectra on a quantum processor with an error-resilient algorithm,
Phys. Rev. X 8, 011021 (2018).



Merkel2013 S. T. Merkel, J. M. Gambetta, J. A. Smolin, S. Poletto, A. D. Córcoles, B. R. Johnson, C. A. Ryan, and M. Steffen,
Self-consistent quantum process tomography,
Phys. Rev. A 87, 062119 (2013).

Greenbaum2015 D. Greenbaum,
Introduction to quantum gate set tomography,
arXiv:1509.02921.



Knill2004 E. Knill,
Fault-tolerant postselected quantum computation: threshold analysis,
arXiv:quant-ph/0404104.

Wallman2015 J. J. Wallman and J. Emerson,
Noise tailoring for scalable quantum computation via randomized compiling,
arXiv:1512.01098.

OGorman2016 J. O'Gorman, N. H. Nickerson, P. Ross, J. J. L. Morton, and S. C. Benjamin,
A silicon-based surface code quantum computer,
npj Quantum Information 2, 15019 (2016).



Ekert2002A. K. Ekert, C. Moura Alves, and D. K. L. Oi,
Direct estimations of linear and nonlinear functionals of a Quantum State, 
Phys. Rev. Lett. 88, 217901(2002).

shallowswap L. Cincio, Y. Subaşi, A. T. Sornborger, P. J. Coles,
Learning the quantum algorithm for state overlap, 
arXiv: 1803.04114. 

AspuruGuzikAutoEncJonathan Romero, Jonathan P Olson, and Alan Aspuru-Guzik,
Quantum autoencoders for efficient compression of quantum data,
Quantum Sci. Technol. 2 045001 (2017).



Song2017 C. Song, K. Xu, W. Liu, C. Yang, S.-B. Zheng, H. Deng, Q. Xie, K. Huang, Q. Guo, L. Zhang, P. Zhang, D. Xu, D. Zheng, X. Zhu, H. Wang, Y.-A. Chen, C.-Y. Lu, S. Han, and J.-W. Pan,
10-qubit entanglement and parallel logic operations with a superconducting circuit,
Phys. Rev. Lett. 119, 180511 (2017).



Bravyi2016 S. Bravyi and D. Gosset,
Improved classical simulation of quantum circuits dominated by Clifford gates,
Phys. Rev. Lett. 116, 250501 (2016).

Bravyi2016PRX S. Bravyi, G. Smith, and J. A. Smolin,
Trading classical and quantum computational resources,
Phys. Rev. X 6, 021043 (2016).

Howard2017 M. Howard and E. Campbell,
Application of a resource theory for magic states to fault-tolerant quantum computing,
Phys. Rev. Lett. 118, 090501 (2017).



Barends2014 R. Barends, J. Kelly, A. Megrant, A. Veitia, D. Sank, E. Jeffrey, T. C. White, J. Mutus, A. G. Fowler, B. Campbell, Y. Chen, Z. Chen, B. Chiaro, A. Dunsworth, C. Neill, P. O'Malley, P. Roushan, A. Vainsencher, J. Wenner, A. N. Korotkov, A. N. Cleland, and J. M. Martinis,
Superconducting quantum circuits at the surface code threshold for fault tolerance,
Nature 508, 500 (2014).

Harty2014 T. P. Harty, D. T. C. Allcock, C. J. Ballance, L. Guidoni, H. A. Janacek, N. M. Linke, D. N. Stacey, and D. M. Lucas,
High-fidelity preparation, gates, memory and readout of a trapped-ion quantum bit,
Phys. Rev. Lett. 113, 220501 (2014).

Lucas C. J. Ballance, T. P. Harty, N. M. Linke, M. A. Sepiol, and D. M. Lucas,
High-fidelity quantum logic gates using trapped-ion hyperfine qubits,
Phys. Rev. Lett. 117, 060504 (2016).

Wineland J. P. Gaebler, T. R. Tan, Y. Lin, Y. Wan, R. Bowler, A. C. Keith, S. Glancy, K. Coakley, E. Knill, D. Leibfried, and D. J. Wineland,
High-fidelity universal gate set for ^9Be^+ ion qubits,
Phys. Rev. Lett. 117, 060505 (2016).



Nielsen2010 M. A. Nielsen and I. L. Chuang,
Quantum Computation and Quantum Information,
Cambridge University Press, Cambridge, (2010).



