Sample-optimal classical shadows for pure states

Daniel Grier
Department of Mathematics and Department of Computer Science and Engineering, University of California, San Diego
Institute for Quantum Computing, University of Waterloo, Canada
Hakop Pashayan
Institute for Quantum Computing, University of Waterloo, Canada
Department of Combinatorics and Optimization, University of Waterloo, Canada
Perimeter Institute for Theoretical Physics, Waterloo, Canada
Dahlem Center for Complex Quantum Systems, Freie Universität Berlin, Germany
Luke Schaeffer
Institute for Quantum Computing, University of Waterloo, Canada
Department of Combinatorics and Optimization, University of Waterloo, Canada
Joint Center for Quantum Information and Computer Science, University of Maryland, College Park

We consider the classical shadows task for pure states in the setting of both joint and independent measurements. The task is to measure few copies of an unknown pure state $\rho$ in order to learn a classical description which suffices to later estimate expectation values of observables. Specifically, the goal is to approximate $\mathrm{Tr}(O \rho)$ for any Hermitian observable $O$ to within additive error $\epsilon$ provided $\mathrm{Tr}(O^2)\leq B$ and $\lVert O \rVert = 1$.
Our main result applies to the joint measurement setting, where we show $\tilde{\Theta}(\sqrt{B}\epsilon^{-1} + \epsilon^{-2})$ samples of $\rho$ are necessary and sufficient to succeed with high probability. The upper bound is a quadratic improvement on the previous best sample complexity known for this problem. For the lower bound, we see that the bottleneck is not how fast we can learn the state but rather how much any classical description of $\rho$ can be compressed for observable estimation.
In the independent measurement setting, we show that $\mathcal O(\sqrt{Bd} \epsilon^{-1} + \epsilon^{-2})$ samples suffice. Notably, this implies that the random Clifford measurements algorithm of Huang, Kueng, and Preskill, which is sample-optimal for mixed states, is not optimal for pure states. Interestingly, our result also uses the same random Clifford measurements but employs a different estimator.

§ INTRODUCTION

How many copies of an unknown state are required to construct a classical description of the state?  The answer to this question will depend on several details:  what constitutes an accurate description; what is already known about the state; and what restrictions are placed on the measurements of the state. Given the fundamental importance of this question, there has been significant prior work in bounding the number of samples of the states required to perform this learning task in a variety of contexts.

The most well-known setting is called quantum state tomography, where the goal is to learn enough about the state to be able to completely reconstruct it—precisely, estimate the unknown $d$-dimensional quantum state to accuracy $\epsilon$ in the Schatten $1$-norm. Tight upper and lower bounds for the number of copies required for this task are known: $\tilde{\Theta}(\epsilon^{-2} d^3)$ copies of the state are needed with independent measurements [11], and $\tilde{\Theta}(\epsilon^{-2} d^2)$ copies are needed when the unknown states can be simultaneously measured in a large joint measurement [18]. Independent measurements are easier to experimentally implement, while the joint measurements explore what is possible with respect to the fundamental limits of quantum mechanics. A key takeaway in the joint measurement setting is that the algorithm for the upper bound is achieving what would naïvely be the best possible result, given that a $d$-dimensional state has $\Theta(d^2)$-many independent parameters, and $\Theta(\epsilon^{-2})$ samples are necessary to estimate any one parameter.

In some sense, the requirements of the quantum tomography question are quite rigid. For many applications, only some properties of the unknown state are important. Can we get away with fewer samples if we relax our notion of approximation? In particular, what if we only wish to learn the expected values of certain Hermitian observables? Aaronson gave a somewhat surprising answer to this question in a joint measurement setting called shadow tomography [1]: given $M$ bounded observables
($\set{O_i}_{i=1}^M$, $\norm{O_i} \le 1$),
estimate $\tr{O_i \rho}$ to within $\epsilon$ additive error.[Aaronson actually stipulates that each observable is positive semi-definite matrix $E_i$ so that $\{E_i, I-E_i\}$ is a 2-outcome POVM.  We note that this is equivalent to the task of estimating expectation values of the (bounded) Hermitian observables $O_i$ via the mapping $E_i=(O_i+I)/2$.] In this setting, Aaronson showed that only $\tilde{\mathcal O}(\epsilon^{-4} \log^4 M \log d)$ samples of the state are needed. Subsequent work by Bădescu and O'Donnell [4] improved this to $\tilde{\mathcal O}(\epsilon^{-4} \log^2 M \log d)$, but there are still no matching lower bounds for this setting. That is, we do not know if we are extracting as much information about the unknown state as we can. In the independent measurement setting, $\tilde{\Theta}(\min\{M, d\}/\epsilon^2)$ samples are necessary and sufficient [5].

One subtlety concerning these observable estimation tasks is whether or not the measurements are allowed to depend on the specific observables $O_i$. In shadow tomography, the measurements can depend on the observables, but an increasingly popular setting (inspired by the work of Huang, Kueng, and Preskill [13]) is one in which the observables $O_i$ are unknown at the time of measurement. That is, the measurements must produce a classical description (called the classical shadow) from which the observable expected values can later be calculated. In their randomized Clifford measurement scheme, Huang, Kueng, and Preskill consider the independent measurement setting and show that $\Theta(B \epsilon^{-2} \log M)$ copies of the unknown state are both necessary and sufficient provided that $\tr{O_i^2}\leq B$ for all $i$ (note that $\tr{O_i^2}\leq d$).

Consider now how the classical shadows setting compares to the quantum state tomography setting with regard to the type of measurements allowed.  In the quantum state tomography setting, we know that joint measurements allow us to extract more information from the state, yielding estimates of the unknown state with provably fewer samples than those required with independent measurements.  In the classical shadows setting, however, it is not known how the type of measurement affects the number of samples required.  Concretely, is it possible to perform the classical shadows task with fewer samples if we switch to joint measurements?  We answer this question affirmatively in the setting of pure states.

Formally, we show the following: $\mathcal O((\sqrt{B}\epsilon^{-1}+\epsilon^{-2})\log M)$ samples of the unknown pure state are sufficient for performing the classical shadows task with constant probability of failure. Compared to [13], this achieves almost a square root reduction in sample complexity.

Remarkably, in analogy with the quantum state tomography setting, our joint measurement procedure is in some sense extracting the maximum amount of information possible.  To see this, consider a simple setting in which $B = d$, $\epsilon$ is constant, and we only wish to estimate a single observable. Our algorithm uses $\mathcal O(\sqrt d)$ samples. However, Gosset and Smolin [9] show that even if you are given the state as an explicit density matrix, you cannot compress your description of the state down to fewer than $\Omega(\sqrt{d})$-many bits of information in order to estimate arbitrary observable expectation values. Notice, however, that to successfully execute the classical shadows task, one would first need to learn such a compressed description through measurement of the unknown state. A priori, the number of measurements required to do this could be much higher than the size of this compressed description.
The fact that we find a matching upper bound implies that accessing the relevant information contained in the state is not the significant bottleneck.

We show that a similar phenomena exists for arbitrary parameters $B$ and $\epsilon$. Namely, we refine the Gosset-Smolin lower bound for compression to $\Omega(\sqrt{B}\epsilon^{-1})$-many bits, which ultimately allows us to show that $\tilde{\Omega}(\sqrt{B}\epsilon^{-1}+\epsilon^{-2})$ samples of the state are required for the classical shadows task. Therefore, our joint-measurement algorithm above is sample-optimal (at least for a single observable and up to log factors).

Finally, we address the classical shadows question with pure states and independent measurements. We show that $\mathcal O((\frac{\sqrt{Bd}}{\epsilon} + \frac{1}{\epsilon^2})\log M)$ copies of the state suffice. It's worth noticing that in certain parameter regimes, this upper bound is smaller than $\Theta(B \epsilon^{-2} \log M)$. In other words, our algorithm uses fewer samples than the classical shadows algorithm of Huang, Kueng, and Preskill which was designed for general mixed states. Indeed, their lower bound methods require the underlying state to be mixed.

§.§ The classical shadows task

We consider the classical shadows task introduced by Huang, Kueng, and Preskill [13]: given several copies of an unknown quantum state, produce a classical description of the state that is sufficiently representative to permit the reliable and accurate estimation of expectation values of some number of observables chosen from a broad class.

To formalise the task, let's begin with the class of observables we will use:

For any $B\in (0,d]$, let

$$\obs(B):=\set{O\in \mathbb C^{d \times d}~|~O=O^{\dagger},\inftynorm{O} = 1, \Tr(O^2)\leq B}.$$

In summary, these observables have been scaled/normalized so that $\inftynorm{O} = 1$ and have a bound of $B$ on their squared Frobenius norm $\Tr(O^2)$. The latter condition is due to the fact that $\Tr(O^2)$ is typically the dominant term in the sample complexity. We could also reasonably upper bound it by the rank of the observable since $\Tr(O^2) \leq \rank{O} \leq d$.

We remark that $\pnorm{O}{2} = \sqrt{\Tr(O^2)}$ and $\inftynorm{O}$ are examples of Schatten $p$-norms where $p = 2$ and $p = \infty$ respectively, but defined in general as $\pnorm{A}{p} := \Tr(|A|^{p})^{1/p}$ for $p \in [1, \infty)$. We will also use the Schatten $1$-norm. Going forward, we write $\pnorm{O}{1}$ for the $1$-norm, $\norm{O}$ for the infinity norm, and prefer $\Tr(O^2)$ over $\pnorm{O}{2}^2$.

The Classical Shadows Task consists of two separate phases—a measurement phase and an observable estimation phase—which are completed by two separate (randomized) algorithms, $\mathcal A_{\mathrm{meas}}$ and $\mathcal A_{\mathrm{est}}$, respectively. In addition to the inputs below, each algorithm also depends on the four parameters $s$, $B$, $\epsilon$, and $\delta$:

Measurement: $\mathcal A_{\mathrm{meas}} \colon \rho^{\otimes s} \to \{0,1\}^*$

Input:  $s$ copies of a state $\rho \in \mathbb C^{d \times d}$.

Output: A bit string called the classical shadow.

Estimation: $\mathcal A_{\mathrm{est}} \colon \obs(B) \times \{0,1\}^* \to \mathbb R$

Input: Observable $O\in \obs(B)$ and a classical shadow.

Output: Estimate $E \in \mathbb R$.

It's worth emphasizing that the input to the measurement algorithm is quantum (the state $\rho^{\otimes s}$) and the output is classical (the classical shadow). This output is computed from measuring the input state with some POVM (with arbitrary post-processing).

We say that $\mathcal A_{\mathrm{meas}}$ and $\mathcal A_{\mathrm{est}}$ constitute a valid protocol for the classical shadows task if their estimate for the expectation of the observable $E := \mathcal A_{\mathrm{est}}(O, \mathcal A_{\mathrm{meas}}(\rho^{\otimes s}))$ is such that

\begin{align}\label{eq:close_in_expectation}
\abs{\Tr(O \rho) - E} < \epsilon
\end{align}

with probability at least $1-\delta$ over the randomness of $\mathcal A_{\mathrm{meas}}$ and $\mathcal A_{\mathrm{est}}$.

Some may find it useful to think about the classical shadows task as a one-way communication protocol where one party (let's call her Melanie) is given copies of an unknown state and another party (say, Esteban) is given an observable. Melanie doesn't know Esteban's observable, and Esteban cannot send hints because we are assuming one way communication from Melanie to Esteban, so there is only one course of action: Melanie must measure her unknown state and send (over a classical channel) a description of the state from which Esteban can estimate the expected value of his given observable.

Throughout this paper, we will focus on the classical shadows task with unknown pure states. This motivates the following definitions:

Let $\task(B,\epsilon,\delta)$ to be the minimum number of samples $s$ required to successfully carry out the classical shadows task on pure states with the set of observables $\obs(B)$, to accuracy $\epsilon$, and failure probability at most $\delta$.

Sometimes we will omit $\delta$ and write $\task(B,\epsilon)$ to denote the minimum number of samples to achieve these tasks with some constant probability of failure, say, $0.001$.

Let $\itask(B,\epsilon,\delta)$ be the sample complexity for the classical shadows task with pure states when the measurement algorithm can only make independent measurements on the input state—that is, the measurement POVM is the tensor product of POVMs on single copies of the state. These POVMs do not have to be identical, but the entire state must be measured at the same time, or in other words, the output from a measurement on one copy of the state cannot influence the measurement on another.

We note that there are many possible variants for the sample complexity of the classical shadows task that we haven't given individual names. Most notably are the settings where the unknown states are mixed states (rather than pure) and/or the measurements are allowed to be adaptive (while still acting on single copies of the state).

§.§ Summary of results

Our main result is to prove matching upper and lower bounds on the sample complexity of performing the classical shadows task with respect to joint measurements and pure states.

$\task(B,\epsilon) = \bigthetalog{\frac{\sqrt{B}}{\epsilon} + \frac{1}{\epsilon^2}}$ provided $B \le \epsilon d$.

Notice that <Ref> consists of separate upper and lower bound results (for constant $\delta$). These match up to logarithmic factors in $B$ and $\epsilon^{-1}$, and the technical relationship between $B$, $\epsilon$, and $d$ is only required for the lower bound. In <Ref>, we will prove the upper bound, where we will also show that the dependence on the failure probability $\delta$ goes as $\log(1/\delta)$. We note that this dependence on $\delta$ implies that there are efficient protocols for the calculation of several observables simultaneously—that is, if the classical shadows task fails with probability at most $\delta$ on a single observable, then it fails with probability at most $M\delta$ on one or more out of $M$ observables by the union bound. In <Ref>, we will prove the lower bound where only the $\epsilon^{-2}$ term will scale with $\log(1/\delta)$.

We also prove an upper bound on the sample complexity of performing the classical shadows task with respect to independent measurements and pure states. Our upper bound can be compared to the matching upper and lower bound of Huang, Kueng, and Preskill [13] which applies to independent measurements and general states. In certain parameter regimes, our upper bound achieves a smaller sample complexity than the lower bound in [13] which implies that in the independent measurement setting, the classical shadows task has smaller sample complexity for pure states.

For all $\epsilon, \delta > 0$,

$$\itask(B, \epsilon, \delta) = \bigo{\min \left\{ \frac{B}{\epsilon^2}, \frac{\sqrt{Bd}}{\epsilon} + \frac{1}{\epsilon^2} \right \} \log(\delta^{-1})}.$$

We discuss and prove <Ref> in <Ref>.

Finally, we note that in all of our algorithms, the estimator $\hat{\rho}$ we use for the unknown state $\rho$ is not itself a proper state. In <Ref>, we show that this is a necessary price for the favorable sample complexity enjoyed by classical shadows schemes. Informally, we show that even for observables in $\obs(1)$, learning an estimate $\hat{\rho}$  that is a proper state to sufficient accuracy to solve the classical shadows task via the formula $\tr{O\hat{\rho}}$, requires a sample complexity that scales linearly in $d$, the dimension of the unknown state.

§ PRELIMINARIES

Here we cover key background material related to Haar random states, their moments, and the symmetric subspace. Throughout, we're working with qudits of dimension $d \ge 2$ unless otherwise specified. Since the unitary group $\mathrm{U}(d)$ acts on the Hilbert space of dimension $d$, it has a corresponding Haar measure which is invariant under the action of the group. Haar random states sampled proportional to this measure are ubiquitous in quantum information, and essential to define our measurement in Section <ref>.

To perform the necessary calculations on Haar random states, we need to discuss their moments, and some ancillary concepts.

For integer $k \geq 1$, $k$-th moment of an ensemble $\mathcal{E}$ of quantum states is

\E_{\ket{\psi} \sim \mathcal{E}}[ \ketbra{\psi}{\psi}^{\otimes k}].

An ensemble $\mathcal{E}$ is a (state) $t$-design if the moments $1 \leq k \leq t$ are identical to those of the Haar distribution (see Lemma <ref>).

Given a permutation $\pi \in \mathrm{S}_s$ (for $s \geq 1$), define a permutation operator $W_{\pi} \in \mathbb C^{d^{s} \times d^{s}}$ such that

W_{\pi} \ket{x_1} \cdots \ket{x_s} = \ket{x_{\pi^{-1}(1)}} \cdots \ket{x_{\pi^{-1}(s)}},

and extend by linearity.
That is, $W_{\pi}$ acts on $(\mathbb C^{d})^{\otimes s}$ by permuting the qudits, sending the qudit in position $i$ to position $\pi(i)$.

The symmetric subspace of an $s$-qudit system $(\mathbb C^{d})^{\otimes s}$ is the subspace invariant under $W_{\pi}$ for all $\pi \in \mathrm{S}_s$. We use $\kappa_s$ to denote its dimension and define $\sym^{(s)}$ to be the projector onto it (notationally omitting the dependence on $d$, the dimension of the qudit).

We have two characterizations of the symmetric subspace.

For all $s$,
$\sym^{(s)} = \frac{1}{s!} \sum_{\pi \in \mathrm{S}_s} W_{\pi}$, and
$\kappa_s = \binom{s+d-1}{d-1}$.

The integral of $\ketbra{\psi}{\psi}$ over the Haar measure is known from, e.g., [20].

\kappa_s \int_{\psi} \left( \ketbra{\psi}{\psi} \right)^{\otimes s} \mathrm{d}\psi = \sym^{(s)} = \frac{1}{s!} \sum_{\pi \in \mathrm{S}_{s}} W_{\pi}

where $\sym^{(s)}$ is the projector onto the symmetric subspace and $W_{\pi}$ is the operator that permutes $s$ qudits by an $s$-element permutation $\pi$.

We will often need to compute the (partial) trace of $(A_1 \otimes A_2 \otimes \cdots \otimes A_s) W_{\pi}$ for some linear operators $A_1, \dots, A_s \in \mathbb C^{d \times d}$. It turns out that there is an extremely useful tensor network based pictorial representation that simplifies these calculations. Let us give a brief introduction to those techniques, though readers may also find more thorough treatments useful [8, 19].

To start, we draw a single $d$ dimensional linear operator $A = \sum_{i,j \in [d]} a_{i,j} \ket{i}\bra{j}$ as a tensor block with a leg for the input and output indices for $A$:

[baseline=(current bounding box.center)]
(A) $A$;
(A.north) edge node[left,very near end] $i$ ++(0,1em);
(A.south) edge node[right,very near end] $j$ ++(0,-1em);

Suppose we have another tensor $B = \sum_{i,j \in [d]} b_{i,j} \ket{i}\bra{j}$. We express composition, tensor product, and trace as the following tensor networks:

[baseline=(current bounding box.center)]
(m) [matrix of math nodes, row sep=1em, column sep=1.5em, text height=0.8em,text depth=.2em]

(m-1-1) – (m-2-1);
(m-2-1.south) – ++(0,-1em);
(m-1-1.north) – ++(0,1em);

[baseline=(current bounding box.center)]
(m) [matrix of math nodes, row sep=, column sep=1em, text height=0.8em,text depth=.2em]

A     B

(m-1-1) – (m-2-1) – (m-3-1);
(m-1-2) – (m-2-2) – (m-3-2);

[baseline=(current bounding box.center)]
(m) [matrix of math nodes, row sep=, column sep=, text height=0.8em,text depth=.2em]

[rounded corners] (m-1-1.south) – ++(0,-1em) – ++(,0) |- ($(m-1-1.north) + (0,1em)$) – ++ (m-1-1);

Composition ($AB$)     Tensor Product ($A \otimes B$)     Trace ($\Tr(A)$)

The reason the tensor network picture is particularly nice for dealing with traces of $W_{\pi}$ terms is because each $W_{\pi}$ term is simply a permutation of wires in the tensor network picture. For example,
for a simple cyclic permutation, we have

\[
\begin{tikzpicture}[baseline=(current bounding box.center)]
\node[draw,rectangle,minimum width=4em] (W) {$W_{(1\,2\,3)}$};
\draw (W.north) ++ (1.5em,0) -- ++(0,1em);
\draw (W.north) ++ (0em,0) -- ++(0,1em);
\draw (W.north) ++ (-1.5em,0) -- ++(0,1em);
\draw (W.south) ++ (1.5em,0) -- ++(0,-1em);
\draw (W.south) ++ (0em,0) -- ++(0,-1em);
\draw (W.south) ++ (-1.5em,0) -- ++(0,-1em);
\end{tikzpicture}
\,\,\, =
\begin{tikzpicture}[baseline=(current bounding box.center)]
\matrix (m) [matrix of math nodes, row sep=\rowspacing, column sep=1em, text height=0.8em,text depth=.2em] {
\phantom{A} & \phantom{B} & \phantom{C} \\
\draw[rounded corners] ($(m-1-1.north)+(0,1em)$) -- (m-1-1.north) -- (m-1-3.south) -- ++(0,-1em);
\draw[rounded corners] ($(m-1-2.north)+(0,1em)$) -- (m-1-2.north) -- (m-1-1.south) -- ++(0,-1em);
\draw[rounded corners] ($(m-1-3.north)+(0,1em)$) -- (m-1-3.north) -- (m-1-2.south) -- ++(0,-1em);
\end{tikzpicture}.
\]

The key feature of tensor networks is that only the topology of the network matters, so we can simplify tensor networks just by moving the elements around. For example, consider a common partial trace  that will arise in this paper: $\Tr_1((A \otimes B) W_{(1 \, 2)})$. Drawing the tensor network, we get

\[
\Tr_1((A \otimes B) W_{(1\,2)}) \quad = \quad
\begin{tikzpicture}[baseline=(current bounding box.center)]
\matrix (m) [matrix of math nodes, row sep=\rowspacing, column sep=\colspacing, text height=0.8em,text depth=.2em] {
A & B \\
& \\

\draw (m-1-2.north) -- ++(0,1em);
\draw[rounded corners=0.5em] (m-1-1) -- (m-2-1.north) -- (m-2-2.south) -- ++(0,-1em);
\draw[rounded corners=0.5em] (m-1-2) -- (m-2-2.north) -- (m-2-1.south) -- ++(0,-1em) -- ++(-\loopwidth,0) |- ($(m-1-1.north) + (0,1em)$) -- ++ (m-1-1);
\end{tikzpicture}
\,=\,\,
\begin{tikzpicture}[baseline=(current bounding box.center)]
\matrix (m) [matrix of math nodes, row sep=\rowspacing, column sep=\colspacing, text height=0.8em,text depth=.2em] {
\phantom{B} & B \\
\phantom{A} & A \\

\draw (m-2-2.south) -- ++(0,-1em);
\draw (m-1-2.north) -- ++(0,1em);
\draw[rounded corners] (m-1-2.south) -- (m-2-1.north) -- (m-2-1.south) -- ++(0,-1em) -- ++ (-1em,0) |- ($(m-1-1.north)+(0,1em)$) -- (m-1-1.south) -- (m-2-2.north);
\end{tikzpicture}
\,=\,
\begin{tikzpicture}[baseline=(current bounding box.center)]
\matrix (m) [matrix of math nodes, row sep=\rowspacing, column sep=\colspacing, text height=0.8em,text depth=.2em] {
B \\
A \\

\draw (m-1-1) -- (m-2-1);
\draw (m-2-1.south) -- ++(0,-1em);
\draw (m-1-1.north) -- ++(0,1em);
\end{tikzpicture}
= \quad BA
\]

where we can push the $B$ tensor through the SWAP and around the trace loop to see that it is composed with $A$. In other words, we have just shown the identity $\Tr_1((A \otimes B) W_{(1 \, 2)}) = BA$.

As a generalization, we have the following useful fact:

Let $n \geq 1$ and $\pi = (1 \, 2 \, \cdots \, n) \in \symm_n$. For any $A_1, \ldots, A_n$, we have

\Tr_{-1}(W_{\pi}(A_1 \otimes A_2 \otimes \cdots \otimes A_n)) = A_n A_{n-1} \cdots A_1,

where $\Tr_{-1}$ indicates the partial trace of all but the first qudit. Thus, $\Tr(W_{\pi}(A_1 \otimes A_2 \otimes \cdots \otimes A_n)) = \Tr(A_n A_{n-1} \cdots A_1)$.

The fact is best seen with a small example. When $n = 5$, for instance, the tensor network diagram is

\[
\begin{tikzpicture}[baseline=(current bounding box.center)]
\matrix (m) [matrix of math nodes, row sep=\rowspacing, column sep=\colspacing, text height=0.8em,text depth=.2em] {
& & & & \\
A_1 & A_2 & A_3 & A_4 & A_5 \\

\draw[rounded corners=0.5em] ($(m-1-1.north) + (0,1em)$) -- (m-1-1.north) -- (m-1-5.south) -- (m-2-5.north);

\draw[rounded corners=0.5em] (m-2-2.south) -- ++(0,-1em) -- ++(\loopwidth,0) |- ($(m-1-2.north) + (0,1em)$) -- (m-1-2.north) -- (m-1-1.south) -- (m-2-1.north);

\draw[rounded corners=0.5em] (m-2-3.south) -- ++(0,-1em) -- ++(\loopwidth,0) |- ($(m-1-3.north) + (0,1em)$) -- (m-1-3.north) -- (m-1-2.south) -- (m-2-2.north);

\draw[rounded corners=0.5em] (m-2-4.south) -- ++(0,-1em) -- ++(\loopwidth,0) |- ($(m-1-4.north) + (0,1em)$) -- (m-1-4.north) -- (m-1-3.south) -- (m-2-3.north);

\draw[rounded corners=0.5em] (m-2-5.south) -- ++(0,-1em) -- ++(\loopwidth,0) |- ($(m-1-5.north) + (0,1em)$) -- (m-1-5.north) -- (m-1-4.south) -- (m-2-4.north);

\draw (m-2-1.south) -- ++ (0,-1em);
\end{tikzpicture}
\]

from which the identity follows.

§ JOINT MEASUREMENT UPPER BOUND

The goal of this section is to prove the following upper bound for the sample complexity of classical shadows for pure states and joint measurements.

$$\task(B,\epsilon,\delta) = \bigo{\left(\frac{\sqrt{B}}{\epsilon}+ \frac{1}{\epsilon^2}\right)\log \frac{1}{\delta}}.$$

The proof of <Ref> is constructive; given $B$, $\epsilon$, $\delta$ and $d$, we specify the number of samples and a pair of algorithms $\mathcal A_{\mathrm{meas}}$ and $\mathcal A_{\mathrm{est}}$ that solve the classical shadows task with that many samples.

In brief, the construction is as follows. We give a measurement $\mathcal{M}_{s}$ on $s$ copies of $\rho$, where the outcome of the measurement is a classical description of a pure state $\Psi$. We apply an affine transformation to the outcome $\Psi$ to produce an unbiased “shadow” estimator: a unital Hermitian matrix $\rhohat$ such that $\E[\rhohat] = \rho$.
Increasing the number of samples, $s$,  suppresses the additive error $\epsilon$ and the failure probability $\delta$ by a factor of $s^{-\bigo{1}}$. To improve this to an inverse exponential suppression in the failure probability,
we repeat the entire procedure $k = \mathcal O(\log(\delta^{-1}))$ times and take the median of the batch estimates akin to the median of means method [17, 16, 13]. A pseudocode description is given in Algorithms <ref> and <ref>.

Algorithm for $\mathcal{A}_{\mathrm{meas}}$ of  Theorem <ref>
[1]
Quantum state $\rho^{\otimes N}$, $B$, $\epsilon$, $\delta$, $d$.
Classical shadow $\set{\hat{\rho}^{(i)}}_{i \in [k]}$.
$s \gets \mathcal O(\sqrt{B} \epsilon^{-1} + \epsilon^{-2})$ Samples per batch
$k \gets \lfloor{N/s}\rfloor$ Number of batches
each batch $i = 1, \ldots, k$
$\psi_i \gets $ Measure new batch of $\rho^{\otimes s}$ with $\mathcal{M}_s$
$\rhohat^{(i)} \gets \frac{(d+s)\psi_i - I}{s}$
$\set{\hat{\rho}^{(i)}}_{i \in [k]}$

Algorithm for $\mathcal{A}_{\mathrm{est}}$ of Theorem <ref> and <Ref>
[1]
Classical shadow $\set{\hat{\rho}^{(i)}}_{i \in [k]}$ and observable $O$.
$E \gets \operatorname{median}(\Tr(O \rhohat^{(1)}), \ldots, \Tr(O \rhohat^{(k)}))$

We now define our measurement $\mathcal{M}_s$.

The standard symmetric joint measurement is a measurement on $s$ qudits. It is defined by the POVM $\mathcal{M}_{s} = \{ A_{\psi} \}_{\psi} \cup \{ I - \sym^{(s)} \}$ with elements

A_{\psi} := \kappa_s \ketbra{\psi}{\psi}^{\otimes s} \mathrm{d} \psi,

for all $d$-dimensional pure states, proportional to the Haar measure, plus a “fail” outcome $I - \sym^{(s)}$ for non-symmetric states.

We will be interested in the setting were $\rho$ is pure[This is the first time we use the purity of $\rho$ in our analysis, but certainly not the last.] and therefore $\rho^{\otimes s}$ is in the symmetric subspace, so we will never see the “fail” outcome—it exists solely to make the POVM sum/integrate to $I$.

The standard symmetric joint measurement is constructed from the Haar measure and can output a continuum of outcomes. This is technically inconsistent with <ref> where the measurement must output a finite length bit string. However, we note that there exist finitely supported ensembles of states known as $t$-designs that can be used to approximate the first $t$ moments of the Haar measure. As we will see, our upper bound uses up to the $s+2^{\mathrm{nd}}$ moment of the Haar measure to prove <ref>, where $s$ is the number of samples jointly measured. Thus any $s+2$ design with finitely many elements can be used in place of the Haar measure to construct a perfectly valid substitute, to the standard symmetric joint measurement, that also has finitely many outcomes.

§.§ Analysis

After defining the measurement, the estimator, and how many samples we need, the only remaining technical component is to bound the probability of failure. This ultimately comes down to Chebyshev's inequality:

$$\Pr\left[|\Tr(O \rhohat) - \E[\Tr(O \rhohat)]| \geq \epsilon \right] \leq \frac{\Var(\Tr(O \rhohat))}{\epsilon^2}.$$

Hence, we need to calculate the mean and variance of the random variable $\Tr(O \rhohat)$. To be precise, let $\rho$ be a pure state and suppose we measure $\rho^{\otimes s}$ with the standard symmetric joint measurement $\mathcal{M}_{s}$. Let $\Psi$ be the density matrix random variable for $\ketbra{\psi}{\psi}$, where $\psi$ is the outcome of the measurement. Let's start with the mean:

For measurement $\mathcal M_s$ on pure state $\rho^{\otimes s}$, we have

\begin{align*}
\E[\Psi] &= \frac{I + s\rho}{d+s}.
\end{align*}

To start, let's express the expectation as a Haar integral using the definition of $\mathcal M_s$:

\E[\Psi] = \int \psi \cdot \Pr[\Psi = \psi]
= \int \psi \cdot \Tr(A_\psi \rho^{\otimes s})
= \kappa_{s} \int \psi \cdot \Tr( \psi^{\otimes s} \rho^{\otimes s}) \, \dd\psi.

Using the identity $A \Tr(B) = \Tr_2 (A \otimes B)$ for all square matrices $A$ and $B$, we can apply <Ref> to compute the integral above:

\E[\Psi] = \kappa_{s} \int \Tr_{-1} \left( \psi^{\otimes s+1} \cdot (I \otimes \rho^{\otimes s}) \right) \, \dd \psi
= \frac{\kappa_{s}}{\kappa_{s+1}} \frac{1}{(s+1)!} \sum_{\pi \in \symm_{s+1}} \Tr_{-1} (W_{\pi} (I \otimes \rho^{\otimes s})).

We attack the right hand side by evaluating $\Tr_{-1}(W_{\pi}(I \otimes \rho^{\otimes s})$ for each $\pi$.
In particular, we will show that

\Tr_{-1}( W_{\pi} (I \otimes \rho^{\otimes s})) = \begin{cases}
I, & \text{if $\pi(1) = 1$,} \\
\rho, & \text{otherwise.}
\end{cases}

To do this, we take the cycle decomposition of $\pi$ and analyze each cycle separately. Notice that any cycle not involving position $1$ is completely traced out and the cycle operator acts on a tensor power of $\rho$ only, so Fact <ref> says the trace is $\Tr( \rho^{k} ) = \Tr( \rho ) = 1$ (since $\rho$ is pure). Thus, only the cycle through position $1$ matters. If $\pi(1) = 1$, then this cycle is trivial, and the result is $I$. Otherwise, the cycle visits $k \geq 1$ copies of $\rho$, leading to the product $\rho^k = \rho$.

There are $s!$ permutations which fix $1$ (i.e., $\pi(1) = 1$) and hence $s \cdot s!$ which do not, so we conclude that

\E[\Psi] = \frac{\kappa_{s}}{\kappa_{s+1}} \frac{1}{(s+1)!} \sum_{\pi \in \symm_{s+1}} \Tr_{-1} (W_{\pi} (I \otimes \rho^{\otimes s})) = \frac{s! \cdot I + s \cdot s! \cdot \rho}{(d+s) s!} = \frac{I + s \rho}{d+s}.

Notice that $\Tr( \E[\Psi] ) = \frac{\Tr(I) + s \Tr(\rho)}{d+s} = 1 = \E[\Tr(\Psi)]$, as a sanity check.

We now turn to the variance calculation, which depends on the second moment of the estimator:

For measurement $\mathcal M_s$ on pure state $\rho^{\otimes s}$, we have

\begin{align*}
\E[\Psi \otimes \Psi] &= \frac{2}{(d+s)(d+s+1)} \left( (I + s \rho)^{\otimes 2} - \frac{s(s+1)}{2} (\rho \otimes \rho) \right) \sym^{(2)}
\end{align*}

As in <Ref>, we evaluate $\E[\Psi \otimes \Psi]$ as

\begin{align*}
\E[\Psi \otimes \Psi]
&= \int (\psi \otimes \psi) \cdot \Pr[\Psi = \psi] \\
&= \int (\psi \otimes \psi) \cdot \kappa_{s} \Tr( \psi^{\otimes s} \rho^{\otimes s}) \, \dd \psi \\
&= \kappa_{s} \int \Tr_{-1,2} \left( \psi^{\otimes s+2} \cdot (I^{\otimes 2} \otimes \rho^{\otimes s}) \right) \, \dd \psi \\
&= \frac{\kappa_{s}}{\kappa_{s+2}} \frac{1}{(s+2)!} \sum_{\pi \in \symm_{s+2}} \Tr_{-1,2} (W_{\pi} (I^{\otimes 2} \otimes \rho^{\otimes s})),
\end{align*}

where the partial trace $\Tr_{-1,2}$ now preserves the first two qudits. In <Ref> and for the special case of $s=1$, we show a complete derivation of how this trace simplifies using the tensor network notation, which may be useful to some readers before proceeding to the more general proof.

Let us evaluate the sum term-by-term. We divide the permutations into two types: those where $1$ and $2$ appear in separate cycles (type A) and those where $1$ and $2$ appear in the same cycle (type B). Consider the type A permutations first:

\Tr_{-1,2}(W_{\pi} (I^{\otimes 2} \otimes \rho^{\otimes s})) = \begin{cases}
I \otimes I & \text{if $\pi(1) = 1$ and $\pi(2) = 2$} \\
I \otimes \rho & \text{if $\pi(1) = 1$ and $\pi(2) \neq 2$} \\
\rho \otimes I & \text{if $\pi(1) \neq 1$ and $\pi(2) = 2$} \\
\rho \otimes \rho & \text{if $\pi(1) \neq 1$ and $\pi(2) \neq 2$}
\end{cases}

Second moment calculation in the special case $s = 1$.

As before, we end up getting $I$ or $\rho$ for each position, depending on whether $1$ and $2$ were fixed by the permutation. The combinatorics is similar but not identical: there are $(s+1)!$ permutations which fix $1$, and of those, $s!$ fix $2$ and $s \cdot s!$ do not. Likewise, $s \cdot s!$ fix $2$ but not $1$. The remainder of the type A permutations fix neither $1$ nor $2$, and to count these we need a fact.

Let $\pi \in \symm_n$ be a permutation, and consider $\pi' := (1 \, 2) \pi$. Then

\text{$1$ and $2$ are in distinct cycles of $\pi$} \iff \text{$1$ and $2$ are in the same cycle of $\pi'$}.

In other words, there is a bijection between A and B permutations, so exactly half of all permutations ($\tfrac{1}{2} (s+2)!$) are type A, and half are type B. It follows that there are $\frac{s(s-1)}{2} \cdot s!$ type A permutations such that $\pi(1) \neq 1$ and $\pi(2) \neq 2$. Therefore, the overall contribution of the type A permutations is equal to $s!$ times

(I \otimes I) + s(I \otimes \rho) + s(\rho \otimes I) + \frac{s(s-1)}{2} (\rho \otimes \rho) = (I + s \rho)^{\otimes 2} - \frac{s(s+1)}{2} \rho^{\otimes 2}.

Fortunately, we do not have to repeat this counting argument for the type B permutations. The bijection from <Ref> decomposes each type B permutation $\pi$ as $(1\,2)\pi'$ where $\pi'$ is type A, and so

\Tr_{-1,2}(W_{\pi}(I \otimes \rho^{\otimes s})) = \Tr_{-1,2}(W_{(1\,2)}W_{\pi'}(I \otimes \rho^{\otimes s})) = W_{(1\,2)} \Tr_{-1,2}(W_{\pi'}(I \otimes \rho^{\otimes s})).

Therefore, we can multiply our result for type A permutations by $W_{(1)(2)} + W_{(1\,2)} = 2 \sym^{(2)}$ to get the total. The result follows from some careful accounting of the scalar factors.

Let $\rhohat = \frac{(d+s) \Psi - I}{s}$. For any observable $O \in \obs(B)$ and $\epsilon > 0$, we bound the probability of failure as

\Pr[|\Tr(O \rhohat) - \Tr(O \rho)| \geq \epsilon] \leq \frac{1}{\epsilon^2 s^2} \left[ \Tr( O^2 ) + 8s \Tr(O^2 \rho) \right].

Our goal will be to compute the mean and variance of the estimate $\Tr(O \rhohat)$ in order to apply Chebyshev's inequality. By <Ref>, we have

\E[\rhohat] = \E \left[\frac{(d+s) \Psi - I}{s} \right] = \rho,

and so the mean of the estimate is correct: $\E[\Tr(O \rhohat)] = \Tr(O \E[\rhohat]) = \Tr(O \rho)$.

To analyze the variance, let us first consider traceless observables $O$, where $\Tr(O) = 0$. As usual, it will be useful to break the variance into a first moment term $\E[\Tr(O \Psi)]$ and a second moment term $\E[\Tr(O \Psi)^2]$:

\begin{align*}
\Var(\Tr(O \rhohat)) = \Var \left( \frac{(d+s)\Tr(O \Psi) - \Tr(O)}{s}\right)
= \left( \frac{d+s}{s} \right)^2 (\E[\Tr(O \Psi)^2] - \E[\Tr(O \Psi)]^2)
\end{align*}

Putting aside the $s^{-2}$ factor for now, the (squared) first moment term is

(d+s)^2 \E[\Tr(O\Psi)]^2 = (d+s)^2 \left(\frac{\Tr(O) + s\Tr(O \rho)}{d+s}\right)^2 = s^2 \Tr(O \rho).

For the second moment term, we use $\E[\Tr(O \Psi)^2] = \E[\Tr((O\otimes O) (\Psi \otimes \Psi))]$ and <Ref> to write

\begin{align*}
(d+s)^2 \E[\Tr(O \Psi)^2]
&= \frac{2(d+s)}{d+s+1} \Tr \left[ O^{\otimes 2}  \left( (I + s \rho)^{\otimes 2} - \frac{s(s+1)}{2} \rho^{\otimes 2} \right) \sym^{(2)} \right] \\
&\le \Tr \left[ O^{\otimes 2}  \left(I^{\otimes 2} + s (I \otimes \rho + \rho \otimes I) - \frac{s(s-1)}{2} \rho^{\otimes 2} \right) (2\sym^{(2)}) \right].
\end{align*}

Recall that $(2\sym^{(2)}) = W_{(1)(2)} + W_{(1\,2)}$, so to simplify, consider the contribution of those two terms:

\begin{align*}
W_{(1)(2)}:&\;\;\Tr(O)^2 + 2s \Tr(O) \Tr(O \rho) + \frac{s(s-1)}{2} \Tr(O \rho)^2 \\
W_{(1\,2)}:&\;\;\Tr(O^2) + 2s \Tr(O^2 \rho) + \frac{s(s-1)}{2} \Tr((O \rho)^2)
\end{align*}

In fact, because $\rho$ is pure, we have[If $\rho := \proj{\psi}$, we get
$\Tr((O\rho)^2) = \Tr(O\proj{\psi}O\proj{\psi}) = \Tr(\bra{\psi}O\ket{\psi}\bra{\psi}O\ket{\psi}) = \bra{\psi}O\ket{\psi}^2 = \Tr((O\rho)^2).$] that $\Tr((O\rho)^2) = \Tr(O \rho)^2$. Combining all of the above, we arrive at a bound for the (scaled) variance of $\Tr(O \Psi)$:

(d+s)^2 \Var(\Tr(O \Psi))
\le \Tr(O^2) + 2s \Tr(O^2 \rho) - s \Tr(O \rho)^2 \\
\le \Tr(O^2) + 2s \norm{O^2}

where the last inequality uses Hölder's inequality.[
Hölder's inequality: For density matrix $\rho$ and observable $O$,
$\Tr(O \rho) \le \sonenorm{O\rho} \leq (\inftynorm{O}\sonenorm{\rho}) \leq \inftynorm{O}.$

It follows that $\Var(\Tr(O \rhohat))$ is

\Var(\Tr(O \rhohat)) \leq \frac{\Tr(O^2) + 2s \| O^2 \|}{s^2},

for traceless observables. Now suppose $O$ has nonzero trace, and let $O_0 := O - \Tr(O) I/d$ be its traceless part. Naturally, we have

\begin{align*}
\Var(\Tr(O \rhohat))
&= \Var(\Tr(O_0 \rhohat)) \leq \frac{\Tr(O_0^2) + 2s \| O_0^2 \|}{s^2}.
\end{align*}

Let's now show that we can bound each of those terms ($\Tr(O_0^2)$ and $\| O_0^2 \|$) using functions of the original observable $O$. First, for the $\Tr(O_0^2)$ term, we have that

\Tr(O_0^2) = \Tr(O^2) - \Tr(O)^2/d \leq \Tr(O^2).

Next consider the $\| O_0^2 \|$ term. We have that $\Tr(O) \le \|O\| d$, and so the largest eigenvalue (in absolute value) of $O - \Tr(O)I/d$ is at most $2\|O\|$. We get

\| O_0^2 \| = \| O_0\|^2 \le (2\| O \|)^2 = 4 \| O^2 \|.

Hence $\Var(\Tr(O \rhohat)) \leq \tfrac{1}{s^2} \left( \Tr(O)^2 + 8s \| O^2 \| \right)$, and the result follows by Chebyshev's theorem.

At last, we can prove the main theorem for this section.

Consider an arbitrary observable $O$, and use Corollary <ref> to bound the probability a single batch estimate $\rhohat^{(i)}$ is wrong by

\Pr[|\Tr(O \rhohat^{(i)}) - \Tr(O \rho)| \geq \epsilon] \leq \frac{1}{\epsilon^2 s^2} \left[ \Tr( O^2 ) + 8s \| O^2 \| \right] \leq \frac{B+8s}{\epsilon^2 s^2}.

Suppose we want this probability to be less than some constant $p < 1/2$; we leave it to the reader to check that at most
$\mathcal O(1/(\epsilon^2 p) + \sqrt{B/(\epsilon^2 p)})$
samples suffice, and note that $s$ is chosen accordingly in Algorithm <ref>.

Recall that our final estimate $E$ is

E := \operatorname{median}(\Tr(O \rhohat^{(1)}), \ldots, \Tr(O \rhohat^{(k)})).

Assume there are an odd number of batches, so the median is actually some $\Tr(O \rhohat^{(i)})$.
If $E$ is a bad estimate, i.e., $|E - \Tr(O \rho)| > \epsilon$ then at least $k/2$ of the batch estimates are wrong: either $E$ and the estimates higher than it, or $E$ and the estimates lower than it.

The batches are independent so Chernoff bounds the chance of seeing $\geq k/2$ failures.

\Pr[|E - \Tr(O \rho)| \geq \epsilon] \leq \Pr \left[\# \{ i : |\Tr(O \rhohat^{(i)}) - \Tr(O \rho)| \geq \epsilon \} \geq \frac{k}{2} \right] \leq \sqrt{4p(1-p)}^{k}

Setting this less than the failure probability $\delta$, we have

k \geq \frac{\log\delta^{-1}}{\log \, (4p(1-p))^{-1/2}}.

Again, we note that $k$ is set accordingly in Algorithm <ref>.

§.§ Discussion

Let us compare this result with the original classical shadows protocol of Huang, Kueng, and Preskill [13]. Their algorithm measures each copy of $\rho$ with $\mathcal{M}_1$[Technically, they use a Clifford unitary instead of a Haar random unitary, presumably for the sake of efficient implementation. However, $\mathcal{M}_1$ will work in place of their measurement, and the analysis is identical since it uses up to third moments of the ensemble of unitaries, which are the same for Clifford vs. Haar random unitaries, i.e., the Cliffords are a $3$-design [21, 23, 15, 24].], producing unbiased single-copy estimates $\rhohat_1, \ldots, \rhohat_s$ for $\rho$, which are then averaged into a batch estimate $\rhohat = \frac{1}{s} \sum_{i=1}^{s} \rhohat_i$. Given the observable $O$, the estimate is then $\Tr(O \rhohat)$, or the median of several batches, if necessary to reduce the probability of failure.

We have just seen that the variance of a single-copy estimate is $\Var(\Tr(O \rhohat_i)) \leq \Tr(O^2) + \Tr(O^2 \rho)$, and averaging $s$ estimates together reduces the variance by a factor of $\frac{1}{s}$. On the other hand, our measurement with $\mathcal{M}_s$ provides an unbiased estimate with variance

\[
\Var(\Tr(O \rhohat)) \leq \frac{\Tr(O^2)}{s^2} + \frac{\Tr(O^2 \rho)}{s}.
\]

Since $\Tr(O^2 \rho)\leq 1$, we see that the quadratic denominator of $\Tr(O^2)$ (which is the dominant term) is making all the difference.

§ JOINT MEASUREMENT LOWER BOUND

$\task(B, \epsilon, \delta) = \Omega \left( \frac{\sqrt{B}}{\epsilon \log(B+1)} + \frac{\log \delta^{-1}}{\epsilon^2} \right)$ provided $B \le \epsilon d$.

Notice that this bound matches the $\mathcal O((\sqrt{B}\epsilon^{-1} + \epsilon^{-2})\log(\delta^{-1}))$ upper bound up to a $\log(B)$ and a $\log(1/\delta)$ factor. We prove this as two separate lower bounds: $\Omega(\frac{\log \delta^{-1}}{\epsilon^2})$ and $\Omega( \frac{\sqrt{B}}{\epsilon \log(B+1)})$.

The first lower bound ($\Omega(\epsilon^{-2}\log(\delta^{-1})$) is derived via a reduction from the problem of distinguishing two pure states, $\rho_0$ and $\rho_1$, at trace distance $2 \epsilon$ from each other. We then use the known performance of the optimal measurement (Helstrom  measurement).

The second lower bound is shown via a reduction from a problem in communication complexity known as Boolean Hidden Matching [7]. We will show that any protocol for the classical shadows task implies a protocol for the Boolean Hidden Matching problem, which has known communication complexity lower bounds. These communication lower bounds will imply that the classical shadow must contain a significant amount of information. However, Holevo's theorem gives an upper bound on the amount of information gained through measurement. Therefore, in order to successfully complete the classical shadows task, many copies of the unknown state are required.

§.§ Omega(epsilon-2) lower bound

The proof of the lower bound uses known results relating the trace distance between two states with our ability to distinguish the states by observables or binary measurements. In particular, the maximum gap for the expectation of a positive semi-definite observable is equal to the trace distance between the states:

For arbitrary states $\rho$ and $\sigma$,

\max_{0 \preceq O \preceq I} |\Tr( O \rho ) - \Tr( O \sigma )| = \tfrac{1}{2} \| \rho - \sigma \|_1.

Furthermore, there is an optimal $O$ satisfying $\Tr(O^2) \leq \tfrac{1}{2} \rank (\rho - \sigma) \leq \tfrac{1}{2}(\rank \rho + \rank \sigma)$.

Diagonalize $\rho - \sigma$ as $\sum_{i} \lambda_i \ketbra{\phi_i}{\phi_i}$. Define two positive semi-definite observables, $O_{+}$ and $O_{-}$.

\begin{align*}
O_{+} &:= \sum_{i : \lambda_i > 0} \ketbra{\phi_i}{\phi_i} &
O_{-} &:= -\sum_{i : \lambda_i < 0} \ketbra{\phi_i}{\phi_i}
\end{align*}

Clearly $O_{+} - O_{-}$ is a projector onto the eigenvectors of $\rho - \sigma$, so $\Tr((O_{+} - O_{-})(\rho - \sigma)) = \Tr(\rho - \sigma) = 0$, and $\rank (O_{+} - O_{-}) = \rank (\rho - \sigma)$. On the other hand,

\Tr((O_{+} + O_{-})(\rho - \sigma)) = \Tr \left( \sum_{i : \lambda_i \neq 0} (\sgn \lambda_i) \lambda_i \ketbra{\phi_i}{\phi_i} \right) = \sum_{i : \lambda_i \neq 0} |\lambda_i| = \| \rho - \sigma \|_{1}.

It follows that $\Tr(O_{+} (\rho - \sigma)) = \Tr(O_{-}(\rho - \sigma)) = \tfrac{1}{2} \| \rho - \sigma \|_1$. Since $O_{+}$ and $O_{-}$ are orthogonal, the rank of their sum, $\rank (\rho - \sigma)$, is the sum of their ranks. Hence, we can take whichever of $O_{+}$ and $O_{-}$ has rank at most $\tfrac{1}{2} \rank (\rho - \sigma)$.

Separately, we know the optimal measurement for distinguishing a uniformly randomly chosen $\rho$ or $\sigma$ is given by:

The optimal measurement for distinguishing states $\rho$ and $\sigma$ succeeds with probability
$\frac{1}{2} + \frac{1}{4} \| \rho - \sigma \|_1$.

$\task(1, \epsilon, \delta) = \Omega(\epsilon^{-2} \log(1/\delta))$.

Let $\rho_0$ and $\rho_1$ be pure states with trace distance $\tfrac{1}{2} \| \rho_0 - \rho_1 \|_1 = 2 \epsilon$.
We claim that a protocol for the classical shadows task to $\epsilon$-approximate the expected values of observables of rank 1 with probability of failure at most $\delta$ can be used to distinguish states $\rho_0, \rho_1$ with probability of failure at most $\delta$. First, apply the measurement subroutine to unknown state $\rho_b^{\otimes s}$ to produce the classical shadow. Then, use the observable $O$ from Lemma <ref> to estimate $\Tr(O \rho_b)$. If the estimate is closer to $\Tr(O \rho_0)$ then guess $b=0$, otherwise guess $b=1$. Since the gap $|\Tr(O \rho_0) - \Tr(O \rho_1)| = \tfrac{1}{2} \| \rho_0 - \rho_1 \| = 2\epsilon$, we succeed whenever the gap between the estimate and $\Tr(O \rho_b)$ is less than $\epsilon$.

We can see this classical shadows protocol as a binary measurement distinguishing $\rho_0$ and $\rho_1$. Since it succeeds with probability $1 - \delta$, the optimal distinguishing measurement from Lemma <ref> must do better, so

\[
\left(1 - 2\delta \right)^2 \leq \frac{1}{4} \| \rho_0^{\otimes s} - \rho_1^{\otimes s} \|_1^2
= 1 - \Tr( \rho_0^{\otimes s} \rho_1^{\otimes s} )
= 1 - \Tr( \rho_0 \rho_1 )^{s}
= 1 - \left(1 - \frac{1}{4} \| \rho_0 - \rho_1 \|_1^{2} \right)^{s}
\!= 1 - (1 - \epsilon^2)^{s}
\]

where we have used the equation $\frac{1}{2} \| \rho_0 - \rho_1 \| = \sqrt{1 - \Tr( \rho_0 \rho_1 )}$ relating trace distance and fidelity for pure states [6]. Rearranging, we have $(1 - \epsilon^2)^{s} \leq 1 - (1 - 2 \delta)^2 = 4 \delta - 4 \delta^2 \leq 4 \delta$. Taking logs and using $1 - \frac{1}{x} \leq \ln x$ we have

-\frac{s\epsilon^2}{1 - \epsilon^2} \leq s \log(1 - \epsilon^2) \leq \log (4\delta),

or equivalently,

s \geq \frac{1-\epsilon^2}{\epsilon^2} \log \left(\frac{1}{4 \delta}\right) = \Omega(\epsilon^{-2} \log(1/\delta)).

§.§ Omega(epsilon-1 B/log(B+1)) lower bound

The lower bound is proved via a reduction from a problem in communication complexity to the classical shadows task. That is, a protocol for the classical shadows task that has low sample complexity will translate to a communication protocol for another task (to be specified) that has low communication complexity. We will then use existing lower bounds on the communication complexity of that problem to prove lower bounds on the classical shadows task. This lower bound technique draws considerable inspiration from that in the work of Gosset and Smolin for compressing classical descriptions of quantum states [9]. Our proof can be seen as generalization of their techniques and also a translation to the setting of sample complexity rather than bit complexity.

§.§.§ One-way communication complexity

Because our proof is based on principles from communication complexity, let's briefly introduce that topic. We are interested in one-way communication protocols where two parties—Alice and Bob—are trying to jointly compute some function $f \colon \mathcal X \times \mathcal Y \to \{0,1\}$. Alice is given some input $x \in \mathcal X$ and Bob is given some input $y \in \mathcal Y$. Alice's goal is to send a single message $m \in \{0,1\}^*$ to Bob, so that he can compute $f(x,y)$. Of course, she could choose to send her entire input $x$, but in many cases it may be possible to communicate fewer bits and still be successful.

To be precise about the size of the message Alice must send, let $X$ and $Y$ be the random variables (possibly correlated) for the inputs of Alice and Bob, respectively, and let $M$ be the random variable for Alice's message to Bob. Notice that implicit in $M$ is Alice's communication strategy, which may be an arbitrary (randomized) function of her input. Let's start with the easiest setting, where Alice and Bob run deterministic algorithms.

$D^{(X,Y)}_\delta(f)$, the bounded-error deterministic one-way communication complexity of $f$, is the minimum number of bits that Alice must send to Bob to compute $f$ with at most $\delta$ probability of error whenever their inputs are chosen according to the distribution $(X,Y)$.

A natural variant of classical one-way protocols is when Alice and Bob are allowed to run randomized algorithms. There are two settings: private-coin protocols, where Alice and Bob each have access to private random strings; and public-coins protocols, where Alice and Bob also have access to a shared random string along with their private strings. It will not be critical to completely understand the nuances of the various types of protocols for our proof, but we define them in order to precisely state the theorems on which the lower bound rests.

$R_\delta(f)$, the bounded-error randomized one-way communication complexity of $f$, is the minimum number of bits that Alice must send to Bob with a public-coin protocol to compute $f$ over all possible inputs with failure probability at most $\delta$.

While randomized strategies may seem more powerful than deterministic strategies, Yao's minimax principle shows that there is always some input distribution for which the randomized and deterministic complexities coincide:

$\max_{(X,Y)} D^{(X,Y)}_\delta (f) = R_\delta(f)$.

It turns out that we will eventually be interested in the amount of information contained in Alice's message $M$, not just the length, which is what is measured by the communication complexity. That said, these two quantities are intuitively related—if the information $I(M : X)$ that Alice's message $M$ reveals about her input $X$ is much lower than the number of bits she is communicating, she should be able to send a smaller message and still be successful. The following theorem formalizes this message compression idea:

$D^{(X,Y)}_{\Delta+\delta}(f) = 2 \Delta^{-1} \left[ \min I(M : X) + O(1) \right]$ where the minimization is over all one-way private-coin protocols for $f$ with input distribution $(X,Y)$ and probability of error at most $\delta$.

In the next section, we will show that classical shadows must also contain a lot of information, which will be the basis of our lower bound.

§.§.§ Classical shadows for Boolean Hidden Matching

Our starting point is a lower bound for the one-way communication complexity for the Boolean Hidden Matching function $\BHM_\alpha \colon \mathcal X \times \mathcal Y \to \{0,1\}$ which was introduced by Gavinsky, Kempe, Kerenidis, Raz, and de Wolf.

$R_\delta(\BHM_\alpha) = \Theta(\sqrt{n/\alpha})$ for any constant $\delta < 1/2$.

Recall that our lower bound technique is to show that a classical shadows strategy with few samples implies a communication protocol for the Boolean Hidden Matching function with low complexity. To do this, let's look more closely at the $\BHM_\alpha$ function.

It will be useful describe the function directly as a communication problem with the inputs $x \in \mathcal X$ for Alice and $y \in \mathcal Y$ for Bob. Alice is given $(0,1)$-assignments for the $n$ vertices of a graph, and Bob is given $(0,1)$-assignments to $\alpha n$ vertex-disjoint edges in the graph. A set of vertex-disjoint edges from a graph is called a matching, hence the name of the function. Importantly, the function is only defined on inputs for Alice and Bob that satisfy the following promise: for each edge in the matching, the parity of the connected vertices (from Alice's input) plus Bob's edge bit assignment is some constant $b\in \set{0,1}$. The output of the function is then defined as this bit $b$.

Formally, for every $n \ge 2$, the Boolean Hidden Matching function $\BHM_\alpha \colon \mathcal X \times \mathcal Y \to \{0,1\}$ with parameter $\alpha \in (0,1/4]$ is defined on inputs $\mathcal X = \{ 0, 1 \}^{n}$ and $\mathcal Y = (\mathbb N, \mathbb N)^{\alpha n} \times \{ 0, 1 \}^{\alpha n}$ as follows:

Alice:     $x \in \mathcal X$.

Bob:     $y = (\mathcal M, w) \in \mathcal Y$, where $\mathcal{M} = \{ (i_1, j_1), \ldots, (i_{\alpha n}, j_{\alpha n})\}$ is a matching and $w \in \{ 0, 1 \}^{\alpha n}$ is a bit assignment.

Promise:     There exists $b \in \{ 0, 1 \}$ such that $b = x_{i_k} \oplus x_{j_k} \oplus w_{k}$ for all $k$.

Output:     $b \in \{0,1\}$

In the setting where $\alpha$ is constant, [7] show that the quantum communication complexity of the Boolean Hidden Matching problem is low—Alice only needs to send a $(\log n)$-qubit state. Gosset and Smolin [9] notice that this implies the existence of a set of states and observables whose expectation values give solutions to the Boolean Hidden Matching function. We generalize this observation to the non-constant $\alpha$ setting below:

Protocol for Boolean Hidden Matching using classical shadows subroutines (shown in blue). From her input $x$, Alice prepares $\rho_x^{\otimes s}$, measures, and sends the classical shadow to Bob. From his input $y$, Bob computes $O_y$ and then estimates $\Tr(O_y \rho_x)$ to accuracy $\alpha$ from the classical shadow that Alice sent to him. He then uses that estimate to answer the Boolean Hidden Matching problem. Correctness follows from the fact that $\Tr(O_y \rho_x) = 2\alpha \cdot \BHM_\alpha(x,y)$. For the details of $\rho_x$ and $O_y$ see <Ref>.

There is an ensemble of states $\{\rho_x \in \mathbb C^n\}_{x \in \mathcal X}$ and observables $\{O_y \in \mathrm{Obs}(\alpha n)\}_{y \in \mathcal Y}$ such that $\Tr(O_y \rho_x) = 2 \alpha \cdot \BHM_\alpha(x,y)$. Furthermore, a protocol for the classical shadows task for observables of squared Frobenius norm $B:= \alpha n$, estimation accuracy $\epsilon := \alpha$, and failure probability $\delta$ implies a one-way private-coin protocol for Boolean Hidden Matching with failure probability $\delta$.

Given valid inputs $x$ and $y = (\mathcal M, w)$ to the $\BHM_\alpha$ function, define the pure state

\ket{\psi_x} := \frac{1}{\sqrt{n}} \sum_{i=1}^n (-1)^{x_i} \ket{i}

and the observable

O_y := \sum_{k=1}^{\alpha n} \frac{1}{2}(\ket{i_k} - (-1)^{w_k} \ket{j_k})(\bra{i_k} - (-1)^{w_k} \bra{j_k}).

Notice that $O_y \in \obs(\alpha n)$ since $O_y$ is a $\alpha n$-rank projector. Letting $\rho_x := \ketbra{\psi_x}{\psi_x}$, we get

\Tr( O_y \rho_x ) = \bra{\psi_x} O_y \ket{\psi_x} = \frac{1}{n} \sum_{k=1}^{\alpha n}(1 - (-1)^{x_{i_k} \oplus x_{j_k} \oplus w_k })= 2 \alpha b = 2\alpha \BHM_\alpha(x,y).

In particular, this implies that if $E$ is an $\alpha$-approximation to $\Tr( O_y \rho_x )$, then

|E - \Tr( O_y \rho_x )| < \alpha \implies \left|\frac{E}{2\alpha} - \BHM_\alpha(x,y)\right| < 1/2,

or in other words, rounding $E/(2\alpha)$ is equal to $\BHM_\alpha(x,y)$.

We now claim that the existence of these states and observables implies a private-coin one-way protocol for the Boolean Hidden Matching problem (see <Ref>): Suppose we want a protocol for $\BHM_\alpha$ with probability of failure at most $\delta$. Let $s = \task(\alpha n, \alpha, \delta)$. On input $x$, Alice prepares the state $\rho_x^{\otimes s}$, measures it with a valid classical shadows strategy, and sends the resulting classical shadow to Bob. On input $y$, Bob computes the observable $O_y$, and then computes an estimate $E$ for $\Tr(O_y \rho_x)$ using the classical shadow sent by Alice. The correctness of the classical shadows strategy implies that $E$ is an $\alpha$-approximation to $\Tr(O_y \rho_x)$ with probability of failure at most $\delta$. As shown above, Bob can then compute $\BHM_\alpha(x,y)$ with failure probability at most $\delta$ by appropriately rounding the estimate.

Let us now note a key property of the one-way protocol in <Ref> for the Boolean Hidden Matching problem. Namely, once Alice prepares $\rho_x^{\otimes s}$, she no longer uses her original input $x$. Her message (the classical shadow) only depends on her measurement of the state $\rho_x^{\otimes s}$. In particular, if her message is to contain a lot of information about her input $x$, then Holevo's theorem stipulates that she must be measuring a state of high dimension, or, in other words, $s$ must be large:

Let $Z$ be the classical outcome of measuring a $d$-dimensional state drawn from an ensemble $\{ \rho_x \}_{x \in \mathcal X}$ according to $x \sim X$. Then, $I(X : Z) \le \log d$.

Naïvely, the states $\rho_{x}^{\otimes s}$ in <Ref> consist of $s$ qudits of dimension $n$, i.e., it is a space of dimension $n^s$. However, since each $\rho_x^{\otimes s}$ is invariant under permutation, it belongs to the symmetric subspace, which has dimension nearly a factor of $s!$ smaller (see <Ref>).

We are now ready to put all of the pieces of the lower bound together:

$\task(B, \epsilon) = \Omega\left(\frac{\sqrt{B}}{\epsilon \log(B+1)}\right)$ provided $B \le \epsilon d$.

Using the communication complexity of $\BHM_\alpha$ as our starting point, we first show that Alice's message to Bob in every successful protocol for the Boolean Hidden Matching problem must contain a significant amount of information. To show this, note that by Yao's minimax principle (Theorem <ref>), there exists a distribution $(X,Y)$ such that $R_\delta(\BHM_{\alpha}) = D^{(X,Y)}_\delta(\BHM_{\alpha})$, for any given $\delta$ is the probability of error. <Ref> lets us upper bound the (deterministic) complexity with mutual information, and <Ref> proves a lower bound. Thus,

\begin{align*}
D^{(X,Y)}_{\delta}(\BHM_\alpha) &= O(\min I(M : X) + 1), \text{ and} \\
D^{(X,Y)}_{\delta}(\BHM_\alpha) &= \Omega(\sqrt{n/\alpha}),
\end{align*}

for any constant $\delta$. It follows that $I(M : X) = \Omega(\sqrt{n/\alpha})$ for any one-way private-coin protocol for $\BHM_\alpha$.

Now consider the classical shadows strategy for solving $\BHM_\alpha$ as described by <Ref>, and suppose that Alice measures $s = \task(\alpha n, \alpha)$ copies of $\rho_x$.[It is possible that Alice could use fewer than $\task(\alpha n, \alpha)$ samples for this specific application to Boolean Hidden Matching, but it will be useful later that $s$ is big enough to estimate a broader class of observables (specifically those used in <Ref>). ]  Recall that Alice's message $M$ depends only on her measurement of $\rho_X^{\otimes s}$ which has classical outcome $Z$. We get

I(X : M) \leq I(X : Z) \leq \log(\dim \sym^{(s)}) \le O(s \log(n/s + 1)).

where we have used (in order) the Data Processing Inequality, Holevo's theorem (<Ref>), the dimension of the symmetric subspace (<Ref>), and the following inequality:

\dim \sym^{(s)} = \binom{n+s-1}{n-1} \leq \binom{n+s}{n} = \binom{n+s}{s} \leq \left( \frac{e(n+s)}{s} \right)^{s}.

Notice that we now have both an upper bound and a lower bound for the mutual information between Alice's input and her message for a one-way protocol for $\BHM_\alpha$ with constant error probability. Setting $\epsilon := \alpha$ and $B := \epsilon n$, we have

\begin{align*}
I(X : M) &= \Omega(\sqrt{B}/\epsilon), \text{ and} \\
I(X : M) &= O(s \log(B/(\epsilon s) + 1)).
\end{align*}

It follows that

\[
s = \Omega\left(\frac{\sqrt{B}}{\epsilon \log(B/(\epsilon s) + 1)}\right).
\]

Notice that if we substitute any lower bound for $s$ in the RHS of the equation above, then we get a new lower bound for $s$ on the LHS. Unfortunately, plugging in the trivial lower bound ($s \ge 1$) is not very tight. Instead, we will use the $s = \Omega(1/\epsilon^2)$ lower bound from the previous section. To justify this, notice that

\[
s = \task(B,\epsilon) \ge \task(1,\epsilon) = \Omega(1/\epsilon^2),
\]

where we have used <Ref> and the fact that $B \ge 1$ since $\| O \| = 1$. Therefore, we can plug $s = \Omega(1/\epsilon^2)$ into the RHS above to arrive at the following:

s = \Omega\left(\frac{\sqrt{B}}{\epsilon \log(B\epsilon + 1)}\right).

Assuming $\epsilon \leq 1$, we simplify this to $s = \Omega( \frac{\sqrt{B}}{\epsilon \log(B+1)})$.[While this simplification apparently makes lower bound weaker for no reason, it doesn't actually effect the overall lower bound. In regimes where $\log(B\epsilon + 1)$ is significantly smaller than $\log(B + 1)$, the additive $1/\epsilon^2$ term in the lower bound becomes dominant.]

Finally, we point out that the construction of <Ref> operates in the regime where the states have dimension $d := n$ and the observables are of rank $B = \epsilon d$. One can extend the lower bound to apply to all observables of rank $B \le \epsilon d$ by embedding the states used in <Ref> into a subspace of dimension $n \leq d$ and keeping the observables the same.

§ INDEPENDENT MEASUREMENT UPPER BOUND

Since the global Clifford group acting on qubits is a $3$-design, the randomized Clifford measurement classical shadows algorithm of Huang, Kueng, and Preskill [13] can be viewed as simulating independent $\mathcal{M}_1$ measurements on all copies of $\rho$ then, constructing an unbiased estimator from the measurement outcome on each copy. Their result is for independent measurements and general mixed states, but it upper bounds pure states as a special case.

For all $\epsilon, \delta > 0$,

$$\itask(B, \epsilon, \delta) = \bigo{ \frac{B\log(\delta^{-1})}{\epsilon^2}}.$$

Huang, Kueng, and Preskill also show a matching lower bound, but the hard instances they construct are with states of full rank. We give an independent measurement classical shadows algorithm for pure states which is better in certain parameter regimes (and is no worse).

For all $\epsilon, \delta > 0$,

$$\itask(B, \epsilon, \delta) = \bigo{\min \left\{ \frac{B}{\epsilon^2}, \frac{\sqrt{Bd}}{\epsilon} + \frac{1}{\epsilon^2} \right \} \log(\delta^{-1})}.$$

Suppose $\delta$ is constant, and consider the case where $B = d^{1/3}$ and $\epsilon = d^{-2/3}$. The Huang, Kueng, Preskill algorithm achieves a sample complexity of $s = \bigo{B/\epsilon^2} = \bigo{d^{5/3}}$, whereas this algorithm uses only $\bigo{d^{4/3}}$ samples. While this regime is not especially natural, it nevertheless demonstrates that the Huang-Kueng-Preskill analysis is not optimal for pure states.

As it turns out, however, our measurement algorithm is the same—on each copy of $\rho$, we make an independent measurement with the POVM $\mathcal{M}_1$, which can be performed with a random Clifford measurement since we only use third moments of the Haar measure. The difference is in how we construct our estimator for the unknown state. To see this, first let $\Psi_1, \ldots, \Psi_s$ be the Hermitian random variables for the measurement outcomes. Using Lemma <ref>, notice that $\rhohat_{i} := (d+1)\Psi_i - I$ is an unbiased estimator for the unknown state, i.e., $\E[\rhohat_i] = \rho.$ The average of the $\rhohat_i$'s, i.e.,

\hat{X} := \frac{1}{s} \sum_{i=1}^{s} \rhohat_i

is effectively the Huang-Kueng-Preskill estimator. Our key observation is that when $\rho$ is pure, $\rhohat_i \rhohat_j$ is also an unbiased estimator of $\rho$:

\E[\rhohat_i \rhohat_j] = \E[\rhohat_i] \E[\rhohat_j] = \rho^2 = \rho.

where we have used the independence of the measurements for the first equality and the purity of $\rho$ for the last.
In light of this, we consider an estimator $\hat{Y}$ defined to be the average of the $s(s-1)$ quadratic terms where $i \neq j$.

\hat{Y} := \frac{1}{s(s-1)} \sum_{i \neq j} \rhohat_i \rhohat_j.

To analyze the accuracy of the estimator $\hat{Y}$, we will once again turn to Chebyshev's inequality:

\begin{align*}
\Pr[|\Tr(O \hat{Y}) - \Tr(O \rho)]| \geq \epsilon] \leq \frac{\Var(\Tr(O \hat{Y}))}{\epsilon^2}.
\end{align*}

Expanding out the variance term using the definition of $\hat Y$, we get

\[
\Var( \Tr(O \hat{Y}) ) = \frac{1}{s^2 (s-1)^2} \sum_{i \neq j} \sum_{k \neq \ell} \Cov( \Tr(O \rhohat_i \rhohat_j), \Tr(O \rhohat_k \rhohat_{\ell})).
\]

We need to bound all of these covariance terms to bound the variance. When all indices $i,j,k,\ell$ are distinct, then the covariance is 0 (by independence). For the other four cases, we rely on corollaries <ref>, <ref>, <ref>, and <ref>, which we summarize in the following lemma (proof in <Ref>):

For each combination of $i,j,k,\ell$, $\Cov(\Tr(O \rhohat_i \rhohat_j), \Tr(O \rhohat_k \rhohat_{\ell}))$ is

* One index matches $(|\{ i, j \} \cap \{ k, \ell \}| = 1)$

* Match in different positions $(i = \ell$ or $j = k)$: $\mathcal O(\|O\|^2)$

* Match in same position $(i = k$ or $j = \ell)$: $\mathcal O(\|O\|^2)$

* Both indices match $(|\{ i, j \} \cap \{ k, \ell \}| = 2)$

* Order swapped $(i = \ell$ and $j = k)$: $\mathcal O(Bd)$

* Same order $(i = j$ and $k = \ell)$: $\mathcal O(Bd)$

Since $\|O\|^2 \le 1$, the contribution from the first two terms is extremely small compared to the last term. This gives us the following bound on the variance:

$\Var(\Tr(O \hat{Y})) = \mathcal O(\frac{Bd}{s^2} + \frac{1}{s})$.

Expand the variance as

\[
\Var( \Tr(O \hat{Y}) ) = \frac{1}{s^2 (s-1)^2} \sum_{i \neq j} \sum_{k \neq \ell} \Cov( \Tr(O \rhohat_i \rhohat_j), \Tr(O \rhohat_k \rhohat_{\ell})).
\]

Using <Ref>, we account for the contribution of each type of covariance term to get

\Var(\Tr(O \hat{Y})) = \bigo{\frac{0\cdot s^4 + \|O\|^2 \cdot s^3 + \|O\|^2 \cdot s^3 + Bd \cdot s^2 + Bd \cdot s^2}{s^4}} = \bigo{Bd/s^2 + 1/s}

where we have used that $\|O\|^2 \le 1$ and that there are $\mathcal O(s^4)$ terms where $i,j,k,\ell$ are distinct; $\mathcal O(s^3)$ terms where exactly one index matches; and $\mathcal O(s^2)$ terms where both indices match.

Putting everything together, we can now prove the claimed sample complexity in Theorem <ref>.

We first point out that the sample complexity of our new estimator is only better (or at least no worse) when $\epsilon \leq \sqrt{B/d}$, so when that does not hold we simply use the original $\hat{X}$ estimator of Huang, Kueng, and Preskill.[Actually, it is better to smoothly transition between the estimators $\hat{X}$ and $\hat{Y}$ using a convex combination rather than use a sharp threshold. However, the improvement is only a constant factor and would require computing covariance of linear vs. quadratic terms (e.g., $\Cov(\Tr(O \rhohat_i \rhohat_j), \Tr(O \rhohat_k))$) to justify rigorously. Hence, we simply use one or the other.]

Otherwise, we use Algorithm <ref> to measure the state and construct several $\hat{Y}$ estimators (line <ref>, constituting the classical shadow. This shadow is then used in Algorithm <ref> for the observable estimation step, which once again uses the median-of-means method where the analysis will be identical to that in the proof of <Ref>.

It suffices to analyze the variance of the estimator constructed within each batch, which is $\bigo{Bd/s^2 + 1/s}$ by <Ref>. To apply Chebyshev's inequality, we need the variance to be at most $\epsilon^2$, which occurs when we have at least $s = \mathcal O(\sqrt{Bd}/\epsilon + 1/\epsilon^2)$ samples.

§ OPEN PROBLEMS

Our upper and lower bounds almost completely settle the question of sample complexity for the classical shadows task with arbitrary measurements and arbitrary observables, but clearly many questions remain. First, can the remaining discrepancies between our upper and lower bounds be removed? That is, can we get the lower bounds to have the correct dependence on $\delta$ (which we conjecture to be $\log(\delta^{-1})$) and remove a $\log B$ factor?

Second, the sample complexity of learning states in the context of tomography is known for all combinations of independent vs. joint measurements, and pure vs. mixed states. Can we characterize the sample complexity of classical shadows as thoroughly? Ref. [13] gives matching bounds for independent measurements and mixed states, and our result gives matching bounds for joint measurements on pure states, but the independent/pure and joint/mixed cases are open. At the very least, we know that in the independent measurement setting that the pure state case is different than the mixed state setting since our upper bound in <Ref> is smaller than the lower bound in Ref. [13] for some regimes.

As a follow up, tomography bounds are sometimes stated as a function of $r$, the rank of the unknown state. For example, the sample complexity is $\tilde{\Theta}(rd \epsilon^{-1})$ for joint measurements, capturing the pure state $(r=1)$ and worst-case mixed state ($r = d$) behaviour simultaneously. We believe the sample complexity of the classical shadows task depends smoothly on $r$, but do not yet have a conjecture.

Third, we do not describe how to concretely implement our large joint measurements. We may replace the continuum of Haar random pure states and elements $A_{\psi}$ in the POVM $\mathcal M_s$ with a concrete $(s+2)$-design, but even then it is not clear how to practically implement the measurement. Alternatively, can we make a similar POVM with a simpler ensemble of states, e.g., $t$-designs for some $t< s+2$, and update the analysis to achieve an equivalent end result? We note that for our exact measurement, a lower bound on $t$ can be computed by using upper bounds on the number of bits required to describe a state $t$-design. To see this, notice that the outcome of the POVM $\mathcal M_s$ suffices as a compressed description of the state for the purpose of observable estimation. Since we proved a state compression lower bound of $\Omega(\sqrt{B}\epsilon^{-1})$-many bits, the number of bits required to specify the state must be at least this large. In particular, this implies that you could not implement our measurement on $n$ qubits with a $3$-design since Clifford states can be specified using $\mathcal O(n^2)$-many bits.

Fourth, there is the question of robustness to error or noise. For any classical shadows protocol, we can ask how it behaves when the samples are not exactly of the form $\rho^{\otimes s}$, due to variation in samples. For our pure state protocols, we are also interested in how fast our algorithms degrade when given mixed states that are close to pure.

Fifth, Ref. [13] introduced two classical shadows protocols known as the random Pauli measurements and the random Clifford measurements schemes. The former schemes targets local observable. The latter scheme, like our protocols, targets observables with low Frobenius norm. These target classes of observables are mutually exclusive and each scheme achieves lower sample complexity with respect to its target class observable. Recent work [3, 2] has also focused on the development of an intermediate scheme that achieves favourable sample complexity scaling in both target classes of observables. All these works focus on general states and independent measurements. Our work identifies an optimal protocol for the low Frobenius norm class of observable in the setting of pure states and joint measurements. So it is natural to consider the pure states and/or joint measurements setting in the context of local observables or the combined class.

Finally, one may consider cubic or higher order generalizations of the quadratic estimator used in the proof of <ref>. We leave the analysis of such estimators to future work.

§ ACKNOWLEDGEMENTS

We thank Anurag Anshu, Stephen Bartlett, Daniel Burgarth, David Gosset, David Gross and Richard Kueng for useful discussions.
Research at Perimeter Institute is supported in part by the Government of Canada through the Department of Innovation, Science and Economic Development Canada and by the Province of Ontario through the Ministry of Colleges and Universities. HP also acknowledges the support of the Natural Sciences and Engineering Research Council of Canada (NSERC) discovery grants [RGPIN-2019-04198] and [RGPIN-2018-05188].

[1]
Scott Aaronson.
Shadow tomography of quantum states.
In Proceedings of the 50th Annual ACM SIGACT Symposium on Theory
of Computing, STOC 2018, page 325–338, New York, NY, USA, 2018.
Association for Computing Machinery.

[2]
A. A Akhtar, H.-Y. Hu, and Y.-Z. You.
Scalable and flexible classical shadow tomography with tensor
arXiv preprint arXiv:2209.02093, 2022.

[3]
Christian Bertoni, Jonas Haferkamp, Marcel Hinsche, Marios Ioannou, Jens
Eisert, and Hakop Pashayan.
Shallow shadows: Expectation estimation using low-depth random
clifford circuits.
arXiv preprint arXiv:2209.12924, 2022.

[4]
Costin Bădescu and Ryan O'Donnell.
Improved quantum data analysis.
In Proceedings of the 53rd Annual ACM SIGACT Symposium on Theory
of Computing, pages 1398–1411, 2021.

[5]
Sitan Chen, Jordan Cotler, Hsin-Yuan Huang, and Jerry Li.
Exponential separations between learning with and without quantum
In 2021 IEEE 62nd Annual Symposium on Foundations of Computer
Science (FOCS), pages 574–585. IEEE, 2022.

[6]
C.A. Fuchs and J. van de Graaf.
Cryptographic distinguishability measures for quantum-mechanical
IEEE Transactions on Information Theory, 45(4):1216–1227,

[7]
Dmitry Gavinsky, Julia Kempe, Iordanis Kerenidis, Ran Raz, and Ronald De
Exponential separations for one-way quantum communication complexity,
with applications to cryptography.
In STOC'07, Proceedings of the Annual ACM Symposium on Theory
of Computing, pages 516–525, 2007.

[8]
David Gross, Felix Krahmer, and Richard Kueng.
A partial derandomization of phaselift using spherical designs.
Journal of Fourier Analysis and Applications, 21(2):229–266,

[9]
David Gosset and John Smolin.
A compressed classical description of quantum states.
In Wim van Dam and Laura Mancinska, editors, 14th Conference on
the Theory of Quantum Computation, Communication and Cryptography (TQC
2019), volume 135 of Leibniz International Proceedings in Informatics
(LIPIcs), pages 8:1–8:9, Dagstuhl, Germany, 2019. Schloss
Dagstuhl–Leibniz-Zentrum fuer Informatik.

[10]
Carl W. Helstrom.
Quantum Detection and Estimation Theory.
Academic Press Inc., New York, NY, USA, 1976.

[11]
Jeongwan Haah, Aram W. Harrow, Zhengfeng Ji, Xiaodi Wu, and Nengkun Yu.
Sample-optimal tomography of quantum states.
In Proceedings of the Forty-Eighth Annual ACM Symposium on
Theory of Computing, STOC '16, page 913–925, New York, NY, USA, 2016.
Association for Computing Machinery.

[12]
Prahladh Harsha, Rahul Jain, David McAllester, and Jaikumar Radhakrishnan.
The communication complexity of correlation.
In Twenty-Second Annual IEEE Conference on Computational
Complexity (CCC'07), pages 10–23. IEEE, 2007.

[13]
Hsin-Yuan Huang, Richard Kueng, and John Preskill.
Predicting many properties of a quantum system from very few
Nature Physics, 16(10):1050–1057, Oct 2020.

[14]
Alexander Semenovich Holevo.
Bounds for the quantity of information transmitted by a quantum
communication channel.
Problemy Peredachi Informatsii, 9(3):3–11, 1973.

[15]
Richard Kueng and David Gross.
Qubit stabilizer states are complex projective 3-designs, 2015.

[16]
M. Lerasle.
Lecture notes: Selected topics on robust statistical learning theory,

[17]
G. Lugosi and S. Mendelson.
Mean estimation and regression under heavy-tailed distributions: A
Found. Comp. Math., 19:1145–1190, 2019.

[18]
Ryan O'Donnell and John Wright.
Efficient quantum tomography.
In Proceedings of the Forty-Eighth Annual ACM Symposium on
Theory of Computing, STOC '16, page 899–912, New York, NY, USA, 2016.
Association for Computing Machinery.

[19]
Daniel A Roberts and Beni Yoshida.
Chaos and complexity by design.
Journal of High Energy Physics, 2017(4):1–64, 2017.

[20]
Andrew J Scott.
Tight informationally complete quantum measurements.
Journal of Physics A: Mathematical and General, 39(43):13507,

[21]
Zak Webb.
The Clifford group forms a unitary 3-design.
Quantum Information and Computation, 16(15-16):1379–1400, Nov

[22]
Andrew Chi-Chin Yao.
Probabilistic computations: Toward a unified measure of complexity.
In 18th Annual Symposium on Foundations of Computer Science
(SFCS 1977), pages 222–227. IEEE Computer Society, 1977.

[23]
Huangjun Zhu.
Multiqubit Clifford groups are unitary 3-designs.
Physical Review A, 96(6), Dec 2017.

[24]
Huangjun Zhu, Richard Kueng, Markus Grassl, and David Gross.
The Clifford group fails gracefully to be a unitary 4-design, 2016.

§ PROPER LEARNING DISCUSSION

The classical shadows task is a learning problem: given many samples of a quantum state, we make measurements with the goal of learning the state well enough to approximate arbitrary observables. A learning problem is said to be proper if it requires the learned representation to be from the same class—in our case, the class of pure states—as the original object. Our classical shadows algorithm fails to be proper on several counts:

* The output of the problem is a real number, not the classical description of a quantum state. By definition, the classical shadows task is not a proper learning task.

* Internally, our algorithm does produce Hermitian matrices $\rhohat^{(i)}$ with trace $1$ (the shadows), which have the potential to represent a quantum state. However, each such estimate is

* high rank, so it cannot represent a pure state, and

* not positive semi-definite, so it cannot represent a mixed state.

* The algorithm uses multiple $\rhohat^{(i)}$ shadows, and takes a median of their estimates on each observable. Even if each $\rhohat^{(i)}$ were a quantum state, there may not exist a state exactly consistent with all our medians on a set of observables.

On the other hand, in the limit where the failure probability is extremely small, we can afford to estimate the expectation of the state on an $\epsilon$-cover of the observables. From these values, we can approximate the original state to accuracy $\epsilon$ in trace distance. In this regime, the problem is equivalent to tomography, which is a proper learning problem.

In this section, we show that any proper learning algorithm for classical shadows would require significantly more samples. Our starting point is a known lower bound for the quantum state tomography question for pure states:

Any quantum algorithm that takes copies of an unknown pure state $\rho$ and outputs a classical estimate $\rhohat$ such that $\pnorm{\rho - \rhohat}{1} \le \epsilon$ with constant failure probability $\delta < 1$ requires $\Omega(d \epsilon^{-2} / \log(d/\epsilon))$ samples.

In particular, we show that a proper classical shadows algorithm implies a state tomography algorithm.

Suppose there exists a quantum learning algorithm that, given $s$ copies of an unknown $d$-dimensional state $\rho$, outputs a classical description of a trace $1$, Hermitian PSD matrix $\hat{\rho}$ such that, for all $O\in \obs(1)$ with failure probability $\delta < 1$:

\begin{align}
\abs{\tr{O\rho}-\tr{O\hat{\rho}}}\leq \epsilon.
\end{align}

Then, $s = \tilde\Omega(d/\epsilon)$.

Run the algorithm and feed it $O = \rho$. We have that $\pnorm{\rho-\hat{\rho}}{1}\leq \sqrt{8\epsilon}$ since

\begin{align*}
\epsilon &\geq |  \Tr(O \rho) - \Tr(O \rhohat) | \\
&= |\Tr( \rho^2 )- \Tr( \rho \rhohat )| \\
&= |1 - \Tr( \rho \rhohat )| & \text{$\rho = \rho^2$ since $\rho$ is pure} \\
&= |1 - F( \rho,  \rhohat )^2| & \text{$\Tr(\rho \sigma) = F(\rho, \sigma)^2$ if either is pure} \\
&\geq 1 - F(\rho, \rhohat)^2 \\
&\geq \tfrac{1}{8} \| \rho - \rhohat \|_1^{2} & \text{Fuchs-van de Graaf inequality}
\end{align*}

where $F(\rho, \sigma) := \Tr(\sqrt{\sqrt{\rho} \sigma \sqrt{\rho}})$ is the fidelity of $\rho$ and $\sigma$.

Hence, if we solve the classical shadows task to error $\epsilon$ on this observable, then we have estimated $\rho$ to within Schatten $1$-norm distance $\mathcal{O}(\sqrt{\epsilon})$. That is, classical shadows learner is also a quantum tomography algorithm, so we can use the known lower bound from <Ref>. It follows that a proper learning algorithm for the classical shadows task requires $\tilde\Omega(d/\epsilon)$ samples.

In the commonly considered regime where $d\gg\epsilon^{-1}$, the above lower bound is significantly more than both our algorithm (from Section <ref>) and the original classical shadows algorithm, which use only  $\bigo{\tfrac{1}{\epsilon^2}}$ samples.

§ COVARIANCE BOUNDS

The goal of this subsection is to prove <Ref>, which gives each of the covariance terms $\Cov(\Tr(O \rhohat_i \rhohat_j), \Tr(O \rhohat_k \rhohat_{\ell}))$ that appears in the expansion of the estimator $\hat Y$. Because $\rhohat_j$ appears twice in all of the covariance terms that are non-zero, it will be convenient to explicitly calculate the second moment of $\rhohat_j$.

For all $j$, the second moment of $\rhohat_j$ is

\E[\rhohat_j^{\otimes 2}] = \left( I \otimes I + I \otimes \rho + \rho \otimes I \right) \left( W_{(1\,2)} - \frac{2}{d+2}\sym^{(2)} \right).

Recall that $\rhohat_j$ is obtained through an independent and identical measurement process, so it suffices to analyze a specific $\rhohat := (d+1)\Psi - I$ term. By <Ref> and <Ref>, we compute the first and second moment of $\Psi$ for the special case $s = 1$ as

\E[\Psi] = \frac{I + \rho}{d+1} \; \text{ and } \;
\E[\Psi \otimes \Psi] = \frac{I \otimes I + I \otimes \rho + \rho \otimes I}{(d+1)(d+2)} \left( W_{(1)(2)} + W_{(1\,2)} \right)

Now, expanding out the second moment, we get

\begin{align*}
\E[\rhohat^{\otimes 2}] &= \E[ ((d+1)\Psi - I)^{\otimes 2}] \\
&= (d+1)^2 \E[ \Psi \otimes \Psi] - (d+1)(\E[\Psi] \otimes I + I \otimes \E[\Psi]) + I \otimes I \\
&= \frac{d+1}{d+2} \left(I \otimes I + I \otimes \rho + \rho \otimes I \right)(W_{(1)(2)} + W_{(1\,2)}) - \left(I \otimes I + I \otimes \rho + \rho \otimes I \right) \\
&= \left(I \otimes I + I \otimes \rho + \rho \otimes I \right) \left( \frac{(d+1) W_{(1\,2)} - W_{(1)(2)}}{d+2} \right),
\end{align*}

where we recall that $W_{(1)(2)} = I \otimes I$ to obtain the last equality. We arrive at the lemma by writing the final line in terms of the symmetric subspace $\sym^{(2)} = (W_{(1\,2)} + W_{(1)(2)})/2$.

Our goal will now be to express all covariance terms in a manner such that we can apply <Ref>. Since these equations can become quite cumbersome to write out fully, we will often drop the “$\otimes$” symbol in expressions with $I$ and $\rho$. For example,

I \otimes I + \rho \otimes I + I \otimes \rho \to (II + I \rho + \rho I)

will be a common abbreviation. We enclose these abbreviations in parentheses when they are multiplied with other terms.

Let's first tackle the covariance terms $\Cov(\Tr(O \rhohat_i \rhohat_j), \Tr(O \rhohat_k \rhohat_{\ell}))$ where there is only 1 index shared, i.e., $(|\{ i, j \} \cap \{ k, \ell \}| = 1)$. There are two subcases: a match in different positions ($i = \ell$ or $j = k$); or a match in same position ($i = k$ or $j = \ell$). While the proofs are quite similar, we break them into two separate corollaries.

For all distinct $i, j, k$,

$$\Cov( \Tr( O\rhohat_i \rhohat_j ), \Tr( O\rhohat_j \rhohat_k ) ) \leq 2\Tr(O \rho)^2 \leq 2\| O \|^2 .$$

First, translate covariance to a second moment calculation:

\begin{align*}
\Cov( \Tr( O\rhohat_i \rhohat_j ), \Tr( O\rhohat_j \rhohat_k ) )&= \E[\Tr( O\rhohat_i \rhohat_j ) \Tr( O\rhohat_j \rhohat_k )^{*}] - \Tr(O \rho) \Tr(O \rho)^{*} \\
&= \E[\Tr( O\rhohat_i \rhohat_j ) \Tr( O\rhohat_k \rhohat_j )] - \Tr(O \rho)^2,
\end{align*}

where the last equality uses that $O$, $\rho_j$, $\rho_k$ are Hermitian. The second moment can be further decomposed using independence of $\rhohat_i$, $\rhohat_j$, and $\rhohat_k$.

\begin{align*}
\E[\Tr( O\rhohat_i \rhohat_j ) \Tr( O\rhohat_k \rhohat_j )]
&= \Tr \left( (O \otimes O) \, \E[\rhohat_i \rhohat_j \otimes \rhohat_k \rhohat_j] \right) \\
&= \Tr \left( (O \otimes O) \, (\E[\rhohat_i] \otimes \E[\rhohat_k]) \, \E[\rhohat_j \otimes \rhohat_j] \right)
\end{align*}

Using $\E[\rhohat_i] = \E[\rhohat_k] = \rho$ and <Ref>, we have

\begin{align*}
(\E[\rhohat_i] \otimes \E[\rhohat_k]) \E[\rhohat_j \otimes \rhohat_j]
&= (\rho \rho) \left( II + \rho I + I \rho \right)( W_{(1\,2)} - 2\sym^{(2)}/(d+2) ) \\
&= 3 (\rho \rho) ( W_{(1\,2)} - 2\sym^{(2)}/(d+2) ),
\end{align*}

where we have once again use the purity of $\rho$. Plugging everything in, we get

\begin{align*}
\Cov( \Tr( O\rhohat_i \rhohat_j ), \Tr( O\rhohat_j \rhohat_k ) ) &= 3\Tr( O^{\otimes 2} \rho^{\otimes 2} W_{(1\,2)}) - \frac{6}{d+2} \Tr( O^{\otimes 2} \rho^{\otimes 2} \sym^{(2)} ) - \Tr(O \rho)^2 \\
&= 3 \Tr(O \rho)^2 - \frac{6}{d+2} \Tr(O \rho)^2 - \Tr(O \rho)^2 \\
&\leq 2 \Tr(O \rho)^2 \leq 2 \| O \|^2.
\end{align*}

For all distinct $i, j, k$,

$$\Cov( \Tr( O\rhohat_i \rhohat_j ), \Tr( O\rhohat_k \rhohat_j ) ) \leq 2 \Tr(O^2 \rho) \leq 2 \| O \|^2.$$

The proof is similar to that of <Ref>. We expand the covariance as

\Cov(\Tr(O \rhohat_i \rhohat_j), \Tr(O \rhohat_k \rhohat_j)) = \E[\Tr(O \rhohat_i \rhohat_j) \Tr(O \rhohat_j \rhohat_k)] - \Tr(O \rho)^2

and compute the second moment term using independence:

\begin{align*}
\E[\Tr(O \rhohat_i \rhohat_j) \Tr(O \rhohat_j \rhohat_k)]
&= \Tr \big((O \otimes O) (\E[\rhohat_i] \otimes I) \E[ \rhohat_j \otimes \rhohat_j] (I \otimes \E[\rhohat_k])] \big) \\
&= \Tr ( (O \otimes O) (\rho I)  ( II + \rho I + I\rho ) ( W_{(1\,2)} - 2\sym^{(2)}/(d+2) ) (I \rho) )
\end{align*}

For the $W_{(1 \, 2)}$ term, we get

\Tr( O^{\otimes 2} (\rho I)(II + I\rho + \rho I) W_{(1\,2)} (I \rho)] = \Tr( O^{\otimes 2} (2 \rho I + \rho \rho) W_{(1\,2)} ) = 2 \Tr(O^2 \rho) + \Tr(O \rho)^2.

The $\sym^{(2)}$ term subtracts a positive quantity, so we drop it to get an upper bound on covariance.

\Cov(\Tr(O \rhohat_i \rhohat_j), \Tr(O \rhohat_k \rhohat_j)) \leq 2 \Tr(O^2 \rho) + \Tr(O \rho)^2 - \Tr(O \rho)^2 = 2 \Tr(O^2 \rho) \leq 2 \| O \|^2.

We now turn to the covariance terms $\Cov(\Tr(O \rhohat_i \rhohat_j), \Tr(O \rhohat_k \rhohat_{\ell}))$ which share two indices, i.e., $(|\{ i, j \} \cap \{ k, \ell \}| = 2)$. Once again, there are two subcases: the order is swapped ($i = \ell$ and $j = k$); or the order is the same ($i = j$ and $k = \ell$).

In both cases, their are terms of the covariance that are proportional to $\Tr(O)$, but interestingly, we cannot assume $O$ is traceless as we have earlier. This is due to the fact that the $\hat Y$ estimator does not necessarily have trace 1. Nevertheless, it will turn out that this cannot affect the overall covariance of $\hat Y$ too much, as we will show in the following corollaries.

For all distinct $i, j$,

$$\Cov( \Tr( O\rhohat_i \rhohat_j ), \Tr( O\rhohat_j \rhohat_i ) ) \leq d \Tr(O^2) + 6 \sqrt{d \Tr(O^2)} + \| O \|^2$$

Using independence, we expand the covariance as

\Cov( \Tr( O\rhohat_i \rhohat_j ), \Tr( O\rhohat_j \rhohat_i ) ) = \Tr( (O \otimes O) \E[ \rhohat_i \otimes \rhohat_i] \E[ \rhohat_j \otimes \rhohat_j] ) - \Tr(O \rho)^2.

Using <Ref>, we get an expression for the second moment terms:

\begin{align*}
\E[ \rhohat_i \otimes \rhohat_i] \E[ \rhohat_j \otimes \rhohat_j]
&= ( II + \rho I + I \rho)^2 ( W_{(1\,2)} - 2\sym^{(2)}/(d+2) )^2 \\
&=\left( II + 3 I \rho + 3 \rho I + 2 \rho \rho \right) ( W_{(1)(2)} - 4\sym^{(2)}/(d+2) + 4\sym^{(2)}/(d+2)^2 ) \\
&=\left( II + 3 I \rho + 3 \rho I + 2 \rho \rho \right) ( W_{(1)(2)} - 4(d+1)\sym^{(2)}/(d+2)^2).
\end{align*}

For the $W_{(1)(2)}$ term in $\Tr( (O \otimes O) \E[ \rhohat_i \otimes \rhohat_i] \E[ \rhohat_j \otimes \rhohat_j])$, we get

\Tr( O^{\otimes 2} (II + 3 I \rho + 3 \rho I + 2 \rho \rho) W_{(1)(2)}) = \Tr(O)^2 + 6\Tr(O) + 2 \Tr(O \rho)^2.

For the $\sym^{(2)}$ term, we get (ignoring the scalar factor)

\Tr( O^{\otimes 2} (II + 3 I \rho + 3 \rho I + 2 \rho \rho) \sym^{(2)}) = \frac{\Tr(O^2) + \Tr(O)^2 + 6 \Tr(O) + 6\Tr(O^2 \rho) + 4 \Tr(O\rho)^2}{2}.

Therefore, the covariance $\Cov( \Tr( O\rhohat_i \rhohat_j ), \Tr( O\rhohat_j \rhohat_i ) )$ is

\begin{align*}
&= \Tr(O)^2 + 6\Tr(O) + \Tr(O \rho)^2 - \frac{2(d+1)}{(d+2)^2} \left( \Tr(O^2) + \Tr(O)^2 + 6 \Tr(O) + 6\Tr(O^2 \rho) + 4 \Tr(O\rho)^2 \right) \\
&\le \Tr(O)^2 + 6|\Tr(O)| + \Tr(O \rho)^2 \\
&\le d \Tr(O^2) + 6 \sqrt{d \Tr(O^2)} + \| O \|^2
\end{align*}

where the first inequality drops the negative terms and the last line comes from Cauchy-Schwarz.

For all distinct $i, j$,

$$\Cov( \Tr( O\rhohat_i \rhohat_j ), \Tr( O\rhohat_i \rhohat_j ) ) = (d+2) \Tr(O^2) + (3d-2) \| O \|^2.$$

We expand the covariance as usual and obtain the following.

\Cov( \Tr( O\rhohat_i \rhohat_j ), \Tr( O\rhohat_i \rhohat_j ) ) = \Tr( (O \otimes O) \E[ \rhohat_i \rhohat_j \otimes \rhohat_j \rhohat_i] ) - \Tr(O \rho)^2

It is not so easy to decompose this into $\E[ \rhohat_i \otimes \rhohat_i ]$ and $\E[ \rhohat_j \otimes \rhohat_j]$. We introduce a third qudit, and use the following identity:

\rhohat_i \rhohat_j \otimes \rhohat_j \rhohat_i = \Tr_{3}( (I \otimes \rhohat_i \otimes \rhohat_i) (\rhohat_j \otimes \rhohat_j \otimes I) W_{(1\,3)(2)} ).

So, plugging back into the expectation, we get

\begin{align*}
\Tr( (O \otimes O) \E[ \rhohat_i \rhohat_j \otimes \rhohat_j \rhohat_i] ) &=  \Tr( (O \otimes O \otimes I) (I \otimes \E[\rhohat_i \otimes \rhohat_i]) (\E[\rhohat_j \otimes \rhohat_j] \otimes I) W_{(1\,3)(2)} ).
\end{align*}

Once again, we use <Ref> to compute

\begin{align*}
I \otimes \E[\rhohat_i \otimes \rhohat_i] &= (III + I\rho I + II\rho)
\left( \tfrac{d+1}{d+2} W_{(1)(2\,3)} - \tfrac{1}{d+2} W_{(1)(2)(3)} \right), \\
\E[\rhohat_j \otimes \rhohat_j] \otimes I &= (III + \rho II + I\rho I)\left( \tfrac{d+1}{d+2} W_{(1\,2)(3)} - \tfrac{1}{d+2} W_{(1)(2)(3)} \right).
\end{align*}

Each expectation is a difference of two $W_{\pi}$ permutation terms. Therefore, computing the product of the two expectations, we get four $W_{\pi}$ terms. It will turn out that the dominant one is the following, where we have taken the $W_{(1)(2\,3)}$ term for $\rhohat_i$ and the $W_{(1\,2)(3)}$ term from $\rhohat_j$ (to visualize the largest contribution from this term, we refer to tensor network picture in <Ref>). Dropping the scalar factor $(d+1)^2/(d+2)^2$ for now, we get

The dominant term in Corollary <ref>, $\Tr( (O \otimes O \otimes I) \, III \, W_{(1)(2\,3)} \, III \, W_{(1\,2)(3)} W_{(1\,3)(2)}) = d \Tr(O^2)$.

\begin{align*}
&\quad \Tr( (O \otimes O \otimes I) (III + I\rho I + II\rho) W_{(1)(2\,3)} (III + \rho II + I\rho I) W_{(1\,2)(3)} W_{(1\,3)(2)} ) \\
&= \Tr( (O \otimes O \otimes I) (III + I\rho I + II\rho)(III + \rho II + II \rho) W_{(1)(2\,3)}  W_{(1\,2)(3)} W_{(1\,3)(2)} ) \\
&= \Tr( (O \otimes O \otimes I) (III + 3II \rho + I \rho I + \rho II + \rho I \rho + I \rho \rho + \rho \rho I) W_{(1\,2)(3)} ) \\
&= \Tr(O^2)(\Tr(I)+3\Tr(\rho)) + \Tr(O^2 \rho)(2\Tr(I) + 2\Tr(\rho)) + \Tr(O \rho)^2 \Tr(I) \\
&=(d+3) \Tr(O^2) + 2(d+1)\Tr(O^2 \rho) + d \Tr(O \rho)^2.
\end{align*}

For completeness, let's also compute the other 3 terms (also without their scalar factors):

\begin{align*}
&\quad \Tr( (O \otimes O \otimes I) (III + I\rho I + II\rho) W_{(1)(2)(3)} (III + \rho II + I\rho I) W_{(1)(2)(3)} W_{(1\,3)(2)} ) \\
&= \Tr( (O \otimes O \otimes I) (III + \rho II + 3I\rho I + II\rho + \rho I \rho + \rho \rho I + I\rho \rho) W_{(1\,3)(2)} ) \\
&= \Tr(O)^2 + 6 \Tr(O) \Tr(O \rho) + 2 \Tr(O \rho)^2 \\
&\le d \Tr(O^2) + 6 \sqrt{d}\Tr(O^2) + 2 \Tr(O \rho)^2 \le 6 d \Tr(O^2) + 2 \Tr(O \rho)^2,
\end{align*}

where inequality comes from Cauchy-Schwarz and the fact that $d \ge 2$. The final two terms turn out to be equal:

\begin{align*}
&\quad \Tr( (O \otimes O \otimes I) (III + I\rho I + II\rho) W_{(1)(2\,3)} (III + \rho II + I\rho I) W_{(1)(2)(3)} W_{(1\,3)(2)} ) \\
% &= \Tr( (O \otimes O \otimes I) (III + I\rho I + II\rho)(III + \rho II + II \rho) W_{(1)(2\,3)}  W_{(1)(2)(3)} W_{(1\,3)(2)} ) \\
&= \Tr( (O \otimes O \otimes I) (III + 3II \rho + I \rho I + \rho II + I \rho \rho  + \rho I \rho + \rho \rho I) W_{(1\,2\,3)} ) \\
&= \Tr(O^2) + 6 \Tr(O^2\rho) + 2 \Tr(O\rho)^2
\end{align*}

\begin{align*}
&\quad \Tr( (O \otimes O \otimes I) (III + I\rho I + II\rho) W_{(1)(2)(3)} (III + \rho II + I\rho I) W_{(1)(2\,3)} W_{(1\,3)(2)} ) \\
&= \Tr( (O \otimes O \otimes I)  (III + II\rho + 3I\rho I + \rho II + I\rho \rho+ \rho I \rho + \rho \rho I ) W_{(1\,2\,3)} ) \\
&= \Tr(O^2) + 6 \Tr(O^2\rho) + 2 \Tr(O\rho)^2
\end{align*}

Notice that these last two terms are non-negative, and so multiplying them by $-(d+1)/(d+2)^2$ makes them non-positive. Since we want to give an upper bound on the covariance, these terms can be dropped. Altogether, and inserting the appropriate constants, we get the following upper bound on the covariance $\Cov( \Tr( O\rhohat_i \rhohat_j ), \Tr( O\rhohat_i \rhohat_j ) )$:

\begin{align*}
=\,& \Tr( (O \otimes O \otimes I) (I \otimes \E[\rhohat_i \otimes \rhohat_i]) (\E[\rhohat_j \otimes \rhohat_j] \otimes I) W_{(1\,3)(2)} ) \\
\leq\,& \tfrac{(d+1)^2}{(d+2)^2} ((d+3) \Tr(O^2) + 2(d+1)\Tr(O^2 \rho) + d \Tr(O \rho)^2) + \tfrac{1}{(d+2)^2} (6 d \Tr(O^2) + 2 \Tr(O \rho)^2 ) - \Tr(O \rho)^2 \\
\leq\,& \frac{1}{(d+2)^2} ( (d^3 + 5d^2 + 13 d + 3) \Tr(O^2) + (3 d^3 + 7 d^2 + 3 d) \| O \|^2 ) \\
\leq\,& (d+2) \Tr(O^2) + (3d-2) \| O \|^2
\end{align*}

where we've used once again that $d \ge 2$.

Algorithm for Theorem <ref>
[1]
Quantum state $\rho^{\otimes N}$, $B$, $\epsilon$, $\delta$, $d$.
Classical shadow $\set{\hat{\rho}^{(i)}}_{i \in [k]}$.
$p \gets \bigo{\log(1/\delta)}$ Number of batches
$s \gets \mathcal O(\sqrt{Bd}/\epsilon + 1/\epsilon^2)$  Samples per batch

$N \gets p s$ Total number of samples
each batch $i = 1, \ldots, p$
$j = 1, \ldots, s$
$\psi_j^{(i)} \gets $ Measure fresh $\rho$ with $\mathcal{M}_1$
$\rhohat_j^{(i)} \gets (d+1) \psi_j^{(i)} - I$
$\rhohat^{(i)} \gets \frac{1}{s(s-1)} \sum_{j \neq k} \rhohat_j^{(i)} \rhohat_k^{(i)}$
$\set{\hat{\rho}^{(i)}}_{i \in [k]}$