Error Mitigation: ZNE & Readout Correction

1. 1
   Step 1:

   $$\ell_i(x)=\prod_{j\ne i}\frac{x-\lambda_j}{\lambda_i-\lambda_j},\qquad\ell_i(\lambda_k)=\delta_{ik}$$

   The Lagrange basis polynomials of the nodes: each has degree $m-1$, vanishes at every node except $\lambda_i$, and equals $1$ there. They exist because the nodes are distinct.
2. 2
   Step 2:

   $$\deg q\le m-1\ \Rightarrow\ q(x)=\sum_{i=1}^mq(\lambda_i)\,\ell_i(x)\quad\text{for all }x$$

   The difference of the two sides is a polynomial of degree at most $m-1$ with the $m$ roots $\lambda_1,\dots,\lambda_m$, so it is the zero polynomial.
3. 3
   Step 3:

   $$\ell_i(0)=\prod_{j\ne i}\frac{0-\lambda_j}{\lambda_i-\lambda_j}=\prod_{j\ne i}\frac{\lambda_j}{\lambda_j-\lambda_i}=\gamma_i$$

   Multiply numerator and denominator of every factor by $-1$. The Richardson weights are the Lagrange basis polynomials evaluated at $x=0$.
4. 4
   Step 4:

   $$q(x)=x^k,\ 0\le k\le m-1:\qquad\sum_{i=1}^m\gamma_i\lambda_i^k=0^k=\delta_{k0}$$

   Evaluate the interpolation identity at $x=0$. For $k=0$ this is $\sum_i\gamma_i=1$, statement (i); for $1\le k\le m-1$ it is statement (ii). In matrix form the weights solve the Vandermonde system $V^\top\gamma=e_0$ with $V_{ik}=\lambda_i^k$, which is nonsingular because $\det V=\prod_{i<j}(\lambda_j-\lambda_i)\ne0$.
5. 5
   Step 5:

   $$x^m-\sum_{i=1}^m\lambda_i^m\,\ell_i(x)=\prod_{i=1}^m(x-\lambda_i)\ \overset{x=0}{\Longrightarrow}\ -\sum_i\gamma_i\lambda_i^m=(-1)^m\prod_i\lambda_i$$

   The left side is monic of degree $m$ (the sum has degree at most $m-1$) and vanishes at all $m$ nodes, so it equals the monic polynomial with these roots. Multiplying by $-1$ gives statement (iii).
6. 6
   Step 6:

   $$\sum_i\gamma_iE(\lambda_i)=E_0\sum_i\gamma_i+\sum_{k=1}^{m-1}a_k\sum_i\gamma_i\lambda_i^k+\sum_i\gamma_iR(\lambda_i)=E_0+\sum_i\gamma_iR(\lambda_i)$$

   Insert the expansion and use (i) and (ii): the constant survives with weight one and every term of order $1,\dots,m-1$ cancels exactly.
7. 7
   Step 7:

   $$\Big\lvert\sum_i\gamma_iR(\lambda_i)\Big\rvert\le C\sum_i\lvert\gamma_i\rvert\lambda_i^m,\qquad R=a_m\lambda^m+O(\lambda^{m+1})\ \Rightarrow\ \sum_i\gamma_iR(\lambda_i)=(-1)^{m-1}a_m\prod_i\lambda_i+O(a_{m+1})$$

   Triangle inequality for the bound; statement (iii) for the leading term. With $a_k=b_k\epsilon^k$ every remaining term carries at least the factor $\epsilon^m$.

∎

1. 1
   Step 1:

   $$\mathcal N_{\text{seq}}=\mathcal U_{\text{seq}}+\epsilon\sum_t\mathcal A_t\,\mathcal L_t\,\mathcal B_t+O(\epsilon^2)$$

   Multiply out any sequence of ideal gates and noisy locations $t$: a first-order term picks $\epsilon\mathcal L_t$ at exactly one location and $\mathcal I$ at all others; $\mathcal B_t$ is the product of the ideal gates before location $t$ and $\mathcal A_t$ the product after it.
2. 2
   Step 2:

   $$U(U^\dagger U)^k=U\ \Rightarrow\ \mathcal U_{\text{fold}}=\mathcal U\ \Rightarrow\ E(\lambda)=E_0+O(\epsilon)$$

   The ideal gates of the folded sequence multiply to $U$ because $U^\dagger U=I$, so the zeroth-order term is the noiseless value for every $k$.
3. 3
   Step 3:

   $$\text{forward block }r\ (0\le r\le k),\ \text{after }G_j:\quad\mathcal B_t=\mathcal U_{\le j}\big(\mathcal U^\dagger\mathcal U\big)^r=\mathcal U_{\le j},\qquad\mathcal A_t=\big(\mathcal U\,\mathcal U^\dagger\big)^{k-r}\,\mathcal U_{>j}=\mathcal U_{>j}$$

   Before the noise of $G_j$ in the $r$-th forward block come $r$ complete pairs $U,U^\dagger$ and the gates $G_1,\dots,G_j$; after it come $G_{j+1},\dots,G_d$ and $k-r$ pairs $U^\dagger,U$. Complete pairs act as the identity, so this location contributes $\epsilon\,c_j$.
4. 4
   Step 4:

   $$\text{inverse block, channel of }G_j\text{ before }G_j^\dagger:\quad\mathcal B_t=\mathcal U_{>j}^\dagger\,\mathcal U=\mathcal U_{\le j},\qquad\mathcal A_t=\mathcal U\,\mathcal U_{\le j}^\dagger=\mathcal U_{>j}$$

   In foldGlobal's order the inverse block runs $\mathcal E_d,G_d^\dagger,\dots,\mathcal E_j,G_j^\dagger,\dots$. Before $\mathcal E_j$ the ideal gates net to $U$ followed by $G_d^\dagger,\dots,G_{j+1}^\dagger$, i.e. $U_{>j}^\dagger U=U_{\le j}$; after it come $G_j^\dagger,\dots,G_1^\dagger$ and gates netting to $U$, i.e. $U\,U_{\le j}^\dagger=U_{>j}$. This location also contributes $\epsilon\,c_j$.
5. 5
   Step 5:

   $$E(\lambda)=E_0+\epsilon\Big[(k+1)\sum_jc_j+k\sum_jc_j\Big]+O(\epsilon^2)=E_0+(2k+1)\,\epsilon\sum_jc_j+O(\epsilon^2)$$

   There are $k+1$ forward blocks and $k$ inverse blocks, each containing the noise of every gate once. This proves (a).
6. 6
   Step 6:

   $$\mathcal E_j\text{ after }G_j^\dagger:\quad\mathcal B_t=\mathcal G_j^\dagger\,\mathcal U_{>j}^\dagger\,\mathcal U=\mathcal U_{<j},\qquad\mathcal A_t=\mathcal U\,\mathcal U_{<j}^\dagger=\mathcal U_{\ge j}\ \Rightarrow\ c_j'$$

   In hardware order the noise of $G_j^\dagger$ acts after it, one gate earlier in the ideal circuit's coordinates: the error sits between $G_{j-1}$ and $G_j$. Counting blocks as before gives (b).
7. 7
   Step 7:

   $$\mathcal L_j\mathcal G_j=\mathcal G_j\mathcal L_j\ \Rightarrow\ c_j'=\mathrm{Tr}\big[O\,\mathcal U_{>j}\mathcal G_j\mathcal L_j\mathcal U_{<j}(\rho_0)\big]=\mathrm{Tr}\big[O\,\mathcal U_{>j}\mathcal L_j\mathcal G_j\mathcal U_{<j}(\rho_0)\big]=c_j$$

   Use $\mathcal U_{\ge j}=\mathcal U_{>j}\mathcal G_j$ and $\mathcal G_j\mathcal U_{<j}=\mathcal U_{\le j}$. The depolarizing channel $\rho\mapsto f\rho+(1-f)\mathrm{Tr}(\rho)\tfrac I2$ commutes with every single-qubit unitary conjugation, so the condition holds for depolarizing noise after single-qubit gates.

∎

1. 1
   Step 1:

   $$\mathbb E\big[\hat E_0\big]=\sum_i\gamma_i\,\mathbb E\big[\hat E_i\big]=\sum_i\gamma_iE(\lambda_i)$$

   Linearity of expectation; each sample mean is unbiased for its node. The Richardson estimate inherits only the extrapolation bias, not a statistical one.
2. 2
   Step 2:

   $$\mathrm{Var}\Big[\sum_i\gamma_i\hat E_i\Big]=\sum_{i,j}\gamma_i\gamma_j\,\mathrm{Cov}\big[\hat E_i,\hat E_j\big]=\sum_i\gamma_i^2\,\mathrm{Var}\big[\hat E_i\big]$$

   Bilinearity of the covariance; estimates from different circuits are independent, so the off-diagonal covariances vanish.
3. 3
   Step 3:

   $$\mathrm{Var}\big[\hat E_i\big]=\mathrm{Var}\Big[\frac1{N_i}\sum_{s=1}^{N_i}o_s\Big]=\frac{\sigma_i^2}{N_i},\qquad o_s\in\{\pm1\}:\ \sigma_i^2=\mathbb E[o^2]-\mathbb E[o]^2=1-E(\lambda_i)^2$$

   The variance of a mean of $N_i$ independent, identically distributed outcomes.
4. 4
   Step 4:

   $$\sum_i\lvert\gamma_i\rvert\sigma_i=\sum_i\frac{\lvert\gamma_i\rvert\sigma_i}{\sqrt{N_i}}\sqrt{N_i}\le\Big(\sum_i\frac{\gamma_i^2\sigma_i^2}{N_i}\Big)^{1/2}\Big(\sum_iN_i\Big)^{1/2}$$

   Cauchy–Schwarz inequality for the vectors $\big(\lvert\gamma_i\rvert\sigma_i/\sqrt{N_i}\big)_i$ and $\big(\sqrt{N_i}\big)_i$.
5. 5
   Step 5:

   $$\mathrm{Var}\big[\hat E_0\big]\ge\frac{\big(\sum_i\lvert\gamma_i\rvert\sigma_i\big)^2}{N_{\text{tot}}},\qquad\text{equality}\iff\frac{\lvert\gamma_i\rvert\sigma_i}{\sqrt{N_i}}\propto\sqrt{N_i}\iff N_i\propto\lvert\gamma_i\rvert\sigma_i$$

   Square the inequality and divide by $N_{\text{tot}}=\sum_iN_i$; Cauchy–Schwarz is an equality exactly when the two vectors are parallel.

∎

1. 1
   Step 1:

   $$\det A=(1-p_{01})(1-p_{10})-p_{10}\,p_{01}=1-p_{01}-p_{10}+p_{01}p_{10}-p_{01}p_{10}=1-p_{01}-p_{10}$$

   The determinant of a $2\times2$ matrix is $ad-bc$; the products $p_{01}p_{10}$ cancel.
2. 2
   Step 2:

   $$\begin{pmatrix}a&b\\c&d\end{pmatrix}^{-1}=\frac1{ad-bc}\begin{pmatrix}d&-b\\-c&a\end{pmatrix}\ \Rightarrow\ A^{-1}=\frac1{1-p_{01}-p_{10}}\begin{pmatrix}1-p_{10}&-p_{10}\\-p_{01}&1-p_{01}\end{pmatrix}$$

   The adjugate formula with $a=1-p_{01}$, $b=p_{10}$, $c=p_{01}$, $d=1-p_{10}$; it exists exactly when $\det A\ne0$.
3. 3
   Step 3:

   $$A\,A^{-1}=\frac1{\det A}\begin{pmatrix}(1-p_{01})(1-p_{10})-p_{10}p_{01}&-(1-p_{01})p_{10}+p_{10}(1-p_{01})\\p_{01}(1-p_{10})-(1-p_{10})p_{01}&-p_{01}p_{10}+(1-p_{10})(1-p_{01})\end{pmatrix}=I$$

   Direct verification: the diagonal entries equal $\det A$ and the off-diagonal entries cancel.
4. 4
   Step 4:

   $$\mathbf 1^\top A=\mathbf 1^\top\ \Rightarrow\ \mathbf 1^\top=\mathbf 1^\top A^{-1}\ \Rightarrow\ \mathbf 1^\top\hat{\mathbf p}=\mathbf 1^\top A^{-1}\tilde{\mathbf p}=\mathbf 1^\top\tilde{\mathbf p}=1$$

   Column-stochasticity is $\mathbf 1^\top A=\mathbf 1^\top$; multiply from the right by $A^{-1}$. Explicitly, the first column of $A^{-1}$ sums to $(1-p_{10}-p_{01})/\det A=1$, and so does the second.
5. 5
   Step 5:

   $$\mathbb E\big[\hat{\mathbf p}\big]=A^{-1}\,\mathbb E\big[\hat{\tilde{\mathbf p}}\big]=A^{-1}A\,\mathbf p=\mathbf p,\qquad\tilde{\mathbf p}=(0.97,0.03)\ \Rightarrow\ \hat{\mathbf p}=(1.0235294117647058,\,-0.02352941176470588)$$

   Linearity of expectation and unbiasedness of measured frequencies. For an individual finite sample, however, $A^{-1}$ can produce negative entries, as the example with $p_{01}=0.05$, $p_{10}=0.10$ shows: a nearly pure $0$ histogram is over-corrected.