Multi-qubit Randomized Benchmarking Using Few Samples

TL;DR

This paper improves the efficiency of multi-qubit randomized benchmarking by reducing the number of samples needed, making it more practical for large quantum systems and providing more reliable error estimates.

Contribution

It introduces a small adaptation to RB that significantly decreases the sample complexity, especially for systems with many qubits and non-unitary noise.

Findings

Sample complexity is essentially independent of the number of qubits.

Sample size scales favorably with the average error rate.

Long sequence lengths require fewer samples for non-unitary noise processes.

Abstract

Randomized benchmarking (RB) is an efficient and robust method to characterize gate errors in quantum circuits. Averaging over random sequences of gates leads to estimates of gate errors in terms of the average fidelity. These estimates are isolated from the state preparation and measurement errors that plague other methods like channel tomography and direct fidelity estimation. A decisive factor in the feasibility of randomized benchmarking is the number of sampled sequences required to obtain rigorous confidence intervals. Previous bounds were either prohibitively loose or required the number of sampled sequences to scale exponentially with the number of qubits in order to obtain a fixed confidence interval at a fixed error rate. Here we show that, with a small adaptation to the randomized benchmarking procedure, the number of sampled sequences required for a fixed confidence interval is dramatically smaller than could previously be justified. In particular, we show that the number of sampled sequences required is essentially independent of the number of qubits and scales favorably with the average error rate of the system under investigation. We also show that the number of samples required for long sequence lengths can be made substantially smaller than previous rigorous results (even for single qubits) as long as the noise process under investigation is not unitary. Our results bring rigorous randomized benchmarking on systems with many qubits into the realm of experimental feasibility.

Equations

$$F_{\mathrm{avg}}(\mathcal{U},\tilde{\mathcal{U}}):=\int \mathrm{d}\varphi \mathrm{Tr}(\mathcal{U}(|\varphi ⟩⟨\varphi |)\tilde{\mathcal{U}}(|\varphi ⟩⟨\varphi |)),$$

$$F_{\mathrm{avg}}(\mathcal{U},\tilde{\mathcal{U}})=F_{\mathrm{avg}}(\mathcal{E},\mathcal{I})$$

$$r(\mathcal{E}):=1-F_{\mathrm{avg}}(\mathcal{E},\mathcal{I})$$

$$f(\mathcal{E}):=\frac{d{F}_{\mathrm{avg}}(\mathcal{E},\mathcal{I})-1}{d-1}$$

$$u(\mathcal{E}):=\frac{d}{d\!-\!1}\!\int\!\mathrm{d}\phi\operatorname{Tr}\!\left(\bigl{|}\mathcal{E}\bigl{(}|\phi\rangle\langle\phi|-\mathbb{1}/d\bigr{)}\bigr{|}^{2}\right).$$

$${\mathbb{E}}_{\vec{G}}(K_m)=|\mathsf{C}{|}^{-m}\sum _{\vec{G}\in {\mathsf{C}}^m}{k}_m(\vec{G})$$

$${\mathbb{E}}_{\vec{G}}(K_m){=}_{\mathrm{fit}}A{f}^m,$$

$$\text{Prob}\left[{\mathbb{E}}_{\vec{G}}(K_m)\in [k_{m,N}-ϵ,k_{m,N}+ϵ]\right]\ge 1-\delta .$$

$$N(\delta ,ϵ,m,r,\chi ,d)=-\log (2/\delta ){\left[\log \left(\frac{1}{1-ϵ}\right)\frac{1-ϵ}{{\mathbb{V}}^2+1}+\log \left(\frac{{\mathbb{V}}^2}{{\mathbb{V}}^2+ϵ}\right)\frac{{\mathbb{V}}^2+ϵ}{{\mathbb{V}}^2+1}\right]}^{-1},$$

$${\mathbb{V}}_m^2\le \frac{d^2-2}{4(d-1{)}^2}{r}^{2}m{f}^{m-1}+\frac{d^2}{(d-1{)}^2}{r}^2{u}^{m-2}\frac{(m-1){\left(\frac{f^2}{u}\right)}^m-m{\left(\frac{f^2}{u}\right)}^{m-1}+1}{(1-f^2/u{)}^2}.$$

$${\mathbb{V}}_m^2\le f^{m-1}\frac{(d^2-2)m}{4(d-1{)}^2}{r}^2+u^{m-2}\frac{d^{2}m(m-1)}{2(d-1{)}^2}{r}^2.$$

$${\mathbb{V}}_{\mathrm{SPAM}}^2\le \frac{d^2-2}{4(d-1{)}^2}{r}^{2}m{f}^{m-1}+\frac{d^2(1+4\eta )r^2}{(d-1{)}^2}\frac{(m-1){\left(\frac{f^2}{u}\right)}^m-m{\left(\frac{f^2}{u}\right)}^{m-1}+1}{(1-f^2/u{)}^2}{u}^{m-2}+\frac{2\eta dmr}{d-1}{f}^{m-1}.$$

$$\frac{(m-1)f^{2m}-m{f}^{2(m-1)}+1}{(1-f^2{)}^2}.$$

$$\frac{1}{(1-f^2{)}^2}\approx O(1/r^2)$$

$$\text{Pr}\left[|k_{m,N}-{\mathbb{E}}_{\vec{G}}(K_m)|\le ϵ\right]\ge 1-\delta .$$

$$\text{Pr}\left[|k_{m,N}-{\mathbb{E}}_{\vec{G}}(K_m)|\ge ϵ\right]$$

$$H({\mathbb{V}}^2,ϵ)={\left(\frac{1}{1-ϵ}\right)}^{\frac{1-ϵ}{{\mathbb{V}}^2+1}}{\left(\frac{{\mathbb{V}}^2}{{\mathbb{V}}^2+ϵ}\right)}^{\frac{{\mathbb{V}}^2+ϵ}{{\mathbb{V}}^2+1}}.$$

$$N=-\frac{\log (2/\delta )}{\log (H({\mathbb{V}}^2,ϵ))}.$$

$${\mathbb{V}}_{\nu }^2={\mathbb{V}}_{(\rho -\hat{\rho })/2}^2=\frac{1}{2^2}({\mathbb{V}}_{\rho }^2+{\mathbb{V}}_{\hat{\rho }}^2)$$

$${\mathbb{V}}_{\nu }^2\le \frac{1}{2}{\mathbb{V}}_{\rho }^2$$

$$\underset{A,f}{\min }\sum _{m\in \mathsf{M}}{w}_m(k_{m,N}-A{f}^m{)}^2$$

$$\mathrm{Tr}(Q\tilde{\mathcal{G}}(\rho ))=\mathrm{Tr}({\mathcal{S}}^{-1}(Q)\mathcal{S}\tilde{\mathcal{G}}{\mathcal{S}}^{-1}(\mathcal{S}(\rho ))).$$

$$\nu =\frac{\mathbf{P}}{2d}$$

$$Q=\frac{1}{2}(1+\mathbf{P})$$

$$\mathrm{Tr}(Q\mathcal{E}(\rho ))=⟨⟨Q|\mathcal{E}(\rho )⟩⟩=⟨⟨Q|\mathcal{E}|\rho ⟩⟩$$

$${\mathcal{D}}_f(\rho )=f\rho +\frac{1-f}{d}{1}_d,$$

$$\sum _{G\in {\mathsf{C}}_q}{\mathcal{G}}^{†}\mathcal{EG}={\mathcal{D}}_f$$

$$F_{\mathrm{avg}}(\mathcal{E},\mathcal{I})=\frac{(d-1)f+1}{d}.$$

$$k_m(\vec{G})$$

$${\mathbb{V}}_m^2$$

$$\displaystyle\mathbb{V}_{m}^{2}=\langle\!\langle Q^{\otimes 2}\big{|}T_{\mathsf{C}}(\mathcal{E}^{\otimes 2})^{m}\!-\!\left[T_{\mathsf{C}}(\mathcal{E})^{m}\right]^{\otimes 2}\!\big{|}\nu^{\otimes 2}\rangle\!\rangle$$

Full text

Here we show that, with a small adaptation to the randomized benchmarking procedure, the number of sampled sequences required for a fixed confidence interval is dramatically smaller than could previously be justified. In particular, we show that the number of sampled sequences required is essentially independent of the number of qubits and scales favorably with the average error rate of the system under investigation. We also investigate the fitting procedure inherent to randomized benchmarking in the light of our results and find that standard methods such as ordinary least squares optimization can give misleading results. We therefore recommend moving to more sophisticated fitting methods such as iteratively reweighted least squares optimization. Our results bring rigorous randomized benchmarking on systems with many qubits into the realm of experimental feasibility.

I Introduction

One of the central problems in the creation of large-scale, functioning quantum computers is the need to accurately and efficiently diagnose the strength and character of the various types of noise affecting quantum operations that arise in experimental implementations. This noise can be due to many factors, such as imperfect manufacturing, suboptimal calibration, or uncontrolled coupling to the external world. Tools that diagnose and quantify these noise sources provide vital feedback on device and control design leading to better quantum devices. They are also used as certification tools, quantifying a device’s ability to e.g. perform successful error correction or implement quantum algorithms. A variety of tools have been developed for this purpose, including state and channel tomography Chuang and Nielsen (1997); Poyatos et al. (1997), direct fidelity estimation (DFE) Flammia and Liu (2011); da Silva et al. (2011), gate set tomography Merkel et al. (2013); Blume-Kohout et al. (2013), and randomized benchmarking (RB) Emerson et al. (2005a); Knill et al. (2008); Magesan et al. (2011a) together with its tomographic extension randomized benchmarking tomography Kimmel et al. (2014). All of these tools have different strengths and weaknesses. State and channel tomography allow the user to get a full characterization of the quantum state or channel of interest but are subject to state preparation and measurement errors (SPAM), which place a noise floor on the accuracy of these characterizations. Moreover these protocols require resources that scale exponentially with the number of qubits even for the more efficient variants using compressed sensing Gross et al. (2010); Flammia et al. (2012), making them prohibitively expensive for use in multi-qubit systems. Randomized benchmarking tomography and gate set tomography remedy the SPAM issue, but require even more resources.

This exponential scaling with the number of qubits is problematic because even though on most quantum computing platforms multi-qubit gates are generally performed as circuits composed of one and two-qubit gates it is still vitally important to obtain aggregate measures of the behavior of multi-qubit quantum circuits. One can in principle gauge the behavior of these circuit by characterizing their component gates but such a characterization will typically give only loose bounds Dugas et al. (2016) on the behavior of the full circuit (even disregarding the possibility of correlated errors inside the circuit McKay et al. (2017)). Therefore their is a need for diagnostic tools that scale efficiently in the number of qubits. Protocols designed with such efficiency in mind, like DFE and RB, do not aspire to a full characterization of the system, but instead aim to estimate a single figure of merit that ideally captures relevant properties of the system under investigation. The figure of merit estimated by both DFE and RB is the average gate fidelity to some target state or gate. However, RB is also robust to SPAM errors (as opposed to DFE). This makes RB the protocol of choice for characterizing many candidate quantum computing platforms Knill et al. (2008); DiCarlo et al. (2009); Gaebler et al. (2012); Barends et al. (2014a); Asaad et al. (2016). Variants of RB that estimate output purity Wallman et al. (2015a), and leakage Wallman et al. (2016, 2015b); Wood and Gambetta (2018) have also been devised.

An important practical problem when using RB is choosing a number of random gate sequences that is sufficiently small to be practical experimentally, and yet gives a good estimate of the gate fidelity. This problem becomes increasingly relevant as error rates improve since estimating small errors accurately ordinarily requires more samples. Early treatments of this problem demanded numbers of sequences that were orders of magnitude larger than were feasible in experiment Magesan et al. (2012a). A more specialized analysis allowed rigorous confidence intervals to be derived for a number of random sequences comparable to the number used in experiments Wallman and Flammia (2014). However, this analysis only provided reasonable bounds on the number of sequences for short sequence lengths and for single qubit experiments while more general multi-qubit bounds had an unfavorable exponential scaling with the number of qubits being benchmarked. The restriction to short sequence lengths is also problematic because long sequences generally lead to better experimental fits Epstein et al. (2014); Granade et al. (2014).

In this paper we propose an adapted version of the standard RB protocol on the set of Clifford gates that requires little experimental overhead. For this protocol we provide a bound on the number of random sequences required to obtain rigorous confidence intervals that is several orders of magnitude sharper than previous multi-qubit bounds. Our result makes rigorous and efficient randomized benchmarking of multi-qubit systems possible using a reasonable amount of experimental resources. In particular, our bounds are approximately independent of the number of qubits being benchmarked As a special case, we also obtain bounds for the single-qubit version of RB that are valid for all sequence lengths and improve on the bounds of Ref. Wallman and Flammia (2014) for long sequence lengths. The key to the analysis of the statistical performance is a novel understanding of the representations of the Clifford group, developed in a companion paper Helsen et al. (2016). Similar representation-theoretic questions have also been studied independently by Zhu et al. Zhu et al. (2016). We also prove a precise sense in which the derived bounds are optimal. Finally we analyze the fitting procedure inherent to randomized benchmarking in light of our results. We conclude that randomized benchmarking yields data that violates the core assumptions of the Ordinary Least Squares fitting procedure, a standard tool for processing randomized benchmarking data Epstein et al. (2014). This means using OLS to analyze RB data can lead to misleading results. As an alternative we propose using the more sophisticated method of iteratively reweighted least squares optimization, which can be guaranteed to lead to correct resultsFedorov (2013); Seber and Wild (1989).

In section II we present an overview of the new contributions of this paper (equations of note here are eqs. 9 and 11) and explain their context. In section III, we discuss the implications of the new bound for experiments, and investigate it in various limits. Finally, in section IV we discuss the derivation of the new bounds and how to apply them in practice, notably with regard to the RB fitting procedure. We also prove that our results are optimal in some well specified sense. We focus on intuition and displace most of the technical proofs to the Supplementary Material. We make heavy use of techniques from group and representation theory, which are of independent interest, but were derived in a more general setting than needed for the purpose of this paper. Readers interested in the details of this part of the derivation are invited to the companion paper Helsen et al. (2016) or the closely related work of Zhu et al. Zhu et al. (2016).

Figure of merit

We begin by introducing the essential quantities we will use to state and derive our results. The central problem that RB addresses is how to efficiently obtain a rigorous figure of merit quantifying how close a physically-performed operation $\mathcal{\tilde{U}}$ (represented by a completely positive, trace preserving (CPTP) map Chuang and Nielsen (1997)) is to an ideal target operation $\mathcal{U}$, which is generally taken to be unitary, that is $\mathcal{U}(\rho)=U\rho U^{\dagger}$ for some unitary $U$ and for all density matrices $\rho$. The quality of a noisy implementation $\mathcal{\tilde{U}}$ relative to its ideal implementation $\mathcal{U}$ is quantified by the average (gate) fidelity,

[TABLE]

where $\mathrm{d}\phi$ is the uniform Haar measure over pure quantum states.

It is convenient (and always possible) to write the physically-performed operation $\mathcal{\tilde{U}}$ as the ideal operation $\mathcal{U}$ up to composition with a “noise operation”, that is we write $\mathcal{\tilde{U}}=\mathcal{E}\circ\mathcal{U}$ where $\mathcal{E}$ is a CPTP map. Note that in general the map $\mathcal{E}$ can depend on the unitary $\mathcal{U}$ being implemented. However, in this paper we shall always consider $\mathcal{E}$ to be the same for all possible unitary operations $\mathcal{U}$. This is called a gate-independent noise model. We will also work with the more general noise model $\mathcal{\tilde{U}}=\mathcal{L}\circ\mathcal{U}\circ\mathcal{R}$ where $\mathcal{R},\mathcal{L}$ are CPTP maps. This ensures compatibility of our results with recent results on RB with gate-dependent noise Wallman (2018); Proctor et al. (2017). However we can always recover the presentation given here by choosing the right gauge. This is explained in section IV.4. Because the map $\mathcal{U}$ is unitary we can also write

[TABLE]

where $\mathcal{I}$ is the identity operation. A useful quantity is the average infidelity $r$ defined as

[TABLE]

We also use the quantity $f=f(\mathcal{E})$ defined as

[TABLE]

where $d$ is the dimension of the state space. One can think of $f$ as the “depolarizing parameter” associated to the quantum channel $\mathcal{E}$. It is this quantity which randomized benchmarking can estimate. In the text, we will often drop the channel $\mathcal{E}$ from the (in)fidelity and depolarizing parameter and simply write $r(\mathcal{E})=r$ because the only channel considered in the text is $\mathcal{E}$ (or equivalently $\mathcal{R}\mathcal{L}$, see section IV.4).

We will also use another quantity associated to quantum channels called the unitarity

[TABLE]

The unitarity has the property that $u(\mathcal{E})=1$ if and only if the quantum channel $\mathcal{E}$ is unitary Wallman et al. (2015a). We will again drop the argument and write $u(\mathcal{E})=u$. Introducing this extra parameter allows us to differentiate between situations where the noise is coherent or incoherent. Randomized benchmarking behaves fundamentally different in each of these situations, as we explain in section IV.7.

The randomized benchmarking protocol

In fig. 1 we lay out our version of the randomized benchmarking protocol as it was analyzed in Knill et al. (2008); Magesan et al. (2012a); Wallman and Flammia (2014). We will perform randomized benchmarking over the Clifford group on $q$ qubits $\mathsf{C}$. This is the group of unitary operations that can be constructed by considering all possible products of CNOT gates, Hadamard gates and $\pi/4$ phase gates on the $q$ qubits Gottesman (1998). We make two essential changes to the standard randomized benchmarking protocol, both of which lead to better guarantees on the precision of randomized benchmarking.

• •

A first modification is to perform each randomized benchmarking sequence twice, but with different input states $\rho,\hat{\rho}$ and then subtracting the result. This is equivalent to performing standard randomized benchmarking with the “input operator” $\nu=\frac{1}{2}(\rho-\hat{\rho})$. A similar idea was suggested in Knill et al. (2008); Granade et al. (2014); Muhonen et al. (2015); Fogarty et al. (2015). The factor $(1/2)$ is not strictly necessary but it allows for a fairer comparison between the original benchmarking protocol and our proposal 111In particular this factor of two ensures that “signal ranges” of the two protocols are equal, that is, standard RB starts at $1$ and decays to $1/2$ for large $m$ and the new protocol starts at $1/2$ and decays to [math]..

• •

Secondly, we do not assume the ideal measurement operator to be the projector on the $|0\cdots 0\rangle$ state. Instead we perform some stabilizer measurement related to a pre-chosen Pauli matrix $\mathbf{P}$. An experimentally good choice would be for instance $\mathbf{P}=Z^{\otimes q}$ but our results hold for any choice of Pauli operator. Correspondingly we pick the input states to be some (impure) states $\rho,\hat{\rho}$ with support on the positive, resp. negative, eigenspaces of the Pauli operator $\mathbf{P}$ That is, we would like to prepare the impure states $\rho=\frac{I+\mathbf{P}}{2d},\hat{\rho}=\frac{I-\mathbf{P}}{2d}$.

Both of these adjustments are done with the purpose of lowering the experimental requirements for rigorous randomized benchmarking. Our first change to the RB protocol, performing randomized benchmarking with a state difference, has two beneficial effects: (1) It changes the regression problem inherent to randomized benchmarking from an exponential fit with a non-zero off-set to an exponential fit (see eq. 7). This eliminates a fitting parameter, lowering experimental requirements. (2) It lowers the statistical fluctuations of randomized benchmarking regardless of what input states are actually used. This improvement is mostly noticeable in the limit of long sequence lengths. We discuss this in more detail in section IV.9.

A much stronger improvement to the statistical fluctuations inherent to randomized benchmarking stems from our second change to the RB protocol; preparing states and performing measurements proportional to $\mathbb{1}+\mathbf{P}$ where $\mathbf{P}$ is a Pauli operator. This change allows us to prove a radically sharper bound on the statistical fluctuations induced by finite sampling relative to preparing other input states. This is discussed in section IV.5 (see in particular eq. 36). In section IV.7 we argue that this behavior is not an artifact of our proof techniques but rather inherent to the statistical behavior of randomized benchmarking. Note that for a single qubit the state $(I\pm\mathbf{P})/2$ is in fact a pure state for any choice of $\mathbf{P}$ (in particular $(I+Z)/2=|0\rangle\langle 0|$). Note that (1) and (2) both reduce the amount of resources needed in a different and independent manner. Using a difference of two input states amounts to effectively “preparing” a traceless input operator. The tracelessness of this operator has two distinct effects. The first effect is that it fixes the constant offset of the decay to be zero, thereby eliminating a fitting parameter. The second effect, which is more subtle, is that it eliminates in the variance expression a representation (which has support on the identity matrix), and hence an extra term in the sequence variance. This means the sequence variance is reduced compared to the sequence variance of standard RB. This effect remains even in the case of imperfect state preparation, as the difference of two density matrices is always traceless (assuming no leakage during the preparation)

As seen in fig. 1 the RB protocol starts by, for a given sequence of Clifford operations $\vec{G}$ of length $m$, computing the expectation value $p_{m}(\vec{G})(\rho)$ of an observable $Q$ for two different input states $\rho$ and $\hat{\rho}$. We subtract these two numbers to obtain a number $k_{m}(\vec{G}):=\frac{1}{2}(p_{m}(\vec{G})(\rho)-p_{m}(\vec{G})(\hat{\rho}))$. Next we obtain an average of this quantity over all possible sequences $\vec{G}$.

[TABLE]

This average over all possible Clifford strings of length $m$ can be fitted for various values $m$ to the exponential decay curve

[TABLE]

with two fitting parameter $A$ and $f$. In the case where all gates performed in the experiment suffer from the same noise, that is $\mathcal{\hat{G}}=\mathcal{E}\circ\mathcal{G}$ for all Clifford operations $\mathcal{G}$ the number $f$ can be interpreted as the depolarizing parameter of the channel $\mathcal{E}$ (as defined in eq. 4) giving an estimate of the average fidelity of the noisy operation $\mathcal{\hat{G}}$ w.r.t. its ideal version $\mathcal{G}$.

In practice the number of possible sequences for a given $m$ is too large to average over completely. Instead one averages over a randomly sampled subset of sequences, which generates an empirical estimate $k_{m,N}$ the validity of which we can interpret using confidence regions. A confidence region, for some set confidence level $1-\delta$ and size $\epsilon$, is an interval $[k_{m,N}-\epsilon,k_{m,N}+\epsilon]$ around the estimate $k_{m,N}$ such that the probability that the (unknown) parameter $\mathbb{E}_{\vec{G}}(K_{m})$ lies in this interval with probability greater than $1-\delta$, i.e.,

[TABLE]

These confidence intervals, obtained for various values of sequence length during the experiment can then be used in the fitting procedure eq. 7 to generate a confidence interval around the empirical estimate $\hat{F}$ for the true channel average fidelity $F_{avg}(\mathcal{E},\mathcal{I})$. This can be done using standard statistical procedures (see e.g. Beck and Arnold (1978)). The number of random sequences $N$ used to obtain $k_{m,N}$ will depend on $\epsilon$ and $\delta$ which are set before the beginning of the experiment, and in general also on some prior estimate of the infidelity $r$ and unitarity $u$. The rest of the paper will be concerned with making this $N$ as small as possible given $\delta$ and $\epsilon$ and (if possible) an a priori bound on the average infidelity $r$.

II Results

In this section we state the main contributions of the paper. We present practical bounds on the number of sequences required to obtain rigorous confidence intervals for randomized benchmarking using the Clifford group under the assumption that the expectation value difference $k_{m}(\vec{G})$ for a given Clifford sequence $\vec{G}$ is estimated easily to a very high precision. This means we assume that any uncertainty on the number $k_{m}$ is mostly due to the fact that we only sample $N$ sequences $\vec{G}$. Epstein et al. (2014); Wallman and Flammia (2014), or equivalently that the uncertainty on the number $k_{m}(\vec{G})$ for a fixed sequence $\vec{G}$ is negligible. In order to construct a $1-\delta$ confidence interval of size $\epsilon$ around a randomized benchmarking sequence average $k_{m,N}$ with sequence length $m$, system dimension $d$ and a prior estimate of the channel infidelity $r$ and unitarity $u$ one needs to average over $N$ random sequences where $N$ is given by Hoeffding (1963):

[TABLE]

where $\mathbb{V}^{2}$ is the variance of the distribution of the samples $k_{m}(\vec{G})$ from a uniform distribution over the Clifford sequences $\vec{G}$. This variance is given below.

The variance of randomized benchmarking

The most important contribution of this paper is a bound on the number of sequences $N$ needed for multi-qubit randomized benchmarking. Previous bounds for multi-qubit RB Wallman and Flammia (2014); Magesan et al. (2012a) are either prohibitively loose or scale exponentially with the number of qubits. Our new bounds, which are derived in detail in theorem 1 of the Supplementary Material, resolve both these issues using techniques from representation theory, enabling multi-qubit RB with practical numbers of random sequences.

Variance bound for SPAM-free multi-qubit RB

For states and measurements that are (very close to) ideal, section IV.5 yields a bound on the variance in terms of the sequence length $m$, the infidelity $r$, the unitarity $u$ and the system size $d$. It is given by

[TABLE]

This bound is asymptotically independent of system size $d$.

To illustrate the improvements due to our bound, consider a single qubit $(d=2)$ RB experiment with sequences of length $m=100$ and average infidelity $r\leq 10^{-4}$. To obtain a rigorous $99\%$ confidence interval of size $\epsilon=10^{-2}$ around $p_{m,N}$, Ref. Wallman and Flammia (2014) reported that $N=145$ random sequences were needed (In the case of perfect state preparation and measurement) while our bounds imply that $N=173$ random sequences are sufficient. However, the new bound has substantially better scaling with $m$. For instance, with $m=5000$, $\epsilon=0.05$ and other parameters as above, our bound only requires $N=470$ compared to the $N=1631$ required by the single qubit bound of Ref. Wallman and Flammia (2014). We illustrate the difference in scaling of the number of sequences needed for a given confidence interval with respect to sequence length $m$ in fig. 2.

A notable upper bound on eq. 9, which is easier to work with, is

[TABLE]

This bound can be further weakened and simplified by setting $u=1$, yielding an upper bound on the variance that is independent of the unitarity. This bound will however rapidly become trivial with increasing sequence length.

Variance bound including SPAM

The above variance bound is sensitive to SPAM errors, which introduce terms into the variance which scale linearly in the infidelity $r$. In theorem 1 of the Supplementary Material, we prove that in the presence of SPAM errors the variance is bounded by

[TABLE]

The correction factor $\eta$ only depends on SPAM. As we show in section IV.8, this SPAM dependence is impossible to avoid if one wants to retain the preferred quadratic scaling in infidelity $r$. This bound is also asymptotically independent of the number of qubits. This means we can perform rigorous randomized benchmarking even in the limit of very many qubits. We illustrate the difference in scaling with respect to system size in fig. 2.

To illustrate the improvements our methods yield we can again compare to Wallman and Flammia (2014). Consider a system with $4$ qubits, that is, $d=16$, with sequence length $m=100$, an a priori estimate of $r\leq 10^{-4}$, and $\eta=0.05$. For a $99\%$ confidence region of size $\epsilon=10^{-2}$ the previous best known bound for multiple qubits Wallman and Flammia (2014) would require $N=3\times 10^{5}$ random sequences, while our dimension independent bound from eq. 11 only requires $N=249$.

Optimality of results

We also prove (see section IV) that for arbitrary SPAM a bound on the variance which is linear in the infidelity $r$ is in fact optimal. This means the result stated above is in some sense the best possible bound on the variance of a randomized benchmarking sequence. It is important to note that this optimality result also holds when RB is performed using a different set of gates than the Clifford group and also when one considers the standard protocol Knill et al. (2008); Magesan et al. (2011a) as opposed to the protocol involving differences of quantum states which we presented in this paper.

Both the SPAM and SPAM-free variance bound also approach a constant independent of the infidelity $r$ in the limit of large sequence length $m$ when the unitarity is one, that is when the noise in the system is purely coherent. In section IV.8 we argue that this behavior is not an artifact of the proof techniques used but is in fact a generic feature of a randomized benchmarking procedure with a unitary noise process.

Fitting procedure

In section IV.3 we discuss the consequences of eqs. 9 and 11 on the fitting procedure used to fit the data $\{k_{m,N}\}$ generated by fig. 1 to the RB fitting relation eq. 7. Our results show that the variance of randomized benchmarking data is strongly heterogeneous with respect to the sequence length $m$. This invalidates the key assumption of homogeneity of variance (homoskedasticity) Seber and Wild (1989) that is necessary for the correct functioning of Ordinary Least Squares (OLS), the standard method used for fitting RB data Epstein et al. (2014). Because of this inferences drawn from can give misleading results when applied to RB data. We recommend switching from OLS to the more sophisticated method of Iteratively Reweighted Least Squares, which can deal with non-homoskedastic data.

III Discussion

In this section we will discuss the behavior of the variance bound eqs. 9 and 11 in various regimes. Of interest are its scaling with respect to the number of qubits in the system, the presence of state preparation and measurement noise and varying amounts of coherence in the noise process.

III.1 Scaling with number of qubits.

We begin by discussing the effect of the number of qubits in the system on the variance and the number of necessary sequences.

As illustrated in fig. 2 (red full) and as can be seen from eq. 9, the derived bound is almost independent of the number of qubits $q$ (where $d=2^{q}$). In fact, the bound on the variance decreases asymptotically to a constant in the limit of many qubits despite the number of possible sequences (that is, $\lvert\mathsf{C}\rvert^{m}$) increasing exponentially with the number of qubits. This constitutes a notable improvement over previous multi-qubit variance bounds with an explicit dependence on the infidelity (dashed green in fig. 2), given in Wallman and Flammia (2014) which had a linear scaling with infidelity but scaled exponentially with the number of qubits. The qualitative behavior of the variance bound in terms of dimension matches a trivial bound on the number of sequences, which can be made by noting that the numbers $k_{m,N}$ are sampled from a distribution bounded on an interval of unit size (and hence has variance at most $1/4$ (dashed blue in fig. 2)) but is much sharper in absolute terms due to its quadratic dependence on the infidelity $r$.

To further illustrate the behavior of the bound, fig. 3(a) shows the number of sequences needed for a $99\%$ confidence interval around $k_{m,N}$ of size $5r$ versus the number of qubits in the system for various values of $r$ ranging from $5\cdot 10^{-3}$ to $10^{-4}$ and sequence length $m=100$. The size of $\epsilon$ was chosen to reflect that for fixed sequence length a smaller infidelity will lead to the need for greater precision around $k_{m.N}$ for a successful fit to the exponential eq. 7 Epstein et al. (2014). This plot was made using the unitarity independent bound in eq. 10 for ideal SPAM, but similar plots can be made for non-negligible SPAM errors using eq. 11. Note also that greater numbers of sequences are needed when the infidelity is small even though the variance in eq. 9 decreases with infidelity. This is due to our setting of the size of the confidence interval and reflects the statistical truism that more samples are in general needed to detect small differences.

III.2 Effects of SPAM terms

In practice it will always be the case that the input state difference $\nu$ and the output measurement POVM element $Q$ are not ideal. This means that in general we must take into account the contributions from non-ideal SPAM when calculating the number of required sequences. These contributions scale linearly in the infidelity $r$ (see eq. 11) rather than quadratically and so will increase the amount of required sequences. The degree to which $\nu$ and $Q$ deviate from the ideal situation is captured by the prefactor $\eta$ (see section IV for more on this factor). To illustrate the effect of the SPAM terms on the variance we plot in fig. 3(b) the variance versus the infidelity $r$ using eq. 11 taking the sequence length $m=100$ and the dimension of the system $d=16$ (four qubits) for SPAM of size $\eta\in\{0,0.01,0.05,0.1,0.5\}$. From this plot we note that for non-zero $\eta$ the variance, and hence the amount of sequences needed increases rapidly, especially in the regime of small $r$. This is due to the fact that increasing the SPAM contribution interpolates the variance between a regime where the terms quadratic in infidelity $r$ are dominant and a regime where the terms linear in infidelity $r$ are dominant. This means that, especially when dealing with systems with very small $r$ it is advantageous to try to suppress SPAM errors. In section IV.7 we show that this type of quadratic-to-linear interpolation behavior is in fact optimal for the variance of randomized benchmarking.

III.3 Scaling with sequence length

Of more immediate relevance is the scaling of the bound with the sequence length. It is easy to see that the variance bound eq. 9 scales quadratically in the sequence length $m$ for any noise process when the sequence length is small (see also eq. 10) but when the sequence length is very long the precise nature of the noise under consideration heavily impacts the variance. If the noise is purely coherent, i.e. the unitarity $u=1$, we see that the scaling of the second term in eq. 9 is set by the factor

[TABLE]

In the limit of $m$ going to infinity this factor goes to

[TABLE]

which means the variance eq. 9 converges to a constant independent of the infidelity $r$. This behavior for unitary noise is strikingly different from the behavior for incoherent noise, that is $u<1$. Here we see that the variance in the limit of long sequences is dominated by the exponential terms $u^{m-2}$ and $f^{2(m-1)}$. Since $f$ and $u$ are strictly less than one by the assumption of incoherence, the variance will decay to zero in the limit of long sequences. As $u\geq f^{2}$ for all possible noise processes Wallman et al. (2015a) the decay rate will be dominated by the size of the unitarity. This is also evident in fig. 4(a). In this figure we see the number of sequences needed (as given by eq. 9) versus sequence length $m$ for fixed infidelity $r=0.1$ and dimension $d=16$, and a fixed confidence interval $\delta=0.99,\epsilon=0.01$ but for different values of the unitarity $u$. Here we have chosen $u=(\kappa+(1-\kappa)f^{2})$ for $\kappa\in\{0.2,0.4,0.6,0.8,1\}$ corresponding to the situations where the noise is relatively incoherent going all the way up to a situation where the unitarity is one. We see that for $u<1$ the number of sequences needed first rises quadratically, tops out and subsequently decays to zero whereas in the case of $u=1$ the number of sequences needed keeps rising with sequence length $m$ until it tops out at some asymptotic value. In section IV.8 we argue that this behavior is not a feature of the variance bound but rather a feature of the variance of randomized benchmarking itself. Therefore, in the case of highly unitary noise, we recommend performing more experiments at shorter sequence lengths rather than trying to map out the entire decay curve.

Another noteworthy feature of the variance bound eq. 9 is the fact that, for non-unitary noise (that is $u<1$) it is in general not monotonically increasing in infidelity $r$. Rather, for a fixed sequence length, the variance increases at first with increasing infidelity but then peaks and decays towards zero. This behavior is illustrated in fig. 4(b). Here we plot a contour plot of the variance with infidelity on the $y$-axis ($r\in[0.01,0.1]$) and sequence length $m$ on the $y$-axis ($m\in[1,100]$) and have set the unitarity to $u=(f^{2}+1)/2$ corresponding to relatively incoherent noise. The take-away from this plot is that it is not enough to have an upper bound on the infidelity to get an upper bound on the variance, rather one must have both an upper and a lower bound on the variance to make full use of the bound eq. 9. Note that the looser upper bound eq. 10 does not share this behavior and always yields an upper bound on the variance given an upper bound on the infidelity $r$.

On the other hand, when the underlying noise process is unitary, that is $u=1$ the variance does increase monotonically with increasing $r$. This strikingly different behavior is illustrated in fig. 4 (c). Here we plot a contour plot of the variance with infidelity on the $y$-axis ($r\in[0.01,0.1]$) and sequence length $m$ on the $y$-axis ($m\in[1,100]$) and have set the unitarity to $u=1$ corresponding to fully coherent noise.

III.4 Future work

An important caveat when applying the confidence bounds is the assumption of gate and time independent noise (this can be relaxed to Markovian, gate independent noise Wallman and Flammia (2014)). This is an assumption that many analyses of RB suffer from to various degrees, hence a major open problem would be to generalize the current bounds to encompass more general noise models. Note however, that since our upper bound captures the correct functional behavior of the RB variance with respect to sequence length (for gate and time independent noise) one could in principle check if these assumptions hold true by computing estimates for the variance at each sequence length (from the measured data) and checking if these estimates deviate significantly from the proposed functional form.

Recent work has also argued that the exponential behavior of randomized benchmarking is robust against Markovian gate-dependent fluctuations Wallman (2018). This however comes at a substantial increase in mathematical complexity. We suspect that similar robustness statements can be made for the variance of randomized benchmarking but new mathematical tools will be needed (perhaps using the Fourier analysis framework proposed recently in Merkel et al. (2018)) to make this suspicion rigorous.

Our work can be straightforwardly extended to interleaved RB Magesan et al. (2012b). However the dominant source of error in the interleaved RB protocol is usually systematic rather than stochastic (due to the fact that the protocol does not yield an estimate of the interleaved gate fidelity but rather provides upper and lower bounds). Interleaved RB essentially consists of two RB experiments: a reference experiment and an interleaved experiment, the latter of which has an extra ‘interleaved gate’ inserted between the random gates of the standard RB protocol. Hence the fidelity extracted from the second experiment corresponds to the fidelity of the composition of the noise due to the random gates and the noise due to the interleaved gate. An estimate of the fidelity of the interleaved gate is then extracted by considering the ratio of the fidelity of the random gates (from the reference experiment) and the fidelity of the above composition. However, the fidelity of a composition of two noise maps is in general not equal to the product of the fidelities of the individual maps and can, depending on the specifics of the noise processes, differ quite radically. Hence in the absence of more knowledge about the underlying noise processes, IRB gives an inaccurate estimate of the fidelity of the interleaved gate. This inaccuracy is not remedied by reducing the imprecision of the fidelity estimates (for a fixed amount of resources), which is what we provide here. And since the inaccuracy due to this lack of fidelity-composition can be much larger than the imprecision for even a modest amount of resources it is less useful to spend significant energy on increasing precision in IRB.

Moreover, it should be noted that while randomized benchmarking is efficient in the complexity theoretical sense, i.e. the amount of resources needed scales polynomially with the number of qubits in the system, the amount of resources required is still significant, and no RB experiment has been performed beyond $3$ qubits so far McKay et al. (2017). Recently several protocols have been devised and implemented that are similar to randomized benchmarking but less resource intensive Helsen et al. (2018); Erhard et al. (2019); Xue et al. (2019), making larger-scale characterization of multi-qubit systems possible. We suspect the bounds derived in this paper can be adapted to these new proposal but we leave this for future work.

Also, successful and rigorous randomized benchmarking not only depends on the number of random sequences needed per sequence length but also on the fitting procedure used to fit the points generated by randomized benchmarking of various lengths to a decay curve in order to extract an estimate of the average gate fidelity. Finding the optimal way to perform this fitting procedure is still an open problem Epstein et al. (2014). Accounting for heteroskedasticity, as we have done here, can be considered a first step in this direction. Performing this accounting is standard practice in statistics but does not seem to be in widespread use in the experimental community. One could also consider directly estimating the variance at each sequence length from obtained data and then using these estimates directly as inputs to a weighted least-squares fitting procedure. We however believe that the parametric model we propose here will be more efficient in terms of data needed for a fixed precision.

Finally, a major theoretical open problem is the extension of the present bounds to non-qubit systems, different varieties of randomized benchmarking Dugas et al. (2015); Cross et al. (2016); Barends et al. (2014b), and to different $2$-designs Dugas et al. (2015); Dankert et al. (2009); Turner and Markham (2016) or even orthogonal $2$-designs Hashagen et al. (2018); Harper and Flammia (2018). If these $2$-designs are assumed to be groups, similar techniques from representation theory might be used Gross et al. (2007) but how this would be done is currently unknown.

IV Methods

In this section, we will discuss the new contributions in detail, and explain how to apply them in an experimental setting. We will give a high level overview of the proof of the bound on the variance of a randomized benchmarking sequence; full details can be found in the Supplementary Material. We will also discuss the behavior of noise terms in the case of non-ideal SPAM and prove that the bounds we obtain are in some sense optimal. Finally we briefly comment on how the variance changes when performing regular randomized benchmarking (using an input state $\rho$ rather than an input state difference $\nu=\frac{1}{2}(\rho-\hat{\rho})$).

IV.1 Estimation theory

In this section, we review confidence intervals and relate the bounding of confidence intervals to the bounding of the variance of a distribution. A first thing of note is that all the variance bounds stated in section II are dependent on the infidelity $r$. The appearance of $r$ in the bound might strike one as odd since this is precisely the quantity one tries to estimate through RB. It is however a general feature of estimation theory that one needs some knowledge of the quantity one tries to estimate in order to use nontrivial estimation methods Beck and Arnold (1978). Note also that while our results are stated in frequentist language, they should also be translatable to Bayesian language, that is, as credible regions on the infidelity given prior beliefs as in Ref. Granade et al. (2014) for example. Bayesian methods are more natural because our bounds depend on prior information about the infidelity, however, a full Bayesian treatment would involve the fitting process, obscuring our primary technical result, i.e. the variance bounds.

Let us now discuss how to use the variance bounds to construct confidence intervals around numbers $k_{m,N}$. We can in general define a $1-\delta$ confidence interval of size $\epsilon$ to be

[TABLE]

Once we have an upper bound on the variance $\mathbb{V}_{m}^{2}$ of an RB distribution we can relate this to an upper bound on number of required sequences through the use of concentration inequalities.

Note that for the case of randomized benchmarking there are two sets of confidence parameters. $(\delta_{N},\epsilon_{N})$ is associated with estimating the average over all possible Clifford sequences, where the relevant parameter is the number of performed sequences $N$ and $(\delta_{L},\epsilon_{L})$ is associated with getting an estimate for the survival probability difference $k_{m}(\vec{G})$ for a given fixed sequence. Here the relevant parameter is $L$, the number of times a single sequence is performed. Since in practice $L<\infty$ there will be some finite $(\delta_{L},\epsilon_{L})$ confidence region around the survival probability difference $k_{m}(\vec{G})$ for a given sequence $\vec{G}$. So in general, when looking at a $\epsilon,\delta$ confidence region for an RB procedure of a given length on should look at $(\epsilon_{N}+\epsilon_{L},\delta_{N}+\delta_{L})$ confidence regions. In what follows we will assume that $L$ is high enough such that $\epsilon_{L},\delta_{L}$ are negligible relative to $(\delta_{N},\epsilon_{N})$. This approach is motivated by experimental realities where it is usually much easier to perform a single string of Cliffords many times quickly than it is to generate, store and implement a large number of random sequences.

For a given variance $\mathbb{V}^{2}$ we can relate the number of sequences $N$ needed to obtain $1-\delta$ confidence intervals of size $\epsilon$ using the following concentration inequality due to Hoeffding Hoeffding (1963):

[TABLE]

with

[TABLE]

We can invert this statement to express the number of necessary sequences $N$ as a function of $\delta,r,\epsilon$ as

[TABLE]

Note that this expression can also be inverted to yield a bound on $\delta,\epsilon$ in terms of a given number of samples $N$. This identity heavily depends on the size of the variance $\mathbb{V}_{m}^{2}$.

IV.2 State preparation and measurement costs

We have argued that our adapted RB protocol allows for a reduction in the number of needed sequences to make rigorous estimates. However implicit in this cost reduction argument is the assumption that estimating the number $k_{m}(\vec{G})$ for a fixed sequence $\vec{G}$ is not more costly than estimating the number $p_{m}(\vec{G})$. Here we justify this assumption for the two changes we made to the randomized benchmarking protocol: using a state difference as input and using an impure input state defined by a single Pauli matrix. In the following we forgo rigor in favor of intuition. We are however only applying standard statistical techniques that can easily be made rigorous.

**State difference

**At first glance one might think that estimating the same sequence twice for difference input states as we propose yields a two-fold overhead in the number of samples per sequence. To see that this is not the case consider the variance $\mathbb{V}_{\rho}^{2}$ associated with estimating the expectation value for a single sequence for a single state $\rho$. From the standard rules of error addition we now have, for the state difference $\nu=(\rho-\hat{\rho})/2$ that

[TABLE]

since the random variables associated to $\rho$ and $\hat{\rho}$ are independently distributed (making the covariance zero). Now assuming that $\rho$ incurs the largest variance, we get

[TABLE]

which means that estimating the expectation value of a single sequence for a difference of states is statistically not harder than estimating it for a single state.

**Optimal input state and measurement

**In our adapted RB procedure we call for preparing the input states $\rho=\frac{\mathbb{1}+\mathbf{P}}{2},\rho=\frac{\mathbb{1}-\mathbf{P}}{2}$ for some Pauli matrix $\mathbf{P}$ and measuring the output operator $\mathbf{P}$. This is different from standard RB where one is asked to prepare and project onto the all zero state $|0\ldots 0\rangle\langle 0\ldots 0|$. We argue that performing RB this way is not more costly than using the standard approach. For concreteness we shall set $\mathbf{P}=Z^{\otimes q}$. Measuring the expectation value of the operator $Z^{\otimes q}$ is trivial; one simply measures all qubits in the standard basis (as one would do in standard RB) and then computes the parity of the outcome. Since standard basis states with even parity precisely span the positive eigenspace of $Z^{\otimes q}$ this amounts to measuring the expectation value of $Z^{\otimes q}$. Preparing the states $\rho=\frac{\mathbb{1}+Z^{\otimes q}}{2},\rho=\frac{\mathbb{1}-Z^{\otimes q}}{2}$ is a little more involved. The state $\rho$ is a probabilistic mixture of all computational basis states $|x\rangle$ of even parity. By the linearity of expectation one could compute (for a fixed Clifford sequence $\vec{G}$ ) the survival probability $p_{m}(\vec{G},|x\rangle)$ and then compute $p_{m}(\vec{G},\rho)=2^{-q/2}\sum_{x}p_{m}(\vec{G},|x\rangle)$. This requires measuring $2^{2/q}$ expectation values $p_{m}(\vec{G},|x\rangle)$, making this approach not scalable. We can remedy this by realizing that we are only interested in a good estimate of the mean $p_{m}(\vec{G},\rho)$. Considering $p_{m}(\vec{G},|x\rangle)$ to be the mean of a Bernoulli random variable with outcomes [math] and $1$, and thus $p_{m}(\vec{G},\rho)$ to be the mean of a normalized binomial distribution we can estimate this mean efficiently by sampling $|x\rangle$ at random (with even parity), estimating $p_{m}(\vec{G},|x\rangle)$ and then computing the empirical mean. Moreover, since we do not need to know the means $p_{m}(\vec{G},|x\rangle)$ very well to get a good estimate of $p_{m}(\vec{G},\rho)$ the about of single data points (clicks) gathered to estimate $p_{m}(\vec{G},\rho)$ is not higher than it would be to accurately estimate $p_{m}(\vec{G},|\psi\rangle)$ for $|\psi\rangle$ some pure state.

IV.3 The fitting procedure

In the previous section we outlined how to use the bound eq. 9 to construct confidence intervals around $k_{m,N}$. However, we have not yet discussed how to integrate the variance bound eq. 9 into the fitting procedure required by eq. 7. A fitting procedure is any method that takes in the set of data points $\{k_{m,N}\}_{m}$, with $m\in\mathsf{M}$, where $\mathsf{M}$ is some set of integers and outputs a tuple $(A^{*},f^{*})$ such that $A^{*}{f^{*}}^{m}$ is a ‘good’ description of the data $\{k_{m,N}\}_{m}$. There are many ways to approach this problem, we refer to Seber and Wild (1989) for a good overview, and finding an optimal procedure is outside the scope of this paper. However we would like to discuss the most commonly used fitting procedure: Ordinary Least Squares (OLS) in the light of the bounds eqs. 9 and 11.

**Ordinary least squares

**Given data $\{k_{m,N}\}_{m}$ and the function $F(A,f)=Af^{m}$ the OLS procedure returns estimates $(\hat{A},\hat{f})$. Through a linearization procedure, as outlined for RB in Epstein et al. (2014), confidence intervals can then be constructed around these estimates. However, for this procedure to yield correct results each data point $k_{m,N}$ must be distributed around $\mathbb{E}_{\vec{G}}(K_{m})$ with the same variance (Seber and Wild, 1989, Chapter 2.8). This assumption, called homoskedasticity in the statistics literature, is not universally valid for randomized benchmarking data $\{k_{m,N}\}_{m}$. This shows in the functional form of the upper bound eq. 9, which strongly depends on the sequence length and from eq. 36 one can see that this is not an artifact of bounding techniques but rather an innate feature of RB data. Moreover OLS assumes that the variance of $k_{m,N}$ is independent of the fitting parameters $A,f$, an assumption which is also explicitly violated in RB data. The violation of these two assumptions (homoskedasticity and independence of fitting parameters) creates problems when performing OLS on the RB data $\{k_{m,N}\}_{m}$. In particular OLS no longer provides an unbiased estimate of the standard error on the fitting parameters $(f,A)$ (Seber and Wild, 1989, Chapter 3.3), which can lead to mis-estimation of confidence intervals around the fitting parameters. Therefore we recommend using a more sophisticated approach.

**Iteratively re-weighted least squares

**Heteroskedasticity (violation of homoskedasticity) and functional dependence of the data distribution on the fitting parameters are well studied problems, and many robust solutions are available. Here we will focus on one particular solution called Iteratively Re-weighted Least Squares (IRLS). For the purposes of this construction we will assume that the data $\{k_{m,N}\}_{m}$ is drawn from a random variable with mean $\mathbb{E}_{\vec{G}}(K_{m})$ and variance $\mathbb{V}_{m}^{2}(m,r)/N$. IRLS constructs estimates for the parameters $(A,f)$ by minimizing the function

[TABLE]

where the weights $w_{m}$ can depend on $f$ and $A$. Under the assumption that eq. 9 is the actual variance $\mathbb{V}_{m}^{2}$ up to a constant factor we can set the weights (Seber and Wild, 1989, Section 2.8.8) to be $w_{m}=w(f,u,m)=1/\sigma(f,u,m)$ where $\sigma$ is the RHS of eq. 9 (if one suspects that $\eta\neq 0$ the eq. 11 can be used instead). We note that this procedure is fairly robust against misspecification of the weights, and moreover that $\sigma$ captures the behavior of $\mathbb{V}_{m}^{2}$ with respect to the sequence length very well (see section IV.5). IRLS now proceeds in the following manner:

It as been shown (Fedorov, 2013, Page 45) (under some mild regularity conditions) that this algorithm converges to estimates $\hat{A},\hat{f}$. If the weights $w_{m}$ are exactly proportional to the variance $\mathbb{V}_{m}^{2}$ then these estimates are asymptotically consistent. In section B.3 in the Supplementary Material we provide a detailed estimate of how close the estimate $\hat{f}$ is to the real depolarizing parameter $f$ in terms of the number of data points in $\{k_{m,N}\}_{m}$ and the number of sequences $N$ sampled per data point.

Finally we would like to note that we have in this procedure kept the number of sequences $N$ constant for varying $N$. It is however possible to let $N$ depend on the sequence length $m$. One choice would be to vary $N$ proportionally to $\mathbb{V}_{0}^{2}$ (assuming a good estimate of $f$ is available). In this scenario, since $k_{m,N}$ is drawn from a distribution with variance $\mathbb{V}_{m}^{2}/N$ this would remedy the issue with heteroskedasticity and OLS could be used to provide reliable fitting.

IV.4 Gate dependent noise and gauge invariance

In recent work Proctor et al. (2017); Wallman (2018) it has been noted that the relation between the parameter estimated by randomized benchmarking and the average fidelity is less than straightforward when the noise channel is allowed to depend on the gate being implemented, that is $\tilde{\mathcal{G}}=\mathcal{E}_{G}\mathcal{G}$. At the heart of the issue is that the only quantities measurable in the lab, probabilities of the form $\operatorname{Tr}(Q\tilde{\mathcal{G}}(\rho))$ for a state $\rho$ and an observable $Q$ are gauge invariant. That is, for any invertible superoperator $\mathcal{S}$ we have that

[TABLE]

This difficulty can be remedied by considering a more general noise model. Instead of choosing $\tilde{\mathcal{G}}=\mathcal{E}\mathcal{G}$ one chooses $\mathcal{\tilde{G}}=\mathcal{L}_{G}\mathcal{G}\mathcal{R}_{G}$ for superoperators $\mathcal{R}_{G},\mathcal{L}_{G}$ Wallman (2018). The individual operators $\mathcal{R}_{G},\mathcal{L}_{G}$ are not gauge invariant but the combined operator $\mathcal{R}_{G}\mathcal{L}_{G}$ is. Since in this paper we deal exclusively with gate-independent noise we can choose the gauge such that $\mathcal{L}=\mathcal{I}$ and $\mathcal{R}=\mathcal{E}$ but our results also hold for the more general choice of gauge with the express caveat that our bounds then work in terms of the infidelity $r$ and unitarity $u$ of the noise in between gates $\mathcal{RL}$. That is we have $r=r(\mathcal{RL})$ and $u=u(\mathcal{RL})$. It is possible to see this explicitly by making the substitution $\mathcal{E}\rightarrow\mathcal{RL}$ in all steps of the derivation of the variance bound in section IV.5 (and theorem 1 in the Supplementary Material).

IV.5 Variance bound

In this section we present a derivation of the multi-qubit variance bound in eq. 9 under the assumption of ideal input difference operator $\nu=\frac{1}{2}(\rho-\hat{\rho})$ and output POVM element $Q$, i.e.

[TABLE]

where $\mathbf{P}$ is some pre-specified target Pauli matrix (fig. 1). Under these ideal conditions we can guarantee that the variance scales quadratically in the infidelity $r$. We will focus on intuition and relegate most technical work to the Supplementary Material. For the remainder of the text we will choose a basis for the space of linear operators $\mathcal{M}_{d}$. This means we can think of density matrices and POVM elements as column and row vectors which we denote with a Dirac-like notation, i.e. $\nu\rightarrow|\nu\rangle\!\rangle$ and $Q\rightarrow\langle\!\langle Q|$. Quantum channels can then be though of as matrices acting on vectors (which represent density matrices). Moreover, in this picture, composition of channels corresponds to matrix multiplication. When measuring the state $\mathcal{E}(\rho)$ using a two component POVM $\{Q,\mathbb{1}-Q\}$ for some quantum channel $\mathcal{E}$ and state $\rho$ and positive operator $Q$ we can write the expectation value $\operatorname{Tr}(Q\mathcal{E}(\rho)$ as a vector inner product

[TABLE]

where we abuse notation by referring to the matrix representation of the quantum channel $\mathcal{E}$ as $\mathcal{E}$ as well. This is variously called the affine or Liouville representation Wolf (2012); Wallman and Flammia (2014).

We assume that every experimental implementation of a Clifford gate $\mathcal{\tilde{G}}$ can be written as $\mathcal{\tilde{G}}=\mathcal{E}\mathcal{G}$ for some fixed CPTP map $\mathcal{E}$ where $G$ is the ideal Clifford gate. That is, we assume the noise is Markovian, constant and independent of the target gate. These assumptions can be relaxed partially Epstein et al. (2014); Wallman and Flammia (2014); Wallman and Emerson (2015); Wallman (2018).

The key to randomized benchmarking is that randomly applying elements of the Clifford group $\mathsf{C}$ and then inverting produces, on average, the depolarizing channel DiVincenzo et al. (2001)

[TABLE]

that is, we have

[TABLE]

with the depolarizing parameter $f$ related to the fidelity by Nielsen (2002)

[TABLE]

Therefore applying a sequence of independently-random gates and then inverting produces $\mathcal{D}_{f^{m}}$ on average. Hence the expectation value of any operator decays as $f^{m}$ on average.

The value of $k_{m}(\vec{G})$ for a fixed sequence of Clifford gates $\vec{G}$ (as defined in fig. 1), and the variance over $\vec{G}\in\mathbb{C}_{q}$ are

[TABLE]

respectively. We can use the identity $a^{2}=a\otimes a$ for $a\in\mathbb{C}$, the distributivity and associativity of the tensor product, and the linearity of quantum channels to write this as Magesan et al. (2011b); Wallman and Flammia (2014)

[TABLE]

where

[TABLE]

The superoperator $T_{\mathsf{C}}(\mathcal{E})$ is often referred to as the twirl of the quantum channel $\mathcal{E}$.

At this point, our analysis diverges from that of Ref. Wallman and Flammia (2014). First, note that for our modified scheme, $\nu^{\otimes 2}$ is traceless and symmetric under the interchange of the tensor factors (we will refer to such a matrix as a traceless symmetric matrix) so

[TABLE]

Furthermore, $T_{\mathsf{C}}(\mathcal{E}^{\otimes 2})$ preserves the trace and symmetry under interchange of tensor factors. Therefore we can define $T_{\rm{TS}}(\mathcal{E}^{\otimes 2})$ to be the restriction of $T_{\mathsf{C}}(\mathcal{E}^{\otimes 2})$ to the space of traceless symmetric matrices. As we prove in lemma 2 and Helsen et al. (2016), the representation $\mathcal{G}^{\otimes 2}$ of the Clifford group restricted to the traceless symmetric subspace decomposes into inequivalent irreducible representations. Therefore by Schur’s lemma (see Supplementary Materials for an explanation of Schur’s lemma),

[TABLE]

where $\mathcal{Z}$ indexes the irreducible subrepresentations of $\mathcal{G}^{\otimes 2}$ on the space of traceless symmetric matrices, $\mathcal{P}_{i}$ are projectors associated to each representation and $\chi_{i}=\chi_{i}(\mathcal{E})\in\mathbb{R}$ are numbers that depend on the quantum channel $\mathcal{E}$. 222When $\nu$ is not traceless, as is the case in regular randomized benchmarking, we can not restrict the two-copy twirl to a twirl over the traceless-symmetric subspace. However the derivation below will still hold, up to the addition of extra terms stemming from equivalent irreducible subrepresentations present in eq. 34. This extra term is discussed in section IV.9 and also Wallman and Flammia (2014) As the $\mathcal{P}_{i}$ are orthogonal projectors that span the space of traceless symmetric matrices, we can write the variance as

[TABLE]

Now we use a telescoping series trick (lemma 7 and in particular 1) on the last factor to write this as

[TABLE]

Here we see that getting a sharp bound on the variance will depend on getting sharp bounds on the difference between the $\chi_{i}$ prefactors and the square of the depolarizing parameter $f^{2}$. Before we start giving upper bounds to eq. 36, we would like to note that the behavior of eq. 36 with respect to the sequence length $m$ is very well matched to that of the final upper bounds given in eqs. 9 and 11. This justifies the use of eqs. 9 and 11 to set the weights in algorithm 1.

Up to this point the derivation has been valid for any input state difference $\nu$ and output positive operator $Q$. However now we will restrict to the case of ideal $Q$ and $\nu$. For the general case of non-ideal $Q$ and $\nu$ see the Supplementary materials. In the case of ideal $Q$ and $\nu$ we can use lemmas 3 and 4 to upper bound

[TABLE]

where $r=1-F_{\mathrm{avg}}(\mathcal{E},I)$ is the infidelity of the quantum channel. We would like to note here that lemma 3 can only be applied if $\nu$ ($Q$) are proportional to $\mathbf{P}$ ($\mathbb{1}+\mathbf{P}$). Moreover, we will see in section IV.7 that without this assumption the variance of RB will scale linearly in infidelity $r$. Continuing the calculation, for $r\leq\frac{1}{3}$, we can say that (lemma 6)

[TABLE]

Hence we can say

[TABLE]

for ideal $Q$ and $\nu$. Now we only need to deal with the $\chi_{i}$ factors in the sum. To do this we will use the fact that every $\chi_{i}$ term is upper bounded by the unitarity $u$ of the quantum channel $\mathcal{E}$. This is derived in lemma 5 in the Supplementary Material. Inserting this we get

[TABLE]

Now we factor $u^{m-2}$ out of the sum over $j$ and use the fact that this sum has a closed form. Using this and lemma 3 to bound the projector inner products we obtain a final bound on the variance

[TABLE]

which is the bound we set out to find. To obtain from this the bound given in eq. 10 we note that $u\geq f^{2}$ and moreover that the fractional term in eq. 42 is monotonically decreasing in $u$ (for fixed $f^{2}$) and reaches a limiting value of $m(m-1)/2$ in the limit of $u\rightarrow f^{2}$ (This can be seen by using l’H$\hat{\text{o}}$pital’s rule).

IV.6 State preparation and measurement

When $Q,\nu$ do not satisfy eq. 22, (which will always happen in practice) the above derivation will not hold exactly and the deviation of $Q,\nu$ from their ideal forms will introduce terms of order $\eta r$ i.e., terms which scale linearly and not quadratically in the infidelity $r$. Deriving an expression of the variance taking into account these these contributions is a little tedious so we will relegate it to the Supplementary Material and instead discuss the form of the prefactor $\eta$. Let $\nu$ be some non-ideal input state difference and let $Q$ be some non-ideal observable. Note from eq. 22 that the ideal input state difference $\nu$ and output POVM $Q$ are related to a pre-chosen “target Pauli matrix” $\mathbf{P}$. We hence have

[TABLE]

the ideal $Q$ and $\nu$. Suppressing some prefactors (the exact expression can be found in eq. 122 in the Supplementary material) we get the following approximate expression for the SPAM factor $\eta$:

[TABLE]

where $\left\lVert\cdot\right\rVert_{2}$ is the Schatten-$2$ norm Wolf (2012) and $Q,\nu$ are the non-ideal operators that are actually implemented. There are several important things to notice here:

• •

$\eta$ goes to zero in the limit of ideal $Q,\nu$. This justifies our choice of the ideal $Q$ and $\nu$ as being defined in terms of a single Pauli matrix rather than preparing and measuring in the $|0\rangle$ state as was the case in the original randomized benchmarking proposal Magesan et al. (2012a)

• •

$\eta$ scales quadratically in the deviation from the ideal of $Q$ and $\nu$. This means that for small deviations $\eta$ is likely to be small.

• •

$\eta$ is non-zero for non-ideal $Q$ even when $\nu$ is ideal and vice versa. This is unfortunate as it means that both state preparation and measurement must be good to ensure small variance. However, as we argue in section IV.7, this is actually optimal.

To get a feel for how the parameter $\eta$ behaves we discuss a particular error model for state preparation and measurement errors, inspired by recent research in superconducting qubits Sagastizabal et al. (2018). Here we see that the dominant error source when preparing states in the computational basis is given by decay to the ground state when in the excited ($|1\rangle$) state and residual excitations when preparing the ground ($|0\rangle$) state. The dominant contribution to measurement errors when measuring in the computational basis are here discrimination errors (mistaking [math] for $1$ and vice versa) as well as errors due to finite sampling. When performing our version of RB, and choosing $\mathbf{P}=Z$, we see that $\nu_{\mathrm{id}}=(|0\rangle\langle 0|-|1\rangle\langle 1|)/2$ and hence we want to ideally prepare the states $|0\rangle\langle 0|$ and $|1\rangle\langle 1|$. Following Sagastizabal et al. (2018) we assume $0.5\%$ residual excitations when preparing the $|0\rangle\langle 0|$ state, $0.8\%$ decay to the ground when preparing $|1\rangle\langle 1|$ and a $1\%$ discrimination error (modeled by a symmetric bit-flip channel) (Here we use the discrimination fidelity given in Riste et al. (2012)). Plugging these numbers into the assumed error models and calculating $\eta$ using eq. 122 in the Supplementary Material we see that in this case $\eta=0.001$. Hence we can say that under realistic scenarios $\eta$ will be quite small. It is possible to make a more fine-grained analysis of the SPAM term $\eta$ as it is defined under eq. 101, as opposed to upper bounding it. However this is likely to be rather involved and given that $\eta$ is already small in realistic scenarios we have opted not to pursue this here.

IV.7 Optimality of maximal variance

In this section we will argue that the bounds on the variance in the case of non-ideal SPAM are optimal in the sense that it is impossible for the variance to scale better that linearly in the infidelity $r$ for arbitrary noise maps when the input POVM element $Q$ is non-ideal even when the input state difference $\nu$ is ideal. The same reasoning will also hold for non-ideal $\nu$ even when $Q$ is ideal. (More generally the reasoning below will also work when randomized benchmarking is performed using a state rather than a state difference but we will not show this explicitly here).

Consider the variance as in eq. 30 for a randomized benchmarking experiment with a quantum channel $\mathcal{E}$ with infidelity $r$ and for simplicity set the sequence length $m=1$ (the argument will work for general $m$). Then we have an expression for the variance

[TABLE]

with the $T_{\mathsf{C}}(\mathcal{E}^{\otimes 2}),T_{\mathsf{C}}(\mathcal{E})^{\otimes 2}$ defined in eq. 31. Now consider setting $\nu=\nu_{\mathrm{id}}$ and maximizing over the POVM element $Q$. That is consider

[TABLE]

Now note that for any unitary $U$ the operator $\mathcal{U}(Q)=UQU^{\dagger}$ is also a POVM element. This means we can write

[TABLE]

where we used the linearity of the inner product and the definition of maximum and the integral is taken over the uniform or Haar measure of the unitary group. Now we use a well known fact from the representation theory of the unitary group which states that the integrated operator $\int dU(\mathcal{U}(Q))$ is precisely proportional to one of the projectors defined in eq. 34. Wallman et al. (2015a). In particular it is proportional to the rank one projector $P_{\mathrm{tr}}=|\Delta\rangle\!\rangle\langle\!\langle\Delta|$ where $\Delta\in\mathcal{M}_{d}$ is some matrix operator (see lemma 2 in the appendix) and $\mathrm{tr}$ is an element of the set $\mathcal{Z}$ which indexed the irreducible representations of the Clifford group in eq. 34. This means we can we can write using eq. 35

[TABLE]

where $\alpha(Q)$ is some positive prefactor function of $Q$. From lemma 5 and Wallman et al. (2015a) it can be seen that $\chi_{\mathrm{tr}}$ is precisely the unitarity $u$ of the quantum channel $\mathcal{E}$. If we now consider $\mathcal{E}$ to be a unitary channel (that is $u=1$), we get (ignoring the prefactors, which can be proven to be strictly positive)

[TABLE]

which is linear in infidelity $r$. Hence when the POVM element $Q$ is allowed to vary freely a linear scaling of the variance with the infidelity $r$ can not be avoided even when the input state difference $\nu$ is ideal. One can perform a similar thought experiment maximizing over $\nu$ while setting $Q=Q_{\mathrm{id}}$ and get the same result. Hence the expression for $\eta$ we discussed in the above section is essentially optimal.

IV.8 Asymptotic behavior of the variance

When looking at the bound on the variance eq. 9 the difference between unitary and non-unitary noise is striking. When the noise is non-unitary, and thus $u<1$ the upper bound on the variance (and hence the variance itself) decays exponentially to zero in the sequence length $m$ but when the noise process is unitary the variance keeps increasing and eventually saturates on a constant that is independent of the infidelity of the noise process. Here we argue that this is not an artifact of the bounding techniques but rather a fundamental feature of performing randomized benchmarking over unitary noise. Moreover this effect is independent of whether RB is performed using a state difference input $\nu$ or a state input $\rho$ (as in standard RB).

Consider a unitary noise process $\mathcal{U}=U\cdot U^{\dagger}$ with infidelity $r>0$ (That is $\mathcal{U}$ is not the identity). Now consider a randomized benchmarking experiment of sequence length $m$. That is, for a random sequence of Clifford unitaries $G_{1},\ldots G_{m}$ we perform the unitary

[TABLE]

Following the reasoning of Magesan et al. (2012a) we can write $V_{m}$ as

[TABLE]

where the unitaries $G^{\prime}_{m},\ldots G^{\prime}_{1}$ are sampled uniformly at random from the Clifford group. We can equally well think of the unitary $U^{\dagger}V_{m}$ as being the product of $m$ uniformly random samples from the set

[TABLE]

Note that this set depends on the unitary $U$. In Emerson et al. (2005b) it was shown that the distribution of the product of $m$ unitaries sampled uniformly at random from a set of unitaries converges to the Haar measure (uniform measure) on the unitary group in the limit of large $m$ as long as this set contains a universal set of gates. Note that this convergence phenomenon is independent of the initial set 333This is similar to how the limiting distribution of a random walk is independent of the initial step-size.

Note now that as long as the unitary $U$ is not a Clifford gate the set $\mathsf{G}_{U}$ will contain a universal gateset Chuang and Nielsen (1997). This means that the distribution from which $V_{m}$ is sampled will converge to the Haar measure in limit of long sequence length (the extra $U^{\dagger}$ factor gets absorbed into the Haar measure). This will happen independently of the unitary $U$ (as long as $U$ is not Clifford). From this we can conclude that the variance of randomized benchmarking with unitary noise must, in the limit of long sequences, converge to the variance of the randomized benchmarking expectation value over the Haar measure independently of what the original unitary noise process is. Note again that the above argument is independent of whether RB is performed using a state difference input or a state input.

IV.9 Relation to regular randomized benchmarking

When performing regular randomized benchmarking, that is using an input state $\rho=\frac{1}{2}(\mathbb{1}+\mathbf{P})$ rather than an input state difference $\nu=\frac{\mathbf{P}}{2}$ the upper bounds on the variance given in eqs. 11 and 9 still hold provided an extra additive term is added to them. This term will stem from the addition of an extra superoperator (that is not a projector) in the sum in eq. 34 which stem from the appearance of two equivalent trivial subrepresentations of the two-copy representation $\mathcal{G}^{\otimes 2}$ of the Clifford group. This term is of the form

[TABLE]

where $\mathcal{E}$ is the noise process under investigation, with infidelity $r$ and unitarity $u$ and system dimension $d$. Here $\left\lVert\mathcal{E}(\mathbb{1}/d)-\mathbb{1}/d\right\rVert_{2}^{2}$ is a measure of how ‘non-unital’ the quantum channel $\mathcal{E}$, that is how far its output deviates from the identity when the identity is the input. This measure can be upper bounded using (Wallman, 2015, Theorem 3) and is already implicitly analyzed in Wallman and Flammia (2014). We will not prove the above explicitly but it can be derived straightforwardly by following the derivation in theorem 1 using $\rho$ as input state. Note however that the upper bound on $T$ does not decay to zero exponentially but rather converges to a non-zero constant even for non-unitary channels. This is not a feature of the upper bound itself but rather of the long sequence behavior of standard randomized benchmarking. It was proven in (Wallman and Flammia, 2014, Theorem 17) that the upper bound $T$ is actually saturated for almost all non-unitary channels. Moreover, for physically relevant noise models such as amplitude damping $T$ can be quite substantial. This very different behavior in the limit of long sequence lengths further motivates the use the state difference $\nu$ for rigorous randomized benchmarking.

Acknowledgements.

We would like to thank Le Phuc Thinh, Jeremy Ribeiro and David Elkouss for enlightening discussions. We would also like to thank Bas Dirkse for pointing out several typos in an earlier version of this paper and diligently checking all proofs in the newest version. We would also like to thank Ramiro Sagastizabal and Adriaan Rol for their suggestions on the dominant SPAM errors in superconducting qubits. JH and SW are funded by STW Netherlands, NWO VIDI and an ERC Starting Grant. This research was supported by the U.S. Army Research Office through grant numbers W911NF-14-1-0098 and W911NF-14-1-0103. STF was supported by the Australian Research Council via EQuS project number CE11001013 and by an Australian Research Council Future Fellowship FT130101744.

Appendix A Preliminaries

A.1 Clifford and Pauli groups

In this section we recall definitions for the Pauli and Clifford groups on $q$ qubits. We begin by defining the Pauli group.

Definition 1** (Pauli group).**

Let $\{v_{0},v_{1}\}$ be an orthonormal basis of $\mathbb{C}^{2}$ and in this basis define the following linear operators by their action on the basis

[TABLE]

for $l\in\{0,1\}$ and addition over indices is taken modulo $2$. Note that $X,Y,Z\in U(2)$. The $q$-qubit Pauli group $\mathcal{P}_{q}$ is now defined as the subgroup of the unitary group $\mathrm{U}(2^{q})$ consisting of all $q$-fold tensor products of $q$ elements of $\mathcal{P}_{1}:=\langle X,Z,i\mathbb{1}_{2}\rangle$.

Elements $P,P^{\prime}$ of the Pauli group have the property that they either commute or anti-commute, that is

[TABLE]

We also define $\hat{\mathcal{P}}_{q}$ as the subset of $\mathcal{P}_{q}$ consisting of all $q$-fold tensor products of element of the set $\{\mathbb{1},X,Y,Z\}$, i.e.$\hat{\mathcal{P}}_{q}=\{\mathbb{1},X,Y,Z\}^{\otimes q}$. Note that the Hermitian subset $\hat{\mathcal{P}}_{q}$ of the Pauli group forms a basis for the Hilbert space $\mathcal{M}_{d}$. We can turn this into an orthonormal basis under the Hilbert-Schmidt inner product which is defined as

[TABLE]

To see this note that $\operatorname{Tr}(P)=0$ for all $P\in\mathcal{P}_{q}/\{\mathbb{1}\}$ and that $\operatorname{Tr}(\mathbb{1})=d$. We introduce the set of normalized Hermitian Pauli matrices.

[TABLE]

where we have given the normalized identity its own symbol for later convenience. We will denote the elements of the set $\bf{\sigma}_{q}$ by Greek letters ($\sigma,\tau,\nu,...$). We also, for later convenience, introduce the normalized matrix product of two normalized Pauli matrices as

[TABLE]

Note that $\sigma\cdot\tau\in\pm\bf{\sigma}_{q}\cup\sigma_{0}$ if $[\sigma,\tau]=0$ and $i\sigma\cdot\tau\in\pm\bf{\sigma}_{q}$ if $\{\sigma,\tau\}=0$. Lastly we define the following parametrized subsets of $\bf{\sigma}_{q}$ .For all $\tau\in\bf{\sigma}_{q}$ we define

[TABLE]

Note that we have $|{\bf{N}}_{\tau}|=\frac{d^{2}}{2}$, $|{\bf{C}}_{\tau}|=\frac{d^{2}}{2}-2$ and ${\bf{C}}_{\tau}$ and ${\bf{N}}_{\tau}$ are disjoint for all $\tau\in\bf{\sigma}_{q}$. We also have for $\sigma,\sigma^{\prime}\in\bf{\sigma}_{q}$ and $\sigma\neq\sigma^{\prime}$ that $|{\bf{C}}_{\sigma}\cap{\bf{C}}_{\sigma^{\prime}}|=\frac{d^{2}}{4}-3$. For a proof of this see (Helsen et al., 2016, Lemma 1).

Next we define the Clifford group. We have

Definition 2**.**

The $q$-qubit Clifford group $\mathsf{C}_{q}$ is the normalizer (up to complex phases) of $\mathcal{P}_{q}$ in $\mathrm{U}(2^{q})$, that is,

[TABLE]

The Clifford group is also often introduced as the group generated by the Hadamard (H), $\pi/4$ phase gate and CNOT gates on all qubits. These are equivalent definitions (up to global phases) Farinholt (2014).

For a more expansive introduction to the Pauli and Clifford groups see e.g. Farinholt (2014) and references therein.

A.2 Representation theory

We recall some useful facts about the representations of finite groups. For a more in depth treatment of this topic we refer to Fulton and Harris (2004); Goodman and Wallach (2009). Let $\mathsf{G}$ be a finite group and let $V$ be some finite dimensional complex vector space. Let also $\mathsf{GL}(V)$ be the group of invertible linear transformations of $V$. We can define a representation $\phi$ of the group $G$ on the space $V$ as a map

[TABLE]

that has the property

[TABLE]

In general we will assume the operators $\phi(g)$ to be unitary. If there is a non-trivial subspace $W$ of $V$ such that

[TABLE]

then the representation $\phi$ is called reducible. The restriction of $\phi$ to the subspace $W$ is also a representation, which we call a subrepresentation of $\phi$. If there are no non-trivial subspaces $W$ such that eq. 61 holds the representation $\phi$ is called irreducible. Two representations $\phi,\phi^{\prime}$ of a group $G$ on spaces $V,V^{\prime}$ are called equivalent if there exists an invertible linear map $T:V\rightarrow V^{\prime}$ such that

[TABLE]

We can also define the twirl $\mathcal{T}_{\phi}(A)$ of a linear map $A:V\to V$ with respect to the representation $\phi$ to be

[TABLE]

The following corollary of Schur’s lemma, an essential result from representation theory. Fulton and Harris (2004); Goodman and Wallach (2009), allows us to evaluate twirls over certain types of representations.

Lemma 1**.**

Let $\mathsf{G}$ be a finite group and let $\phi$ be a representation of $\mathsf{G}$ on a complex vector space $V$ with decomposition

[TABLE]

into inequivalent irreducible subrepresentations $\phi_{i}$. Then for any linear operator $A$ from $V$ to $V$, the twirl of $A$ over $\mathsf{G}$ takes the form

[TABLE]

where $P_{i}$ is the projector onto the subspace carrying the irreducible subrepresentation $\phi_{i}$. In the rest of the text we will often denote the prefactor $\operatorname{Tr}(AP_{i})/\operatorname{Tr}(P_{i})$ by $\chi_{i}$.

A.3 Liouville representation of quantum channels

Quantum channels Chuang and Nielsen (1997) are completely positive and trace-preserving (CPTP) linear maps $\mathcal{E}:\mathcal{M}_{d}\rightarrow\mathcal{M}_{d}$. We will denote quantum channels by calligraphic font throughout. The canonical example of a quantum channel is conjugation by a unitary $U$, which we denote by the corresponding calligraphic letter, i.e., $\mathcal{U}(\rho)=U\rho U^{\dagger}$ for all density matrices $\rho$. We will denote the noisy implementation of a channel by an overset tilde, e.g., $\tilde{\mathcal{G}}$ denotes a noisy implementation some ideal quantum channel channel $\mathcal{G}$.

It is often useful to think of quantum channels as matrices acting on vectors. This is variously known as the Liouville Wallman and Flammia (2014) or affine Wolf (2012) representation. This representation corresponds to fixing an orthonormal basis for $\mathcal{M}_{d}$ according to the Hilbert-Schmidt or trace-inner product and then expressing elements of $\mathcal{M}_{d}$ as vectors in $\mathbb{C}^{d^{2}}$. The Hilbert-Schmidt inner product is again defined as

[TABLE]

Now let $\{B_{j}\}_{j}$ for $j\in\mathbb{Z}_{d^{2}}$ be an orthonormal basis for $\mathbb{C}^{d\times d}$ with respect to the Hilbert-Schmidt inner product. We can construct a map $|.\rangle\!\rangle:\mathcal{M}_{d}\to\mathbb{C}^{d^{2}}$ by setting $|B_{j}\rangle\!\rangle=e_{j}$ where $e_{j}$ is the $j$th canonical basis vector for $\mathbb{C}^{d^{2}}$. Linearly extending the map $|\cdot\rangle\!\rangle$ to all elements $M\in\mathcal{M}_{d}$ we get

[TABLE]

Defining $\langle\!\langle M|=|M\rangle\!\rangle^{\dagger}$, we then have

[TABLE]

so that the Hilbert-Schmidt inner product is equivalent to the standard vector inner product.

We will generally construct the Liouville representation using the basis spanned by the normalized (with respect to the Hilbert-Schmidt inner product) Pauli matrices $\{\sigma_{0}\}\cup\boldsymbol{\sigma}_{q}$ where $\sigma_{0}:=I_{d}/\sqrt{d}$ with $d=2^{q}$ is the normalized identity matrix and

[TABLE]

is the set of normalized Hermitian Pauli matrices excluding the identity.

As any quantum channel $\mathcal{E}$ is a linear map from $\mathcal{M}_{d}$ to $\mathcal{M}_{d}$ we have

[TABLE]

so that we can represent $\mathcal{E}$ by the matrix

[TABLE]

where we abuse notation by using the same symbol to refer to an abstract channel and its matrix representation. The action of a channel $\mathcal{E}$ on a density matrix $\rho$ now corresponds to the standard matrix action on the vector $|\rho\rangle\!\rangle$, hence for a density matrix $\rho$ and a POVM element $Q$ in $\mathcal{M}_{d}$ we have

[TABLE]

The Liouville representation has the nice properties (as can be easily checked) that the composition of quantum channels is equivalent to matrix multiplication of their Liouville matrices and that tensor products of channels correspond to tensor products of the corresponding Liouville matrices, that is, for all channels $\mathcal{E}_{1}$ and $\mathcal{E}_{2}$ and all $A\in\mathcal{M}_{d}$,

[TABLE]

In the Liouville picture the depolarizing parameter and the unitarity Wallman et al. (2015a) of a quantum channel $\mathcal{E}$ are

[TABLE]

and the Liouville representation of a depolarizing channel with depolarizing parameter $f$ is given by Wallman and Flammia (2014)

[TABLE]

A.4 Traceless-Symmetric representation

In the rest of the text we will often work with quantum channels which have a tensor product structure. That is we will often be dealing with channels of the form

[TABLE]

where $\mathcal{E}_{i}$ is a CPTP map for all $i$ and $\lambda_{i}\in\mathbb{C}$ is some abstract parameter. Note that $\mathcal{W}$ is now a linear map from $\mathcal{M}_{d}^{\otimes 2}$ to $\mathcal{M}_{d}^{\otimes 2}$. Maps of these form have a number of useful properties which we will now consider. We begin by defining the traceless-symmetric subspace $V_{\mathrm{TS}}$ which is a subspace of $\mathcal{M}_{d}^{\otimes 2}$ of the form

[TABLE]

where we have suppressed the tensor product (that is $\sigma\tau:=\sigma\otimes\tau$). The traceless-symmetric subspace has several desirable properties which we note here. First let $\rho,\hat{\rho}\in\mathcal{M}_{d}$ be density matrices and call their difference $\nu:=\rho-\hat{\rho}$, then we have that

[TABLE]

Moreover, for any quantum channel $\mathcal{W}$ of the form defined in eq. 78 we have that

[TABLE]

or equivalently we have that

[TABLE]

where $\mathcal{P}_{\mathrm{TS}}$ is the projector onto the space $V_{\mathrm{TS}}$ (note that $\mathcal{P}_{\mathrm{TS}}$ is a linear map from $\mathcal{M}_{d}^{\otimes 2}$ to $\mathcal{M}_{d}^{\otimes 2}$). This observation follows from the fact than $\mathcal{W}$ is a linear combination of two-fold tensor products of quantum channels (which preserve the trace and map operators that are symmetric under interchange of the two copies of $\mathcal{M}_{d}^{\otimes 2}$ to operators that are symmetric under interchange of the two copies of $\mathcal{M}_{d}^{\otimes 2}$).

We will in particular be interested in how a representation of of the Clifford group $\mathsf{C}$ behaves on the traceless symmetric subspace. Define the two-fold tensor product representation of the Clifford group on $\mathcal{M}_{d}^{\otimes 2}$ as

[TABLE]

for all where $\mathcal{G}$ is the Liouville representation of $G$ for all $G\in\mathsf{C}$. This representation has a natural restriction to the subspace $V_{\mathrm{TS}}$ since $\mathcal{G}^{\otimes 2}$ is of the form described in eq. 78. We can define the subrepresentation $\phi_{\mathrm{TS}}$ of $\phi_{2}$ as

[TABLE]

for all $G\in\mathsf{C}$. This representation is in general not irreducible but decomposes further into a collection of irreducible subrepresentations. In Helsen et al. (2016) we derived these irreducible subrepresentations of $\phi_{\mathrm{TS}}$ and studied their properties. In the following lemma we will quote several results from Helsen et al. (2016) which will be useful for our purposes.

Lemma 2**.**

Let $\mathsf{C}$ be the Clifford group and let $\phi_{\mathrm{TS}}$ be the traceless symmetric representation. This representation is a direct sum of three subrepresentations $\phi_{\mathrm{d}}$ (diagonal), $\phi_{\mathrm{[S]}}$ (symmetric commuting) and $\phi_{\mathrm{\{S\}}}$ (symmetric anti-commuting) acting on the spaces

[TABLE]

The diagonal subrepresentation $\phi_{\mathrm{d}}$ decomposes into three subrepresentations denoted by $\phi_{\mathrm{tr}},\phi_{1},\phi_{2}$ with $\phi_{\mathrm{tr}}$ the trivial representation spanned by

[TABLE]

We will index these representations by the set $\mathcal{Z}_{\mathrm{d}}:=\{\mathrm{tr},1,2\}$.

The symmetric commuting representation $\phi_{[S]}$ decomposes into $3$ irreducible subrepresentations denoted as $\phi_{[\mathrm{adj}]},\phi_{[1]},\phi_{[2]}$. We will index these representations by the set $\mathcal{Z}_{[S]}:=\{[\mathrm{adj}],[1],[2]\}$. The spaces carrying these representations can be written as a direct sum of subspaces in the following way

[TABLE]

where $V_{i}^{\tau}\subset V^{[\tau]}$ with

[TABLE]

The symmetric anti-commuting representation $\phi_{\{S\}}$ decomposes into $2$ irreducible subrepresentations denoted as $\phi_{\{1\}},\phi_{\{2\}}$. We will index these representations by the set $\mathcal{Z}_{\{S\}}:=\{\{1\},\{2\}\}$. The spaces carrying these representations can be written as a direct sum of subspaces in the following way

[TABLE]

where $V_{i}^{\tau}\subset V^{\{\tau\}}$ with

[TABLE]

Finally we denote the set indexing all irreducible subrepresentations of $\phi_{\mathrm{TS}}$ as $\mathcal{Z}=\mathcal{Z}_{\mathrm{d}}\cup\mathcal{Z}_{[S]}\cup\mathcal{Z}_{\{S\}}$ and we note that all irreducible representations indexed by $\mathcal{Z}$ are mutually inequivalent.

Note that we have only given an explicit basis for the space on which the representation $\phi_{\mathrm{tr}}$ acts. It is possible to write down explicit bases for all relevant vector spaces but we will not need to do see (see however Helsen et al. (2016)).

Appendix B Randomized benchmarking

B.1 Variance bound

In this section we prove the main theorem of the paper. Concretely we prove the following.

Theorem 1**.**

Let $Q$ be an observable and $\rho,\hat{\rho}$ density matrices and set $\nu=\frac{1}{2}(\rho-\hat{\rho})$. Consider a randomized benchmarking experiment using the Clifford group $\mathsf{C}$ with noisy implementation $\mathcal{\tilde{G}}=\mathcal{E}\mathcal{G}$ for all $G\in\mathsf{C}$. Then the variance $\mathbb{V}^{2}_{m}$ of this experiment is upper bounded by

[TABLE]

where $u=u(\mathcal{E})$ is the unitarity, $r=r(\mathcal{E})$ is the infidelity, $d$ is the system dimension, $m$ is the sequence length, $f=1-\frac{dr}{d-1}$ is the depolarizing parameter and $\eta$ is a function capturing the deviation from the ideal $Q$ and $\nu$. This bound is valid for $r\leq\frac{1}{3}$.

Proof.

We begin from an exact expression of the variance expressed in the Liouville representation eq. 30:

[TABLE]

where $\mathcal{T}_{\phi_{2}}$ is the twirl over the two-copy representation of the Clifford group as defined in eq. 83 and $\mathcal{T}_{\phi}$ is the twirl over the (single copy) Liouville representation. Note now that $|\nu^{\otimes 2}\rangle\!\rangle\in V_{\mathrm{TS}}$ and that both $\mathcal{T}_{\phi_{2}}(\mathcal{E}^{\otimes 2})$ and $\mathcal{T}_{\phi}(\mathcal{E})^{\otimes 2}$ are CPTP maps of the form described in eq. 78. This means we can restrict both twirls to the traceless symmetric subspace. In this subspace we have from lemma 1 and lemma 2 that $\mathcal{T}_{\phi_{2}}(\mathcal{E}^{\otimes 2})$ and $\mathcal{T}_{\phi}(\mathcal{E})^{\otimes 2}$ are of the form

[TABLE]

where $\mathcal{Z}$ (as defined in lemma 2) indexes the irreducible subrepresentations of the traceless symmetric representation of the Clifford group and $\chi_{i}=\operatorname{Tr}(\mathcal{P}_{i}\mathcal{E}^{\otimes 2})/\operatorname{Tr}(\mathcal{P}_{i})$ are the prefactors associated to the different subrepresentations. We also used that $\mathcal{T}_{\phi}(\mathcal{E})$ is a depolarizing channel with depolarizing parameter $f$ Wallman and Flammia (2014). Using that $\mathcal{P}_{i}^{2}=\mathcal{P}_{i}$ and $\mathcal{P}_{i}\mathcal{P}_{j}=0$ for $i,j\in\mathcal{Z},i\neq j$ we can rewrite the variance as

[TABLE]

We now apply a telescoping series identity, which is proven in 1 of lemma 7, to the factor $\chi_{i}^{m}-f^{2m}$ in the above equation (for all $i\in\mathcal{Z}$). This gives

[TABLE]

This equation contains two terms, eq. 94a and eq. 94b which we will bound separately. We now proceed to upper bound the first term, that is eq. 94a. For this we will split the the input and output operators $Q,\nu$ into their ideal parts (that is, the Pauli operator $\sigma_{\mathbf{P}}:=\mathbf{P}/\sqrt{d}$) and deviations from that ideal. We define the functions

[TABLE]

for all $i\in\mathcal{Z}$ where $Q_{\mathbf{P}}=\operatorname{Tr}(Q\sigma_{\mathbf{P}})$ and similarly for $\nu_{\mathbf{P}}$. Using this we can write eq. 94a as

[TABLE]

Now consider the first term of the RHS, eq. 96a. First note from lemma 3 that for $i\not\in\mathcal{Z}_{\mathrm{d}}=\{\mathrm{tr},1,2\}$ we have $\mathcal{P}_{i}|\sigma_{\mathbf{P}}^{\otimes 2}\rangle\!\rangle=0$. Hence we have

[TABLE]

where we used lemma 3 in the first and second equalities and the fact that

[TABLE]

in the last equality (this can be seen from lemma 2). Now we use lemma 4 and the fact that $Q_{\mathbf{P}}\nu_{\mathbf{P}}\leq 1/4$ to obtain an upper bound

[TABLE]

This leaves us with the second term in the RHS, eq. 96b. Here we cannot attain a bound that is quadratic in $r$. Instead we will attempt a bound that is linear in $r$ using lemma 6. We can write

[TABLE]

subject to the condition $r\leq\frac{1}{3}$. Writing $\eta(Q,\nu):=\sum_{i\in\mathcal{Z}}|H_{i}(Q,\nu)|$ we have a bound on eq. 94a.

We continue by upper bounding the second term in the variance, that is eq. 94b. We again split off the ideal components of $Q$ and $\nu$ and write

[TABLE]

where we have used the definition of the function $H_{i}(Q,\nu)$, lemma 3 and the triangle inequality. Now we use lemma 6 to upper bound this quantity as

[TABLE]

where we have used the fact that $\sum_{i\in\mathcal{Z}_{\mathrm{d}}}\operatorname{Tr}(\mathcal{P}_{i})=d^{2}-1$. It remains to deal with the last factor. This we do by using lemma 5 which states that $\chi_{i}\leq u$ for all $i\in\mathcal{Z}$, where $u$ is the unitarity of the channel $\mathcal{E}$. Writing again $\eta(Q,\nu):=\sum_{i\in\mathcal{Z}}|H_{i}(Q,\nu)|$ we then have

[TABLE]

We can further make sense of this quantity by using the well known series identity

[TABLE]

Factoring out a factor of $u^{m-2}$ and setting $x=f^{2}/u$ we obtain the following

[TABLE]

This finishes the upper bounding of eq. 94b. Gathering all terms we come to a final bound

[TABLE]

which is the bound we set out to find. ∎

Noting that $f^{2}\leq u$ and that the factor

[TABLE]

is monotonically decreasing in $u$ we can upper bound this factor by taking the limit $u\rightarrow f^{2}$. This gives

[TABLE]

which can be confirmed by an application of l’H$\hat{\text{o}}$pital’s rule. Plugging this in to eq. 106 we obtain eq. 10.

B.2 State preparation and measurement (SPAM) terms

In the central bound on the variance ( theorem 1) we had to account for the fact that the variance can depend on how well the input states $\rho,\hat{\rho}$ and the output POVM $Q$ can be implemented. The ideal behavior of $\nu=\frac{1}{2}(\rho-\hat{\rho})$ and $Q$ are given by

[TABLE]

where $\mathbf{P}$ is a pre-specified element of the Pauli group (see fig. 1). The deviation of $Q$ and $\nu$ from this ideal can be captured by writing

[TABLE]

where $\langle Q_{\mathrm{id}},Q_{\mathrm{spam}}\rangle=\langle\nu_{\mathrm{id}},\nu_{\mathrm{spam}}\rangle=0$.

In the variance bound the deviation from the ideal has an effect which is measured by the parameter $\eta(Q,\nu)$. This parameter $\eta(Q,\nu)$ was defined as

[TABLE]

where $\mathcal{Z}$ indexes the irreducible representations of the traceless symmetric representation of the Clifford group and the $\mathcal{P}_{i}$ are projectors onto the spaces carrying these subrepresentations (lemma 2). Let us now analyze these terms further. For $i\in\mathcal{Z}_{\mathrm{d}}$ we have

[TABLE]

where we have used that $\langle Q_{\mathrm{id}},Q_{\mathrm{spam}}\rangle=\langle\nu_{\mathrm{id}},\nu_{\mathrm{spam}}\rangle=0$ which implies that $\langle\!\langle Q_{\mathrm{id}}\otimes Q_{\mathrm{spam}}|\mathcal{P}_{i}=\mathcal{P}_{i}|\nu_{\mathrm{id}}\otimes\nu_{\mathrm{spam}}\rangle\!\rangle=0$ for $i\in\mathcal{Z}_{\mathrm{d}}$. Using the triangle inequality and the Cauchy-Schwarz inequality we can get

[TABLE]

where $\left\lVert\mathcal{P}_{i}\right\rVert_{2\rightarrow 2}$ is the induced $2$-norm of the superoperator $\mathcal{P}_{i}$. It is well known that this norm is equal to the largest singular value of the Liouville representation of $\mathcal{P}_{i}$ Wallman and Flammia (2014), which since the Liouville representation of $\mathcal{P}_{i}$ is an orthonormal projection, is equal to one. This means we have for $i\in\mathcal{Z}_{\mathrm{d}}$ that

[TABLE]

Note that this expression is zero when both $Q$ and $\nu$ are ideally implemented but is non-zero when either of them is not. This behavior is in general unavoidable as we argue in the main text (section IV.7). But first we will consider the functions $H_{i}(Q,\nu)$ for $i\in Z_{[S]}\cup Z_{\{S\}}$. Note first that since $\mathrm{supp}(\mathcal{P}_{i})\subset\mathrm{Span}\{S_{\sigma,\sigma^{\prime}}\;\|\;\;\sigma,\sigma^{\prime}\in\bf{\sigma}_{q},\;\sigma\neq\sigma^{\prime}\}$ we must have that $\mathcal{P}_{i}|\nu_{\mathrm{id}}^{\otimes 2}\rangle\!\rangle=\langle\!\langle Q_{\mathrm{id}}^{\otimes 2}|\mathcal{P}_{i}=0$. This means we can write

[TABLE]

which we can rewrite as

[TABLE]

which makes manifest that $H_{i}(Q,\nu)=0$ if $Q$ and $\nu$ are ideal and moreover that this term actually scales with the product of the deviations in $Q$ and $\nu$ (as measured in the $2$-norm). Hence we see that to lowest order in $Q_{\mathrm{spam}}$ and $\nu_{\mathrm{spam}}$ the SPAM parameter $\eta(Q,\nu)$ is proportional to

[TABLE]

with the exact expression being

[TABLE]

where the factors $3$ and $5$ arise from the fact that $|\mathcal{Z}_{\mathrm{d}}|=3$ and $|Z_{[S]}\cup Z_{\{S\}}|=5$ respectively (this is for $q\geq 3$, for $q=1$ we get the significantly better $|\mathcal{Z}_{\mathrm{d}}|=2$ and $|Z_{[S]}\cup Z_{\{S\}}|=1$ instead Helsen et al. (2016)).

B.3 Sample complexity of iteratively reweighted least squares

In this section we analyze the sample complexity of the RB fitting procedure using iteratively reweighted least squares, as outlined in section IV.3. Given a set of sequence lengths $\mathsf{M}$ we will assume that $N$ random sequences are sampled for each sequence length. It is possible to let $N$ be a function of the sequence length $m$ and prove a more general version of the theorem presented here but we will not pursue this here. We will also only be interested in the uncertainty around the estimate for the depolarizing parameter $f$, it is straightforward to extend our analysis to also include the uncertainty around estimate for the pre-factor $A$. The methods we use are all standard and can be found in Fedorov (2013); Seber and Wild (1989). See also Epstein et al. (2014) for an earlier calculation of this form in the context of randomized benchmarking (not taking into account the heteroskedasticity of randomized benchmarking data).

Theorem 2**.**

Let $\mathsf{M}$ be a set of integers denoting sequence lengths and let $\{k_{m,N}\}_{m\in\mathsf{M}}$ be a set of RB data points obtained by sampling $N$ random sequences for each sequence length $m\in\mathsf{M}$. Denote by $f^{*},A^{*}$ the true values for the RB fitting parameters and denote by $f_{\rm est},A_{\rm est}$ their estimates as obtained by the iteratively reweighted least squares procedure outlined in algorithm 1. We then have that

[TABLE]

where $\delta$ is upper bounded by

[TABLE]

with $H$ defined in eq. 15 and

[TABLE]

and

[TABLE]

and $J_{1}$ is the first entry of this vector.

Proof.

The starting off point for this proof is given by Eq. 1.6.19 in (Fedorov, 2013, Page 45) which states that the outcome of the IRLS procedure satisfies the following equality

[TABLE]

where $w(f,m)$ is the weight function given by the inverse of eq. 30 (we suppress the dependency on the unitarity here for notational simplicity). We can rewrite eq. 128 as

[TABLE]

We can now think of the LHS of eq. 130 as a function of the vector $[f_{\rm est},A_{\rm est}]$. Assuming $[f_{\rm est},A_{\rm est}]$ is close to $[f^{*},A^{*}]$ we can expand the LHS of eq. 130 to first order to get

[TABLE]

where $J$ is the Jacobian associated to the LHS of eq. 130, that is:

[TABLE]

Taking the Moore-Penrose inverse $J^{\rm MP}=(J^{T}J)^{-1}J^{T}$ of $J$ and inserting this in the first entry of eq. 131 we can say that

[TABLE]

where $J_{1}$ is the first entry of $J$. Now we can say that

[TABLE]

Now note that $k_{m,N}$ can be seen as a number drawn from a random variable $K_{m}$ with mean $A^{*}(f^{*})^{m}$ and variance $\mathbb{V}_{m}(f^{*})/N^{2}$ where $N$ is the number of random sequences drawn for each data-point $k_{m,N}$. Moreover $k_{m,N}$ and $k_{N,m^{\prime}}$ for $m\neq m^{\prime}$ are drawn from independent random variables $K_{m}$ and $K_{m^{\prime}}$. Hence we can apply the concentration inequality given in eq. 15 to eq. 135 to get

[TABLE]

with $\mathbb{V}_{\mathrm{fit}},\epsilon_{\mathrm{fit}}$ given by

[TABLE]

which completes the proof. ∎

Using eq. 9 or eq. 11 then gives an upper bound on total amount of data that needs to be gathered for rigorous RB.

Appendix C Technical lemmas

In this section we give proofs of all technical lemmas used in the main result theorem 1.

C.1 Projectors in the traceless symmetric subspace

In lemma 3 we prove a series of useful upper bounds on the trace overlap between the superoperator-projectors associated to the traceless-symmetric representation of the Clifford group and the normalized Pauli matrices. The saturated versions of these inequalities are critical to establishing the quadratic scaling with infidelity of the variance bound in the case of SPAM-free RB.

Lemma 3**.**

Let $\mathcal{E}:\mathcal{M}_{d}\to\mathcal{M}_{d}$ be a quantum channel and consider the twirled operator $\mathcal{T}_{\phi_{\mathrm{TS}}}(\mathcal{E}^{\otimes 2})$ with respect to the traceless-symmetric representation. This operator can then be written as (lemmas 2 and 1)

[TABLE]

with $\mathcal{Z}=\{\mathrm{tr},1,2,[1],[2],[3],\{1\},\{2\}\}$ and $\mathcal{P}_{i}$ the projector onto the spaces $V_{i}\subset\mathcal{M}_{d}^{\otimes 2}$. Let $I(x\in A)$ be the indicator function for the set $A$ (that is $I(x\in A)=1$ if $x\in A$ and $I(x\in A)=0$ otherwise). We have the following statements

• •

For $i\in\mathcal{Z}$ and $\sigma,\sigma^{\prime}\in\bf{\sigma}_{q}$ we have that

[TABLE]

with equality when $\sigma=\sigma^{\prime}$ .

• •

For $i\in\mathcal{Z}$, $\tau,\tau^{\prime}\in\bf{\sigma}_{q}$ and $\sigma\in\bf{C}_{\tau},\sigma^{\prime}\in\bf{C}_{\tau^{\prime}}$ we have that

[TABLE]

with equality when $\sigma=\sigma^{\prime}$.

• •

For $i\in\mathcal{Z}$, $\tau,\tau^{\prime}\in\bf{\sigma}_{q}$ and $\sigma\in\bf{N}_{\tau},\sigma^{\prime}\in\bf{N}_{\tau^{\prime}}$ we have that

[TABLE]

with equality when $\sigma=\sigma^{\prime}$.

where the sets $Z_{\mathrm{d}},\mathcal{Z}_{[S]},\mathcal{Z}_{\{S\}}$ are defined in lemma 2.

Proof.

We begin by proving the first claim. Let $\mathcal{P}_{i}$ be a projector as defined in the lemma statement with $i\in\mathcal{Z}$ and take $\sigma,{\sigma^{\prime}}\in\bf{\sigma}_{q}$. From lemma 2 we have immediately that

[TABLE]

Now consider $i\in\mathcal{Z}_{\mathrm{d}}$. Note that since $\mathcal{P}_{i}$ is a projector it is a real matrix and we have that $\mathcal{P}_{i}\geq 0$, that is $\mathcal{P}_{i}$ is a positive semidefinite matrix. This means that we have, by the Sylvester principal minor conditions, that

[TABLE]

for all $\sigma,\sigma^{\prime}\in\bf{\sigma}_{q}$. Now consider the case $\sigma=\sigma^{\prime}$. Note that for all $\tau,\sigma\in\bf{\sigma}_{q}$ there is a $G_{\tau}^{\sigma}\in\mathsf{C}$ such that $\mathcal{G}_{\tau}^{\sigma}(\tau)=\pm\sigma$. That is, the Clifford group acts transitively on $\bf{\sigma}_{q}$ Zhu (2015). This means we can write

[TABLE]

where we used the fact that $\mathcal{P}_{i}$ commutes with $\mathcal{G}^{\otimes 2}$ for all $G\in\mathsf{C}$ and the fact that $V_{i}\subset V_{\mathrm{d}}$ (where $V_{\mathrm{d}}$ is defined in lemma 2). This proves the first claim of the lemma.

Next we consider the second claim of the lemma. Let $\tau,\tau^{\prime}\in\bf{\sigma}_{q}$ and take $\sigma\in\bf{C}_{\tau}$ and $\sigma^{\prime}\in\bf{C}_{\tau^{\prime}}$. Again from lemma 2 we have immediately that

[TABLE]

Now consider $i\in\mathcal{Z}_{[S]}$. From lemma 2 we have that we can write

[TABLE]

where $\mathcal{P}^{\tau}_{i}$ has support in the space

[TABLE]

From this we immediately get

[TABLE]

Now consider $\tau=\tau^{\prime}$. Again from the Sylvester minor conditions we get for all $\sigma,\sigma^{\prime}\in\bf{C}_{\tau}$ that

[TABLE]

Now consider the case $\sigma=\sigma^{\prime}$. From Zhu (2015) we can see that the action of the Clifford group on the set $A=\{(\sigma,\sigma\cdot\tau)\;\|\;\;\tau\in\bf{\sigma}_{q},\;\sigma\in\bf{C}_{\tau}\}$ is $2$-transitive. That is, for all pairs $(\nu,\mu)\in A$ there is a $G_{\nu,\mu}^{\sigma,\tau}\in\mathsf{C}$ such that

[TABLE]

This implies we can make essentially the same argument as before, that is

[TABLE]

where we have used the fact that $\mathcal{G}^{\otimes 2}$ commutes with $\mathcal{P}_{i}$ for all $G\in\mathsf{C}$ and also the definition of the space $V_{[S]}$ (given in lemma 2). The factor of two appears from the fact that the set $A$ counts the basis of $V_{[S]}$ twice since $S_{\nu,\nu\cdot\mu}=S_{\nu\cdot\mu,\nu}$ for all $(\mu,\nu\cdot\mu)\in A$. We have also used that $|A|=|{\bf{\sigma}_{q}}||{\bf{C}_{\tau}}|=(d^{2}-1)(d^{2}/2-2)$. This proves the second claim of the lemma.

The proof of the third claim of the lemma proceeds in the same way as the proof of the second claim with the difference that anti-commuting, rather than commuting elements of the Pauli group must considered. We will not write it down explicitly. ∎

C.2 Bound on sum of squares of the diagonal elements of a quantum channel

This lemma (lemma 4) proves that the diagonal elements of a CPTP map are generically quite close to their mean. The key technique used here is the fact that the diagonal elements of a CPTP map are invariant under Pauli twirling. This is a structural result about quantum channels on arbitrarily many qubits and might be of independent interest. We use it to establish the quadratic scaling of the variance in the infidelity in the case of SPAM-free RB.

Lemma 4**.**

Let $\mathcal{E}:\mathcal{M}_{d}\to\mathcal{M}_{d}$ be a quantum channel with infidelity $r$ and depolarizing parameter $f=1-\frac{dr}{d-1}$. The quantity

[TABLE]

where $\mathcal{E}_{\tau,\tau}=\langle\tau,\mathcal{E}(\tau)\rangle$, has the following upper and lower bounds in terms of the infidelity $r$

[TABLE]

Proof.

We begin by noting that upper and lower bounds of the quantity eq. 153 can be found by maximizing and minimizing respectively the following optimization

[TABLE]

Here we maximize (minimize) the quantity eq. 153 over all possible CPTP maps which have depolarizing parameter $f$ . Solving this optimization problem is not easy since it not clear how to express the CP condition in terms of the optimization parameters $\mathcal{E}_{\tau\tau}$. We will therefore relax this problem to an easier one which we can solve. We begin by noting that the optimization variables $\mathcal{E}_{\tau\tau}$ are invariant under the action of a Pauli channel, i.e. for all $G\in\mathsf{P}$ with $\mathsf{P}$ the Pauli group, we have that

[TABLE]

for all $\tau\in\boldsymbol{\sigma}_{q}\cup\sigma_{0}$ where $\text{sgn}(\tau,G)$ is defined as

[TABLE]

which, since $\tau\in\boldsymbol{\sigma}_{q}\cup\sigma_{0}$ is a normalized element of the Pauli group, is well defined because elements of the multi-qubit Pauli group can either commute $([.,.])$ or anti-commute $(\{.,.\})$ with each other Gottesman (1998). By eq. 156 and linearity we can now note that the optimization variables in the optimization eq. 155 are invariant under twirling over the Pauli group $\mathsf{P}$, i.e.

[TABLE]

Note also that the “twirl” operation, for any group, preserves complete positivity Wolf (2012). This means we can relax the optimization eq. 155 to

[TABLE]

Note that this is a relaxation of the previous optimization because while the twirl of a CP map will always be CP the opposite need no be true. Now we use the following result due to Holevo Holevo (2005) which states that any CPTP map $\mathcal{E}$, twirled over the Pauli group, is of the form

[TABLE]

where $\{p_{G}\}_{G}$ is a probability distribution, i.e. $p_{G}\geq 0,\forall G\in\mathsf{P}$ and $\sum_{G\in\mathsf{P}}p_{G}=1$. Let us now rewrite the optimization eq. 159 in terms of this probability distribution. We begin by noting that since $\mathcal{E}$ is TP we have that $\mathcal{E}_{\sigma_{0}\sigma_{0}}=1$ and hence we can write the depolarizing constraint in eq. 159 as

[TABLE]

Now, using the form of the Pauli-twirled channel, we can write the RHS of this equivalence as

[TABLE]

where in the last line we used that the identity Pauli element $I$ commutes with all Pauli matrices $\tau\in\boldsymbol{\sigma}_{q}\cup\sigma_{0}$, whereas every non-identity Pauli $G$ commutes with exactly of the elements of $\boldsymbol{\sigma}_{q}\cup\sigma_{0}$ and anti-commutes with the other half. We also used that $|\boldsymbol{\sigma}_{q}\cup\sigma_{0}|=d^{2}$. We can make a similar calculation for the objective of eq. 159 which gives

[TABLE]

where we have used that $\text{sgn}(\tau,G)\text{sgn}(\tau,\hat{G})=\text{sgn}(\tau,G\hat{G})$, that $GG^{\dagger}=I,\forall G\in\mathsf{P}$ and again that the Pauli identity $I$ commutes with all elements of $\boldsymbol{\sigma}_{q}\cup\sigma_{0}$ while every non-identity Pauli $G\hat{G}^{\dagger},G\neq G\hat{G}$ commutes with exactly half of the elements of $\boldsymbol{\sigma}_{q}\cup\sigma_{0}$ and anti-commutes with the other half. We have now rewritten the optimization eq. 159 completely in terms of the probability distribution $\{p_{G}\}_{G}$. This becomes

[TABLE]

Noting that the element $p_{I}$ is essentially fixed we can eliminate this element from the optimization and obtain an even simpler optimization

[TABLE]

The above optimization is a well studied instance of a class of optimization problems called quadratic programs Boyd and Vandenberghe (2004). This problem has the minimum (Boyd and Vandenberghe, 2004, Chapter 4, Section 4):

[TABLE]

and has $d^{2}-1$ degenerate maxima indexed by the non-identity Pauli elements $\tilde{G}$ of the form

[TABLE]

This means we can lower bound the quantity eq. 153, for any CPTP map $\mathcal{E}$, by:

[TABLE]

By now using the relation $f=1-\frac{dr}{d-1}$ we can rewrite this lower bound in terms of the infidelity $r$. This process is straightforward but rather tedious so we will not write it down. At the end of the calculation we obtain

[TABLE]

Similarly we can write for the upper bound

[TABLE]

which, by essentially the same tedious but straightforward calculation yields

[TABLE]

which completes the lemma.

∎

C.3 Eigenvalues of twirled quantum channels

Lemma 5 proves that the unitarity upper bounds the eigenvalues of the twirled superoperator $\mathcal{T}_{\phi_{\mathrm{TS}}}(\mathcal{E}^{\otimes 2})$. This resolves an open question posed in Wallman and Flammia (2014) and allows us to establish the long sequence length behavior of the variance of RB.

Lemma 5**.**

Let $\mathcal{E}:\mathcal{M}_{d}\to\mathcal{M}_{d}$ be a quantum channel with unitarity $u$ and consider the twirled operator $\mathcal{T}_{\phi_{\mathrm{TS}}}(\mathcal{E}^{\otimes 2})$ with respect to the traceless-symmetric representation. This operator can then be written as (lemmas 2 and 1)

[TABLE]

with $\mathcal{Z}=\{\mathrm{tr},1,2,[1],[2],[3],\{1\},\{2\}\}$, $\mathcal{P}_{i}$ the projector onto the spaces $V_{i}\subset\mathcal{M}_{d}^{\otimes 2}$ and

[TABLE]

where the trace is taken over superoperators. We now have for all $i\in\mathcal{Z}$ that

[TABLE]

Proof.

We begin by considering $i\in\mathcal{Z}_{\mathrm{d}}$. Note first that for $i=\mathrm{tr}$ we have that

[TABLE]

where we have used the definition of $\mathcal{P}_{\mathrm{tr}}$ (lemma 2). We can calculate

[TABLE]

where we have used the definition of the unitarity. Now consider $i\in\mathcal{Z}_{\mathrm{d}}$. We have

[TABLE]

Where we have used that the support of $\mathcal{P}_{i}$ lies in $V_{\mathrm{d}}$ (defined in lemma 2). Now we can use lemma 3 to upper bound this quantity. We have

[TABLE]

where we have again used the definition of the unitarity.

Next we consider the case of $i\in\mathcal{Z}_{[S]}$. We have

[TABLE]

where we have used that the support of $\mathcal{P}_{i}$ lies in $V_{[S]}$ (defined in lemma 2) and the factor of $1/4$ accounts for the fact that we are double counting the basis of $V_{[S]}$ since $S_{\sigma,\sigma\cdot\tau}=S_{\sigma\cdot\tau,\sigma}$ (we double count twice: once in the definition of the trace and once in the resolution of the identity on $V_{[S]}$). From lemma 3 we can lose one of the sums and get

[TABLE]

We can further use lemma 3 to upper bound this quantity as

[TABLE]

where we have also used the triangle inequality for the absolute value. Using the triangle inequality again together with the fact that $2|ab|\leq a^{2}+b^{2}$ for all $a,b\in\mathbb{R}$ we can write

[TABLE]

Now since $\sigma\in\bf{C}_{\tau}\iff\sigma\cdot\tau\in\bf{C}_{\tau}$ we can roll the four sums in the above expression into one, that is

[TABLE]

where we used the fact that $\sigma\in\bf{C}_{\tau}\iff\tau\in{\bf C}_{\sigma}$, the fact that $|{\bf{C}_{\sigma}}\cap{\bf{C}_{\sigma^{\prime}}}|\leq|{\bf{C}_{\sigma}}|=d^{2}/2-2$ and the definition of the unitarity. This means we have $\chi_{i}\leq u$ for all $i\in\mathcal{Z}_{[S]}$. The argument for $i\in\mathcal{Z}_{\{S\}}$ is conceptually the same as that for $i\in\mathcal{Z}_{[S]}$ so we will not write it down. ∎

Lemma 6 proves that the eigenvalues of the twirled superoperator $\mathcal{T}_{\phi_{\mathrm{TS}}}(\mathcal{E}^{\otimes 2})$ are close to the depolarizing parameter $f$. This fact is key in our analysis of the variance of RB in the presence of SPAM.

Lemma 6**.**

Let $\mathcal{E}:\mathcal{M}_{d}\to\mathcal{M}_{d}$ be a quantum channel with infidelity $r$ and depolarizing parameter $f=1-\frac{dr}{d-1}$ and consider the twirled operator $\mathcal{T}_{\phi_{\mathrm{TS}}}(\mathcal{E}^{\otimes 2})$ with respect to the traceless-symmetric representation. This operator can then be written as (lemmas 2 and 1)

[TABLE]

with $\mathcal{Z}=\{\mathrm{tr},1,2,[1],[2],[3],\{1\},\{2\}\}$, $\mathcal{P}_{i}$ the projector onto the spaces $V_{i}\subset\mathcal{M}_{d}^{\otimes 2}$ and

[TABLE]

where the trace is taken over superoperators. We now have for all $i\in\mathcal{Z}_{\mathrm{d}}$

[TABLE]

and for all $i\in\mathcal{Z}_{[S]}\cup\mathcal{Z}_{\{S\}}$

[TABLE]

subject to the constraint $r\leq\frac{1}{3}$

Proof.

From lemma 5 we have that $\chi_{i}\leq u$ for all $i\in\mathcal{Z}$. And since $u\leq 1$ for all quantum channels Wallman et al. (2015a) we certainly have that

[TABLE]

Hence we are only interested in upper bounding $f^{2}-\chi_{i}$, and thus lower bounding $\chi_{i}$ for all $i\in\mathcal{Z}$. First consider $i\in\mathcal{Z}_{\mathrm{d}}$. We proceed in much the same way as lemma 5. We have

[TABLE]

We begin by considering the first term in eq. 189. Using lemma 3 we can say

[TABLE]

where we have also used the lower bound from lemma 4. Now let us consider the second term in eq. 189. We have

[TABLE]

where we have again used lemma 3, the lower bound from lemma 4 and the definition of unitarity. We can now see that for $i\in\mathcal{Z}_{\mathrm{d}}$ we have

[TABLE]

Now consider $i\in\mathcal{Z}_{[S]}$ (note that we are implicitly taking $d\geq 4$ for this part of the proof, this is justified as the set $\mathcal{Z}_{[S]}$ is empty for $q=1$). From lemma 5 and in particular eq. 180 we get

[TABLE]

We can rewrite this a little bit as follows

[TABLE]

where we used that $\sigma^{\prime}\in{\bf{C}_{\tau}}\iff\sigma^{\prime}\cdot\tau\in{\bf{C}_{\tau}}$, that $(\sigma^{\prime}\cdot\tau)\cdot\tau=\sigma^{\prime}$ and that $S_{\sigma^{\prime},\sigma^{\prime}\cdot\tau}=S_{\sigma^{\prime}\cdot\tau,\sigma^{\prime}}$. We can again separate off the ‘diagonal’ terms to get

[TABLE]

We will analyze the terms eq. 198a and eq. 198b separately. We begin with eq. 198a. We can use lemma 3 to get

[TABLE]

Now we use the generic statement $2ab=a^{2}+b^{2}-(a-b)^{2}$ for all $a,b\in\mathbb{R}$ to write

[TABLE]

where we again used that $\sigma\in{\bf{C}_{\tau}}\iff\sigma\cdot\tau\in{\bf{C}_{\tau}}$ and that $\sigma\in{\bf{C}_{\tau}}\iff\tau\in{\bf{C}_{\sigma}}$ and also the lower bound from lemma 4. It remains to bound the second term in eq. 204. To do this we will maximize the quantity $(\mathcal{W}_{\nu,\nu}-\mathcal{W}_{\nu\cdot\mu,\nu\cdot\mu})^{2}$ for $\mu\in\bf{\sigma}_{q}$ and $\nu\in\bf{C}_{\mu}$ subject to the constraint that $\mathcal{W}$ is a CPTP map with depolarizing parameter $f$. That is, we will try to solve the maximization problem

[TABLE]

As in lemma 4 we can restrict ourselves to $\mathcal{W}$ being a Pauli channel (since the optimization function is a function of only the diagonal elements of $\mathcal{W}$). That is we can consider $\mathcal{W}(X)=\sum_{G\in\mathsf{P}}p_{G}GXG^{\dagger}$ where $\{p_{G}\}_{G}$ is a probability distribution over the Pauli group. We can write the optimization objective as

[TABLE]

where the $\text{sgn}(\nu,G)$ (as defined in eq. 157) encodes the commutation relations of the elements of the Pauli group. Note that the above quantity does not depended on $p_{\mathbb{1}}$ (the weight associated with the Pauli identity) since $\text{sgn}(\mu,\mathbb{1})=1$ for all $\mu\in\bf{\sigma}_{q}$. Hence we can solve the optimization problem

[TABLE]

This problem has an easily spotted maximum in that we want to put all probability weight on a single $G\in C_{\nu}\cap N_{\mu}$ and set all other $p_{G}$ to zero (subject to the constraint that the overall channel must have depolarizing parameter $f$, which is encoded in the first constraint of eq. 207 ). Hence we have

[TABLE]

We can feed this back into eq. 204 to obtain

[TABLE]

This is a suitable lower bound on eq. 198a. Next we consider eq. 198b. We have

[TABLE]

where we used an array of steps that have been used before: the triangle inequality and lemma 3 for the first inequality, the fact that $2|ab|\leq a^{2}+b^{2}$ for all $a,b\in\mathbb{R}$ for the second inequality, the fact that $\sigma\in\bf{C}_{\tau}\iff\tau\in\bf{C}_{\sigma}$ for the third equality, the fact that $|{\bf{C}_{\sigma}\cap\bf{C}_{\sigma}}|=d^{2}/4-3$ for $\sigma\neq\sigma^{\prime}$ Helsen et al. (2016) for the fourth equality and lemma 4 and the definition of unitarity for the last equality. This is a good lower bound on eq. 198b. We can now combine the lower bounds on eq. 198a and eq. 198b to get

[TABLE]

for $i\in\mathcal{Z}_{[S]}$. This gives a final bound (using $u\leq 1$)

[TABLE]

which we can rewrite to yield

[TABLE]

Setting $\left(1-\frac{1}{2}\frac{dr}{d-1}\right)\leq 1$ and working out we get

[TABLE]

for

[TABLE]

This completes the proof for $i\in\mathcal{Z}_{[S]}$. The proof for $i\in\mathcal{Z}_{\{S\}}$ is conceptually the same as that of $i\in\mathcal{Z}_{[S]}$ and yields the same bound so we will not write it down here. The only notable difference is the difference in size for the sets ${\bf{N}}_{\tau}$ and ${\bf{N}}_{\tau}\cap{\bf{N}}_{\tau^{\prime}}$ for $\tau,\tau^{\prime}\in\bf{\sigma}_{q}$ which gives a different area of validity for the bound, namely

[TABLE]

Choosing $r\leq 1/3$ satisfies both constraints for all $d$ and thus completes the proof. ∎

C.4 Telescoping series

Lemmas 7 and 1 provide us with a powerful tool to break up the analysis of the variance of RB into manageable pieces.

Lemma 7**.**

For two arbitrary ordered lists of $m$ elements $\{a_{1},\ldots,a_{m}\}$ and $\{b_{1},\ldots,b_{m}\}$ of an algebra with associative and distributed addition and multiplication we have,

[TABLE]

where $a_{j:k}$ with $j\geq k$ is defined with respect to the list $\{a_{1},\ldots,a_{m}\}$ as

[TABLE]

Proof.

We will prove this by induction. For $m=1$ the statement is trivial. For $m+1$, we have

[TABLE]

by induction hypothesis. This proves the lemma. ∎

Corollary 1**.**

For $a,b,c\in\mathbb{C}$ with $c\geq a$, we have

[TABLE]

Proof.

Note first that the statement is trivial if $a=b$. Therefore assume $a\neq b$. We begin by applying lemma 7 to $a^{m}-b^{m}$. This gives

[TABLE]

We now perform the following manipulation

[TABLE]

Note that be have used the fact that $a,b\in\mathbb{C}$ are commutative. Now we can apply lemma 7 again to the factors $(a^{m-j}-b^{m-j})$ in the second term in the above to obtain

[TABLE]

Performing the substitution $s=j+t$ and working out we obtain

[TABLE]

Now we can factor out $a^{m-1}$ from the second term to obtain

[TABLE]

We can further rewrite this using the standard series identity

[TABLE]

The upper bound follows by upper bounding each term in the sum. ∎