developed by Brandão _et al._ [brandao2019quantum] was shown to greatly
improve these costs. Primarily, the ability to solve SDPs on a quantum
computer bypasses the need to manipulate an exponentially large classical
representation. That being said, while this makes the space complexity
polynomial, the total number of quantum gates (hence time complexity) remains
exponential, as $\tilde{\mathcal{O}}(\sqrt{L}G)+\operatorname{poly}(2^{n})$.
These costs can be somewhat improved by placing further restrictions, for
instance by only considering POVMs that can be implemented in
$G=\operatorname{poly}(n)$ gates.

### Classical shadows

While a breakthrough in sample complexity for learning quantum states, it is
clear that the computational cost of shadow tomography is still impractical.
Even ignoring gate complexity, the problem involves coherently manipulating
the state $\rho^{\otimes N}$, which requires a massive register of
$\tilde{\mathcal{O}}(n^{2}\epsilon^{-4}\log^{4}L\log\delta^{-1})$ qubits. A
pressing question is therefore: can the philosophy of shadow tomography be
adapted into a practical, near-term-friendly protocol?

In their seminal 2020 paper, Huang, Kueng, and Preskill [huang2020predicting]
affirmed this question, introducing the method of _classical
shadows_.252525Paini _et al._ [paini2019approximate_superseded,
paini2021estimating] also introduced essentially the same idea independently,
albeit specializing to $\mathcal{D}=\mathrm{SU}(2)^{\otimes n}$ rather than
providing a completely general formulation. Importing ideas from quantum state
tomography, in particular linear-inversion estimation [guta2020fast], into the
“shadow” framework, they developed a protocol which converts single copies of
a quantum state into a minimal classical sketch, called its classical shadow.
This efficient description is designed such that different properties
statistically converge to the ground truth at different rates; for example, a
classical shadow built out of $\tilde{\mathcal{O}}(\epsilon^{-2})$ samples may
already predict local observables to within $\epsilon$ error, while
nonetheless being exponentially far from the true density matrix. Importantly,
the required gates can be as simple as single-qubit rotations.

The central idea is to craft a distribution of POVMs such that the desired
properties are the ones converged quickly. Thus the task is somewhat of an
inverse to the PAC learning problem: instead of being given $\mathcal{D}$, we
seek such a distribution which yields nice prediction features with respect to
the target set of observables. The setting can be framed as follows. Given
single-copy measurement access to the $n$-qubit state $\rho$, we wish to
predict the values of $\operatorname{tr}(O_{j}\rho)$ for a collection of
observables $O_{1},\ldots,O_{L}$. Specifically, for accuracy and confidence
parameters $\epsilon,\delta\in(0,1)$, output estimates
$\hat{o}_{1},\ldots,\hat{o}_{L}\in\mathbb{R}$ such that

$|\hat{o}_{j}-\operatorname{tr}(O_{j}\rho)|\leq\epsilon\quad\forall j\in[L],$
(119)

with success probability at least $1-\delta$. The difference from shadow
tomography is that, rather than demanding the protocol be efficient over any
possible collection of observables, we accept that it may only be efficient
over certain classes of observables (for example, local observables or low-
rank observables).

#### Warm-up example: random Pauli measurements

To build some intuition for classical shadows, we first examine a specific but
highly relevant protocol of _random Pauli measurements_. Such measurements are
important because they are arguably the simplest informationally complete set
of POVMs, requiring only local, single-qubit control to implement. Note that
random Pauli measurements have been considerably studied prior to the advent
of classical shadows, especially in the context of partial state learning
[yu2020sample, cotler2020quantum, paini2019approximate_superseded,
evans2019scalable]. The intuition we develop here draws from those works as
well.

##### A single qubit

We begin with a single qubit. A Pauli measurement is achieved by mapping the
computational basis, defined by $Z=|0\rangle\\!\langle 0|-|1\rangle\\!\langle
1|$, to one of either $X=\mathsf{H}^{\dagger}Z\mathsf{H}$,
$Y=\mathsf{G}^{\dagger}Z\mathsf{G}$, or itself,
$Z=\mathbb{I}^{\dagger}Z\mathbb{I}$. These transformations are achieved by the
unitaries

$\mathsf{H}\coloneqq\frac{1}{\sqrt{2}}\begin{pmatrix}1&1\\\
1&-1\end{pmatrix},\quad\mathsf{G}\coloneqq\mathsf{H}S^{\dagger}\text{ where
}S\coloneqq\begin{pmatrix}1&0\\\ 0&\mathrm{i}\end{pmatrix}.$ (120)

Let us first review how Pauli measurements are typically performed. As our
native POVM is limited to the $Z$ basis, predicting expectation values of an
arbitrary Pauli observable $W\in\\{X,Y,Z\\}$ makes use of the identity

$\operatorname{tr}(W\rho)=\operatorname{tr}(U^{\dagger}ZU\rho)=\operatorname{tr}(ZU\rho
U^{\dagger})$ (121)

for an appropriate choice of $U\in\\{\mathbb{I},\mathsf{H},\mathsf{G}\\}$.
Thus, by applying $U$ to the state immediately before measurement in the
native $Z$-basis, we obtain statistics of the original state in the $W$ basis.

By the postulates of quantum mechanics, each measurement outcome
$x\in\\{0,1\\}$ is obtained with probability $p_{x}=\langle x|U\rho
U^{\dagger}|x\rangle$. Using $Z|x\rangle=(-1)^{x}|x\rangle$, we can construct
an estimator for $\operatorname{tr}(W\rho)$ from

$\displaystyle\operatorname{tr}(W\rho)=\operatorname{tr}(ZU\rho U^{\dagger})$
$\displaystyle=\sum_{x\in\\{0,1\\}}p_{x}(-1)^{x}$ (122)
$\displaystyle\approx\frac{1}{N}\sum_{\ell=1}^{N}(-1)^{x_{\ell}},$ (123)

where $x_{1},\ldots,x_{N}\in\\{0,1\\}$ are the measurement outcomes of $N$
rounds of this experiment. We want the absolute error of this estimate to be
no larger than $\epsilon$:

$\mathopen{}\left|\operatorname{tr}(W\rho)-\frac{1}{N}\sum_{\ell=1}^{N}(-1)^{x_{\ell}}\right|\mathclose{}\leq\epsilon.$
(124)

To show how many samples $N$ we need for this bound to hold, one can appeal to
Chebyshev’s inequality. For a random variable $\mathcal{X}$ with finite first
and second moments, $|{\operatorname*{\mathbb{E}}[\mathcal{X}]}|<\infty$ and
$\operatorname{Var}[\mathcal{X}]<\infty$, Chebyshev’s inequality states that
the probability that $\mathcal{X}$ deviates from its mean value
$\operatorname*{\mathbb{E}}[\mathcal{X}]$ by more than $\epsilon$ obeys

$\Pr[|\mathcal{X}-\operatorname*{\mathbb{E}}[\mathcal{X}]|\geq\epsilon]\leq\frac{\operatorname{Var}[\mathcal{X}]}{\epsilon^{2}}.$
(125)

This is the probability that our experiment has failed to produce
$\epsilon$-accurate estimates, so we want it to be small. In our setting, the
random variable is simply the average of the measurement outcomes:

$\mathcal{X}=\frac{1}{N}\sum_{\ell=1}^{N}(-1)^{x_{\ell}},$ (126)

which obeys
$\operatorname*{\mathbb{E}}[\mathcal{X}]=\operatorname{tr}(W\rho)$. Because
the samples are independent and identically distributed (i.i.d.), we can
compute its variance as

$\begin{split}\operatorname{Var}[\mathcal{X}]=\operatorname{Var}\mathopen{}\left[\frac{1}{N}\sum_{\ell=1}^{N}(-1)^{x_{\ell}}\right]\mathclose{}=\frac{1}{N^{2}}\sum_{\ell=1}^{N}\operatorname{Var}\mathopen{}\left[(-1)^{x_{\ell}}\right]\mathclose{}&=\frac{1}{N}\operatorname{Var}\mathopen{}\left[(-1)^{x}\right]\mathclose{}\\\
&=\frac{1}{N}\mathopen{}\left[1-\operatorname{tr}(W\rho)^{2}\right]\mathclose{}.\end{split}$
(127)

Assuming no prior information about $\rho$, we can apply the trivial bound
$\operatorname{tr}(W\rho)^{2}\geq 0$ to bound the variance as
$\operatorname{Var}[\mathcal{X}]\leq 1/N$. Chebyshev’s inequality thus implies
that the failure probability is at most $\delta=1/(N\epsilon^{2})$.
Rearranging this statement informs us that taking $N=1/(\delta\epsilon^{2})$
samples suffices262626More sophisticated tail bounds, such as Bernstein’s
inequality, can exponentially improve the dependence on $\delta$ to the
optimal scaling of $\log(1/\delta)$. We stick with Chebyshev’s inequality here
for simplicity, and besides we take $\delta$ to be constant. to learn
$\operatorname{tr}(W\rho)$ to error $\leq\epsilon$, with success probability
$\geq 1-\delta$. Doing this for all three Pauli matrices means we need to a
total of $N\to 3/(\delta\epsilon^{2})$ measurements. Taking $\delta$ to be
some small constant, say $0.01$, implies that the statistical error scales as
$\epsilon\sim\sqrt{3/N}$.

Now we invert the approach: rather than choosing the Pauli matrix
$W\in\\{X,Y,Z\\}$ that we want to measure and selecting the appropriate
$U\in\\{\mathbb{I},\mathsf{H},\mathsf{G}\\}$, we _first_ select the unitary
and then ask which Pauli basis it rotates $Z$ into. This perspective is
perfectly valid if we are interested in learning all three Pauli observables,
as we will eventually need to choose all three unitaries anyways. Since we
need to perform many rounds of the experiment to collect ample statistics in
the first place, we can allocate the three rotations uniformly among the $N$
rounds.

Suppose instead that we have decision anxiety about which unitary to implement
at any given round of the experiment.272727For example, perhaps we can only
run the experiment for an unknown amount of time before funding runs out and
power to the lab gets cut off. Fortunately, we can appeal to classical
randomness to make the decision for us: for each $\ell\in[N]$, flip a fair
three-sided coin282828We leave the construction of such an object as an
exercise to the reader. to decide which $U$ to apply. This promotes $U$ to a
random variable $U_{\ell}$, $\ell=1,\ldots,N$, with uniform probabilities
$\Pr[U_{\ell}=\mathbb{I}]=\Pr[U_{\ell}=\mathsf{H}]=\Pr[U_{\ell}=\mathsf{G}]=1/3$.
This defines our distribution $\mathcal{D}$ of POVMs.

For clarity of exposition, suppose we want to estimate
$\operatorname{tr}(X\rho)$ from the $N$ samples of this _randomized_
measurement scheme. (By uniformity, the following argument will apply to $Y$
and $Z$ as well.) On average, only $N/3$ samples will contain any information
about $X$, which we assign a value of $(-1)^{x_{\ell}}$ as usual. The other
$2N/3$ experiments do not measure in the $X$ basis (on average), and so are
assigned a value of $0$. Thus in expectation,

$\operatorname*{\mathbb{E}}_{U_{\ell}\sim\mathcal{D}}[\operatorname{tr}(U_{\ell}^{\dagger}ZU_{\ell}\rho)]=\frac{1}{3}\operatorname{tr}(\mathsf{H}^{\dagger}Z\mathsf{H}\rho)+\frac{2}{3}\cdot
0=\frac{1}{3}\operatorname{tr}(X\rho).$ (128)

To correct for this factor of $1/3$, we merely need to redefine the estimator
as

$\hat{o}_{\ell}=\begin{cases}3(-1)^{x_{\ell}}&\text{if
}U_{\ell}=\mathsf{H},\\\ 0&\text{else},\end{cases}$ (129)

which obeys
$\operatorname*{\mathbb{E}}[\hat{o}_{\ell}]=\operatorname{tr}(X\rho)$. To get
an accurate estimate, average over all $N$ samples,

$\hat{o}(N)\coloneqq\frac{1}{N}\sum_{\ell=1}^{N}\hat{o}_{\ell},$ (130)

which by linearity still obeys
$\operatorname*{\mathbb{E}}[\hat{o}(N)]=\operatorname{tr}(X\rho)$. The
convergence of this estimator also scales as $\epsilon\sim\sqrt{3/N}$; we can
show this using the same Chebyshev argument from before. The only modification
we need to make is to compute the variance of $\hat{o}(N)$. As the samples are
i.i.d. (including the random draws of $U_{\ell}$), we have

$\begin{split}\operatorname{Var}[\hat{o}(N)]&=\frac{1}{N}\operatorname{Var}[\hat{o}(1)]\\\
&=\frac{1}{N}\mathopen{}\left(\frac{1}{3}\cdot 3^{2}+\frac{2}{3}\cdot
0^{2}\right)\mathclose{}-\frac{1}{N}\operatorname{tr}(X\rho)^{2}\\\
&\leq\frac{3}{N}.\end{split}$ (131)

Note that because all $N$ samples are used to estimate all three Pauli matrix
expectation values, this is the final variance bound. Hence, we have
sidestepped our decision anxiety through the use of randomization, while
enjoying the same statistical guarantees as the deterministic approach.

##### Many qubits

It turns out that this sidestepping of decision making has fundamental
consequences for constructing an efficient protocol. To see this, we must
generalize to $n\geq 1$ qubits. Fortunately, it is simple to generalize the
distribution $\mathcal{D}$ from the single-qubit setting: simply perform
random Pauli measurements on each qubit independently. Estimating any single-
qubit Pauli observable, for example $X\otimes\mathbb{I}^{\otimes(n-1)}$,
proceeds as described above; we merely have to keep track of qubit labels for
the random unitaries $U=U^{(1)}\otimes\cdots\otimes U^{(n)}$ and measurement
outcomes $b=b^{(1)}\cdots b^{(n)}\in\\{0,1\\}^{n}$ that we acquire from each
experimental round. We can estimate such observables in parallel over the $n$
qubits.

Since we have a composite system, we would like to learn not only single-body,
but also many-body properties. For example, consider the set of two-local
operators, $\mathbb{I}^{\otimes(i-1)}\otimes
W^{(i)}\otimes\mathbb{I}^{\otimes(j-i-1)}\otimes
W^{(j)}\otimes\mathbb{I}^{\otimes(n-j)}$, where $1\leq i<j\leq n$ and
$W^{(i)},W^{(j)}\in\\{X,Y,Z\\}$. In this scenario, a randomly selected
$n$-qubit Pauli measurement “hits” $W^{(i)}$ and $W^{(j)}$ with probability
$1/9$. Then, the appropriate estimator that we should define will be analogous
to Eq. 129, but with a factor of $9$ instead of $3$:

$\hat{o}_{\ell}=\begin{cases}9(-1)^{b_{\ell}^{(i)}+b_{\ell}^{(j)}}&\text{if
}(U_{\ell}^{(i)})^{\dagger}ZU_{\ell}^{(i)}=W^{(i)}\text{ and
}(U_{\ell}^{(j)})^{\dagger}ZU_{\ell}^{(j)}=W^{(j)},\\\
0&\text{else}.\end{cases}$ (132)

By the same analysis, the error of this estimate scales as $\epsilon\sim
3/\sqrt{N}$.

Now we make the key observation that reveals the power of randomization. For
each $n$-qubit Pauli basis, there are $\binom{n}{2}$ two-local Pauli
observables that are “hit,” i.e., acquire a nonzero estimate according to Eq.
132. On the other hand, there are $9\binom{n}{2}$ two-local Pauli operators in
total. Thus we can conceptually lift our probabilistic interpretation of
hitting individual operators to the hitting of collections of operators: with
probability $1/9$, a whole collection of $\binom{n}{2}$ operators are learned
by a single choice of random measurement basis. By a coupon-collector
argument, one may expect that on the order of ${\sim}(9\log 9)$ random draws
of measurement bases will result in every two-local Pauli operator being hit
at least once.292929This argument is not entirely valid, because commutation
relations only allow certain collections of operators to be grouped together.
Nonetheless, the Pauli matrices are sufficiently uniform such that this
intuition at least morally holds. Thus, up to some logarithmic corrections,
taking $N\sim 9/\epsilon^{2}$ random Pauli measurements suffices to obtain
$\epsilon$-accurate estimates of all $9\binom{n}{2}$ two-local Pauli
observables. Very crucially, this sample complexity is independent of $n$,
despite being able to learn quadratically many observables!

In the absence of randomization, if one were to learn each operator
individually (say, by allocating each experimental round to only estimating a
single two-local Pauli observable), this would require taking $N\sim
9\binom{n}{2}/\epsilon^{2}$ measurements. But this is far too naive: as we
have already pointed out, each measurement basis hits many operators. The
question then is, can one construct a deterministic set of Pauli measurements
which covers all two-local operators with only around $9$ Pauli POVMs? This is
certainly possible in principle; however, we need to construct this cover
efficiently for any $n$. What about covering the set of three-local Paulis
instead—do we need an entirely new construction algorithm? Four-local? What if
we had performed the experiment with our measurement scheme optimized for two-
local operators, but then in retrospect we now want to extract some three-body
information from that data without running more experiments?

The elegance of randomized measurement schemes, and classical shadows in
particular, is that all these challenges and concerns are sidestepped.
Instead, an appropriately chosen distribution of random measurements does much
of the heavy lifting for us, and we are left in the comfortable position of
“measure first, ask questions later” [elben2023randomized]. Furthermore, this
approach is fully rigorous: using statistical analysis, we can control the
estimation accuracy and bound our confidence intervals.

#### Description of the general theory

This section provides an overview of full theory of classical shadows,
generalizing beyond Pauli measurements [huang2020predicting]. The randomized
measurement primitive can described as follows: fix an ensemble of $n$-qubit
unitaries $\mathcal{D}\subseteq\mathrm{U}(2^{n})$ equipped with a probability
distribution. For each copy of $\rho$, draw an independently random
$U\sim\mathcal{D}$ and apply it to $\rho$. Then measure the state $U\rho
U^{\dagger}$, obtaining a bit string $b\in\\{0,1\\}^{n}$ with probability
$\langle b|U\rho U^{\dagger}|b\rangle$. Finally, store a classical description
of the postmeasurement state with the inverse unitary applied,
$U^{\dagger}|b\rangle\\!\langle b|U$.

This measurement primitive implements the POVM
$\\{p(U)U^{\dagger}|b\rangle\\!\langle b|U\mid
U\in\mathcal{D},b\in\\{0,1\\}^{n}\\}$, where $p(U)$ is the probability measure
on $\mathcal{D}$. In expectation, the quantum process is a channel (completely
positive trace-preserving map):

$\mathcal{M}(\rho)=\operatorname*{\mathbb{E}}_{U\sim\mathcal{D}}\sum_{b\in\\{0,1\\}^{n}}\langle
b|U\rho U^{\dagger}|b\rangle U^{\dagger}|b\rangle\\!\langle b|U.$ (133)

One can also view this as the partial trace of a two-fold twirl by
$\mathcal{D}$:

$\mathcal{M}(\rho)=\operatorname{tr}_{1}\mathopen{}\left[\sum_{b\in\\{0,1\\}^{n}}\mathcal{T}_{2,\mathcal{D}}(|b\rangle\\!\langle
b|^{\otimes 2})(\rho\otimes\mathbb{I})\right]\mathclose{},$ (134)

where the $t$-fold twirl is defined as

$\mathcal{T}_{t,\mathcal{D}}(A)\coloneqq\operatorname*{\mathbb{E}}_{U\sim\mathcal{D}}(U^{\dagger})^{\otimes
t}AU^{\otimes t}$ (135)

for any $A\in(\mathbb{C}^{2^{n}\times 2^{n}})^{\otimes t}$. Alternatively,
from a representation of superoperator composition, one can write

$\mathcal{M}=\operatorname*{\mathbb{E}}_{U\sim\mathcal{D}}\mathcal{U}^{\dagger}\mathcal{M}_{Z}\mathcal{U},$
(136)

where $\mathcal{U}(\cdot)\coloneqq U(\cdot)U^{\dagger}$ and
$\mathcal{M}_{Z}(\cdot)\coloneqq\sum_{b\in\\{0,1\\}^{n}}\langle
b|(\cdot)|b\rangle|b\rangle\\!\langle b|$.

If the effective POVM is informationally complete, then it is possible to
invert $\mathcal{M}$,303030If not, then one can still define a pseudoinverse
which inverts over the span of the POVM elements and annihilates the kernel.
which leads to the seemingly trivial statement

$\rho=\mathcal{M}^{-1}\mathopen{}\left(\mathcal{M}(\rho)\right)\mathclose{}.$
(137)

However, from the sampling perspective of $\mathcal{M}$ according to Eq. 133,
we can use linearity to write

$\rho=\operatorname*{\mathbb{E}}_{U\sim\mathcal{D}}\operatorname*{\mathbb{E}}_{b\sim
U\rho U^{\dagger}}\mathcal{M}^{-1}(U^{\dagger}|b\rangle\\!\langle b|U),$ (138)

where $b\sim U\rho U^{\dagger}$ denotes draws of $|b\rangle$ according to the
Born rule. Thus each sample $(U,b)$ of the measurement primitive is converted
to a _classical shadow snapshot_
$\hat{\rho}_{U,b}\coloneqq\mathcal{M}^{-1}(U^{\dagger}|b\rangle\\!\langle
b|U)$, called such because it is an unbiased estimator for $\rho$. That is,
$\hat{\rho}_{U,b}$ obeys
$\operatorname*{\mathbb{E}}_{U,b}[\hat{\rho}_{U,b}]=\rho$. While
$\mathcal{M}^{-1}$ is not a physical process (it is not completely positive),
the postmeasurement state $U^{\dagger}|b\rangle\\!\langle b|U$ is already
stored as some classical description, so applying $\mathcal{M}^{-1}$ amounts
to further classical computation. This can be achieved, for example, by
computing an analytic expression for $\mathcal{M}$ using tools from
representation theory whenever $\mathcal{D}$ is a compact group with the Haar
measure. See Chapter Group-theoretic Error Mitigation Enabled by Classical
Shadows and Symmetries for further details.

Because classical shadows reproduce the density matrix in expectation, by
linearity they can be used to predict expectation values:

$\operatorname{tr}(O_{j}\rho)=\operatorname*{\mathbb{E}}_{U,b}\operatorname{tr}(O_{j}\hat{\rho}_{U,b}).$
(139)

To study how quickly these estimates converge to the ground truth, we can
analyze their variance. Let $T=KN$ be the total number of samples,
$\hat{\rho}_{1},\ldots,\hat{\rho}_{T}$, partitioned into $K$ groups of $N$
snapshots each. Huang _et al._ [huang2020predicting] employ the median-of-
means estimator

$\hat{o}_{j}(T,K)\coloneqq\mathrm{median}\mathopen{}\left\\{\operatorname{tr}(O_{j}\hat{\rho}_{(1)}),\ldots,\operatorname{tr}(O_{j}\hat{\rho}_{(K)})\right\\}\mathclose{},$
(140)

where each $\hat{\rho}_{(k)}$ is defined as the average classical shadow of
the $k$th batch,

$\hat{\rho}_{(k)}\coloneqq\frac{1}{N}\sum_{\ell=(k-1)N+1}^{kN}\hat{\rho}_{\ell}.$
(141)

An application of the Chebyshev and Hoeffding bounds implies that a sample
complexity of

$T=\mathcal{O}\mathopen{}\left(\epsilon^{-2}\log(L/\delta)\max_{j\in[L]}\operatorname{Var}[\hat{o}_{j}]\right)\mathclose{}$
(142)

suffices to solve the problem stated in Eq. 119. The variance appearing here
is the single-shot variance of $\hat{o}_{j}\equiv\hat{o}_{j}(1,1)$, defined in
the standard way:

$\begin{split}\operatorname{Var}[\hat{o}_{j}]&=\operatorname*{\mathbb{E}}_{U,b}\mathopen{}\left[\operatorname{tr}(O_{j}\hat{\rho}_{U,b})^{2}\right]\mathclose{}-\operatorname{tr}(O_{j}\rho)^{2}\\\
&=\operatorname*{\mathbb{E}}_{U\sim\mathcal{D}}\sum_{b\in\\{0,1\\}^{n}}\langle
b|U\rho U^{\dagger}|b\rangle\langle
b|U\mathcal{M}^{-1}(O_{j})U^{\dagger}|b\rangle^{2}-\operatorname{tr}(O_{j}\rho)^{2}\\\
&=\operatorname{tr}\mathopen{}\left[\sum_{b\in\\{0,1\\}^{n}}\mathcal{T}_{3,\mathcal{D}}(|b\rangle\\!\langle
b|^{\otimes 3})\mathopen{}\left(\rho\otimes\mathcal{M}^{-1}(O_{j})^{\otimes
2}\right)\mathclose{}\right]\mathclose{}-\operatorname{tr}(O_{j}\rho)^{2}.\end{split}$
(143)

Above, we used the fact that $\mathcal{M}$ and its inverse are self-adjoint,
i.e.,
$\operatorname{tr}(A\mathcal{M}^{-1}(B))=\operatorname{tr}(\mathcal{M}^{-1}(A)B)$.
Observe that the variance depends on the third moment of the distribution
$\mathcal{D}$, as captured by the three-fold twirl
$\mathcal{T}_{3,\mathcal{D}}$.

Using the universal bound $\operatorname{tr}(O_{j}\rho)^{2}\geq 0$ and taking
the supremum over all quantum states yields a state-independent variance bound
called the _shadow norm_ :

$\begin{split}\|O_{j}\|_{\mathrm{shadow}}&\coloneqq\max_{\text{states
}\sigma}\sqrt{\operatorname*{\mathbb{E}}_{U\sim\mathcal{D}}\sum_{b\in\\{0,1\\}^{n}}\langle
b|U\rho U^{\dagger}|b\rangle\langle
b|U\mathcal{M}^{-1}(O_{j})U^{\dagger}|b\rangle^{2}}\\\
&\geq\operatorname{Var}[\hat{o}_{j}].\end{split}$ (144)

Indeed, the shadow norm obeys all the properties of a norm on the space of
linear operators. This quantity is particularly useful for evaluating the
efficiency of the classical shadows protocol as it is valid for all quantum
states, and it can be (in principle) computed before running the experiment,
given only knowledge of the user-specified inputs to the protocol (namely, the
unitary ensemble $\mathcal{D}$ and the observables $O_{j}$ of interest).

The two examples considered by [huang2020predicting] were the local Clifford
group $\mathrm{Cl}(1)^{\otimes n}$ and the global $n$-qubit Clifford group
$\mathrm{Cl}(n)$. For the former, they showed that $k$-local observables
$O_{j}$ exhibit a shadow norm $\|O_{j}\|_{\mathrm{shadow}}^{2}\leq
4^{k}\|O_{j}\|_{\infty}$, which simplifies to
$\|O_{j}\|_{\mathrm{shadow}}^{2}=3^{k}$ whenever $O_{j}$ is a single $k$-local
Pauli operator. Thus this protocol, essentially equivalent to the random Pauli
measurements highlighted in Warm-up example: random Pauli measurements, is
highly efficient for estimating local qubit observables. For the latter
distribution, the variance instead scales with the Frobenius norm of the
operator, $\|O_{j}\|_{\mathrm{shadow}}^{2}\leq 3\operatorname{tr}(O_{j}^{2})$,
which is favorable for low-rank observables such as
$|\psi\rangle\\!\langle\psi|$ (e.g., as in fidelity estimation).

In Chapter Fermionic Partial Tomography via Classical Shadows we will
demonstrate another distribution, essentially of fermionic Gaussian unitaries,
which achieves optimal shadow norm
$\|O_{j}\|_{\mathrm{shadow}}^{2}=\binom{2n}{2k}/\binom{n}{k}$ for $k$-body
fermionic (Majorana) observables [zhao2021fermionic].

Altogether we see that, unlike in shadow tomography, the scaling with respect
to $1/\epsilon$ is only quadratic, which is much more manageable (and is in
fact optimal without additional quantum resources or special assumptions). The
scaling with the number of target observables $L$ is also only logarithmic, as
is with $1/\delta$, both of which are optimal.

It is a somewhat subtle point that median-of-means estimation was necessary in
order to achieve this optimal scaling with $L/\delta$, as otherwise (say, with
a standard mean estimator) a union bound over all failure probabilities would
result in a sampling bound linear in $L/\delta$, not logarithmic. This would
be unacceptably loose, as there are many applications for which we would like
to take $L=\operatorname{poly}(n)$ (for instance, local observable
estimation). In Chapter Fermionic Partial Tomography via Classical Shadows,
Performance guarantees without median-of-means estimation, we touch on this
subject, showing that under certain scenarios one can in fact apply a
Bernstein inequality to the mean estimator to nonetheless achieve the desired
$\log(L/\delta)$ scaling. [helsen2022thrifty] also examined this point in some
more detail, showing that median-of-means estimation can also be successfully
replaced with the mean estimator when $\mathcal{D}=\mathrm{U}(2^{n})$, the
full unitary group, but not when $\mathcal{D}=\mathrm{Cl}(n)$, the $n$-qubit
Clifford group. This is a particularly surprising result because the Clifford
group is a unitary $3$-design, and as such has identical estimation formulas
and variances to $\mathrm{U}(2^{n})$ itself. This discrepancy is due to the
difference in the higher moments ($t>3$) of the distributions, which causes
the tails of $\hat{o}_{j}$ to be exponentially bounded under
$\mathrm{U}(2^{n})$ but heavy under $\mathrm{Cl}(n)$. Unfortunately, while
Clifford circuits can be efficiently constructed from $\mathcal{O}(n^{2})$
gates [aaronson2004improved, bravyi2021hadamard], Haar-random unitaries
require exponentially many gates [knill1995approximation]. Thus it seems that
the nicer tail properties of $\mathrm{U}(2^{n})$ cannot be efficiently taken
advantage of.

#### Learning nonlinear functions

As an aside, it turns out that classical shadows yield more than what was
bargained for: they can also be used to predict nonlinear functions of $\rho$.
This is accomplished by manipulating the classical shadow $\hat{\rho}$
appropriately. For example, any quadratic function, such as
$\operatorname{tr}(O\rho^{2})$, can be expressed as the expectation value of
some $2n$-qubit operator $O^{\prime}$ on two copies of the state:
$\operatorname{tr}(O\rho^{2})=\operatorname{tr}(O^{\prime}\rho^{\otimes 2})$.
Under this equivalence, the $T$ classical shadows can be used to construct
estimates for $\rho^{\otimes 2}$, and hence the quadratic function, as

$\frac{1}{T(T-1)}\sum_{\ell\neq\ell^{\prime}\in[T]}\operatorname{tr}(O^{\prime}(\hat{\rho}_{\ell}\otimes\hat{\rho}_{\ell^{\prime}}))\approx\operatorname{tr}(O\rho^{2}).$
(145)

The estimator is unbiased due to linearity and independence of each
$\hat{\rho}_{\ell}$.

This aspect of classical shadows further broadens its appeal, as quantities
such as the $\alpha$-Rényi entanglement entropy
[$\log(\operatorname{tr}(\rho^{\alpha}))/(1-\alpha)$] are both of high
importance to many-body physics and also experimentally difficult to learn
with standard techniques [elben2020mixed]. Although we will only be interested
in predicting linear functions within this dissertation, further analysis of
the nonlinear scenario can be found in the original paper
[huang2020predicting].

## Noisy Quantum Computation

One of the greatest challenges facing the development and deployment of large-
scale quantum processors is their susceptibility to noise. By virtue, quantum
computers manipulate delicate states of Nature which are fragile to the
effects of their harsh environment. Thus while quantum computation is
theoretically believed to be more powerful than classical computation
(expressed as the complexity-theoretic conjecture that the inclusion
$\mathsf{BPP}\subset\mathsf{BQP}$ is strict), the task of building a quantum
machine that is fully robust to its environment is both a necessary yet
daunting one.

In light of such challenges, the clearest path forward is a paradigm known as
quantum error correction [lidar2013quantum]. Shor first codified the idea in
1995 [shor1995scheme], which was soon followed by his celebrated threshold
theorem [shor1996fault]. These results established the ability of a quantum
computer, with access only to noisy components, to achieve _fault tolerance_
—the ability to correct errors faster than the rate at which they occur. His
theorem set an upper bound on the required physical error rate, known as the
_threshold_ , for fault tolerance; above this threshold, errors would occur
too frequently for a quantum error-correction code to keep up with. Initially,
Shor’s threshold required that each noisy gate have error rate at most
$\mathcal{O}(1/\operatorname{poly}(\log t))$ to fault-tolerantly run a quantum
circuit of $t$ gates. Within about a year, others improved this threshold to
$\mathcal{O}(1)$, i.e., an error rate per gate independent of the size of the
circuit [aharonov1997fault, knill1998resilient_rspa, kitaev2003fault]. This
strengthened threshold solidified the promise of quantum error correction to
enable universal, scalable quantum computation [gottesman1998theory].

Unfortunately, current quantum technology is not yet below the required
threshold for promising error-correction codes, such as the surface code
[fowler2012surface]. At the time of this writing, substantial experimental
progress has been made in demonstrating the essential components for a fault-
tolerant, logically encoded qubit [egan2021fault, postler2022demonstration,
zhao2022realization, sundaresan2023demonstrating, google2023suppressing,
sivak2023real, ni2023beating]. However, device imperfections and limited qubit
counts ultimately impede the realization of fault-tolerant quantum computation
in the near term. In the face of this chasm, researchers began to investigate
the capabilities of noisy devices in the absence of quantum error correction,
a time frame coined by Preskill as the noisy intermediate-scale quantum (NISQ)
era [preskill2018quantum]. The broad motivation is that, already at a modest
number of qubits, say $n=100$, the corresponding Hilbert space of complex
dimension $2^{n}\approx 10^{30}$ is well beyond the reach of being simulated
by the world’s most powerful supercomputers. For example, strongly correlated
quantum systems can be studied using classical algorithms on up to $50\sim 70$
qubits, beyond which one must typically resort to approximate methods which
may fail to provide sufficient accuracy [wecker2014gate]. The hope then is
that quantum devices of such size, despite the imperfections and errors, might
be able to outperform the best-known classical methods.

### Quantum supremacy and advantage

There are two tall asks associated with demonstrating such a quantum
advantage. First, because we can always compare against existing classical
approximations, even if the the quantum algorithm could solve the problem
exactly, we must assess how much accuracy is retained in the presence of
decoherence and errors. Second, we are ultimately seeking a meaningful
application of NISQ computers, i.e., the target problem should be one of
potential scientific or industrial relevance.

To date, there have been few experimental demonstrations of so-called _quantum
supremacy_ —a computational task that would take up to multiple millennia for
the fastest classical supercomputers to solve, but at most a few hours for the
quantum machine. The most salient example is _random circuit sampling_ (RCS),
first experimentally demonstrated in 2019 by Google [arute2019quantum] and
then in 2021 by USTC [wu2021strong]. Both experiments ran circuits of up to 20
cycles (layers of random single- and two-qubit gates) with 53 and 56 qubits,
respectively. Both groups also recently reran their RCS experiments at larger
sizes with improved devices (24 cycles with 70 and 60 qubits, respectively)
[morvan2023phase, zhu2022quantum].

The idea of RCS is that its hardness is derived from the hardness of
simulating the quantum computer itself [aaronson2017complexity,
boixo2018characterizing, bouland2019complexity]. The goal is essentially to
output samples consistent with a distribution of random quantum circuits. Thus
the task is designed such that the quantum computer can address it natively.
On the other hand, the distribution is crafted so that the best-known
classical algorithm is simply to simulate it with exponential cost. Even with
such an inherent head start, the largest and most recent iteration of this
experiment in 2023 (70 qubits at 24 cycles) [morvan2023phase] exhibited a
fidelity of only $F\sim 0.0017$. That paper also further improved the best-
known classical algorithms for simulating such low-fidelity RCS experiments
[pan2022solving], from which they estimated that it would take roughly
$4.4\times 10^{27}$ floating-point operations to classically simulate RCS with
the same $F$. On the Frontier supercomputer at full capacity (ignoring issues
of insufficient memory), they translated this to a wall-clock time of 3.3
millennia.

Regardless of its ultimate fate as a demonstration of quantum supremacy, RCS
is is widely believed to be too artificially constructed to be relevant to
meaningful problems. With that in consideration, researchers have instead
looked toward problems such as combinatorial optimization, machine learning,
and quantum simulation as potential avenues [bharti2022noisy]. The overarching
strategy to tackle such problems is the use of short, parametrized quantum
circuits to variationally optimize a cost function. This heuristic is broadly
classified under variational quantum algorithms, which we discuss below.

### Variational quantum algorithms

Noisy-device fidelities rapidly degrade with the length of the quantum
circuit. Thus, there is a strong restriction on the quantum complexity on a
noisy machine before its output becomes entirely meaningless. At the same
time, many challenging computational tasks can be framed as optimization
problems, e.g., searching for the minimum eigenvalue of some Hermitian
operator $H$:

$E^{\star}=\min_{\text{states }|\psi\rangle}\langle\psi|H|\psi\rangle.$ (146)

For example, when $H$ is the Hamiltonian of a quantum system, then $E^{\star}$
is its ground-state energy. If $H$ is a classical (i.e., diagonal)
Hamiltonian, then the problem corresponds to finding an optimal combinatorial
solution.

The space of all quantum states is exponentially large. _Variational quantum
algorithms_ (VQAs) therefore restrict this search space via a parametrized
quantum circuit (PQC) $U(\bm{\theta})$, where
$\bm{\theta}\in\Theta\subseteq\mathbb{R}^{m}$ is a vector of $m$ real
parameters. For instance, the parameters are typically encoded as the angles
of the gates making up the PQC. This restricts the problem to finding

$E(\bm{\theta}^{\star})=\min_{\bm{\theta}\in\Theta}\langle\psi(\bm{\theta})|H|\psi(\bm{\theta})\rangle,$
(147)

where, without loss of generality, $|\psi(\bm{\theta})\rangle\coloneqq
U(\bm{\theta})|0^{n}\rangle$.

To offload arithmetic complexity from the quantum device, VQAs perform the
minimization via a quantum–classical feedback loop. This involves calls to the
quantum computer only to evaluate (measure) the cost function $E(\bm{\theta})$
by executing the PQC. Meanwhile, the navigation through parameter space via
estimated energy gradients is performed on a classical computer (e.g., by the
Nelder–Mead simplex method). The philosophy behind VQAs is to identify the
essential “quantum hardness” of a problem, which the PQC handles; the rest of
the algorithm is handed to a reliable classical processor which we assume can
do arbitrary amounts of polynomial-time computation. By the variational
principle, we always have that $E(\bm{\theta}^{\star})\geq E^{\star}$, with
equality only when the PQC is expressible enough to capture a ground state
within its parameter space.

While it is known that the ground-state problem is $\mathsf{QMA}$-hard in
general [kempe2006complexity] (hence intractable even for a quantum computer),
the hope is that either some Hamiltonian of interest has enough structure to
avoid that complexity, or otherwise PQCs can achieve higher-quality
approximations than competing classical algorithms in practice. Indeed, this
was the motivation behind the first proposed VQA, the _variational quantum
eigensolver_ (VQE) for approximating ground-state energies of quantum-chemical
systems [peruzzo2014variational]. To date, a variety of experiments have been
performed to investigate the performance of VQE in state-of-the-art quantum
processors [omalley2016scalable, kandala2017hardware, colless2018computation,
dumitrescu2018cloud, hempel2018quantum, kandala2019error, kokail2019self,
nam2020ground, arute2020hartree, stanisic2022observing, motta2023quantum,
obrien2023purification]; however, neither the system sizes nor hardware
performance are yet sufficient to demonstrate a quantum advantage in this
domain.

The framework of VQAs has been imported into a variety of other domains, such
as the quantum approximate optimization algorithm (QAOA313131This acronym was
later co-opted and reworked into an generalized framework called the “quantum
alternating operator ansatz” [hadfield2019quantum].) [farhi2014quantum] for
combinatorial optimization (which is $\mathsf{NP}$-hard in general) and
quantum neural networks for machine-learning tasks [mitarai2018quantum]. While
differing in scope and precise methodology, all such variational approaches
face a number of universal challenges:

  1. 1.

Ansatz expressibility. Clearly, we desire PQCs whose parameter space actually
contains a good approximation to the problem solution. This is essentially an
educated guess, or _ansatz_ , based on perhaps some physical intuition or
symmetry. Expressibility, however, is hindered by the fact that we desire
compact PQCs to avoid the overwhelming accumulation of device errors. Broadly
speaking, $\operatorname{poly}(\log n)$-depth circuits may already be
sufficient to demonstrate quantum advantage, since current classical
algorithms for simulating arbitrary quantum circuits are only efficient at
constant depth [bravyi2021classical]). This is however a subtle issue to
consider, because certain classes of quantum circuits may have
$\operatorname{poly}(n)$ gates yet be simulable by special techniques (e.g.,
Clifford [gottesman1998heisenberg] or matchgate [valiant2001quantum]
circuits).

  2. 2.

Barren plateaus. It was first realized in [mcclean2018barren] that the
parameter landscape for common cost functions (e.g., local Hamiltonians)
exhibits _barren plateaus_ : large regions of parameter space wherein the
gradients in any direction are exponentially small. This effectively renders
the search for the optimal parameters intractable. The problem of barren
plateaus has been investigated in a wide number of scenarios, revealing close
relations to a number of other PQC properties such as circuit depth and cost-
function locality [cerezo2021cost], the presence of noise [wang2021noise], and
trade-offs with ansatz expressibility [holmes2022connecting]. On the other
hand, the barren-plateaus phenomenon is a statement made on average over the
parameter landscape, and it has been proposed that good initial guesses might
circumvent the problem [zhang2022escaping]. Alternatively, diagnostics can be
probed to detect and potentially avoid the encounter of a barren plateau
[sack2022avoiding].

  3. 3.

Number of measurements. The actual value of any $E(\bm{\theta})$ cannot be
extracted exactly from a quantum computer. Instead, it must be learned from
sampling the quantum circuit, as discussed in Chapter Learning from Quantum
Systems. This can be “astronomically large” [wecker2015progress]; for example,
the electronic-structure Hamiltonian naively decomposes into
$\mathcal{O}(n^{4})$ simple-to-measure terms. Estimating gradients of the cost
function also has this issue; the analytic-gradient method
[schuld2019evaluating] is a numerically stable approach (in contrast to finite
difference), but requires measuring even more properties from the quantum
state. Note that these costs are further magnified by the fact that we require
this large number of measurements at each step of the variational
minimization.

  4. 4.

Noisy quantum hardware. Even ignoring all the above challenges, the effects of
noise in non-error-corrected devices impose a floor to the ultimate accuracy
of the VQA. While one can always hope to build hardware with less
imperfections and more precise control, one question to ask is what kinds of
_algorithmic_ protocols can be applied, to bootstrap more accurate quantum
calculations from an inherently noisy machine.

Quantum error mitigation was developed to address this final point, and below
we discuss this topic in closer detail.

### Quantum error mitigation

Quantum error correction aims to detect and, as the name suggests, correct
errors throughout the course of the computation in a coherent fashion. In
contrast, _quantum error mitigation_ (QEM) works by running the noisy device
at or above its inherent error rate, over an ensemble of various circuits. The
errors are not coherently corrected, but rather the effects of the errors
(say, in the estimation of observables) are approximately countered by
appropriate classical postprocessing of the outcomes from this ensemble of
circuit runs.

Different QEM strategies have been developed, although they share a common
penalty: the variance of error-mitigated estimation increases, thereby
requiring more samples compared to a noiseless protocol. Generally speaking,
the goal is to decrease the bias of estimation at the cost of higher variance,
so that although it takes longer to converge to the mean, this mean is closer
to the output of an ideal quantum computer than what the noisy device can
produce on its own. How this is accomplished, however, varies from strategy to
strategy. The best way to get a sense of what QEM looks like is by examining a
few of the most popular ideas. This is by no means a comprehensive survey, and
we recommend the interested reader to see [cai2022quantum] for an extensive
review article on QEM.

In all the examples below, we take the task to be the estimation of the
expectation value $\operatorname{tr}(O\rho)$ of some observable $O$ with
respect to an ideally prepared state $\rho$ (although we do not have the
ability to prepare $\rho$ noise-free). The sampling overhead factor will be
denoted by $C$.

#### Zero-noise extrapolation

Zero-noise extrapolation (ZNE) [li2017efficient, temme2017error] operates
under the assumption that the accumulation of errors is sufficiently well-
behaved. Suppose that the noise can be effectively modeled in terms of a
“circuit fault rate” $\lambda\geq 0$, such that the noisy quantum circuit
prepares the state $\rho(\lambda)$. The ideal, noiseless state is
$\rho\equiv\rho(0)$, although at best we only have access to
$\rho(\lambda_{\min})$ for some $\lambda_{\min}>0$. Then, rather than reduce
the fault rate, one can consider _boosting_ the error rate within a sequence
$\lambda_{\min}\leq\lambda_{1}<\lambda_{2}<\cdots$. This might be achieved,
for example, by inserting “noisy identities” within the circuit (e.g.,
$G^{\dagger}G=\mathbb{I}$ for any gate $G$) [giurgica2020digital].

The expectation value $\operatorname{tr}(O\rho(\lambda))$ can then be modeled
as a function $f(\lambda;\bm{p})$ with some free model parameters
$\bm{p}\in\mathbb{R}^{M}$ to be fit. If $\lambda$ is small, one can make the
following approximation:

$\operatorname{tr}(O\rho(\lambda))\approx
f(\lambda;\bm{p})=\sum_{\ell=0}^{M-1}p_{\ell}\frac{\lambda^{\ell}}{\ell!}.$
(148)

Supposing we evaluate the noisy device at $M$ different fault rates
$\lambda_{1},\ldots,\lambda_{M}$, we can find the optimal parameters
$\bm{p}^{\star}$ such that the model best fits the measured expectation
values: $\operatorname{tr}(O\rho(\lambda_{m}))\approx f(\lambda_{m};\bm{p})$.
Then we extrapolate to the zero-noise value by evaluating the model
$f(\lambda;\bm{p}^{\star})$ at $\lambda=0$. For example, a popular option is
Richardson extrapolation, for which the ZNE prediction is equal to the zeroth
component of the optimized parameters [giurgica2020digital]:

$\operatorname{tr}(O\rho(0))\approx
p_{0}^{\star}=\sum_{m=1}^{M}\operatorname{tr}(O\rho(\lambda_{m}))\prod_{k\neq
m}\frac{\lambda_{k}}{\lambda_{k}-\lambda_{m}}.$ (149)

Assuming that the noise behaves under this circuit-fault-rate model, the bias
in this estimate is of order $\lambda^{M}$. The sampling overhead scales by a
factor of

$C\sim\mathopen{}\left(\sum_{m=1}^{M}\prod_{k\neq
m}\frac{\lambda_{k}}{\lambda_{k}-\lambda_{m}}\right)^{2}\mathclose{},$ (150)

which can be shown from a variance analysis of $p_{0}^{\star}$.

#### Probabilistic error cancellation

Probabilistic error cancellation (PEC) [temme2017error] is an approach using
Monte Carlo sampling techniques to, on average, cancel the effects of noisy
gates. The general theory can be described as follows. Let $U=U_{L}\cdots
U_{1}$ be the ideal quantum circuit desired, decomposed into a sequence of
elementary gates $U_{j}$. The channel representation of a unitary $U$ will be
denoted by $\mathcal{U}$. In the experiment, we only have access to noisy
physical gate operations $\\{\widetilde{\mathcal{G}}_{k}\\}_{k}$. Suppose each
ideal gate $\mathcal{U}_{j}$ can be decomposed into a linear combination over
the noisy gate set,

$\mathcal{U}_{j}=\sum_{k}\gamma_{k}^{(j)}\widetilde{\mathcal{G}}_{k},\quad\gamma_{k}^{(j)}\in\mathbb{R}.$
(151)

This is called a _quasiprobability decomposition_ of $\mathcal{U}_{j}$, since
it can be implemented in expectation by treating the coefficients $\gamma_{k}$
as a quasiprobability distribution. That is, for each gate $j\in[L]$, we
sample $k$ from the probability distribution
$|\gamma_{k}^{(j)}|/\|\bm{\gamma}^{(j)}\|_{1}$, implement the corresponding
$\widetilde{\mathcal{G}}_{k}$, and store the sign of $\gamma_{k}^{(j)}$. Then
after performing all the experiments, when averaging over the measurement data
we reweight this particular sample by a factor of
$\|\bm{\gamma}^{(j)}\|_{1}\operatorname{sign}(\gamma_{k}^{(j)})$.

As a simple demonstrative example of how this works, suppose the noisy
operation $\widetilde{\mathcal{U}}$ is the ideal unitary $\mathcal{U}$,
followed by a depolarizing channel of strength $\lambda$. The form of this
channel can be written as
$\widetilde{\mathcal{U}}=(1-\lambda)\mathcal{U}+\lambda\mathcal{E}$, where
$\mathcal{E}(\rho)\coloneqq\operatorname{tr}(\rho)\mathbb{I}/2^{n}$ sends
every input to the maximally mixed state. Rewriting this expression reveals

$\mathcal{U}=\frac{1}{1-\lambda}\widetilde{\mathcal{U}}-\frac{\lambda}{1-\lambda}\mathcal{E}.$
(152)

The channel $\mathcal{E}$ must be further written in terms of the physical
gate set. For example, we can use its Kraus decomposition in terms of all
Pauli operators:

$\mathcal{E}=\frac{1}{4^{n}}\sum_{P\in\\{\mathbb{I},X,Y,Z\\}^{\otimes
n}}\mathcal{P}.$ (153)

Assuming we can perform perfect Pauli gates,323232This is typically a
reasonable assumption, because in most platforms single-qubit gates have
infidelities at least an order of magnitude lower than those of two-qubit
gates. this gives a physically implementable quasiprobability decomposition as
desired. With probability $1/(1+\lambda)$, we implement the noisy unitary
$\widetilde{\mathcal{U}}$ as usual. But with probability
$\lambda/(1+\lambda)$, we implement a random Pauli matrix on each qubit
instead. Then after a sufficient number of rounds of this randomized protocol,
each measurement outcome wherein we performed $\widetilde{\mathcal{U}}$ is
scaled by $(1+\lambda)/(1-\lambda)$, while the outcomes associated with the
random Pauli gates is scaled by $-(1+\lambda)/(1-\lambda)$ instead.

Returning to the general framework with $L$ gates, the overall rescaling
factor is a product of each gate’s quasiprobability $1$-norm, and in a
straightforward manner the variance of estimation is augmented by the square
of that factor:

$C\sim\prod_{j=1}^{L}\|\bm{\gamma}^{(j)}\|_{1}^{2}.$ (154)

An alternative approach to PEC is to assume that each noisy gate takes the
form $\widetilde{\mathcal{U}}_{j}=\mathcal{E}_{j}\mathcal{U}_{j}$. After each
gate, we would like to apply the channel inverse $\mathcal{E}_{j}^{-1}$.
However, this inverse is typically not a physically implementable map (usually
failing to be completely positive). Instead, we find a quasiprobability
decomposition for
$\mathcal{E}_{j}^{-1}=\sum_{k}\gamma_{k}^{(j)}\widetilde{\mathcal{G}}_{k}$ and
effectively implement it using the Monte Carlo technique outlined in the
previous approach. The sampling overhead for this approach is of the same
form, although the $\gamma_{k}^{(j)}$’s will be different and therefore might
be advantageous depending on the details of the noise.

Essential to either of these PEC methods is an accurate characterization of
the error channels, as otherwise one does not know how to define the correct
quasiprobability decomposition. While in principle this can be achieved by
careful calibration of the quantum device, in practice this is both costly and
might not accurately capture error effects if they do not appear in
calibration (e.g., crosstalk or time-dependent phenomena). Therefore, a number
of learning-based approaches have been developed, wherein effective models for
the quasiprobability decomposition are fit to via training data
[strikis2021learning, van2023probabilistic]. For example, the training data
might consist of circuits which are structurally similar to the desired
circuit of interest, but can be efficiently validated by classical simulation.

#### Training error-mitigation models

Ideas from machine learning can be imported into QEM more broadly as well.
Here we will describe two similar methods that train a simple regression model
using noisy experimental data from classically simulable quantum circuits:
Clifford data regression (CDR) [czarnik2021error] and training with fermionic
linear optics (TFLO) [montanaro2021error]. As their names imply, they train
models from (nearly) Clifford and free-fermion circuits, respectively, both of
which can be simulated efficiently with classical computation. The structure
of the circuits, however, are chosen to be close to that of the target
circuits (which are ideally not classically simulable), such that we expect
the generalization error of this model to be small.

Both ideas use a linear regression model, based on the fact that many forms of
decohering errors corrupt the quantum state as
$\rho\mapsto\widetilde{\rho}=(1-\lambda)\rho+\lambda\sigma$, where $\sigma$ is
some other state and $0\leq\lambda\leq 1$. The quantum state itself is
produced by some PQC with parameters $\bm{\theta}\in\Theta$, so we write the
ideal states as $\rho(\bm{\theta})$. Assuming that neither the overall fault
rate $\lambda$ nor the erroneous state $\sigma$ change when the parameters are
varied, the noisy expectation value can be modeled as

$\operatorname{tr}(O\widetilde{\rho}(\bm{\theta}))=(1-\lambda)\operatorname{tr}(O\rho(\bm{\theta}))+\lambda\operatorname{tr}(O\sigma).$
(155)

The noisy values $\operatorname{tr}(O\widetilde{\rho}(\bm{\theta}))$ are
measured from the quantum computer, while the noise-free values
$\operatorname{tr}(O\rho(\bm{\theta}))$ can be classically computed at special
values of $\bm{\theta}\in\mathcal{T}\subset\Theta$, corresponding to the set
of training circuits. Then $a=1/(1-\lambda)$ and
$b=-\lambda\operatorname{tr}(O\sigma)/(1-\lambda)$ are the unknown model
parameters that we optimize via least-squares minimization:

$(a^{\star},b^{\star})=\operatorname*{arg\,min}_{a,b}\sum_{\bm{\theta}\in\mathcal{T}}\mathopen{}\left|a\operatorname{tr}(O\widetilde{\rho}(\bm{\theta}))+b-\operatorname{tr}(O\rho(\bm{\theta}))\right|^{2}\mathclose{}.$
(156)

The linear-regression model is therefore

$f(\bm{\theta};a^{\star},b^{\star})=a^{\star}\operatorname{tr}(O\widetilde{\rho}(\bm{\theta}))+b^{\star},$
(157)

which we hope will be close to the noise-free expectation value
$\operatorname{tr}(O\rho(\bm{\theta}))$ at any $\bm{\theta}\in\Theta$, not
just within the training set.

In the case of CDR, one sets the training parameters such that most of the
gates are Clifford and only at most $\mathcal{O}(\log n)$ gates are non-
Clifford; such nearly Clifford circuits can be classically simulated in
polynomial time [bravyi2016improved]. For example, the parametrized gate
$e^{-\mathrm{i}\theta P/2}$, where $P$ is any Pauli operator, is Clifford
whenever $\theta\in\\{-\pi/2,0,\pi/2,\pi\\}$.

For TFLO, the circuits are usually based on elementary one- and two-body
fermionic unitaries. For example,
$e^{-\mathrm{i}\theta(a_{p}^{\dagger}a_{q}+a_{q}^{\dagger}a_{p})/2}$ is one-
body, while $e^{-\mathrm{i}\theta a_{p}^{\dagger}a_{p}a_{q}^{\dagger}a_{q}}$
is two-body. Then, given some circuit structure of such gates, we simply
remove all the two-body gates333333Or set their angles to be a full $2\pi$
rotation such that they implement a noisy identity. such that the resultant
circuit is purely one-body, or a fermionic linear-optical network. Such
circuits can be efficiently classically simulated [terhal2002classical] (see
also Chapter Theory of Fermions).

Unfortunately, the sampling overhead of these approaches is difficult to
analytically study. In part, this is because the assumption that the noise
obeys the linear regression model is virtually never observed in practice,
although it is a reasonable approximation. It is also difficult to rigorously
ascertain that the training data is representative of the full parameter
space, which may require, e.g., out-of-distribution techniques [caro2023out]
to properly analyze. However, numerical and experimental evidence suggests
that the resources required are reasonable to achieve substantial error
mitigation in practice [czarnik2021error, montanaro2021error,
stanisic2022observing].

#### Symmetry verification

Symmetry-based QEM is a conceptually simple approach that uses the ubiquitous
presence of symmetries in physical systems. If we know that the ideal state
$\rho$ should lie in a symmetry sector (eigenspace) of some symmetry operator
$S$, then we can detect and potentially mitigate against errors which leak the
state into the wrong sectors by measuring that symmetry. Note that this is the
core idea behind much of quantum error correction [lidar2013quantum]; for
example in stabilizer codes, the joint $+1$-eigenspace of a set of Pauli
operators defines a subspace of quantum states called a codespace, in which
the encoded, error-free quantum information should lie. Correctable errors
then take the state outside of that codespace, which we detect and correct
through syndrome measurements of the code.

In the case of QEM, we do not aim to correct the errors coherently, but rather
use classical postprocessing (and perhaps low-overhead detection circuitry) to
counteract the effects of symmetry-violating noise. Besides, it is too costly
to build a quantum code which can detect a large variety of errors; instead,
we typically restrict ourselves to whichever symmetries are inherently
provided by the physics of the system. The simplest method is postselection on
measurement outcomes, wherein we discard all the data that do not obey the
symmetry. This is broadly known as symmetry verification (SV)
[mcardle2019error, bonet2018low]. However, this is only possible, without
additional circuitry, when the symmetry operator commutes with the observable
$O$ (or more generally, its expansion terms). For example, suppose the system
obeys a parity symmetry, corresponding to the operator $S=Z^{\otimes n}$. This
means the ideal state only has support on either all even- or odd-Hamming-
weight basis states. Then only observables which are also diagonal in the $Z$
basis can be estimated from the data that measured $S$.

If we allow for additional circuitry, then it is possible to develop SV
techniques with broader scope. For example, one can introduce an ancilla qubit
and apply $n$ $\mathrm{CNOT}$ gates to it, where each $\mathrm{CNOT}$ is
controlled on a qubit from the system register. If the system did not
experience a parity-violating error, then the measurement outcome of the
ancilla is predictable (either $|0\rangle$ or $|1\rangle$, depending on the
correct parity sector). Otherwise, if such an error did occur, then the
ancilla will reveal this fact by returning the wrong output. The system
itself, meanwhile, is free to be measured in whichever basis one desires. The
sampling overhead of performing postselection in this manner is inversely
proportional to the “pass rate,” or the percentage of outcomes that pass the
symmetry check. Specifically, let $\widetilde{\rho}$ be the noisy state and
$\Pi$ the projector onto the correct symmetry sector. That is,

$C\sim\operatorname{tr}(\Pi\widetilde{\rho})^{-1}$ (158)

for this form of coherent SV. While this overhead is more favorable compared
to some other approaches (for example, the incoherent SV that we will describe
below), parity verification already requires an additional $n$ two-qubit
($\mathrm{CNOT}$) gates. The verification circuit for more complex symmetries,
for example particle number $\eta$, requires $n\lceil\log\eta\rceil$
controlled-phase gates [mcardle2019error].

An alternative to the coherent detection of errors is to perform a virtual
projection via classical postprocessing on the measurement record
[bonet2018low].343434This idea is closely related to techniques called quantum
subspace expansion [mcclean2017hybrid, mcclean2020decoding]. This is an
incoherent approach to SV. Let $\Pi$ be the projector as before; for example,
in our simple parity example, $\Pi=(\mathbb{I}+S)/2$ projects onto the even-
parity subspace. Using the fact that the postselected state is
$\rho^{\prime}=\Pi\widetilde{\rho}\Pi/\operatorname{tr}(\Pi\widetilde{\rho})$,
we can write expectation values for $\rho^{\prime}$ as

$\operatorname{tr}(O\rho^{\prime})=\frac{\operatorname{tr}(O^{\prime}\widetilde{\rho})}{\operatorname{tr}(\Pi\widetilde{\rho})},$
(159)

where $O^{\prime}=\Pi O\Pi=(O+SO+OS+SOS)/4$. Then rather than actually
postselect the measurements, we instead estimate the noisy expectation values
these four operators and $S$. Using Eq. 159, we can construct estimates for
$\operatorname{tr}(O\rho^{\prime})$ without needing the additional controlled
operations and ancilla. The cost of this virtual postselection is that the
sampling overhead now scales as

$C\sim\operatorname{tr}(\Pi\widetilde{\rho})^{-2},$ (160)

because the appearance of the pass rate in the denominator of Eq. 159
translates to an inverse-quadratically scaling factor in the variance
[huggins2021efficient, cai2021quantum].

This idea of postprocessing to effect postselection has been extended in a
number of ways, for example, recombining the expansion operators with
different weights to reduce bias [cai2021quantum] or using classical shadows
to efficiently estimate the multiple expansion operators with a single
protocol [jnane2023quantum]. However, these SV approaches suffer from only
being blind to symmetry-respecting errors. Thus projecting the noisy state
into the correct symmetry sector has no guarantee for the closeness of the
projected state to the ideal state. In Chapter Group-theoretic Error
Mitigation Enabled by Classical Shadows and Symmetries, we introduce a
technique called _symmetry-adjusted classical shadows_ [zhao2023group], which
combines classical shadows with symmetry information in a paradigm beyond that
of SV. Under simplifying noise assumptions, one can guarantee the correctness
of the error-mitigated estimates with this protocol for a wide range of noise
models. The key idea is that the twirling aspect of classical shadows
scrambles the noise uniformly throughout the quantum state, so that (1) a much
wider class of errors become scrambled into symmetry-violating form, and (2)
detectable errors from just the symmetry become valid diagnostics for how the
errors affect the rest of the state.

#### Virtual distillation

Virtual distillation (VD) [huggins2021virtual], also introduced
contemporaneously as error suppression by derangement (ESD)
[koczor2021exponential], takes advantage of the fact that the ideal state is
typically pure, $\rho=|\psi\rangle\\!\langle\psi|$. It furthermore assumes
that the dominant errors in real hardware are mostly incoherent, or can be
made as such [wallman2016noise]. Consider the spectral decomposition of a
noisy density matrix,

$\widetilde{\rho}=\sum_{i=0}^{2^{n}-1}p_{i}|\chi_{i}\rangle\\!\langle\chi_{i}|,$
(161)

where the labels are such that $p_{0}\geq\cdots\geq p_{2^{n}-1}$. Assuming
that $p_{0}>p_{1}$ holds strictly, then taking powers of the density matrix

$\rho^{(M)}=\frac{\widetilde{\rho}^{M}}{\operatorname{tr}(\widetilde{\rho}^{M})}$
(162)

converges to $\lim_{M\to\infty}\rho^{(M)}=|\chi_{0}\rangle\\!\langle\chi_{0}|$
exponentially quickly in $M$.

In order to access the state $\rho^{(M)}$, VD/ESD proposes to “virtually
distill” the dominant eigenvector by preparing $M$ copies of
$\widetilde{\rho}$ in parallel and performing entangled measurements across
the copies. Rather than coherently purify the state, it estimates properties
of $\rho^{(M)}$ by employing the identity

$\operatorname{tr}(O\rho^{M})=\operatorname{tr}(\overline{O}S_{M}\widetilde{\rho}^{\otimes
M}),$ (163)

where $S_{M}$ is the cyclic shift operator on the $M$ systems and

$\overline{O}=\frac{1}{M}\sum_{j=1}^{M}\mathbb{I}\otimes\cdots\otimes\mathbb{I}\otimes\underbrace{O}_{j\text{th
system}}\otimes\mathbb{I}\otimes\cdots\otimes\mathbb{I}$ (164)

is the target observable, symmetrized over the $M$ systems. Then the mitigated
expectation value can be calculated as

$\operatorname{tr}(O\rho^{(M)})=\frac{\operatorname{tr}(O\rho^{M})}{\operatorname{tr}(\rho^{M})}=\frac{\operatorname{tr}(\overline{O}S_{M}\widetilde{\rho}^{\otimes
M})}{\operatorname{tr}(S_{M}\widetilde{\rho}^{\otimes M})},$ (165)

which involves measuring the observables $\overline{O}S_{M}$ and $S_{M}$ on
the $M$-system state. Is it worth noting that this technique can be viewed as
another use of symmetry to perform QEM; in this case, the symmetry is
artificially introduced through the preparation of $M$ copies of the same
$n$-qubit state.

While a powerful technique for suppressing device errors, VD/ESD involves a
number of important considerations. First, note that the pure state (the
dominant eigenvector) that is being distilled is $|\chi_{0}\rangle$, which has
no guarantee to be the ideal state $|\psi\rangle$. As a trivial
counterexample, a purely coherent error will unitarily rotate $|\psi\rangle$
to some other pure state, against which VD/ESD has no hope of mitigating. In
general, this “coherent mismatch,” defined as the infidelity
$c=1-|\langle\psi|\chi_{0}\rangle|^{2}$, was studied in depth in
[koczor2021dominant]. One important result from that work is that, while the
ultimate noise floor is bounded as
$\Delta=|\langle\psi|O|\psi\rangle-\langle\chi_{0}|O|\chi_{0}\rangle|\leq
2\sqrt{c}\|O\|_{\infty}$ for arbitrary states, a quadratically smaller bound
of $\Delta\leq 2c\|O\|_{\infty}$ holds whenever $|\psi\rangle$ is an
eigenstate of $O$, which is the setting of most VQAs.

Another challenge is the measurement of the highly global observables
$\overline{O}S_{M}$ and $S_{M}$. The presence of the cyclic shift operator
particularly complicates matters. For $M=2$ and $O$ being a single-qubit
observable, circuits were given in [huggins2021virtual] that require only a
single layer of transversal two-qubit gates between identically labeled qubits
across the copies. However, already for $M=3$ the measurement circuit needed
to be numerically approximated and optimized by [huggins2021virtual]. For
larger $M$, ancilla qubits are needed to perform Hadamard tests for the
measurement [huggins2021virtual, koczor2021exponential]. Furthermore, if $O$
is not a single-qubit observable then even more circuitry and ancilla are
required. The use of classical shadows for VD/ESD has been considered to
circumvent these challenges [seif2023shadow, hu2022logical], requiring only a
single copy and using the shadows to predict the nonlinear functions of the
state. However, such approaches feature an exponentially large sample
complexity upfront.

For $M=2$, the sampling overhead of the standard VD/ESD technique vanishes
quadratically with the impurity of the noisy state:

$C\sim\operatorname{tr}(\widetilde{\rho}^{2})^{-2}.$ (166)

One may observe that this behavior is similar to that of SV via
postprocessing, as a consequence of the “virtual” (incoherent) aspect. For
general $M$, $\operatorname{tr}(\widetilde{\rho}^{M})$ is exponentially small
in $M$ (unless $\widetilde{\rho}$ is pure) and so the sampling overhead will
grow exponentially in $M$.

#### Verified phase estimation

The last example we review here is verified phase estimation (VPE)
[obrien2021error], which falls broadly under a category of techniques known as
echo verification [cai2021resource, huo2022dual]. In a sense, VPE is similar
to VD/ESD in that it uses two copies of a(n ideally pure) state
$\rho=|\psi\rangle\\!\langle\psi|$ to perform error mitigation; however,
whereas in VD those copies are spatially separated, in VPE they are
_temporally_ separated. The basic premise is that, after preparing
$|\psi\rangle=U|0^{n}\rangle$ from some state-preparation circuit $U$,
applying the inverse $U^{\dagger}$ (hence the “echo”) and measuring the system
should always return a deterministic output, $|0^{n}\rangle$, in the absence
of noise. Of course, in the presence of errors corrupting both $U$ and
$U^{\dagger}$, the final noisy state will support some nontrivial distribution
in the computational basis. Thus we can use this measurement as a verification
of whether or not any errors have occurred—if any output other than
$|0^{n}\rangle$ is seen, we throw that run out (technically, we do not
postselect but rather assign such unverified runs an estimate of $0$ while
retaining that sample in our statistics). It can be shown that this is
equivalent to measuring the state
$(\widetilde{\rho}\rho+\rho\widetilde{\rho})/(2\operatorname{tr}(\widetilde{\rho}\rho))$
[huo2022dual], hence the similarity to second-order VD/ESD.

Because we are measuring the system qubits to perform the echo verification,
we need an alternative way to extract the information about observables from
the state. This is achieved using ideas from phase estimation
[somma2019quantum]. First, decompose the target observable $O$ into a linear
combination of terms $H_{s}$,

$O=\sum_{s=1}^{L}h_{s}H_{s}.$ (167)

such that the unitary $V_{s}(t)=e^{-\mathrm{i}H_{s}t}$ can be implemented for
any $t$ in a NISQ-friendly manner. For example, a local Pauli decomposition
satisfies this condition. Then for each $s\in[L]$, we perform a single-ancilla
phase estimation experiment with $V_{s}(t)$: apply the unitary c-$V_{s}(t)$ on
$|\psi\rangle$, controlled on an ancilla qubit prepared in
$|+\rangle=(|0\rangle+|1\rangle)/\sqrt{2}$. This yields the $(n+1)$-qubit
state

$|\Psi(t)\rangle=\sum_{j}a_{s,j}\mathopen{}\left(\frac{|0\rangle+e^{-\mathrm{i}E_{s,j}t}|1\rangle}{\sqrt{2}}\right)_{\mathrm{anc}}\mathclose{}\otimes|E_{s,j}\rangle_{\mathrm{sys}},$
(168)

where $H_{s}=\sum_{j}E_{s,j}|E_{s,j}\rangle\\!\langle E_{s,j}|$ is the
spectral decomposition of the term and $a_{s,j}=\langle E_{s,j}|\psi\rangle$.
The state of the control qubit is then

$\rho_{\mathrm{anc}}(t)=\operatorname{tr}_{\mathrm{sys}}(|\Psi(t)\rangle\\!\langle\Psi(t)|)=\frac{1}{2}\begin{pmatrix}1&g_{s}(t)\\\
g_{s}^{*}(t)&1\end{pmatrix},$ (169)

where

$g_{s}(t)=\sum_{j}|a_{s,j}|^{2}e^{-\mathrm{i}E_{s,j}t}$ (170)

is the phase function containing information about the observable $H_{s}$
within $|\psi\rangle$. Note that performing the inverse circuit $U^{\dagger}$
on the system does not change the state of the ancilla qubit. Thus we can both
perform verification and learn properties from the system by probing the
ancilla.

Specifically, one can learn this function by implementing c-$V_{s}(t_{k})$ at
various time points $t_{1},t_{2},\ldots$ and measuring the control qubit in
the $X$ and $Y$ bases. This reveals the signal because we can express
$g_{s}(t)=\operatorname{tr}(X\rho_{\mathrm{anc}}(t))+\mathrm{i}\operatorname{tr}(Y\rho_{\mathrm{anc}}(t))$.
Classical signal-analysis methods can then be employed to extract the
amplitudes $|a_{s,j}|^{2}$. Finally, the expectation value can be estimated
from these amplitudes using
$\langle\psi|H_{s}|\psi\rangle=\sum_{j}|a_{s,j}|^{2}E_{s,j}$. This requires
knowledge of the eigenvalues of $H_{s}$, which is generally available because
they are our basis operators for the decomposition of $O$; for example, the
spectrum of Pauli operators is simply $\\{\pm 1\\}$. For such bounded
operators, the number of time points required to resolve the signal is also
bounded, so the number of shots required is the usual
$\mathcal{O}(\epsilon^{-2})$ incurred from sampling the quantum circuits.

Repeating this for each $s=1,\ldots,L$ yields $\epsilon$-accurate estimates
for each $\langle\psi|H_{s}|\psi\rangle$, and so by linearity we have an
$\mathcal{O}(\|\bm{h}\|_{1}\epsilon)$-accurate estimate for
$\langle\psi|O|\psi\rangle$. Furthermore, a variance analysis reveals that the
overhead due to noise is

$C\sim p^{-1}_{\mathrm{NE}},$ (171)

where $p_{\mathrm{NE}}$ is the probability of no error occurring in the entire
quantum circuit (recall that this is $U^{\dagger}(\text{c-}V_{s}(t))U$, which
has essentially double the size of $U$ alone). Thus to obtain an
$\epsilon$-accurate estimate for $\langle\psi|O|\psi\rangle$, VPE naively
requires taking
$\mathcal{O}(L\|\bm{h}\|_{1}^{2}\epsilon^{-2}p_{\mathrm{NE}}^{-1})$ samples,
although this can be ameliorated, e.g., by parallelizing over groups of
commuting terms (although this requires a separate control qubit for each such
term). After paying these costs, however, VPE is notable in that it has the
potential to mitigate errors down to order of $(1-p_{\mathrm{NE}})^{2}$ or
even $(1-p_{\mathrm{NE}})^{3}$, whereas other QEM techniques typically only
suppress errors to first order. (For comparison, VD/ESD suppresses errors to
$M$th order.)

#### Comment on sampling overhead

As we have seen from these examples, there is always some overhead $C$ in the
number of samples required to resolve the same accuracy in error-mitigated
estimates. Indeed, this is a necessary feature of QEM [takagi2022fundamental,
takagi2022universal, tsubouchi2022universal, quek2022exponentially].
Furthermore, the accumulation of noise generally implies that this factor will
scale exponentially in the circuit size. As a simple example, let $D$ be the
depth of the circuit and consider the VPE method. Assuming a simple
depolarizing noise model wherein each qubit has error probability
$p=\mathcal{O}(1)$ at each gate, then $p_{\mathrm{NE}}=(1-p)^{nD}$, hence
generically [obrien2021error]

$C\sim\mathopen{}\left(\frac{1}{1-p}\right)^{nD}\mathclose{}=2^{\operatorname{poly}(n)}.$
(172)

Qualitatively speaking, the overheads for other QEM methods behave similarly
[cai2022quantum].

## Fermionic Partial Tomography via Classical Shadows

### Preface

This chapter is based on [zhao2021fermionic], coauthored by the author of this
dissertation, Nicholas C. Rubin, and Akimasa Miyake. New material has been
added to this chapter in Theorem 0.4, proven in Information-theoretic lower
bounds on predicting local fermionic observables, which establishes an
information-theoretic lower bound matching the performance of the protocol
constructed here.

### Introduction

One of the most promising applications of quantum computation is the study of
strongly correlated systems such as interacting fermions. While quantum
algorithms such as phase estimation [nielsen2010quantum, kitaev2002classical]
allow for directly computing important quantities such as ground-state
energies with quantum speedup [abrams1999quantum, somma2002simulating,
aspuru2005simulated], current hardware limitations [preskill2018quantum] have
directed much attention toward variational methods. Of note is the variational
quantum eigensolver (VQE) [peruzzo2014variational, mcclean2016theory], where
short-depth quantum circuits are repeatedly executed in order to estimate
observable expectation values.

Initial bounds on the number of these circuit repetitions associated with
fermionic two-body Hamiltonians were prohibitively high [wecker2015progress],
spurring on much recent work addressing this problem. We roughly classify
these strategies into two categories: those that specifically target energy
estimates [mcclean2014exploiting, mcclean2016theory, kandala2017hardware,
babbush2018low, rubin2018application, izmaylov2019revising,
izmaylov2019unitary, huggins2021efficient, crawford2019efficient,
zhao2020measurement, torlai2020precise, arrasmith2020operator,
paini2019approximate_superseded, hadfield2022measurements, yen2020cartan,
gonthier2022measurements, huang2021efficient, garcia2021learning,
hillmich2021decision, hadfield2021adaptive, wu2023overlapped], referred to as
Hamiltonian averaging, and more general techniques that can learn the $k$-body
reduced density matrices ($k$-RDMs) of a quantum state [aaronson2020shadow,
aaronson2018online, aaronson2019gentle, yu2019quantum, yu2020sample,
verteletskyi2020measurement, jena2019pauli, yen2020measuring, gokhale2019on3,
cotler2020quantum, bonet2020nearly, hamamura2020efficient, garcia2020pairwise,
jiang2020optimal, evans2019scalable, huang2020predicting, smart2020lowering,
tilly2021reduced]. (Not all works fit neatly into this dichotomy, e.g., Refs.
[harrow2019low, wang2019accelerated, kubler2020adaptive, sweke2020stochastic,
van2020measurement, wang2021minimizing].) Hamiltonian averaging is ultimately
interested in a single observable, allowing for heavy exploitation in its
structure. In contrast, reconstructing an RDM requires estimating all the
observables that parametrize it.

Though generally more expensive than Hamiltonian averaging, calculating the
$k$-RDM allows one to determine the expectation value of any $k$-body
observable [coleman1980reduced]. For example, the electronic energy of
chemical systems is a linear functional of the 2-RDM, while in condensed-
matter systems, effective models for electrons can require knowledge of the
3-RDM [tsuneyuki2008transcorrelated, PhysRevB.87.245129]. Beyond the energy,
other important physical properties include pair-correlation functions and
various order parameters [mazziotti2012two, jensen2017introduction]. The 2-RDM
is also required for a host of error-mitigation techniques for near-term
quantum algorithms [mcclean2017hybrid, rubin2018application,
takeshita2020increasing], which have been experimentally demonstrated to be
crucial in obtaining accurate results [colless2018computation,
sagastizabal2019experimental, mccaskey2019quantum, arute2020hartree].
Additionally, promising extensions to VQE such as adaptive ansatz construction
[grimsley2019adaptive, ryabinkin2020iterative, tang2019qubit,
wang2020resource] and multireference- and excited-state calculations
[mcclean2017hybrid, parrish2019quantum, takeshita2020increasing,
huggins2020non, stair2020multireference, urbanek2020chemistry] can require up
to the 4-RDM.

Motivated by these considerations, in this work we focus on partial tomography
for fermionic RDMs. While numerous works have demonstrated essentially optimal
sample complexity for estimating qubit RDMs [cotler2020quantum,
bonet2020nearly, jiang2020optimal, evans2019scalable, huang2020predicting],
such approaches necessarily underperform in the fermionic setting. Recognizing
this fundamental distinction, Bonet-Monroig _et al._ [bonet2020nearly]
developed a measurement scheme that achieves optimal scaling for
fermions.353535Ref. [jiang2020optimal] also claims an optimally scaling
protocol (which requires a specific fermion-to-qubit mapping and doubling the
number of qubits); unfortunately, the claim does not hold due to an error in
their conversion of a derived “attenuation factor” to sample complexity.
However, the construction is not readily generalizable for $k>2$, in part due
to its highly complicated design.

In this chapter, we propose a randomized scheme that is free from these
obstacles. It is based on the theory of classical shadows
[huang2020predicting]: a protocol of randomly distributed measurements from
which one acquires a partial classical representation of an unknown quantum
state (its “shadow”). Classical shadows are sufficient for learning a limited
collection of observables, making this framework ideal for partial state
tomography. Our key results identify efficient choices for the ensemble of
random measurements, suitable for the structure of fermionic RDMs.

### Fermionic RDMs

Consider a fixed-particle state $\rho$ represented in second quantization on
$n$ fermion modes. The $k$-RDM of $\rho$, obtained by tracing out all but $k$
particles, is typically represented as a $2k$-index tensor,

${}^{k}D_{q_{1}\cdots q_{k}}^{p_{1}\cdots
p_{k}}\coloneqq\operatorname{tr}(a_{p_{1}}^{\dagger}\cdots
a_{p_{k}}^{\dagger}a_{q_{k}}\cdots a_{q_{1}}\rho),$ (173)

where $a_{p}^{\dagger},a_{p}$ are fermionic creation and annihilation
operators, $p\in\\{0,\ldots,n-1\\}$. By linearity, these matrix elements may
be equivalently expressed using Majorana operators, starting with the
definitions

$\gamma_{2p}\coloneqq
a_{p}+a_{p}^{\dagger},\quad\gamma_{2p+1}\coloneqq-\mathrm{i}(a_{p}-a_{p}^{\dagger}).$
(174)

Then for each $2k$-combination $\bm{\mu}\equiv(\mu_{1},\ldots,\mu_{2k})$,
where $0\leq\mu_{1}<\cdots<\mu_{2k}\leq 2n-1$, we define a $2k$-degree
Majorana operator

$\Gamma_{\bm{\mu}}\coloneqq(-\mathrm{i})^{k}\,\gamma_{\mu_{1}}\cdots\gamma_{\mu_{2k}}.$
(175)

All unique $2k$-degree Majorana operators are indexed by the set of all
$2k$-combinations of $\\{0,\ldots,2n-1\\}$, which we shall denote by
$\mathcal{C}_{2n,2k}$. Because Majorana operators possess the same algebraic
properties as Pauli operators (Hermitian, self-inverse, and Hilbert–Schmidt
orthogonal), any fermion-to-qubit encoding maps between the two in a one-to-
one correspondence.

The commutativity structure inherited onto $\mathcal{C}_{2n,2k}$ constrains
the maximum number of mutually commuting (hence simultaneously measurable)
operators to be $\mathcal{O}(n^{k})$ [bonet2020nearly]. As there are
$\mathcal{O}(n^{2k})$ independent $k$-RDM elements, this implies an optimal
scaling of $\mathcal{O}(n^{k})$ measurement settings to account for all matrix
elements.

### Classical shadows and randomized measurements

We briefly review the framework of classical shadows introduced by Huang _et
al._ [huang2020predicting], upon which we build our fermionic extension and
prove sampling bounds. Let $\rho$ be an $n$-qubit state and
$\\{O_{1},\ldots,O_{L}\\}$ a set of $L$ traceless observables for which we
wish to learn
$\operatorname{tr}(O_{1}\rho),\ldots,\operatorname{tr}(O_{L}\rho)$. Classical
shadows require a simple measurement primitive: for each preparation of
$\rho$, apply the unitary map $\rho\mapsto U\rho U^{\dagger}$, where $U$ is
randomly drawn from some ensemble $\mathcal{U}$; then perform a projective
measurement in the computational basis, $\\{|z\rangle\mid
z\in\\{0,1\\}^{n}\\}$.

Suppose we have an efficient classical representation for inverting the
unitary map on postmeasurement states, yielding
$U^{\dagger}|z\rangle\\!\langle z|U$. Then the process of repeatedly applying
the measurement primitive and classically inverting the unitary may be viewed,
in expectation, as the quantum channel

$\mathcal{M}_{\mathcal{U}}(\rho)\coloneqq\operatorname*{\mathbb{E}}_{U\sim\mathcal{U},|z\rangle\sim
U\rho U^{\dagger}}\mathopen{}\left[U^{\dagger}|z\rangle\\!\langle
z|U\right]\mathclose{},$ (176)

where $|z\rangle\sim U\rho U^{\dagger}$ is defined by the usual probability
distribution from Born’s rule, $\Pr[|z\rangle\mid U\rho U^{\dagger}]=\langle
z|U\rho U^{\dagger}|z\rangle$. Informational completeness of $\mathcal{U}$
ensures that this channel is invertible, which allows us to define the
classical shadow

$\hat{\rho}_{U,z}\coloneqq\mathcal{M}_{\mathcal{U}}^{-1}\mathopen{}\left(U^{\dagger}|z\rangle\\!\langle
z|U\right)\mathclose{}$ (177)

associated with the particular copy of $\rho$ for which $U$ was applied and
$|z\rangle$ was obtained. Classical shadows form an unbiased estimator for
$\rho$, and so they can be used to estimate the expectation value of any
observable $O$:

$\operatorname*{\mathbb{E}}_{U\sim\mathcal{U},|z\rangle\sim U\rho
U^{\dagger}}\mathopen{}\left[\operatorname{tr}(O\hat{\rho}_{U,z})\right]\mathclose{}=\operatorname{tr}(O\rho).$
(178)

The number of repetitions $M$ required to obtain an accurate estimate for each
$\operatorname{tr}(O_{j}\rho)$ is controlled by the estimator’s variance,
which may be upper bounded by

$\max_{\text{states }\sigma}\mathbb{E}_{\begin{subarray}{c}U\sim\mathcal{U}\\\
|z\rangle\sim U\sigma U^{\dagger}\end{subarray}}\mathopen{}\left[\langle
z|U\mathcal{M}_{\mathcal{U}}^{-1}(O_{j})U^{\dagger}|z\rangle^{2}\right]\mathclose{}\eqqcolon\|O_{j}\|_{\mathcal{U}}^{2}.$
(179)

This quantity is referred to as the (squared) shadow norm. Then by median-of-
means estimation, one may show that

$M=\mathcal{O}\mathopen{}\left(\frac{\log L}{\varepsilon^{2}}\max_{1\leq j\leq
L}\|O_{j}\|_{\mathcal{U}}^{2}\right)\mathclose{}$ (180)

samples suffice to estimate all expectation values to within additive error
$\varepsilon$. To minimize Eq. (180) for a fixed collection of observables,
the only available freedom is in $\mathcal{U}$. One must therefore properly
choose the ensemble of unitaries, with respect to the target observables.

### Naive application to fermionic observables

A natural ensemble for near-term considerations is the group of single-qubit
Clifford gates, $\mathrm{Cl}(1)^{\otimes n}$ (i.e., Pauli measurements). For
an $\ell$-local Pauli observable $P$, Huang _et al._ [huang2020predicting]
showed that $\|P\|_{\mathrm{Cl}(1)^{\otimes n}}^{2}=3^{\ell}$, similar to the
results of prior approaches also based on Pauli measurements
[cotler2020quantum, bonet2020nearly, jiang2020optimal, evans2019scalable].
Although optimal for qubit $\ell$-RDMs, such strategies cannot achieve the
desired $\mathcal{O}(n^{k})$ scaling in the fermionic setting due to the
inherent nonlocality of fermion-to-qubit mappings. Indeed, assuming that the
$n$ fermion modes are encoded into $n$ qubits, the 1-degree Majorana operators
necessarily possess an average qubit locality of at least $\log_{3}(2n)$
[jiang2020optimal]. This implies that, under random Pauli measurements, the
squared shadow norm maximized over all $2k$-degree Majorana operators cannot
do better than $3^{2k\log_{3}(2n)}=4^{k}n^{2k}$. In fact, for commonly used
mappings such as the Jordan–Wigner [jordanwigner] or Bravyi–Kitaev
[bravyi2002fermionic, seeley2012bravyi, tranter2015bravyi,
havlivcek2017operator] transformations, the scalings are poorer ($3^{n}$ and
${\sim}\,9^{k}n^{3.2k}$, respectively).

### Randomized measurements with fermionic Gaussian unitaries

To obtain optimal scaling in the shadow norm for fermionic observables, we
propose randomizing over a different ensemble: the group of fermionic Gaussian
Clifford unitaries. First, the group of fermionic Gaussian unitaries
$\mathrm{FGU}(n)$ comprises all unitaries of the form

$U(e^{A})\coloneqq\exp\mathopen{}\left(\frac{1}{4}\sum_{\mu,\nu=0}^{2n-1}A_{\mu\nu}\gamma_{\mu}\gamma_{\nu}\right)\mathclose{},$
(181)

where $A=-A^{\mathsf{T}}\in\mathbb{R}^{2n\times 2n}$. This condition implies
that $\mathrm{FGU}(n)$ is fully characterized by the Lie group
$\mathrm{SO}(2n)$ [sattinger1986lie]. In particular, the adjoint action

$U(Q)^{\dagger}\gamma_{\mu}U(Q)=\sum_{\nu=0}^{2n-1}Q_{\mu\nu}\gamma_{\nu}\quad\forall
Q\in\mathrm{SO}(2n)$ (182)

allows for efficient classical simulation of this group [knill2001fermionic,
terhal2002classical, bravyi2004lagrangian, divincenzo2005fermionic,
jozsa2008matchgates]. Second, the Clifford group $\mathrm{Cl}(n)$ is the set
of all unitary transformations that permute $n$-qubit Pauli operators among
themselves. It also admits an efficient classical representation
[gottesman1998heisenberg, aaronson2004improved].

Because Majorana operators are equivalent to Pauli operators, we may deduce
from Eq. (182) that a unitary that is both Gaussian and Clifford corresponds
to $Q$ being a signed permutation matrix. Note that this defines the full
group of Majorana swap circuits [bonet2020nearly]. As the signs are irrelevant
for our purpose, we simply consider the group of $2n\times 2n$ permutation
matrices with determinant 1, known as (the faithful matrix representation of)
the alternating group, $\mathrm{Alt}(2n)$.

Concretely, we set

$\mathcal{U}_{\mathrm{FGU}}\coloneqq\\{U(Q)\in\mathrm{FGU}(n)\mid
Q\in\mathrm{Alt}(2n)\\}.$ (183)

Given the context of fermionic tomography, the motivation for studying
$\mathrm{FGU}(n)$ is clear, as it preserves the degree of Majorana operators.
On the other hand, the restriction to the discrete Clifford elements is
valuable for practical considerations. As we show in Section Computations with
the fermionic Gaussian Clifford ensemble of the Supplemental Material (SM),
the permutational property of Clifford transformations necessarily implies
that $\mathcal{M}_{\mathrm{FGU}}$, as a linear map on the algebra of fermionic
observables, is diagonal in the Majorana-operator basis,

$\mathcal{M}_{\mathrm{FGU}}(\Gamma_{\bm{\mu}})=\lambda_{\bm{\mu}}\Gamma_{\bm{\mu}}\quad\forall\bm{\mu}\in\mathcal{C}_{2n,2k},$
(184)

with eigenvalues

$\lambda_{\bm{\mu}}=\mathopen{}\left.\binom{n}{k}\middle/\binom{2n}{2k}\right.\mathclose{}\equiv\lambda_{n,k}.$
(185)

In this diagonal form, the channel is readily invertible. Thus one may obtain
closed-form expressions for the classical shadows $\hat{\rho}_{Q,z}$, and,
importantly, their corresponding estimators for
$\operatorname{tr}(\Gamma_{\bm{\mu}}\rho)$:

$\operatorname{tr}(\Gamma_{\bm{\mu}}\hat{\rho}_{Q,z})=\lambda_{n,k}^{-1}\sum_{\bm{\nu}\in\mathcal{C}_{2n,2k}}\langle
z|\Gamma_{\bm{\nu}}|z\rangle\det[Q_{\bm{\nu},\bm{\mu}}].$ (186)

Here, $Q_{\bm{\nu},\bm{\mu}}$ denotes the submatrix of $Q$ formed from its
rows and columns indexed by $\bm{\nu}$ and $\bm{\mu}$, respectively
[chapman2018classical]. Because $Q$ is a permutation matrix, for each
$\bm{\mu}$ there is exactly one $\bm{\nu}^{\prime}$ such that
$\det[Q_{\bm{\nu}^{\prime},\bm{\mu}}]\neq 0$. Thus Eq. (186) is nonzero if and
only if that $\Gamma_{\bm{\nu}^{\prime}}$ is diagonal (i.e., maps to a
Pauli-$Z$ operator under a fermion-to-qubit transformation). In other words,
the Clifford operation $U(Q)$ sends $\Gamma_{\bm{\mu}}$ to
$\pm\Gamma_{\bm{\nu}^{\prime}}$, which can be estimated only if it is diagonal
in the computational basis.

From Eq. (179), the eigenvalues $\lambda_{n,k}^{-1}$ of the inverse channel
$\mathcal{M}^{-1}_{\mathrm{FGU}}$ determine the shadow norm. The sample
complexity of our approach then follows from Eq. (180). We summarize this
first key result with the following theorem.

######  Theorem 0.3.

Consider all $2k$-degree Majorana operators $\Gamma_{\bm{\mu}}$ on $n$
fermionic modes, labeled by $\bm{\mu}\in\mathcal{C}_{2n,2k}$. Under the
ensemble $\mathcal{U}_{\mathrm{FGU}}$ defined in Eq. (183), the shadow norm
satisfies

$\|\Gamma_{\bm{\mu}}\|_{\mathrm{FGU}}^{2}=\mathopen{}\left.\binom{2n}{2k}\middle/\binom{n}{k}\right.\mathclose{}\approx\binom{n}{k}\sqrt{\pi
k}$ (187)

for all $\bm{\mu}\in\mathcal{C}_{2n,2k}$. Thus the method of classical shadows
estimates the fermionic $k$-RDM of any state $\rho$, i.e.,
$\operatorname{tr}(\Gamma_{\bm{\mu}}\rho)\
\forall\bm{\mu}\in\bigcup_{j=1}^{k}\mathcal{C}_{2n,2j}$, to additive error
$\varepsilon$, given

$M=\mathcal{O}\mathopen{}\left[\binom{n}{k}\frac{k^{3/2}\log
n}{\varepsilon^{2}}\right]\mathclose{}$ (188)

copies of $\rho$. Additionally, there is no subgroup
$G\subset\mathrm{FGU}(n)\cap\mathrm{Cl}(n)$ for which
$\|\Gamma_{\bm{\mu}}\|_{G}<\|\Gamma_{\bm{\mu}}\|_{\mathrm{FGU}}$.

The proof is presented in the SM, Section Computations with the fermionic
Gaussian Clifford ensemble. Furthermore, noting from Eq. (186) that
$\mathopen{}\left|\operatorname{tr}(\Gamma_{\bm{\mu}}\hat{\rho}_{Q,z})\right|\mathclose{}\leq\lambda_{n,k}^{-1}$,
we also show in the SM that Bernstein’s inequality
[boucheron2013concentration] guarantees the above sample complexity via
standard sample-mean estimation, rather than requiring the median-of-means
technique proposed in the original work on classical shadows
[huang2020predicting].

This result has an intuitive conceptual interpretation. In the computational
basis, there are precisely $\binom{n}{k}$ diagonal Majorana operators within
$\mathcal{C}_{2n,2k}$, corresponding to the unique $k$-fold products of
occupation-number operators (e.g.,
$\prod_{j=1}^{k}a_{p_{j}}^{\dagger}a_{p_{j}}$) on $n$ modes. As a permutation
on $\mathcal{C}_{2n,2k}$, each element of $\mathcal{U}_{\mathrm{FGU}}$ defines
a different basis in which some other subset of $\binom{n}{k}$ operators are
diagonal. Then, one may expect to account for all
$|\mathcal{C}_{2n,2k}|=\binom{2n}{2k}$ Majorana operators by randomly
selecting on the order of $\binom{2n}{2k}/\binom{n}{k}$ such bases; Theorem
0.3 makes this claim rigorous.

We also establish a matching lower bound on the sample complexity for the task
of estimating multiple $2k$-degree Majorana-operator expectation values with
single-copy, local fermionic measurements. Such measurements encompass any
POVM that can be implemented with fermionic Gaussian operations, such as our
Gaussian Clifford ensemble. This result implies that, under single-copy
measurements, any protocol that can uniformly learn all $k$-body fermionic
observables necessarily requires the same number of samples as our
constructive protocol described above, up to constant factors.

######  Theorem 0.4 (Informal).

Any learning protocol based on single-copy, local-fermionic measurements that
can estimate the expectation values of any collection of $L\leq\binom{2n}{2k}$
$k$-body fermionic operators to additive error $\varepsilon$ requires at least

$M\geq\Omega\mathopen{}\left(\frac{\lambda_{n,k}^{-1}\log
L}{\varepsilon^{2}}\right)\mathclose{}$ (189)

copies of $\rho$.

The formal statement of the theorem and its proof are provided in the SM,
Information-theoretic lower bounds on predicting local fermionic observables.
The proof technique is a fermionic generalization of the local-qubit argument
presented in Huang _et al._ [huang2020predicting], which is based on
information-theoretic limits of communicating classical information through
quantum states. Our result is complementary to the lower bound proven by
Bonet-Monroig _et al._ [bonet2020nearly], which addresses nearly the same
task. Their bound differs from ours by a logarithmic factor, and their proof
method is also drastically different from ours. Also, they considered the
scenario of general Clifford measurements, whereas our result concerns local
fermionic (Gaussian) measurements. It interesting that both our and their
constructive protocols employ measurements that are at the intersection of
Gaussian and Clifford transformations, hence the two lower bounds nicely
complement each other. Finally, we note that our result applies to arbitrary
states $\rho$, so it does not apply to, e.g., the protocol of Low
[low2022classical] which targets only fixed-particle-number states.

Fermionic Gaussian circuits have a well-studied compilation scheme based on a
Givens-rotation decomposition [wecker2015solving, kivlichan2018quantum,
jiang2018quantum]. For a general element of $\mathcal{U}_{\mathrm{FGU}}$, we
require a circuit depth of at most $2n$ with respect to this decomposition
[jiang2018quantum]. Additionally, as pointed out in Ref. [bonet2020nearly],
Gaussian unitaries commute with the global parity operator
$\Gamma_{(0,\ldots,2n-1)}$, allowing for error mitigation via symmetry
verification [bonet2018low, mcardle2019error].

Such compilation schemes make use of a group homomorphism property,
$U(Q_{1})U(Q_{2})=U(Q_{1}Q_{2})$. Therefore, if the circuit preparing $\rho$
itself features fermionic Gaussian operations at the end, then we may further
compile the measurement unitary into the state-preparation circuit
[takeshita2020increasing]. In the case of indefinite particle number, this
concatenation is essentially free. However, rotations with particle-number
symmetry have depth at most $n$ [kivlichan2018quantum, jiang2018quantum], so
they must be embedded into the larger Gaussian unitary of depth $2n$. This
observation motivates us to explore classical shadows over the number-
conserving (NC) subgroup of $\mathrm{FGU}(n)$.

### Modification based on particle-number symmetry

Fermionic Gaussian unitaries that preserve particle number are naturally
parametrized by $\mathrm{U}(n)$. We express an element of this NC subgroup as

$U(e^{\kappa})\coloneqq\exp\mathopen{}\left(\sum_{p,q=0}^{n-1}\kappa_{pq}a_{p}^{\dagger}a_{q}\right)\mathclose{},$
(190)

where $\kappa=-\kappa^{\dagger}\in\mathbb{C}^{n\times n}$, hence
$e^{\kappa}\in\mathrm{U}(n)$. Because the particle-number symmetry manifests
as a global phase factor $e^{\operatorname{tr}\kappa/2}\in\mathrm{U}(1)$,
without loss of generality we may consider $\operatorname{tr}\kappa=0$, or
equivalently, $e^{\kappa}\in\mathrm{SU}(n)$. Such unitaries are also called
orbital-basis rotations, owing to their adjoint action,

$U(u)^{\dagger}a_{p}U(u)=\sum_{q=0}^{n-1}u_{pq}a_{q}\quad\forall
u\in\mathrm{SU}(n).$ (191)

This action on Majorana operators follows by linear extension.

Taking the intersection with the Clifford group requires that $u$ be an
$n\times n$ generalized permutation matrix, with nonzero elements taking
values in $\\{\pm 1,\pm i\\}$. This corresponds to the group of fermionic swap
circuits [bravyi2002fermionic, kivlichan2018quantum]. Again, the phase factors
on the matrix elements are irrelevant, so we shall restrict to
$u\in\mathrm{Alt}(n)$. By itself, this ensemble is insufficient to perform
tomography. To see this, consider an arbitrary reduced density operator
$A_{\bm{p}}^{\dagger}A_{\bm{q}}\coloneqq a_{p_{1}}^{\dagger}\cdots
a_{p_{k}}^{\dagger}a_{q_{k}}\cdots a_{q_{1}}$, where
$\bm{p},\bm{q}\in\mathcal{C}_{n,k}$. Such operators are diagonal in the
computational basis only if $\bm{p}=\bm{q}$. Informational completeness thus
requires that there exists some $U(u)$ that maps
$A_{\bm{p}}^{\dagger}A_{\bm{q}}$ to $A_{\bm{r}}^{\dagger}A_{\bm{r}}$, for some
$\bm{r}\in\mathcal{C}_{n,k}$. Because $u\in\mathrm{Alt}(n)$, conjugation by
$U(u)$ simply permutes $\bm{p}$ and $\bm{q}$ independently. However, as
permutations are bijective, it is not possible to permute both $\bm{p}$ and
$\bm{q}$ to the same $\bm{r}$ if $\bm{p}\neq\bm{q}$.

Therefore, this ensemble will necessarily require operations beyond either the
NC or Gaussian constraints. The simplest option for maintaining the low-depth
structure of the basis rotations is to append Pauli measurements at the end of
the circuit. Although the resulting circuit no longer preserves particle
number, this addition incurs only a single layer of single-qubit gates.
Specifically, we define the ensemble

$\mathcal{U}_{\mathrm{NC}}\coloneqq\\{V\circ U(u)\mid
V\in\mathrm{Cl}(1)^{\otimes n},\,u\in\mathrm{Alt}(n)\\}.$ (192)

By virtue of introducing the notion of “single-qubit” gates, this method is
dependent on the choice of fermion-to-qubit mapping. Let
$\mathrm{loc}(\Gamma_{\bm{\mu}})$ denote the qubit locality of
$\Gamma_{\bm{\mu}}$ under some chosen mapping. While Pauli measurements incur
a factor of $3^{\mathrm{loc}(\Gamma_{\bm{\mu}})}$ in the variance, the
randomization over fermionic swap circuits effectively averages this quantity
over all same-degree Majorana operators (rather than depending solely on the
most nonlocal operator). Formally, we find that the shadow norm here is

$\|\Gamma_{\bm{\mu}}\|_{\mathrm{NC}}^{2}=\operatorname*{\mathbb{E}}_{u\sim\mathrm{Alt}(n)}\mathopen{}\left[3^{-\mathrm{loc}[U(u)^{\dagger}\Gamma_{\bm{\mu}}U(u)]}\right]^{-1}\mathclose{}.$
(193)

Although this expression does not possess a closed form, the following theorem
provides a universal upper bound.

######  Theorem 0.5.

Under the ensemble $\mathcal{U}_{\mathrm{NC}}$ defined in Eq. (192), the
shadow norm obeys

$\max_{\bm{\mu}\in\mathcal{C}_{2n,2k}}\|\Gamma_{\bm{\mu}}\|_{\mathrm{NC}}^{2}\leq\mathopen{}\left.9^{k}\binom{n}{2k}\middle/\binom{n-k}{k}\right.\mathclose{}=\mathcal{O}(n^{k})$
(194)

for a fixed integer $k$ and for all fermion-to-qubit mappings. Thus the method
of classical shadows with $\mathcal{U}_{\mathrm{NC}}$ estimates the $k$-RDM to
additive error $\varepsilon$ with sample complexity

$M=\mathcal{O}\mathopen{}\left(\frac{n^{k}\log
n}{\varepsilon^{2}}\right)\mathclose{}.$ (195)

We provide derivations for the above results in the SM, Section Computations
with the number-conserving modification. Note that we have fixed $k$ as a
constant here, so the asymptotic notation may hide potentially large
prefactors depending on $k$. To understand such details, we turn to numerical
studies.

### Numerical calculations

Instead of drawing a new circuit for each repetition, here we employ a
simplification more amenable to practical implementation. Fixing some integer
$r\geq 1$, we generate a random collection
$\\{U^{(j)}\sim\mathcal{U}\\}_{j=1}^{K_{r}}$ of $K_{r}$ unitaries such that
all target observables are covered at least $r$ times. We say a Majorana
operator $\Gamma_{\bm{\mu}}$ is covered by the measurement unitary $U$ if
$U\Gamma_{\bm{\mu}}U^{\dagger}$ is diagonal in the computational basis.
Because the ensembles considered here consist of Gaussian and Clifford
unitaries, we can determine all covered operators efficiently. Additionally,
for the $\mathcal{U}_{\mathrm{NC}}$ calculations, the qubit mappings were
automated through OpenFermion [openfermion].

To achieve precision corresponding to $S=\mathcal{O}(1/\varepsilon^{2})$
samples per observable, one repeats each circuit $\lceil S/r\rceil$ times. The
total number of circuit repetitions for our randomized protocols is then
$\lceil S/r\rceil K_{r}$. For practical purposes, we fix $r=50$ in this work
(see Section Hyperparameter tuning of the SM for further details). To compare
against prior deterministic strategies, we compute $S\times C$ for each such
strategy, where $C$ is the number of sets of commuting observables constructed
by a given strategy.

Figure 1: (Left) Summary of the methods compared here, cataloging their
required circuit types and scalings in the number of measurement settings.
Because graph-based methods [jena2019pauli, yen2020measuring, gokhale2019on3,
hamamura2020efficient] are resource intensive, we employ sorted insertion
[crawford2019efficient] as a more tractable alternative. The Majorana clique
cover [bonet2020nearly], which employs the same class of fermionic Gaussian
Clifford circuits as our classical shadows (FGU) unitaries, possesses optimal
asymptotic scaling; however, it exhibits jumps at powers of 2 due to a divide-
and-conquer approach. Furthermore, the construction exists only for $k\leq 2$.
The measurement strategy using fermionic swap networks is a generalization of
the optimal 1-RDM strategy introduced in Ref. [arute2020hartree], which we
describe in Section Fermionic swap network bounds of the SM. (Right) Numerical
performances (log–log scale). Note that sorted insertion and the Majorana
clique cover are equivalent for $k=1$. Because our scheme uses randomization,
we include error bars of 1 standard deviation, averaged over 10 instances.
However, they are not visible at the scale of the plots, indicating the
consistency of our method.

For the comparisons presented in Fig. 1, we focus on the most competitive
prior strategies applicable to fermionic RDM tomography. Because the 1-RDM has
a relatively simple structure, optimal strategies are known [bonet2020nearly,
arute2020hartree], and so randomization underperforms for $k=1$. However, the
advantage of our $\mathcal{U}_{\mathrm{FGU}}$-based method becomes clear for
$k\geq 2$. When comparing against the Majorana clique cover, which features
asymptotically optimal $\mathcal{O}(n^{2})$ scaling for the 2-RDM
[bonet2020nearly], we find a roughly twofold factor improvement by our
approach.

For the $\mathcal{U}_{\mathrm{NC}}$ case, we observe a trade-off between
circuit size and measurement efficiency. As expected, the choice of fermion-
to-qubit mapping matters here; the Jordan–Wigner (JW) mapping performs worse
than Bravyi–Kitaev (BK), as the former possesses more qubit nonlocality.
Although more measurement settings are required compared to the
$\mathcal{U}_{\mathrm{FGU}}$ ensemble (e.g., a factor of ${\sim}\,\text{2--5}$
under BK, depending on $k$), each circuit itself requires only half the depth
of general fermionic Gaussian circuits. Notably, however,
$\mathcal{U}_{\mathrm{NC}}$ classical shadows for the 2-RDM under the BK
mapping is closely comparable to the Majorana clique cover.

### Conclusions

We have adapted the framework of classical shadows to the efficient tomography
of fermionic $k$-RDMs, applicable for all $k$. Numerical calculations
demonstrate that our approach consistently outperforms prior strategies using
measurement circuits of comparable sizes when $k\geq 2$, despite the
logarithmic factor in the sample complexity (a consequence of rigorously
bounding the worst-case probabilistic instances). The power of randomization
here lies in avoiding the $\mathsf{NP}$-hard problem of partitioning
observables into commuting cliques [verteletskyi2020measurement,
yen2020measuring, jena2019pauli, gokhale2019on3]. Instead, we show that a
highly overlapping cover of the observables suffices to perform partial
tomography efficiently, as a factor of $\mathcal{O}(1/\varepsilon^{2})$
repetitions is already required for this task.

An outlook for further applications is to adapt these ensembles, e.g., for
Hamiltonian averaging. As expected, our method is less efficient in this
context than those tailored for the task (see Section Hamiltonian averaging of
the SM for preliminary numerical calculations). Possible modifications may
include biasing the distribution of unitaries [hadfield2022measurements,
hillmich2021decision, hadfield2021adaptive, wu2023overlapped], or
derandomization techniques [huang2021efficient].

### Supplemental Material

### Additional notation

First, we define some notation and preliminary concepts not discussed in the
main text. For completeness, we generalize the definition of Majorana
operators to include odd-degrees:

$\Gamma_{\bm{\mu}}\coloneqq(-\mathrm{i})^{\binom{k}{2}}\gamma_{\mu_{1}}\cdots\gamma_{\mu_{k}},$
(196)

where $\bm{\mu}\in\mathcal{C}_{2n,k}$. Where required, we may define the empty
product ($k=0$) as the identity, i.e.,
$\Gamma_{\varnothing}\coloneqq\mathbb{I}$. Then for each $0\leq k\leq 2n$, we
define the $\mathbb{R}$-linear span of $k$-degree Majorana operators on $n$
fermion modes,

$\mathcal{A}_{k}^{(n)}\coloneqq\operatorname{span}\\{\Gamma_{\bm{\mu}}\mid\bm{\mu}\in\mathcal{C}_{2n,k}\\},$
(197)

which is isomorphic (as a vector space) to $\mathbb{R}^{\binom{2n}{k}}$. The
algebra of all even-degree fermionic observables shall be denoted by

$\mathcal{A}_{\mathrm{even}}^{(n)}\coloneqq\bigoplus_{k=0}^{n}\mathcal{A}_{2k}^{(n)}.$
(198)

Due to the parity superselection rule [streater2000pct], physical fermionic
operators lie in $\mathcal{A}_{\mathrm{even}}^{(n)}$; thus the notions of
informational completeness for fermionic tomography are understood with
respect to this algebra. To see this, suppose we interpret $|z\rangle$ as a
Fock basis state. Only products of occupation-number operators,
$\prod_{p}a_{p}^{\dagger}a_{p}$ (modulo anticommutation relations), have
nonzero diagonal matrix elements in the Fock basis. Thus if $U$ is a unitary
generated by a fermionic Hamiltonian (hence respecting parity supersymmetry),
only those $O\in\mathcal{A}_{\mathrm{even}}^{(n)}$ are able to satisfy
$\langle z|UOU^{\dagger}|z\rangle\neq 0$. This argument also holds if one
insists on viewing $|z\rangle$ as a qubit computational basis state, since
they are necessarily mapped to Fock basis states under a fermion-to-qubit
encoding [steudtner2018fermion].

The group of fermionic Gaussian unitaries $\mathrm{FGU}(n)$ is generated by
$\mathrm{i}\mathcal{A}_{2}^{(n)}$, and its adjoint action on any $k$-degree
Majorana operator straightforwardly generalizes as

$U(Q)^{\dagger}\Gamma_{\bm{\mu}}U(Q)=\sum_{\bm{\nu}\in\mathcal{C}_{2n,k}}\det[Q_{\bm{\mu},\bm{\nu}}]\Gamma_{\bm{\nu}},$
(199)

where $Q_{\bm{\mu},\bm{\nu}}$ denotes the submatrix formed by taking the rows
and columns of $Q$ indexed by $\bm{\mu}$ and $\bm{\nu}$, respectively
[chapman2018classical, Appendix A]. This defines an adjoint representation
$\mathrm{FGU}(n)\to\mathrm{SO}[\mathfrak{su}(2^{n})]\cong\mathrm{SO}(4^{n}-1)$,
understood in the sense that $\mathrm{FGU}(n)\subseteq\mathrm{SU}(2^{n})$. We
will rather be interested in the orthogonal representation $\Phi$ of
$\mathrm{SO}(2n)$ induced by this adjoint representation, through the group
homomorphism $U\colon\mathrm{SO}(2n)\to\mathrm{FGU}(n)$. That is, we define
$\Phi\colon\mathrm{SO}(2n)\to\mathrm{SO}(4^{n}-1)$ by the matrix elements

$[\Phi(Q)]_{\bm{\mu}\bm{\nu}}\coloneqq\det[Q_{\bm{\mu},\bm{\nu}}].$ (200)

By the Cauchy–Binet formula one may verify that this indeed defines an
orthogonal matrix [chapman2018classical]. Additionally, $\Phi$ inherits the
homomorphism property from $U$; hence it is an orthogonal representation.

Since determinants are defined only for square matrices, $\Phi$ possesses a
natural decomposition as $\Phi=\bigoplus_{k=0}^{2n}\phi_{k}$, where each
$\phi_{k}\colon\mathrm{SO}(2n)\to\mathrm{SO}(\mathcal{A}_{k}^{(n)})\cong\mathrm{SO}\mathopen{}\left[\binom{2n}{k}\right]\mathclose{}$
is defined just as in Eq. 200, restricted a particular $k$. These
subrepresentations will be the main focus of our analysis.

Finally, because of their relation to the Clifford group, we make frequent use
of permutations and their generalizations. We establish the relevant notation
here. Let $d,m>0$ be integers. We denote the symmetric group on $d$ symbols by
$\mathrm{Sym}(d)$, which is faithfully represented by the group of $d\times d$
permutation matrices. The alternating group $\mathrm{Alt}(d)$ is the subgroup
of all even parity permutations. The generalized symmetric group of cyclic
order $m$ over $d$ symbols is defined via the wreath product,
$\mathrm{Sym}(m,d)\coloneqq\mathbb{Z}_{m}\wr\mathrm{Sym}(d)\equiv\mathbb{Z}_{m}^{d}\rtimes\mathrm{Sym}(d)$.
Its faithful matrix representation is the group of $d\times d$ generalized
permutation matrices, wherein each nonzero matrix element can take on values
from the $m$th roots of unity. The determinant of such matrices is the sign of
the underlying permutation, multiplied by all $d$ nonzero elements. In
particular, we shall refer to $\mathrm{Sym}(2,d)$ as the group of signed
permutation matrices, and denote its $({+1})$-determinant subgroup by
$\mathrm{Sym}^{+}(2,d)$.

### Computations with the fermionic Gaussian Clifford ensemble

We now derive the main results leading to Theorem 0.3 of the main text. In The
classical shadows channel we find an expression for the channel
$\mathcal{M}_{\mathrm{FGU}}$, showing that the permutational property of the
Clifford group necessarily implies that the channel is diagonalized by the
basis of Majorana operators [Eq. 211]. Then in The shadow norm we compute the
corresponding eigenvalues of $\mathcal{M}_{\mathrm{FGU}}$, which are directly
related to the ensemble’s shadow norm (Lemma 0.9). In particular, Lemma 0.9
provides necessary and sufficient conditions for saturating the minimum value
of the shadow norm under Gaussian Clifford measurements. In Theorem 0.14 we
explicitly calculate this minimum value, and finally with Theorem 0.16 we
prove that both $\mathrm{FGU}(n)\cap\mathrm{Cl}(n)$ and its
$\mathrm{Alt}(2n)$-generated subgroup $\mathcal{U}_{\mathrm{FGU}}$ (the one
presented in the main text) satisfy the necessary and sufficient conditions.

After proving these main results, in The classical shadow estimator we
explicitly derive the expression for the classical shadows $\hat{\rho}_{Q,z}$.
For completeness, we show how to compute the shadow norm for an arbitrary
observable in Variance bounds for arbitrary observables. Finally, in
Performance guarantees without median-of-means estimation we show how the
boundedness of classical shadow estimators and Bernstein’s inequality allow
one to avoid requiring median-of-means estimation, as mentioned in the main
text.

#### The classical shadows channel

Our goal is to find an analytic expression for the channel
[huang2020predicting]
$\mathcal{M}_{\mathrm{FGU}}\colon\mathcal{A}_{\mathrm{even}}^{(n)}\to\mathcal{A}_{\mathrm{even}}^{(n)}$,

$\begin{split}\mathcal{M}_{\mathrm{FGU}}(O)&=\mathbb{E}_{U\sim\mathrm{FGU}(n)\cap\mathrm{Cl}(n)}\mathopen{}\left[\sum_{z\in\\{0,1\\}^{n}}\langle
z|UOU^{\dagger}|z\rangle U^{\dagger}|z\rangle\\!\langle
z|U\right]\mathclose{}\\\
&=\mathbb{E}_{Q\sim\mathrm{Sym}^{+}(2,2n)}\mathopen{}\left[\sum_{z\in\\{0,1\\}^{n}}\langle
z|U(Q)OU(Q)^{\dagger}|z\rangle U(Q)^{\dagger}|z\rangle\\!\langle
z|U(Q)\right]\mathclose{}.\end{split}$ (201)

Since this is a linear map, we need only to evaluate it on a basis of
$\mathcal{A}_{\mathrm{even}}^{(n)}$, the most natural choice being the
Majorana operators. Distinguishing the Majorana operators which are diagonal
(with respect to the computational basis) is highly important. For each $1\leq
k\leq n,$ we shall define the subset
$\mathcal{D}_{2n,2k}\subseteq\mathcal{C}_{2n,2k}$ of $2k$-combinations
corresponding to the diagonal $2k$-degree Majorana operators. Formally,

$\mathcal{D}_{2n,2k}\coloneqq\\{\bm{\mu}\in\mathcal{C}_{2n,2k}\mid\langle
z|\Gamma_{\bm{\mu}}|z\rangle\neq 0\ \forall z\in\\{0,1\\}^{n}\\}.$ (202)

Since there are $\binom{n}{k}$ independent $k$-fold products of occupation-
number operators, each set has cardinality
$|\mathcal{D}_{2n,2k}|=\binom{n}{k}$, so that

$\bigg{|}\bigcup_{k=1}^{n}\mathcal{D}_{2n,2k}\bigg{|}=\sum_{k=1}^{n}\binom{n}{k}=2^{n}-1,$
(203)

which indeed matches the maximal number of simultaneously commuting Pauli
operators [lawrence2002mutually] (e.g., all Pauli-$Z$ operators of locality
$1$ to $n$). For instance, under the Jordan–Wigner transformation, the
corresponding sets are

$\begin{split}\mathcal{D}_{2n,2}&\coloneqq\\{(p,p+1)\mid 0\leq p\leq
2n-2:p\text{ even}\\},\\\ \mathcal{D}_{2n,4}&\coloneqq\\{(p,p+1,q,q+1)\mid
0\leq p<q\leq 2n-2:p,q\text{ even}\\},\\\
\mathcal{D}_{2n,6}&\coloneqq\\{(p,p+1,q,q+1,r,r+1)\mid 0\leq p<q<r\leq
2n-2:p,q,r\text{ even}\\},\end{split}$ (204)

and so forth up to $\mathcal{D}_{2n,2n}$.

With this formalism, we can express basis states as

$|z\rangle\\!\langle
z|=\frac{1}{2^{n}}\mathopen{}\left(\mathbb{I}+\sum_{j=1}^{n}\sum_{\bm{\mu}\in\mathcal{D}_{2n,2j}}\langle
z|\Gamma_{\bm{\mu}}|z\rangle\Gamma_{\bm{\mu}}\right)\mathclose{},$ (205)

which undergo fermionic Gaussian evolution as

$\begin{split}U(Q)^{\dagger}|z\rangle\\!\langle
z|U(Q)&=\frac{1}{2^{n}}\mathopen{}\left(\mathbb{I}+\sum_{j=1}^{n}\sum_{\bm{\mu}\in\mathcal{D}_{2n,2j}}\langle
z|\Gamma_{\bm{\mu}}|z\rangle
U(Q)^{\dagger}\Gamma_{\bm{\mu}}U(Q)\right)\mathclose{}\\\
&=\frac{1}{2^{n}}\mathopen{}\left(\mathbb{I}+\sum_{j=1}^{n}\sum_{\bm{\mu}\in\mathcal{D}_{2n,2j}}\langle
z|\Gamma_{\bm{\mu}}|z\rangle\sum_{\bm{\nu}\in\mathcal{C}_{2n,2j}}\det\mathopen{}\left[Q_{\bm{\mu},\bm{\nu}}\right]\mathclose{}\Gamma_{\bm{\nu}}\right)\mathclose{}.\end{split}$
(206)

For the corresponding Born-rule probability factor, we use the fact that,
since $U$ is a group homomorphism, inverses are preserved
[$U(Q)^{\dagger}=U(Q^{\mathsf{T}})$], and so for any
$\bm{\tau}\in\mathcal{C}_{2n,2k}$ we have

$\begin{split}\langle
z|U(Q)\Gamma_{\bm{\tau}}U(Q)^{\dagger}|z\rangle&=\sum_{\bm{\sigma}\in\mathcal{C}_{2n,2k}}\det\mathopen{}\left[\mathopen{}\left(Q^{\mathsf{T}}\right)\mathclose{}_{\bm{\tau},\bm{\sigma}}\right]\mathclose{}\langle
z|\Gamma_{\bm{\sigma}}|z\rangle\\\
&=\sum_{\bm{\sigma}\in\mathcal{D}_{2n,2k}}\det\mathopen{}\left[\mathopen{}\left(Q^{\mathsf{T}}\right)\mathclose{}_{\bm{\tau},\bm{\sigma}}\right]\mathclose{}\langle
z|\Gamma_{\bm{\sigma}}|z\rangle.\end{split}$ (207)

Now recall some basic properties of Majorana (equiv. Pauli) operators. They
are traceless,

$\sum_{z\in\\{0,1\\}^{n}}\langle
z|\Gamma_{\bm{\sigma}}|z\rangle=\operatorname{tr}\Gamma_{\bm{\sigma}}=0,$
(208)

and more generally they are Hilbert–Schmidt (trace) orthogonal,

$\operatorname{tr}\mathopen{}\left(\Gamma_{\bm{\sigma}}\Gamma_{\bm{\mu}}\right)\mathclose{}=2^{n}\delta_{\bm{\sigma}\bm{\mu}}.$
(209)

In the case that $\Gamma_{\bm{\sigma}}$ and $\Gamma_{\bm{\mu}}$ are both
diagonal, we have that

$\begin{split}\sum_{z\in\\{0,1\\}^{n}}\langle
z|\Gamma_{\bm{\sigma}}|z\rangle\langle
z|\Gamma_{\bm{\mu}}|z\rangle&=\sum_{z\in\\{0,1\\}^{n}}\langle
z|\Gamma_{\bm{\sigma}}\Gamma_{\bm{\mu}}|z\rangle\\\
&=\operatorname{tr}\mathopen{}\left(\Gamma_{\bm{\sigma}}\Gamma_{\bm{\mu}}\right)\mathclose{}=2^{n}\delta_{\bm{\sigma}\bm{\mu}}.\end{split}$
(210)

With these relations in hand, we multiply Eqs. 206 and 207 and sum over all
$z\in\\{0,1\\}^{n}$ to obtain

$\begin{split}\sum_{z\in\\{0,1\\}^{n}}\langle
z|U(Q)\Gamma_{\bm{\tau}}U(Q)^{\dagger}|z\rangle
U(Q)^{\dagger}|z\rangle\\!\langle
z|U(Q)&=\frac{1}{2^{n}}\sum_{\bm{\sigma}\in\mathcal{D}_{2n,2k}}\det\mathopen{}\left[\mathopen{}\left(Q^{\mathsf{T}}\right)\mathclose{}_{\bm{\tau},\bm{\sigma}}\right]\mathclose{}\Bigg{(}\operatorname{tr}(\Gamma_{\bm{\sigma}})\mathbb{I}\\\
&\quad+\sum_{j=1}^{n}\sum_{\bm{\mu}\in\mathcal{D}_{2n,2j}}2^{n}\delta_{\bm{\sigma}\bm{\mu}}\sum_{\bm{\nu}\in\mathcal{C}_{2n,2j}}\det\mathopen{}\left[Q_{\bm{\mu},\bm{\nu}}\right]\mathclose{}\Gamma_{\bm{\nu}}\Bigg{)}\\\
&=\sum_{\begin{subarray}{c}\bm{\sigma}\in\mathcal{D}_{2n,2k}\\\
\bm{\nu}\in\mathcal{C}_{2n,2k}\end{subarray}}\det\mathopen{}\left[Q_{\bm{\sigma},\bm{\tau}}\right]\mathclose{}\det\mathopen{}\left[Q_{\bm{\sigma},\bm{\nu}}\right]\mathclose{}\Gamma_{\bm{\nu}}.\end{split}$
(211)

One may be tempted to use the Cauchy–Binet formula to evaluate the sum over
$\bm{\sigma}$; crucially, however, the sum is restricted to
$\mathcal{D}_{2n,2k}$, so the identity does not apply here. Instead, we will
first evaluate the sum over $\bm{\nu}$ by formalizing our notion of Gaussian
Clifford transformations as degree-preserving permutations of Majorana
operators. For generality, we state the following lemma with regards to any
generalized permutation matrix.

######  Lemma 0.6.

Let $Q\in\mathrm{Sym}(m,d)$ and fix $\bm{\tau}\in\mathcal{C}_{d,j}$, $1\leq
j\leq d$. Then there exists exactly one $\bm{\sigma}\in\mathcal{C}_{d,j}$ for
which $|{\det[Q_{\bm{\sigma},\bm{\tau}}]}|=1$; otherwise,
$\det[Q_{\bm{\sigma},\bm{\tau}}]=0$. In particular, we have

$\det[(Q^{\dagger})_{\bm{\tau},\bm{\sigma}}]\det[Q_{\bm{\sigma},\bm{\nu}}]=|{\det[Q_{\bm{\sigma},\bm{\tau}}]}|\,\delta_{\bm{\tau}\bm{\nu}}$
(212)

for all $\bm{\nu}\in\mathcal{C}_{d,j}$.

###### Proof.

By definition of permutation matrices, for each column $q$ of $Q$ there is
exactly one row $p$ for which $Q_{pq}\neq 0$. Furthermore, this row is unique
to each column. This property generalizes from matrix elements to
subdeterminants: for each set of columns indexed by $\bm{\tau}$, there is
exactly one unique set of rows $\bm{\sigma}$ for which
$Q_{\bm{\sigma},\bm{\tau}}$ has a nonzero element in each row. In other words,
$Q_{\bm{\sigma},\bm{\tau}}\in\mathrm{Sym}(m,j)$ and hence
$|{\det[Q_{\bm{\sigma},\bm{\tau}}]}|=1$ (recall that the determinant of these
matrices is the sign of the underlying permutation multipled by $m$th roots of
unity). Otherwise, for all $\bm{\sigma}^{\prime}\neq\bm{\sigma}$,
$Q_{\bm{\sigma}^{\prime},\bm{\tau}}$ possesses at least one row or column of
all zeros, and hence has determinant $0$.

Next, we show that submatrices behave under conjugate transposition as
$(Q_{\bm{\sigma},\bm{\tau}})^{\dagger}=(Q^{\dagger})_{\bm{\tau},\bm{\sigma}}$,
which can be seen by examining their matrix elements:

$\begin{split}[Q_{\bm{\sigma},\bm{\tau}}]_{pq}&=Q_{\sigma_{p}\tau_{q}}=[Q^{\dagger}]_{\tau_{q}\sigma_{p}}^{*}=[(Q^{\dagger})_{\bm{\tau},\bm{\sigma}}]_{qp}^{*}=[([Q^{\dagger}]_{\bm{\tau},\bm{\sigma}})^{\dagger}]_{pq}.\end{split}$
(213)

Since the determinant is invariant under transposition and preserves complex
conjugation,
$\det[(Q^{\dagger})_{\bm{\tau},\bm{\sigma}}]=\det[(Q_{\bm{\sigma},\bm{\tau}})]^{*}$.
This gives us

$\det[(Q^{\dagger})_{\bm{\tau},\bm{\sigma}}]\det[Q_{\bm{\sigma},\bm{\nu}}]=\det[(Q_{\bm{\sigma},\bm{\tau}})]^{*}\det[Q_{\bm{\sigma},\bm{\nu}}].$
(214)

But since $\bm{\tau}$ is the unique $j$-combination for which
$\det[(Q_{\bm{\sigma},\bm{\tau}})]\neq 0$ for a fixed $\bm{\sigma}$, the above
expression can only be nonzero when $\bm{\nu}=\bm{\tau}$. Equation (212) thus
follows. ∎

######  Remark 0.7.

This uniqueness property of nonzero subdeterminants implies that the image of
the orthogonal representation
$\Phi\colon\mathrm{SO}(2n)\to\mathrm{SO}(4^{n}-1)$, when restricted to a
subgroup $G\subseteq\mathrm{Sym}^{+}(2,2n)\subset\mathrm{SO}(2n)$, lies in
$\mathrm{Sym}^{+}(2,4^{n}-1)$. Its subrepresentations also satisfy
$\phi_{k}(G)\subseteq\mathrm{Sym}\mathopen{}\left[2,\binom{2n}{k}\right]\mathclose{}$.
From a physical perspective, these are all straightforward consequences of the
Clifford property imposed on our unitary ensemble.

Applying Lemma 0.6 to Eq. 211 reveals that the Majorana operators are in fact
the eigenbasis of $\mathcal{M}_{\mathrm{FGU}}$:

$\begin{split}\mathcal{M}_{\mathrm{FGU}}(\Gamma_{\bm{\tau}})&=\operatorname*{\mathbb{E}}_{Q\sim\mathrm{Sym}^{+}(2,2n)}\mathopen{}\left[\sum_{z\in\\{0,1\\}^{n}}\langle
z|U(Q)\Gamma_{\bm{\tau}}U(Q)^{\dagger}|z\rangle
U(Q)^{\dagger}|z\rangle\\!\langle z|U(Q)\right]\mathclose{}\\\
&=\operatorname*{\mathbb{E}}_{Q\sim\mathrm{Sym}^{+}(2,2n)}\mathopen{}\left[\sum_{\bm{\sigma}\in\mathcal{D}_{2n,2k}}|{\det[Q_{\bm{\sigma},\bm{\tau}}]}|\right]\mathclose{}\Gamma_{\bm{\tau}}.\end{split}$
(215)

#### The shadow norm

To evaluate the eigenvalues of Eq. 215, we invoke the theory of finite frames
[han2000frames, waldron2018introduction]. We begin with the definition of a
frame.

######  Definition 0.8.

Let $V$ be a Hilbert space with inner product $\langle\cdot,\cdot\rangle$. A
_frame_ is a sequence $\\{x_{j}\\}_{j}\subset V$ which satisfies

$\alpha\|v\|_{V}^{2}\leq\sum_{j}|\langle
x_{j},v\rangle|^{2}\leq\beta\|v\|_{V}^{2}\quad\forall v\in V,$ (216)

for some real constants $\alpha,\beta>0$ (called the _frame bounds_). Here,
$\|v\|_{V}\coloneqq\sqrt{\langle v,v\rangle}$. A frame is called _tight_ if
$\alpha=\beta$.

Recall that $\det[Q_{\bm{\sigma},\bm{\tau}}]$ defines the matrix elements of
the orthogonal representation
$\phi_{2k}\colon\mathrm{SO}(2n)\to\mathrm{SO}\mathopen{}\left[\binom{2n}{2k}\right]\mathclose{}$.
In order to demonstrate the optimality of our ensemble, we shall generalize
Eq. 215 to take the average over any subgroup
$G\subseteq\mathrm{Sym}^{+}(2,2n)$, in which case we consider the restricted
representations $\phi_{2k}\big{|}_{G}\colon
G\to\mathrm{SO}\mathopen{}\left[\binom{2n}{2k}\right]\mathclose{}$. (Whenever
the context is clear, we shall simply write $\phi_{2k}$.) To simplify
notation, we take $\mathbb{R}^{\binom{2n}{2k}}$ as our representation space,
spanned by the standard basis
$\\{e_{\bm{\mu}}\mid\bm{\mu}\in\mathcal{C}_{2n,2k}\\}$.

Consider the group orbit

$\phi_{2k}(G)e_{\bm{\sigma}}\coloneqq\\{\phi_{2k}(Q)e_{\bm{\sigma}}\mid Q\in
G\\}.$ (217)

A sufficient condition for the eigenvalues of $\mathcal{M}_{G}$ to be nonzero
(hence guaranteeing the existence of $\mathcal{M}_{G}^{-1}$) is that
$\phi_{2k}(G)e_{\bm{\sigma}}$ be a frame for some
$\bm{\sigma}\in\mathcal{D}_{2n,2k}$. Indeed, assuming the frame condition, we
have

$\begin{split}\operatorname*{\mathbb{E}}_{Q\sim
G}|{\det[Q_{\bm{\sigma},\bm{\tau}}]}|&=\frac{1}{|G|}\sum_{Q\in
G}|{\det[Q_{\bm{\sigma},\bm{\tau}}]}|^{2}\\\ &=\frac{1}{|G|}\sum_{Q\in
G}|\langle\phi_{2k}(Q)^{\mathsf{T}}e_{\bm{\sigma}},e_{\bm{\tau}}\rangle|_{2}^{2}\\\
&\geq\frac{\alpha}{|G|}\|e_{\bm{\tau}}\|_{2}^{2}>0\end{split}$ (218)

for all $\bm{\tau}\in\mathcal{C}_{2n,2k}$. Note that the transpose is
irrelevant since $G$ is a group and $\phi_{2k}$ an orthogonal representation.
While $\mathcal{M}_{G}$ may be made positive even without taking $G$ to be a
group, the group structure has desirable implications for the measurement
complexity (shadow norm) of our scheme. (It is also generally easier to sample
from a well-established group like the symmetric group, rather than some ad
hoc subset of $\mathrm{Sym}^{+}(2,2n)$.) To see this, we shall focus attention
exclusively to tight frames, which is motivated by the following observation.

######  Lemma 0.9.

Let $G\subseteq\mathrm{Sym}^{+}(2,2n)$ be a subset such that
$\phi_{2k}(G)e_{\bm{\sigma}}$ is a frame for all
$\bm{\sigma}\in\mathcal{D}_{2n,2k}$, $1\leq k\leq n$. Let
$\alpha_{k},\beta_{k}>0$ be the cumulative frame bounds over all $\bm{\sigma}$
for each $k$, i.e.,

$\alpha_{k}\leq\frac{1}{|G|}\sum_{\bm{\sigma}\in\mathcal{D}_{2n,2k}}\sum_{Q\in
G}|{\det[Q_{\bm{\sigma},\bm{\tau}}]}|^{2}\leq\beta_{k}$ (219)

for all $\bm{\tau}\in\mathcal{C}_{2n,2k}$. Then the shadow norm associated
with the ensemble $G$ satisfies

$\alpha_{k}/\beta_{k}^{2}\leq\|\Gamma_{\bm{\tau}}\|_{G}^{2}\leq\beta_{k}/\alpha_{k}^{2}.$
(220)

These bounds are saturated if and only if the cumulative frame is tight, i.e.,
$\|\Gamma_{\bm{\tau}}\|_{G}^{2}=\alpha_{k}^{-1}$.

######  Remark 0.10.

In Eq. 219, and hence in Eq. 215, the actual frame being used is
$\phi_{2k}(G)(e_{\bm{\sigma}}|G|^{-1/2})$, i.e., the “orbit” of the basis
vector scaled by a factor of $|G|^{-1/2}$. (We use scare quotes to indicate
that, since we do not require $G$ to be a group here, the corresponding set of
vectors may not be an orbit proper.)

###### Proof.

First, we recall the definition of the shadow norm:

$\|\Gamma_{\bm{\tau}}\|_{G}^{2}=\max_{\text{states
}\rho}\mathopen{}\left(\mathbb{E}_{Q\sim
G}\mathopen{}\left[\sum_{z\in\\{0,1\\}^{n}}\langle z|U(Q)\rho
U(Q)^{\dagger}|z\rangle\langle
z|U(Q)\mathcal{M}^{-1}(\Gamma_{\bm{\tau}})U(Q)^{\dagger}|z\rangle^{2}\right]\mathclose{}\right)\mathclose{}.$
(221)

Adapting Eq. 215 to the language of frames yields the inequalities

$\beta_{k}^{-1}\leq\|\mathcal{M}_{G}^{-1}(\Gamma_{\bm{\tau}})\|\leq\alpha_{k}^{-1},$
(222)

where we have used the fact that the spectral norm $\|\cdot\|$ of the Majorana
operators is $1$. Thus

$\begin{split}\|\Gamma_{\bm{\tau}}\|_{G}^{2}&\leq\alpha_{k}^{-2}\max_{\text{states
}\rho}\mathopen{}\left(\mathbb{E}_{Q\sim
G}\mathopen{}\left[\sum_{z\in\\{0,1\\}^{n}}\langle z|U(Q)\rho
U(Q)^{\dagger}|z\rangle\langle
z|U(Q)\Gamma_{\bm{\tau}}U(Q)^{\dagger}|z\rangle^{2}\right]\mathclose{}\right)\mathclose{}\\\
&=\alpha_{k}^{-2}\max_{\text{states }\rho}\mathopen{}\left(\mathbb{E}_{Q\sim
G}\mathopen{}\left[\sum_{z\in\\{0,1\\}^{n}}\langle z|U(Q)\rho
U(Q)^{\dagger}|z\rangle\mathopen{}\left(\sum_{\bm{\mu}\in\mathcal{C}_{2n,2k}}\det[(Q^{\mathsf{T}})_{\bm{\tau},\bm{\mu}}]\langle
z|\Gamma_{\bm{\mu}}|z\rangle\right)^{2}\mathclose{}\right]\mathclose{}\right)\mathclose{}.\end{split}$
(223)

Let us examine the innermost bracketed term. Since $\langle
z|\Gamma_{\bm{\mu}}|z\rangle=\pm 1$ when $\bm{\mu}\in\mathcal{D}_{2n,2k}$ and
$0$ otherwise, we can restrict the sum to run over $\mathcal{D}_{2n,2k}$.
Furthermore, by Lemma 0.6, we know that
$\det[(Q^{\mathsf{T}})_{\bm{\tau},\bm{\mu}}]$ is nonzero only for one
particular $\bm{\mu}$, which may or may not lie in $\mathcal{D}_{2n,2k}$. This
means that the sum contains at most one nonzero term, which takes the value
$\pm 1$, allowing us to “transform” the square into an absolute value:

$\mathopen{}\left(\sum_{\bm{\mu}\in\mathcal{C}_{2n,2k}}\det[(Q^{\mathsf{T}})_{\bm{\tau},\bm{\mu}}]\langle
z|\Gamma_{\bm{\mu}}|z\rangle\right)^{2}\mathclose{}=\sum_{\bm{\mu}\in\mathcal{D}_{2n,2k}}|{\det[Q_{\bm{\mu},\bm{\tau}}]}|.$
(224)

Importantly, this quantity does not depend on $z$, which makes the sum over
$z$ trivial (as well as the maximum over all states):

$\sum_{z\in\\{0,1\\}^{n}}\langle z|U(Q)\rho
U(Q)^{\dagger}|z\rangle=\operatorname{tr}\rho=1.$ (225)

We are then left with taking the average over $Q\sim G$ in Eq. 224, which is
simply the original quantity of interest obeying the frame bounds [Eq. 219].
Therefore

$\begin{split}\|\Gamma_{\bm{\tau}}\|_{G}^{2}&\leq\alpha_{k}^{-2}\operatorname*{\mathbb{E}}_{Q\sim
G}\mathopen{}\left[\sum_{\bm{\mu}\in\mathcal{D}_{2n,2k}}|{\det[Q_{\bm{\mu},\bm{\tau}}]}|\right]\mathclose{}\\\
&\leq\alpha_{k}^{-2}\beta_{k}.\end{split}$ (226)

To obtain the lower bound, simply reverse the roles of $\alpha_{k}$ and
$\beta_{k}$. Since the only inequalities invoked were those of the frame
bounds, it follows that equality holds in both directions if and only if the
frame is tight. ∎

Lemma 0.9 is useful in two ways: first, it gives an estimate on the shadow
norm for any valid subset of fermionic Gaussian Clifford unitaries, assuming
one has bounds on the eigenvalues of $\mathcal{M}_{G}$. Second, it tells us
that those subsets which give rise to _tight_ frames exhibit optimal sample
complexity in the sense of the shadow norm. Tight frames generated by the
action of finite groups have been completely characterized through
representation theory [vale2004tight, waldron2018introduction]. Below we
restate the primary results relevant to our context, which will motivate us to
restrict $G$ to a group.

######  Proposition 0.11 ([vale2004tight, Theorems 6.3 and 6.5]).

Let $H$ be a finite group, $V$ a Hilbert space, and $\varphi\colon
H\to\mathrm{U}(V)$ a unitary representation. Then:

  1. 1.

Every orbit $\varphi(H)v$, $v\in V\setminus\\{0\\}$, is a tight frame if and
only if every orbit spans $V$ [i.e., $\varphi$ is an irreducible
representation (irrep)].

  2. 2.

There exists $v\in V\setminus\\{0\\}$ for which $\varphi(H)v$ is a tight frame
if and only if there exists $w\in V\setminus\\{0\\}$ such that
$\operatorname{span}(\varphi(H)w)=V$.

In our context, this means that if $\phi_{2k}\big{|}_{G}$ is irreducible, then
every orbit is tight, hence saturating the bound of Lemma 0.9. Alternatively,
if $\phi_{2k}\big{|}_{G}$ is not an irrep, there may be only specific orbits
which form tight frames. Fortunately, this complication does not arise in our
setting due to the fact that we are dealing exclusively with (signed)
permutation matrices. We formalize this notion with the following lemma.
(Although we are only interested in the even-degree representations, we
formulate the statement to apply to all $1\leq k\leq 2n$ for completeness.)

######  Lemma 0.12.

Let $G\subseteq\mathrm{Sym}^{+}(2,2n)$ be a group such that