# Informationally complete POVM-based shadow tomography

Atithi Acharya  Department of Physics and Astronomy, Rutgers University,
Piscataway, NJ 08854, USA Siddhartha Saha  Department of Physics and
Astronomy, Rutgers University, Piscataway, NJ 08854, USA Anirvan M. Sengupta
Department of Physics and Astronomy, Rutgers University, Piscataway, NJ 08854,
USA Center for Computational Mathematics and Center for Computational Quantum
Physics, Flatiron Institute, New York, NY 10010 USA

###### Abstract

Abstract: Recently introduced shadow tomography protocols use ‘classical
shadows’ of quantum states to predict many target functions of an unknown
quantum state. Unlike full quantum state tomography, shadow tomography does
not insist on accurate recovery of the density matrix for high rank mixed
states. Yet, such a protocol makes multiple accurate predictions with high
confidence, based on a moderate number of quantum measurements. One particular
influential algorithm, proposed by Huang, Kueng, and Preskill [10], requires
additional circuits for performing certain random unitary transformations. In
this paper, we avoid these transformations but employ an arbitrary
informationally complete POVM and show that such a procedure can compute
$k$-bit correlation functions for quantum states reliably. We also show that,
for this application, we do not need the median of means procedure of Huang et
al. Finally, we discuss the contrast between the computation of correlation
functions and fidelity of reconstruction of low rank density matrices.

POVM based shadows

††preprint: arXiv

##  I Introduction

Recent advances in quantum information processing often require characterizing
quantum states prepared during various stages of a procedure. As a result, the
problem of characterising a quantum state, more specifically, a density
matrix, from measurements on an ensemble of identical states, known as quantum
state tomography (QST), has seen a surge of interest [20, 5, 10]. One of the
key challenges is that, for $n$-qubit quantum systems, the density matrix is
of size $2^{n}\times 2^{n}$. As the number of qubits become large, inferring
the density matrix from a limited number of measurements becomes difficult.

Can we get away without fully characterizing the quantum state, but by
constructing an approximate classical description that predicts many different
functions of the state accurately? Shadow Tomography [2] precisely aims to do
this, namely, predict a power law number of observations in number of qubits,
$n$, from $O(n)$ copies of the density matrix $\rho$. This idea was taken
further by Huang et al. [10] who have constructed such a description of low
sample complexity via classical shadows ($\hat{\rho}$), related to states
without any entanglement in the appropriate basis, corresponding to each copy
of $\rho$ .

Quantum measurement requires specifying a set of Positive Operator Valued
Measures (POVMs) [14] which is a generalization of a complete set of
projection operators. The work by Huang et al. [10] involves measurements via
projection operators. Since projection operators are not informationally
complete (see Sec. II.1), Huang et al. employ a set of random unitary
transformations before taking measurements. In this work, we directly employ a
complete or overcomplete POVM system and perform shadow tomography. We also
gain some insight into how the choice of POVM affects the efficiency of the
method.

##  II Generalized Measurements

A projective measurement is described by an observable, $A$, a Hermitian
operator on the state space of the system being observed. The observable has a
spectral decomposition, $A=\sum_{a}aP_{a}$ where $P_{a}$ is the projector onto
the eigenspace of $A$ with eigenvalue $a$. The possible outcomes of the
measurement corresponding to the eigenvalues, $a$, of the observable and the
outcome probability is $p(a)=\bra{\psi}P_{a}\ket{\psi}$. Projection Valued
Measures (PVMs) are a special case of general measurements, where the
measurement operators are Hermitian and orthogonal projectors. A set of
Positive Operator Valued Measures (POVMs) [14] forms a generalization of PVMs.
The index $a$ in the POVM element $M_{a}$ refers to the measurement outcomes
that may occur in the experiment. The probability of the measurement outcome
is given by $p(a)=tr(\rho M_{a})$ and the post measurement density matrix can
be written as $\frac{K_{a}\rho K^{\dagger}_{a}}{tr(K_{a}\rho
K^{\dagger}_{a})}$, where $\\{K_{a}\\}$ are the Kraus operators [14]
corresponding to the POVM, with $K^{\dagger}_{a}K_{a}=M_{a}$. The operators
$\\{M_{a}\\}$ form a complete set of Hermitian non-negative operators. Namely,
they satisfy $Hermiticity:M_{a}=M^{\dagger}_{a}$,
$Positivity:\bra{\psi}M_{a}\ket{\psi}\geq 0$ for any vector $\ket{\psi}$ and
$Completeness:\sum_{a}M_{a}=\mathbb{I}$. Such a POVM could be thought of as a
partition of unity by non-negative operators.

###  II.1 Informational completeness

The density matrix, $\rho$, is a Hermitian and unit trace operator. If we have
a $d$-dimensional system, $\rho$ will be a complex square matrix represented
by $d^{2}-1$ real parameters. The operator space for this $d$-dimensional
operator will be however spanned by $d^{2}$ linearly independent basis
operators. Note that PVMs only have $d$ projection operators. They are capable
of providing only the diagonal elements of $\rho$ in a particular orthonormal
basis, leaving out potential entanglement-related information from the off-
diagonal elements. Thus, PVMs are examples of POVMs that are informationally
undercomplete.

If the number of outcomes $k$ satisfies $k\geq d^{2}$, and we can form exactly
$d^{2}$ linearly independent operators by linearly combining the set of POVMs,
such POVMs will be called informationally complete (IC). However, in the most
common terminology, informationally complete actually refers to the minimally
complete POVM ($k=d^{2}$). If we proceed to reconstruct the density matrix for
an informationally complete POVM, we can expand $\rho$ as

$\rho=\sum^{d^{2}-1}_{a=0}\xi_{a}M_{a}.$ (1)

If $k=d^{2}$ we have a informationally complete (minimally complete) basis
set. However, if we have $k>d^{2}$ $\implies$ it forms an informationally
overcomplete set [16]. IC POVM has been used for entanglement detection [8]
and for individual elements of the density matrix [13].

We start out by giving the example of a rather simple overcomplete set in the
single-qubit Hilbert space, $\bm{M}_{\text{Pauli}-6}$. Pauli-6 POVM has 6
outcomes
$\bm{M}_{\text{Pauli}-6}=\big{\\{}M_{0}=\frac{1}{3}\times\outerproduct{0}{0},M_{1}=\frac{1}{3}\times\outerproduct{1}{1},M_{2}=\frac{1}{3}\times\outerproduct{+}{+},M_{3}=\frac{1}{3}\times\outerproduct{-}{-},M_{4}=\frac{1}{3}\times\outerproduct{l}{l},M_{5}=\frac{1}{3}\times\outerproduct{r}{r}$
where $\\{\ket{0},\ket{1}\\}$, $\\{\ket{+},\ket{-}\\}$, and
$\\{\ket{r},\ket{l}\\}$ stand for the eigenbases of the Pauli operators
$\sigma^{z}$, $\sigma^{x}$, and $\sigma^{y}$, respectively. Experimentally, it
can be implemented directly by first randomly choosing $x$, $y$, or $z$, and
then measuring the respective Pauli operator, which justifies the $1/3$
factor. However, other probabilities will also be valid for this example of an
overcomplete POVM.

Now, let us give an example of a minimally complete POVM, the Pauli-4
POVM:$\,\,\bm{M}_{\text{Pauli}-4}=\big{\\{}M_{0}=\frac{1}{3}\times\outerproduct{0}{0},M_{1}=\frac{1}{3}\times\outerproduct{+}{+},M_{2}=\frac{1}{3}\times\outerproduct{l}{l},M_{3}=\frac{1}{3}\times\big{(}\outerproduct{1}{1}+\outerproduct{-}{-}+\outerproduct{r}{r}\big{)}\big{\\}}$.
As a sanity check for the completeness relation, one can see
$\sum_{a}M_{a}=1/3(\outerproduct{0}{0}+\outerproduct{1}{1}+\outerproduct{+}{+}+\outerproduct{-}{-}+\outerproduct{l}{l}+\outerproduct{r}{r})=\mathbb{I}$.
The experimental procedure will be similar to that of Pauli-6 POVM, with an
additional step where three different outcomes of Pauli-6 are identified as
the single element of Pauli-4, $M_{3}$. Thus, this set contains an element
which is not a rank-1 projector.
The third one is the tetrahedral POVM
$\bm{M}_{\text{tetra}}=\big{\\{}M_{a}=\frac{1}{4}(\mathbb{I}+\bm{s}_{a}\cdot\bm{\sigma})\big{\\}}_{a\in\\{0,1,2,3\\}}$,
whose outcomes correspond to sub-normalized rank-1 projectors along the
directions $\bm{s}_{0}=(0,0,1)$,
$\bm{s}_{1}=(\frac{2\sqrt{2}}{3},0,-\frac{1}{3})$,
$\bm{s}_{2}=(-\frac{\sqrt{2}}{3},\sqrt{\frac{2}{3}},-\frac{1}{3})$, and
$\bm{s}_{3}=(-\frac{\sqrt{2}}{3},-\sqrt{\frac{2}{3}},-\frac{1}{3})$ in the
Bloch sphere. Since the tetrahedron formed is regular, it forms an example of
a symmetric informationally complete (SIC) POVM. The experimental
implementation of $\bm{M}_{\text{tetra}}$ relies on Neumark’s dilation
theorem. The theorem implies that $\bm{M}_{\text{tetra}}$ can be physically
realized by coupling the system qubit to an ancillary qubit and performing a
von Neumann measurement on the two qubits (see Ref. [5, 19] for explicit
constructions).

##  III Classical Shadows with POVMs

Aaronson introduced the idea of “pretty good tomography”[1], with the focus on
predicting many observations accurately, based on $N$ copies of the density
matrix. This idea parallels the “learnability” of quantum states in a Probably
Approximately Correct (PAC) sense [21]. Proceeding along this line, he later
introduced the concept of Shadow Tomography [2], where from $N$ copies of the
density matrix $\rho$, we want to predict $L$ different linear target
functions $tr(O_{1}\rho),tr(O_{2}\rho).\ldots tr(O_{L}\rho)$ up to an additive
error less than $\epsilon$.

Huang et al. [10] build their methods on the idea of Shadow Tomography [2].
They repeatedly perform a measurement procedure, i.e. apply a random unitary
to rotate the state ($\rho\mapsto U\rho U^{\dagger}$) and perform a
computational-basis measurement. Then, after the measurement, they apply the
inverse of $U$ to the resulting computational basis state. This procedure
collapses $\rho$ to a snapshot
$U^{\dagger}|\hat{b}\rangle\\!\langle\hat{b}|U$, producing a quantum channel
$\mathcal{M}$, which depends on the ensemble of (random) unitary
transformations.

If the collection of unitaries is defined to be tomographically complete,
namely, if the condition i.e. for each $\sigma\neq\rho$, there exist
$U\in\mathcal{U}$ and $b$ such that $\langle b|U\sigma
U^{\dagger}|b\rangle\neq\langle b|U\rho U^{\dagger}|b\rangle$ is met, then
$\mathcal{M}$ — viewed as a linear map — has a unique inverse
$\mathcal{M}^{-1}$. Huang et al. [10] set

$\displaystyle\hat{\rho}=\mathcal{M}^{-1}\left(U^{\dagger}|\hat{b}\rangle\\!\langle\hat{b}|U\right)$
$\displaystyle\text{(classical shadow)}.$ (2)

Although the inverted channel $\mathcal{M}^{-1}$ is not physical (it is not
completely positive), one can still apply $\mathcal{M}^{-1}$ to the
(classically stored) measurement outcome
$U^{\dagger}|\hat{b}\rangle\\!\langle\hat{b}|U$ in a completely classical
post-processing step. Even if an individual sample of $\hat{\rho}$ is not a
density matrix, the expectation of $\hat{\rho}$’s is the original density
matrix $\rho$. One can use this property to get a good prediction of
measurements performed on $\rho$.

If, instead of working with the computational basis measurements, we decide to
use an IC POVM (Sec. II.1), we can avoid dealing with particular random
unitary ensembles. The only thing we need to make sure is that the resulting
channel $\mathcal{M}$ is invertible.

###  III.1 POVMs for the the $n$-qubit system

From single qubit POVMs $\\{M_{a}\\}$, we introduce $k^{n}$ operators by
taking tensor products and form POVMs for the $n$-qubit system:
$\bm{M}=\big{\\{}M_{a_{1}}\otimes
M_{a_{2}}\otimes..M_{a_{n}}\big{\\}}_{a_{1},\ldots a_{n}}$. The outcomes of
this measurements in this system are of the form
$\vec{\mathbf{a}}=(a_{1},a_{2},....a_{n})$. Now, we discuss how to form
shadows from such an observation.

###  III.2  A synthetic measurement channel

Let the POVM elements be diagonalised as follows:
$M_{a}=\sum_{i}\lambda^{a}_{i}\outerproduct{i,a}{i,a},\quad\forall
i,a,\lambda^{a}_{i}\geq 0$, since $M_{a}\succeq 0$. Let
$f:\mathcal{R}^{+}_{0}\to\mathcal{R}^{+}_{0}$ be a strictly monotonic function
which will be applied to the eigenvalues of the POVM elements. The function
$f$ is defined on $\mathcal{R}^{+}_{0}$ since the eigenvalues are non-
negative. The probability outcome ‘$a$’ is given as

$p(a)=tr(\rho M_{a}).$ (3)

Each time we perform a measurement and get an outcome ‘$a$’, we construct a
pure output state $\outerproduct{i,a}{i,a}$ with probability
$p(i|a)=\frac{f(\lambda^{a}_{i})}{\sum_{j}f(\lambda^{a}_{j})}$. We assume each
$M_{a}$ to be non-zero, guaranteeing that the denominator
$\sum_{j}f(\lambda^{a}_{j})\neq 0$. Although this is a synthetic channel, we
will refer to it as the measurement channel, in analogy with the case where
$\\{M_{a}\\}$ are projections.

The measurement channel, for a single qubit, can be defined as

$\tilde{\rho}=\mathcal{M}(\rho)=\sum_{a}p(a)\sum_{i}p(i|a)\outerproduct{i,a}{i,a}.$
(4)

For simplicity, in the following discussion, we consider the case where the
highest eigenvalue of each $M_{a}$ is non-degenerate. The modifications needed
for the general case are obvious. If a particular POVM element is not a rank
one projector and the function $f$ is very steeply increasing, then the
overwhelmingly likely output is $\ket{\psi_{a}}\bra{\psi_{a}}$, where
$\ket{\psi_{a}}$ is the eigenvector corresponding to the highest eigenvalue of
$M_{a}$. An example of such a function is $f(\lambda)=\lambda^{m}$ in the
large $m$ limit. In the large $m$ limit, as we perform a measurement, the
output is $\outerproduct{\psi_{a}}{\psi_{a}}$ (snapshots) with probability
$tr(\rho M_{a})$. The measurement channel can be defined as

$\tilde{\rho}=\mathcal{M}(\rho)=\sum_{a}tr(\rho
M_{a})\outerproduct{\psi_{a}}{\psi_{a}}.$ (5)

In a more general scheme, like the one mentioned in the beginning of the
subsection, $\ket{\psi_{a}}$ is a random vector chosen according to a
probability distribution. For example, in the current scheme, if the largest
eigenvalue of $M_{a}$ is degenerate, we choose any one of the corresponding
eigenvectors with equal probability.

In the formalism developed in [10], the channel and its inversion were related
to the ensemble of (random) unitary transformations (e.g. Clifford unitary
ensemble). The condition of tomographical completeness depended on the
existence of a unitary transformation in the chosen ensemble to distinguish
different density matrices [10]. However, with our reformulation of the
measurement channel, we need to use an informationally complete set POVMs
(e.g. Pauli-6, see Sec. II).

In the example of a single qubit measured using the 6 projectors coming from
the 3 Pauli matrices i.e. Pauli-6 POVM, the channel and its inverse can be
explicitly computed. Similar to the classical shadows built out of random
Pauli measurements [10], we get a depolarizing channel i.e. a channel that
contracts a pure state (lying on the surface of the Bloch sphere) towards the
‘center’ of the Bloch sphere, namely, the maximally mixed state
$\rho=\mathbb{I}_{2}/{2}$. The inverse (a non-physical map) can be computed,
which can map a point inside the Bloch ball to the outside.

Figure 1: The convex region in the figure is the set of admissible density
matrices. We schematically describe the process of forming classical shadows
from $N$ copies of $\rho$. For the $i$-th observation with outcome ‘$a$’, the
inverse of the channel, $\mathcal{M}^{-1}$, acts on the projectors
$\outerproduct{\psi_{a}}{\psi_{a}}$ to construct the shadow $\hat{\rho}^{i}$.
The sample mean of the shadows cast by $\rho$ i.e.
$\frac{1}{N}\sum_{i}\hat{\rho}^{i}$ fluctuates around the true $\rho$ and
could be outside the convex region. However, while measuring $L$ $k$-local
observables [10] $\vec{O}=(O_{1},\ldots O_{L})$, the convergence of the sample
averages $\vec{\hat{o}}(N)$ to the true expected values can be guaranteed with
a number of samples $O(\log L)$. See Theorem 1.

Multi-qubit system: For local measurements (not necessarily the depolarizing
channel), the inverse channel for the $n$-qubit system can be written as

$\mathcal{M}_{n}^{-1}=\bigotimes_{j=1}^{n}\mathcal{M}_{1}^{-1}.$ (6)

We can now reformulate the shadows with our overcomplete POVM set and its
corresponding channel. For instance, when we work with Pauli-6 POVM, we will
get

$\hat{\rho}=\bigotimes_{j=1}^{n}\mathcal{M}_{1}^{-1}(\outerproduct{\psi_{a,j}}{\psi_{a,j}})\qquad(classical\
shadow),$ (7)

where $\mathcal{M}_{1}^{-1}(X)=3X-tr(X)\mathbb{I}$ (see Sec. A.2). Note that
the $2^{n}\times 2^{n}$ matrix $\hat{\rho}$ need not be constructed
explicitly. We just need to store $\ket{\psi_{a,j}}$ for each qubit $j$.

Since the inverted channel $\mathcal{M}^{-1}$ is not physical (it is not
completely positive), the $\hat{\rho}$ in Eq. (7) need not be physical. In
other words, there is no guarantee the output of the inverse channel is
positive semidefinite. See Fig. 1 for a schematic description. We recover the
true density matrix only in expectation. However, if the shadow matrix is
forced to be positive semidefinite, we can see how the observations such as
fidelity changes (see Sec. IV.1).

###  III.3 Noisy shadow

Earlier, we defined our measurement channel, Eq. (5). However, we can also let
each of our qubits pass through a previously characterized noise channel
$\mathcal{E}_{1}$ and then take the measurements [11]. The combined channel
${M}_{\mathcal{E},1}$ is given by

$\tilde{\rho}=\mathcal{M}_{\mathcal{E},1}(\rho)=\sum_{a}tr(\mathcal{E}_{1}(\rho)M_{a})\outerproduct{\psi_{a}}{\psi_{a}}.$
(8)

We used IC POVMs to ensure that the measurement channel was invertible. As
long as the action of the noise channel $\mathcal{E}_{1}$ itself is
invertible, $\mathcal{M}_{\mathcal{E},1}$ is also invertible. We will work
with an $n$-qubit noise channel of the form
$\mathcal{E}_{n}=\bigotimes_{j=1}^{n}\mathcal{E}_{1}$. Thus, we can still
write the inverse of the new noisy measurement channel for the $n$-qubit
system in terms of the single qubit inverse shadow channel
$\mathcal{M}_{\mathcal{E},1}^{-1}$:

$\mathcal{M}_{\mathcal{E},n}^{-1}=\bigotimes_{j=1}^{n}\mathcal{M}_{\mathcal{E},1}^{-1}.$
(9)

If we choose an amplitude damping channel with damping parameter $\gamma$, one
of the Kraus operator representations can be given as

$\mathcal{E}_{AD}(\rho)=K_{0}\rho K^{\dagger}_{0}+K_{1}\rho K^{\dagger}_{1},$
(10)

where $K_{0}=\begin{bmatrix}1&&0\\\ 0&&\sqrt{(1-\gamma)}\end{bmatrix}$,
$K_{1}=\begin{bmatrix}0&&\sqrt{\gamma}\\\ 0&&0\end{bmatrix}$.

The inverse of the noisy shadow channel $\mathcal{M}_{{AD}}^{-1}(X)$ is given
in Eq. (24). Its action on $\mathbb{I}$ and $\sigma_{x,y,z}$ is given as
$\mathcal{M}_{{AD}}^{-1}(\mathbb{I})=\mathbb{I}-\frac{\gamma}{1-\gamma}\sigma_{z}$,
$\mathcal{M}_{{AD}}^{-1}(\sigma_{x,y})=\frac{3}{\sqrt{1-\gamma}}\sigma_{x,y}$
and $\mathcal{M}_{{AD}}^{-1}(\sigma_{z})=\frac{3}{1-\gamma}\sigma_{z}$. See
Sec.A.2 for a general description on the inversion of a noisy shadow channel.
Here, we will construct the shadows (noisy) with the following definition:

$\hat{\rho}=\bigotimes_{j=1}^{n}\mathcal{M}_{AD}^{-1}(\outerproduct{\psi_{a,j}}{\psi_{a,j}}).$
(11)

###  III.4 Predicting linear functions with classical shadows

Using the statistical properties of a single shadow, we can predict linear
functions in the unknown state $\rho$ as

$o=tr(O\rho)=\mathbb{E}[\hat{o}],\qquad\textrm{where}\quad\hat{o}=tr(O\hat{\rho}).$
(12)

In practice, using an array of shadows (i.e. $N$ snapshots), we can estimate
the expectation $o$. Given an array of $N$ independent classical snapshots
(each defined as in Eq. (7)) :

$\mathsf{S}(\rho;N)=\left\\{\hat{\rho}^{(1)},\hat{\rho}^{(1)},\ldots,\hat{\rho}^{(N)}\right\\}.$
(13)

The sample mean is
$\bar{o}=\frac{1}{N}\sum_{j=1}^{N}\mathrm{tr}\left(O\hat{\rho}^{(j)}\right).$
This sample mean will fluctuate around the true prediction, with
$\mathbb{E}(\bar{o})=o$.

###  III.5 The algorithm and the guarantee of performance

We want to predict the expected value of multiple $k$-local observables
$O_{1},\ldots O_{L}$ based on shadows using the two algorithms below.

Input: IC POVM with $k$ outcomes, $\rho\in\mathcal{C}^{2^{n}}$ (N copies of
the unknown density matrix)

1 Compute the measurement channel $\mathcal{M}_{1}$ and its inverse
$\mathcal{M}^{-1}_{1}$ for the chosen IC POVM. (See Sec. A.2) ;

2

3for  _$i=1,\ldots N$_ do

4 Perform measurements using the POVM elements $M_{a}$ to get outcomes
$a_{ji}\in\\{1,\ldots,k\\}$ ;

5 Construct shadows
$\hat{\rho}_{i}=\bigotimes_{j=1}^{n}\mathcal{M}_{1}^{-1}(\outerproduct{\psi_{a_{ji}}}{\psi_{a_{ji}}})$
(See Sec. III.2, A for the general version)

6 end for

7

Output: $\hat{\rho}_{1},\hat{\rho}_{2}\ldots\hat{\rho}_{N}$

8

Algorithm 1 Generating Shadows with POVMs

1

2

Input: A POVM set, N copies of unknown density matrix $\rho$, L different
$k$-local Pauli observables $O_{1},O_{2},\ldots O_{L}$ and error parameters
$\epsilon,\delta$

3

4Find bounds on the local observables $B(\\{O\\},\mathcal{M})$. (See Sec.A.3
for details).

5

6Using algorithm.1, collect
$N\geq\frac{B(\\{O\\},\mathcal{M})\log(\frac{2L}{\delta})}{2\epsilon^{2}}$
shadows.

7Compute means
$\hat{o}_{i}=\frac{1}{N}\sum_{j=1}^{N}\mathrm{tr}\left(O_{i}\hat{\rho}^{(j)}\right)$

8

Output: $\bar{o}_{1},\bar{o}_{2}\ldots\bar{o}_{L}$

9

10

Algorithm 2 Predicting many properties using mean as an estimate

The existence of the bound is guaranteed by the following theorem.

######  Theorem 1.

With
$N\geq\frac{B(\\{O\\},\mathcal{M})\log(\frac{2L}{\delta})}{2\epsilon^{2}}$
samples of $\rho$, we can predict $L$ different linear target functions
$tr(O_{1}\rho),\ tr(O_{2}\rho),\ldots,tr(O_{L}\rho)$ up to additive error
$\epsilon$ with maximum failure probability $\delta$.

The constant bound $B(\\{O\\},\mathcal{M})$ will depend on the measurement
channel $\mathcal{M}$ (which depends on the choice of POVM) and on the
operator set $\\{O\\}=\\{O_{1},O_{2},..\hskip 1.0pt,O_{L}\\}$). The important
thing is that $B(\\{O\\},\mathcal{M}$ is bounded for so called $k$-local
operators, as defined in [10].

For instance, if we choose Pauli-6, the bound is given as
$\hat{o}^{(j)}_{i}\in[-3^{k},3^{k}]$, in which case
$B(k\text{-local},\text{Pauli-}6)=4\times 9^{k}$ (See Sec.A.3). See section
III.2 for the algorithm, including the construction of the measurement
channel. In the appendix section on sample complexity (Sec. A.3), the details
of the proof is provided.

##  IV Numerical Results

For many quantum systems in Condensed Matter Physics, one of the objects of
interest is the two-point correlation function. Two-point correlators could be
efficiently estimated using classical shadows based on Pauli-6 POVM. The
predictions of two-point functions
$\langle\sigma^{Z}_{i}\sigma^{Z}_{j}\rangle$ for the GHZ states with varying
degree of noise is shown in Fig. 2.

Figure 2: Prediction of two-point correlations,
$\langle\sigma^{z}_{i}\sigma^{z}_{j}\rangle$, for noisy GHZ target states
using classical shadows for Pauli-6, with 1-standard deviation band. The
standard deviations are estimated over ten independent runs, each of which
involved $N=5000$ samples. The parameter $p$, representing the local
depolarizing noise strength, is described in Eq. (14).

We can write the action of the single qubit depolarizing noise [14] on an
arbitrary $\rho$ written in the Bloch sphere representation:

$\rho^{\prime}=(1-\frac{4}{3}p)\rho+\frac{4}{3}p\mathbb{I}.$ (14)

Applying this channel to every qubit, we generate a noise GHZ state [9] from a
pure one. The expected two-point correlations
$\langle\sigma^{Z}_{i}\sigma^{Z}_{j}\rangle$ varies as $(1-4p/3)^{2}$ with the
noise parameter $p$.
While predicting multiple $1,\ldots,L$, two-point or $k$-point correlations,
we monitor the maximum possible error among all the observables. This measure
of error is expected to go down with increasing number of samples. This
scaling, as seen in Fig. 3, gives us some idea of the appropriateness of a
POVM set for a particular task.

(a) GHZ

(b) Spin down

Figure 3: Maximum error in two-point correlators. (a) Scaling of maximum error
among all two-point correlations in 30 qubit pure GHZ state, plotted against
different number of samples for different choice of POVMs: Pauli-6, Pauli-4
and tetrahedral. (b) Scaling of maximum error for all spin down state with
Pauli-4 and Pauli-6 POVM. Pauli-4 ensures a much better scaling. See Sec. A.3
for details.

####  IV.0.1 1D Transverse Field Ising Model

We take antiferromagnetic ($J>0$ in Eq. (15) ) transverse field Ising model in
1D:

$H=J\sum_{<ij>}\sigma^{z}_{i}\sigma^{z}_{j}+h\sum_{i}\sigma^{x}.$ (15)

The quantum critical point at $h/J=1$ will be exhibited by the power-law decay
of the correlations. See Fig. 4 for results in the three regimes: critical,
ordered and paramagnetic. The exact numerical correlations are plotted using
the matrix product representations of the ground states. [15]. In [5] and
[12], POVM-based measurements, followed by a neural-network-centric approach
for constructing the ground state, and computing the resulting two-point
correlations were presented for the same system.

(a)  $J=h=1$

(b)  $J>h$, $J=1$ and $h=0.5$

(c)  $J<h$, $J=0.5$ and $h=1$

Figure 4: Two-point functions $\langle\sigma^{z}_{0}\sigma^{z}_{i}\rangle$ for
ground states of antiferromagnetic 1D tranverse field Ising model using
Pauli-6 POVM based shadows and the true value, as computed using matrix
product states. The correlations are plotted against the lattice separation.
The lattice size is 30 and the number of samples used is 5000. (a) Critical
($J=h$) antiferromagnetic 1D TFIM, showing signatures of power-law
correlation. (b) Ordered state ($J>h$), where correlations saturates with
increasing lattice separation. (c) The paramagnetic state ($J<h)$, displaying
exponential decay of the correlations.

####  IV.0.2 1D Disordered Heisenberg Model

The Hamiltonian for the 1D disordered Heisenberg model is given by,

$\hat{H}=-\frac{1}{2}\sum_{j=1}^{N}(J^{x}_{j}\sigma_{j}^{x}\sigma_{j+1}^{x}+J^{y}_{j}\sigma_{j}^{y}\sigma_{j+1}^{y}+J^{z}_{j}\sigma_{j}^{z}\sigma_{j+1}^{z}+h\sigma_{j}^{z}).$
(16)

The properties of spin-$\frac{1}{2}$ antiferromagnetic chains with various
types of random exchange coupling has been studied in an exact decimation
renormalization-group (strong-disorder) schemes, some of which involve
generalization or modifications of the scheme introduced by Dasgupta and Ma
[7]. The numerical studies done by R.N. Bhatt and P.A. Lee [4] indicate that
the system could be in a random-singlet phase. In such a phase, each spin is
paired with another spin that may be far away on the lattice. We perform exact
diagonalization, obtain the ground state and then compute two-point quantum
correlations. The 2d plot of the correlation matrix will also inform us about
the locations of the singlet formations in the chain. We can also reconstruct
these behavior of a ground-state corresponding to one particular disorder
realization of the XXZ-Heisenberg model Eq. (16) ($J_{x}=J_{y}=2J_{z}$, $h=0$)
with sufficient number of shadows. See Fig. 5, where the singlet formations
are indicated by the schematics drawn on the axes of the matrix visualization
plots and the results from the two methods are compared.

(a)  Exact Diagonalization

(b) Reconstruction with Shadows

Figure 5: Two-point functions $\langle\sigma^{z}_{i}\sigma^{z}_{j}\rangle$ for
ground states of disordered 1D Heisenberg spin chain with length=10 and open
boundary conditions. (a) Exact diagonalization results. (b) Results from using
Pauli-6 POVM based shadows using 5000 samples.

###  IV.1 Exploring quantum fidelity

In our approach to construct shadows using local POVMs, we ensure prediction
of local observables. However, we can also explore non-local observables such
as fidelity. Using sample mean as an estimator, we can construct a hypothesis
state ($\sigma$):

$\sigma=\frac{1}{N}\sum^{N}_{i}\hat{\rho}.$ (17)

When our target state is pure, we can rewrite quantum fidelity as a linear
prediction with our target observable given as $O=\outerproduct{\psi}{\psi}$.
Starting from this definition of quantum fidelity i.e.
$F_{Q}(\rho,\sigma)=(tr(\sqrt{\sqrt{\rho}\sigma\sqrt{\rho}}))^{2}$, using
$\rho=\outerproduct{\psi}{\psi}$ we get
$(tr\sqrt{(\bra{\psi}\sigma\ket{\psi}\outerproduct{\psi}{\psi})})^{2}$.
Further simplification of the fidelity gives us

$F_{Q}(\rho,\sigma)=\bra{\psi}\sigma\ket{\psi}(tr(\sqrt{\outerproduct{\psi}{\psi}}))^{2}=\bra{\psi}\sigma\ket{\psi}=tr(\sigma
O).$

Figure 6: Quantum fidelity predicted for the pure GHZ state using sample mean
of shadows constructed on $10^{4}$ samples. The shaded regions are the
standard deviation over ten independent runs. The inset shows the scaling of
the variance of fidelity which grows exponentially with number of qubits.

The measure $tr(\sigma O)$ is equivalent to fidelity, only when the latter is
defined, i.e. when $\sigma\succeq 0$. That property is likely to hold only
when the number of samples is large. We expect $tr(\sigma O)$ to fluctuate
around its mean value 1, as seen in Fig. 6, even when the typical $\sigma$ is
not a physical state, meaning it is not positive semidefinite. Also, the
fluctuation around this mean keeps on growing exponentially with the number of
qubits (see Fig. 6). This growth cannot be dealt with even by the median of
means (MoM) procedure [10] within the shadow formalism. Numerical computations
using MoM also show no advantage over sample means here.

Hence, we need a procedure to find the ‘closest’ physical state to $\sigma$.
The trace condition $tr(\hat{\sigma})=1$ ensures that once $\hat{\sigma}\nsucc
0$, some of the eigenvalues will be greater than 1 to compensate for the
negative eigenvalues. Thus, we cannot just throw away the negative
eigenvalues, as would be done for projecting a Hermitian matrix to the space
of positive semidefinite matrices.

We define the the convex set of physical states to be $C=\\{\rho|\rho\succeq
0,tr(\rho)=1\\}$. Our nonlinear projection to $C$ is

$\Pi_{C}(\sigma)=\operatorname*{arg\,min}_{\rho\in C}tr((\rho-\sigma)^{2}).$
(18)

We achieve this by diagonalizing $\sigma$, projecting the eigenvalues
$\lambda_{i}$ of $\sigma$ onto a canonical simplex
$\Delta=\\{(\lambda^{p}_{1},\dots,\lambda^{p}_{D})\mid\lambda^{p}_{i}\geq
0,\sum_{i=1}^{D}\lambda^{p}_{i}=1\\}$, using the recipe from Ref. [23], while
leaving the eigenvectors untouched. Here, $D=2^{n}$ where $n$ is the total
number of qubits. The projected state is a biased estimator. We can hope that
the price paid by accepting some bias comes with the benefit of reduced
variance. This expectation seems to be born out in Fig. 7. However, as number
of qubits increase, the bias itself reduces fidelity. To compensate this
effect, we need larger sample sizes ($N$). Fig. 7 shows all these trends.

Figure 7: Quantum fidelity of the projected shadows (onto the physical
positive definite space) with the noiseless GHZ state. As we increase the
number of samples, from $10^{3}$ to $10^{4}$, the quantum fidelity improves.
The shaded regions indicate 1-standard deviation bands, estimated over ten
independent runs.

##  V Discussions

We provide an approach to predict expectations of local observables without
having to apply random unitary transformations, which sometimes require
complex circuits of its own, and can become a practical bottleneck. We show
that this can rather be done using an IC POVM. For illustrations, we show
faithful reconstruction properties of low energy states coming from different
many body Hamiltonians relevant to near-term applications of quantum devices.
When we have additional information about the possible noisy channels we also
adapt the shadow channel as a composition of the noise channel and the
measurement channel. The invertibility becomes straightforward in the proposed
framework. We also comment on why the mean as an estimator is sufficient
throughout our discussion. And as long as we are dealing with local
observables, we can provide efficient sample complexity using Hoeffding’s
inequality directly.

We provided instances where the choice of POVM impacts the sample complexity
for predicting 2-point correlators in certain quantum states for fixed maximum
error. We noted that the different POVMs work better for different states. It
is an exciting endeavour to understand which sets of POVM would be ideal for
different classes of quantum states and observables.

Although, an exploration, we attempt to reconstruct fidelity using the locally
built shadows and show that we cannot benefit from median of means as an
estimator, since variance of fidelity becomes exponential in number of qubits.
Additionally, when presented with few samples we raise the issue of unphysical
i.e. not positive semidefinite $\hat{\rho}$ and then provide a projection
tecnique, similar to [17], to estimate fidelity. Unfortunately, the estimator
no longer remains unbiased. Addressing this issue would require methods to
deal with non-local observables.

We did not provide an effective analog of the global Clifford unitary
transformation-based method in [10]. There has been work which provides
description of global alternatives using stabilizer states [17]. Whether there
can be a scheme based on such states that is competitive with the classical
shadows method [10] remains to be seen.

The use of generalized measurement to unambiguously discriminate non-
orthogonal states with lower failure probability is well known [3, 14, 6].
Efficient prediction of expectations of local observables combined with the
generalized measurement scheme to obtain the shadows can be used as an optimal
framework in the discrimination of non-orthogonal states. In the future, it is
a promising direction of exploration.

## Acknowledgement

We would like to thank Shagesh Sridharan, James Stokes, Miles Stoudenmire for
insightful discussions.

##  Appendix A Appendix

###  A.1 The Measurement Channel for Pauli-6

We can take the simple rank-1 Pauli-6 POVMs to see the action of a measurement
channel:

$\tilde{\rho}=\mathcal{M}(\rho)=\sum_{a}tr\big{(}\frac{1}{2}(\mathbb{I}+\mathbf{r}.\pmb{\sigma})M_{a}\big{)}\outerproduct{\psi_{a}}{\psi_{a}},$
(19)

where we use the Bloch representation
$\rho=\tfrac{1}{2}(\mathbb{I}+\mathbf{r}.\pmb{\sigma})$.

The contribution of the first two POVM elements of Pauli-6 only gets
contribution from $\mathbb{I}$ and $r_{z}\sigma_{z}$, generating

$tr(\frac{1}{2}(\mathbb{I}+r_{z}\sigma_{z})M_{0})\outerproduct{0}{0}+tr(\frac{1}{2}(\mathbb{I}+r_{z}\sigma_{z})M_{1})\outerproduct{1}{1}.$

Using $M_{0}=\frac{1}{3}\times\outerproduct{0}{0}$ and
$M_{1}=\frac{1}{3}\times\outerproduct{1}{1}$ this expression becomes:
$\frac{1}{6}\mathbb{I}+\frac{1}{3}r_{z}\sigma_{z}$ Following similar steps for
pairs $M_{2},M_{3}$ and $M_{4},M_{5}$, we get that:

$\mathcal{M}(\rho)=\mathcal{M}\big{(}\frac{1}{2}(\mathbb{I}+\mathbf{r}.\pmb{\sigma})\big{)}=\frac{1}{2}(\mathbb{I}+\frac{1}{3}\mathbf{r}.\pmb{\sigma}),$

making $\mathcal{M}$ a depolarizing channel.

###  A.2 Inverse of the measurement channel

Given any single qubit channel, the inverse can be easily computed using the
Bloch-sphere representation. We can write any 2 dimensional (single qubit)
quantum operation ($X$) as $X=(x_{0}\mathbb{I}+\vec{r}.\vec{\sigma})$. Any
arbitrary trace-preserving quantum operation is given as
$\mathcal{E}(X)=(x_{0}\mathbb{I}+\vec{r}^{\hskip 2.0pt\prime}.\vec{\sigma})$.
The map $\vec{r}\xrightarrow{\mathcal{E}}\vec{r}^{\hskip 2.0pt\prime}$ is
equivalent to,

$\vec{r}^{\hskip 2.0pt\prime}=T\vec{r}+x_{0}\vec{c},\quad
T_{i,j}=\frac{1}{2}tr(\sigma_{i}\mathcal{E}(\sigma_{j})).$ (20)

The components of displacement ($\vec{c}$) is given as
$c_{i}=\frac{1}{2}tr(\sigma_{i}\mathcal{E}(\mathbb{I}))$. The affine map
between the Bloch sphere and itself is given by $T$, and its meaning is
understood better by doing a singular value decomposition i.e.
$T=O_{1}DO^{\prime}_{2}$ where $O_{1},O_{2}$ are orthogonal matrices. The
singular values capture the deformation of the Bloch sphere about its
principal axes. A superoperator $\hat{T}_{4\times 4}$ can be defined as

$\displaystyle\hat{T}\begin{bmatrix}x_{0}\\\
\vec{r}\end{bmatrix}=\begin{bmatrix}x_{0}\\\ \vec{r}^{\hskip
2.0pt\prime}\end{bmatrix},\quad\hat{T}_{4\times 4}$
$\displaystyle=\begin{bmatrix}1&&0\\\ \vec{c}&&\small{T_{3\times
3}}\end{bmatrix}.$ (21)

Computing inverse of the channel is equivalent to writing $(x_{0},\vec{r})$
from $(x_{0},\vec{r}^{\hskip 2.0pt\prime})$ i.e. computing $\hat{T}^{-1}$.

In the main text, the Pauli measurement channel
($\mathcal{E}=\mathcal{M}_{1}$) turns out to be a depolarizing channel, and
its inverse that acts on the local qubit is given as

$\mathcal{M}_{1}^{-1}(X)=3X-tr(X)\mathbb{I}.$

We take a more general example following our definition of a measurement
channel:

$\tilde{\rho}=\mathcal{M}(\rho)=\sum_{a}tr(\rho
M_{a})\outerproduct{\psi_{a}}{\psi_{a}}.$

If a particular POVM element is not rank one, $\ket{\psi_{a}}$ can be taken as
the eigenvector corresponding to the highest eigenvalue of $M_{a}$. For
Pauli-4, except for
$M_{3}=\frac{1}{3}(\outerproduct{1}{1}+\outerproduct{-}{-}+\outerproduct{r}{r})=\mathbb{I}$,
all other elements are rank-1. Since $M_{3}$ is rank-2, when the outcome is
$3$, we take the eigenvector $\ket{t}$ corresponding to eigenvalue
$\frac{1}{2}(1+\frac{1}{\sqrt{3}})$ instead of the other corresponding to
$\frac{1}{2}(1-\frac{1}{\sqrt{3}})$. Rewriting Eq. (A.2), we get

$\begin{split}\tilde{\rho}=\mathcal{M}_{1}(\rho)=tr(\rho
M_{0})\outerproduct{0}{0}+tr(\rho M_{1})\outerproduct{+}{+}\\\ +tr(\rho
M_{2})\outerproduct{l}{l}+M_{3})\outerproduct{t}{t}.\end{split}$ (22)

The inverse of the channel can be written as

$\begin{split}\mathcal{M}_{1}^{-1}(X)=6X-(\frac{x_{s}}{\sqrt{3}}-x_{0}(\sqrt{3}-1))(\sum_{i}\sigma_{i})-5x_{0}\mathbb{I}\\\
\quad x_{0}=\frac{tr(X)}{2},\hskip
3.0ptx_{s}=\frac{3(tr(X(\sum_{i}\sigma_{i}))+(\sqrt{3}-1)x_{0})}{\sqrt{3}+1}.\end{split}$
(23)

When we are working with a known noise channel $\mathcal{E}$, the inverse is
given as
$\hat{T}_{\mathcal{M}_{\mathcal{E}}}^{-1}=\hat{T}_{\mathcal{E}}^{-1}\hat{T}_{\mathcal{M}}^{-1}$.
If we choose an amplitude damping channel with a damping parameter $\gamma$,
the inverse can be given as

$\begin{split}\mathcal{M}_{\mathcal{E}}^{-1}(X)=\frac{3(1-\frac{1}{2}tr(X\sigma_{z}))}{\sqrt{(1-\gamma)}}X+\frac{3tr(\sigma_{z}X)}{2(1-\gamma)}X\\\
+(\frac{1}{2}-\frac{3}{2\sqrt{(1-\gamma)}})tr(X)\mathbb{I}+\frac{\gamma
tr(X)}{2(\gamma-1)}\sigma_{z}.\end{split}$ (24)

###  A.3 Sample complexity

####  A.3.1 Variance of the Estimate for a Single Observable

Given an array of $N$ independent, classical snapshots (each defined as Eq.
(7)) :

$\mathsf{S}(\rho;N)=\left\\{\hat{\rho}^{(1)},\hat{\rho}^{(1)},\ldots,\hat{\rho}^{(N)}\right\\}.$
(25)

The sample mean is
$\hat{o}=\frac{1}{N}\sum_{j=1}^{N}\mathrm{tr}\left(O\hat{\rho}^{(j)}\right).$
The bound on probability of deviation of the sample mean is given by
Chebyshev’s inequality:

$Pr(\absolutevalue{\hat{o}-\mathbb{E}\left[\hat{o}\right]}\geq\epsilon)\leq\frac{Var[\hat{o}]}{\epsilon^{2}}$
(26)

$\mathbb{E}\left[\hat{o}\right]=\mathrm{tr}\left(O\rho\right)$ where $\rho$ is
the true density matrix. Fluctuations of $\hat{o}$ around this desired
expectation are controlled by the variance.
$Var[\hat{o}]=\frac{1}{N}Var[tr(O\hat{\rho}^{(j)})]=\frac{Var[o^{(j)}]}{N}$.
However, since the classical shadows are unit trace by construction, the
variance depends only on the trace-less part of the observable i.e.
$O_{0}=O-\frac{tr(O)}{2^{n}}\mathbb{I}$. The minimum number of samples needed
to assure a maximum failure probability ($\delta$) using Eq. (26) is

$N\geq\frac{Var[o^{(j)}]}{\epsilon^{2}\delta}.$ (27)

####  A.3.2 Dependence on POVM

Given a measurement channel and an observable, we can bound the variance of
its estimator, using familiar maneuvers with superoperators [10],

$\displaystyle
Var[o^{(j)}]=\mathbb{E}((o^{(j)})^{2})-(\mathbb{E}(o^{(j)}))^{2}\leq\mathbb{E}((o^{(j)})^{2})$
$\displaystyle=\sum_{a_{1},.,a_{n}}Pr(a_{1}..a_{n})\bra{a_{1},.,a_{n}}[\mathcal{M}^{-1}_{n}]^{\dagger}(O_{0})\ket{a_{1},.,a_{n}}^{2}$
$\displaystyle\text{where}\quad Pr(a_{1},..,a_{n})=tr(\rho M^{a_{1}}\otimes
M^{a_{1}}..\otimes M^{a_{n}}).$

We broadly define a $k$-local Pauli-observable as an operator which acts
nontrivially only on $k$ qubits. Traceless $k$ local operators can be
expressed as linear conbination of tensor products of indentity matrices and
$k$ or less Pauli matrices. Hence, we need to focus only on special class of
$k$-local operators. Denoting $P_{i}$ as one of the Pauli matrices acting on
the $i$th qubit, we focus on of tenor products like $O_{0}=P_{1}\otimes
P_{2}\otimes..\otimes P_{k}\otimes\mathbb{I}^{\otimes(n-k)}$, where, without
loss of generality, we assume that the operator acts non-trivially on only the
first $k$ qubits.

For Pauli-6 POVM, the inverse of the measurement channel is a self-adjoint
map, and thus one can verify its action as:

$\displaystyle[\mathcal{M}_{1}^{-1}]^{\dagger}(P_{\alpha})=\mathcal{M}_{1}^{-1}(P_{\alpha})=3P_{\alpha},$

where $P_{\alpha}$ denotes a Pauli matrix and
$[\mathcal{M}_{1}^{-1}]^{\dagger}(\mathbb{I})=\mathcal{M}_{1}^{-1}(\mathbb{I})=\mathbb{I}$.
Given a $k$-local observable, we can further compute the bound on variance:

$\displaystyle
Var[o^{(j)}]\leq\sum_{a_{1},.,a_{n}}Pr(a_{1}..a_{n})\prod^{k}_{i=1}\bra{a_{i}}3P_{i}\ket{a_{i}}^{2}$
$\displaystyle=\sum_{a_{1},.,a_{n}}tr(\rho M^{a_{1}}\otimes M^{a_{2}}..\otimes
M^{a_{n}})\prod^{k}_{i=1}\bra{a_{i}}3P_{i}\ket{a_{i}}^{2}$
$\displaystyle=tr[\rho\sum_{a_{1},.,a_{n}}(M^{a_{1}}\otimes M^{a_{1}}..\otimes
M^{a_{n}}\prod^{k}_{i=1}\bra{a_{i}}3P_{i}\ket{a_{i}})^{2}]$
$\displaystyle=tr[\rho\bigotimes^{k}_{i=1}\underbrace{\sum_{a_{i}}M^{a_{i}}\bra{a_{i}}3P_{i}\ket{a_{i}}^{2}}_{3\mathbb{I}}\bigotimes\mathbb{I}^{\otimes(n-k)}]=3^{k}.$

Now, we take up Pauli-4 POVM. One can verify the action of
$[\mathcal{M}_{1}^{-1}]^{\dagger}$ as:

$\displaystyle[\mathcal{M}_{1}^{-1}]^{\dagger}(P_{\alpha})=(2-\sqrt{3})\mathbb{I}+(3+\sqrt{3})P_{\alpha}-\sum_{\beta\neq\alpha}(3-\sqrt{3})P_{\beta},$

where $P_{\alpha}$ denotes a Pauli matrix. Using the fact that
$\mathcal{M}_{1}^{-1}$ is a trace preserving map, one can say its adjoint has
to be unital $i.e.$ $[\mathcal{M}_{1}^{-1}]^{\dagger}(\mathbb{I})=\mathbb{I}$.
Given a $k$-local Pauli-observable, one can again compute the bound on
variance:

$\displaystyle
Var[o^{(j)}]\leq\sum_{a_{1},.,a_{n}}Pr(a_{1}..a_{n})\prod^{k}_{i=1}\bra{a_{i}}[\mathcal{M}_{1}^{-1}]^{\dagger}(P_{i})\ket{a_{i}}^{2}$
$\displaystyle=\sum_{a_{1},.,a_{n}}tr(\rho M^{a_{1}}\otimes M^{a_{2}}..\otimes
M^{a_{n}})\prod^{k}_{i=1}\bra{a_{i}}[\mathcal{M}_{1}^{-1}]^{\dagger}(P_{i})\ket{a_{i}}^{2}$
$\displaystyle=tr[\rho\sum_{a_{1},.,a_{n}}(M^{a_{1}}\otimes M^{a_{1}}..\otimes
M^{a_{n}}\prod^{k}_{i=1}\bra{a_{i}}[\mathcal{M}_{1}^{-1}]^{\dagger}(P_{i})\ket{a_{i}}^{2}]$
$\displaystyle=tr[\rho\bigotimes^{k}_{i=1}\underbrace{\sum_{a_{i}}M^{a_{i}}\bra{a_{i}}[\mathcal{M}_{1}^{-1}]^{\dagger}(P_{i})\ket{a_{i}}^{2}}_{5\mathbb{I}+4P_{i}}\bigotimes\mathbb{I}^{\otimes(n-k)}]$
$\displaystyle=tr[\rho\bigotimes^{k}_{i=1}{(5\mathbb{I}+4P_{i})}\bigotimes\mathbb{I}^{\otimes(n-k)}].$

Clearly, the above bound on variance is dependent on the state $\rho$, unlike
the bound we obtained using Pauli-6 POVM. Since $\rho$ is a density matrix and
the operator $5\mathbb{I}+4P_{i}$ is a PSD operator, one gets the minimum
value for the bound when $\rho$ is of the form:

$\rho=\left(\bigotimes^{k}_{i=1}{\outerproduct{p_{i}}{p_{i}}}\right)\bigotimes\tilde{\rho}_{n-k},$
(28)

where $\outerproduct{p_{i}}{p_{i}}$ is the projector into the eigenvector
corresponding to the lowest eigenvalue of the operator $5\mathbb{I}+4P_{i}$,
and $\tilde{\rho}_{n-k}$ is a valid density matrix in the Hilbert space of
$n-k$ qubits on which the $k$-local Pauli-observable acts trivially. For the
above $\rho$, it is simple to verify that the value of the variance bound is 1
(independent of $k$). Thus, for example, if the unknown state $\rho$ is the
all spin down state, then Pauli-4 POVM works better than Pauli-6 POVM in
predicting two-point correlators $\langle\sigma^{Z}_{i}\sigma^{Z}_{j}\rangle$,
since the variance is higher in the latter.

####  A.3.3 Improved Bound Using Hoeffding’s Inequality

Furthermore, we can use Hoeffding’s inequality to provide theoretical bounds
when we are dealing with $k$-local Pauli observable, since we are working with
bounded random variables. If
$\hat{o}=\frac{1}{N}\sum_{j=1}^{N}\mathrm{tr}\left(O\hat{\rho}^{(j)}\right)$,
$\hat{o}^{(j)}\in[a,b]$ for all $j$, where $-\infty<a\leq b\leq\infty$, we can
write,

$Pr(\absolutevalue{\hat{o}-\mathbb{E}\left[\hat{o}\right]}\geq\epsilon)\leq
2e^{(\frac{-2N\large{\epsilon}^{2}}{(b-a)^{2}})}.$ (29)

The minimum number of samples needed to assure a maximum failure probability
($\delta$) among all the observables using Eq. (31) is

$N\geq\log(\frac{2}{\delta})\frac{(b-a)^{2}}{2\epsilon^{2}}.$ (30)

$B(\\{O\\},\mathcal{M})=(b-a)^{2}$ depends on the locality of the observable
and the maximum eigenvalue $\lambda_{max}$ of the inverse channel acting on
the observable. The bound on random variable $\hat{o}^{(j)}$ can be found as
the range of the Rayleigh quotient of the inverse of the measurement channel,
acting on the observable over all possible states. For instance, if we choose
Pauli-6, the bounds can be shown to lie within
$\hat{o}^{(j)}\in[-3^{k},3^{k}]$ in which case
$B(k\text{-local},\text{Pauli-}6)=4\times 9^{k}$. Using the action of
$[\mathcal{M}^{-1}]^{\dagger}$, one can verify that the value of the random
variable $tr([\mathcal{M}^{-1}]^{\dagger}(P_{i})\outerproduct{a_{i}}{a_{i}})$
belongs to the set $\\{5,-1\\}$ for any Pauli matrix $P_{i}$ when
$\outerproduct{a_{i}}{a_{i}})$ is the inferred state for Pauli-4. Thus, the
random variable $\hat{o}$ is contained in the range $\\{-5^{k},5^{k}\\}$ which
is exponential on the locality rather than the number of qubits.

####  A.3.4 The Guarantee of Performance for Multiple Observables

If we have $L$ different $k$-local Pauli observables
$\hat{o}_{1},\ldots,\hat{o}_{i},\ldots,\hat{o}_{L}$ with the sample mean
corresponding to the observable $i$ defined as
$\hat{o}_{i}=\frac{1}{N}\sum_{j=1}^{N}\mathrm{tr}\left(O_{i}\hat{\rho}^{(j)}\right)$.
If
$\hat{o}_{i}=\frac{1}{N}\sum_{j=1}^{N}\mathrm{tr}\left(O_{i}\hat{\rho}^{(j)}\right)$,
$\hat{o}^{(j)}_{i}\in[a,b]$ for all $j$, where $-\infty<a\leq b<\infty$, we
can combine the union bound with Hoeffding’s inequality to write

$Pr(\max\limits_{\begin{subarray}{c}1\leq i\leq
L\end{subarray}}\absolutevalue{\hat{o}_{i}-\mathbb{E}\left[\hat{o}_{i}\right]}\geq\epsilon)\leq
2Le^{(\frac{-2N\large{\epsilon}^{2}}{(b-a)^{2}})}.$ (31)

The minimum number of samples needed to assure a maximum failure probability
($\delta$) among all the observables using Eq. (31) is

$N\geq\log(\frac{2L}{\delta})\frac{(b-a)^{2}}{2\epsilon^{2}}.$ (32)

The scaling is logarithmic in the number of observables $L$, instead of linear
behavior we get using Chebyshev’s inequality. We do not need to use MoM
procedure [10], which would have been necessary if we were dealing with
estimate distributions with long tails (unlike the bounded estimates for
$k$-local Pauli observables).

###  A.4 Numerical computations

The computations for GHZ states has been carried out using MPS matrix product
state (MPS) representations for noiseless states and matrix product operator
(MPO) representations for the noisy states. The details on the simulation of
mixed states using MPOs has been shown in [22, 5]. The data-sets corresponding
to the ground states of spin Hamiltonians such as the Transverse-field Ising
model have been generated using the density matrix renormalization group
(DMRG). The library used is mpnum (a matrix product represenation library for
Python). [18]. Given a particular spin model, with the Hamiltonian expressed
as a MPO, the DMRG algorithm attempts to find the optimal MPS with the lowest
energy. However, for the visualization of correlations in the disordered 1d
Heisenberg spin chain, the computations are made without using the DMRG
framework. This is because the number of sites were low for this particular
numerical experiment.

## References

  * [1] Scott Aaronson.  The learnability of quantum states.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 463(2088):3089–3114, Sep 2007.
  * [2] Scott Aaronson.  Shadow tomography of quantum states, 2018.
  * [3] Stephen M. Barnett and Sarah Croke.  Quantum state discrimination.  Adv. Opt. Photon., 1(2):238–278, Apr 2009.
  * [4] R. N. Bhatt and P. A. Lee.  Scaling studies of highly disordered spin-½ antiferromagnetic systems.  Phys. Rev. Lett., 48:344–347, Feb 1982.
  * [5] Juan Carrasquilla, Giacomo Torlai, Roger G. Melko, and Leandro Aolita.  Reconstructing quantum states with generative models.  Nature Machine Intelligence, 1(3):155–161, Mar 2019.
  * [6] Anthony Chefles.  Quantum state discrimination.  Contemporary Physics, 41(6):401–424, Nov 2000.
  * [7] Chandan Dasgupta and Shang-Keng Ma.  Low-temperature properties of the random Heisenberg antiferromagnetic chain.  Phys. Rev. B, 22(3):1305–1319, August 1980.
  * [8] Aleksandra Dimić and Borivoje Dakić.  Single-copy entanglement detection.  npj Quantum Information, 4(1):1–8, 2018.
  * [9] Daniel M. Greenberger, Michael A. Horne, and Anton Zeilinger.  Going beyond Bell’s theorem.  eprint arXiv: 0712.0921, 2007.
  * [10] Hsin-Yuan Huang, Richard Kueng, and John Preskill.  Predicting many properties of a quantum system from very few measurements.  Nature Physics, Jun 2020.
  * [11] Dax Enshan Koh and Sabee Grewal.  Classical shadows with noise.  e-print arXiv:2011.11580, 2020.
  * [12] Ilia A. Luchnikov, Alexander Ryzhov, Pieter-Jan Stas, Sergey N. Filippov, and Henni Ouerdane.  Variational autoencoder reconstruction of complex many-body physics.  Entropy, 21(11):1091, Nov 2019.
  * [13] Joshua Morris and Borivoje Dakić.  Selective quantum state tomography.  arXiv preprint arXiv:1909.05880, 2019.
  * [14] M. A Nielsen and I. L. Chuang.  Quantum Computation and Quantum Information: 10th Anniversary Edition.  Cambridge University Press, 10 edition, 2011.
  * [15] Román Orús.  A practical introduction to tensor networks: Matrix product states and projected entangled pair states.  Annals of Physics, 349:117–158, Oct 2014.
  * [16] Joseph M. Renes, Robin Blume-Kohout, A. J. Scott, and Carlton M. Caves.  Symmetric informationally complete quantum measurements.  Journal of Mathematical Physics, 45(6):2171–2180, Jun 2004.
  * [17] G.I. Struchalin, Ya. A. Zagorovskii, E.V. Kovlakov, S.S. Straupe, and S.P. Kulik.  Experimental estimation of quantum state properties from classical shadows.  PRX Quantum, 2(1), Jan 2021.
  * [18] Daniel Suess and Milan Holzäpfel.  mpnum: A matrix product representation library for Python.  Journal of Open Source Software, 2(20):465, 2017.
  * [19] Gelo Noel M. Tabia.  Experimental scheme for qubit and qutrit symmetric informationally complete positive operator-valued measurements using multiport devices.  Phys. Rev. A, 86:062107, Dec 2012.
  * [20] Giacomo Torlai, Guglielmo Mazzola, Juan Carrasquilla, Matthias Troyer, Roger Melko, and Giuseppe Carleo.  Neural-network quantum state tomography.  Nature Physics, 14(5):447–450, Feb 2018.
  * [21] L. G. Valiant.  A theory of the learnable.  Communications of the ACM, 1984.
  * [22] F. Verstraete, J. J. García-Ripoll, and J. I. Cirac.  Matrix product density operators: Simulation of finite-temperature and dissipative systems.  Phys. Rev. Lett., 93:207204, Nov 2004.
  * [23] Weiran Wang and Miguel Á. Carreira-Perpiñán.  Projection onto the probability simplex: An efficient algorithm with a simple proof, and an application.  CoRR, abs/1309.1541, 2013.