Efficient estimation of Pauli channels

TL;DR

This paper presents efficient methods for learning and estimating Pauli channels in quantum systems, enabling accurate characterization of noise and errors with fewer measurements, which is crucial for quantum computing reliability.

Contribution

The paper introduces new measurement-efficient algorithms for estimating Pauli channels and their error rates, including channels with local correlations, advancing quantum noise characterization.

Findings

Achieves high-probability estimation of n-qubit Pauli channels with O(ε^{-2} n 2^n) measurements.

Provides measurement bounds for estimating s Pauli error rates with O(ε^{-4} log s log s/ε).

Learns channels with k-local correlations using O_k(ε^{-2} n^2 log n) measurements, efficient in qubit number.

Abstract

Pauli channels are ubiquitous in quantum information, both as a dominant noise source in many computing architectures and as a practical model for analyzing error correction and fault tolerance. Here we prove several results on efficiently learning Pauli channels, and more generally the Pauli projection of a quantum channel. We first derive a procedure for learning a Pauli channel on $n$ qubits with high probability to a relative precision $\epsilon$ using $O\bigl(\epsilon^{-2} n 2^n\bigr)$ measurements, which is efficient in the Hilbert space dimension. The estimate is robust to state preparation and measurement errors which, together with the relative precision, makes it especially appropriate for applications involving characterization of high-accuracy quantum gates. Next we show that the error rates for an arbitrary set of $s$ Pauli errors can be estimated to a relative precision $\epsilon$ using $O\bigl(\epsilon^{-4} \log s\log s/\epsilon\bigr)$ measurements. Finally, we show that when the Pauli channel is given by a Markov field with at most $k$-local correlations, we can learn an entire $n$-qubit Pauli channel to relative precision $\epsilon$ with only $O_k\bigl(\epsilon^{-2} n^2 \log n \bigr)$ measurements, which is efficient in the number of qubits. These results enable a host of applications beyond just characterizing noise in a large-scale quantum system: they pave the way to tailoring quantum codes, optimizing decoders, and customizing fault tolerance procedures to suit a particular device.

Equations

$$\Vert \hat{\mathbf{p}}-\mathbf{p}{\Vert }_2\le O(ϵ)(1-p_0),$$

$$\Vert \hat{\mathbf{p}}-\mathbf{p}{\Vert }_{\infty }$$

$$\Vert \hat{\mathbf{p}}-\mathbf{p}{\Vert }_1\le O(ϵ)\Vert 1_I-\mathbf{p}{\Vert }_{\infty }$$

$$P_a=P_{(a_x,a_z)}=i^{a_x\cdot a_z}X[a_x]Z[a_z]$$

$$⟨a,b⟩=a_x\cdot b_z+a_z\cdot b_x\mathrm{mod}2$$

$${\mathsf{C}}_{\mathsf{G}}=\{a\in {\mathsf{P}}^n:\forall g\in \mathsf{G},⟨a,g⟩=0\}.$$

$$b=a+c$$

$$1[x]=\{\begin{array}{ll}1& \text{if x is true},\\ 0& \text{if x is false}.\end{array}$$

$$\frac{1}{|\mathsf{G}|}\sum _{b\in \mathsf{G}}(-1{)}^{⟨a,b⟩}=1[a\in {\mathsf{C}}_{\mathsf{G}}].$$

$$\displaystyle|{\sf{O}}|\leq\min\bigl{\{}|{\sf{X}}|,\sqrt{|\langle{\sf{X}}\rangle/{\sf{S_{X}}}|}+1\bigr{\}}\,,$$

$${\rho }_{\mathsf{S},a}=\frac{1}{d}\sum _{s\in \mathsf{S}}(-1{)}^{⟨s,a⟩}{P}_s$$

$$E_{\mathsf{S}}=\left\{E_{\mathsf{S},e}=\frac{d}{|{\mathsf{A}}_{\mathsf{S}}|}{\rho }_{\mathsf{S},e}:e\in {\mathsf{A}}_{\mathsf{s}}\right\}$$

$$\mathcal{L}(M)=\sum _k{A}_{k}M{B}_k^{†}.$$

$${\mathcal{P}}_a(P_b)=P_a{P}_b{P}_a^{†}=(-1{)}^{⟨a,b⟩}{P}_b.$$

$$\mathcal{E}(\rho )=\sum _{a\in {\mathsf{P}}^n}{p}_a{P}_a\rho P_a$$

$${\mathcal{L}}^{\mathbb{T}}=\frac{1}{|\mathbb{T}|}\sum _{T\in \mathbb{T}}\mathcal{T}\mathcal{L}{\mathcal{T}}^{†}.$$

$${\mathcal{L}}^{\mathbb{T}}(M)=\frac{1}{|\mathbb{T}|}\sum _{T\in \mathbb{T}}\sum _k\mathcal{T}(A_k)\rho \mathcal{T}(B_k{)}^{†}.$$

$$p_a=\sum _k{l}_{k,a}{r}_{k,a}^{\ast }.$$

$${\mathcal{L}}^{{\mathsf{P}}^n}(M)=|{\mathsf{P}}^n{|}^{-1}\sum _{c\in {\mathsf{P}}^n}\sum _k{\mathcal{P}}_c(A_k)M{\mathcal{P}}_c(B_k{)}^{†}.$$

$$\begin{array}{rl}{\mathcal{L}}^{{\mathsf{P}}^n}(M)& =|{\mathsf{P}}^n{|}^{-1}\sum _k\sum _{a,b,c\in {\mathsf{P}}^n}{l}_{k,a}{r}_{k,b}^{\ast }(-1{)}^{⟨c,a+b⟩}{P}_{a}M{P}_b\\ & =\sum _k\sum _{a,b\in {\mathsf{P}}^n}{l}_{k,a}{r}_{k,b}^{\ast }1[a=b]P_{a}M{P}_b\\ & =\sum _{a\in {\mathsf{P}}^n}\left(\sum _k{l}_{k,a}{r}_{k,a}^{\ast }\right)P_{a}M{P}_a,\end{array}$$

$$\mathcal{E}=\frac{1}{d}\sum _{a,b\in {\mathsf{P}}^n}(-1{)}^{⟨a,b⟩}{p}_b|P_a)(P_a|.$$

$$|a)=\frac{1}{\sqrt{d}}|P_a)$$

$${\mathbf{W}}_{\mathsf{A},\mathsf{B}}=\sum _{a\in \mathsf{A},b\in \mathsf{B}}(-1{)}^{⟨a,b⟩}|a)(b|.$$

$$\mathbf{f}=\mathbf{W}\mathbf{p},$$

$$\mathcal{E}=\sum _{a\in {\mathsf{P}}^n}{f}_a|a)(a|.$$

$$r_h=1-f_h,{\mathbf{r}}_{\mathsf{A}}=\sum _{a\in \mathsf{A}}{r}_a|a),1_{\mathsf{A}}=\sum _{a\in \mathsf{A}}|a).$$

$${\mathbf{W}}_{\mathsf{A},\mathsf{B}}^{†}{\mathbf{W}}_{\mathsf{A},\mathsf{B}}=|\mathsf{A}|\sum _{b,b^{\prime }\in \mathsf{B}:b+b^{\prime }\in {\mathsf{C}}_{\mathsf{A}}}|b^{\prime })(b|.$$

$${\mathbf{W}}_{\mathsf{A},\mathsf{B}}^{†}{\mathbf{W}}_{\mathsf{A},\mathsf{B}}$$

$$=\sum _{b,b^{\prime }\in \mathsf{B}}\left(\sum _{a\in \mathsf{A}}(-1{)}^{⟨a,b+b^{\prime }⟩}\right)|b^{\prime })(b|$$

$$=|\mathsf{A}|\sum _{b,b^{\prime }\in \mathsf{B}}1[b^{\prime }+b\in {\mathsf{C}}_{\mathsf{A}}]|b^{\prime })(b|,$$

$$\frac{1}{2}\Vert |0)-\mathbf{p}{\Vert }_1={ϵ}_{\diamond }(\mathcal{E})=(1+1/d)r(\mathcal{E}),$$

Full text

I Introduction

Pauli channels are among the most basic noise channels in quantum information science. Nearly all theoretical studies of quantum error correction and fault tolerance, including most threshold and overhead estimates, rely on modeling noise as a Pauli channel [1].

Originally, the theoretical focus on Pauli channels was primarily motivated by the ease with which they can be theoretically analyzed and simulated on a classical computer. However, Pauli channels are now justified by the technique of randomized compiling [2], which maps general quantum noise to its Pauli projection, a Pauli channel having the same average fidelity to an ideal channel. This addresses the concern that coherent noise, which is in general as hard to model as full quantum computation, may create distorted comparisons with threshold error rates that were computed using Pauli noise [3, 4]. Moreover, quantum error correction of coherent noise will, under reasonable assumptions, lead to less coherent noise at the logical level, meaning that noise at that level is better approximated by a Pauli channel [5, 6]. Finally, recent experiments [7] have shown that enforcing the Pauli projection by randomized compiling works extremely well in practice, further solidifying the importance of Pauli channels.

Despite the central role played by Pauli channels, to date there have been few systematic studies of how to estimate them efficiently, meaning with a complexity that improves over what follows from a naive application of full channel tomography. We will review the relevant literature below. Filling this gap becomes even more pressing in light of work showing just how much the threshold changes under biased or correlated noise models [8, 9, 10, 11], in some cases by more than a factor of 4 for a code capacity threshold [12, 13, 14]. Such factors can be substantial because the logical error rate depends exponentially on the distance below the threshold. Accurate estimation of the Pauli noise in an architecture would allow for many applications, including tailoring codes and decoders to match the noise [15], customizing fault-tolerance schemes [16], and accurate estimation of thresholds and overheads [17].

I.1 Summary of Results

In this paper, we give a comprehensive treatment of the sample complexity of Pauli channel estimation. We will present three main results corresponding to efficient estimation procedures with error guarantees for three separate and complementary regimes: estimation of a complete Pauli channel, estimation of error rates for an arbitrary set of Pauli errors, and estimation of a Pauli channel that factorizes over a bounded-degree factor graph. The latter is equivalent to considering the Pauli error probability as a tensor network state comprised of tensors with a bounded number of indices, but no restriction on the topology of the connections. This section contains only informal statements of our results; the precise versions of these statements along with rigorous proofs are given in the subsequent sections.

The basic procedure that we develop and analyze is a variant of randomized benchmarking [18, 19, 20] and its recently introduced cousins, character benchmarking [21] and cycle benchmarking [22]. The procedure, defined in Section III, inherits the robustness to errors in state preparation and measurement (SPAM) enjoyed by all benchmarking variants. It uses preparations and measurements in a stabilizer basis (or equivalently, in the computational basis together with a single Clifford group element), and repeated rounds of random Pauli gates to average over the noise. Our results apply when the noise on the random Pauli gates, the state preparations, and the measurement all obey certain mild regularity assumptions such as being gate-independent, time-stationary, Markovian, and not too far from ideal. These assumptions can be relaxed still further, particularly the gate-independence assumption [23, 24, 25, 26]. The precise conditions on the noise for which our proofs hold are given in Definitions 2 and 3. In this section, we informally refer to these restrictions as “nice” noise.

Our first main result is that our procedure for estimating a complete Pauli channel on $n$ qubits requires only O\bigl{(}\epsilon^{-2}n2^{n}\bigr{)} measurements to estimate the channel to constant relative precision $\epsilon$ and with a constant success probability. This result forms a core subroutine in our subsequent two results. The output is the vector $\bm{p}$ of all $4^{n}$ Pauli error probabilities, which we achieve by using the full power of $n$-bit measurements. Here the relative precision and robustness to SPAM are crucial for applications, since average error rates in quantum gates are now routinely below 1%, and in some cases as low as $10^{-6}$ [27]. In this regime, a meaningful additive error approximation would require at least $10^{12}$ samples, making it far outside the realm of practicality even if it is still technically “efficient”.

Result 1** (Informal summary of Propositions 8 and 9).**

For a nice Pauli channel on $n$ qubits, the Pauli error rates $\bm{p}$ can be estimated using O\bigl{(}\epsilon^{-2}n2^{n}\bigr{)} measurements by $\hat{\bm{p}}$ such that

[TABLE]

holds with high probability, where $p_{0}$ is the probability of no error.

We can apply a randomized sampling routine to estimate a subnormalized probability distribution over an arbitrary set ${\sf{E}}$ of size $|{\sf{E}}|=s$. We note that the $\epsilon^{-4}$ scaling is an artifact of the proof technique and expect it can be improved to $\epsilon^{-2}$ by a more careful analysis of the bias in the output of the subroutine ${\bf Ratio}$ described in section IV.

Result 2** (Informal statement of Theorem 11).**

For a nice Pauli channel on $n$ qubits, a subnormalized distribution over any set ${\sf{E}}$ of $s$ Pauli errors can be estimated using $O(\epsilon^{-4}\log(s)\log(s/\epsilon^{2}))$ measurements by $\hat{\bm{p}}$ such that

[TABLE]

holds with constant probability, where $p_{0}$ is the probability of no error.

The above sampling protocol can be applied directly to efficiently estimate the probability of all low-weight Pauli errors, which will account for the majority of the distribution under realistic physical assumptions in near-term high-performance quantum systems. We also provide a tree-based search heuristic that can be used to identify sets of errors of interest. Two natural candidates for interesting sets of errors are the best $s$-sparse approximation without assuming local errors, and high-weight Pauli errors that occur with probabilities that differ substantially from model predictions. Although we do not give a formal analysis of the probability with which the search heuristic identifies a correct set, when it finds a sparse set with large measure it is certifiably correct as a description of the Pauli channel. This follows because Theorem 11 can be applied to the final output without accounting for failure probabilities from intermediate steps.

Our final main result is the efficient reconstruction of any nice $n$-qubit Pauli channel in polynomial time in $n$ whenever the channel has bounded-degree correlations and the local marginals are positive, conditions which we quantify in section VII. We assume throughout that the correlations are modeled by a factor graph with fixed known topology.

Result 3** (Informal statement of Proposition 18.).**

Let $\bm{p}$ be a nice $n$-qubit Pauli channel with a $k$-degree factor graph and positive marginals. Then an estimate $\hat{\bm{p}}$ as a tensor network can be obtained using O_{k}\bigl{(}\epsilon^{-2}n^{2}\log n\bigr{)} measurements such that

[TABLE]

holds with high probability, where $O_{k}$ means the implicit big-$O$ constant depends exponentially on $k$. Moreover, an estimate proportional to $\hat{\bm{p}}$ can be found in time $\mathrm{poly}(n)$.

Here we use the notation $\bm{1}_{I}$ to mean a probability distribution supported only on the identity Pauli operator, which corresponds to a noiseless channel.

I.2 Proof techniques

As mentioned above, we make use of a variant of randomized benchmarking (RB). RB was originally applied to the full unitary group and the Clifford group [18, 19], but subsequent work has developed it for many different groups [28, 29, 30, 31, 32, 33, 34, 35]; see Refs. [36, 21] for a general treatment. Our work focuses on RB over the Pauli group, and introduces several important advances over prior art.

As shown by Harper et al. [37], RB over the Clifford group provably gives relative precision estimates of a single parameter associated to a quantum channel, the average error rate. Our proof extends this to the case of RB over the Pauli group. The original proof of Ref. [37] isolates a single exponential decay from an RB signal using the method from Ref. [38]. However, for the case of interest here, the Pauli group, there are $4^{n}$ parameters that must be learned instead of just one. While a general solution exists to the problem of isolating exponential decays [21], prior work estimates the individual decay rates serially using one-bit measurements. Our work greatly improves this by showing how to take the natural $n$-bit measurements arising from measuring the $n$ individual qubits and using these measurements to estimate up to $2^{n}$ parameters in parallel. When sampled in this way, a Hadamard transformation can be done to isolate the exponentials on all data collected from the $n$-bit measurements simultaneously, avoiding the need for serial exponential fits. Our proof for the case of estimating probabilities for arbitrary sets of Pauli errors uses the essential ideas above with random sampling techniques.

Finally, we assume that the Pauli channel error rates are given by a Markov random field over a known factor graph. We can again apply the core subroutine for estimating Pauli channels to the marginal channel on each of the factors and then round these estimated marginals into a coherent global probability distribution using the Hammersley-Clifford theorem [39].

I.3 Comparison with Prior Work

These results comprise the first proofs in the literature of recovery guarantees for quantum channels to relative precision while avoiding bias from SPAM errors. As noted above, this is absolutely essential for applications. Furthermore, Results 2 & 3 are the first recovery guarantees known that apply to broad classes of non-trivial quantum channels on $n$ qubits with a scaling that is efficient in $n$ (when $s=\mathrm{poly}(n)$ or $k=O(1)$, respectively). Lastly, although numerical practicality is not a focus of the present paper (it will be explored elsewhere [40]) we wish to stress that these results are practical and immediately applicable to characterizing error rates in near-term devices, such as the 50- and 72-qubit processors announced by IBM and Google respectively. To emphasize this last point, the results of our numerical simulations of these procedures on up to 100 qubits show excellent performance, and we have successfully implemented a variant of the method on a publicly accessible 14-qubit quantum device [40]. To put this in perspective, the previous record for quantum channel tomography is of a 3-qubit quantum gate [41], which was only made possible by employing sophisticated methods from compressed sensing [42, 43, 44]. It is clear from these examples that our results represent a qualitative shift in the characterization of quantum devices.

Let us compare the scaling of our procedures to the best previously known results from the literature. A naive, non-adaptive application of channel tomography with (single-copy) two-outcome measurements [45] would require \tilde{O}\bigl{(}\tfrac{d^{4}}{\epsilon^{2}}\bigr{)} state preparations and measurements to achieve an additive precision $\epsilon$ in, say, the average fidelity of the channel, where $d=2^{n}$ for $n$ qubits. Applying the best known procedures for state tomography [46, 47, 48] to the Choi state of the channel using a collective measurement requires \tilde{\Theta}\bigl{(}\tfrac{d^{4}}{\epsilon^{2}}\bigr{)} copies, so this is roughly tight. For sparse channels in a known superoperator basis (such as the Pauli basis), Ref. [44] argues from the theory of compressed sensing that $O(s\log d)$ random product basis measurements should suffice to reconstruct an $s$-sparse channel in the Pauli basis, but it is not clear that the technical “incoherence” or restricted isometry properties that are used for efficiently reconstructing sparse vectors are satisfied in this case [43]. The recently introduced procedure of shadow tomography [49] could achieve an estimate of the Pauli error rates with a small ($\mathrm{poly}\log d$) number of copies of the Choi matrix for the channel, but at the cost of a large ($\mathrm{poly}(d)$) amount of nontrivial quantum computation and a copy complexity of $~{}1/\epsilon^{5}$ with respect to the additive error. Our procedure uses much more trivial quantum resources, and with a circuit depth that is independent of the dimension. The idea of tomography of matrix product states [50] has been extended to matrix product operators that are unitary in the numerical work of Ref. [51]. There is as yet no systematic mathematical treatment of matrix product operator tomography, or tomography of quantum channels with bounded correlations of the kind discussed in this paper.

The above tomographic techniques only provide additive precision estimates rather than the relative precision estimates achieved here. Furthermore, these tomographic methods are not robust against SPAM errors, which create an unknown systematic error, or, equivalently, places an absolute floor on the achievable additive precision that can be obtained using these tomographic techniques [52]. The magnitude and complexity of SPAM errors scales with the system size, so that the achievable precision decreases for larger numbers of qubits. These tomographic techniques could be made robust to SPAM errors using gate-set tomography [52, 53], at the cost of substantially increasing the resources required. However, it is unclear how to obtain relative-precision estimates from any of these methods.

We achieve multiplicative precision and robustness to SPAM errors by using techniques originating from randomized benchmarking. Randomized benchmarking tomography [54] also employs randomized benchmarking techniques to achieve robustness to SPAM, and can obtain tomographic reconstructions using O\bigl{(}d^{2}\log(d)\bigr{)} randomized benchmarking experiments with an additive reconstruction error that scales as $O(d/\sqrt{\log d})$ [55]. Note that the additive error is essentially fundamental to the approach used in randomized benchmarking tomography because the decay rates are the fidelities between the noisy process and distinct Clifford gates. These fidelities are significantly less than 1—even for ideal noise processes—resulting in rapid decays that are difficult to estimate precisely.

Besides the prior work on quantum tomography, there has also been a large amount of prior work on learning Markov random fields, both in terms of parameter estimation (including hardness results [56]) and in terms of learning the structure of the factor graph [57, 58, 59, 60, 61, 62]. These results are certainly relevant, but to make a direct comparison to the present results is complicated by the fact that our probability distribution contains variables that might be called “quasi-latent”, in the sense that they are not observable to all queries, only to some, and queries that probe one variable must necessarily hide others. This is due to the symplectic structure of the Pauli group: we can only simultaneously measure observables that commute. Therefore learning a “symplectic Markov field” seems to be an inherently different task from what has been previously considered.

II Mathematical Preliminaries

Given a set of $n$ qubits with Hilbert space dimension $d=2^{n}$, let us introduce the following notation. Let $\mathbb{P}^{n}$ denote the group of Pauli operators acting on $n$ qubits, and ${\sf{P}}^{n}=\mathbb{P}^{n}/\langle i\rangle$ be the quotient of $\mathbb{P}^{n}$ with its center. The group ${\sf{P}}^{n}$ is Abelian and isomorphic to $\mathbb{Z}_{2}^{2n}$, so we will label elements in ${\sf{P}}^{n}$ by bit strings of length $2n$. We will abuse notation and talk about $a\in{\sf{P}}^{n}$ when we mean a matrix $P_{a}\in{\sf{P}}^{n}$ up to an overall phase. For concreteness, for a bit string $a$ and a single-qubit Pauli operator $A$ we write $A[a]=\bigotimes_{a_{i}\in a}A^{a_{i}}$ and choose the elements of ${\sf{P}}^{n}$ to be

[TABLE]

where $a\in\mathbb{Z}_{2}^{2n}$ and $X$ and $Z$ are the standard single-qubit Pauli matrices. The choice of phase is arbitrary but this choice of phase ensures that every $P_{a}$ is Hermitian. For any two Pauli matrices $P_{a}$ and $P_{b}$ with respective images $a$ and $b$ in $\mathbb{Z}_{2}^{2n}$, we have $P_{a}P_{b}=(-1)^{\langle a,b\rangle}P_{b}P_{a}$ where

[TABLE]

is a binary symplectic form, and so it is symmetric and linear in each argument. We will typically omit the mod 2 because it is clear from context.

We define a stabilizer group to be a linear subspace of $\mathbb{Z}_{2}^{2n}$ such that $\langle a,b\rangle=0$ for all $a,b\in{\sf{S}}$. Therefore every stabilizer group is a group ${\sf{S}}\subset\sf{P}^{n}$ whose elements all commute as matrices. Note that stabilizer groups are often defined in terms of the Pauli matrices with specific phase conventions [63].

The commutant of a set ${\sf{G}}\subseteq{\sf{P}}^{n}$ is the group ${\sf{C_{G}}}\subseteq{\sf{P}}^{n}$ that is orthogonal to ${\sf{G}}$ according to the binary symplectic form. That is,

[TABLE]

We refer to ${\sf{C_{G}}}$ as the commutant because every element of ${\sf{C_{G}}}$ commutes with every element of ${\sf{G}}$ as matrices in $\mathbb{P}^{n}$. When ${\sf{G}}$ is a group (i.e., closed under multiplication), ${\sf{C_{C_{G}}}}={\sf{G}}$ by the double commutant theorem. We define the anti-commutant ${\sf{A_{G}}}$ of ${\sf{G}}\subseteq{\sf{P}}^{n}$ to be the quotient group ${\sf{P}}^{n}/{\sf{C_{G}}}$. Any $b\in{\sf{P}}^{n}$ can be uniquely decomposed as

[TABLE]

for some $a\in{\sf{A_{G}}}$ and $c\in{\sf{C_{G}}}$.

Any element of the Pauli group (except those proportional to the identity) commutes with exactly half of the other elements in the group. A similar statement is true for any subgroup, and we have the following lemma. In what follows, we use the Iverson bracket notation, where if $x$ is any logical proposition (such as a set inclusion), then

[TABLE]

Lemma 1**.**

For any Pauli $a\in{\sf{P}}^{n}$ and any group ${\sf{G}}\subseteq{\sf{P}}^{n}$,

[TABLE]

Our Pauli estimation procedure varies over states and measurements whose stabilizer groups have elements that cover some set ${\sf{X}}$ of errors. This motivates the following definition.

Definition 1** (Stabilizer coverings).**

A stabilizer covering ${\sf{O}}$ of a set ${\sf{X}}\subseteq{\sf{P}}^{n}$ is a set of stabilizer groups ${\sf{O}}=\{{\sf{S}}_{j}\}$ such that ${\sf{X}}\subseteq\bigcup_{j}{\sf{S}}_{j}$.

For important applications, ${\sf{X}}$ will not necessarily be a group. For example, ${\sf{X}}$ could be the set of Pauli errors with weight $\leq w$, which is not generally a group. The next lemma bounds the size of a stabilizer covering for any set ${\sf{X}}$, and in some cases we can improve over the trivial bound.

Lemma 2**.**

For any set ${\sf{X}}\in{\sf{P}}^{n}$ there exists a stabilizer covering ${\sf{O}}$ of ${\sf{X}}$ with cardinality

[TABLE]

where $\langle{\sf{X}}\rangle$ is the group generated by the elements of ${\sf{X}}$ and ${\sf{S_{X}}}$ is any maximal stabilizer subgroup of $\langle{\sf{X}}\rangle$. In particular, the bound is tight when ${\sf{X}}$ is a group, and $|{\sf{O}}|\leq 2^{n}+1$ unconditionally.

Proof.

For any set ${\sf{X}}$, a simple construction of a stabilizer covering is the set of all two-element stabilizer groups generated by each nontrivial element of ${\sf{X}}$. That is, ${\sf{O}}=\{\langle r\rangle:r\in{\sf{X}}\backslash\{0\}\}$. Clearly this set has $|{\sf{O}}|\leq|{\sf{X}}|$.

Now consider a stabilizer covering of $\langle{\sf{X}}\rangle$, which obviously will also cover ${\sf{X}}$. As a subgroup of ${\sf{P}}^{n}$, the generators can be partitioned into $s$ generators of a stabilizer group ${\sf{S_{X}}}$ and $2k$ generators of the quotient $\langle{\sf{X}}\rangle/{\sf{S_{X}}}$. The quotient comprises $k$ logical qubits, so it is isomorphic to ${\sf{P}}^{k}$ and has $4^{k}$ elements. Then we can construct ${\sf{O}}$ to be the stabilizer groups of a complete set of $|{\sf{O}}|=2^{k}+1$ mutually unbiased bases (MUBs) [64] acting on the $k$ logical qubits, where each MUB is taken in direct sum with ${\sf{S_{X}}}$. The cardinality of this choice for ${\sf{O}}$ is just $2^{k}+1=\sqrt{4^{k}}+1$, as claimed.

The MUBs give a minimal covering when ${\sf{X}}$ is a group, as there are at most $2^{k}-1$ non-identity Pauli elements in a stabilizer group on $k$ qubits, so at least $2^{k}+1$ stabilizer groups are required to cover all $4^{k}-1$ non-identity elements of ${\sf{P}}^{k}$. ∎

II.1 States and measurements

We will use the following states and measurements constructed from stabilizer groups. Recall that $d=2^{n}$. For every $a\in{\sf{P}}^{n}$ and every stabilizer group ${\sf{S}}$,

[TABLE]

is a valid quantum state, known as a stabilizer state. Note that for any $a\in{\sf{P}}^{n}$ and $b\in{\sf{C_{S}}}$, $\rho_{{\sf{S}},a+b}=\rho_{{\sf{S}},a}$ and so we can uniquely label stabilizer states by the error syndromes $e\in{\sf{A_{S}}}$. Moreover,

[TABLE]

is a valid measurement, otherwise known as a syndrome measurement. The normalization is such that each $E_{{\sf{S}},e}$ is an orthogonal projector. If $\dim{\sf{S}}=n$, then $\rho_{{\sf{S}},e}$ is a pure state and $E_{{\sf{S}}}$ is a rank-1 projective measurement.

II.2 Quantum channels

We wish to partially characterize a general linear map $\mathcal{L}:\mathbb{C}^{d\times d}\to\mathbb{C}^{d\times d}$ which has a Kraus operator representation

[TABLE]

We will generally be interested in completely positive (CP) maps, where we can choose $A_{k}=B_{k}$ for all $k$. A specific class of channels that we frequently use are CP maps with a single unitary Kraus operator $U$. For any unitary matrix $U\in U(d)$, we implicitly define the corresponding ideal channel $\mathcal{U}$ that maps an arbitrary matrix $M$ to $UMU^{\dagger}$. Note that the implicit function is many to one, as it maps phase multiples of a unitary matrix to the same channel. For any $a,b\in\sf{P}^{n}$,

[TABLE]

We will partially characterize a general linear map by performing a Pauli twirl of the channel to reduce the effective channel to a Pauli channel, that is, to a linear map with a Kraus representation

[TABLE]

where the $p_{a}$ are Pauli error rates with $p_{a}\geq 0$ for CP maps. For any set $\mathbb{T}\subset U(d)$, we define the $\mathbb{T}$-twirl of a channel $\mathcal{L}$ to be

[TABLE]

In Kraus form, we have

[TABLE]

As we now show, the Pauli-twirled channel $\mathcal{L}^{{\sf{P}}^{n}}$ is a Pauli channel where the error rates are directly related to the Kraus operators of the untwirled channel $\mathcal{L}$. This lemma has previously appeared in [65, Lemma 5.2.4.], but we provide a proof here for completeness.

Lemma 3** (Pauli error rates).**

For a linear map $\mathcal{L}(M)$ as in eq. 9 with Kraus operators $A_{k}=\sum_{a\in{\sf{P}}^{n}}l_{k,a}P_{a}$ and $B_{k}=\sum_{b\in{\sf{P}}^{n}}r_{k,b}P_{b}$, the Pauli error rates $p_{a}$ of the twirled channel $\mathcal{L}^{{\sf{P}}^{n}}$ are given by

[TABLE]

Proof.

From eq. 13,

[TABLE]

As the Pauli matrices are a Hermitian orthogonal basis for $\mathbb{C}^{d\times d}$, we can expand the Kraus operators as $A_{k}=\sum_{a\in{\sf{P}}^{n}}l_{k,a}P_{a}$ and $B_{k}=\sum_{b\in{\sf{P}}^{n}}r_{k,b}P_{b}$, so that

[TABLE]

where in the second line we have used lemma 1 with the fact that ${\sf{C}}_{{\sf{P}}^{n}}=0$. ∎

In the most important case of a CP map, $l_{k,a}=r_{k,a}$ and the error rates are manifestly nonnegative. However, effective negative error rates are possible in the presence of initial correlations, and this would correspond to the inner product in eq. 14 being negative. It is also noteworthy that this definition is independent of the freedom in the Kraus operators of a CP map, since the unitary freedom on the index $k$ above does not change the inner product.

We can simplify the representation of Pauli channels using the vectorization map, a linear map $|\,\cdot\,):X\to|X)$ that acts by stacking the columns of $X$ in a given basis and satisfies the identity $|AXB^{\dagger})=B^{T}\otimes A|X)$. Note that a consequence of this is that the vectorization map can be implemented as $|X)=I\otimes X|\Gamma)$, where $|\Gamma)=\sum_{j}|j\rangle\otimes|j\rangle$ is a maximally entangled state and the $|j\rangle$ form an orthonormal basis. In this notation, the dual vectors map to scalars via the Hilbert-Schmidt inner product, $(A|B)=\mathrm{Tr}(A^{\dagger}B)$. Then we can write the action of a channel $\mathcal{L}$ as a superoperator acting via left multiplication so that \mathcal{L}|\rho)=\bigl{|}\mathcal{L}(\rho)\bigr{)}. Note that we overload our notation by using the same symbol for the abstract channel and the superoperator matrix.

Any Pauli channel can be expressed as a superoperator in the Pauli basis as

[TABLE]

It is convenient to make the following two simplifications. First, we note that $|P_{a})$ is not a normalized vector, so we adopt the convention that

[TABLE]

so that $(a|b)=\tfrac{1}{d}\mathrm{Tr}(P_{a}^{\dagger}P_{b})=\delta_{a,b}$. This notation is unambiguous as long as we distinguish carefully between $a$ and $P_{a}$ from the given context inside the vectorization map. Second, for any two sets of Pauli matrices ${\sf{A}},{\sf{B}}\subseteq{\sf{P}}^{n}$, we define a Walsh-Hadamard transform $\bm{W}_{{\sf{A}},{\sf{B}}}:\mathbb{C}^{|{\sf{B}}|}\to\mathbb{C}^{|{\sf{A}}|}$ whose matrix representation is

[TABLE]

As we prove in lemma 4 below, this transformation is proportional to an isometry for specific groups ${\sf{A}}$ and ${\sf{B}}$.

Given these definitions, the vector of Pauli error rates $\bm{p}$ is related to the superoperator representation in terms of the Pauli channel eigenvalues, or just Pauli eigenvalues for short, which are defined as,

[TABLE]

where we have simplified the notation by using $\bm{W}=\bm{W}_{{\sf{P,P}}}$. After this change of variables we have simply

[TABLE]

With this convention, $f_{0}=\sum_{a}p_{a}\leq 1$ with equality if the channel is trace preserving, and all other channel eigenvalues lie inside the interval $[-f_{0},f_{0}]$ if the channel is completely positive. We are most interested in channels that are close to the identity in the sense that the eigenvalues are in some small, strictly positive interval $[\epsilon,1]$, which is appropriate when the noise is weak. Consequently, for any $h\in{\sf{P}}^{n}$ and any set ${\sf{A}}\subseteq{\sf{P}}^{n}$ we define

[TABLE]

We will omit the subscripts from the vector notation when ${\sf{A}}={\sf{P}}^{n}$.

In what follows, our estimation strategy will be to infer both $\bm{r}$ and, by then doing the inverse transform, $\bm{p}=\bm{W}^{-1}(\bm{1}-\bm{r})$, where we can evaluate the inverse using the following lemma.

Lemma 4**.**

For any subgroups ${\sf{A}},{\sf{B}}\subseteq{\sf{P}}^{n}$, $\bm{W}_{\sf{A,B}}$ satisfies

[TABLE]

In particular, when ${\sf{C_{A}\cap B}}$ is trivial, we have that $\bm{W}_{\sf{A,B}}$ is proportional to an isometry, with $\bm{W}_{\sf{A,B}}^{\dagger}\bm{W}_{\sf{A,B}}=|{\sf{A}}|\Pi_{\sf{B}}$ where $\Pi_{\sf{B}}=\sum_{b\in{\sf{B}}}|b)(b|$ is the projector onto ${\sf{B}}$.

Proof.

From eq. 18,

[TABLE]

where we have used lemma 1 to obtain the third line.

When ${\sf{B}}$ is a group and ${\sf{C_{A}\cap B}}$ is trivial, $b+b^{\prime}\in{\sf{C_{A}}}$ if and only if $b^{\prime}=-b=b$. ∎

Finally, we can evaluate the quality of a reconstructed set of errors $\hat{\bm{p}}$ either in terms of a norm $\|\hat{\bm{p}}-\bm{p}\|$ or a figure of merit on the space of linear maps. The natural norm would be the 1-norm, which is related to the average gate infidelity $r(\mathcal{E})$ and the diamond distance to the identity $\epsilon_{\diamond}(\mathcal{E})$ by

[TABLE]

which holds for Pauli channels [66]. Here we have r(\mathcal{E})=1-\int\mathrm{d}\psi\operatorname{Tr}\bigl{[}\psi\mathcal{E}(\psi)\bigr{]}, where $\mathrm{d}\psi$ is the uniform (Haar) measure over all pure state projectors $\psi$ and $\epsilon_{\diamond}(\mathcal{E})=\sup_{\rho}\frac{1}{2}\|\rho-\mathcal{E}\otimes\mathcal{I}(\rho)\|_{1}$, where the supremum is over all valid quantum states on an enlarged space with an arbitrary-sized ancilla.

However, motivated by physical considerations we will project reconstructions into the nearest point in the set of probability distributions, or potentially to the set of subnormalized probability distributions. The nearest point is not uniquely defined in general according to some norms, such as the 1-norm, and so in some cases we will instead use the 2-norm for which the nearest point is uniquely defined.

II.3 Model assumptions

Our procedure makes use of the following primitives: preparations of stabilizer states, Pauli transformations, and syndrome measurements. An experimental implementation of our procedure will necessarily involve noisy versions of all primitives. For clarity, we always denote the noisy version of a primitive with an overset $\sim$, e.g., $\tilde{A}$ is a noisy implementation of the ideal operation $A$ (whether $A$ is a state preparation, a measurement operator, or a channel).

For ease of analysis, we assume that the noise on the primitives is independent, time-stationary, and Markovian. In particular, we assume that the noisy preparation of a state $\rho$ and a noisy syndrome measurement of a stabilizer group ${\sf{H}}$ are independent of the circuit to be applied, so that they can be written as a density operator $\tilde{\rho}$ and a positive-operator-valued measure (POVM) $\{\tilde{E}_{{\sf{H}},e}:e\in{\sf{E_{H}}}\}$ respectively.

Similarly, we require that the noise in the implementations of a set of twirling channels $\mathcal{T}$ are independent of the particular twirl being implemented and the remainder of the circuit and so can be written as $\tilde{\mathcal{T}}=\mathcal{T}\Lambda$ for some fixed completely positive (CP) map $\Lambda$. We codify this in the following definition.

Definition 2** (GTM noise).**

A noise model is time-stationary if the noisy implementation $\tilde{\mathcal{U}}(t)$ of a gate $\mathcal{U}$ at time $t$ is a linear map that is independent of $t$ and if state preparations and measurements are respectively described by fixed density operators and POVMs. A noise model for the Pauli group is called GTM (gate-independent, time-stationary, Markovian) if it is time-stationary and there exists a completely positive trace-preserving map $\Lambda$ such that $\tilde{\mathcal{P}}=\mathcal{P}\Lambda$ for all $\mathcal{P}\in{\sf{P}}^{n}$.

The above assumption is routinely assumed in analyses of RB for groups $\sf{X}$ that contain complex circuits of multi-qubit gates and can be relaxed with sufficient effort [23, 24, 25]. In contrast, we only make the assumption for groups consisting of tensor products of channels acting on individual qubits.

The other assumption we will make on the noise is that it is sufficiently weak that we can ignore certain algebraic limitations in the use of our specific estimators. In particular, we will assume that our noise and the state preparation and measurements (SPAM) satisfy the following definitions for choices of a parameter $c$ that we will specify later.

Definition 3** (Weak, stable).**

A noise map $\mathcal{L}$ is $\bm{c}$-weak if the Pauli twirl $\mathcal{L}^{{\sf{P}}^{n}}$ is close to the identity channel in the operator norm, $\|\mathcal{I}-\mathcal{L}^{{\sf{P}}^{n}}\|\leq c$. A SPAM parameter $A$ is called $\bm{c}$-stable if $A\geq 1-c$.

The first definition is equivalent to saying that the Pauli twirl $\mathcal{L}^{{\sf{P}}^{n}}$ has Pauli eigenvalues $\bm{f}$ that all lie in the interval $[1-c,1]$. It is likely that this assumption too can be relaxed, but to do so would require different estimators than the ratio estimator that we use below. We note that a simple sufficient condition to ensure that a CP map is $c$-weak is that the Pauli error rate for the identity, $p_{0}$, obeys $p_{0}\geq 1-\tfrac{c}{2}$. This can be verified from the inequality $f_{i}\geq 2p_{0}-1$ [22]. The second definition will become clearer once we define (in eq. 24 below) how the SPAM parameters $A_{j}$ in the definition depend on the noisy state preparations and measurements. Intuitively, it just ensures that the measurements are sufficiently good to give a reasonable signal, so a $c$-stable SPAM parameter is $c$-close to being ideal.

III Estimation procedure and sample complexity

We now specify a procedure for estimating the Pauli eigenvalues of a noisy implementation of the Pauli group ${\sf{P}}^{n}$, where we prepare states using a fixed stabilizer group ${\sf{G}}$ and we make syndrome measurements of a second, possibly different, stabilizer group ${\sf{H}}$. As we will see below, for this procedure to yield useful information about the noise we will choose ${\sf{H}}\subseteq{\sf{G}}$, but we do not require this choice in our proofs. The backbone of the procedure is running a generalized cycle benchmarking sequence of length $m$ [22] and recording the result of a syndrome measurement. One of the core innovations of the present procedure is to consider measurements of general stabilizer groups, which enables us to obtain more information from each experiment. Character benchmarking [21] over the Pauli group, direct RB [67] with the Pauli group as the generator distribution, and cycle benchmarking with no interleaved gate [22] can all be rephrased in terms of the following more general procedure by choosing the measurement to be of a two-element stabilizer group, that is, a group generated by a single Pauli operator.

${\bf RunCB}({\sf{G}},{\sf{H}},m)$

Prepare the approximate stabilizer state $\tilde{\rho}_{{\sf{G}},0}$. 2. 2.

For each $i=0,\ldots,m$, apply a Pauli gate $P_{a_{i}}\in{\sf{P}}^{n}$ uniformly at random. 3. 3.

Perform a syndrome measurement of ${\sf{H}}$ and record the outcome $b\in{\sf{A_{H}}}$ from the noisy POVM element $\tilde{E}_{b}$. 4. 4.

Return $z\in{\sf{A_{H}}}$ such that $b+\sum_{i}a_{i}\equiv z$.

We will use the shorthand ${\bf RunCB}({\sf{G}},m)={\bf RunCB}({\sf{G}},{\sf{G}},m)$ for the special case ${\sf{G}}={\sf{H}}$. The following proposition exactly characterizes the output of ${\bf RunCB}({\sf{G}},{\sf{H}},m)$ under gate-independent and time-stationary noise. We note that ${\bf RunCB}({\sf{G}},{\sf{H}},m)$ can be straightforwardly generalized to include interleaved gates [22]; this will be explored in future work.

Proposition 5**.**

Let $Z=Z({\sf{G}},{\sf{H}},m)$ be the random variable taking values $z\in{\sf{E_{H}}}$ that results from one call to the subroutine ${\bf RunCB}({\sf{G}},{\sf{H}},m)$. Under GTM noise with $\tilde{\mathcal{P}}=\mathcal{P}\Lambda$ for all $P\in{\sf{P}}^{n}$, $Z$ has the probability distribution function

[TABLE]

where the $f_{h}^{m}=(h|\Lambda|h)^{m}$ are the $m$th powers of the Pauli eigenvalues of $\Lambda$, and

[TABLE]

Proof.

By the linearity of the Born rule for density matrices, the probability of a given instantiation $Z=z$ is the sum of all the probabilities of random noisy Pauli gates and measurement outcomes that would lead to that result. This is equivalent to first writing the entire joint distribution over any possible sequence, and then marginalizing over the specific sequences of Pauli gates and measurement outcomes that give the same value of $z$.

Under the assumption of GTM noise, the probability of choosing $a_{0},a_{1},\ldots,a_{m}$ and obtaining the outcome $b$ is

[TABLE]

where for non-commutative products we use the convention $\prod_{i=a}^{b}x_{i}=x_{a}\ldots x_{b}$ (i.e., the lower index is on the left of the product). We can rewrite the probability as

[TABLE]

where we set $a^{\prime}_{0}=a_{0}$ and recursively define $a^{\prime}_{i}=a_{i}+a^{\prime}_{i-1}$ for $i=1,\ldots,m-1$, where $a^{\prime}_{m}=\sum_{i}a_{i}$ (recall that addition is mod 2). Averaging uniformly and independently over the $a_{i}$ is equivalent to averaging uniformly and independently over the $a^{\prime}_{i}$ because ${\sf{P}}^{n}$ is a group. Therefore the marginal probability of $a^{\prime}_{m}\in{\sf{P}}^{n}$ and $b\in{\sf{E_{H}}}$ is

[TABLE]

by eqs. 20 and 10, where in particular we use $(h|\Lambda^{{\sf{P}}^{n}}|h)=(h|\Lambda|h)=f_{h}$ and we define $A_{b,h}=(\tilde{E}_{{\sf{H}},b}|h)(h|\Lambda|\tilde{\rho}_{{\sf{G}},0})$.

Splitting $a^{\prime}_{m}=a+c$ where $a\in{\sf{A_{H}}}$ and $c\in{\sf{C_{H}}}$ and averaging over $c$ gives the marginal probability

[TABLE]

by lemma 1. Now by the double commutant theorem we have ${\sf{C_{C_{H}}}}={\sf{H}}$ for any group ${\sf{H}}\subseteq{\sf{P}}^{n}$, so

[TABLE]

The probability of obtaining the outcome $z\in{\sf{A_{H}}}$ is then

[TABLE]

which completes the proof. ∎

We remark that the coefficients $A_{h}$ appearing in proposition 5 are $1$ for ideal states, measurements, and transformations. Indeed, setting $\tilde{\rho}_{{\sf{G}},0}=\rho_{{\sf{G}},0}$, $\Lambda=1$ and $\tilde{E}_{{\sf{H}},b}=E_{{\sf{H}},b}$ for any fixed outcome $b\in{\sf{A_{H}}}$ in eq. 24 and substituting in eqs. 7 and 8, we have that

[TABLE]

where we have used the normalization from eq. 17. Plugging this into eq. 24, we find that in the ideal case,

[TABLE]

We refer to the coefficients $A_{h}$ as SPAM coefficients (despite having a contribution from edge noise terms) because the coefficients are the only model parameters affected by errors in the state preparations and measurements. We see that the information contained in the likelihood function is maximized when ${\sf{H}}={\sf{G}}$.

Proposition 5 exactly quantifies the output of one call to ${\bf RunCB}({\sf{G}},{\sf{H}},m)$. One could naively obtain many samples of ${\bf RunCB}({\sf{G}},{\sf{H}},m)$ for different values of $m$ and perform a multi-exponential fit to estimate the model parameters $\{A_{h},f_{h}:h\in{\sf{H}}\}$. However, there are exponentially many such parameters, making the fit and a rigorous theoretical treatment nigh impossible. We can resolve this problem by transforming the output using the following lemma to obtain a new random variable with a single exponential decay as in [21]. We note that considering syndrome measurements instead of the two-outcome measurements considered in Refs. [21, 22] allows us to estimate multiple random variables from the same data using the following procedure, which can significantly reduce the sample complexity.

${{\bf V}}({\sf{X}},{\sf{G}},t,m)$: Estimator function for vector of SPAM-dependent Pauli eigenvalues.

Set $\bm{V}\coloneqq 0$. 2. 2.

For $k=1,\ldots,t$, Do

• •

Set $z\coloneqq{\bf RunCB}({\sf{G}},m)$,

• •

Set $\bm{V}\mathrel{+}=\sum_{x\in{\sf{X}}}(-1)^{\langle x,z\rangle}|x)$. 3. 3.

Return $\hat{\bm{V}}=\frac{1}{t}\bm{V}$.

The elements in the output of ${\bf V}({\sf{X}},{\sf{G}},t,m)$ are correlated binomial variables, but their covariance can also be exactly computed using the following lemma.

Lemma 6**.**

Let $\bm{V}$ be the output of ${\bf V}({\sf{X}},{\sf{G}},t,m)$. Under GTM noise, we have

[TABLE]

where $x^{\bot}$ is the component of $x$ in ${\sf{P}}^{n}/{\sf{G}}$ and the $A_{x}$ are as in proposition 5, and, in particular, are independent of $m$.

Proof.

As the calls to ${\bf RunCB}({\sf{G}},m)$ are independent, the expectation value of $\bm{V}$ is

[TABLE]

For any fixed $x\in{\sf{X}}$, substituting in proposition 5 with ${\sf{H}}={\sf{G}}$, noting we ${\bf RunCB}({\sf{G}},m)$), and using lemma 1 gives

[TABLE]

Note that for a stabilizer group ${\sf{G}}$, ${\sf{C_{A_{G}}}}={\sf{P}}^{n}/{\sf{G}}$. The result then follows as ${\sf{G}}$ is a group and ${\sf{P}}^{n}$ splits into ${\sf{G}}\oplus{\sf{P}}^{n}/{\sf{G}}$.

Similarly, for fixed $x,x^{\prime}\in{\sf{X}}$ we have

[TABLE]

which completes the proof. ∎

Since the elements of $\bm{V}$ are bounded, an average of sufficiently many independent samples will converge quickly to the mean. The next proposition provides a simple and conservative tail bound on the probability of estimates obtained using $T$ samples will deviate from their mean by more than $\epsilon$. Specifically, we bound the failure probability of the following procedure, where we allow ${\sf{X}}\subseteq{\sf{G}}$ both for later analysis and to illustrate the fundamental scaling.

Proposition 7**.**

Let ${\sf{G}}$ be a stabilizer group, ${\sf{X}}\subseteq{\sf{G}}$, and $\hat{\bm{V}}$ be the output of $\bm{V}({\sf{X}},{\sf{G}},t,m)$ for some fixed positive integers $t$ and $m$. Then for any $\epsilon>0$,

[TABLE]

Proof.

The proof follows trivially from the union bound and Hoeffding’s inequality [68] where each element of $\hat{\bm{V}}$ is in the interval $[-1,1]$. ∎

Note that the above bound is very loose because, as shown in lemma 6, there are correlations between the elements of $\bm{\hat{V}}$ that are neglected when applying the union bound. However, we still see that the sample complexity above in estimating $A_{j}f_{j}^{m}$ is independent of the Pauli eigenvalues. The precision in estimating the eigenvalues $f_{j}$ with a fixed failure probability will, however, depend on the bare precision $\epsilon$ and on the specific values of $m$ that we choose. In the next section we will analyze a specific procedure and estimator to get bounds on this dependency.

IV Data fitting and error analysis

To understand how the sequence lengths should be chosen to get a precise estimate of a Pauli channel, we must introduce a specific estimation strategy. We use an estimator that samples $Af^{m}$ at exponentially increasing values of $m$, and then we use only the data at the endpoints to obtain a ratio estimator. We also ignore repeated elements from distinct stabilizer groups in a stabilizer covering of a set. Our estimation strategy throws away a lot of data, but it has two advantages. First, it is easy to analyze explicitly, and our proof is a straightforward adaption of the theorem presented in Ref. [37]. Second, the endpoints are where the data are most sensitive, so discarding the intermediate data is not as damaging as one might expect since not all points contribute equally to the variance in the estimator. In what follows we use the notation $V_{i}(g)=(g|V_{i}$ for the $g$th element of $V_{i}$ since the subscript real estate is occupied by the iteration index $i$.

${\bf Ratio}({\sf{O}},{\sf{X}},t)$ Ratio estimator for exponential regression of a stabilizer covering ${\sf{O}}$ of a set ${\sf{X}}$.

Set $\hat{r}_{x}=\texttt{NaN}$ for all $x\in{\sf{X}}$. 2. 2.

Set ${\sf{A}}\coloneqq\emptyset$. 3. 3.

For all ${\sf{G}}\in{\sf{O}}$, Do

1. (a)

Set ${\sf{X}}^{\prime}\coloneqq{\sf{G}}\backslash{\sf{A}}$. 2. (b)

Set ${\sf{A}}\coloneqq{\sf{A}}\cup{\sf{G}}$. 3. (c)

Set $\hat{V}\coloneqq{\bf V}({\sf{X}}^{\prime},{\sf{G}},t,0)$. 4. (d)

Set $m\coloneqq 1$. 5. (e)

While $\exists x\in{\sf{X}}^{\prime}$ such that $\hat{r}_{x}=\texttt{NaN}$, Do

• •

Set $\hat{W}\coloneqq{\bf V}({\sf{X}}^{\prime},{\sf{G}},t,m)$,

• •

For all $x\in{\sf{X}}^{\prime}$ such that $\hat{r}_{x}=\texttt{NaN}$, Do

• –

Set $v\coloneqq\hat{V}(x)$ and $w\coloneqq\hat{W}(x)$.

• –

If $w\leq v/3$ and $w,v>0$, Set \hat{r}_{x}\coloneqq 1-\bigl{(}w/v\bigr{)}^{1/m};

• –

Else If $w\leq 0$ or $v\leq 0$, Set $\hat{r}_{x}\coloneqq 1$.

• •

Set $m\coloneqq 2m$ 4. 4.

Return $\hat{\bm{r}}_{\sf{X}}=\sum_{x\in{\sf{X}}}\hat{r}_{x}|x)$.

This ratio estimator guarantees a pointwise multiplicative precision estimate of the true vector of infidelities $\bm{r}_{\sf{X}}$ (recall eq. 21), at least when certain mild assumptions hold that place us close to the regime of interest. This is where the notions of $c$-weak and $c$-stable from section II.3 come into play. We will assume that the noise map is $\tfrac{1}{2}$-weak, which in particular implies that $f_{h}\geq\tfrac{1}{2}$, or equivalently, that $r_{h}\leq\tfrac{1}{2}$. We will also assume that each SPAM parameter $A_{h}$ from eq. 24 is $\tfrac{1}{2}$-stable, so that $A_{h}\geq\tfrac{1}{2}$. However, the only thing that our proof really requires is that both $A$ and $f$ are bounded from below by a positive constant, and the specific choice of $\tfrac{1}{2}$ is motivated only by convenience.

We note that an equivalent theorem to the following can be proven if we remove the $\tfrac{1}{2}$-stable assumption and instead assume that $\epsilon$ is chosen to give a relative precision proportional to $\min_{h}A_{h}f_{h}$ at each step instead of only a constant. However, this makes the statements about sample complexity dependent on the SPAM parameters $A_{h}$. This obfuscates the actual sample complexity of the procedure, so we prefer to add the stability assumption.

To understand how the set of sequence lengths $\kappa$ should be chosen to get a precise estimate of a Pauli channel, we must introduce the spectral gap of the channel, or of a subset of the channel corresponding to the Paulis in the set ${\sf{X}}$. Since a Pauli channel is already diagonal in the Pauli basis, as in Eq. (20), and every trace-preserving channel has a largest absolute eigenvalue of 1, the spectral gap $\Delta_{{\sf{X}}}=\Delta_{{\sf{X}}}(\Lambda)$ of the Pauli channel $\Lambda$ over the set ${\sf{X}}$ is given by

[TABLE]

In the most interesting regime, $\Lambda$ is $\epsilon$-weak for some $\epsilon<1$, and the absolute value is not necessary because then all eigenvalues are positive. In that case, we have $\Delta_{\sf{X}}=\min_{x\in{\sf{X}}}r_{x}$. When ${\sf{X}}={\sf{P}}^{n}$, this is exactly the spectral gap of the superoperator for the channel. Note that $\Delta_{\sf{X}}\geq\Delta_{{\sf{P}}^{n}}$.

As we will now show, one can completely learn all the Pauli eigenvalues in ${\sf{X}}$ with relative precision $\epsilon$ using a set $\kappa$ of sequence lengths with |\kappa|=O\bigl{(}\log\tfrac{1}{\Delta_{\sf{X}}}\bigr{)} and m_{\max}=\max\kappa=O\bigl{(}\tfrac{1}{\Delta_{\sf{X}}}\bigr{)}.

Proposition 8**.**

Let ${\sf{X}}\subseteq{\sf{P}}^{n}$ and ${\sf{O}}$ be a stabilizer covering of ${\sf{X}}$. For any sufficiently small $\epsilon,\delta>0$ and assuming $\tfrac{1}{2}$-weak, $\tfrac{1}{2}$-stable, GTM noise, the following holds with probability $1-\delta$. Running ${\bf Ratio}({\sf{O}},{\sf{X}},t)$ with t=\frac{2}{\epsilon^{2}}\log\bigl{(}\frac{2|{\sf{X}}||\kappa|}{\delta}\bigr{)} uses a set $\kappa$ of at most O\bigl{(}\log\tfrac{1}{\Delta_{\sf{X}}}\bigr{)} sequence lengths $m$ with m_{\max}=O\bigl{(}\tfrac{1}{\Delta_{\sf{X}}}\bigr{)}. Moreover, the output $\hat{\bm{r}}_{\sf{X}}$ satisfies

[TABLE]

Proof.

The following proof of multiplicative precision is a modification of a similar result proven recently in Ref. [37].

We begin by considering each iteration of the loop in step 3. Let $t$ be a fixed integer, ${\sf{G}}\in{\sf{O}}$ be a fixed group, and ${\sf{X}}^{\prime}$ be the value in step 3a in the corresponding iteration of ${\bf Ratio}({\sf{O}},{\sf{X}},t)$. The inner loop (e) yields estimators $\hat{\bm{V}}_{m}={\bf V}({\sf{X}}^{\prime},{\sf{G}},t,m)$ for each $m\in\kappa=\{0,1,2,4,\ldots,m_{\rm max}\}$ for some yet-unspecified terminating sequence length $m_{\rm max}$. Note that we do not store all these intermediate estimators in the procedure as we only use them to prove correctness.

By proposition 7 and the union bound, for any fixed $\epsilon>0$,

[TABLE]

with probability at least $1-\delta_{{\sf{X}}^{\prime}}$ where $\delta_{{\sf{X}}^{\prime}}=2|\kappa||{\sf{X}}^{\prime}|\exp(-t\epsilon^{2}/2)$. To bound the terminating value $m_{\rm max}$ for each ${\sf{G}}\in{\sf{O}}$ and hence prove the stated sample complexity, fix $g\in{\sf{G}}$, let $m$ be the corresponding value in the procedure when $\hat{r}_{x}$ is assigned and $w=\hat{V}_{m}({\sf{X^{\prime}}},{\sf{G}},t,m)$ and $u=\hat{V}_{m}({\sf{X^{\prime}}},{\sf{G}},t,\lfloor m/2\rfloor)$. Then $u,v,w$ satisfy the following inequalities

[TABLE]

For simplicity, we now use the notation $f=f_{g}$ and $A=A_{g}$. We also assume that $m>1$ as otherwise the proof is trivial. Let us define the deviations of $u,v,w$ from their mean,

[TABLE]

By proposition 7, $|\epsilon_{x}|\leq\epsilon$ for $x=u,v,w$. Using the fact that $A\geq 1/2$ for $\tfrac{1}{2}$-stable noise,

[TABLE]

Similarly,

[TABLE]

Substituting eqs. 27 and 28 into eq. 26, we have

[TABLE]

We have that $r=1-f\geq\Delta_{{\sf{X}}}$, and note that $f^{m}=(1-r)^{m}$ will satisfy eq. 30 for fixed $\epsilon$ if and only if $m=\Theta(1/r)$.

To show that $|\hat{r}-r|\leq O(\epsilon)r$, we have to analyze the accuracy of the estimator $(w/v)^{1/m}$. Using eqs. 27 and 29, we have

[TABLE]

We can factor out the $f^{m}$ (since $f>0$ as we have assumed the noise is $1/2$-weak) and use the lower bound from eq. 30 to find

[TABLE]

For sufficiently small $\epsilon$, this implies

[TABLE]

Now we Taylor expand around $\epsilon=0$ for sufficiently small $\epsilon$ and use $\tfrac{1}{m}=\Theta(r)$ to obtain the stated accuracy.

Now let $g\in{\sf{X}}$ be the element that requires the largest value of $m_{\rm max}$. Then the set of sequence lengths $\kappa$ used to estimate all the $f_{a}$ to within the stated precision with probability at least $1-\delta_{{\sf{X}}^{\prime}}$ is $\kappa=\{2^{i}:i=0,\ldots,\ell=\log_{2}(m_{\rm max})\}$. As $m_{\rm max}=\Theta(1/r_{g})$, we have |\kappa|-1=\log_{2}m_{\rm max}=O\bigl{(}\log\frac{1}{\Delta_{\sf{X}}}\bigr{)}, which is well defined as the noise is $1/2$-weak by assumption. The total failure probability is then at most $\delta$ where $\delta=2|\kappa||{\sf{X}}|\exp(-2t\epsilon^{2}/2)$ by the union bound. Rearranging gives the stated sample complexity. ∎

V Reconstructing error rates for a group

We have seen how we can estimate each of the Pauli eigenvalues $f_{a}$ of a Pauli channel for all $a\in{\sf{X}}$ with a number of measurements that scales like O\bigl{(}|{\sf{O}}|\log|{\sf{X}}|\bigr{)} for any stabilizer covering ${\sf{O}}$ of ${\sf{X}}$, assuming that $\epsilon$, $\delta$ and the channel spectral gap $\Delta_{\sf{X}}$ are fixed. When ${\sf{X}}$ is a group, we can choose ${\sf{O}}$ such that $|{\sf{O}}|\leq\sqrt{|{\sf{X}}|}+1$ using the construction in lemma 2, so the total number of samples is at most $O(\sqrt{|{\sf{X}}|}\log|{\sf{X}}|)$ (or even less if ${\sf{X}}$ has a nontrivial stabilizer subgroup). Moreover, the precision of these estimates is multiplicative. We now show how the estimates of Pauli eigenvalues can be used to estimate Pauli error rates.

Suppose we estimate all the Pauli eigenvalues of some group ${\sf{G}}$ (which could be the full Pauli group). Then eq. 19 becomes

[TABLE]

where we use the shorthand $\bm{v}_{{\sf{G}}}$ to denote the vector with entries $\{v_{g}:g\in{\sf{G}}\}$. The columns of $\bm{W}_{{\sf{G}},{\sf{P}}^{n}}$ for Paulis that differ by an element of ${\sf{C_{G}}}$ are identical and so the corresponding Pauli error rates cannot be distinguished using only the estimated Pauli eigenvalues. That is, we can only reconstruct a marginal or coarse-grained probability distribution $\bm{p}_{{\sf{A_{G}}}}$ over ${\sf{A_{G}}}$ with marginal probabilities

[TABLE]

for $a\in{\sf{A_{G}}}$ via

[TABLE]

When ${\sf{G}}$ is the group generated by the stabilizers and logical operators of an error-correcting code, the marginal probabilities give exactly the distribution of logical errors for each syndrome. This marginal distribution can therefore be used to construct the maximum likelihood decoder for the actual noise afflicting a device.

In practice, inverting eq. 36 to find the marginal error rates might be solved using, for example, a non-negative least squares solver, being careful to take the nontrivial covariance into account as well when estimating the errors. However, to prove concrete theorems, we will use the estimator

[TABLE]

using lemma 4. Here we use the notation $[\bm{v}]_{\mathcal{D}}$ to mean taking a $k$-dimensional vector $\bm{v}$ and projecting it to the nearest point (according to a Euclidean metric) in the $k$-point simplex, which we denote $\mathcal{D}$, with the dimension being inferred from the context. The informal statement in result 1 follows as $\|\bm{r}_{\sf{G}}\|_{\infty}\leq 2(1-p_{0})$ by [22, Lemma 4].

Proposition 9**.**

Let ${\sf{G}}$ be a group of Pauli matrices with a maximal stabilizer subgroup ${\sf{S}}$. Then, under the conditions of proposition 8, an estimator $\hat{\bm{p}}_{\sf{A_{G}}}$ of $\bm{p}_{\sf{A_{G}}}$ satisfying

[TABLE]

with probability at least $1-\delta$ can be obtained using t=\frac{2}{\epsilon^{2}}\log\bigl{(}\frac{2|\kappa||{\sf{G}}|}{\delta}\bigr{)} measurements per round for $|\kappa|(\sqrt{|{\sf{G}}/{\sf{S}}|}+1)$ rounds, where |\kappa|=O\bigl{(}\log\tfrac{1}{\Delta_{\sf{G}}}\bigr{)} and m_{\rm max}=O\bigl{(}\tfrac{1}{\Delta_{\sf{G}}}\bigr{)} is the longest sequence length used.

Proof.

Suppose we have an estimator $\hat{\bm{r}}_{\sf{G}}$ satisfying

[TABLE]

First note that projecting to the simplex can only decrease the error, that is,

[TABLE]

By lemma 4 and using ${\sf{C_{A_{G}}}}\cap{\sf{G}}=0$, $|{\sf{G}}|^{-1/2}\bm{W}_{{\sf{A_{G}}},{\sf{G}}}$ is a unitary transformation on the space $\bm{v}_{\sf{G}}$, so, using a standard norm equivalence,

[TABLE]

Under the conditions of proposition 8, $\epsilon_{0}\leq O(\epsilon)\|\bm{r}_{{\sf{G}}}\|_{\infty}$. The bound on the number of rounds and the total number of measurements follows from the bound on the size of a stabilizer covering of ${\sf{G}}$ from lemma 2. ∎

In the special case that the group ${\sf{G}}={\sf{P}}^{n}$, proposition 9 provides a direct method for estimating the Pauli projection of a noise channel and gives pointwise precision to within $O(\epsilon)(1-p_{0})$, where $p_{0}$ is the probability of no error. Because of its frequent use in applications, it may be desirable to state the following corollary of proposition 9 in terms of the diamond distance.

Corollary 10**.**

Let \Lambda=\mathcal{E}^{{\sf{P}}^{n}}\bigr{|}_{\sf{G}} be the restriction to a subgroup ${\sf{G}}\subseteq{\sf{P}}^{n}$ of the Pauli projection of a channel $\mathcal{E}$. For any sufficiently small $\epsilon,\delta\geq 0$, an estimate $\hat{\Lambda}$ of any $\tfrac{1}{2}$-weak, $\tfrac{1}{2}$-stable, GTM noise model $\Lambda$ can be reconstructed with

[TABLE]

using O\Bigl{[}\tfrac{|\kappa||{\sf{G}}|^{1/2}}{\epsilon^{2}}\log\bigl{(}\tfrac{|\kappa||{\sf{G}}|}{\delta}\bigr{)}\Bigr{]} samples with |\kappa|=O\bigl{(}\log(1/\Delta)\bigr{)} and sequence lengths at most $O(1/\Delta)$.

Proof.

This follows from eq. 38 in proposition 9 using the norm equivalence

[TABLE]

with $\|\hat{\bm{r}}_{\sf{G}}-\bm{r}_{\sf{G}}\|_{2}\leq O(\epsilon)\|\bm{r}_{\sf{G}}\|_{2}$ from proposition 8 and $\|\bm{r}_{\sf{G}}\|_{2}=\|\mathcal{I}_{{\sf{G}}}-\Lambda\|_{F}$, where $\mathcal{I}_{{\sf{G}}}$ is the restriction of the identity channel to ${\sf{G}}$. ∎

VI Reconstructing a subset of errors

We have shown that all Pauli error rates for a marginal model can be reliably reconstructed with substantially fewer resources than might have been anticipated. However, reconstructing an exponential number of probabilities is manifestly inefficient. We now show how to reconstruct the dominant error rates in an approximately sparse noise model. The motivating application is that errors can be divided into a “background” noise process due to, e.g., independent dephasing and amplitude damping errors, and some additional errors that arise from unknown couplings. The goal is to learn these couplings and either determine the physical mechanism causing them and eliminate them via re-engineering the device or to introduce compensating pulses.

To set the scale for the number of dominant error rates, consider a background local depolarizing process on $n$ qubits where independent single-qubit errors occur with probability $p$. Then a weight $w$ error (that is, an error that acts nontrivially on $w$ qubits) occurs with probability $\binom{n}{w}(1-p)^{n-w}p^{w}$. For $np\ll 1$ (where, e.g., 50 qubits with $p=0.001$ is approximately the current state of the art), the probability of an error with more than weight 1 is approximately $(np)^{2}/2$. Now suppose an additional noise process is mixed in that, with probability $p^{\prime}$ uniformly at random applies one of $t$ additional errors with no a priori structure, representing, e.g., unknown correlated errors. Then $3n+t$ errors account for almost all of the errors in the device.

We first show how to use the ${\bf Ratio}$ subroutine to estimate the probabilities of all errors within an arbitrary set ${\sf{E}}$. Formally, we will prove pointwise precision to within $O(\epsilon)(1-p_{0})$, where $p_{0}$ is the probability of no error. By eq. 22 we have $1-p_{0}=(1+1/d)r$ where $r$ is the so-called average error rate. We will then present a divide-and-conquer routine to efficiently identify a set ${\sf{E}}$ of fixed size that can be used to identify correlated errors. The informal version of the following theorem stated as result 2 in the introduction follows from using the trivial stabilizer covering ${\sf{O}}$ of ${\sf{X}}$ with $|{\sf{O}}|=|{\sf{X}}|$.

Theorem 11**.**

For any set ${\sf{E}}\subseteq{\sf{P}}^{n}$, an estimator $\hat{\bm{p}}_{\sf{E}}$ of $\bm{p}_{\sf{E}}$ satisfying

[TABLE]

with probability at least $1-\delta$ can be obtained from the output of ${\bf Ratio}({\sf{O}},{\sf{X}},t)$ using a random set ${\sf{X}}\subseteq{\sf{P}}^{n}$ with $|{\sf{X}}|=\tfrac{1}{\epsilon^{2}}\log(4|{\sf{E}}|/\delta)$, any stabilizer covering ${\sf{O}}$ of ${\sf{X}}$, and $t=\tfrac{1}{\epsilon^{2}}\log(4|{\sf{X}}||\kappa|/\delta)$ measurements per round for $|{\sf{O}}|$ rounds.

Proof.

By eqs. 19, 4 and 1, we have

[TABLE]

for any $a\in{\sf{P}}^{n}$. Let ${\sf{X}}\subseteq{\sf{P}}^{n}$ be a set with $|{\sf{X}}|$ elements sampled independently and uniformly at random without replacement and

[TABLE]

By [22, Lemma 4], we have $r_{b}\in[0,2(1-p_{0})]$ for all $b\in{\sf{P}}^{n}$. Applying Serfling’s inequality [69] to the first term gives

[TABLE]

Now suppose we have estimates $\hat{r}_{b}$ of $r_{b}$ for all $b\in{\sf{X}}$ and let

[TABLE]

By the triangle inequality and eq. 39, we have that with probability at least $1-\delta/2$,

[TABLE]

By proposition 8, $|\hat{r}_{b}-r_{b}|\leq O(\epsilon)r_{b}$ with probability at least $1-\delta/2$ using t=\frac{2}{\epsilon^{2}}\log\bigl{(}\frac{4|{\sf{X}}||\kappa|}{\delta}\bigr{)} measurements per round. Therefore by the union bound, we have that with probability at least $1-\delta$,

[TABLE]

for any $a\in{\sf{P}}^{n}$. The final result holds for any set ${\sf{E}}\subseteq{\sf{P}}^{n}$ by the union bound, redefining $\delta\to\delta/|{\sf{E}}|$. ∎

We now provide a search heuristic to identify sets ${\sf{E}}$ of interest using a function ${\bf Select}$ that returns the indices (as Pauli operators) of entries in a vector $\bm{p}$ of probabilities satisfying a desired condition. We use three additional subroutines in what follows. ${\bf Choose}({\sf{A}},b)$ returns a set of $\min\{b,|{\sf{A}}|\}$ elements of ${\sf{A}}$ chosen uniformly at random without replacement. ${\bf Cover}({\sf{F}})$ returns any valid stabilizer covering of a set ${\sf{F}}\subseteq{\sf{P}}^{n}$. ${\bf Ratio}$ is as described in section IV. For any ${\sf{A}}\subset\mathbb{N}$, let ${\sf{P}}_{\sf{A}}\cong{\sf{P}}^{|{\sf{A}}|}$ be the set of Pauli operators that act trivially on all qubits not in ${\sf{A}}$ and the support of a set $\mathrm{supp}({\sf{X}})\subseteq{\sf{P}}^{n}$ be the set of qubits that are acted on nontrivially by some element of ${\sf{X}}$.

${\bf TreeReconstruction}(t,{\bf Select},u,n)$ Reconstruct an $n$-qubit error model using a function ${\bf Select}$ to select error terms at each stage of the iteration.

Set $E^{\prime}\coloneqq\{{\sf{P}}_{\{j\}}:j\in\mathbb{Z}_{n}\}$. 2. 2.

For $j=0,1,\ldots,\lceil\log_{2}n\rceil$, Do

1. (a)

Set $E\coloneqq E^{\prime}$. 2. (b)

Set $F\coloneqq\emptyset$. 3. (c)

For each $e\in E$, Do

• •

Set $F\coloneqq F\cup{\bf Choose}({\sf{P}}_{\mathrm{supp}(e)},u)$. 4. (d)

Set {\sf{O}}\coloneqq{\bf Cover}\bigl{(}\cup_{f\in F}f\bigr{)}. 5. (e)

Set $\bm{r}\coloneqq{\bf Ratio}({\sf{O}},\cup_{f\in F}f,t)$. 6. (f)

Set $E^{\prime}\coloneqq\emptyset$. 7. (g)

For $k=0,1,\ldots,|E|-1$, Do

• •

Set $\displaystyle\bm{p}\coloneqq\sum_{e\in E_{k}}\mathbb{E}_{x\in F_{k}}(-1)^{\langle e,x\rangle}r_{x}|e)$.

• •

Set $E_{k}\coloneqq{\bf Select}(\bm{p})$.

• •

If $k$ is odd, Set $E^{\prime}\coloneqq E^{\prime}\cup\{E_{k-1}\otimes E_{k}\}$. 8. (h)

If $|E|$ is odd, set $E^{\prime}\coloneqq E^{\prime}\cup\{E_{|E|-1}\}$. 3. 3.

Return $\displaystyle\bm{p}\coloneqq\sum_{e\in E_{0}}\mathbb{E}_{x\in F_{0}}(-1)^{\langle e,x\rangle}r_{x}|e)$.

Theorem 12**.**

Let ${\bf Select}$ be any function, $u$ be any positive integer and ${\bf Cover}({\sf{X}})$ return the trivial cover from lemma 2. Then for sufficiently large $t$, ${\bf TreeReconstruction}(t,{\bf Select},u,n)$ uses O\bigl{(}utn\log(1/\Delta)\bigr{)} measurements.

Proof.

At step 2(c) in the $j$th iteration, we are guaranteed that each of the at most $2^{j}$ subsets of $F$ each contain at most $u$ terms. Therefore $|\cup_{f\in{\sf{F}}}f|=|{\sf{O}}|\leq 2^{\lceil\log_{2}(n)\rceil-j}u$ in each iteration. For sufficiently large $t$, ${\bf Ratio}({\sf{O}},\cup_{f\in F}f,t)$ uses $t$ measurements for each of $|{\sf{O}}||\kappa|$ rounds where |\kappa|=O\bigl{(}\log\tfrac{1}{\Delta}\bigr{)}. Therefore the $j$th iteration uses O\bigl{(}2^{\lceil\log_{2}(n)\rceil-j}ut\log\tfrac{1}{\Delta}\bigr{)} measurements. Summing over $j=0,\ldots,\lceil\log_{2}n\rceil$ completes the proof. ∎

A natural choice of ${\bf Select}$ is the function that returns the indices of the $s$ largest elements of the input $\bm{p}$. With this function, ${\bf TreeReconstruction}$ will heuristically return a set ${\sf{E}}$ that is a good $s$-sparse approximation of the Pauli channel. We note that we do not attempt to prove that the resulting set is a best $s$-sparse approximation, that is, that it is a maximizer (or near-maximizer) of $\|\bm{p}_{\sf{E}}\|_{1}$ over all ${\sf{E}}$ of size $s$. It would be possible to prove an upper bound on how close it gets to an optimal approximation by estimating the precision of each error in the set to high precision and backtracking if $\sum_{e\in E_{k}}\bm{p}_{e}$ is below some threshold value at any iteration. However, a rigorous statement along these lines would likely only hold with impractical amounts of resources or a tighter analysis of the distribution of the output of ${\bf Ratio}$. We leave the question of provable approximation ratios to future work.

We also note that a better run-time may be achieved in this setting using methods for reconstructing a sparse signal from few measurements [70, 71, 72]. Specifically, the present variant can be cast in a similar form to sparse Fourier and Hadamard transforms [73, 74, 75, 76]. However, it is an open question if the methods in those references can be adapted to the symplectic structure of the Pauli group to achieve near-optimal sparse reconstruction. Moreover, while finding such a set is desirable, we do not require any property of the set ${\sf{E}}$ for the reconstruction of $\bm{p}_{\sf{E}}$ to be pointwise convergent. Thus, when the heuristic converges to a sparse set with large measure it is a certificate of correctness since each estimate in the set is pointwise accurate.

Another interesting choice of ${\bf Select}$ is a function that returns the indices of errors that do not fit a background model. The most natural background model is for independent single-qubit errors, which can be directly estimated using theorem 11.

VII Bounded degree graphical models

In this section we show how the reconstruction procedures above for complete Pauli channels can be used to learn a Pauli channel on $n$ qubits with bounded degree correlations. It will be convenient to adopt a slightly different notation than what was used previously that treats the probability distribution $p$ over Pauli errors as a function of a random variable rather than a vector. That is, we write $p(\bm{x})$ for the probability of the string of Paulis $\bm{x}$, $p(\bm{x}|\bm{y})$ for the conditional probability of $\bm{x}$ given $\bm{y}$, and so on. In particular, if $p(\bm{x},\bm{y})$ is a joint distribution, then we write $p(\bm{x})$ for the marginal. We will now review some concepts from the field of probabilistic graphical models; see Ref. [77] for an introduction.

Probabilistic graphical models are, at the most basic level, probability distributions over collections of random variables where the dependency structure between the variables is specified by a bipartite graph called a factor graph. Consider a collection of random variables with $x_{j}$ denoting the $j$th random variable. The variable $\bm{x}$ will be a random $n$-qubit Pauli, with $x_{j}$ being the $j$th tensor factor (i.e., the single-qubit Pauli acting on qubit $j$), but we leave the discussion more general for the moment. Suppose that there exist factors $C_{k}$ that are subsets of the random variables such that the joint probability distribution $p$ over all the $x_{j}$ obeys

[TABLE]

where $\phi_{k}:\bm{x}_{C_{k}}\to\mathbb{R}^{+}$ are strictly positive functions called factor potentials supported on the factors; the argument of $\phi_{k}$ is shorthand notation for the subset of variables contained in a factor, $\bm{x}_{C}=\{x_{j}:j\in C\}$; and $Z$ is a normalization constant called the partition function. (Note that some authors use the term factor potential to refer instead to $\log\phi_{k}$.) Then the factor graph for this model is the bipartite graph with one set of nodes labeled by the variable labels $j$ and the other set of nodes labeled by the factors $C_{k}$, and an edge between $j$ and $C_{k}$ if and only if $x_{j}\in C_{k}$. Such a factorized strictly positive probability distribution is known as a Gibbs random field since the factor potentials play the role of exponentiated energy potentials in the Gibbs distribution of a statistical mechanical model. Note that there always exists a Gibbs random field for any probability distribution, although the factor graph may be trivial, that is, contain only one factor.

The structure of a nontrivial factor graph implies that not all variables in $\bm{x}$ can be arbitrarily correlated. In particular, the correlation is controlled by the Markov blanket $\bm{x}_{\partial S}$ of a subset of variables $\bm{x}_{S}$, where the set $\partial S$ is defined as

[TABLE]

That is, it is the set of variables that are adjacent to $S$ via a factor node in the factor graph, but excluding $S$ itself. We will refer to both $\partial S$ and $\bm{x}_{\partial S}$ as the Markov blankets of both $S$ and $\bm{x}_{S}$ with context resolving any ambiguity. The closure $\bar{S}$ is the union of a set and its Markov blanket, $\bar{S}=S\cup\partial S$.

The celebrated Hammersley-Clifford theorem [39, 78] describes how the Markov blanket controls the dependency between variables in a Gibbs random field. The variables $\bm{x}_{A}$ and $\bm{x}_{A^{c}}$, where $A^{c}$ is the complement of $A$, obey the following conditional independence relation in terms of the Markov blanket of $A$,

[TABLE]

This relation, known as the local Markov property, states that the variables in a set $A$ have only bounded correlation in the sense that they are conditionally independent of any other variables outside of their Markov blanket $\partial A$.

To illustrate some of these concepts, consider the example factor graph given in fig. 1. Here the factors are the neighborhoods of the black squares, so for example $C_{1}=\{x_{1},x_{2},x_{5},x_{6}\}$. The Markov blanket of $x_{1}$ is given by $\{x_{2},x_{5},x_{6}\}$, and the Markov blanket of $\{x_{1},x_{5}\}$ is $\{x_{2},x_{6},x_{9},x_{10}\}$. The closure of $\{x_{1},x_{2}\}$ is $\{x_{1},x_{2},x_{3},x_{5},x_{6},x_{7}\}$. The variable $x_{1}$ is conditionally independent of every variable except $\{x_{2},x_{5},x_{6}\}$ and the set of random variables $\{x_{6},x_{7},x_{10},x_{11}\}$ depends conditionally on every other variable.

We will provide an estimate of a Gibbs random field by estimating the factor potentials. We will not provide estimates for the partition function because estimating general partition functions is believed to be computationally hard [79, 80], and so properly normalizing our probability estimate will in general be intractable. However, there are still several cases of practical interest where this is provably not the case. For example, when the factor graph is a tree, and when $p(\bm{0})$ is not too small to be estimated via sampling. Even without being able to efficiently compute the normalization exactly, either heuristics could be used, or ratios of probabilities can be used unconditionally. Ratios of probabilities are all that is needed in, for example, most Monte Carlo methods. However, as we now show, the error on the renormalized distribution can be bounded by the error on the factor potentials without needing to obtain the partition functions.

Lemma 13**.**

Let $p$ and $q$ be Gibbs random fields with the same factor graph and factorizations $p(\bm{x})=\tfrac{1}{Z_{p}}\prod_{C}\phi_{p,C}(\bm{x}_{C})$ and $q(\bm{x})=\tfrac{1}{Z_{q}}\prod_{C}\phi_{q,C}(\bm{x}_{C})$ respectively. Then

[TABLE]

Proof.

We begin by recalling Pinsker’s inequality [81], one form of which is given by

[TABLE]

where $D(p\|q)$ is the relative entropy, defined by

[TABLE]

Thus it suffices to bound the symmetric relative entropy, which simplifies to

[TABLE]

for any probability distributions $p$ and $q$. Substituting the definitions of the Gibbs random fields into the logarithms, we have

[TABLE]

where to obtain the second line we have used the fact that $Z_{p}$ and $Z_{q}$ are independent of $\bm{x}$ and that $\sum_{\bm{x}}p(x)-q(x)=0$, which holds as $p$ and $q$ are probability distributions. By applying Holder’s inequality separately for each $C$, we obtain

[TABLE]

where the second line follows from the bounded dependence of the factors. Substituting the bound on the symmetric relative entropy into eq. 44 and canceling the common factor of $\|p-q\|_{1}$ gives the desired result. ∎

Given a probability distribution $p(\bm{x})$ that is a Gibbs random field, the choice of factor potentials and partition function are not unique. However, the Hammersley-Clifford theorem also gives an explicit description of a certain canonical choice for the factor potentials as a function of the marginal probability distributions on the factors and their Markov blankets, as well as an explicit partition function. It is this description that we will use for our estimators below. We use the parameterization from Abbeel et al. [58] that applies directly to factor graphs (as opposed to less general formulations like Bayesian networks or Markov random fields) and makes explicit use of the local Markov property. We first fix a fiducial assignment to $\bm{x}$ which we call $\bm{0}$. This is simply a reference value and can be any arbitrary fixed assignment to $\bm{x}$; however we will see below that a convenient choice for our purposes will be to choose as a reference value the “identity Pauli” outcome for each variable. Next we will augment the factors in our graph and consider $2^{C_{k}}$, the set of all subsets of $C_{k}$. The union of all of these, minus the empty set, defines a new, larger collection of factors ${\sf{C}}^{*}$, given explicitly by

[TABLE]

Next, given a factor $C^{*}_{k}\in{\sf{C}}^{*}$, define the canonical factor potential for $C^{*}_{k}$ by

[TABLE]

Let $C_{j}$ be a factor such that $C^{*}_{k}\in 2^{C_{j}}$. Then every $S\subset C^{*}_{k}$ is also a subset of $C_{j}$ and so every probability on the right-hand-side of eq. 46 depends only upon $\bar{C}_{j}$, the closure of $C_{j}$ with its Markov blanket. That is, for all $S$ in the sum we have

[TABLE]

Therefore each $\phi_{k}$ is determined by $p(\bm{x}_{\bar{C}_{j}})$ for some factor $C_{j}$. If each $\bar{C}_{j}$ has constant size, each $\log\phi_{k}$ is completely determined by a small amount of data.

The Gibbs random field can now be expressed in terms of the canonical factor potentials as

[TABLE]

We note that $1/p(\bm{0})$ plays the role of the partition function $Z$ and can be independently estimated via theorem 11.

The error bounds below are derived assuming that an independent estimate is obtained for each of the raw factors in the canonical factor decomposition. This adds a constant factor overhead to the overall sample complexity (when the degree of the factor graph is bounded). This overhead could be avoided by using estimates from a covering set of marginals and then directly computing estimates of the marginals on the subfactors, that is, on the subsets of the factors. However, this might bias the overall estimate somewhat and would complicate the analysis presented below, so we leave an understanding of this more efficient protocol to future work.

Importantly, the estimates for the individual raw factors will not generally be self-consistent, and expressing this requires careful notation. Let $\hat{p}_{S}(\bm{x}_{S})$ denote an empirical estimate of the marginal distribution $p(\bm{x}_{S})$ for a set $S$, which in our case will be obtained from proposition 9. Then for two sets $A,B$ such that $S\subset A$, $S\subset B$ and $A\neq B$, the two marginal distributions obtained over $S$ from the empirical estimates do not need to agree, that is,

[TABLE]

Indeed, they will generically not agree in the presence of sampling errors. Moreover, neither will generally agree with the corresponding marginal of the global reconstruction $q$ obtained by substituting the $\hat{p}_{S}$ into eq. 46. Thus, the canonical factor potentials are essential to round each of the local empirical marginals into a global coherent probability distribution with self-consistent local marginals.

As a final difficulty, we cannot use an estimate of $p(\bm{0})$ as the partition function to exactly normalize the empirically reconstructed distribution. Fortunately, as shown in lemma 13, we do not need to know the exact value of the partition function for the empirical distribution to bound the error between our empirical estimate and the true distribution.

Our estimator of the global Gibbs distribution is as follows. We will simply use the expressions in eqs. 46 and 48 where empirical estimates of the local conditional probabilities are used in place of the exact distributions. These can be obtained by using proposition 9 to estimate the complete marginal on the closure of each factor, $\bar{C}$. When an estimate $q$ is obtained in this way, we say that $q$ is a canonical estimator of $p$.

Definition 4** (Canonical estimator).**

Given a Gibbs random field $p$ with known factor graph ${\sf{C}}$, a distribution q=q\bigl{(}\{\hat{p}_{\bar{C}}:C\in{\sf{C}}^{\star}\}\bigr{)} obtained by substituting empirical estimates $\{\hat{p}_{\bar{C}}:C\in{\sf{C}}^{\star}\}$ into eqs. 46 and 48 is called a canonical estimator of $p$.

We now bound the error on our empirical estimate assuming it has been properly normalized. In order to obtain a nontrivial bound, we require that the empirical estimates of the marginal probabilities used to estimate the factor potentials are sufficiently close to their true values. We also require that the true and estimated values are strictly positive, a point that we will quantify below. Finally, we will also use the bound on the number of variables that are in any given factor, defining

[TABLE]

to be the maximum degree of all the factors.

Lemma 14**.**

Let $p$ be a Gibbs random field with a factor graph ${\sf{C}}$ and factorization $\tfrac{1}{Z_{p}}\prod_{C\in{\sf{C}}}\phi_{p,C}$, let $N$ be the number of factors in ${\sf{C}}$, and let $\nu$ be the maximum degree of all the factors. Then any canonical estimator $q$ satisfies

[TABLE]

where the maximum is over all variables $\bm{x}_{C}$ in either a factor or a subfactor in ${\sf{C}}$.

Proof.

Applying lemma 13 to the canonical factor decomposition gives

[TABLE]

Using the canonical factor potentials from eq. 46, we have

[TABLE]

Taking the absolute value and using the triangle inequality over the $2^{|S|}$ subsets of $S$, we obtain

[TABLE]

where in the last line we used the bound on the number of factors $N$, the bound on the degree of the factors $|C_{k}|\leq\nu$, and the identity $\sum_{A\subseteq B\subseteq C}1=3^{|C|}$. The maximization in the last line is over all variables $\bm{x}_{C}$ where $C$ is a factor or a subfactor, and the result is immediate. ∎

The next lemmas let us translate the previous bound in terms of the logarithms of the conditional marginals and their empirical estimates into a bound in terms of the quantities that we naturally have control over, namely the (non-conditional) marginals and their empirical estimates.

Lemma 15**.**

For $ab>0$, \bigl{|}\log\tfrac{a}{b}\bigr{|}\leq\tfrac{|a-b|}{\sqrt{ab}}.

Proof.

This result is a consequence of a simpler inequality in terms of a single variable. For $x\geq 1$ we can use the integral representation of $\log x$ and the Cauchy-Schwarz inequality to obtain

[TABLE]

Multiplying by $-1$, we find

[TABLE]

and by letting $0<z=1/x\leq 1$, we see that the reverse inequality is true in the interval $x\in(0,1]$. This establishes that for all $x>0$,

[TABLE]

The lemma follows by letting $x=a/b$ for $ab>0$ and using some basic algebra. ∎

Lemma 16**.**

For any two strictly positive probability distributions $p(x,y)$, $q(x,y)$ on alphabets $X,Y$ and for any fixed elements $x\in X$ and $y\in Y$ we have

[TABLE]

Proof.

By lemma 15, we have

[TABLE]

From the definition of conditional probability, the triangle inequality, and the fact that $p$ and $q$ are probability distributions, we then have

[TABLE]

Combining these two inequalities yields the desired bound. We remark that this proof works equally well with $p\leftrightarrow q$, so the $p(y)$ in the denominator could also be replaced with $\max\{p(y),q(y)\}$ to get a stronger bound. ∎

The following definitions help facilitate a direct application of proposition 8 to the problem of learning the parameters of a Pauli channel with a bounded-degree factor graph. Here, as before, the minimizations in these definitions are over every factor or subfactor in the factor graph of $p$ and $q$. We first introduce a quantitative notion of positivity given by the minimum geometric mean marginal probability over a factor $C$, defined by

[TABLE]

for the fixed but arbitrary outcome $\bm{0}_{\partial C}$. We will also quantify the deviation of the local marginals from their empirical estimates by introducing the following error parameters. Let

[TABLE]

For Pauli channels with maximum degree $\nu$, the triangle inequality shows that $\epsilon_{2}\leq 4^{\nu}\epsilon_{1}$. Finally, we introduce

[TABLE]

From these definitions, we have the following lemma.

Lemma 17**.**

Let $p$ be a Gibbs random field with a factor graph ${\sf{C}}$ and factorization $\tfrac{1}{Z_{p}}\prod_{C\in{\sf{C}}}\phi_{p,C}$, let $N$ be the number of factors in ${\sf{C}}$, and let $\nu$ be the maximum degree of all the factors. Then any canonical estimator $q$ satisfies

[TABLE]

where $G$ is as in eq. 52, $\epsilon_{1}$, $\epsilon_{2}$ are as in eq. 53, and $\gamma$ is as in eq. 54.

Proof.

The result follows immediately by combining lemma 14 with lemma 16 and the definitions in eqs. 52, 53 and 54. ∎

Lemma 17 can be used together with our procedure for learning local marginal distributions (that is, complete Pauli noise models on a subset of qubits) to obtain a global guarantee on the 1-norm error between our estimate and the true distribution. As we see in the lemma, the quality of the estimate will depend on two quantities related to the local marginals as well as the precision of our local estimate, $\epsilon$ from proposition 9. The first additional quantity is the local marginal that is furthest from the noiseless case, as quantified by the 2-norm,

[TABLE]

Note that by eq. 22 we have $\|\bm{1}_{I}-p(\bm{x}_{\bar{C}})\|_{\infty}\leq\|\bm{1}_{I}-p(\bm{x}_{\bar{C}})\|_{1}=2\tfrac{d+1}{d}r\leq 3r$ where $r$ is the average error rate of the local noise channel supported on $\bar{C}$, so we expect that this term is small. The second is the minimum geometric mean marginal probability from eq. 52. Finally, the geometry of the factor graph is also important. We will need to use the bounded degree assumption to get a bound on $|\bar{C}|$ that is independent of $n$ and $N$. We will quantify the geometry of the factor graph by defining, in analogy with the definition of $\nu$ from eq. 50,

[TABLE]

In general $\bar{\nu}$ can depend on both $n$ and $N$, but it is $O(1)$ for factor graphs having factors of bounded size and where each variable participates in a bounded number of factors. From these quantities, we have the following global guarantee.

Proposition 18**.**

Let $p$ be a Gibbs distribution with a known factor graph ${\sf{C}}$ having $N$ factors, $n$ variables, and $\bar{\nu}=O(1)$ such that $p$ corresponds to the Pauli error rates for a quantum channel. Suppose that this noisy channel is GTM, $\tfrac{1}{2}$-weak, $\tfrac{1}{2}$-stable, and that $\Delta_{\bar{C}}=O(1)$. Then there exists a canonical estimator $q$ such that for all sufficiently small $\epsilon>0$, using $T$ samples with

[TABLE]

we have, with probability at least $1-\delta$,

[TABLE]

Moreover, an estimate proportional to $q$ can be found in time $\mathrm{poly}(n)$.

Proof.

We construct the following canonical estimator $q$. Using RunCB and Ratio, we obtain estimates of the marginal distribution $\hat{p}_{\bar{C}}(\bm{x}_{\bar{C}})$ for every factor $C$ of the graph. By proposition 9, each of these estimates can be obtained such that

[TABLE]

with probability $\delta$ using T_{\bar{C}}=O\bigl{(}\tfrac{1}{\epsilon^{2}}\log\tfrac{1}{\delta}\bigr{)} samples per marginal. Here the big-$O$ notation hides factors of \kappa=O\bigl{(}\log\tfrac{1}{\Delta_{\bar{C}}}\bigr{)} where $\Delta_{\bar{C}}$ is the spectral gap of the marginal channel as defined in eq. 25, but these are $O(1)$ by assumption. By the triangle inequality, these bounds apply to any marginals that are computed from these estimates as well, so they can also be applied to the subfactors in the canonical factor decomposition.

The above bounds for the empirical marginals imply that $\epsilon_{1}$ and $\epsilon_{2}$ from eq. 53 are also both bounded as $\epsilon_{1},\epsilon_{2}\leq O(1)\epsilon r^{\star}$ with probability at least $1-N\delta$. Because the noise is assumed to be 1/2-weak, we also have that $\gamma$ from eq. 54 is bounded as $\gamma\geq 1/2$. Now applying lemma 17, we have a canonical estimator $q$ such that with probability at least $1-\delta$

[TABLE]

using T=O\bigl{(}\tfrac{N}{\epsilon^{2}}\log\tfrac{N}{\delta}\bigr{)} samples.

This sample complexity can be improved by a factor of $N$ as follows. Consider the graph whose vertices are the closures of each factor, $\bar{C}$, and whose edges are present if and only if $\bar{C}_{1}\cap\bar{C}_{2}\not=\emptyset$. Because the factor graph is assumed to have constant degree, this adjacency graph also has constant degree, and so the $\bar{C}$ vertices can be partitioned into a constant number of independent sets. By construction, nonadjacent factors are conditionally independent, and so they can be estimated from the same sample. That is to say, each independent set can be chosen to comprise $\Omega(N)$ vertices, and the constituent marginals associated to these factors can be estimated in parallel. Such a choice can be found in polynomial time in $N$ using, e.g., a greedy algorithm. The estimation is repeated across the $O(1)$ distinct independent sets, which reduces the sample complexity by a factor of $O(N)$, as claimed.

The final result is obtained by rescaling $\epsilon$ to $\epsilon/N$. As all the factors are assumed to have size $O(1)$, the canonical estimator $q$ can be constructed in polynomial time in $n$ (or equivalently $N$) except for the normalizing factor. ∎

Finally, we demonstrate that significantly better bounds that scale as $\sqrt{n}$ can be obtained in terms of the differences between the marginal distributions of $p$ and $q$. The following bound is difficult to apply theoretically because we cannot directly compute or bound the marginals of the empirical distribution without first computing or bounding the partition function for the reconstructed distribution. However, the marginals of the true probability distribution can be accurately estimated by proposition 8 and the marginals of the empirical distribution can be estimated heuristically by Monte Carlo sampling from the factor potentials.

Lemma 19**.**

Given two Gibbs random fields $p$ and $q$ on $n$ variables with the same factor graph, we have the bound

[TABLE]

where $G^{\star}_{u}=\min_{\bm{x}_{\bar{u}}}\sqrt{p(\bm{x}_{\bar{u}})q(\bm{x}_{\bar{u}})}$ and $\bar{u}$ denotes the closure $u\cup\partial u$ of $u$.

Proof.

We will bound the relative entropy and then use Pinsker’s inequality. Recall the chain rule for relative entropy [82],

[TABLE]

By iteratively applying the chain rule to each variable, we have

[TABLE]

where in the second step we use the local Markov property to restrict the conditioning to the remaining part of the Markov blanket of $u$, which we denote $\bm{x}_{\partial^{+}u}$. Because the relative entropy is positive, the chain rule also proves that marginalizing or conditioning on random variables only decreases the relative entropy. Therefore, we can add back the variables from the Markov blanket of $u$ to get the upper bound

[TABLE]

Now we note that from lemma 15 we have

[TABLE]

The result then follows by applying this bound term by term, taking a maximum, and then substituting into eq. 44. ∎

VIII Conclusion

We have shown that Pauli channels can be learned with high precision using far fewer resources than previous methods. As noise can be engineered to be effectively Pauli noise using randomized compiling [2], the methods described here should be broadly applicable to learning the residual noise under randomized compiling for large-scale devices, and, as such, represent a paradigm shift for characterizing quantum hardware.

We envision hardware developers using the description of the noise to efficiently discover previously unknown properties of devices beyond regimes that have been explored at present. This should enable engineering effects to be concentrated on removing the most relevant residual errors. Moreover, we envision theoretical efforts that will result in error mitigation and correction techniques that are tailored to the specific noise afflicting a device as reconstructed using the methods in this paper. Such techniques will result in substantially better device performance by enabling optimized design of codes and decoders, bespoke fault tolerance schemes, and error-aware compiling for quantum simulation.

While we have proven rigorous upper bounds, many interesting open questions remain. For example, we have not attempted to prove lower bounds on the problems that we study in this work. It would be interesting to prove such bounds or to find procedures with even better asymptotic scaling.

Other natural avenues for future work are to explicitly prove robustness to gate-dependent noise and to generalize the method to estimate the noise on an interleaved gate [83]. With such a generalization, a natural question would be to compare the predicted circuit performance under the reconstructed noise model with assessments of circuit-level performance obtained through other protocols, such as the accreditation protocol of Ref. [84].

We have also proven our convergence guarantees in the idealized “single-shot” regime of Ref. [85]. It should be straightforward to generalize our arguments to handle the reuse of individual sequences instead of using a fresh random sequence each time. We have also made exclusive use of the qubit Pauli group, but it would be interesting to extend this work to $d$-level quantum systems as well and more formally incorporate techniques from compressed sensing.

For the reconstruction of Pauli channels as Markov random fields, it should also be possible to incorporate structure learning techniques [57, 58, 59, 60, 61] to learn the dependency structure. If the correlations are promised to be bounded degree, then we expect that structure learning can be also be done efficiently, in time $\mathrm{poly}(n)$ in the number of qubits. The sample complexity can be improved by making use of parallel samples on the same system so long as the measurements being used in parallel lie in disjoint neighborhoods of the factor graph. This can be used to improve the efficiency still further. Finally, a natural generalization of this idea is to learn a matrix product operator quantum channel [86, 87]. The state and unitary analogs of this idea have been explored before [50, 51], but not yet in a way that is robust to SPAM errors. Perhaps combining the ideas in the present paper with those of Kimmel et al. [54] could lead to a more general procedure.

Finally, there is the practical application of our methods to real world quantum devices. A first exploration along these lines can be found in a companion paper by R. Harper and the present authors [40].

Acknowledgements.

We thank Robin Harper and Wenjun Yu for comments on an earlier draft. This work was supported in part by the US Army Research Office grant numbers W911NF-14-1-0098 and W911NF-14-1-0103, the Australian Research Council Centre of Excellence for Engineered Quantum Systems (EQUS) CE170100009, the Government of Ontario, and the Government of Canada through the Canada First Research Excellence Fund (CFREF) and Transformative Quantum Technologies (TQT), the Natural Sciences and Engineering Research Council (NSERC), Industry Canada.