# Provably efficient machine learning for quantum many-body problems

Hsin-Yuan Huang  Institute for Quantum Information and Matter and
Department of Computing and Mathematical Sciences, Caltech, Pasadena, CA, USA
Richard Kueng  Institute for Integrated Circuits, Johannes Kepler University
Linz, Austria  Giacomo Torlai  AWS Center for Quantum Computing, Pasadena, CA,
USA  Victor V. Albert  Joint Center for Quantum Information and Computer
Science, National Institute of Standards and Technology and University of
Maryland, College Park, MD, USA  John Preskill  Institute for Quantum
Information and Matter and
Department of Computing and Mathematical Sciences, Caltech, Pasadena, CA, USA
AWS Center for Quantum Computing, Pasadena, CA, USA

(August 27, 2024)

###### Abstract

Classical machine learning (ML) provides a potentially powerful approach to
solving challenging quantum many-body problems in physics and chemistry.
However, the advantages of ML over more traditional methods have not been
firmly established. In this work, we prove that classical ML algorithms can
efficiently predict ground state properties of gapped Hamiltonians in finite
spatial dimensions, after learning from data obtained by measuring other
Hamiltonians in the same quantum phase of matter. In contrast, under widely
accepted complexity theory assumptions, classical algorithms that do not learn
from data cannot achieve the same guarantee. We also prove that classical ML
algorithms can efficiently classify a wide range of quantum phases of matter.
Our arguments are based on the concept of a classical shadow, a succinct
classical description of a many-body quantum state that can be constructed in
feasible quantum experiments and be used to predict many properties of the
state. Extensive numerical experiments corroborate our theoretical results in
a variety of scenarios, including Rydberg atom systems, 2D random Heisenberg
models, symmetry-protected topological phases, and topologically ordered
phases.

##  I Introduction

Figure 1:  (a) Efficient quantum-to-classical conversion. The classical shadow
of a quantum state, constructed by measuring very few copies of the state, can
be used to predict many properties of the state with a rigorous performance
guarantee. (b) Predicting ground state properties. After training on data
obtained in quantum experiments, a classical ML model predicts a classical
representation of the ground state $\rho(x)$ of the Hamiltonian $H(x)$ for
parameters $x$ spanning the entire phase. This representation yields estimates
of the properties of $\rho(x)$, avoiding the need to run exhaustive classical
computations or quantum experiments. (c) Classifying quantum phases. After
training, a classical ML receives a classical representation of a quantum
state, and predicts the phase from which the state was drawn. (d) Training
data. For predicting ground states, the classical ML receives a classical
representation of $\rho(x)$ for each value of $x$ sampled during training. For
predicting quantum phases of matter, the training data consists of classical
representations of quantum states accompanied by labels identifying the phase
to which each state belongs.

Solving quantum many-body problems, such as finding ground states of quantum
systems, has far-reaching consequences for physics, materials science, and
chemistry. While classical computers have facilitated many profound advances
in science and technology, they often struggle to solve such problems.
Powerful methods, such as density functional theory HohenbergKohn ; NobelKohn
, quantum Monte Carlo CEPERLEY555 ; SandvikSSE ; becca_sorella_2017 and
density-matrix renormalization group DMRG1 ; DMRG2 , have enabled solutions to
certain restricted instances of many-body problems, but many general classes
of problems remain outside the reach of even the most advanced classical
algorithms.

Scalable fault-tolerant quantum computers will be able to solve a broad array
of quantum problems, but are unlikely to be available for years to come.
Meanwhile, how can we best exploit our powerful classical computers to advance
our understanding of complex quantum systems? Recently, classical machine
learning (ML) techniques have been adapted to investigate problems in quantum
many-body physics CarleoRMP ; APXReview , with promising results dassarma2017
; carrasquilla2017nature ; Carleo_2017 ; torlai_learning_2016 ; Nomura2017 ;
evert2017nature ; leiwang2016 ; gilmer2017neural ; torlai_Tomo ;
vargas2018extrapolating ; schutt2019unifying ; Glasser2018 ;
rodriguez2019identifying ; qiao2020orbnet ; choo_fermionicnqs2020 ;
kawai2020predicting ; moreno2020deep ; Kottmann2021 . So far these approaches
are mostly heuristic, reflecting the general paucity of rigorous theory in ML.
While shown to be effective in some intermediate-size experiments Bohrdt2018 ;
Rem2018 ; torlai_rydberg19 , these methods are generally not backed by
convincing theoretical arguments to ensure good performance, particularly for
problem instances where traditional classical algorithms falter.

In general, simulating quantum many-body physics is hard for classical
computers, because accurately describing an $n$-qubit quantum system may
require an amount of classical data that is exponential in $n$. In prior work,
some of us addressed this bottleneck using _classical shadows_ — succinct
classical descriptions of quantum many-body states that can be used to
accurately predict a wide range of properties with rigorous performance
guarantees huang2020predicting ; paini2019approximate . Furthermore, this
quantum-to-classical conversion technique can be readily implemented in
various existing quantum experiments Struchalin2021Shadows ;
Elben2020Entanglement ; choi2021emergent . Classical shadows open new
opportunities for addressing quantum problems using classical methods such as
ML. In this paper, we build on the classical shadow formalism and devise
polynomial-time classical ML algorithms for quantum many-body problems which
are supported by rigorous theory.

We consider two applications of classical ML, indicated in Figure 1. The first
application we examine is learning to predict classical representations of
quantum many-body ground states. We consider a family of Hamiltonians, where
the Hamiltonian $H(x)$ depends smoothly on $m$ real parameters (denoted by
$x$). The ML algorithm is trained on a set of training data consisting of
sampled values of $x$, each accompanied by the corresponding classical shadow
for the ground state $\rho(x)$ of $H(x)$. This training data could be obtained
from either classical simulations or quantum experiments. During the
prediction phase, the ML algorithm predicts a classical representation of
$\rho(x)$ for new values of $x$ different from those in the training data.
Ground state properties can then be estimated using the predicted classical
representation.

This learning algorithm is efficient, provided that the ground state
properties to be predicted do not vary too rapidly as a function of $x$.
Indeed, sufficient upper bounds on the gradient can be derived for any family
of gapped geometrically-local Hamiltonians in any finite spatial dimension, if
the property of interest is the expectation value of a sum of few-body
observables. The conclusion is that any such property can be predicted with a
small average error, where the amount of training data and the classical
computation time are polynomial in $m$ and linear in the system size.
Furthermore, we show that classical algorithms that do not learn from data
cannot provide the same rigorous guarantee without violating widely accepted
complexity-theoretic conjectures. This is a manifestation of the advantage of
ML algorithms with data over those without data huang2020power .

If the training data are obtained from quantum experiments, then one might
choose to learn about properties of $\rho(x)$ for a new input $x$ by
conducting new experiments rather than by using the classical ML to generalize
from the training data. However, ML could be far more convenient in some
cases, especially when changing some parameters may even require costly re-
engineering of the entire experiment. ML algorithms open up the possibility of
efficiently and accurately predicting properties of quantum states that are
extremely challenging to prepare and measure in the laboratory.

Classical ML could be used to generalize from training data that are obtained
from either quantum experiments or classical simulations; the same rigorous
performance guarantees apply in either case. Even if the training data are
generated classically, it could be more efficient and more accurate to use ML
to predict properties for new values of the input $x$, rather than doing new
simulations which could be computationally very demanding and of unverified
reliability. Promising insights into quantum many-body physics are already
being obtained using classical ML based on classical simulation data
dassarma2017 ; Nomura2017 ; Carleo_2017 ; zhang2017machine ;
zhang2020interpreting ; gilmer2017neural ; vargas2018extrapolating ;
schutt2019unifying ; qiao2020orbnet ; choo_fermionicnqs2020 ;
kawai2020predicting . Our rigorous analysis identifies general conditions that
guarantee the success of classical ML models, and elucidates the advantages of
classical ML models over non-ML algorithms. These results enhance the
prospects for interpretable ML techniques ribeiro2016should ;
arrieta2020explainable ; zhang2020interpreting to further shed light on
quantum many-body physics.

In the second application we examine, the goal is to classify quantum states
of matter into phases Read2012 in a supervised learning scenario. Suppose that
during training we are provided with sample quantum states which carry labels
indicating whether each state belongs to phase $A$ or phase $B$. Our goal is
to predict the phase label for new quantum states that were not encountered
during training. We assume that, during both the learning and prediction
stages, each quantum state is represented by its classical shadow, which could
be obtained either from a classical computation or from an experiment on a
quantum device. The classical ML, then, trains on labeled classical shadows,
and learns to predict labels for new classical shadows.

We assume that the $A$ and $B$ phases can be distinguished by a nonlinear
function of marginal density operators of subsystems of constant size. This
assumption is reasonable because we expect the phase to be revealed in
subsystems which are larger than the correlation length, but independent of
the total system size. We show that if such a function exists, a classical ML
can learn to distinguish the phases using an amount of training data and
classical processing which are polynomial in the system size. We do not need
to know anything about this nonlinear function in advance, apart from its
existence.

In what follows, we review the classical shadow formalism huang2020predicting
, and use this formalism to derive rigorous guarantees for ML algorithms in
predicting ground state properties and classifying quantum phases of matter.
We also describe numerical experiments in a wide range of physical systems to
support our theoretical results.

##  II Constructing efficient classical representations of quantum systems

We begin with an overview of the randomized measurement toolbox Ohliger_2013 ;
VanEnk2012 ; Elben2018 ; evans2019scalable ; huang2020predicting ;
paini2019approximate , relegating further details to Appendix A. We
approximate an $n$-qubit quantum state $\rho$ by performing randomized single-
qubit Pauli measurements on $T$ copies of $\rho$. That is, we measure every
qubit of the unknown quantum state $\rho$ in a random Pauli basis $X,Y$ or $Z$
to yield a measurement outcome of $\pm 1$. Collapse of the wavefunction
implies that this measurement procedure transforms $\rho$ into a random pure
product state $\ket{s^{(t)}}=\bigotimes_{i=1}^{n}\ket{s^{(t)}_{i}}$, where
$\ket{s^{(t)}_{i}}\in\\{\ket{0},\ket{1},\ket{+},\ket{-},\ket{\mathrm{i}+},\ket{\mathrm{i}-}\\}$
are eigenstates of the selected Pauli matrices. Performing one randomized
measurement grants us classical access to one such snapshot. Performing a
total of $T$ randomized measurements grants us access to an entire collection
$S_{T}(\rho)=\big{\\{}|s_{i}^{(t)}\rangle:\;i\in\\{1,\ldots,n\\},\;t\in\\{1,\ldots,T\\}\big{\\}}$.
Each element is a highly structured single-qubit pure state, and there are
$nT$ of them in total. So, $3nT$ bits suffice to store the entire collection
in classical memory. The randomized measurements can be performed in actual
physical experiments or through classical simulations. Resulting data can then
be used to approximate the underlying $n$-qubit state $\rho$:

$\rho\approx\sigma_{T}(\rho)=\frac{1}{T}\sum_{t=1}^{T}\sigma_{1}^{(t)}\otimes\dots\otimes\sigma_{n}^{(t)}\quad\text{where}\quad\sigma_{i}^{(t)}=3\lvert
s_{i}^{(t)}\rangle\\!\langle s_{i}^{(t)}\rvert-\mathbb{I},$ (1)

and $\mathbb{I}$ denotes the $2\times 2$ identity matrix. This _classical
shadow_ representation huang2020predicting ; paini2019approximate exactly
reproduces the global density matrix in the limit $T\to\infty$, but
$T=\mathcal{O}(\mathrm{const}^{r}\log(n)/\epsilon^{2})$ already provides an
$\epsilon$-accurate approximation of _all_ reduced $r$-body density matrices
(in trace distance). This, in turn, implies that we can use $\sigma_{T}(\rho)$
to predict any function that depends on only reduced density matrices, such as
expectation values of (sums of) local observables and (sums of) entanglement
entropies of small subsystems. Classical storage and postprocessing cost also
remain tractable in this regime. To summarize, the classical shadow formalism
equips us with an efficient quantum-to-classical converter that allows
classical machines to efficiently and reliably estimate subsystem properties
of any quantum state $\rho$.

##  III Predicting ground states of quantum many-body systems

We consider the task of predicting ground state representations of quantum
many-body Hamiltonians in finite spatial dimensions. Suppose that a family of
geometrically local, $n$-qubit Hamiltonians $\\{H(x):\;x\in[-1,1]^{m}\\}$ is
parametrized by a classical variable $x$. That is, $H(x)$ smoothly maps a
bounded $m$-dimensional vector $x$ (parametrization) to a Hermitian matrix of
size $2^{n}\times 2^{n}$ ($n$-qubit Hamiltonian). We do not impose any
additional structure on this mapping; in particular, we do not assume
knowledge about how the physical Hamiltonian depends on the parameterization.
The goal is to learn a model $\hat{\sigma}(x)$ that can predict properties of
the ground state $\rho(x)$ associated with Hamiltonian. This problem arises in
many practical scenarios. Suppose diligent experimental effort has produced
experimental data for ground state properties of various physical systems. We
would like to use this data to train an ML model that predicts ground state
representations of hitherto unexplored physical systems.

###  III.1 An ML algorithm with rigorous guarantee

We will prove that a classical ML algorithm can predict classical
representations of ground states after training on data belonging to the same
quantum phase of matter. Formally, we consider a smooth family of Hamiltonians
$H(x)$ with a constant spectral gap. During the training phase of the ML
algorithm, many values of $x$ are randomly sampled, and for each sampled $x$,
the classical shadow of the corresponding ground state $\rho(x)$ of $H(x)$ is
provided, either by classical simulations or quantum experiments. The full
training data of size $N$ is given by
$\big{\\{}x_{\ell}\rightarrow\sigma_{T}(\rho(x_{\ell}))\big{\\}}_{\ell=1}^{N}$,
where $T$ is the number of randomized measurements in the construction of the
classical shadows at each value of $x_{\ell}$.

We train classical ML models using the size-$N$ training data, such that when
given the input $x_{\ell}$, the ML model can produce a classical
representation $\hat{\sigma}(x)$ that approximates
$\sigma_{T}(\rho(x_{\ell}))$. During prediction, the classical ML produces
$\hat{\sigma}(x)$ for new values of $x$ different from those in the training
data. While $\hat{\sigma}(x)$ and $\sigma_{T}(\rho(x_{\ell}))$ classically
represent exponentially large density matrices, the training and prediction
can be done efficiently on a classical computer using various existing
classical ML models, such as neural networks with large hidden layers
jacot2018neural ; li2019enhanced ; du2019graph ; neuraltangents2020 and kernel
methods cortes1995support ; CC01a . In particular, the predicted output of the
trained classical ML models can be written as the extrapolation of the
training data using a learned metric $\kappa(x,x_{\ell})\in\mathbb{R}$,

$\hat{\sigma}(x)=\frac{1}{N}\sum_{\ell=1}^{N}\kappa(x,x_{\ell})\sigma_{T}(\rho(x_{\ell})).$
(2)

For example, prediction using a trained neural network with large hidden
layers jacot2018neural is equivalent to using the metric
$\kappa(x,x_{\ell})=\sum_{\ell^{\prime}=1}^{N}f^{\mathrm{(NTK)}}(x,x_{\ell^{\prime}})(F^{-1})_{\ell^{\prime}\ell}$,
where $f^{\mathrm{(NTK)}}(x,x^{\prime})$ is the neural tangent kernel
jacot2018neural corresponding to the neural network and
$F_{\ell^{\prime}\ell}=f^{\mathrm{(NTK)}}(x_{\ell^{\prime}},x_{\ell})$; see
Appendix C for more discussion. The ground state properties are then estimated
using these predicted classical representations $\hat{\sigma}(x)$.
Specifically, $f_{O}(x)=\mathrm{tr}\left(O\rho(x)\right)$ can be predicted
efficiently whenever $O$ is a sum of few-body operators.

To derive a provable guarantee, we consider the simple metric
$\kappa(x,x_{\ell})=\sum_{k\in\mathbb{Z}^{m},\left\lVert
k\right\rVert_{2}\leq\Lambda}\cos(\pi k\cdot(x-x_{\ell}))$ with cutoff
$\Lambda$, which we refer to as the $l_{2}$-Dirichlet kernel. We prove that
the prediction will be accurate and efficient if the function $f_{O}(x)$ does
not vary too rapidly when $x$ changes in any direction. Indeed, sufficient
upper bounds on the gradient magnitude of $f_{O}(x)$ can be derived using
quasi-adiabatic continuation hastings2005quasiadiabatic ;
bachmann2012automorphic .

Under the $l_{2}$-Dirichlet kernel, the classical ML model is equivalent to
learning a truncated Fourier series to approximate the function $f_{O}(x)$.
The parameter $\Lambda$ is a cutoff for the wavenumber $k$ that depends on
(upper bounds on) the gradient of $f_{O}(x)$. Using statistical analysis, one
can guarantee that
$\operatorname*{{\mathbb{E}}}_{x}|\operatorname{tr}(O\hat{\sigma}(x))-f_{O}(x)|^{2}\leq\epsilon$
as long as the amount of training data obeys $N=m^{\mathcal{O}(1/\epsilon)}$
in the $m\to\infty$ limit. The conclusion is that any such $f_{O}(x)$ can be
predicted with a small _constant_ average error, where the amount of training
data and the classical computation time are polynomial in $m$ and at most
linear in the system size $n$. Moreover, the training data need only contain a
single classical shadow snapshot at each point $x_{\ell}$ in the parameter
space (i.e., $T=1$). An informal statement of the theorem is given below; we
explain the proof strategy in Appendix E, and provide more details in Appendix
F.  We also discuss how one could generalize the proof to long-range
interacting systems, electronic Hamiltonians, and other settings in Appendix
E.2.

######  Theorem 1 (Learning to predict ground state representations;
informal).

For any smooth family of Hamiltonians $\\{H(x):\;x\in[-1,1]^{m}\\}$ in a
finite spatial dimension with a constant spectral gap, the classical machine
learning algorithm can learn to predict a classical representation of the
ground state $\rho(x)$ of $H(x)$ that approximates few-body reduced density
matrices up to a constant error $\epsilon$ when averaged over $x$. The
required training data size $N$ and computation time are polynomial in $m$ and
linear in the system size $n$.

Though formally “efficient” in the sense that $N$ scales polynomially with $m$
for any fixed approximation error $\epsilon$, the required amount of training
data scales badly with $\epsilon$. This unfortunate scaling is not a
shortcoming of the considered ML algorithm, but a necessary feature. In
Appendix G, we show that the data size and time complexity cannot be improved
further without making stronger assumptions about the class of gapped local
Hamiltonians. However, in cases of practical interest, the Hamiltonian may
obey restrictions such as translational invariance or graph structure that can
be exploited to obtain better results. Incorporating these restrictions can be
achieved by using a suitable $\kappa(x,x_{\ell})$, such as one that
corresponds to a large-width convolutional neural network li2019enhanced or a
graph neural network du2019graph . Rigorously establishing that neural-
network-based ML algorithms can achieve improved prediction performance and
efficiency for particular classes of Hamiltonians is a goal for future work.

###  III.2 Computational hardness for non-ML algorithms

In the following proposition, we show that a classical algorithm that does not
learn from data cannot achieve the same guarantee in estimating ground state
properties without violating the widely believed conjecture that NP-complete
problems cannot be solved in randomized polynomial time. This proposition is a
corollary of standard complexity-theoretic results Lichtenstein1982PlanarFA ;
VALIANT198685 . See Appendix H for the detailed statement and proof.

######  Proposition 1 (Informal).

Consider a randomized polynomial-time classical algorithm $\mathcal{A}$ that
does not learn from data. Suppose for any smooth family of two-dimensional
Hamiltonians $\\{H(x):\;x\in[-1,1]^{m}\\}$ with a constant spectral gap,
$\mathcal{A}$ can efficiently compute expectation values of one-body
observables in the ground state $\rho(x)$ of $H(x)$ up to a constant error
when averaged over $x$. Then there is a randomized classical algorithm that
can solve NP-complete problems in polynomial time.

It is instructive to observe that a classical ML algorithm with access to data
can perform tasks that cannot be achieved by classical algorithms which do not
have access to data. This phenomenon is studied in huang2020power , where it
is shown that the complexity class defined by classical algorithms that can
learn from data is strictly larger than the class of classical algorithms that
do not learn from data. (The data can be regarded as a restricted form of
randomized advice string.) We caution that obtaining the data to train the
classical ML model could be challenging. However, if we focus only on data
that could be efficiently generated by quantum-mechanical processes, it is
still possible that a classical ML that learns from data could be more
powerful than classical computers. In Appendix H we present a contrived family
of Hamiltonians that establishes this claim, based on the (classical)
computational hardness of factoring.

##  IV Classifying quantum phases of matter

Classifying quantum phases of matter is another important application of
machine learning to physics. We will consider this classification problem in
the case where quantum states are succinctly represented by their classical
shadows. For simplicity, we consider the classification of two phases (denoted
$A$ and $B$), but the analysis naturally generalizes to classifying any number
of phases.

###  IV.1 ML algorithms

We envision training a classical ML with classical shadows, where each
classical shadow carries a label $y$ indicating whether it represents a
quantum state $\rho$ from phase $A$ $(y(\rho)=1)$ or phase $B$ $(y(\rho)=-1)$.
We want to show that a suitably chosen classical ML can learn to efficiently
predict the phase for new classical shadows beyond those encountered during
training. Following a strategy which is standard in learning theory, we
consider a classical ML that maps each classical shadow to a corresponding
feature vector in a high-dimensional feature space, and then attempts to find
a hyperplane that separates feature vectors in the $A$ phase from feature
vectors in the $B$ phase. The learning is efficient if the geometry of the
feature space is efficiently computable, and if the feature map is
sufficiently expressive. Thus, our task is to construct a feature map with the
desired properties.

In the simpler task of classifying symmetry-breaking phases, there is
typically a local order parameter $O=\sum_{i}O_{i}$ given as a sum of $r$-body
observables for some $r>0$ that satisfies

$\operatorname{tr}(O\rho)\geq 1,\forall\rho\in\mbox{phase
$A$},\quad\quad\operatorname{tr}(O\rho)\leq-1,\forall\rho\in\mbox{phase $B$}.$
(3)

Under this criterion, the classification function may be chosen to be
$y(\rho)=\operatorname*{sign}(\operatorname{tr}(O\rho))$. Hence, classifying
symmetry-breaking phases can be achieved by finding a hyperplane that
separates the two phases in the high-dimensional feature space that subsumes
all $r$-body reduced density matrices of the quantum state $\rho$. The feature
vector consisting of all $r$-body reduced density matrices of the quantum
state $\rho$ can be accurately reconstructed from the classical shadow
representation $S_{T}(\rho)$ when $T$ is sufficiently large.

Finding a suitable choice of hyperplane in the feature space can be cast as a
convex optimization problem known as the soft-margin support vector machine,
discussed in more detail in Appendix J.1. With a sufficient amount of training
data, the hyperplane found by the classical ML model will generalize so the
phase $y(\rho)$ can be predicted accurately for a previously unseen quantum
state $\rho$. The classical ML is not merely a black box; it exhibits the
order parameter (encoded by the hyperplane), guiding physicists toward a
deeper understanding of the phase structure.

For more exotic quantum phases of matter, such as topologically ordered
phases, the above classical ML model no longer suffices. The topological phase
of a state is invariant under a constant-depth quantum circuit, and a phase
containing the product state $\ket{0}^{\otimes n}$ is called the trivial
phase. Using these notions, we can prove that no observable — not even one
that acts on the entire system — can be used to distinguish between two
topological phases. The proof, given in Appendix I, uses the observation that
random single-qubit unitaries can confuse any global or local order parameter.

######  Proposition 2.

Consider two distinct topological phases $A$ and $B$ (one of the phases could
be the trivial phase). No observable $O$ exists such that

$\operatorname{tr}(O\rho)>0,\forall\rho\in\mbox{phase
$A$},\quad\quad\operatorname{tr}(O\rho)\leq 0,\forall\rho\in\mbox{phase $B$}.$
(4)

While this proposition implies that no linear function
$\operatorname{tr}(O\rho)$ can be used to classify topologically ordered
phases, it does not exclude nonlinear functions, such as quadratic functions
$\operatorname{tr}(O\rho\otimes\rho)$, degree-$d$ polynomials
$\operatorname{tr}(O\rho^{\otimes d})$ and more general analytic functions.
For example, it is known that the topological entanglement entropy
kitaev2006topological ; levin2006detecting , a nonlinear function of $\rho$,
can be used to classify a wide variety of topologically ordered phases. For
this purpose, it suffices to consider a subsystem whose size is large compared
to the correlation length of the state, but is independent of the total size
of the system. The correlation length in the ground state of a local
Hamiltonian increases when the spectral gap between the ground state and the
first excited state becomes smaller hastings2006spectral . On the other hand,
a linear function on the full system will fail even with constant correlation
length.

To learn nonlinear functions, we need a more expressive ML model. For this
purpose we devise a powerful feature map that takes the classical shadow
$S_{T}(\rho)$ of the quantum state $\rho$ to a feature vector that includes
arbitrarily-large $r$-body reduced density matrices, as well as an
arbitrarily-high-degree polynomial expansion,

$\phi^{\mathrm{(shadow)}}(S_{T}(\rho))=\lim_{D,R\to\infty}\bigoplus_{d=0}^{D}\sqrt{\frac{\tau^{d}}{d!}}\left(\bigoplus_{r=0}^{R}\sqrt{\frac{1}{r!}\left(\frac{\gamma}{n}\right)^{r}}\bigoplus_{i_{1}=1}^{n}\ldots\bigoplus_{i_{r}=1}^{n}\mathrm{vec}\left[\frac{1}{T}\sum_{t=1}^{T}\bigotimes_{\ell=1}^{r}\sigma_{i_{\ell}}^{(t)}\right]\right)^{\otimes
d},$ (5)

where $\tau,\gamma>0$ are hyper-parameters. The direct sum
$\bigoplus_{r=0}^{R}$ is a concatenation of all $r$-body reduced density
matrices, and the other direct sum $\bigoplus_{d=0}^{D}$ subsumes all
degree-$d$ polynomial expansions. The computational cost of finding a
hyperplane in feature space that separates the training data into two classes
is dominated by the cost of computing inner products between feature vectors.
The inner product
$\langle\phi^{\mathrm{(shadow)}}(S_{T}(\rho)),\phi^{\mathrm{(shadow)}}(S_{T}(\tilde{\rho}))\rangle$
can be analytically computed by reorganizing the direct sums, writing it as a
double series, and wrapping both series into an exponential, which gives

$k^{\mathrm{(shadow)}}\left(S_{T}(\rho),S_{T}(\tilde{\rho})\right)=\exp\left(\frac{\tau}{T^{2}}\sum_{t,t^{\prime}=1}^{T}\exp\left(\frac{\gamma}{n}\sum_{i=1}^{n}\operatorname{tr}\left(\sigma_{i}^{(t)}\tilde{\sigma}_{i}^{(t^{\prime})}\right)\right)\right),$
(6)

where $S_{T}(\rho)$ and $S_{T}(\tilde{\rho})$ are classical shadow
representations of $\rho$ and $\tilde{\rho}$, respectively. The computation
time for the inner product is $\mathcal{O}(nT^{2})$, linear in the system size
$n$ and quadratic in $T$, the number of copies of each quantum state which are
measured to construct the classical shadow.

###  IV.2 Rigorous guarantee

By statistical analysis, we can establish a rigorous guarantee for the
classical ML model
$\langle\alpha,\phi^{\mathrm{(shadow)}}(S_{T}(\rho))\rangle$, where $\alpha$
is the trainable vector defining the classifying hyperplane. The result is the
following theorem proven in Appendix J.

######  Theorem 2 (Classifying quantum phases of matter; informal).

If there is a nonlinear function of few-body reduced density matrices that
classifies phases, then the classical algorithm can learn to classify these
phases accurately. The required amount of training data and computation time
scale polynomially in system size.

If there is an efficient procedure based on _few-body reduced density
matrices_ for classifying phases, the proposed ML algorithm is guaranteed to
find the procedure efficiently. This includes local order parameters for
classifying symmetry breaking phases, and topological entanglement entropy in
a sufficiently large local region for partially classifying topological phases
kitaev2006topological ; levin2006detecting . We expect that, to classify
topological phases accurately, the classical ML will need access to local
regions that are sufficiently large compared to the correlation length, and as
we approach the phase boundary, the correlation length increases. As a result,
the classifying function for topological phases may depend on $r$-body
subsystems with a larger $r$, and the amount of training data and computation
time required would increase accordingly. Note that the classical ML not only
classifies phases accurately, but also constructs a classifying function
explicitly.

Our classical ML model may also be useful for classifying and understanding
symmetry-protected topological (SPT) phases. SPT phases are characterized much
like topological phases, but with the additional constraint that all
structures involved (states, Hamiltonians, and quantum circuits) respect a
particular symmetry. It is reasonable to expect that an SPT phase can be
identified by examining reduced density matrices on constant-size regions
Li2008 ; Pollmann2010 ; pollmann2012detection ; Haegeman2012 ; Shapourian2017
; Dehghani2021 , where the size of the region is large compared to the
correlation length. The existence of classifying functions based on reduced
matrices have been rigorously established in some cases Kitaev2006 ;
zhang2017quantum ; Hastings2015 ; Kapustin2020 ; Bachmann2020 ; Bachmann2014a
; tasakiPRL2018 ; Tasaki2020 . In Appendix K, we prove that the ML algorithm
is guaranteed to efficiently classify a class of gapped spin-$1$ chains in one
dimension. For more general SPT phases, the ML algorithm should be able to
corroborate known classification schemes, determine new and potentially more
compact classifiers, and shed light on interacting SPT phases in two or more
dimensions for which complete classification schemes have not yet been firmly
established.

The hypothesis of Theorem 2, stating that phases can be recognized by
inspecting regions of constant size independent of the total system size, is
particularly plausible for gapped phases, but might apply to some gapless
phases as well. Our classical ML model would be able to efficiently classify
such gapless phases. On the other hand, the contrapositive of Theorem 2
asserts that if the classical ML is not able to distinguish between two
distinct gapless phases, then nonlocal data is required to characterize at
least one of those phases.

##  V Numerical experiments

Figure 2: Numerical experiment for predicting ground-state properties in a 1D
Rydberg atom system with 51 atoms. (a) Hamiltonian. Illustration of the
Rydberg array geometry, Hamiltonian, and phases. (b) Phase diagram. The
system’s three distinct phases bernien2017probing are characterized by two
order parameters (for $Z_{2}$ and $Z_{3}$ orders). Training data are enclosed
by gray circles, and three specific testing points are indicated by the star,
diamond, and cross, respectively. (c) Local expectation values. We use
classical ML (the best model is selected from a set of ML models) to predict
the expectation values of Pauli operators $X_{i}$ and $Z_{i}$ for each atom at
the three testing points. We compare with “predictions” obtained from the
training data nearest to the testing points. The markers denote predicted
values, while the solid lines denote exact values obtained from DMRG.
Additional predictions are shown in Appendix D.1.

We have conducted numerical experiments assessing the performance of classical
ML algorithms in some practical settings. The results demonstrate that our
theoretical claims carry over to practice, with the results sometimes turning
out even better than our guarantees suggest.

###  V.1 Predicting ground state properties

For predicting ground states, we consider classical ML models encompassed by
Eq. (2). We examine various metrics $\kappa(x,x_{\ell})$ equivalent to
training neural networks with large hidden layers jacot2018neural ;
neuraltangents2020 or training kernel methods cortes1995support ;
murphy2012machine . We find the best ML model and the hyperparameters using a
validation set to minimize root-mean-square error (RMSE) and report the
predictions on a test set. The full details of the models and hyperparameters,
as well as their comparison, are given in Appendix D.2 and D.3.

_Rydberg atom chain_ — Our first example is a system of trapped Rydberg atoms
Fendley2004 ; browaeys_many-body_2020 , a programmable and highly-controlled
platform for Ising-type quantum simulations Schauss1455 ; Endres2016 ;
bernien2017probing ; Labuhn ; Misha256 ; 2020arXiv201212268S . Following
bernien2017probing , we consider a one-dimensional array of $n=51$ atoms, with
each atom effectively described as a two-level system composed of a ground
state $|g\rangle$ and a highly-excited Rydberg state $|r\rangle$. The atomic
chain is characterized by a Hamiltonian $H(x)$ (given in Figure 2(a)) whose
parameters are the laser detuning $x_{1}=\Delta/\Omega$ and the interaction
range $x_{2}=R_{b}/a$. The phase diagram (shown in Figure 2(b)) features a
disordered phase and several broken-symmetry phases, stemming from the
competition between the detuning and the Rydberg blockade (arising from the
repulsive Van der Waals interactions).

We trained a classical ML model using $20$ randomly chosen values of the
parameter $x=(x_{1},x_{2})$; these values are indicated by gray circles in
Figure 2(b). For each such $x$, an approximation to the exact ground state was
found using DMRG DMRG1 based on the formalism of matrix product states (MPS)
SCHOLLWOCK201196 . For each MPS, we performed $T=500$ randomized Pauli
measurements to construct a classical shadow. The classical ML then predicted
classical representations at the testing points in the parameter space, and
these predicted classical representations were used to estimate expectation
values of local observables at the testing points.

Predictions for expectation values of Pauli operators $Z_{i}$ and $X_{i}$ at
the testing points are shown in Figure 2(c), and found to agree well with
exact values obtained from the DMRG computation of the ground state at the
testing points. Additional predictions can be found in Appendix D.1. Also
shown are results from a more naive procedure, in which properties are
predicted using only the data at the point in the training set which is
closest to the testing point. The naive procedure predicts poorly,
illustrating that the considered classical ML model effectively leverages the
data from multiple points in the training set.

This example corroborates our expectation that classical machines can learn to
efficiently predict ground state representations. An important caveat is that
the rigorous guarantee in Theorem 1 applies only when the training points and
the testing points are sampled from the same phase, while in this example the
training data includes values of $x$ from three different phases.
Nevertheless, our numerics show that classical machines can still learn to
predict well.

Figure 3: Numerical experiment for predicting ground state properties in the
2D antiferromagnetic Heisenberg model. (a) Hamiltonian. Illustration of the
Heisenberg model geometry and Hamiltonian. We consider random couplings
$J_{ij}$, sampled uniformly from $[0,2]$. A particular instance is shown, with
coupling strength indicated by the thickness of the edges connecting lattice
points. (b) Two-point correlator. Exact values and ML predictions of the
expectation value of the correlation function
$C_{ij}=\tfrac{1}{3}(X_{i}X_{j}+Y_{i}Y_{j}+Z_{i}Z_{j})$ for all spin pairs
$(ij)$ in the lattice, for the Hamiltonian instance shown in (a). The absolute
value of $C_{ij}$ is represented by the size of each circle, while the
circle’s color indicates the actual value. (c) Prediction error. Each blue
point indicates the root-mean-square error (averaged over Heisenberg model
instances) of the correlation function for a particular pair $(ij)$, where the
estimate of $C_{ij}$ is obtained using a classical shadow with $T=500$
randomized Pauli measurements of the true ground state. Red points indicate
errors in ML predictions for $C_{ij}$.

_2D antiferromagnetic Heisenberg model_ — Our next example is the two-
dimensional antiferromagnetic Heisenberg model. Spin-$\tfrac{1}{2}$ particles
(i.e. qubits) occupy sites on a square lattice, and for each pair $(ij)$ of
neighboring sites the Hamiltonian contains a term
$J_{ij}\left(X_{i}X_{j}+Y_{i}Y_{j}+Z_{i}Z_{j}\right)$ where the couplings
$\\{J_{ij}\\}$ are uniformly sampled from the unit interval $[0,2]$. The
parameter $x$ is a list of all $J_{ij}$ couplings; hence in this case the
dimension of the parameter space is $m=O(n)$, where $n$ is the number of
qubits. The Hamiltonian $H(x)$ on a $5\times 5$ lattice is shown in Figure
3(a).

We trained a classical ML model using 90 randomly chosen values of the
parameter $x=\\{J_{ij}\\}$. For each such $x$, the exact ground state was
found using DMRG, and we simulated $T=500$ randomized Pauli measurements to
construct a classical shadow. The classical ML predicted the classical
representation at new values of $x$, and we used the predicted classical
representation to estimate a two-body correlation function, the expectation
value of $C_{ij}=\frac{1}{3}\left(X_{i}X_{j}+Y_{i}Y_{j}+Z_{i}Z_{j}\right)$,
for each pair of qubits $(ij)$. In Figure 3(b), the predicted and actual
values of the correlation function are displayed for a particular value of
$x$, showing reasonable agreement.

Figure 3(c) shows the prediction performance for all pairs of spins and for
variable system size. Each red point in the plot represents the RMSE in the
correlation function estimated using our predicted classical representation,
for a particular pair of spins and averaged over sampled values of $x$. For
comparison, each blue point is the RMSE when the correlation function is
predicted using the classical shadow obtained by measuring the actual ground
state $T=500$ times. For most correlation functions, the prediction error
achieved by the best classical ML model is comparable to the error achieved by
measuring the actual ground state.

###  V.2 Classifying quantum phases of matter

For classifying quantum phases of matter, we consider an unsupervised
classical ML model that constructs an  infinite-dimensional _nonlinear_
feature vector for each quantum state $\rho$ by applying the map
$\phi^{\mathrm{(shadow)}}$ in Eq. (5) with $\tau,\gamma=1$ to the classical
shadow $S_{T}(\rho)$ of the quantum state $\rho$. We then perform a principal
component analysis (PCA) pearson1901liii in the  infinite-dimensional
_nonlinear_ feature space.  The low-dimensional subspace found by PCA in the
nonlinear feature space corresponds to a nonlinear low-dimensional manifold in
the original quantum state space. This method is efficient using the shadow
kernel $k^{\mathrm{(shadow)}}$ given in Eq. (6) and the kernel PCA procedure
scholkopf1998nonlinear . Details are given in Appendix D.4 and D.5.

Bond-alternating XXZ model — We begin by considering the bond-alternating XXZ
model with $n=300$ spins. The Hamiltonian is given in Figure 4(a); it
encompasses the bond-alternating Heisenberg model ($\delta=1$) and the bosonic
version of the Su-Schrieffer-Heeger model ($\delta=0$) su1979solitons . The
phase diagram in Figure 4(b) is obtained by evaluating the partial reflection
many-body topological invariant pollmann2012detection ; elben2020many . There
are three different phases: trivial, symmetry-protected topological, and
symmetry broken.

For each value of $J$ and $\delta$ considered, we construct the exact ground
state using DMRG, and find its classical shadow by performing randomized Pauli
measurement $T=500$ times. We then consider a two-dimensional principal
subspace of the  infinite-dimensional nonlinear feature space found by the
unsupervised ML based on the shadow kernel, which is visualized in Figure 4(c,
d). We can clearly see that the different phases are well separated in the
principal subspace. This shows that even without any phase labels on the
training data, the ML model can already classify the phases accurately. Hence,
when trained with only a small amount of labeled data, the ML model will be
able to correctly classify the phases as guaranteed by Theorem 2.

Figure 4: Numerical experiments for classifying quantum phases in the bond-
alternating XXZ model. (a) Hamiltonian. Illustration of the model — a one-
dimensional qubit chain, where the coefficient of
$(X_{i}X_{i+1}+Y_{i}Y_{i+1}+\delta Z_{i}Z_{i+1})$ alternates between $J$ and
$J^{\prime}$. (b) Phase diagram. The system’s three distinct phases are
characterized by the many-body topological invariant $\tilde{Z}_{R}$ discussed
in Refs. pollmann2012detection ; elben2020many . Blue denotes
$\tilde{Z}_{R}=1$, red denotes $\tilde{Z}_{R}=-1$, and gray denotes
$\tilde{Z}_{R}\approx 0$. (c, d) Unsupervised phase classification. Bottom
panels: $\tilde{Z}_{R}$ vs. $J^{\prime}/J$ at cross sections (c) $\delta=0.5$
and (d) $\delta=3.0$ of the phase diagram. Top panels: visualization of the
quantum states projected to two dimensions using the unsupervised ML (
nonlinear PCA with shadow kernel). In all panels, colors of the points
indicate the value of $J^{\prime}/J$, indicating that the two phases naturally
cluster in the expressive feature space.

_Distinguishing a topological phase from a trivial phase_ — We consider the
task of distinguishing the toric code topological phase from the trivial phase
in a system of $n=200$ qubits. Figure 5(a) illustrates the sampled topological
and trivial states. We generate representatives of the nontrival topological
phase by applying low-depth geometrically local random quantum circuits to
Kitaev’s toric code state kitaev2003fault with code distance 10, and we
generate representatives of the trivial phase by applying random circuits to a
product state.

Randomized Pauli measurements are performed $T=500$ times to convert the
states to their classical shadows, and these classical shadows are mapped to
feature vectors in the high-dimensional feature space using the feature map
$\phi^{\mathrm{(shadow)}}$. Figure 5(b) displays a one-dimensional projection
of the feature space using the unsupervised classical ML for various values of
the circuit depth, indicating that the phases become harder to distinguish as
the circuit depth increases. In Figure 5(c), we show the classification
accuracy of the unsupervised classical ML model. We also compare to training
convolutional neural networks (CNN) that use measurement outcomes from the
Pauli-$6$ POVM carrasquilla2019reconstructing as input to learn an observable
for classifying the phases. Since Proposition 2 establishes that no observable
(even a global one) can classify topological phases, this CNN approach is
doomed to fail. On the other hand, if the CNN takes classical shadow
representations as input, then it can learn nonlinear functions and
successfully classify the phases.

##  VI Outlook

We have rigorously established that classical machine learning (ML)
algorithms, informed by data collected in physical experiments, can
effectively address  some quantum many-body problems. These results boost our
hopes that classical ML trained on experimental data can solve practical
problems in chemistry and materials science that would be too hard to solve
using classical processing alone.

Our arguments build on the concept of a classical shadow derived from
randomized Pauli measurements. We expect, though, that other succinct
classical representations of quantum states could be exploited by classical ML
with similarly powerful results. For example, some currently available quantum
simulators are highly programmable, but lack the local control needed to
perform arbitrary single-qubit Pauli measurements. Instead, after preparing a
many-body quantum state of interest, one might switch rapidly to a different
Hamiltonian and then allow the state to evolve for a short time before
performing a computational basis measurement. How can we make use of that
measurement data to predict properties reliably? Answering such questions, and
thereby expanding the reach of near-term programmable quantum platforms, will
be an important goal for future research cotler2021emergent ;
hu2021hamiltonian .

Viewed from a broader perspective, by illustrating how experimental data can
be exploited to make accurate predictions about features of quantum systems
that have never been studied directly, our work exemplifies a potentially
powerful methodology for advancing the physical sciences. With further
theoretical developments, perhaps we can learn how to use experimental data
that is already routinely available to accelerate the discovery of new
chemical compounds and materials with remarkable properties that could benefit
humanity.

Figure 5: Numerical experiments for distinguishing between trivial and
topological phases. (a) State generation. Trivial or topological states are
generated by applying local random quantum circuits of some circuit depth to a
product state or exactly-solved topological state, respectively. (b)
Unsupervised phase classification. visualization of the quantum states
projected to one dimension using the unsupervised ML ( nonlinear PCA with
shadow kernel), shown for varying circuit depth (divided by the “code
distance” 10, which quantifies the depth at which the topological properties
are washed out). The feature space is sufficiently expressive to resolve the
phases for a small enough depth without training, with classification becoming
more difficult as the depth increases. (c) Classification accuracy for three
ML algorithms described in Section V.2.

### Acknowledgments:

The authors thank Nir Bar-Gill, Juan Carrasquilla, Sitan Chen, Yifan Chen,
Andreas Elben, Matthew Fishman, Martin Fraas, Scott Glancy, Jeongwan Haah,
Felix Kueng, Jarrod McClean, Spiros Michalakis, Jacob Taylor, Yuan Su, and
Thomas Vidick for valuable input and inspiring discussions. HH thanks Andreas
Elben for providing the code on bond-alternating XXZ model. The numerical
simulations were performed on AWS EC2 computing infra-structure, using the
software packages ITensors itensor and PastaQ pastaq . HH is supported by the
J. Yang & Family Foundation. JP acknowledges funding from the U.S. Department
of Energy Office of Science, Office of Advanced Scientific Computing Research,
(DE-NA0003525, DE-SC0020290), and the National Science Foundation
(PHY-1733907). The Institute for Quantum Information and Matter is an NSF
Physics Frontiers Center. Contributions to this work by NIST, an agency of the
US government, are not subject to US copyright. Any mention of commercial
products does not indicate endorsement by NIST. VVA thanks Olga Albert, Halina
and Ryhor Kandratsenia, as well as Tatyana and Thomas Albert for providing
daycare support throughout this work.

Appendices: contents & navigation guide

A. Background on classical shadows A

B. Related work B

Estimating ground state properties $\bullet$ Classifying quantum phases of
matter $\bullet$ Classical representations of quantum systems

C. Neural networks with classical shadow for quantum many-body problems C

Predicting ground state representation $\bullet$ Classifying phases of matter

D. Details regarding numerical experiments D

Additional numerical experiments $\bullet$ Ground state properties of the
Rydberg atom chain $\bullet$ Ground state properties of the 2D
antiferromagnetic Heisenberg model $\bullet$ Classifying phases of the bond-
alternating XXZ model $\bullet$ Distinguishing a topological phase from a
trivial phase

E. Proof idea for the efficiency in predicting ground states E

Main result $\bullet$ Generalization to other systems and settings

F. Proof of efficiency for predicting ground states F

Overview for sample complexity upper bound $\bullet$ Controlling the
truncation error $\bullet$ Controlling generalization errors from using the
training data $\bullet$ Computational time for training and prediction
$\bullet$ Spectral gap implies smooth parametrizations

G. Sample complexity lower bound for predicting ground states G

Learning problem formulation $\bullet$ Communication protocol $\bullet$
Information-theoretic analysis

H. Computational hardness for non-ML algorithms to predict ground state
properties H

NP-hardness for estimating one-body observables in the ground state of 2D
Hamiltonians $\bullet$ Computational hardness for a class of Hamiltonians
based on factoring

I. No observable can classify topological phases I

J. Proof of efficiency for classifying phases of matter J

Training support vector machines $\bullet$ Prediction using support vector
machines $\bullet$ Kernel functions for classical shadows $\bullet$ Physical
assumptions about classifying quantum phases of matter $\bullet$ Training with
shadow kernels $\bullet$ Prediction based on shadow kernels

K. Classifying SPT phases with $O(2)$ symmetry using a few-body observable K

Symmetry-protected topological phases $\bullet$ $O(2)$-symmetric qutrit spin
chains

For those unfamiliar with the classical shadow formalism huang2020predicting ,
Appendix A provides a concise introduction and contains all the necessary
information on classical shadows to follow this paper. Discussion of related
literature on methods in many-body physics, in particular works for training
machine learning models to solve quantum many-body problems, is given in
Appendix B.

Readers interested in the numerical experiments can jump directly to Appendix
C, which provides a detailed discussion on practical approaches for combining
neural networks with classical shadows. The reader could then continue to
Appendix D, which gives additional numerical experiments, details of numerical
experiments, and example codes (in Python) for using the machine learning
models.

The rest of the appendices are dedicated to providing a rigorous understanding
in using classical machine learning models to solve quantum many-body
problems.

##### Predicting ground states

: We recommend that readers start with Appendix E, which provides the idea for
why classical machine learning models can be trained to predict ground state
representations of quantum systems. In Appendix E, we also provide the proof
idea for generalizing to other settings (such as Fermionic systems, long-range
interacting systems, etc). The detailed proof is given in Appendix F. In
Appendix G, we give a fundamental lower bound in the required data size for
learning to predict ground state properties for general classes of
Hamiltonians. In Appendix H, we show why non-ML algorithms cannot achieve a
similar guarantee as ML algorithms in predicting ground state representations.

##### Classifying phases of matter

: The reader could begin with the basic proposition given in Appendix I, which
shows that no (local or global) observable $\operatorname{tr}(O\rho)$ can be
used to classify topological phases. This motivates the need to consider
stronger machine learning models that can learn nonlinear functions in the
quantum state $\rho$. The readers could then proceed to Appendix J, which
provides a general theory for establishing provable guarantees in training
machine learning models based on classical shadows to classify quantum phases
of matter. In Appendix K, we briefly introduce symmetry-protected topological
phases and prove that the proposed machine learning model can classify a
particular subset of such phases.

##  APPENDIX A Background on classical shadows

The classical shadows formalism uses randomized (single-shot) measurements to
predict many properties of an unknown quantum state $\rho$ at once
huang2020predicting , see also paini2019approximate ; evans2019scalable . The
underlying idea dates back to Ohliger_2013 and also features prominently in
VanEnk2012 ; Elben2018 ; Vermersch2018 . In particular, the classical shadows
formalism comes with rigorous performance guarantees in terms of approximation
accuracy, classical storage, as well as data processing. Here, we focus on
randomized single-qubit Pauli measurements and repeat the following procedure
a total of $T$ times: (i) prepare an independent copy of $\rho$; (ii) select
$n$ single qubit Pauli measurements uniformly at random ($Z$, $X$ and $Y$
occur with probability $1/3$ each) and (iii), perform the associated
measurement to obtain $n$ classical bits ($+1$ if we measure ‘up’ and $-1$ if
we measure ‘down’). Subsequently, we store the associated post-measurement
state

$|s_{1}^{(t)}\rangle\otimes\cdots\otimes|s_{n}^{(t)}\rangle\quad\text{with}\quad|s_{1}^{(t)}\rangle,\ldots,|s_{n}^{(t)}\rangle\in\left\\{|0\rangle,|1\rangle,|+\rangle,|-\rangle,|\mathrm{i}+\rangle,|\mathrm{i}-\rangle\right\\}\subset\mathbb{C}^{2}$
(7)

in classical memory. This is very cheap, because there are only six
possibilities for each qubit. After $T$ repetitions, we obtain an entire
collection of $nT$ single-qubit states that we arrange in a two-dimensional
array:

$S_{T}(\rho)=\left\\{|s_{i}^{(t)}\rangle:\;i\in\\{1,\ldots,n\\},t\in\\{1,\ldots,T\\}\right\\}\in\left\\{|0\rangle,|1\rangle,|+\rangle,|-\rangle,|\mathrm{i}+\rangle,|\mathrm{i}-\rangle\right\\}^{n\times
T}$ (8)

The distribution of product states contains valuable information about the
underlying $n$-qubit density matrix $\rho$. In fact, we can use $S_{T}(\rho)$
to approximate $\rho$ via

$\sigma_{T}(\rho)=\frac{1}{T}\sum_{t=1}^{T}\left(3|s_{1}^{(t)}\rangle\\!\langle
s_{1}^{(t)}|-\mathbb{I}\right)\otimes\cdots\otimes\left(3|s_{n}^{(t)}\rangle\\!\langle
s_{n}^{(t)}|-\mathbb{I}\right),$ (9)

where $\mathbb{I}$ denotes the identity matrix (here, a 2-by-2 identity). It
is instructive to view this as the empirical average of $T$ independent and
identically (iid) random matrices. Each random matrix is an iid copy of
$\sigma_{1}(\rho)=\big{(}3|s_{1}\rangle\\!\langle
s_{1}|-\mathbb{I}\big{)}\otimes\cdots\otimes\big{(}3|s_{n}\rangle\\!\langle
s_{n}|-\mathbb{I}\big{)}$. Each tensor factor is guaranteed to have
eigenvalues $\lambda_{+}=2$ and $\lambda_{-}=-1$. This ensures that

$\displaystyle\mathrm{tr}\left(\sigma_{1}(\rho)\right)=$
$\displaystyle\mathrm{tr}\left(|s_{1}\rangle\\!\langle
s_{1}|-\mathbb{I}\right)\cdots\mathrm{tr}\left(|s_{n}\rangle\\!\langle
s_{n}|-\mathbb{I}\right)=1\quad\text{and}$ (10a)
$\displaystyle\left\|\sigma_{1}(\rho)\right\|_{p}=$
$\displaystyle\left\|3|s_{1}\rangle\\!\langle
s_{1}|-\mathbb{I}\right\|_{p}\cdots\left\|3|s_{n}\rangle\\!\langle
s_{n}|-\mathbb{I}\right\|_{p}=\left(|\lambda_{+}|^{p}+|\lambda_{-}|^{p}\right)^{n/p}=\left(2^{p}+1^{p}\right)^{n/p},$
(10b)

regardless of the concrete realization (and the underlying quantum state
$\rho$). The most relevant Schatten-$p$ norms are
$\|\sigma_{1}(\rho)\|_{1}=3^{n}$, $\|\sigma_{1}(\rho)\|_{2}=5^{n/2}$ and
$\|\sigma_{1}(\rho)\|_{\infty}=2^{n}$. Note, however, that the matrix
$\sigma_{1}(\rho)$ is never positive semidefinite.

The random matrix $\sigma_{1}(\rho)$ is a highly structured tensor product
that can assume a total of $6^{n}$ values. Each of them reflects the outcome
of performing randomly selected single-qubit Pauli measurements on the
$n$-qubit state $\rho$. Let us denote these Pauli matrices by
$W_{1},\ldots,W_{n}\in\left\\{X,Y,Z\right\\}$ and let
$o_{1},\ldots,o_{n}\in\left\\{\pm 1\right\\}$ be the observed outcomes ($+1$
if we measure ‘spin up’ and $-1$ if we measure ‘spin down’). Elementary
reformulations and Born’s rule then imply

$\displaystyle\sigma_{1}(\rho)=\frac{1}{2}\left(\mathbb{I}+3o_{1}W_{1}\right)\otimes\cdots\otimes\frac{1}{2}\left(\mathbb{I}+3o_{n}W_{n}\right)\quad\text{with
prob.}\quad\frac{1}{3^{n}}\mathrm{tr}\left(\frac{1}{2}(\mathbb{I}+o_{1}W_{1})\otimes\cdots\otimes\frac{1}{2}(\mathbb{I}+o_{n}W_{n})\rho\right).$
(11)

This construction ensures that $\sigma_{1}(\rho)$ exactly reproduces the
underlying quantum state $\rho$ in expectation. That is, if we average over
all $3^{n}$ choices of Pauli measurements and the associated (single-shot)
outcomes $o_{i}\in\left\\{\pm 1\right\\}$, we obtain

$\displaystyle\operatorname*{{\mathbb{E}}}_{s_{1},\ldots,s_{n}}\left[\sigma_{1}(\rho)\right]$
$\displaystyle=\operatorname*{{\mathbb{E}}}_{s_{1},\ldots,s_{n}}\left[\frac{1}{2}\left(\mathbb{I}+3o_{1}W_{1}\right)\otimes\cdots\otimes\frac{1}{2}\left(\mathbb{I}+3o_{n}3W_{n}\right)\right]$
(12a)
$\displaystyle=\sum_{W_{1},\ldots,W_{n}=X,Y,Z}\sum_{o_{1},\ldots,o_{n}=\pm
1}\frac{1}{3^{n}}\mathrm{tr}\left(\frac{1}{2}(\mathbb{I}+o_{1}W_{1})\otimes\cdots\otimes\frac{1}{2}(\mathbb{I}+o_{n}W_{n})\rho\right)$
(12b) $\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\times\frac{1}{2}\left(\mathbb{I}+o_{1}3W_{1}\right)\otimes\cdots\otimes\frac{1}{2}\left(\mathbb{I}+o_{n}3W_{n}\right)$
(12c) $\displaystyle=\rho.$ (12d)

We refer to Ref. huang2020predicting for a more detailed derivation and
context.

The classical shadow (9) attempts to approximate this expectation value by an
empirical average over $T$ independent samples, much like Monte Carlo sampling
approximates an integral. The accuracy of the approximation increases with
$T$, but insisting on accurate approximations of the global state $\rho$ is
prohibitively expensive. Known fundamental lower bounds flammia2012quantum ;
haah2017sample state that classical shadows of exponential size (at least)
$T=\Omega\left(2^{n}/\epsilon^{2}\right)$ are required to
$\epsilon$-approximate $\rho$ in trace distance. This quickly becomes
intractable in terms of both measurement budget, as well as classical storage
and processing.

This bleak picture lightens up considerably if we restrict our attention to
subsystem approximations. The classical shadow size required to accurately
approximate _all_ reduced $r$-body density matrices scales exponentially in
subsystem size $r$, but is independent of the total number of qubits $n$.

######  Lemma 1.

Fix $\epsilon,\delta\in(0,1)$, a subsystem size $r\leq n$ and let
$\sigma_{T}(\rho)$ be a classical shadow (9) of an $n$-qubit quantum state
$\rho$ with size

$T=(8/3)12^{r}\left(r\left(\log(n)+\log(12)\right)+\log(1/\delta)\right)/\epsilon^{2}=\mathcal{O}\left(r12^{r}\log(n/\delta)/\epsilon^{2}\right).$
(13)

Then, with probability at least $1-\delta$,

$\left\|\mathrm{tr}_{\neg A}\left(\sigma_{T}(\rho)\right)-\mathrm{tr}_{\neg
A}\left(\rho\right)\right\|_{1}\leq\epsilon\quad\text{for all subsystems
$A\subset\\{1,\ldots,n\\}$ with size $|A|\leq r$.}$ (14)

###### Proof.

Let us start by considering a fixed subsystem $A=\\{i_{1},\ldots,i_{r}\\}$
comprised of (at most) $r$ qubits. Use linearity to exchange partial trace
with expectation value to obtain

$\displaystyle\operatorname*{{\mathbb{E}}}_{s_{i_{1}}^{(t)},\ldots,s_{i_{r}}^{(t)}}\left[\big{(}3|s_{i_{1}}^{(t)}\rangle\\!\langle
s_{i_{1}}^{(t)}|-\mathbb{I}\big{)}\otimes\cdots\otimes\big{(}3|s_{i_{r}}^{(t)}\rangle\\!\langle
s_{i_{r}}^{(t)}|-\mathbb{I}\big{)}\right]$ (15a)
$\displaystyle=\mathrm{tr}_{\neg
A}\left(\operatorname*{{\mathbb{E}}}_{s_{1}^{(t)},\ldots,s_{n}^{(t)}}\left[\big{(}3|s_{1}^{(t)}\rangle\\!\langle
s_{1}^{(t)}|-\mathbb{I}\big{)}\otimes\cdots\otimes\big{(}3|s_{n}^{(t)}\rangle\\!\langle
s_{n}^{(t)}|-\mathbb{I}\big{)}\right]\right)$ (15b)
$\displaystyle=\mathrm{tr}_{\neg A}(\rho),$ (15c)

according to Eq. (12). In words, each reduced tensor product is an independent
random matrix that reproduces the $r$-qubit state $\mathrm{tr}_{\neg A}(\rho)$
exactly in expectation. Empirical averages of $T$ such independent and
identically distributed (iid) random matrices tend to concentrate sharply
around this expectation value. The matrix Bernstein inequality, see e.g.
Tro12:User-Friendly , provides powerful tail bounds in terms of operator norm
deviation. Let $X_{1},\ldots,X_{T}$ be iid random $D$-dimensional matrices
that obey $\|X_{t}-\operatorname*{{\mathbb{E}}}X_{t}\|_{\infty}\leq R$ almost
surely. Then, for $\tilde{\epsilon}>0$

$\mathrm{Pr}\left[\left\|\tfrac{1}{T}\sum_{t=1}^{T}\left(X_{t}-\operatorname*{{\mathbb{E}}}X_{t}\right)\right\|_{\infty}\geq\tilde{\epsilon}\right]\leq
2D\exp\left(-\frac{T\tilde{\epsilon}^{2}/2}{\sigma^{2}+R\tilde{\epsilon}/3}\right)\quad\text{where}\quad\sigma^{2}=\left\|\tfrac{1}{T}\sum_{t}\operatorname*{{\mathbb{E}}}X_{t}^{2}\right\|_{\infty}.$
(16)

Let us apply this tail bound to classical shadow concentration. We have $D\leq
2^{r}$ (at most $r$ qubits) and set
$X_{t}=\big{(}3|s_{i_{1}}^{(t)}\rangle\\!\langle
s_{i_{1}}^{(t)}-\mathbb{I}\big{)}\otimes\cdots\otimes\big{(}3|s_{i_{r}}^{(t)}\rangle\\!\langle
s_{i_{r}}^{(t)}|-\mathbb{I}\big{)}$, such that
$\operatorname*{{\mathbb{E}}}X_{t}=\mathrm{tr}_{\neg A}(\rho)$. Eq. (10) then
implies
$\|X_{t}-\operatorname*{{\mathbb{E}}}X_{t}\|_{\infty}\leq\|X_{t}\|+\|\operatorname*{{\mathbb{E}}}X_{t}\|_{\infty}\leq
2^{r}+1=:R$. Accurately bounding $\sigma^{2}$ is somewhat more involved, and
we turn to existing literature. A computation detailed in (guta2020fast, ,
Appendix C.3) yields $\sigma^{2}=3^{r}$. We are now ready to apply the matrix
Bernstein inequality. For $\tilde{\epsilon}>0$,

$\displaystyle\mathrm{Pr}\left[\left\|\mathrm{tr}_{\neg
A}\left(\sigma_{T}(\rho)\right)-\mathrm{tr}_{\neg
A}(\rho)\right\|_{\infty}\geq\tilde{\epsilon}\right]\leq
2^{r+1}\exp\left(-\frac{T\tilde{\epsilon}^{2}/2}{3^{r}+(2^{r}+1)\tilde{\epsilon}/3}\right)\leq
2^{r+1}\exp\left(-\frac{3T\tilde{\epsilon}^{2}}{8\times 3^{r}}\right),$ (17)

for $\tilde{\epsilon}\in(0,1)$. This is a powerful concentration statement in
operator norm. We can use the equivalence relation between trace- and operator
norm, $\|X\|_{\infty}\leq\|X\|_{1}\leq D\|X\|_{\infty}$, to obtain a tail
bound for trace norm deviations:

$\mathrm{Pr}\left[\left\|\mathrm{tr}_{\neg
A}\left(\sigma_{T}(\rho)\right)-\mathrm{tr}_{\neg
A}(\rho)\right\|_{\infty}\geq\epsilon\right]\leq\mathrm{Pr}\left[\left\|\mathrm{tr}_{\neg
A}\left(\sigma_{T}(\rho)\right)-\mathrm{tr}_{\neg
A}(\rho)\right\|_{1}\geq\epsilon/2^{r}\right]\leq
2^{r+1}\exp\left(-\frac{3T\epsilon^{2}}{8\times 12^{r}}\right).$ (18)

We see that, for a fixed subsystem $A=\\{i_{1},\ldots,i_{r}\\}$, the
probability of an $\epsilon$-deviation in trace distance is exponentially
suppressed in the size $T$ of the classical shadow. A union bound allows us to
extend this assertion to _all_ subsystems comprised of (at most) $r$ qubits:

$\displaystyle\mathrm{Pr}\left[\max_{A\subset\\{1,\ldots,n\\},|A|\leq
r}\left\|\mathrm{tr}_{\neg A}\left(\sigma_{T}(\rho)\right)-\mathrm{tr}_{\neg
A}(\rho)\right\|_{1}\geq\epsilon\right]$
$\displaystyle\leq\sum_{A\subset\\{1,\ldots,n\\},|A|\leq
r}\mathrm{Pr}\left[\left\|\mathrm{tr}_{\neg
A}\left(\sigma_{T}(\rho)\right)-\mathrm{tr}_{\neg
A}(\rho)\right\|_{1}\geq\epsilon\right]$ (19a) $\displaystyle\leq
n^{r}2^{r+1}\exp\left(-\frac{3T\epsilon^{2}}{8\times 12^{r}}\right).$ (19b)

Setting
$T=(8/3)12^{r}\left(\log(n^{r}12^{r})+\log(1/\delta)\right)/\epsilon^{2}=(8/3)12^{r}\left(r\left(\log(n)+\log(12)\right)+\log(1/\delta)\right)/\epsilon^{2}$
ensures that this upper bound on failure probability does not exceed $\delta$.
∎

##  APPENDIX B Related work

###  B.1 Estimating ground state properties

Determining the ground state of a system governed by a known many-body
Hamiltonian is a long-standing problem in quantum science. Despite having
several well-established and practically successful algorithms at our
disposal, we are typically faced with either a runtime that scales
exponentially with system size, and/or a lack of rigorous performance
guarantees. The literature is vast, and surveying it is beyond the scope of
this article. Instead, we review a few families of established algorithms in
order to put our work into proper context.

Density functional theory (DFT) has been a workhorse for determining
properties of interacting electronic systems in quantum chemistry and solid-
state physics. DFT recasts the problem of finding the many-body state with
minimal energy into finding a few-body energy functional. While the “true”
functional corresponding to the ground state is known to exist in theory
HohenbergKohn ; NobelKohn , determining it to polynomial accuracy in the
number of electrons is QMA-hard Schuch2009 . Various efficient approximations
to the true functional have seen much practical success, but they are
difficult to justify rigorously (except so far for some special cases
Lewin2020 ). These limitations present an opportunity for ML approaches to be
used instead of or to supplement DFT methods Schleder2019 ; moreno2020deep .
In moreno2020deep , they utilized the Hohenberg-Kohn theorem in DFT which
implies the existence of a mapping from single-particle densities to
correlation function. Then moreno2020deep proposed to train deep learning
models to approximately learn the mapping between correlation functions and
the densities. The authors in moreno2020deep found that across the phase
boundary, the mapping becomes hard to learn. This empirical observation could
be seen as a manifestation of the fact that Theorem 1 only holds when we do
not consider the presence of phase transitions. However, our numerical
experiments on Rydberg atoms suggest that ML models could still predict well
when there are multiple phase boundaries. Understanding when predictions would
be accurate in the presence of phase transitions is an interesting and
important question to study further.

The family of algorithms known as Quantum Monte Carlo (QMC) CEPERLEY555 ;
SandvikSSE utilizes probabilistic sampling techniques to estimate observable
properties at either finite or zero temperature. For ground states,
expectation values can be obtained using an imaginary-time evolution projector
or a high-power of the model Hamiltonian KaulQMC . The efficiency of QMC
methods depends on the structure of the Hamiltonian, specifically on whether
all of its off-diagonal matrix elements are negative (i.e. it is stoquastic
BravyiStoquastic ). In this case, the ground state wavefunction is real-valued
and positive, and the algorithmic complexity of the QMC estimators scales
polynomially with the number of particles. For non-stoquastic Hamiltonians,
the QMC suffers from the so-called sign problem, which makes evaluation of
statistical properties of the system NP-hard and renders QMC intractable for
large systems or low temperatures TroyerSignProblem . It is important to note
that, even for stoquastic Hamiltonians with polynomial computational
complexity, the success of QMC simulations heavily relies on the existence of
efficient sampling schemes (e.g. cluster updates) which are sufficiently
ergodic, and leading to small auto-correlation time EvertzLoop . In general,
it is not possible to prove the existence of such update algorithms, nor their
ergodicity.

An alternative approach to solve for the ground state properties of a many-
body Hamiltonian is based on the variational principle in quantum mechanics,
which states that the expectation value of the energy on any valid
wavefunction is always greater or equal than the ground state energy. It is
then possible to design classical parametric representations of the many-body
wavefunction, and update their parameters to minimize the corresponding energy
estimator. A notable example is the density-matrix renormalization group DMRG1
; DMRG2 (DMRG). This algorithm can be interpreted as a variational
optimization of a Matrix Product State (MPS) Garcia07 ; Vestraete2008 ;
SCHOLLWOCK201196 , which is a local decomposition of a wavefunction as a one-
dimensional tensor network. These parametrized wavefunctions display area law
of entanglement and exponentially-decaying correlations ORUS2014117 , which
lend themselves most effective for systems described by one-dimensional gapped
Hamiltonians. Furthermore, a standout feature of DMRG is that modifications of
the original procedure, such as rigorous renormalization group algorithms
arad2017rigorous , are guaranteed to find the ground state of one-dimensional
geometrically-local gapped Hamiltonians in polynomial time
landau2015polynomial ; arad2017rigorous ; abrahamsen2020polynomialtime . In
two spatial dimensions, MPS-based DMRG can still be applied to solve for
ground states Stoudenmire12 ; Wu2019 ; Moore2020 , though it suffers an
exponential scaling in one of the two linear dimensions of the system.
Projected entangled pair states (PEPS) Jordan2008 ; Vanderstraeten2016 ;
Corboz2016 , the two-dimensional generalization of MPSs, are instead a more
suitable ansatz in this context. However, while improved algorithms for PEPS
optimization are routinely put forward Zalatel2020 ; Haghshenas2019 ;
Hyatt2019 , the same level of performance achieved by DMRG in 1d systems is
still out of reach.

Another class of variational wavefunctions that has recently received a lot of
attention are neural-network quantum states Carleo_2017 . In this framework, a
neural network is used as a parametric function approximator of a many-body
wavefunction $\psi_{\lambda}(\sigma)=\langle\sigma|\psi\rangle$, where the
classical state $\sigma$ is interpreted as the neural-network input, and
$\lambda$ is a set of neural-network parameters (i.e. weights and biases). In
a variational setting, these parameters are iteratively optimized to lower the
total energy becca_sorella_2017 , or additionally the energy variance. Neural-
network quantum states have been explored in a variety of setups, including
topological phases dassarma2017 , Fermi-Hubbard models Nomura2017 , molecular
ground states choo_fermionicnqs2020 , frustrated magnetism Ferrari2019 , and
more Glasser2018 ; Choo2018 ; rbm_su2 ; morawetz2020 ; Luo2021 . The auto-
regressive property of some types of neural networks (e.g. recurrent neural
networks, transformers, etc.) has also been leveraged to improve convergence
of variational Monte Carlo hibat2020recurrent . In contrast to tensor-network
states, this class of wavefunctions can more easily display non-local
correlations, allowing in principle to capture quantum states with higher
entanglement Deng2017 .

Another class of machine learning methods gilmer2017neural ;
vargas2018extrapolating ; schutt2019unifying ; qiao2020orbnet ;
kawai2020predicting train neural networks to predict the ground state or
excited state properties directly. The input to the neural network is a
description of the Hamiltonian, and the output is a ground state property of
interest. The training data is a set of different Hamiltonians (inputs) and
their corresponding ground state properties (outputs). This class of ML
methods is closest to the setting considered in this paper. Methods in this
class lack rigorous guarantees, so it is not clear when such approaches could
outperform non-ML algorithms. In one of our main contributions, given in
Theorem 1 and Proposition 1, we introduce an ML model that, when trained with
experimental data, can accurately predict ground state representations better
than any classical algorithm that does not learn from data. Our model is
relatively basic, utilizing the well-known $l_{2}$-Dirichlet kernel, but it is
already enough to establish a rigorous guarantee. Similarly determining when
other ML models yield an advantage over non-ML algorithms is an interesting
topic for future work.

###  B.2 Classifying quantum phases of matter

Proposals for classifying quantum phases of matter abound. These include
quantum neural networks cong2019quantum , classical neural networks
carrasquilla2017nature ; evert2017nature ; neupeurt2017 ; Beach2017 ;
Greplova_2020 , or other classical machine learning models leiwang2016 ;
Wetzel2017 ; rodriguez2019identifying ; che2020topological. Since these models
do not come with rigorous guarantees, relying on them too much can lead one
astray. For example, some deep neural networks can misclassify the original
phase if the corresponding state is distorted by noise, even if the distortion
is very slight jiang2019vulnerability .

In this work, we provide rigorous machine learning approaches that are
guaranteed to classify accurately under the specified conditions given in
Theorem 2. We believe similar analyses can be performed on other machine
learning models to understand their limitation and potential, which will be an
important future direction. The ML model used in rodriguez2019identifying ,
which is based on defining diffusion maps over classical spins systems, is the
ML approach that seems most similar to that used in our work. Hence, it is
very likely that the models considered in rodriguez2019identifying can be
rigorously analyzed via similar techniques. Neural network approaches for
classifying phases of matter will be harder to analyze, but one should be able
to study neural network with large hidden layers using the theory developed in
this work and the theory of neural tangent kernels jacot2018neural ;
neuraltangents2020 .

###  B.3 Classical representations of quantum systems

One of the most important ingredients in designing classical ML procedures for
understanding quantum spin systems is the construction of efficient classical
representations of the underlying quantum systems. The properties of the
quantum system retained by the classical representation directly determine the
set of functions the classical ML procedure can learn. The classical shadow
formalism huang2020predicting ; paini2019approximate , developed by some of us
and others and used throughout this work, is a versatile framework for this
purpose. It has been extended to fermionic systems zhao2021fermionic ;
hadfield2020measurements , suggesting that our ML approaches may be extendable
as well. Classical shadows have also been shown to allow sample-efficient
reconstruction of Hamiltonian from thermal states anurag2021sample , although
such an algorithm is not yet time efficient. However, classical shadows
provide only one of many promising and actively studied approaches for
efficient representation cotler2020quantum ; paini2019approximate ;
huang2020predicting ; evans2019scalable ; koh2020classical ; chen2020robust ;
huang2021efficient .

While the curse of dimensionality prevents one from representing general
quantum spin systems both exactly and efficiently, simplifying assumptions can
lift the curse and drastically reduce both the overhead and complexity of
representing and characterizing the system. A prime example is a classical
system, whose Hamiltonian is diagonal in the computational basis. ML methods
for such systems, such as those in Ref. rodriguez2019identifying (discussed
above), do not require an additional quantum-to-classical compression. Another
example, relevant to electronic material characterization, is the family of
solid-state band insulators topobook — gapped two- and three-dimensional non-
interacing fermionic systems with various crystalline symmetries. Their myriad
topological phases can typically be characterized by data at a discrete set of
high-symmetry points in the Brillouin zone Bacry ; Bradlyn2017 ; Po2017 ;
Kruthoff2017 ; Po2020 ; Cano2021 . The techniques developed here should pave
the way for certifying accuracy of current ML methods for band insulator
characterization (e.g., Sun2018 ; Ming2019 ; Balabanov2021 ; Scheurer2020 ;
Claussen2020 ; Andrejevic2020 ; Peano2021 ) as well as developing new ones.

##  APPENDIX C Neural networks with classical shadow for quantum many-body
problems

Imposing inductive biases in the ML model is a common technique for boosting
the prediction performance of ML models. One approach is to enhance the
proposed ML algorithms with neural networks, such as convolutional or graph
neural networks. These neural networks could better capture structure of the
underlying function we are trying to learn and hence may require significantly
less data than the very expressive ML model given in the main text. We leave
the proof that neural network enhancements can lead to better prediction
performance as a goal for future work.

There are multiple ways of combining classical shadows and neural networks.
Here, we will only showcase one such approach by utilizing the theory of
neural tangent kernels jacot2018neural . Remarkably, this theory allows us to
efficiently train various types of neural networks (convolutional/graph/etc.)
with an infinite number of neurons in each hidden layer (infinite width). As
such, this line of work has gained a lot of attention du2019graph ;
arora2019exact ; neuraltangents2020 in recent years. In the limit of infinite
width, one can analytically solve for the neural network after training on a
set of data $\\{x_{\ell},y_{\ell}\\}_{\ell=1}^{N}$, where $x_{\ell}$ and
$y_{\ell}$ are vectors of some size. For example, consider training a neural
network that takes in a vector $x$ and produces a vector
$f^{\mathrm{NN}}_{\theta}(x)$ through the following optimization problem using
gradient descent,

$\min_{\theta}\sum_{\ell=1}^{N}\left\lVert
f^{\mathrm{NN}}_{\theta}(x_{\ell})-{y}_{\ell}\right\rVert_{2}^{2},$ (20)

where we begin on a randomly initialized $\theta$. Note that due to the
infinite number of neurons, $\theta$ is a vector of infinite dimension. The
trained neural network $f^{\mathrm{NN}}_{\theta^{*}}({x})$ can always be
written in the following form

$f^{\mathrm{NN}}_{\theta^{*}}(x)=\sum_{\ell=1}^{N}\sum_{\ell^{\prime}=1}^{N}k^{\mathrm{(NTK)}}(x,x_{\ell})(K^{-1})_{\ell\ell^{\prime}}y_{\ell^{\prime}},$
(21)

where $k^{\mathrm{(NTK)}}(x,x^{\prime})$ is a function called the neural
tangent kernel jacot2018neural , and
$K_{\ell,\ell^{\prime}}=k^{\mathrm{(NTK)}}({x}_{\ell},{x}_{\ell^{\prime}})$ is
the kernel matrix of the neural tangent kernel. One can see that the infinite-
dimensional vector $\theta^{*}$ does not appear on the right hand side of Eq.
(21). And as long as we can efficiently evaluate the neural tangent kernel
$k^{\mathrm{(NTK)}}(x,x^{\prime})$, we can evaluate the infinite-dimensional
neural network in polynomial time. This is the main contribution of
jacot2018neural , which enables one to efficiently train infinite-width neural
networks. For a given neural network architecture, one can compute
$k^{\mathrm{(NTK)}}({x},{x}^{\prime})$ efficiently using open-source software,
such as neuraltangents2020 . In Appendix D.2, we give the code for training
infinite-width neural networks using the open-source software: Neural Tangents
neuraltangents2020 .

###  C.1 Predicting ground state representation

For the task of predicting ground state representation, we consider the
training data to be

$\big{\\{}{x}_{\ell}\rightarrow\sigma_{T}(\rho({x}_{\ell}))\big{\\}}_{\ell=1}^{N}\leavevmode\nobreak\
,$ (22)

where $\sigma_{T}(\rho({x}_{\ell}))$ is the classical shadow representation of
$\rho({x}_{\ell})$ given in Eqs. (1) and (9) based on $T$ randomized Pauli
measurements. Recall that $\sigma_{T}(\rho({x}_{\ell}))$ is a $2^{n}\times
2^{n}$ matrix that reproduces $\rho({x}_{\ell})$ in expectation over the
randomized Pauli measurements. Suppose we now train an infinite-width neural
network parameterized by $\theta$ that takes in an input ${x}$ and produces an
exponential-size matrix $\sigma^{\mathrm{NN}}_{\theta}({x})$, by solving the
optimization problem

$\displaystyle\min_{\theta}$
$\displaystyle\quad\sum_{\ell=1}^{N}\left\lVert\sigma^{\mathrm{NN}}_{\theta}({x}_{\ell})-\sigma_{T}(\rho({x}_{\ell}))\right\rVert_{F}^{2}.$
(23)

The squared Frobenius difference between two matrices is equal to the squared
Euclidean norm of their vectorizations (flattenings). In turn, the theory of
infinite-width neural networks jacot2018neural shows that the trained neural
network $\sigma^{\mathrm{NN}}_{\theta^{*}}({x})$ could be written in the form

$\sigma^{\mathrm{NN}}_{\theta^{*}}({x})=\sum_{\ell=1}^{N}\sum_{\ell^{\prime}=1}^{N}k^{\mathrm{(NTK)}}({x},{x}_{\ell})(K^{-1})_{\ell\ell^{\prime}}\sigma_{T}(\rho({x}_{\ell^{\prime}})).$
(24)

The kernel function $k^{\mathrm{(NTK)}}({x},{x}^{\prime})$ depends on the
neural network architecture and could be calculated utilizing existing open-
source software neuraltangents2020 . This also falls into the general form
shown in the main text; see Eq. (2). Hence, training an infinite-width neural
networks to predict an exponentially large density matrix can be done
efficiently on a classical computer. For a given neural network architecture,
all one has to do is compute the kernel function
$k^{\mathrm{(NTK)}}({x},{x}^{\prime})$. Then the neural network optimized
using the training data could be analytically solved as given in Eq. (24). To
estimate a property on the predicted ground state using the neural network is
as simple as evaluating

$\operatorname{tr}(O\sigma^{\mathrm{NN}}_{\theta^{*}}({x}))=\sum_{\ell=1}^{N}\sum_{\ell^{\prime}=1}^{N}k^{\mathrm{(NTK)}}({x},{x}_{\ell})(K^{-1})_{\ell\ell^{\prime}}\operatorname{tr}(O\sigma_{T}(\rho({x}_{\ell^{\prime}}))),$
(25)

which can be done by first computing
$\operatorname{tr}(O\sigma_{T}(\rho({x}_{\ell}))),\forall\ell=1,\ldots,N$ and
compute the linear interpolation.

###  C.2 Classifying phases of matter

We want to learn how to classify two phases of $n$-qubit states. A fully
classical training set would simply consist of $N$ labeled classical
representations of quantum states $\\{\rho_{\ell}\rightarrow
y_{\ell}\\}_{\ell=1}^{N}$, where $y_{\ell}=+1$ ($-1$) if $\rho_{\ell}$ belongs
to phase $A$ ($B$). However, insisting on perfect knowledge of each
$\rho_{\ell}$ is impractical for a variety of reasons. Instead, we assume that
we have access to classical shadows of $\rho_{\ell}$. The raw data
$S_{T}(\rho_{\ell})$ behind each classical shadow is a 2-dimensional array,

$\displaystyle
S_{T}(\rho_{\ell})=\left\\{|s_{i}^{(t)}\rangle:\;i\in\\{1,\ldots,n\\},t\in\\{1,\ldots,T\\}\right\\}\quad\text{where}\quad|s_{i}^{(t)}\rangle\in\\{\ket{0},\ket{1},\ket{+},\ket{-},\ket{\mathrm{i}+},\ket{\mathrm{i}-}\\}.$
(26)

In the main text, we propose to use this data to train a support vector
machine based on the shadow kernel

$k^{\mathrm{(shadow)}}\left(S_{T}(\rho_{\ell}),\tilde{S}_{T}(\rho_{\ell^{\prime}})\right)=\exp\left(\frac{\tau}{T^{2}}\sum_{t,t^{\prime}=1}^{T}\exp\left(\frac{\gamma}{n}\sum_{i=1}^{n}\mathrm{tr}\left(\big{(}3|s_{i}^{(t)}\rangle\\!\langle
s_{i}^{(t)}|-\mathbb{I}\big{)}\big{(}3|\tilde{s}_{i}^{(t)}\rangle\\!\langle\tilde{s}_{i}^{(t)}|-\mathbb{I}\big{)}\right)\right)\right).$
(27)

This specific choice of (deterministic) kernel function allows us to carry out
a thorough theoretical analysis of the entire learning procedure; see Appendix
J.

But there are other sensible kernels that may perform even better in practice.
For instance, we could feed the two-dimensional data array (26) into a neural
network architecture, e.g. a convolutional neural network. In the limit of an
infinite number of neurons in each hidden layer, this produces the neural
tangent kernel
$k^{\mathrm{(NTK)}}\left(S_{T}(\rho_{\ell}),\tilde{S}_{T}(\rho_{\ell^{\prime}})\right)$
jacot2018neural . This kernel is positive-semidefinite and should be viewed as
a measure of similarity induced by the trained neural network. Mercer’s
theorem mercer1909functions allows us to make this intuition precise by
reformulating the neural tangent kernel as a Gram matrix in feature space:

$k^{\mathrm{(NTK)}}\left(S_{T}(\rho_{\ell}),\tilde{S}_{T}(\rho_{\ell^{\prime}})\right)=\left\langle\phi^{\mathrm{(NTK)}}\left(S_{T}(\rho_{\ell})\right),\phi^{\mathrm{(NTK)}}\left(\tilde{S}_{T}(\rho_{\ell^{\prime}})\right)\right\rangle\leavevmode\nobreak\
.$ (28)

Hence, any infinite-width neural network with input array $S_{T}(\rho)$
induces a feature map $\phi^{\mathrm{(NTK)}}$ that can be used instead of the
doubly-infinite feature map $\phi^{\mathrm{(shadow)}}$ (5) that is associated
with the shadow kernel (27).

##  APPENDIX D Details regarding numerical experiments

In this appendix, we provide additional numerical experiments as well as more
details about the numerical experiments described in the main text.

###  D.1 Additional numerical experiments

Figure 6: Numerical experiment for predicting ground state properties
(Pauli-$Z$ in each atom) in a 1D Rydberg atom system with 51 atoms. We use the
classical ML to predict the ground state properties at the three testing
points. Also shown are “predictions” obtained from the training data nearest
to the testing points. The markers denote predicted values, while the solid
lines denote exact values obtained from DMRG.

Figure 7: Numerical experiment for predicting ground state properties
(Pauli-$X$ in each atom) in a 1D Rydberg atom system with 51 atoms. We use the
classical ML to predict the ground state properties at the three testing
points. Also shown are “predictions” obtained from the training data nearest
to the testing points. The markers denote predicted values, while the solid
lines denote exact values obtained from DMRG.

Figure 8: “Predictions” obtained by performing bivariate B-spline
interpolation using the training data. The markers denote interpolated values,
while the solid lines denote exact values obtained from DMRG.  Figure 9:
Numerical experiments for distinguishing trivial and topological phases.
Trivial or topological states are generated by applying low-depth local random
quantum circuits to a product state or exactly solved topological state
respectively. (a) Kernel matrix for SPT/trivial phases The exactly solved
topological state is the cluster state. The $(i,j)$-entry denotes the inner
product of the $i$-th and $j$-th feature vectors in the infinite-dimensional
feature space defined by the classical shadow representation. To the left,
states from the two phases are randomly mixed. To the right, the two phases
are ordered. (b) Kernel matrix for topologically-ordered/trivial phases. The
exactly solved topological state is the toric code ground state.

_Rydberg atom chain_ — In the main text, we have provided partial prediction
outcomes for a one-dimensional chain of $n=51$ Rydberg atoms; see Figure 2.
Here, we supply predictions of expectation values of Pauli operators $Z_{i}$
and $X_{i}$ on all $51$ atoms at the testing points marked in Figure 2(b).
These are shown in Figure 8 and Figure 8, respectively. These extend the more
restricted presentation in the main text to all qubits. In Figure 8, we show a
different baseline considering bivariate B-spline interpolation from the
training data.

_Distinguishing an SPT phase from a trivial phase_ — We consider a one-
dimensional chain of $n=50$ qubits with $Z_{2}\times Z_{2}$ symmetry. The 1D
cluster state is in the nontrivial SPT phase. We generate other
representatives of the nontrivial SPT phase by applying symmetric depth-$3$
geometrically local random quantum circuits to the cluster state, and we
generate representatives of the trivial phase by applying symmetric depth-$3$
random circuits to a product state.

Randomized Pauli measurements are performed $T=500$ times to convert the
states to their classical shadows, and these classical shadows are mapped to
feature vectors in the infinite-dimensional feature space using the feature
map $\phi^{\mathrm{(shadow)}}$ (5). In Figure 9(a), inner products of feature
vectors (matrix elements of the shadow kernel) are displayed. Figure 9(b)
shows the feature vectors projected onto a two-dimensional subspace using
_nonlinear_ principal component analysis (PCA) based on the shadow kernel
$k^{\mathrm{(shadow)}}$. Both figures show that feature vectors representing
distinct phases can be distinguished easily. Correspondingly, the classical ML
efficiently learns how to classify phases accurately, even if the training
data is unlabeled.

_Distinguishing a topologically-ordered phase from a trivial phase_ — We
consider the task of distinguishing the toric code kitaev2003fault
topologically-ordered phase from the trivial phase in a system of $n=200$
qubits. We generate other representatives of the topologically-ordered phase
by applying two-dimensional depth-$3$ geometrically local random quantum
circuits to the toric code state, and we generate representatives of the
trivial phase by applying two-dimensional depth-$3$ random circuits to a
product state.

Randomized Pauli measurements are performed $T=500$ times to convert the
states to their classical shadows, and these classical shadows are mapped to
feature vectors in the infinite-dimensional feature space using the feature
map $\phi^{\mathrm{(shadow)}}$. In Figure 9(c, d), inner products of feature
vectors (matrix elements of the shadow kernel) and the projection of feature
space data onto the two-dimensional subspace spanned by the largest principal
components is shown. Once more, one can clearly see that feature vectors
representing distinct phases can be distinguished easily. Correspondingly, the
classical ML efficiently learns how to classify phases accurately, even if the
training data is unlabeled.

###  D.2 Ground state properties of the Rydberg atom chain

Our first example is a one-dimensional chain of $n=51$ Rydberg atoms
Fendley2004 ; browaeys_many-body_2020 ; bernien2017probing . Each atom can be
in either its ground state or a highly excited Rydberg state. Such systems can
effectively be regarded as a qubit, where the basis state $|0\rangle$ is the
ground state $|g\rangle$ and the basis state $|1\rangle$ is the Rydberg state
$|r\rangle$. The Hamiltonian of the atomic chain is

$H=\frac{\Omega}{2}\sum_{i}X_{i}-\Delta\sum_{i}N_{i}+\Omega\sum_{i<j}\left(\frac{R_{b}}{a|i-j|}\right)^{6}N_{i}N_{j}\leavevmode\nobreak\
,$ (29)

where $\Omega$ is the (fixed) Rabi frequency, $\Delta$ is the laser detuning,
$N_{i}$ is the Rydberg occupation number operator, $a$ is the separations of
the atoms, and $R_{b}$ is the so called Rydberg blockade radius. For large and
negative $\Delta$, the ground state of $H$ is a vacuum state, where all atoms
are in the ground state $|g\rangle$. In contrast, for large and positive
$\Delta$, different broken-symmetry ground states can be engineered depending
on the value of $R_{b}$.

Approximations of the exact ground states of the Rydberg chain were found
using the density-matrix renormalization group (DMRG) based on matrix product
states (MPS). Starting from a random MPS with bond dimension $\chi=10$, we
variationally optimize the MPS using a singular value decomposition (SVD)
cutoff of $10^{-9}$. We perform a number of DMRG sweeps until the change in
energy is below $\epsilon=10^{-6}$. Upon convergence, we perform randomized
Pauli measurements simply by performing local rotations into the corresponding
Pauli bases, and sampling the resulting state FerrisSampling .

In Figure 2(b), the color in the phase diagram corresponds to the phase
obtained by two order parameters for characterizing $Z_{2}$ and $Z_{3}$ order.
For $Z_{2}$ order, where the atoms are in $\ket{rgrgrg\ldots}$ or
$\ket{grgrgr\ldots}$, we consider the order parameter,

$O_{Z_{2}}=\frac{1}{n-1}\sum_{i=1}^{n-1}\left(\lvert
r_{i}g_{i+1}\rangle\\!\langle r_{i}g_{i+1}\rvert+\lvert
g_{i}r_{i+1}\rangle\\!\langle g_{i}r_{i+1}\rvert\right).$ (30)

For $Z_{3}$ order, where the atoms are in $\ket{rggrgg\ldots}$ or
$\ket{grggrg\ldots}$ or $\ket{ggrggr\ldots}$, we consider the order parameter,

$O_{Z_{3}}=\frac{1}{n-2}\sum_{i=1}^{n-2}\left(\lvert
r_{i}g_{i+1}g_{i+2}\rangle\\!\langle r_{i}g_{i+1}g_{i+2}\rvert+\lvert
g_{i}r_{i+1}g_{i+2}\rangle\\!\langle g_{i}r_{i+1}g_{i+2}\rvert+\lvert
g_{i}g_{i+1}r_{i+2}\rangle\\!\langle g_{i}g_{i+1}r_{i+2}\rvert\right).$ (31)

We estimate the two order parameters of the ground state $\rho$. First we
check which order parameter ($O_{Z_{2}}$ or $O_{Z_{3}}$) yields a larger
expectation value. Then, we check if that expectation value is larger than the
threshold value $0.8$. If $O_{Z_{2}}>O_{Z_{3}}$ and $O_{Z_{2}}>0.8$, we
associate the state with the $Z_{2}$-order phase (red color). Else if
$O_{Z_{3}}>O_{Z_{2}}$ and $O_{Z_{3}}>0.8$, we say that the state is in the
$Z_{3}$-order phase (vanilla color). If neither of these conditions is
satisfied (both expectation values are less than $0.8$), we assign the
disordered phase (blue color) to this state.

For the Rydberg atom experiment, the input parameter vector $x$ is two-
dimensional. We first normalize the values to lie within a square
$[-1,1]^{2}$. Then we consider classical machine learning models given by

$\hat{\sigma}_{N}(x)=\sum_{\ell=1}^{N}\kappa(x,x_{\ell})\sigma_{T}(x_{\ell})=\sum_{\ell=1}^{N}\underbrace{\left(\sum_{\ell^{\prime}=1}^{N}k(x,x_{\ell^{\prime}})(K+\lambda
I)^{-1}_{\ell^{\prime}\ell}\right)}_{\kappa(x,x_{\ell})}\sigma_{T}(x_{\ell}),$
(32)

where $\lambda>0$ is a parameter to regularize the model when $K$ is not
invertible, $\sigma_{T}(x_{\ell})$ is shorthand for
$\sigma_{T}\left(\rho_{\ell}\right)$ and denotes the classical shadow
representation of the ground state $\rho_{\ell}=\rho(x_{\ell})$ under $T$
randomized Pauli measurements. Moreover, $K_{ij}=k(x_{i},x_{j})$ is the kernel
matrix, $k(x,x^{\prime})$ is a kernel function, and $\kappa(x,x_{\ell})$ is a
function that depends on the kernel function, the kernel matrix $K$, and
$\lambda$. We consider a set of different regularization parameters,

$\lambda\in\\{0.0125,0.025,0.05,0.125,0.25,0.5,1.0,2.0,4.0,8.0\\},$ (33)

and we also consider a set of different kernel functions
$k(x,x^{\prime})=\tilde{k}(x,x^{\prime})/\sqrt{\tilde{k}(x,x)\tilde{k}(x^{\prime},x^{\prime})}$,
where

$\displaystyle\tilde{k}(x,x^{\prime})$ $\displaystyle=\exp(-\gamma\left\lVert
x-x^{\prime}\right\rVert_{2}^{2}),$ $\displaystyle\quad\mbox{(Gaussian
kernel)},$ (34a) $\displaystyle\tilde{k}(x,x^{\prime})$
$\displaystyle=\sum_{k_{1}=-3}^{3}\sum_{k_{2}=-3}^{3}\cos\left(\pi(k_{1}(x_{1}-x_{1}^{\prime})+k_{2}(x_{2}-x_{2}^{\prime}))\right),$
$\displaystyle\quad\mbox{(Dirichlet kernel)},$ (34b)
$\displaystyle\tilde{k}(x,x^{\prime})$
$\displaystyle=k^{\mathrm{(NTK)}}(x,x^{\prime}),$
$\displaystyle\quad\mbox{(Neural tangent kernel)}.$ (34c)

The hyperparameter $\gamma>0$ in the Gaussian kernel is chosen to be equal to
$N^{2}/\sum_{i=1}^{N}\sum_{j=1}^{N}\left\lVert
x_{i}-x_{j}\right\rVert_{2}^{2}$, the inverse of the average distance between
$x_{i}$ and $x_{j}$. We consider the neural tangent kernel
$k^{\mathrm{(NTK)}}(x,x^{\prime})$ jacot2018neural ; neuraltangents2020 that
is equivalent to an infinite-width feed-forward neural network with $2,3,4,5$
hidden layers and that uses the rectified linear unit (ReLU) as the activation
function. Computing the neural tangent kernel can be implemented easily using
the open-source software Neural Tangents neuraltangents2020 . Suppose that the
input data $\\{x_{\ell}\\}_{\ell=1}^{N}$ is stored in a numpy array of size
$N\times m$, denoted as dataX in the following code. We can use then use
following code to generate the neural tangent kernel matrix. The imported
package neural_tangents can be downloaded from
https://github.com/google/neural-tangents.

⬇

import numpy as np

from neural_tangents import stax

init_fn, apply_fn, kernel_fn = stax.serial(

stax.Dense(32), stax.Relu(),

stax.Dense(32), stax.Relu(),

stax.Dense(1)

)

kernel_NN2 = kernel_fn(dataX, dataX, ’ntk’)

init_fn, apply_fn, kernel_fn = stax.serial(

stax.Dense(32), stax.Relu(),

stax.Dense(32), stax.Relu(),

stax.Dense(32), stax.Relu(),

stax.Dense(1)

)

kernel_NN3 = kernel_fn(dataX, dataX, ’ntk’)

init_fn, apply_fn, kernel_fn = stax.serial(

stax.Dense(32), stax.Relu(),

stax.Dense(32), stax.Relu(),

stax.Dense(32), stax.Relu(),

stax.Dense(32), stax.Relu(),

stax.Dense(1)

)

kernel_NN4 = kernel_fn(dataX, dataX, ’ntk’)

init_fn, apply_fn, kernel_fn = stax.serial(

stax.Dense(32), stax.Relu(),

stax.Dense(32), stax.Relu(),

stax.Dense(32), stax.Relu(),

stax.Dense(32), stax.Relu(),

stax.Dense(32), stax.Relu(),

stax.Dense(1)

)

kernel_NN5 = kernel_fn(dataX, dataX, ’ntk’)

list_kernel_NN = [kernel_NN2, kernel_NN3, kernel_NN4, kernel_NN5]

# Normalization of the kernel matrix

for r in range(len(list_kernel_NN)):

for i in range(len(list_kernel_NN[r])):

for j in range(len(list_kernel_NN[r])):

list_kernel_NN[r][i][j] /= (list_kernel_NN[r][i][i] \

* list_kernel_NN[r][j][j]) ** 0.5

In order to predict the expectation value
$\operatorname{tr}(O\hat{\sigma}_{N}(x))$ of an observable $O$ for a new
ground state $\hat{\sigma}_{N}(x)$, we utilize the following property of
expectation values,

$\operatorname{tr}(O\hat{\sigma}_{N}(x))=\sum_{\ell=1}^{N}\kappa(x,x_{\ell})\operatorname{tr}(O\sigma_{T}(x_{\ell})).$
(35)

Hence, we first compute $\operatorname{tr}(O\sigma_{T}(x_{\ell}))$, which can
be done efficiently for $r$-body observables that factorize nicely into tensor
products. Indeed, an $O=O_{i_{1}}\otimes\ldots\otimes O_{i_{r}}$ ensures

$\operatorname{tr}(O\sigma_{T}(x_{\ell}))=\frac{1}{T}\sum_{t=1}^{T}\operatorname{tr}\left(O\sigma_{1}^{(t)}(x_{\ell})\otimes\dots\otimes\sigma_{n}^{(t)}(x_{\ell})\right)=\frac{1}{T}\sum_{t=1}^{T}\operatorname{tr}\left(O_{i_{1}}\sigma_{i_{1}}^{(t)}(x_{\ell})\right)\ldots\operatorname{tr}\left(O_{i_{r}}\sigma_{i_{r}}^{(t)}(x_{\ell})\right),$
(36)

and the right hand side can be computed with $\mathcal{O}(Tn)$ arithmetic
operations. Then, we can compute $\operatorname{tr}(O\hat{\sigma}_{N}(x))$ by
extrapolating $\operatorname{tr}(O\sigma_{T}(x_{\ell}))$ using
$\kappa(x,x_{\ell})$. We utilize scikit-learn, a Python package scikit-learn ,
for the training of these machine learning models.

Due to the different classical ML models one could consider (corresponding to
different regularization parameters $\lambda$ and kernel functions
$k(x,x^{\prime})$), we have to perform model selection to find an appropriate
ML model. Typically, the prediction performance will be quite sensitive to
these parameters, so one has to select them carefully. To evaluate the ML
models, we consider $100$ different points $x\in[-1,1]^{2}$ in parameter
space. Among these $100$ points, we select $N=20$ to be training data. These
are the circled points in Figure 2(b). For each property we would like to
predict, we choose one of the the three kernels and the different values of
$\lambda$ such that the prediction error is minimized on a validation set
containing $80-3$ inputs of $x$. The validation set is disjoint from the $20$
training points and the $3$ testing points for evaluating the prediction
performances (special markers in Figure 2(a)). Their purpose is to perform
model selection. Finally, we test on the three input $x$’s shown by the
special markers (cross, diamond and star) in Figure 2(b).

We found that for each property we would like to predict, the prediction
performance for different classical ML model varies moderately. When we have
sufficiently large training data size $N$, most choices of $\lambda$ and the
kernel function should yield good prediction performance. However, we are
using a very small number of training data in our experiments, hence the
choice of these options becomes more important. In particular, the best choice
of $\lambda$ can differ quite significantly over the different properties we
would like to predict.

Figure 10: Numerical experiment for predicting ground state properties
(Pauli-$X$ and $Z$ in each atom) in a 1D Rydberg atom system with 51 atoms
under different hyperparameters. (Left) The prediction error (root-mean-square
error) over different training sizes $N$ with a fixed number $T=10$ of
randomized Pauli measurements, also referred to as the shadow size. (Right)
The prediction error over different shadow sizes $T$ with a fixed training
data size $N=31$.

For completeness, we include a set of experiments where we vary the training
data size $N$ or the classical shadow size $T$, where by “shadow size” we mean
the number of randomized Pauli measurements used to approximate each state.
The result in given in Figure 10.  For this set of experiments, we consider a
fixed set of $70$ validation points in the phase space. Recall that we are
using the ML model to predict ground state properties. Here, we consider the
properties to be the expectation values of single-site Pauli-$X$ and Pauli-$Z$
operators. Because there are a total of $51$ atoms, there are a total of
$51\times 2=102$ properties. For each property, we randomly draw $10$
different points in the phase space (not in the training set or the validation
set). Therefore, the test set is of size $1020$, where each instance in the
test set corresponds to a property of a point in the phase space. The
prediction error is given by the root-mean-square error over the $1020$
instances in the test set. We can see that as training set size $N$ increases,
the prediction becomes better. However, we see that as the training size
increases, the slope of the prediction error (RMSE) over $N$ flattens. This is
expected from the theorem we established showing that
$N=m^{\mathcal{O}(1/\epsilon)}$, where $N$ is the training set size, $m$ is
the number of parameters, and $\epsilon$ is the prediction error. While the
theorem only provides an upper bound on $N$, if we assume the upper bound is
saturated, then we can use elementary calculus to derive

$\frac{d\epsilon}{dN}\mbox{ is proportional to
}-\frac{\epsilon^{2}}{N\log(m)}.$ (37)

Hence, the analysis is compatible with the observation that the slope of RMSE
over $N$ flattens as $N$ becomes larger.  While we proved a rigorous result
using the Dirichlet kernel, other commonly used ML models may yield a better
prediction performance in practice. Proving rigorous prediction guarantees and
understanding the limitations and strengths for other more commonly used ML
models are important future directions.

###  D.3 Ground state properties of the 2D antiferromagnetic Heisenberg model

Our next example is the two-dimensional antiferromagnetic Heisenberg model.
Spin-$\tfrac{1}{2}$ particles (i.e. qubits) occupy sites on a square lattice,
and for each pair $(ij)$ of neighboring sites the Hamiltonian contains a term
$J_{ij}\left(X_{i}X_{j}+Y_{i}Y_{j}+Z_{i}Z_{j}\right)$ where the couplings
$\\{J_{ij}\\}$ are uniformly sampled from the interval $[0,2]$. The parameter
$x$ is a list of all $J_{ij}$ couplings; hence in this case the dimension of
the parameter space is $m=O(n)$, where $n$ is the number of qubits. The
Hamiltonian $H(x)$ on a $5\times 5$ lattice is shown in Figure 3(a). The exact
ground state was found using DMRG. Analogously to the Rydberg atoms
experiments, we fixed the SVD cutoff to $10^{-8}$ and stopped the DMRG runs
when the difference in energy was below $10^{-4}$.

The classical ML models we considered are the same as the Rydberg atom chain
experiment. The only difference is that we slightly modify the Dirichlet
kernel (34b) to

$k(x,x^{\prime})=\sum_{i\neq
j}\sum_{k_{i}=-3}^{3}\sum_{k_{j}=-3}^{3}\cos\left(\pi(k_{i}(x_{i}-x_{i}^{\prime})+k_{j}(x_{j}-x_{j}^{\prime}))\right),\quad\mbox{(Dirichlet
kernel)}.$ (38)

We trained the classical ML model using a training set containing $N=90$
randomly chosen values of the parameter $x=\\{J_{ij}\\}$. Then, for each
property we would like to predict, we find the top-performing ML model setting
(out of all $\lambda$ parameters and kernel functions $k(x,x^{\prime})$) on a
validation set containing $100$ parameters $x$ distinct from the training set.
Finally, we test on $10$ newly sampled parameters $x$ to estimate the
prediction error. Figure 3(b) shows the prediction outcome from one of the
input parameter $x$. Figure 3(c) shows the RMSE from all $10$ input
parameters.

Figure 11: Numerical experiment for predicting ground state properties (two-
point correlation functions) in a 2D antiferromagnetic Heisenberg model with
$5\times 5$ spins under different hyperparameters. (Left) The predict error
(root-mean-square error) over different training size $N$ with a fixed number
of randomized Pauli measurements $T=10$, also referred to as the shadow size.
(Right) The prediction error (root-mean-square error) over different shadow
size $T$ with a fixed training data size $N=90$.

Similar to the Rydberg atom experiment, the best-performing ML model setting
differs across the properties we would like to predict. The three kernels
perform similarly at larger training data size $N$ and larger number of
randomized Pauli measurements $T$. But neural networks and Gaussian kernel
methods tend to perform better in most cases. The best choice of $\lambda$
differs substantially across the different properties: there is no single
choice of $\lambda$ that performs uniformly better than the other choices.

To showcase these effects, we also include a set of experiments where we vary
the training data size $N$ or the classical shadow size $T$, that is, the
number of randomized Pauli measurements used to approximate each state. The
numerical results are summarized in Figure 11. For this set of experiments, we
consider fixed sets of $100$ validation points.  For this set of experiments,
we consider a fixed set of $70$ validation points in the phase space. Recall
that we are using the ML model to predict ground state properties. Here, we
consider the properties to be the two-point correlation functions over every
pair of the $25$ spins. This results in a total of $25\times 25=625$
properties. For each property, we randomly draw $10$ testing points in the
$m=\mathcal{O}(n)$ dimensional parameter space (not in the training set or the
validation set). Therefore, the test set is of size $6250$, where each
instance in the test set corresponds to a property of a point in the parameter
space. The prediction error is given by the root-mean-square error over the
$6250$ instances in the test set. The results resemble what was found in the
Rydberg atom experiments, but with one notable difference — in the Rydberg
experiments, but not for the 2D antiferromagnet, the Dirichlet kernel has the
best performance for the largest shadow size $T$ we considered. This may be
because the dimension $m$ of the parameter space is much lower in the Rydberg
case.

###  D.4 Classifying phases of the bond-alternating XXZ model

To illustrate our classical ML for classifying quantum phases of matter, we
consider the bond-alternating XXZ model with $n=300$ spin-$\frac{1}{2}$
particles (i.e. qubits). The Hamiltonian is given by

$\sum_{i:\mathrm{odd}}J(X_{i}X_{i+1}+Y_{i}Y_{i+1}+\delta
Z_{i}Z_{i+1})+\sum_{i:\mathrm{even}}J^{\prime}(X_{i}X_{i+1}+Y_{i}Y_{i+1}+\delta
Z_{i}Z_{i+1}),$ (39)

and encompasses the bond-alternating Heisenberg model ($\delta=1$), as well as
the bosonic version of the Su-Schrieffer-Heeger model su1979solitons
($\delta=0$). The phase diagram in Figure 4(b) is obtained by evaluating the
partial reflection many-body topological invariant pollmann2012detection ;
elben2020many . It is given by

$\tilde{\mathcal{Z}}_{\mathcal{R}}=\frac{\mathcal{Z}_{\mathcal{R}}}{\sqrt{[\operatorname{tr}(\rho_{I_{1}}^{2})+\operatorname{tr}(\rho_{I_{2}}^{2})]/2}},\quad\mbox{where}\quad\mathcal{Z}_{\mathcal{R}}=\operatorname{tr}(\rho_{I_{1}\cup
I_{2}}\mathcal{R}_{I_{1}\cup I_{2}}),$ (40)

and we consider $I_{1}$ with $6$ spins: the $145$-th spin to the $150$-th
spin. Likewise, we fix $I_{2}$ to also contain $6$ spins: the $151$-th spin to
the $156$-th spin. Hence, the union $I_{1}\cup I_{2}$ contains $12$ spins. The
symbols $\rho_{I_{1}},\rho_{I_{2}}$ and $\rho_{I_{1}\cup I_{2}}$ denote the
reduced density matrices associated with each local region. The reflection
operator $\mathcal{R}_{I_{1}\cup I_{2}}$ acts on the local region $I_{1}\cup
I_{2}$ and is given by

$\mathcal{R}_{I_{1}\cup I_{2}}\ket{s_{1},\ldots,s_{|I_{1}\cup
I_{2}|}}=\ket{s_{|I_{1}\cup I_{2}|},\ldots,s_{1}},\quad\forall
s_{1},\ldots,s_{|I_{1}\cup I_{2}|}\in\\{0,1\\}.$ (41)

The partial reflection many body-topological invariant can resolve three
phases: trivial ($\tilde{\mathcal{Z}}_{\mathcal{R}}=+1$), symmetry-protected
topological (SPT) ($\tilde{\mathcal{Z}}_{\mathcal{R}}=-1$) and symmetry broken
($\tilde{\mathcal{Z}}_{\mathcal{R}}=0$). In Figure 4(b), we use the colors
blue (trivial), red (SPT) and gray (symmetry broken) to visualize these
different types of phases.

For each value of $J^{\prime}/J$ and $\delta$ considered, we construct the
exact ground state using DMRG, and find its classical shadow by performing
randomized single-qubit Pauli measurements a total of $T=500$ times. To
simulate this experiment, we follow the same setting for DMRG used in
elben2020many . We limit the maximum number of sweeps to $100$ and set the
DMRG cutoff to $10^{-9}$. We initialize the state to be the Néel state
$\ket{0101\ldots}$. To pin one of the degenerate ground state in the symmetry
broken phase, we include a penalty term given by $0.1JZ_{1}$ in the
Hamiltonian.

After obtaining the classical shadow representation $S_{T}(\rho_{\ell})$ for
each quantum state $\rho_{\ell}$, we compute the kernel matrix
$K\in\mathbb{R}^{N\times N}$, where each entry is given by the shadow kernel
$k^{\mathrm{(shadow)}}(S_{T}(\rho_{\ell}),S_{T}(\rho_{\ell^{\prime}}))$.
Recall that the shadow kernel is defined as

$k^{\mathrm{(shadow)}}(S_{T}(\rho),S_{T}(\tilde{\rho}))=\exp\left(\frac{1}{T^{2}}\sum_{t,t^{\prime}=1}^{T}\exp\left(\frac{1}{n}\sum_{i=1}^{n}\operatorname{tr}\left(\sigma_{i}^{(t)}\tilde{\sigma}_{i}^{(t^{\prime})}\right)\right)\right),\,\,\text{where}\,\,\sigma_{i}^{(t)}=3\lvert
s_{i}^{(t)}\rangle\\!\langle s_{i}^{(t)}\rvert-\mathbb{I},$ (42)

and the classical shadow representation is given by

$S_{T}(\rho)=\left\\{\ket{s_{i}^{(t)}}:\;i\in\left\\{1,\ldots,n\right\\},\;t\in\left\\{1,\ldots,T\right\\}\right\\},\quad\text{where}\quad\ket{s_{i}^{(t)}}\in\\{\ket{0},\ket{1},\ket{+},\ket{-},\ket{\mathrm{i}+},\ket{\mathrm{i}-}\\}\leavevmode\nobreak\
.$ (43)

Care should be taken when computing diagonal elements of the kernel matrix
$K$. The problem is that for $\rho=\tilde{\rho}$ and $t=t^{\prime}$, we
necessarily have
$\mathrm{tr}\left(\sigma_{i}^{(t)}\tilde{\sigma}_{i}^{(t)}\right)=5$ for all
$1\leq i\leq n$. And the double exponential will amplify this already
substantial contribution enormously. We found that counteracting this blow-up
improves the numerical stability of the kernel method substantially. When
$\ell=\ell^{\prime}$, when we compute
$k^{\mathrm{(shadow)}}(S_{T}(\rho_{\ell}),S_{T}(\rho_{\ell}))$, we sum over
$t\neq t^{\prime}$ instead of all $t,t^{\prime}$. In particular, when
$\rho=\tilde{\rho}$, we consider a slight modification to the kernel
definition,

$k^{\mathrm{(shadow)}}(S_{T}(\rho),S_{T}(\rho))=\exp\left(\frac{1}{T(T-1)}\sum_{t\neq
t^{\prime}}\exp\left(\frac{1}{n}\sum_{i=1}^{n}\operatorname{tr}\left(\sigma_{i}^{(t)}\sigma_{i}^{(t^{\prime})}\right)\right)\right),$
(44)

This modification also seems to slightly improve the classification
performance.

After evaluating the kernel matrix $K$, we renormalize the entries to obtain
the standardized kernel matrix

$\overline{K}_{\ell\ell^{\prime}}=\frac{K_{\ell\ell^{\prime}}}{\sqrt{K_{\ell\ell}K_{\ell^{\prime}\ell^{\prime}}}}\quad\text{for}\quad\ell,\ell^{\prime}\in\left\\{1,\ldots,N\right\\}.$
(45)

Subsequently, we perform kernel principal component analysis (PCA) on
$\overline{K}$. The implementation we used for kernel PCA is based on scikit-
learn sklearn_api . The output of kernel PCA is a list of low-dimensional
vectors (the dimension can be chosen arbitrarily, but we choose two dimensions
for this experiment). Each low-dimensional vector corresponds to a quantum
state. In Figure 4(c, d), we can see that the low-dimensional vectors are
clustered into different quantum phases of matter.

_Distinguishing an SPT phase from a trivial phase_ — We consider a one-
dimensional chain of $n=50$ qubits with $Z_{2}\times Z_{2}$ symmetry. The 1D
cluster state is in the nontrivial SPT phase. We generate other
representatives of the nontrivial SPT phase by applying symmetric low-depth
geometrically local random quantum circuits to the cluster state, and we
generate representatives of the trivial phase by applying symmetric random
circuits to a product state. We simulate the application of symmetric low-
depth geometrically local random quantum circuits to the cluster state through
matrix product states (MPS). Each circuit layer consists of patterns of random
two-qubit gates acting on next-to-nearest neighbors sites. We generate the
random gates in a block-sparse structure in the parity symmetry sectors. This
choice, together with the choice of connectivity, guarantees that the
$Z_{2}\times Z_{2}$ symmetry is conserved during the circuit evolution.

Randomized Pauli measurements are performed $T=500$ times to convert the
states to their classical shadows. We perform kernel PCA to find low-
dimensional representation for the quantum states using exactly the same
method as the experiment on bond-alternating XXZ model.

###  D.5 Distinguishing a topological phase from a trivial phase

We consider the task of distinguishing the toric code topological phase from
the trivial phase in a system of $n=200$ qubits. Kitaev’s toric code state
kitaev2003fault is in the nontrivial topologically-ordered phase, while a
product state represents the trivial phase. To populate both phases, we apply
low-depth geometrically local random Clifford circuits aaronson2004improved to
Kitaev’s toric code state kitaev2003fault with code distance 10, and we
generate representatives of the trivial phase by applying random Clifford
circuits to a product state. We utilize Clifford circuits to ensure efficient
simulation of in total $n=200$ qubits (and with a depth up to $9$) by means of
the Gottesman-Knill theorem. We again perform kernel PCA to find low-
dimensional representations for the quantum states using exactly the same
method as the experiment on bond-alternating XXZ model. This is used to
generate the plot in Figure 5(b) for a one-dimensional projection of the
feature space, as well as the plot in Figure 9(d) for a two-dimensional
projection.

For the unsupervised ML model shown in Figure 5(c), we consider a combination
of kernel PCA and randomized projections pmlr-v23-karnin12 . First we perform
kernel PCA to map the data to a six-dimensional subspace of the infinite-
dimensional feature space. Then we repeat the following procedure $500$ times.
We select a one-dimensional subspace uniformly at random in the six-
dimensional subspace. We project all the quantum states to the one-dimensional
subspace. Then, we find the center point (according to median instead of mean)
to split up the quantum states into two phases. We also record the sum of the
absolute values from all points to the center point in the one-dimensional
subspace. Finally, we consider the classification obtained from the random
one-dimensional projection that results in the largest sum of the absolute
values.

For the convolutional neural network (CNN) approach shown in Figure 5(c), we
consider the following CNN built from Keras chollet2015keras .

⬇

from tensorflow.keras import datasets, layers, models

model = models.Sequential()

model.add(layers.Conv2D(32, (2, 2), activation=’relu’, padding=’same’,

input_shape=(2*L, L, 6)))

model.add(layers.MaxPooling2D((2, 2)))

model.add(layers.Conv2D(32, (2, 2), activation=’relu’, padding=’same’))

model.add(layers.MaxPooling2D((2, 2)))

model.add(layers.Conv2D(32, (2, 2), activation=’relu’, padding=’same’))

model.add(layers.Flatten())

model.add(layers.Dense(32, activation=’relu’))

model.add(layers.Dense(2))

In the above code, $L$ is the code distance for the toric code and is equal to
$10$ in this experiment (recall that toric code ground state has $n=2L^{2}$
qubits). This CNN model is supervised and requires a training data with a
corresponding label for indicating which phase the training data point is in.
We first perform the Pauli-$6$ POVM on each qubit
carrasquilla2019reconstructing to transform the quantum state into a array of
size $n$ where each entry has six outcomes. We perform one-hot encoding to
yield a classical vector of size $6n$, where each entry in the classical
vector is either $0$ or $1$. Because the toric code ground state is two-
dimensional ($2L\times L$), we restructure the classical vector into a three-
dimensional tensor of size $2L\times L\times 6$. The first two dimensions
corresponds to the spatial dimension of the toric code ground state. The last
dimension corresponds to the one-hot encoded vector for the six-outcome POVM.
We then train the above model using the Adam optimizer kingma2014adam with the
categorical cross entropy as the loss function. The code is given below.

⬇

loss=tf.keras.losses.SparseCategoricalCrossentropy(from_logits=True),

metrics=[’accuracy’])

We train the convolutional neural network using $100$ training points (half
are topologically-ordered states, and the other half are trivial states). Then
we use a validation set of $100$ points to perform early stopping. This is
because the longer we train, the more likely the neural network is going to
overfit. Hence, it is a good practice to perform model selection by choosing
which model to use at different time points (during the training process). We
choose the model that performs the best on the validation set. Then we test
the classification accuracy (the percentage that the prediction of the phases
is correct) on a testing set consisting of $100$ points.

The performance of the above ML model is not substantially different from
random guessing. Hence, we also consider a very simple CNN enhanced with
classical shadow under $T=500$ randomized Pauli measurements. In particular,
we compute the local reduced density matrix using the classical shadow. Then
for each qubit, we represent it with the local reduced density matrix. For
simplicity, we consider the $i$-th qubit to be represented by a vector of size
$16$, which includes the $2$-body reduced density matrix for the subsystem
consisting of the $i$-th and the $i+1$-th qubit. Hence, each quantum state is
now represented by a classical vector of dimension $2L^{2}\times 16$. We
reshape the classical vector into a three-dimensional tensor of size $2L\times
L\times 16$. The classical vector is feed into the convolutional neural
network structured as follows. We also apply the Adam optimizer kingma2014adam
with the categorical cross entropy as the loss function. The evaluation
process is exactly the same as the CNN approach based on the Pauli-$6$ POVM.

⬇

from tensorflow.keras import datasets, layers, models

model = models.Sequential()

model.add(layers.Conv2D(16, (1, 1), activation=’relu’,\

padding=’same’, input_shape=(2*L, L, 16)))

model.add(layers.MaxPooling2D((2, 2)))

model.add(layers.Conv2D(16, (2, 2), activation=’relu’,\

padding=’same’))

model.add(layers.MaxPooling2D((2, 2)))

model.add(layers.Conv2D(16, (2, 2), activation=’relu’,\

padding=’same’))

model.add(layers.MaxPooling2D((2, 2)))

model.add(layers.Flatten())

model.add(layers.Dense(32, activation=’relu’))

model.add(layers.Dense(2))

model.compile(optimizer=’adam’,

loss=tf.keras.losses.SparseCategoricalCrossentropy(

from_logits=True),

metrics=[’accuracy’])

##  APPENDIX E Proof idea for the efficiency in predicting ground states

###  E.1 Main result

In order to illustrate the proof of Theorem 1, let us begin by looking at a
simpler task: training a machine learning model to predict a specified ground
state property instead of the classical representation of the ground state.
Consider the property $\operatorname{tr}(O\rho)$, where $\rho$ is the ground
state and $O$ is a local observable. In this simpler task, we consider the
training data to be

$\big{\\{}x_{1}\rightarrow\operatorname{tr}(O\rho(x_{1})),\quad\ldots,\quad
x_{N}\rightarrow\operatorname{tr}(O\rho(x_{N}))\big{\\}},$ (46)

where $x_{\ell}\in[-1,1]^{m}$ is a classical description of the Hamiltonian
$H(x_{\ell})$ and $\rho(x_{\ell})$ is the ground state of $H(x)$. Intuitively,
in a quantum phase of matter, the ground state property
$\operatorname{tr}(O\rho(x))$ changes smoothly as a function of the input
parameter $x$. The smoothness condition can be rigorously established as an
upper bound on the average magnitude of the gradient of
$\operatorname{tr}(O\rho(x))$ using quasi-adiabatic evolution
hastings2005quasiadiabatic ; bachmann2012automorphic , assuming that the
spectral gap of $H(x)$ is bounded below by a nonzero constant throughout the
parameter space. The upper bound on the average gradient magnitude enables us
to design a simple classical ML model based on an $l_{2}$-Dirichlet kernel for
generalizing from the training set to a new input $x\in\left[-1,1\right]^{m}$:

$\hat{O}_{N}(x)=\frac{1}{N}\sum_{\ell=1}^{N}\kappa(x,x_{\ell})\operatorname{tr}(O\rho(x_{\ell}))\,\,\,\mbox{with}\,\,\,\kappa(x,x_{\ell})=\sum_{k\in\mathbb{Z}^{m},\left\lVert
k\right\rVert_{2}\leq\Lambda}\cos(\pi k\cdot(x-x_{\ell}))\in\mathbb{R}.$ (47)

The $l_{2}$-Dirichlet kernel is often used in the study of high-dimensional
Fourier series weisz2012summability and the proposed ML model is equivalent to
learning a truncated Fourier series to approximate the function
$\operatorname{tr}(O\rho(x))$, where the parameter $\Lambda$ is a cutoff on
the wavenumber $k$ that depends on the upper bound on the gradient of
$\operatorname{tr}(O\rho(x))$. Using statistical analysis, one can guarantee
that
$\operatorname*{{\mathbb{E}}}_{x}|\hat{O}_{N}(x)-\operatorname{tr}(O\rho(x))|^{2}\leq\epsilon$
as long as the amount of training data $N=m^{\mathcal{O}(1/\epsilon)}$ where
our big-$\cal{O}$ notation is with respect to the $m\to\infty$ limit. Hence,
we can achieve a small _constant_ prediction error with an amount of training
data and computational time that are both polynomial in the number $m$ of
input parameters. The training is efficient because the number of modes needed
for the truncated Fourier series to provide an accurate approximation to
$\operatorname{tr}(O\rho)$ scales polynomially with $m$.

The key to the statistical analysis is to bound the model complexity of the
above machine learning model. In particular, the model complexity depends on
the number of wave vectors we consider in the $l_{2}$-Dirichlet kernel. The
more wave vectors $k$ we include, the higher the model complexity; and we
would have to use more data to train the ML model to achieve good
generalization performance. Furthermore, one could show that the amount of
data is proportional to the number of wave vectors we consider. In order to
achieve a prediction error
$\operatorname*{{\mathbb{E}}}_{x}|\hat{O}_{N}(x)-\operatorname{tr}(O\rho(x))|^{2}\leq\epsilon$,
we would need to select $\Lambda$ to be of order $\sqrt{1/\epsilon}$. Hence,
the number of wave vectors is proportional to the number of lattice points in
an $m$-dimensional $l_{2}$ ball of radius $\Lambda$. The volume of an
$m$-dimensional $l_{2}$ ball with radius $\Lambda$ is proportional to
$\Lambda^{m}=(1/\epsilon)^{m/2}$. If the number of lattices points is
proportional to the volume, then this would imply an exponential scaling in
the number of parameters $m$. However, through a proper combinatorial
analysis, we show that the number of lattices points is actually proportional
to $m^{\mathcal{O}(\Lambda^{2})}=m^{\mathcal{O}(1/\epsilon)}$, which is only
polynomial in the number of parameters $m$.

We can build on this idea to address the task of predicting ground state
representations. Now instead of predicting $\operatorname{tr}(O\rho)$ for a
new input $x$, the goal is to predict the classical shadow of the ground state
$\rho(x)$. We consider the training data to be
$\big{\\{}x_{\ell}\rightarrow\sigma_{1}(\rho(x_{\ell}))\big{\\}}_{\ell=1}^{N}$,
where $\sigma_{1}(\rho(x_{\ell}))$ is the classical shadow representation of
$\rho(x_{\ell})$ obtained from just a single randomized Pauli measurement of
the state (the $T=1$ case of Eq. (1)). Following the same approach as outlined
above for the case of predicting a single property, the predicted ground state
representation is now given by

$\hat{\sigma}_{N}(x)=\frac{1}{N}\sum_{\ell=1}^{N}\kappa(x,x_{\ell})\sigma_{1}(\rho(x_{\ell}))\,\,\,\mbox{with}\,\,\,\kappa(x,x_{\ell})=\sum_{k\in\mathbb{Z}^{m},\left\lVert
k\right\rVert_{2}\leq\Lambda}\cos(\pi k\cdot(x-x_{\ell}))\in\mathbb{R}.$ (48)

One can then guarantee that this representation accurately predicts
expectation values for a wide range of observables.

The fact that only a single snapshot $\sigma_{1}$ per parameter point is