Statistical analysis of randomized benchmarking

TL;DR

This paper improves randomized benchmarking for quantum computers by proposing a modification that enhances accuracy, simplifies experimental procedures, and provides reliable error estimates with credible regions.

Contribution

It introduces a simple modification to RB+ that eliminates a nuisance parameter and improves the efficiency and reliability of error rate estimation.

Findings

Eliminates a nuisance parameter in RB+

Provides a method for credible regions of parameters

Achieves error estimates with multiplicative precision

Abstract

Randomized benchmarking and variants thereof, which we collectively call RB+, are widely used to characterize the performance of quantum computers because they are simple, scalable, and robust to state-preparation and measurement errors. However, experimental implementations of RB+ allocate resources suboptimally and make ad-hoc assumptions that undermine the reliability of the data analysis. In this paper, we propose a simple modification of RB+ which rigorously eliminates a nuisance parameter and simplifies the experimental design. We then show that, with this modification and specific experimental choices, RB+ efficiently provides estimates of error rates with multiplicative precision. Finally, we provide a simplified rigorous method for obtaining credible regions for parameters of interest and a heuristic approximation for these intervals that performs well in currently relevant regimes.

Equations

$$\bar{q}(m)=A{p}^m+B$$

$$\begin{array}{rl}\hat{a}(m|P)& =\sum _s\hat{q}(m,s|P)/k_m\\ \hat{a}(m)& =\frac{1}{4^n-1}\sum _P\hat{a}(m|P)\\ \hat{b}(m)& =\frac{1}{k_m}\sum _{P,s}\hat{q}(m,s|P{)}^2-\hat{a}(m|P{)}^2.\end{array}$$

$$\begin{array}{rl}\bar{a}(m)& =A{l}^m\\ \bar{b}(m)& =A^{\prime }{u}^m\end{array}$$

$$\displaystyle B:=\operatorname{Tr}\bigl{[}E\Lambda(\mathbbm{1}/2^{n})\bigr{]}$$

$$\hat{q}(m|b)=\frac{1}{|{\mathbb{S}}_{m,b}|}\sum _{s\in {\mathbb{S}}_{m,b}}q(m,s).$$

$$\bar{q}(m)=\bar{q}(m|0)-\bar{q}(m|1)=A{p}^m$$

$$\displaystyle B=\frac{1}{k}\sum_{j=1}^{k}\operatorname{Tr}\bigl{[}E_{j}\Lambda(\tfrac{\mathbbm{1}}{2^{n}})\bigr{]}=\frac{1}{k}\operatorname{Tr}\bigl{[}\mathbbm{1}\Lambda(\tfrac{\mathbbm{1}}{2^{n}})\bigr{]}=\frac{1}{k}\,.$$

$$\begin{array}{rl}A& ={\left[\bar{q}(m_1)-B\right]}^{m_2/\delta m}{\left[\bar{q}(m_2)-B\right]}^{-m_1/\delta m},\\ p& ={\left[\bar{q}(m_1)-B\right]}^{-1/\delta m}{\left[\bar{q}(m_2)-B\right]}^{1/\delta m}.\end{array}$$

$$\mathbb{E}{\left[\hat{q}(m_j)-B\right]}^{\alpha }$$

$$=x_j^{\alpha }\left[1+\frac{1}{2}\alpha (\alpha -1)\mathbb{E}{ϵ}_j^2+O(\mathbb{E}\alpha {ϵ}^3)\right]$$

$$\approx x_j^{\alpha }+\frac{1}{2}\alpha (\alpha -1)x_j^{\alpha -2}\mathbb{V}\hat{q}(m_j).$$

$$\mathbb{V}{\left[\hat{q}(m_j)-B\right]}^{\alpha }$$

$$\displaystyle\Pr\Bigl{(}\lvert\hat{p}-p\rvert>k\sqrt{\mathbb{V}(\hat{p})}\Bigr{)}\leq k^{-2},$$

$$\displaystyle\lvert\hat{p}-p\rvert\leq\frac{3}{\delta\!m}\biggl{(}\sum_{j}\mathbb{V}[\hat{q}(m_{j})]\biggr{)}^{1/2}+O\bigl{(}\mathbb{E}(\epsilon_{j}^{3/2})/\delta\!m^{3/2}\bigr{)}.$$

$$\begin{array}{rl}\log \frac{\hat{q}(m_2)}{\hat{q}(m_1)}& \sim \mathcal{N}\left(\log (\frac{\bar{q}(m_2)}{\bar{q}(m_1)}),{\sigma }^2\right)\\ {\sigma }^2& =\sum _j\frac{\overline{q}(m_j)\left(1-\overline{q}(m_j)\right)}{k_j(\overline{q}(m_j)-B{)}^2}.\end{array}$$

$$\log \hat{p}=\frac{1}{\delta m}\log \frac{\hat{q}(m_2)}{\hat{q}(m_1)}\sim \mathcal{N}(\log p,{\sigma }^2/\delta m^2).$$

$$\underset{m_2}{\mathrm{argmin}}\left(\log \left[p^{-2m_1}\overline{q}(m_1)\left(1-\overline{q}(m_1)\right)+p^{-2m_2}\overline{q}(m_2)\left(1-\overline{q}(m_2)\right)\right]-2\log (m_2-m_1)\right).$$

$$\displaystyle M=O\biggl{(}\frac{1}{\epsilon^{2}}\log\bigl{(}\tfrac{1}{r}\bigr{)}\log\Bigl{[}\tfrac{1}{\delta}\log\bigl{(}\tfrac{1}{r}\bigr{)}\Bigr{]}\biggr{)}$$

$$\frac{(1-4ϵ{)}^2}{9}<p^m\le \frac{1+4ϵ}{3}.$$

$$\displaystyle\hat{q}_{\ell}=Ap(p^{2^{\ell}}+\epsilon_{\ell})\leq\tfrac{1}{3}\hat{q}_{1}=\frac{Ap}{3}\bigl{(}1+\epsilon_{1}\bigr{)}.$$

$$\displaystyle\hat{q}_{\ell-1}=Ap(p^{2^{\ell-1}}+\epsilon_{\ell-1})>\tfrac{1}{3}\hat{q}_{1}=\frac{Ap}{3}\bigl{(}1+\epsilon_{1}\bigr{)}.$$

$$\frac{(1-3{ϵ}_{\ell -1}+{ϵ}_1{)}^2}{9}<p^m\le \frac{1-3{ϵ}_{\ell }+{ϵ}_1}{3}.$$

$${\hat{p}}_{\pm }:={\left(\frac{q_{\ell }\pm Apϵ}{q_1\mp Apϵ}\right)}^{1/m}=p{\left(\frac{1\pm ϵ/p^m}{1\mp ϵ}\right)}^{1/m}.$$

$${\hat{p}}_{+}\le p{\left(1+\frac{2ϵ\left(5-4ϵ+8{ϵ}^2\right)}{(1-4ϵ{)}^2(1-ϵ)}\right)}^{1/m}\text{and }{\hat{p}}_{-}>p{\left(1-\frac{2ϵ\left(5-4ϵ+8{ϵ}^2\right)}{(1-4ϵ{)}^2(1+ϵ)}\right)}^{1/m}.$$

$$\displaystyle\hat{p}_{+}<p\bigl{(}1+O(\epsilon)\bigr{)}^{1/m}\quad\text{and }\quad\hat{p}_{-}>p\bigl{(}1-O(\epsilon)\bigr{)}^{1/m},$$

$$\displaystyle\hat{p}_{+}<p\bigl{(}1+O(\epsilon q)\bigr{)}\quad\text{and}\quad\hat{p}_{-}>p\bigl{(}1-O(\epsilon r)\bigr{)}.$$

$$|\hat{p}-p|=|\hat{r}-r|\le O(ϵr).$$

Full text

I Introduction

Characterizing large scale quantum devices is a prerequisite to optimizing their performance and being able to reliably perform useful information processing tasks. Full characterization is manifestly not scalable for general errors, so that scalable methods can only partially characterize the noise. Currently, the only fully scalable protocols that partially characterize quantum devices are randomized benchmarking Knill et al. (2008); Magesan et al. (2011, 2012a) and variants thereof (RB+) Magesan et al. (2012b); Carignan-Dugas et al. (2015); Wallman et al. (2015a, b); Sheldon et al. (2016); Cross et al. (2016); Harper and Flammia (2017). This family of protocols can provide a wide variety of information about noise parameters, including the average error rate Knill et al. (2008); Magesan et al. (2011, 2012a); Wallman (2018), error rates for specific gates Magesan et al. (2012b); Carignan-Dugas et al. (2015); Cross et al. (2016); Harper and Flammia (2017), leakage rates Wallman et al. (2016), loss rates Wallman et al. (2015b); Wallman (2018), and the amount of residual unitary (calibration) errors Wallman et al. (2015a); Sheldon et al. (2016); Yang et al. (2018).

RB+ provides estimates of noise parameters by applying long sequences of random gates to amplify errors in the implementation of gates and estimate them independently from state preparation and measurement errors (SPAM). Typically, descriptions of RB+ state that experiments should be repeated to obtain a desired precision without necessarily specifying (or recommending) any of the following: (1) estimators for finite data; (2) how many repetitions should be performed; or (3) how finite and heteroscedastic data should be fit to a specified model. This last point is important because RB+ data are generally heteroscedastic, meaning that the variance across the data is non-uniform, since the variance over random sequences increases with the sequence length Wallman and Flammia (2014). Abstaining from specifics on these points was perhaps warranted by the fact that particular choices are difficult to derive or justify as being optimal, or nearly optimal. Obtaining a fully general and optimal specification is confounded by the unknown distribution of errors over the random sequence of gates Wallman and Flammia (2014); Helsen et al. (2017). However, in one experimental regime, Bayesian techniques can be applied to obtain rigorous credible intervals for the model parameters as well as efficient allocation of experimental measurements Granade et al. (2014, 2017). We discuss and utilize this work by Granade et al. in section IV.

In this paper, we present a minor modification of RB+ that improves the efficiency of the method by eliminating a nuisance model parameter, yet it adds no experimental overhead. Similar methods have been presented previously in the literature for the case of single-qubit RB Muhonen et al. (2015); Fogarty et al. (2015). For the two remaining model parameters, we then provide estimators which do not have to be weighted to correct for heteroscedasticity because they can be estimated from two independent sequence lengths. Without our modification, at least 3 sequence lengths are required, which in turn require a weighted fit where the correct weights are not generally inferable from the data. We then study the distribution of the parameter estimators and show how to obtain simple and rigorous credible intervals in the regime studied in Ref. Granade et al. (2014), that is, when each random sequence is repeated once. Finally, we also provide a simple proof that certain experimental design choices enable RB+ to efficiently provide estimates of error rates that have multiplicative precision. By showing that the estimates of such error rates have multiplicative precision, we confirm that RB+ will continue to allow efficient estimation of the model parameters as gate fidelity rates improve through the simple expedient of increased sequence lengths.

In what follows, we use the notation that $\hat{x}$ is an estimator of a quantity $\bar{x}$, where the bar denotes that either an expected value or a sample average has been taken over realizations of a random variable $x$.

II RB+ protocol

We begin by providing a general framework that describes all existing RB+ protocols except the leakage protocol of Ref. Wood and Gambetta (2018) and the unitarity protocol of Ref. Wallman et al. (2015a). We will also present a modified version of the unitarity protocol. We exclude the leakage protocol of Ref. Wood and Gambetta (2018) because it has a more complicated fit model that is not robust to SPAM errors.

RB+ protocols are of the following form.

Choose a positive integer $m$.

Choose a random sequence of gates $s$ from a set $\mathbb{S}_{m}$, typically of Clifford gates. Note that these gates are often chosen to leave a state invariant. However, as discussed below, uniformly choosing $s$ to either leave the state invariant or map it to an orthogonal state eliminates a nuisance model parameter.

Obtain an estimate $\hat{q}(m,s)$ of the expectation value $q(m,s)$ of an observable $E$ after preparing a state $\rho$ and applying the gates in $s$. Typically $\rho$ should be close to an ideal computational basis state and $E$ should be close to a projector onto a pure state in the computational basis.

Repeat steps 2–3 $k_{m}$ times to obtain an estimate $\hat{q}(m)$ of $\bar{q}(m)=\lvert\mathbb{S}_{m}\rvert^{-1}\sum_{s\in\mathbb{S}_{m}}q(m,s)$.

Repeat steps 1–4 and fit to the model

$\displaystyle\bar{q}(m)=Ap^{m}+B$

(1)

where $p$ is related to some parameter of interest (e.g., the average gate fidelity to the identity) and $A$ and $B$ are SPAM-dependent constants.

We assume throughout that $A\gg 0$, which holds in current regions of interest, as otherwise it is unclear how to efficiently gather useful statistics.

II.1 Unitarity

We now introduce a slight variant of the RB+ protocol to handle the special case of the unitarity protocol from Ref. Wallman et al. (2015a). The variant enables an independent estimate of the unitarity (which quantifies how coherent the errors are) and the leakage rate Wallman et al. (2015b, 2016). Note that the following protocol is not strictly scalable as it involves sampling every Pauli matrix and also assumes that there is no (or minimal) state-dependent loss. The scalability could be improved by, for example, performing importance sampling of the Pauli matrices conditioned on the sequence Flammia and Liu (2011); however, we leave this as an open problem.

Choose a positive integer $m$. 2. 2.

Choose a random sequence of $m$ $n$-qubit Clifford gates $s$. 3. 3.

For each $n$-qubit Pauli matrix $P$, obtain an estimate $\hat{q}(m,s|P)$ of the expectation value of the observable $P$ after preparing a fixed state $\rho$ and applying the gates in $s$. 4. 4.

Repeat steps 2–3 $k_{m}$ times. For each $n$-qubit Pauli matrix $P\neq I$, set

[TABLE] 5. 5.

Repeat steps 1–4 and fit to the models

[TABLE]

where $l$ and $u$ are the leakage rate Wallman et al. (2016) and the unitarity Wallman et al. (2015a) respectively.

Unlike other protocols, the combined unitarity/loss protocol does not require any truncation of $\bar{q}(m)$ to avoid negative values. We also note that recently an alternative protocol has been proposed which proposes a method for efficient unitarity benchmarking in the regime of few-qubit Clifford gates Dirkse et al. (2018).

II.2 Eliminating the offset

While superficially benign, the variable offset $B$ in eq. 1 can severely increase the marginal uncertainty in $p$, the parameter of interest Muhonen et al. (2015). We now present a method of rigorously and exactly eliminating this constant offset without having to estimate its value.

First note that for gate independent noise $\Lambda$, the constant

[TABLE]

and the decay parameter $p$ Magesan et al. (2012a); Wallman (2018) do not change if we compile any gate into the sequence. (In fact, strictly this holds even for non-trace-preserving noise where $B$ is multiplied by a second exponential). In particular, let $X$ be any gate that maps the input state to an orthogonal state, as the single-qubit Pauli $X$ operator does for states in the computational basis. Let $\mathbb{S}_{m,b}$ be the set of sequences obtained from $\mathbb{S}_{m}$ by compiling $X^{b}$ into the sequence and let

[TABLE]

Then we have

[TABLE]

where now $A\in[0,1]$. A similar idea was suggested for single qubits in Muhonen et al. (2015); Fogarty et al. (2015). One disadvantage of this approach is that the remaining $A$ coefficient in eq. 1 may be small for some values of $b$, especially for multiple qubits, thus reducing the signal from some experiments.

Alternatively, consider an $n$-qubit POVM $\{E_{1},\ldots,E_{k}\}$ and suppose that a set of gates $\{X_{1},\ldots,X_{k}\}$ are such that $E_{j}\approx X_{j}E_{1}X_{j}^{\dagger}$. Then by compiling $X_{j}$ into the sequence uniformly at random and recording the probability of observing the corresponding $E_{j}$ and averaging over $j$, the average value of $B$ becomes

[TABLE]

III Estimating the decay rate

With a known value of $B$, eqs. 1 and 3 have two unknown parameters and so we need at least two values of $m$, denoted $m_{1}<m_{2}$, to estimate either (or both) parameters. Alternatively, we could use a single value of $m$ if we are content to accept a lower bound on $p$ (by assuming that $A=1$, and that other model assumptions are respected), and this may be sufficient for certain purposes.

From eq. 1,

[TABLE]

Each of these terms is of the form $x_{1}^{\alpha_{1}}x_{2}^{\alpha_{2}}$ where $x_{j}^{\alpha}=\left[\bar{q}(m_{j})-B\right]^{\alpha}$. A natural approach would be to estimate $x_{j}^{\alpha}$ by $\left[\hat{q}(m_{j})-B\right]^{\alpha}$, where $\hat{q}(m_{j})$ is an unbiased estimator for $\bar{q}(m_{j})$ (e.g., the sample mean). However, this approach has two issues. First, the estimator is complex or undefined if $\hat{q}(m_{j})\leq B$. To address this issue, we can truncate $\bar{q}(m_{j})$ to $B+\delta$ for some fixed $0<\delta\ll 1$, which will occur with negligible probability provided enough sequences are taken and the $\hat{q}(m_{i})$ are sufficiently far from [math].

A second issue is that if $\hat{q}(m_{j})$ is an unbiased estimator for $\bar{q}(m_{j})$, then $\left[\hat{q}(m_{j})-B\right]^{\alpha}$ is a biased estimator of $x_{j}^{\alpha}$ for any $\alpha\neq 1$. To address this second issue, we can estimate and then substract the bias if necessary. Repeating steps 2–4 for a fixed $m$ and using final gates compiled in to remove the offset $B$ yields an estimate $\hat{q}(m_{j})-B=x_{j}(1+\epsilon_{j})$ for some random variable $\epsilon_{j}$ with zero mean. We then have

[TABLE]

Therefore $\left[\hat{q}(m_{j})-B\right]^{\alpha}$ is biased but consistent and so the bias can be neglected when sufficiently many sequences are sampled at each $m_{j}$. Moreover, the bias is modulated by $\alpha$, which as we prove below in section V, is $O(r)$ in the optimal regime for $p$. For small numbers of sampled sequences, the bias term can be subtracted using sample estimates of $x_{j}^{\alpha}$ and $\mathbb{V}\hat{q}(m_{j})$ on the right-hand-side of section III. As our numerics will show in section IV, the bias is negligible for intermediate numbers of measurements but is noticeable for very small numbers of sequences.

To determine approximately optimal values of $m_{1}$ and $m_{2}$ assuming that the bias in section III is negligible, note that

[TABLE]

Now we note that by Chebyshev’s inequality,

[TABLE]

so that with probability 8/9 (for example) and using the trivial bound $p\leq 1$, we have

[TABLE]

The only unknowns in eq. 12 are the variances at the two sequence lengths. Choosing $\delta\!m\approx 1/(1-p)$ therefore gives a multiplicative precision estimate of the error rate $1-p$, as claimed.

We would like to make the error term as small as possible. To achieve a large $\delta\!m$, and hence a small error, we want $m_{1}$ as small as possible and $m_{2}$ as large as possible. However, eq. 1 is only typically accurate for $m\geq 4$ Wallman (2018), so we henceforth set $m_{1}=4$. Furthermore, the number of sequences required to make the truncation probability negligible increases with $m_{2}$, and eq. 10 is inversely proportional to $x_{j}$ for $\alpha\leq 1$, so $m_{2}$ can only be increased to some fixed value. Thus, when sampling to some fixed constant accuracy we must not choose $m_{2}$ to be too large.

Ref. Granade et al. (2014) recommended $m_{2}=\lceil 1/(1-p)\rceil$ as the optimal choice of $m_{2}$, however, this was for $B=0$ (i.e., the infinite-dimensional limit). Numerically, we observe that $m_{2}=\lceil 1/[2(1-p)]\rceil$ results in a more precise estimate.

Consequently, provided gate lengths ($m$) can be increased as specified above, the number of sequences required to determine $1-p$ to within a specified factor remain approximately independent of $p$. This demonstrates that RB+ protocols scale favorably with the error rate.

This derivation is not quite rigorous only for a very trivial reason, namely because our use of Taylor’s theorem requires that we control the smoothness of the functions being expanded over some region, and this region should also be appropriately defined. These expressions are nonetheless useful for practical analysis of RB+. By contrast, the derivation in section V is completely rigorous, but the proof is not meant to provide anything more than coarse guidance about how to choose parameter settings in practical situations.

IV Accelerated RB

We now analyze accelerated randomized benchmarking (ARB) Granade et al. (2014) using the modified protocol discussed in section II.2. In ARB, each individual estimate $\hat{q}(m,s)$ is in $\{0,1\}$, that is, each random sequence is measured once. Therefore, $\hat{q}(m)\sim\mathcal{B}(k_{m},\bar{q}(m))/k_{m}$, where $\mathcal{B}(n,p)$ denotes the binomial distribution with $n$ trials and probability $p$.

For sufficiently many samples, log ratios of binomial variables are approximately normally distributed Katz et al. (1978), so that

[TABLE]

Therefore

[TABLE]

In fig. 1, we illustrate that the normal approximation is sufficiently accurate by comparing the exact cumulative density function for the estimator $\hat{p}$ from binomial statistics with the normal approximation of eq. 14 for multiple values of $A$, $B$, and $p$. Note in particular that the shape of the cumulative density function for the estimator $\hat{p}$ is essentially independent of $p$, but has heavier tails for smaller values of $A$.

Under the log-normal approximation for $\hat{p}$, the value of $m_{2}$ that minimizes the variance of the estimate is given by

[TABLE]

While this minimizes the variance of $\log\hat{p}$ rather than $\hat{p}$, since $p$ near $1$, we have $\mathbb{V}(\log\hat{p})\approx\mathbb{V}(\hat{p})$. This minimization can be performed numerically using an initial value of $-1/\log p$. In fig. 2, we plot the variance and optimal $m_{2}$ values in several relevant parameter regimes. The optimal value of $m_{2}$ depends on the true values of $A$, $B$, and $p$. In fig. 2(a-b) we see that when choosing a future experiment based on present knowledge with multiplicative uncertainty in $p$, it is best to err on the side of $m_{2}$ that is short with respect to the optimal value.

V Rigorous proof of multiplicative precision for the ratio estimator

In this section we give a rigorous proof that RB converges to an estimate with multiplicative precision using the ratio estimator in eq. 8. Here the focus is not on obtaining tight answers, but on having a simple and clear statement of the scaling of the precision that can be achieved assuming the decay model in eq. 6, and achieving a given sample complexity. We therefore largely neglect to track estimation errors closely, focusing instead on the simplest proof possible and a big-$O$ estimate of the resources required.

Let $r=1-p$ and fix some small $1/16>\epsilon_{0}>0$. We are most interested in the regime where $r$ is small, or equivalently $p$ is close to 1. We assume that we can estimate the quantities $q_{i}=Ap^{m_{i}}$ with an unbiased estimator $\hat{q}_{i}=q_{i}+Ap\epsilon_{i}$ where the estimation error $\epsilon_{i}$ is a random variable with zero mean. Note that this is multiplicative precision for the case $m_{i}=1$, but for larger values of $m_{i}$ we have just rescaled an additive precision by $Ap$ for algebraic convenience. If the estimator is the sample mean of $t$ Bernoulli random variables with mean $q_{i}$, then $\epsilon_{i}=O(1/\sqrt{t})$ with high probability.

Under these conditions, we can estimate $\hat{p}$ using the following algorithm.

Algorithm 1

Ratio estimator for exponential regression.

1. 1.

Set $i:=1$ and $m_{1}:=1$. 2. 2.

Estimate $\hat{q}_{1}:=q_{1}(1+\epsilon_{1})$ using $t$ samples. 3. 3.

While $\hat{q}_{i}>\frac{1}{3}\hat{q}_{1}$, Do

• •

Set $i:=i+1$,

• •

Set $m_{i}:=2^{i}+1$,

• •

Estimate $\hat{q}_{i}:=q_{i}+Ap\epsilon_{i}$ using $t$ samples. 4. 4.

Set $\ell:=i$ and $m=2^{\ell}$. 5. 5.

Return \hat{p}:=\Bigl{(}\frac{\hat{q}_{\ell}}{\hat{q}_{1}}\Bigr{)}^{1/m} and $\hat{r}:=1-\hat{p}$.

We now rigorously prove that the above algorithm returns an estimator with multiplicative precision of $\hat{r}$.

Theorem 1**.**

For any sufficiently small $\epsilon_{0}>0$, the algorithm above returns estimates $\hat{r}$ such that $|\hat{r}-r|\leq O(\epsilon r)$ with probability $1-\delta$ using

[TABLE]

measurements.

The proof relies on a few simple lemmas, which we now state and prove.

Lemma 2**.**

Given a set of $\ell$ independent estimates $\hat{q}_{i}$ obtained from sampling each $t$ times as described above, the probability that $|\hat{q}_{i}-q_{i}|\geq Ap\epsilon$ for any $i>1$ or $|\hat{q}_{1}-q_{1}|\geq q_{1}\epsilon$ is at most $\delta$ if we choose t=O\bigl{(}\frac{1}{\epsilon^{2}}\log\frac{\ell}{\delta}\bigr{)}, where the implied constant depends on $Ap$.

Proof.

The proof is an elementary application of the Chernoff bound and the union bound. We omit the details. ∎

Thus, we can assume that each of the random estimates $\epsilon_{i}$ satisfies $|\epsilon_{i}|\leq\epsilon$ in the algorithm, and we will fail with probability at most $\delta$. Next, we will see that the algorithm converges in a small number of steps $\ell$, and with $m$ taking a value that scales like $1/r$.

Lemma 3**.**

With probability at least $1-\delta$, the above algorithm converges with \ell=\Theta\bigl{(}\log\tfrac{1}{r}\bigr{)} using O\bigl{(}\frac{\ell}{\epsilon^{2}}\log\frac{\ell}{\delta}\bigr{)} total samples, and with $m$ such that

[TABLE]

Proof.

The algorithm exits the while loop when

[TABLE]

As the algorithm did not exit for $i=\ell-1$,

[TABLE]

Squaring the latter inequality and then rearranging both to be in terms of $p^{m}$ with $m=2^{\ell}$, we have

[TABLE]

Supposing that for some fixed $\epsilon>0$, $\lvert\epsilon_{i}\rvert\leq\epsilon$ for all $i$ with with probability $1-\delta$, the claim about $p^{m}$ follows by taking the worst-case choices of the $\epsilon_{i}$. Taking logarithms of this and using $p=1-r$, we find that for any sufficiently small $\epsilon$ we have \ell=\Theta\bigl{[}-\log\bigl{(}-\log(1-r)\bigr{)}\bigr{]}. As long as $r$ is bounded away from $1$ then $-\log(1-r)=\Theta(r)$, and this is equivalent to \ell=\Theta\bigl{(}\log\tfrac{1}{r}\bigr{)}.

If we sample as per 2, then $|\epsilon_{i}|\leq\epsilon$ for all $i$ with probability $1-\delta$. Therefore the claim about the total number of samples follows immediately from 2. ∎

Now we are ready to prove the main theorem.

Proof of Theorem.

From the above lemmas, we know that with probability $1-\delta$ the ratio estimator converges with m=O\bigl{(}\log\tfrac{1}{r}\bigr{)} and errors bounded by $\epsilon$ in the numerator and denominator. We have the bounds $\hat{p}_{-}\leq\hat{p}\leq\hat{p}_{+}$, where

[TABLE]

From the inequality in eq. 17, we have

[TABLE]

Now we choose any $\epsilon_{0}<1/16$ so that the $\epsilon$ dependent terms above are $O(\epsilon)$ and the term for $\hat{p}_{-}$ remains less than 1. Explicitly evaluating the $\epsilon$ dependent terms, we have the bounds

[TABLE]

where the implied constant decreases with $\epsilon_{0}$. Now Taylor expanding in $1/m$ and using the result from 3 that \ell=\log_{2}m=\Theta\bigl{(}\log\tfrac{1}{r}\bigr{)}, we find that

[TABLE]

Adopting the bounds $\hat{r}_{\pm}=1-\hat{p}_{\mp}$ gives the analogous result for $\hat{r}$. This establishes that the estimator has multiplicative precision,

[TABLE]

The statement about complexity follows directly from the lemmas, and the theorem is proven. ∎

VI Conclusion

We have provided a modification to RB+ and concrete recommendations for how to obtain precise estimates of error rates in practical regimes. We have rigorously shown that the precision is multiplicative, and our derivations and numerics demonstrate the utility of the heuristics that we use. Our recommendations are based upon the assumption that the model in eq. 1 is correct. For standard randomized benchmarking, there are only two factors that can cause a deviation from eq. 1 for sequence lengths $m\geq 4$, namely, noise that is time-dependent or non-Markovian Epstein et al. (2014); Wallman and Flammia (2014); Ball et al. (2016); Wallman (2018). Both types of noise are ubiquitous in experiments and neither can be detected using only two sequence lengths. A standard approach is to take data from more sequence lengths, perform a joint fit and then use the goodness-of-fit as an indicator for non-Markovian noise or drift. However, fitting more sequence lengths is nontrivial as the data are heteroscedastic. Furthermore, adding more sequence lengths does not significantly increase the quality of the error estimates when eq. 1 is correct, and so performing a joint fit provides little extra information and introduces correlations between model estimation and model validation. We instead recommend fitting data using only two sequence lengths and then using hypothesis testing to determine if data taken at other sequence lengths are consistent with the hypothesis that the noise is static and Markovian.

Acknowledgements.

RH and STF were supported by the Australian Research Council through the Centre of Excellence in Engineered Quantum Systems CE170100009. This research was supported by the US Army Research Office through grant numbers W911NF-14-1-0098 and W911NF-14-1-0103. IH and JJW gratefully acknowledge contributions from the Canada First Research Excellence Fund, Industry Canada, the Province of Ontario, and Quantum Benchmark Inc.