plainlemlemma#1 plainlemLemma#1 plainthmtheorem#1 plainthmTheorem#1
plaincorcorollary#1 plaincorCorollary#1 plaindefidefinition#1
plaindefiDefinition#1 plainalgalgorithm#1plainalgAlgorithm#1
plainappappendix#1plainappAppendix#1

# Learning quantum many-body systems from a few copies

Cambyse Rouzé<EMAIL_ADDRESS>Department of Mathematics, Technische
Universität München, 85748 Garching, Germany  Daniel Stilck França
<EMAIL_ADDRESS>QMATH, Department of Mathematical Sciences, University of
Copenhagen, Denmark

###### Abstract

Estimating physical properties of quantum states from measurements is one of
the most fundamental tasks in quantum science. In this work, we identify
conditions on states under which it is possible to infer the expectation value
of all quasi-local observables of a given locality up to a relative error from
a number of samples that grows polylogarithmically with system size and
polynomially on the locality of the target observables. This constitutes an
exponential improvement over known tomography methods in some regimes. We
achieve our results by combining one of the most well-established techniques
to learn quantum states, the maximum entropy method, with techniques from the
emerging fields of quantum optimal transport and classical shadows. We
conjecture that our condition holds for all states exhibiting some form of
decay of correlations and establish it for several subsets thereof. These
include widely studied classes of states such as one-dimensional thermal and
gapped ground states and high-temperature Gibbs states of local commuting
Hamiltonians on arbitrary hypergraphs. Moreover, we show improvements of the
maximum entropy method beyond the sample complexity of independent interest.
These include identifying regimes in which it is possible to perform the
postprocessing efficiently and novel bounds on the condition number of
covariance matrices of many-body states.

##  I Introduction

The subject of quantum tomography has as its goal devising methods for the
efficient classical description of a quantum system from access to
experimental data. However, all tomographic methods for general quantum states
inevitably require resources that scale exponentially in the system size
Wright2016 ; Haah2017 , be it in terms of the number of samples required; or
the post-processing needed to perform the task.

Fortunately, most of the physically relevant quantum systems can be described
in terms of a (quasi)-local structure. These range from that of a local
interaction Hamiltonian corresponding to a finite temperature Gibbs state to
that of a shallow quantum circuit. Hence, locality is a physically motivated
requirement that brings the number of parameters describing the system to a
tractable number. Thus, effective tomographic procedures should be able to
incorporate this information. And, indeed, many protocols in the literature
achieve a good recovery guarantee in trace distance from a number of copies
that scales _polynomially_ with system size Cramer2010 ; Baumgratz_2013 ;
2004.07266 ; Torlai2018 ; PhysRevA.101.032321 ; 2004.07266 ; Eisert2020 .

Furthermore, in many scenarios one is interested in learning only physical
properties of the state on which tomography is being performed on. These are
often encoded into the expectation value of quasi-local observables that only
depend on reduced density matrices of subregions of the system. By Helstrom’s
theorem, obtaining a good recovery guarantee in trace distance is equivalent
to demanding that the expectation value of _all bounded observables_ are close
for the two states.

It is in turn natural to design tomographic procedures that can take advantage
of the fact that we only wish to approximate such observables, instead of
demanding a recovery in trace distance. For instance, the overlapping
tomography or classical shadows methods of Huang2020 ; cotler2020quantum ;
jena_pauli_2019 ; crawford_efficient_2020  allow for approximately learning
all $k$-local reduced density matrices of an $n$-qubit state with failure
probability $\delta$ using
$\mathcal{O}(e^{ck}k\log(n\delta^{-1})\epsilon^{-2})$ copies without imposing
any assumptions on the locality of the underlying state. This constitutes an
exponential improvement in the system size compared to the previously
mentioned many-body setting at the expense of an undesirable exponential
dependency in the locality of the observables. Unfortunately, without making
any assumptions on the underlying states, this exponential scaling is also
unavoidable Huang2020 .

In light of the previous discussion, it is natural to ask if, combining the
stronger structural assumptions required for many-body tomography with
techniques like classical shadows, it is possible to obtain the best of the
two worlds: a sample complexity that is logarithmic in system size and
polynomial in the locality of the observables we wish to estimate.

In this work, we achieve precisely this demand for a large class of physically
motivated states by combining recent techniques from quantum optimal transport
Carlen_2019 ; Rouz2019 ; gao_fisher_2020 ; palma_optimal ; kiani2021quantum ;
giacomo_cambyse , the well-established maximum entropy method
jaynes_information_1957  to learn quantum states and the method of classical
shadows. This results in an exponential improvement over known methods of
many-body tomography Cramer2010 ; Baumgratz_2013 ; 2004.07266 ; Torlai2018 ;
PhysRevA.101.032321 ; 2004.07266 ; Eisert2020  and recent shadow tomography or
overlapping tomography techniques Huang2020 ; cotler2020quantum ;
jena_pauli_2019 ; crawford_efficient_2020 , as summarized in Table 1.

Structure | Assumptions on State | Assumptions on Observable | Samples
---|---|---|---
Many-body tomography | Many-body | none | $\textrm{poly}(n,\epsilon^{-1})$
Classical Shadows | none |  $k$-local | $\textrm{poly}(e^{ck},\log(n),\epsilon^{-1})$
This work | Many body+transportation cost |  $k$-local | $\textrm{poly}(k,\log(n),\epsilon^{-1})$
Table 1: Summary of underlying assumptions and sample complexity of other
approaches to perform tomography on quantum many-body states.

Our revisit of the maximum entropy principle is further motivated by recent
breakthroughs in Hamiltonian learning 2004.07266 , shadow tomography Huang2020
, the understanding of correlations and computational complexity of quantum
Gibbs states kliesch_locality_2014 ; PhysRevLett.124.220601 ;
harrow2020classical ; kuwahara2020gaussian ; kuwahara_improved_2021  and
quantum functional inequalities capel2020modified ; giacomo_cambyse  that shed
new light on this seasoned technique. Our examples where we obtain exponential
improvements include thermal states, gapped ground states of $1D$ systems and
high-temperature thermal states of commuting Hamiltonians on arbitrary
hypergraphs. Furthermore, based on results by kuwahara2020gaussian , we
conjecture that our results should hold for any high-temperature Gibbs state,
even with long-range interactions, and give evidence in this direction.

The main ingredient to obtain our improvements are so-called transportation
cost (TC) inequalities Talagrand_1996 . They allow us to relate by how much
the expectation values of two states on Lipschitz observables, a concept we
will review shortly, differ to their relative entropy. Such inequalities
constitute a powerful tool from the theory of optimal transport
villani_optimal_2009  and are typical tools to prove sharp concentration
inequalities (raginsky_concentration_2014, , Chapter 3). Moreover, they have
been recently extended to quantum states palma_optimal ; Rouz2019 ;
giacomo_cambyse . By combining such inequalities with the maximum entropy
principle, we are able to easily control the relative entropy between the
states and, thus, the difference of expectation values of Lipschitz
observables.

Before we summarize our contributions in more detail, we will first define and
revise the main concepts required for our results, namely Lipschitz
observables, transportation cost inequalities and the basics of the maximum
entropy principle.

##### Lipschitz observables

In the classical setting, given a metric $d$ on a sample space $\mathcal{S}$,
the regularity of a function $f:\mathcal{S}\to\mathbb{R}$ can be quantified by
its _Lipschitz constant_ (raginsky_concentration_2014, , Chapter 3)

$\displaystyle\|f\|_{\operatorname{Lip}}=\sup_{x,y\in\mathcal{S}}|f(x)-f(y)|/d(x,y)\,.$
(1)

For instance, if we consider functions on the $n$-dimensional hypercube
$\\{-1,1\\}^{n}$ endowed with the Hamming distance, the Lipschitz constant
quantifies by how much a function can change per flipped spin. It should then
be clear that physical quantities like average magnetization have a small
Lipschitz constant. Some recent works Rouz2019 ; palma_optimal  extended this
notion to the noncommutative setting. Although the definition put forth in
Rouz2019  and the approach in palma_optimal  are not equivalent, they both
recover the intuition discussed previously. E.g., for the approach followed in
palma_optimal , the Lipschitz constant of an observable on $n$ qudits is
defined as 111the definition of palma_optimal  has a different normalization
and does not have the $\sqrt{n}$ term. This normalization will be convenient
to treat the constants of Rouz2019  and palma_optimal  on an equal footing.

$\displaystyle\frac{\|O\|_{\operatorname{Lip},\square}}{\sqrt{n}}:=\max\limits_{1\leq
i\leq n}\underset{\begin{subarray}{c}\rho,\sigma\in\mathcal{D}_{d^{n}}\\\
\operatorname{tr}_{i}[\rho]=\operatorname{tr}_{i}[\sigma]\end{subarray}}{\max}\,\operatorname{tr}\left[O(\rho-\sigma)\right]\,,$
(2)

where $\mathcal{D}_{d^{n}}$ denotes the set of $n$-qudit states. That is,
$\|O\|_{\operatorname{Lip},\square}$ quantifies the amount by which the
expectation value of $O$ changes for states that are equal when tracing out
one site. It is clear that $\|O\|_{\operatorname{Lip}}\leq
2\sqrt{n}\|O\|_{\infty}$ always holds by Hölder’s inequality, but it can be
the case that $\|O\|_{\operatorname{Lip},\square}\ll\sqrt{n}\|O\|_{\infty}$.
For instance, consider for some $k>0$ the $n$-qubit observable
$O=\sum_{i=1}^{n}\otimes_{j=i}^{i+k}Z_{j}$, where for each site $j$, $Z_{j}$
denotes the Pauli observable $Z$ acting on site $j$ and we take addition
modulo $n$. It is not difficult to see that
$\|O\|_{\operatorname{Lip},\square}=2k\sqrt{n}$, while $\|O\|_{\infty}=n$. See
the discussion in Fig. 1 for another example.

Moreover, one can show that shallow local circuits or short-time local
dynamics satisfying a Lieb-Robinson bound cannot substantially increase the
Lipschitz constant of an observable. That is, if we have that $\Phi_{t}^{*}$
is the quantum channel that describes some local dynamics at time/depth $t$ in
the Heisenberg picture and it satisfies a Lieb-Robinson bound, then we have:

$\displaystyle\|\Phi_{t}^{*}(O)\|_{\operatorname{Lip},\square}=\mathcal{O}(e^{\mu
t}\|O\|_{\operatorname{Lip},\square}).$

This result is discussed in more detail in Section B.1 of the supplemental
material. Thus, averages over local observables and short-time evolutions
thereof all belong to the class of observables that have a small Lipschitz
constant.

Another notion of a Lipschitz constant was put forth in Rouz2019 . Their
definition is of a more geometric flavour and is connected to so-called
differential structures related to a quantum state Carlen_2019 . Let
$a=\ket{1}\bra{0}$ and $a_{i}$ act like $a$ on qubit $i$ and identity else.
Then a differential Lipschitz constant can be defined as

$\displaystyle\|O\|_{\operatorname{Lip},\nabla}^{2}:=\sum_{i}\|[a_{i},O]\|_{\infty}^{2}+\sum_{i}\|[a_{i}^{\dagger},O]\|_{\infty}^{2}.$

It is not difficult to see that this definition exhibits a similar scaling for
local observables as the definition in Eq. (2).

The previous definitions allow us to quantify the regularity of all
observables, not only the strictly local ones considered so far. However, it
will also be fruitful to explicitly look at $(k,g)$-quasi local observables
kuwahara2020gaussian . Given a graph $V$ of cardinality $n$ corresponding to
the system $\mathcal{H}^{\otimes n}$, an observable $X\in\mathcal{M}_{d^{n}}$
is said to be $(k,g)$ quasi-local for $k\in\mathbbm{N}$ and $g>0$ if it can be
decomposed as

$\displaystyle X=\sum_{|A|\leq k}\,X_{A}\text{ with }\sum_{A\ni
v}\|X_{A}\|_{\infty}\leq g\text{ for }v\in V\,,$ (3)

i.e., each $X_{A}$ acts on at most $k$ systems and the norms of terms acting
on each system is bounded.

Once we are given a Lipschitz constant on observables or a set of quasi-local
observables, we can define a Wasserstein-1 distance on states by duality
Rouz2019 ; palma_optimal . The latter quantifies how well we can distinguish
two states by their action on certain regular or local observables
$\mathcal{L}_{\Lambda}$ and is given by

$\displaystyle{W}_{1,\Lambda}(\rho,\sigma):=\sup_{O\in\mathcal{L}_{\Lambda}}\operatorname{tr}\left[O(\rho-\sigma)\right].$
(4)

Here we will consider
$\mathcal{L}_{\Lambda}=\\{O=O^{\dagger},\|O\|_{\operatorname{Lip},\Lambda}\leq
1\\}$ for
$\|.\|_{\operatorname{Lip},\Lambda}\in\\{\|.\|_{\operatorname{Lip},\square},\|.\|_{\operatorname{Lip},\nabla}\\}$
and $\mathcal{L}_{\operatorname{Loc},(k,g)}$ as the set of $(k,g)$ quasi-local
observables as in Eq. (3).

Definition 4 is in direct analogy with the definition of the trace distance,
where we substituted the optimization over bounded observables by a more
structured set. We refer to Sec. B.1 for a thorough discussion of all these
definitions. We remark that $\mathcal{W}_{1,\operatorname{Loc}}$ corresponds
to the weakest of them, as it can only capture information contained in
$k$-body reduced density matrices, while the other definitions also capture
global information. On the other hand, $\mathcal{W}_{1,\operatorname{Loc}}$
allows for a more straightforward analysis, which will allow us to obtain
results for high-temperature Hamiltonians even with long-range interactions.

##### Transportation cost inequalities

The paragraphs before motivated the idea that observables with a small
Lipschitz constant capture most physically relevant observables and, thus,
that controlling the Wasserstein distance between two states gives rise to a
more physically motivated distance measure than the trace distance. However,
it is a priori not clear how to effectively control the Wasserstein distance
between states.

In this work, we will achieve this by relating Wasserstein distances to the
relative entropy between two states,
$D(\rho\|\sigma):=\operatorname{tr}\left[\rho(\log(\rho)-\log(\sigma))\right]$,
for $\sigma$ full-rank. This can be achieved through the notion of a
transportation cost inequality: a state $\sigma$ is said to satisfy a
transportation cost inequality with constant $\alpha>0$ if the Wasserstein
distance of $\sigma$ to any other state $\rho$ can be controlled by their
relative entropy, i.e.

$\displaystyle{W}_{1,\Lambda}(\rho,\sigma)\leq\sqrt{\frac{D(\rho\|\sigma)}{2\alpha}}$
(5)

holds for all states $\rho\in\mathcal{D}_{d^{n}}$.

Such inequalities are particularly powerful whenever the constant $\alpha$
does not depend on the system size $n$ or does so at most inverse
polylogarithmically, and can be thought of as a strengthening of Pinkser’s
inequality.

Transportation cost inequalities are closely related to the notion of Gaussian
concentration bobkov_exponential_1999 ; Rouz2019 ; palma_optimal , i.e.
Lipschitz functions strongly concentrating around their mean. Establishing
analogues of such concentration inequalities for quantum-many body systems has
also been a fruitful line of research in the last years and they are related
to fundamental questions in statistical physics, see e.g.
brandao_equivalence_2015 ; Anshu_2016 ; kuwahara2020gaussian ;
tasaki_local_2018 ; kuwahara_eigenstate_2020 . Although we are certain that
inequalities like Eq. (5) also shed new light on this matter, here we will
focus on their application to learning a classical description of a state
through maximum entropy methods. We refer to Table 2 for a summary of classes
of states known to satisfy it, as discussed in more detail below.

Recent works have established transportation cost inequalities with $\alpha$
either constant or logarithmic in system size for several classes of Gibbs
states of commuting Hamiltonians capel2020modified ; palma_optimal ;
giacomo_cambyse . They essentially are known to hold for local Hamiltonians on
arbitrary hypergraphs at high enough temperature or in $1D$. Moreover, these
works obtain the inequality for the stronger
$\|.\|_{\operatorname{Lip},\square}$ and $\|.\|_{\operatorname{Lip},\nabla}$
notions of Lipschitz constant. In this work we enlarge the class of examples
by showing them for outputs of short-depth circuits or evolutions. Moreover,
building upon the results of kuwahara2020gaussian , we show them also to hold
for high-temperature Gibbs states, even in the presence of long-range
interactions. We refer to Section B.1.3 for a precise definition and proof.
However, we only achieve this result for the
$\|\cdot\|_{\operatorname{Lip},\operatorname{Loc},(k,g)}$ definition and for
temperatures $\beta=\mathcal{O}((kg)^{-1})$. Although nontrivial, such an
inequality still has the undesired properties of only applying to strictly
local observables and the inverse temperature depending on the locality.
Despite these shortcomings, we believe that this preliminary result allows for
conjecturing that a transportation cost inequality might hold for such systems
for the stronger definitions of the Lipschitz constant and for all Gibbs
states above a threshold temperature that is independent of the locality of
the observables. Establishing such an inequality would allow us to apply our
result to a very broad range of physical systems.

##### Maximum-entropy methods

Let us now show how transportation cost inequalities can be combined with
maximum entropy methods. Such methods depart from the assumption that we are
given a set of self-adjoint, linearly independent Hermitian observables over
an $n$-qudit system,
$E_{1},\ldots,E_{m}\in\mathcal{M}_{d^{n}}^{\operatorname{sa}}$, with
$\|E_{i}\|_{\infty}\leq 1$, a maximal inverse temperature $\beta$ and the
promise that the state we wish to learn can be expressed as:

$\displaystyle\sigma\equiv\sigma(\lambda)=\frac{\operatorname{exp}\left(-\beta\sum\limits_{i=1}^{m}\lambda_{i}E_{i}\right)}{\mathcal{Z}(\lambda)}\,,$
(6)

where $\lambda\in\mathbbm{R}^{m}$ with sup norm
$\|\lambda\|_{\ell_{\infty}}\leq 1$ and

$\displaystyle\mathcal{Z}(\lambda)=\operatorname{tr}\left[\operatorname{exp}\left(-\beta\sum\limits_{i=1}^{m}\lambda_{i}E_{i}\right)\right]\,.$

Denoting by $e(\lambda)$ the vector of components
$e_{i}(\lambda)=\operatorname{tr}\left[E_{i}\,\sigma(\lambda)\right]$, we then
have that $\lambda$ is the unique optimizer of

$\displaystyle\min_{\begin{subarray}{c}\mu\in\mathbbm{R}^{m}\\\
\,\|\mu\|_{\ell_{\infty}}\leq
1\end{subarray}}\log(\mathcal{Z}(\mu))+\beta\sum\limits_{i=1}^{m}\mu_{i}e_{i}(\lambda)\,,$
(7)

which gives us a variational principle to learn the state given $e(\lambda)$.

Typical examples of observables $E_{i}$ are e.g. all $2-$local Pauli
observables corresponding to edges of a known graph. This models the situation
in which we know a state to be the thermal state of a Hamiltonian with a known
locality structure. More generally, for most of the examples discussed here
the $E_{i}$ are given as follows: we depart from a hypergraph $G=(V,E)$ and
assume that there is a maximum radius $r_{0}\in\mathbb{N}$ such that, for any
hyperedge $A\in E$, there exists a vertex $v\in V$ such that the ball
$B(v,r_{0})$ centered at $v$ and of radius $r_{0}$ includes $A$. The $E_{i}$
are then given by a basis of the traceless matrices on each hyperedge $A$.
This definition captures the notion of a local Hamiltonian w.r.t. the
hypergraph.

This framework also encompasses pure states after making an appropriate
approximation. For instance, in this article we will also consider the outputs
of shallow quantum circuits and unique ground states of gapped Hamiltonians,
which are pure. To see how to perform the approximation, suppose that
$\ket{\psi}=\mathcal{U}\ket{0}^{\otimes n}$ is the output of an unknown
shallow circuit $\mathcal{U}$ of depth $L$ with respect to a known interaction
graph $(V,E)$ on $n$ vertices. That is,

$\displaystyle\mathcal{U}=\prod_{\ell\in[L]}\,\bigotimes_{e\in
E_{\ell}}\,\mathcal{U}_{\ell,e},$

where each $E_{\ell}\subset E$ are a subset of the edges. That is, in this
setting the locality of the circuit is known, but the underlying unitary is
not. As we then show in Theorem B.3, the state $\ket{\psi}$ is $\epsilon$
close in Wasserstein distance to the state
$\sigma=\textrm{exp}(-\beta\sum\limits_{i=1}^{n}\mathcal{U}Z_{i}\mathcal{U}^{\dagger})/\mathcal{Z}$
for $\beta=\log(\epsilon^{-1})$. By a simple light-cone argument we can bound
the support of each $\mathcal{U}Z_{i}\mathcal{U}^{\dagger}$, as we know the
underlying structure of the circuit. Thus, to learn $\ket{\psi}$ up to a small
error, it is sufficient to learn $\sigma$ at temperature
$\mathcal{O}(\log(\epsilon^{-1}))$. A slightly different bound holds for
ground states of gapped Hamiltonians. They are well approximated in trace
distance by states whose inverse temperature scales like
$\mathcal{O}(\log(n\epsilon^{-1}))$ hastings_spectral_2006 . However, in light
of results like dalzell_locally_2019 , which assert that there exists a state
that approximates all $k$-local reduced density matrices of the ground state
up to and is the output of a ciricuit whose depth depends ony on $k$, we
believe that inverse temperatures scaling like
$\mathcal{O}(\log(\epsilon^{-1}))$ suffice to obtain a good approximation in
Wasserstein distance. Thus, when we say that a pure state satisfies a TC
inequality, we mean that a corresponding full-rank approximation of high
enough temperature does so, as otherwise the inequality is trivial.

We see from Eq. (7) that the expectation values of the $E_{i}$ completely
characterize the state $\sigma(\lambda)$. But it is possible to obtain a more
quantitative version of this statement through the following identity, also
observed in 2004.07266 :

$\displaystyle
D(\sigma(\mu)\|\sigma(\lambda))+D(\sigma(\lambda)\|\sigma(\mu))$
$\displaystyle\qquad\qquad=-\beta\,\langle\lambda-\mu|e(\lambda)-e(\mu)\rangle.$
(8)

In addition to showing that if $e(\lambda)=e(\mu)$ then
$\sigma(\mu)=\sigma(\lambda)$, Eq. (8) implies that by controlling how well
the local expectations values of one state approximate another, we can also
easily control their relative entropies. In particular, if $m=\mathcal{O}(n)$
and $\|e(\mu)-e(\lambda)\|_{\ell_{1}}=\mathcal{O}(\epsilon n)$ for some
$\epsilon>0$, we obtain from an application of Hölder’s inequality that:

$\displaystyle D(\sigma(\mu)\|\sigma(\lambda))$ $\displaystyle\leq
D(\sigma(\mu)\|\sigma(\lambda))+D(\sigma(\lambda)\|\sigma(\mu))$
$\displaystyle=\mathcal{O}(\beta\epsilon n)\,.$ (9)

We refer to Section A for more details. Thus, if we can find a state that
approximates the expectation value of each $E_{i}$ up to $\epsilon$, we are
guaranteed to have a $\mathcal{O}(\beta\epsilon)$ relative entropy density.
This observation is vital to ensure that the maximum entropy principle still
yields a good estimate of the state even under some noise in the vector of
expectation values $e(\lambda)$. Indeed, the variational principle of Eq. (7)
would allow us to recover the state exactly if we had access to the exact
values of $e(\lambda)$. However, it turns out that solving Eq. (7) with some
estimate $\hat{e}(\lambda)$ such that
$\|\hat{e}(\lambda)-e(\lambda)\|_{\ell_{\infty}}\leq\epsilon$ still yields a
Gibbs state $\sigma(\mu)$ satisfying Eq. (I). At this stage we will not
discuss the computational complexity of solving the maximum entropy problem,
which is deferred to the next section, but focus rather on the sample
complexity and recovery guarantees we obtain from the previous observations.

##### Combining TC with the maximum entropy principle

Suppose now that we have that each of the $E_{i}$ acts on at most $l_{0}$
qudits. Then, by using e.g. the method of classical shadows, we can estimate
the expectation values of all $E_{i}$ up to $\epsilon$ with failure
probability at most $\delta>0$ with
$\mathcal{O}(4^{l_{0}}\epsilon^{-2}\log(m\delta^{-1}))$ samples. For the rest
of this article we will set the probability of failure to
$\delta=\textrm{poly}(n)^{-1}$ and suppress it. From our discussion above, we
see that this is enough to obtain a state $\sigma(\mu)$ satisfying Eq. (I).
Further assuming that we have a TC with some constant $\alpha>0$ for
$\sigma(\lambda)$ we conclude that:

$\displaystyle\left|\operatorname{tr}\left[O(\sigma(\lambda)-\sigma(\mu))\right]\right|$
$\displaystyle\leq\|O\|_{\operatorname{Lip}}W_{1}(\sigma(\lambda),\sigma(\mu))$
$\displaystyle\leq\|O\|_{\operatorname{Lip}}\sqrt{\frac{D(\sigma(\mu)\|\sigma(\lambda))}{2\alpha}}$
$\displaystyle=\mathcal{O}(\sqrt{\beta\epsilon n}\|O\|_{\operatorname{Lip}}).$

Finally, recall that for sums of $k$ local operators like in Fig. 1, where we
have a $k=L^{2}$, the Lipschitz constant satisfies
$\|O\|_{\operatorname{Lip}}=\mathcal{O}(\sqrt{n})$ and we require a precision
of $\mathcal{O}(\epsilon n/k)$ to obtain a relative error of $\epsilon$.
Putting all these elements together, we conclude that by setting
$\epsilon=\beta\tilde{\epsilon}^{2}/(k^{2})$ for some $\tilde{\epsilon}>0$ we
arrive at

$\displaystyle\left|\operatorname{tr}\left[O(\sigma(\lambda)-\sigma(\mu))\right]\right|=\mathcal{O}(\tilde{\epsilon}k^{-1}n),$
(10)

which constitutes a relative error for the expectation value. In particular,
we see that the sample complexity required to obtain this error was

$\displaystyle\mathcal{O}(4^{r_{0}}k^{4}\beta^{2}\tilde{\epsilon}^{-4}\log(m)).$
(11)

We then obtain:

######  Theorem I.1 (Learning Gibbs states).

Let $\sigma(\lambda)$ be a Gibbs state as defined in Eq. (6) and such that
each $E_{i}$ acts on at most $l_{0}$ qudits. Moreover, suppose that
$\sigma(\lambda)$ satisfies a TC inequality with $\alpha$ depending at most
inverse logarithmically with system size. Then we can obtain a state
$\sigma(\mu)$ such that

$\displaystyle\left|\operatorname{tr}\left[O(\sigma(\lambda)-\sigma(\mu))\right]\right|\leq\mathcal{O}(\epsilon\sqrt{n}\,\|O\|_{\operatorname{Lip}})\,$
(12)

for all observables $O\in\mathcal{M}_{d^{n}}$ from
$\mathcal{O}(4^{l_{0}}\beta^{2}\textrm{poly}(\epsilon^{-1},log(m)))$ samples
of $\sigma(\lambda)$.

This gives an exponential speedup over known results for the tomography of
many-body states Cramer2010 ; Baumgratz_2013 ; Lanyon2017 ;
PhysRevA.101.032321 ; 2004.07266 , which have a polynomial sample complexity
to learn all $k$-body density matrices. Theorem I.1 also provides an
exponential improvement over shadow techniques in the locality of the
observables. However, unlike our methods, shadow techniques do not need to
make any assumptions on the underlying states.

We also remark that it is possible to improve the scaling in accuracy in Eq.
(11) from $\epsilon^{-4}$ to the expected $\epsilon^{-2}$. In order to do
that, it is important to bound the condition number of the Hessian of the log-
partition function. That is, we need a bound of the form

$\displaystyle LI\leq\nabla^{2}\log(\mathcal{Z}(\mu))\leq UI$ (13)

for constants $L,U>0$ and all $\mu\in B_{\ell_{\infty}}(0,1)$. We refer to
Sec. A of the supplemental material for a thorough discussion. We note that in
2004.07266  the authors show such results in a more general setting, although
with $U,L$ polynomial in $n$, which is not sufficient for our purposes. For us
it will be important to ensure that the condition number of the log-partition
function is at most polylogarithmic in system size (i.e.
$L^{-1}U=\tilde{\mathcal{O}}(1)$). It then follows that

$\displaystyle\|\lambda-\mu\|_{\ell_{2}}=\tilde{\mathcal{O}}(\|e(\lambda)-e(\mu)\|_{\ell_{2}}).$

That is, whenever the expectation values are close, the underlying parameters
must be close as well. In this case, we have from
$\|e(\lambda)-e(\mu)\|_{\ell_{2}}=\mathcal{O}(\epsilon\sqrt{m})$ that:

$\displaystyle-\beta\,\langle\lambda-\mu|e(\lambda)-e(\mu)\rangle$
$\displaystyle\qquad\leq\beta\|\lambda-\mu\|_{\ell_{2}}\|e(\lambda)-e(\mu)\|_{\ell_{2}}=\mathcal{O}(\beta\epsilon^{2}m).$

This quadratic improvement then yields the expected $\epsilon^{-2}$ scaling in
the accuracy. For commuting Hamiltonians we further show in Prop. E.1 and
Lemma E.1 of the supplementary material that:

######  Proposition I.1 (Strengthened strong convexity constant for commuting
Hamiltonians).

For each $\mu\in B_{\ell_{\infty}}(0,1)$, let $\sigma(\mu)$ be a Gibbs state
at inverse temperature $\beta$ corresponding to the commuting Hamiltonian
$H(\mu)=\sum_{j}\,\mu_{j}\,E_{j}$ on the hypergraph $G=(V,E)$, where
$\operatorname{tr}\left[E_{i}E_{j}\right]=0$ for all $i,j$ and each local
operator $E_{j}$ is traceless on its support. Then for $\beta$ such that the
states $\sigma(\mu)$ satisfy exponential decay of correlations, the Hessian of
the log-partition function is bounded by

$\displaystyle\mathcal{O}(\beta^{2})I\geq\nabla^{2}\,\log\,\mathcal{Z}(\mu)\geq\Omega(\beta^{2}\,e^{-c\beta})I.$

Structure | Samples | Postprocessing | Lipschitz-constant
---|---|---|---
Depth $l$ $2D$ circuit | $4^{2l}\epsilon^{-4}$ | $\operatorname{exp}(3n)$ | $\nabla,\square$
| Commuting Gibbs
---
for $\beta<\beta_{c}$
$\epsilon^{-2}$ | $\textrm{poly}(n)$ | $\nabla,\square$
$1D$ commuting | $\min\\{\epsilon^{-2}e^{\mathcal{O}(\beta)},\epsilon^{-4}\\}$ | $\textrm{poly}(n)$ | $\nabla,\square$
Gibbs for $\beta<\beta_{c}=\mathcal{O}((kg)^{-1})$ | $\epsilon^{-4}$ | $\textrm{poly}(n)$ | $\operatorname{loc}(k,g)$
Table 2: Performance of the algorithm under various assumptions of the
underlying state. All the estimates are up to $\textrm{polylog}(n)$ factors.
We refer to Sec. B.3 for proofs of the TC used and C for how to combine them
with maximum entropy methods. Figure 1: example of observable
$O=\sum_{i}O_{i}$ for $2D$ lattice system of size $n$. Each $O_{i}$ is
supported on a $L\times L$ square ($L=3$ in the figure). We have
$\|O\|_{\operatorname{Lip}}=\mathcal{O}(\sqrt{n})$ and $\|O\|=n/L^{2}$. Thus
our methods require $\textrm{poly}(L,\log(n),\epsilon^{-1})$ samples to
estimate the expectation value of all such observables. Shadow-like methods
require $\textrm{poly}(e^{cL},\log(n),\epsilon^{-1})$ samples, an
exponentially worse dependency in $L$. Even for moderate values of $L$, say
$L=5$, this can lead to $10^{7}$ factor savings sample complexity and gives an
exponential speedup for $L=\textrm{poly}(\log(n))$. Other many-body methods
have a $\textrm{poly}(L,n,\epsilon^{-1})$ Cramer2010 ; Baumgratz_2013 ;
Lanyon2017 ; PhysRevA.101.032321 ; 2004.07266  scaling, which in turn is
exponentially worse in the system size.

##  II Summary of the maximum entropy procedure and contributions

Now that we have discussed how our results yield better sample complexities
for some classes of states, we will discuss the maximum entropy algorithm in a
more detail and comment on how our results equip it with better performance
guarantees.

The flowchart in Figure 2 gives the general scheme behind the maximum entropy
method. Besides the exponential improvements in sample complexity laid out in
Table 2, we also provide structural improvements we elaborate on while also
explaining the general scheme:

##### Input:

The input consists of $m$ linearly independent operators $E_{i}$ of operator
norm $1$, some upper bound $\beta>0$, precision parameter $\epsilon>0$ and
step size $\eta$. Moreover, we are given the promise that the state of
interest satisfies (6). Although we will be mostly concerned with the case in
which these are local, we show the convergence of the algorithm in general in
Sec. A. The step size should be picked as $\eta=\mathcal{O}(U^{-1})$, as
explained in Sec. A.1.

##### Require:

We assume that we have access to copies of $\sigma(\lambda)$ and that we can
perform measurements to estimate the expectation values of the observables
$E_{i}$ up to precision $\epsilon>0$. For most of the examples considered
here, this will mostly require implementing simple, few-qudit measurements.

##### Output:

The output is in the form of a vector of parameters $\mu$ of a Gibbs state
$\sigma(\mu)$ as in Eq. (6). Note that unlike 2004.07266 , our goal is not to
estimate the vector of parameters $\lambda$, but rather obtain an
approximation of the state satisfying $\sigma(\lambda)\simeq\sigma(\mu)$. Here
we will focus on quantifying the output’s quality in relative entropy. More
precisely, the output of the algorithm is guaranteed to satisfy
$D(\sigma(\mu)\|\sigma(\lambda))=\mathcal{O}(\epsilon n)$.

##### Step 1:

In this step, we estimate the expectation values of each observable $E_{i}$ on
the state $\sigma(\lambda)$ up to an error $\epsilon$. The resources to be
optimized here are the number of samples of $\sigma(\lambda)$ we require and
the complexity to implement them. Using shadow tomography or Pauli grouping
methods Huang2020 ; cotler2020quantum ; bonet2020nearly  we can do so
requiring $\mathcal{O}(4^{r_{0}}\epsilon^{-2}\operatorname{polylog}(m))$
samples and Pauli or $1-$qubit measurements. This is discussed in more detail
in Sec. C.

##### Step 2:

The maximum entropy problem in Eq. (7) can be solved with gradient descent, as
it corresponds to a strictly convex optimization problem boyd2004convex . At
this step we simply initialize the algorithm to start at the maximally mixed
state.

##### Step 3:

It turns out that the gradient of the maximal entropy problem target function
at $\mu_{t}$ is proportional to $e(\mu_{t})-e(\lambda)$. Thus, to implement an
iteration of gradient descent, it is necessary to compute $e(\mu_{t})$, as we
assumed we obtained an approximation of $e(\lambda)$ in Step 1. Moreover, it
is imperative to ensure that the algorithm also converges with access to
approximations to $e(\mu)$ and $e(\lambda)$. This is because most algorithms
to compute $e(\mu_{t})$ only provide approximate values hastings_quantum_2007
; kliesch_locality_2014 ; PhysRevLett.124.220601 ; harrow2020classical  . In
addition, they usually have a poor scaling in the accuracy
PhysRevLett.124.220601 , making it necessary to show that the process
converges with rough approximations to $e(\mu_{t})$ and $e(\lambda)$. Here we
show that it is indeed possible to perform this step with only approximate
computations of expectation values. This allows us to identify classes of
states for which the postprocessing can be done efficiently. These results are
discussed in more detail in Sec. A.2.

##### Convergence loop:

Now that we have seen how to compute one iteration of gradient descent, the
next natural question is how many iterations are required to reach the
stopping criterium. As this is a strongly convex problem, the convergence
speed depends on the eigenvalues of the Hessian of the function being
optimized (boyd2004convex, , Section 9.1.2). For max-entropy, this corresponds
to bounding the eigenvalues of a generalized covariance matrix. In 2004.07266
the authors already showed such bounds for local Hamiltonians implying the
convergence of the algorithm in a number of steps depending polynomially in
$m$ and logarithmically on the tolerance $\epsilon$ for a fixed $\beta$. Here
we improve their bound in several directions. First, we show that the
algorithm converges after a polynomial in $m$ number of iterations for
arbitrary $E_{i}$, albeit with a polynomial dependence on the error, as
discussed in Sec. A.2. We then specialize to certain classes of states to
obtain various improvements. For high-temperature, commuting Hamiltonians we
provide a complete picture and show that the condition number of the Hessian
is constant in Prop. I.1. This implies that gradient descent converges in a
number of iterations that scales logarithmically in system size and error.

Figure 2: Flowchart for general maximum entropy algorithms.

##### Stopping condition and recovery guarantees:

the stopping condition,

$\displaystyle\|e(\lambda)-e(\mu_{t})\|_{\ell_{2}}\leq\epsilon\sqrt{m}\,,$

can be immediately converted to a relative entropy bound between the target
state and the current iterate by the identity (8). This justifies its choice
as a stopping criterion.

Since we already discussed the sample complexity of the maximum entropy
method, let us now discuss some of its computational aspects. There are two
quantities that govern the complexity of the algorithm: how many iterations we
need to perform until we converge and how expensive each iteration is.

As the maximum entropy problem is strongly convex, one can show that
$\mathcal{O}(UL^{-1}\log(n\epsilon^{-1}))$ suffice to converge. Here again
$U,L$ are bounds on the Hessian as in Eq. (13). Nevertheless, we also show how
to bypass requiring such bounds in Sec. A.1 and obtain that the maximum
entropy algorithm converges after $\mathcal{O}(m\epsilon^{-2})$ iterations
without any locality assumptions on the $E_{i}$ or strong convexity
guarantees. That is, the number of iterations is at most polynomial in $m$.

Let us now discuss the cost of implementing each iteration of the algorithm on
a classical computer. This boils down to estimating $e(\mu_{t})$ for the
current iterate, which can be achieved in various ways. In the worst case, it
is possible to just compute the matrix exponential and the expectation values
directly, which yields a complexity of $\mathcal{O}(d^{3n}m)$. However, for
many of the classes considered here it is possible to do this computation in
polynomial time. For instance, in PhysRevLett.124.220601  the authors show
that for high-temperature Gibbs states it is possible to approximate the
partition function efficiently. Thus, for the case of high-temperature Gibbs
states, not necessarily commuting ones, we can do the postprocessing
efficiently. It is also worth mentioning tensor network techniques to estimate
$e(\mu_{t})$. As we only require computing the expectation value of local
observables, recent works show that it is possible to find tensor network
states of constant bond dimension that approximate all expectation values of a
given locality well dalzell_locally_2019 ; alhambra_locally_2021 ;
huang_locally_2021 . From such a representation it is then possible to compute
$e(\mu_{t})$ efficiently in the $1D$ case by contracting the underlying tensor
network. Unfortunately, however, in higher dimensions the contraction still
takes exponential time. Table 2 provides a summary of the complexity of the
postprocessing for various classes.

It is also worth considering the complexity of the postprocessing with access
to a quantum computer, specially for commuting Hamiltonians. As all high-
temperature Gibbs states satisfy exponential decay of correlations, the
results of brandao_finite_2019  imply that high-temperature commuting Gibbs
states can be prepared with a circuit of depth that is logarithmic in system
size. Thus, by using the same method we used to estimate $e(\lambda)$ we can
also estimate $e(\mu_{t})$ by using the copies provided by the quantum
computer. The complexity of the postprocessing for shadows is linear in system
size. Thus, with access to a quantum computer we can perform the post-
processing for each iteration in time $\tilde{O}(m\epsilon^{-2})$. As in this
case we showed that the number of iterations is $\tilde{\mathcal{O}}(1)$, we
conclude that we can perform the postprocessing in time
$\tilde{\mathcal{O}}(m\epsilon^{-2})$. That is, for this class of systems our
results give an arguably complete picture regarding the postprocessing, as it
can be performed in a time comparable to writing down the vector of
parameters, up to polylogarithmic factors. Furthermore, given that the
underlying Gibbs states are known to satisfy TC and Prop. I.1 gives
essentially optimal bounds on the covariance matrices, we believe that the
present work essentially settles the question of how efficiently we can learn
such Hamiltonians and corresponding Gibbs states. We discuss this in more
detail in Sec. E.

Finally, an example of our bounds is illustrated in Fig. 3, where we show that
the number of samples required to estimate a local observable to relative
precision is essentially system-size independent.

Figure 3: Error in estimating a Lipschitz observable after performing the
maximum entropy reconstruction method. The underlying state is a classical
$1D$-Gibbs state with randomly chosen nearest-neighbor interactions and
inverse temperature $\beta=1$. We estimated all the $Z_{i}Z_{i+1}$ expectation
values from the original state based from $10^{3}$ samples of the original
state. We then computed the upper bound on the trace distance predicted by Eq.
8 and Pinsker’s inequality and compared it to the actual of discrepancy for a
Lipschitz observable on the reconstructed and actual state. The Lipschitz
observable was chosen as $\sum_{i}n^{-1}UZ_{i}Z_{i+2}U^{\dagger}$, where we
picked $U$ as a depth $3$ quantum circuit. We observe that the error incurred
is essentially independent of system size, and we get good predictions even
when the number of samples is smaller than it.

##  III Conclusion

In this article we have demonstrated that ideas from quantum optimal transport
can yield exponential improvements in the sample complexity required to learn
a quantum state when compared to state of the art methods. More precisely, we
showcased how the interplay between maximum entropy methods and the
Wasserstein distance, which is mediated by a transportation cost inequality,
allows for fine-tuning the complexity of observables whose expectation value
we wish to estimate and the number of samples required for that. Through our
techniques we essentially settled most questions related to how efficiently it
is possible to learn a commuting Gibbs state.

We believe that the framework we began to develop here will find applications
in other areas of quantum information theory such as in process tomography,
quantum machine learning and in quantum many-body physics. Indeed, in many
practical settings demanding a trace distance bound might be too restrictive,
and replacing it by a Wasserstein distance bound can lead to substantial
gains, as in this article.

Some of the outstanding open questions raised by this article are establishing
that a suitable notion of exponential decay of correlations in general implies
a transportation cost inequality and showing that TC holds for a larger class
of systems. Moreover, it would also be interesting to investigate other
applications of Gaussian concentration in many-body physics
brandao_equivalence_2015 ; Anshu_2016 ; kuwahara2020gaussian ;
tasaki_local_2018 ; kuwahara_eigenstate_2020  from the angle of transportation
cost inequalities.

##  IV Acknowledgements

DSF was supported by VILLUM FONDEN via the QMATH Centre of Excellence under
Grant No. 10059. The research of CR has been supported by project QTraj
(ANR-20-CE40-0024-01) of the French National Research Agency (ANR) and by a
Junior Researcher START Fellowship from the MCQST. DSF and CR are grateful to
Richard Kueng, Fernando Brandão and Giacomo De Palma for interesting
discussions.

## References

  * [1] Ryan O’Donnell and John Wright.  Efficient quantum tomography.  In Daniel Wichs and Yishay Mansour, editors, Proceedings of the 48th Annual ACM SIGACT Symposium on Theory of Computing, STOC 2016, Cambridge, MA, USA, June 18-21, 2016, pages 899–912. ACM, 2016.
  * [2] Jeongwan Haah, Aram Wettroth Harrow, Zhengfeng Ji, Xiaodi Wu, and Nengkun Yu.  Sample-optimal tomography of quantum states.  IEEE Trans. Inf. Theory, 63(9):5628–5641, 2017.
  * [3] Marcus Cramer, Martin B. Plenio, Steven T. Flammia, Rolando Somma, David Gross, Stephen D. Bartlett, Olivier Landon-Cardinal, David Poulin, and Yi-Kai Liu.  Efficient quantum state tomography.  Nature Communications, 1(1):149, dec 2010.
  * [4] T Baumgratz, A Nüßeler, M Cramer, and M B Plenio.  A scalable maximum likelihood method for quantum state tomography.  New Journal of Physics, 15(12):125004, dec 2013.
  * [5] Anurag Anshu, Srinivasan Arunachalam, Tomotaka Kuwahara, and Mehdi Soleimanifar.  Sample-efficient learning of interacting quantum systems.  Nature Physics, may 2021.
  * [6] Giacomo Torlai, Guglielmo Mazzola, Juan Carrasquilla, Matthias Troyer, Roger Melko, and Giuseppe Carleo.  Neural-network quantum state tomography.  Nature Physics, 14(5):447–450, may 2018.
  * [7] Jun Wang, Zhao-Yu Han, Song-Bo Wang, Zeyang Li, Liang-Zhu Mu, Heng Fan, and Lei Wang.  Scalable quantum tomography with fidelity estimation.  Phys. Rev. A, 101:032321, Mar 2020.
  * [8] Jens Eisert, Dominik Hangleiter, Nathan Walk, Ingo Roth, Damian Markham, Rhea Parekh, Ulysse Chabaud, and Elham Kashefi.  Quantum certification and benchmarking.  Nature Reviews Physics, 2(7):382–390, jul 2020.
  * [9] Hsin-Yuan Huang, Richard Kueng, and John Preskill.  Predicting many properties of a quantum system from very few measurements.  Nature Physics, June 2020.
  * [10] Jordan Cotler and Frank Wilczek.  Quantum overlapping tomography.  Physical Review Letters, 124(10):100401, 2020.
  * [11] Andrew Jena, Scott Genin, and Michele Mosca.  Pauli Partitioning with Respect to Gate Sets.  arXiv:1907.07859 [quant-ph], July 2019.  arXiv: 1907.07859.
  * [12] Ophelia Crawford, Barnaby van Straaten, Daochen Wang, Thomas Parks, Earl Campbell, and Stephen Brierley.  Efficient quantum measurement of Pauli operators in the presence of finite sampling error.  arXiv:1908.06942 [quant-ph], April 2020.  arXiv: 1908.06942.
  * [13] Eric A. Carlen and Jan Maas.  Non-commutative calculus, optimal transport and functional inequalities in dissipative quantum systems.  Journal of Statistical Physics, 178(2):319–378, nov 2019.
  * [14] Cambyse Rouzé and Nilanjana Datta.  Concentration of quantum states from quantum functional and transportation cost inequalities.  Journal of Mathematical Physics, 60(1):012202, January 2019.
  * [15] Li Gao, Marius Junge, and Nicholas LaRacuente.  Fisher Information and Logarithmic Sobolev Inequality for Matrix-Valued Functions.  Annales Henri Poincaré, 21(11):3409–3478, November 2020.
  * [16] Giacomo De Palma, Milad Marvian, Dario Trevisan, and Seth Lloyd.  The quantum Wasserstein distance of order 1, 2020.  arXiv:2009.04469v1.
  * [17] Bobak Toussi Kiani, Giacomo De Palma, Milad Marvian, Zi-Wen Liu, and Seth Lloyd.  Quantum earth mover’s distance: A new approach to learning quantum data.  arXiv preprint arXiv:2101.03037, 2021.
  * [18] Giacomo De Palma and Cambyse Rouzé.  Quantum concentration inequalities.  arXiv:2106.15819 [math-ph, physics:quant-ph], June 2021.  arXiv: 2106.15819.
  * [19] E. T. Jaynes.  Information Theory and Statistical Mechanics.  Physical Review, 106(4):620–630, May 1957.
  * [20] M. Kliesch, C. Gogolin, M.J. Kastoryano, A. Riera, and J. Eisert.  Locality of Temperature.  Physical Review X, 4(3):031019, July 2014.
  * [21] Tomotaka Kuwahara, Kohtaro Kato, and Fernando G. S. L. Brandão.  Clustering of conditional mutual information for quantum Gibbs states above a threshold temperature.  Phys. Rev. Lett., 124:220601, Jun 2020.
  * [22] Aram W Harrow, Saeed Mehraban, and Mehdi Soleimanifar.  Classical algorithms, correlation decay, and complex zeros of partition functions of quantum many-body systems.  In Proceedings of the 52nd Annual ACM SIGACT Symposium on Theory of Computing, pages 378–386, 2020.
  * [23] Tomotaka Kuwahara and Keiji Saito.  Gaussian concentration bound and ensemble equivalence in generic quantum many-body systems including long-range interactions.  Annals of Physics, 421:168278, 2020.
  * [24] Tomotaka Kuwahara, Álvaro M. Alhambra, and Anurag Anshu.  Improved Thermal Area Law and Quasilinear Time Algorithm for Quantum Gibbs States.  Physical Review X, 11(1):011047, March 2021.
  * [25] Ángela Capel, Cambyse Rouzé, and Daniel Stilck França.  The modified logarithmic sobolev inequality for quantum spin systems: classical and commuting nearest neighbour interactions.  arXiv preprint arXiv:2009.11817, 2020.
  * [26] M. Talagrand.  Transportation cost for gaussian and other product measures.  Geometric and Functional Analysis, 6(3):587–600, may 1996.
  * [27] Cédric Villani.  Optimal Transport, volume 338 of Grundlehren der mathematischen Wissenschaften.  Springer Berlin Heidelberg, Berlin, Heidelberg, 2009.
  * [28] Maxim Raginsky and Igal Sason.  Concentration of Measure Inequalities in Information Theory, Communications, and Coding.  Now Publishers, Norwell, MA, 2014.  OCLC: 1193114287.
  * [29] the definition of [16] has a different normalization and does not have the $\sqrt{n}$ term. This normalization will be convenient to treat the constants of [14] and [16] on an equal footing.
  * [30] S.G Bobkov and F Götze.  Exponential Integrability and Transportation Cost Related to Logarithmic Sobolev Inequalities.  Journal of Functional Analysis, 163(1):1–28, April 1999.
  * [31] Fernando G. S. L. Brandao and Marcus Cramer.  Equivalence of Statistical Mechanical Ensembles for Non-Critical Quantum Systems.  arXiv:1502.03263 [cond-mat, physics:quant-ph], February 2015.  arXiv: 1502.03263.
  * [32] Anurag Anshu.  Concentration bounds for quantum states with finite correlation length on quantum spin lattice systems.  New Journal of Physics, 18(8):083011, aug 2016.
  * [33] Hal Tasaki.  On the Local Equivalence Between the Canonical and the Microcanonical Ensembles for Quantum Spin Systems.  Journal of Statistical Physics, 172(4):905–926, August 2018.
  * [34] Tomotaka Kuwahara and Keiji Saito.  Eigenstate Thermalization from the Clustering Property of Correlation.  Physical Review Letters, 124(20):200604, May 2020.
  * [35] Matthew B. Hastings and Tohru Koma.  Spectral Gap and Exponential Decay of Correlations.  Communications in Mathematical Physics, 265(3):781–804, August 2006\.
  * [36] Alexander M. Dalzell and Fernando G. S. L. Brandão.  Locally accurate MPS approximations for ground states of one-dimensional gapped local Hamiltonians.  Quantum, 3:187, September 2019.
  * [37] B. P. Lanyon, C. Maier, M. Holzäpfel, T. Baumgratz, C. Hempel, P. Jurcevic, I. Dhand, A. S. Buyskikh, A. J. Daley, M. Cramer, M. B. Plenio, R. Blatt, and C. F. Roos.  Efficient tomography of a quantum many-body system.  Nature Physics, 13(12):1158–1162, September 2017.
  * [38] Xavier Bonet-Monroig, Ryan Babbush, and Thomas E O’Brien.  Nearly optimal measurement scheduling for partial tomography of quantum states.  Physical Review X, 10(3):031064, 2020.
  * [39] Stephen Boyd, Stephen P Boyd, and Lieven Vandenberghe.  Convex optimization.  Cambridge university press, 2004.
  * [40] M. B. Hastings.  Quantum belief propagation: An algorithm for thermal quantum systems.  Physical Review B, 76(20):201102, November 2007.
  * [41] Álvaro M. Alhambra and J. Ignacio Cirac.  Locally accurate tensor networks for thermal states and time evolution.  arXiv:2106.00710 [cond-mat, physics:quant-ph], June 2021.  arXiv: 2106.00710.
  * [42] Yichen Huang.  Locally accurate matrix product approximation to thermal states.  arXiv:2106.03854 [cond-mat, physics:math-ph, physics:quant-ph], June 2021.  arXiv: 2106.03854.
  * [43] Fernando G. S. L. Brandão and Michael J. Kastoryano.  Finite Correlation Length Implies Efficient Preparation of Quantum Thermal States.  Communications in Mathematical Physics, 365(1):1–16, January 2019\.
  * [44] Brian Swingle and Isaac H. Kim.  Reconstructing quantum states from local data.  Physical Review Letters, 113(26):260501, December 2014.  arXiv: 1407.2658.
  * [45] Christian Kokail, Rick van Bijnen, Andreas Elben, Benoît Vermersch, and Peter Zoller.  Entanglement Hamiltonian Tomography in Quantum Simulation, 2020\.  arXiv:2009.09000v1.
  * [46] M. B. Hastings.  Quantum belief propagation: An algorithm for thermal quantum systems.  Physical Review B, 76(20), nov 2007.
  * [47] Scott Aaronson.  Shadow tomography of quantum states.  In STOC’18—Proceedings of the 50th Annual ACM SIGACT Symposium on Theory of Computing, pages 325–338. ACM, New York, 2018.
  * [48] Akram Youssry, Christopher Ferrie, and Marco Tomamichel.  Efficient online quantum state estimation using a matrix-exponentiated gradient method.  New J. Phys., 21(3):033006, 2019.
  * [49] Fernando G. S. L. Brandão, Amir Kalev, Tongyang Li, Cedric Yen-Yu Lin, Krysta M. Svore, and Xiaodi Wu.  Quantum SDP solvers: large speed-ups, optimality, and applications to quantum learning.  In 46th International Colloquium on Automata, Languages, and Programming, volume 132 of LIPIcs. Leibniz Int. Proc. Inform., pages Art. No. 27, 14. Schloss Dagstuhl. Leibniz-Zent. Inform., Wadern, 2019.
  * [50] Fernando G. S. L. Brandão, Richard Kueng, and Daniel Stilck França.  Fast and robust quantum state tomography from few basis measurements.  arXiv:2009.08216 [quant-ph], September 2020.  arXiv: 2009.08216.
  * [51] Eric A Carlen and Jan Maas.  An analog of the 2-wasserstein metric in non-commutative probability under which the fermionic fokker–planck equation is gradient flow for the entropy.  Communications in mathematical physics, 331(3):887–926, 2014.
  * [52] Eric A Carlen and Jan Maas.  Gradient flow and entropy inequalities for quantum Markov semigroups with detailed balance.  Journal of Functional Analysis, 273(5):1810–1869, 2017.
  * [53] Maxim Raginsky and Igal Sason.  Concentration of measure inequalities in information theory, communications and coding.  Foundations and Trends in Communications and Information Theory; NOW Publishers: Boston, MA, USA, 2018.
  * [54] M. B. Hastings.  Locality in Quantum Systems.  arXiv:1008.5137 [math-ph, physics:quant-ph], August 2010.  arXiv: 1008.5137.
  * [55] David Poulin.  Lieb-Robinson Bound and Locality for General Markovian Quantum Dynamics.  Physical Review Letters, 104(19):190401, May 2010.
  * [56] Thomas Barthel and Martin Kliesch.  Quasilocality and efficient simulation of markovian quantum dynamics.  Physical review letters, 108(23):230504, 2012.
  * [57] Martin Kliesch, Christian Gogolin, and Jens Eisert.  Lieb-Robinson Bounds and the Simulation of Time-Evolution of Local Observables in Lattice Systems.  In Volker Bach and Luigi Delle Site, editors, Many-Electron Approaches in Physics, Chemistry and Mathematics, pages 301–318. Springer International Publishing, Cham, 2014.  Series Title: Mathematical Physics Studies.
  * [58] Nathael Gozlan and Christian Léonard.  Transport inequalities. a survey.  arXiv preprint arXiv:1003.3852, 2010.
  * [59] Salman Beigi, Nilanjana Datta, and Cambyse Rouzé.  Quantum reverse hypercontractivity: its tensorization and application to strong converses.  Communications in Mathematical Physics, 376(2):753–794, 2020.
  * [60] Roland L Dobrushin and Senya B Shlosman.  Completely analytical interactions: constructive description.  Journal of Statistical Physics, 46(5):983–1014, 1987.
  * [61] Huzihiro Araki.  Gibbs states of a one dimensional quantum lattice.  Communications in Mathematical Physics, 14(2):120–157, 1969.
  * [62] Marc Vuffray, Sidhant Misra, Andrey Lokhov, and Michael Chertkov.  Interaction screening: Efficient and sample-optimal learning of ising models.  In Advances in Neural Information Processing Systems, pages 2595–2603, 2016.
  * [63] Andrea Montanari et al.  Computational implications of reducing data to sufficient statistics.  Electronic Journal of Statistics, 9(2):2370–2390, 2015.
  * [64] Michael J Kastoryano and Fernando GSL Brandao.  Quantum Gibbs samplers: the commuting case.  Communications in Mathematical Physics, 344(3):915–957, 2016.
  * [65] Ivan Bardet, Angela Capel, and Cambyse Rouzé.  Approximate tensorization of the relative entropy for noncommuting conditional expectations.  arXiv preprint arXiv:2001.07981, 2020.

## Supplemental Material

This is the supplemental material to ”Learning many-body states from very few
copies”. We will start in Sec. A with a review of the basic properties of the
maximum entropy principle to learn quantum states. This is followed by a
discussion of Lipschitz constants, Wasserstein distances and transport cost
inequalities in Sec. B. After that, in Sec. C we more explicitely discuss the
interplay between the maximum entropy method and transportation cost
inequalities. We then briefly discuss scenarios in which the postprocessing
required for the maximum entropy method can be performed efficiently in Sec.
D. Finally, in Sec. E we discuss a class of examples were we show that all
technical results required to obtain the strongest guarantees of our work
hold: commuting Hamiltoniians at high enough temperature and 1D commuting
Hamiltonians.

We start by setting some basic notation. Throughout this article, we denote by
$\mathcal{M}_{k}$ the algebra of $k\times k$ matrices on $\mathbb{C}^{k}$,
whereas $\mathcal{M}_{k}^{\operatorname{sa}}$ denotes the subspace of self-
adjoint matrices. The set of quantum states over $\mathbb{C}^{k}$ is denoted
by $\mathcal{D}_{k}$. Typically, $k$ will be taken as $d^{n}$ for $n$ qudit
systems. The trace on $\mathcal{M}_{k}$ is denoted by $\operatorname{tr}$.
Given two quantum states $\rho,\sigma$, we denote by
$S(\rho)=-\operatorname{tr}\left[\rho\log(\rho)\right]$ the von Neumann
entropy of $\rho$, and by $D\left(\rho\|\sigma\right)$ be the relative entropy
between $\rho$ and $\sigma$, i.e.
$D(\rho\|\sigma)=\operatorname{tr}\left[\rho(\log(\rho)-\log(\sigma))\right]$
whenever the support of $\rho$ is contained in that of $\sigma$ and $+\infty$
otherwise. The trace distance will be denoted by
$\|\rho-\sigma\|_{tr}:=\operatorname{tr}\left[|\rho-\sigma|\right]$ and the
operator norm of an observable by $\|O\|_{\infty}$. Scalar products are
denoted by $\langle\cdot|\cdot\rangle$. Moreover, we will denote the
$\ell_{p}$ norm of vectors by $\|\cdot\|_{\ell_{p}}$ and for
$x\in\mathbbm{R}^{m},r\in\mathbbm{R}$ $B_{\ell_{p}}(x,r)$ will denote the
radius $r$ ball in the norm around $x$. The identity matrix is denoted by $I$.
The adjoint of an operator $A$ is denoted by $A^{\dagger}$ and that of a
channel $\Phi$ by $\Phi^{*}$.

##  Appendix A Maximum entropy principle for quantum Gibbs states

This work concerns the effectiveness of the maximum entropy method for the
tomography of quantum Gibbs states in various settings and regimes. Thus, we
start by recalling some basic properties of the maximum entropy method. Our
starting assumption is that the target state is well-described by a quantum
Gibbs state with respect to a known set of operators and that we are given an
upper bound on the inverse temperature:

######  Definition A.1 (Gibbs state with respect to observables).

Given a set of observables $\mathcal{E}=\\{E_{i}\\}_{i=1}^{m}$,
$E_{1},\ldots,E_{m}\in\mathcal{M}^{\operatorname{sa}}_{d^{n}}$ being linearly
independent with $\|E_{i}\|_{\infty}\leq 1$, we call a state
$\sigma\in\mathcal{D}_{d^{n}}$ a Gibbs state at inverse temperature $\beta>0$
if there exists a vector $\lambda\in\mathbbm{R}^{m}$ with
$\|\lambda\|_{\ell_{\infty}}\leq 1$ such that:

$\displaystyle\sigma=\operatorname{exp}\left(-\beta\sum\limits_{i=1}^{m}\lambda_{i}E_{i}\right)/\mathcal{Z}(\lambda),\qquad\text{where}\qquad\mathcal{Z}(\lambda)=\operatorname{tr}\left[\operatorname{exp}\left(-\beta\sum\limits_{i=1}^{m}\lambda_{i}E_{i}\right)\right]$
(14)

denotes the partition function. In what follows, we will denote $\sigma$ by
$\sigma(\lambda)$ and $\sum_{i}\lambda_{i}E_{i}=H(\lambda)$, where the
dependence of $\sigma(\lambda)$ on $\beta$ is implicitly assumed.

We are mostly interested in the regime where $m\ll d^{n}$. Then the above
condition can be interpreted as imposing that the matrix $\log(\sigma)$ is
sparse with respect to a known basis $\mathcal{E}$. A canonical example for
such states are Gibbs states of local Hamiltonians on a lattice , for which
$m=\mathcal{O}(n)$ and the observables $E_{i}$ are taken as tensor products of
Pauli matrices acting on neighboring sites. But we could also consider a basis
consisting of quasi-local operators or some subspace of Pauli strings.

Next, we review some basic facts about quantum Gibbs states. One of their main
properties is that they satisfy a well-known maximum entropy principle [19].
This allows us to simultaneously show that the expectation values of the
observables $\mathcal{E}$ completely characterize the state $\sigma(\lambda)$
and further provides us with a variational principle to learn a description
from which we can infer an approximation of other expectation values. Let us
start with the standard formulation of the maximum entropy principle:

######  Proposition A.1 (Maximum entropy principle).

Let $\sigma(\lambda)\in\mathcal{D}_{d^{n}}$ be a quantum Gibbs state (14) with
respect to the basis $\mathcal{E}$ at inverse temperature $\beta$ and
introduce $e_{i}(\lambda):=\operatorname{tr}\left[\sigma(\lambda)E_{i}\right]$
for $i=1,\ldots,m$. Then $\sigma(\lambda)$ is the unique optimizer of the
maximum entropy problem:

$\displaystyle\underset{\rho\in\mathcal{D}_{d^{n}}}{\operatorname{maximize}}$
$\displaystyle\quad S(\rho)$ (15) $\displaystyle\mathrm{subject~{}to}$
$\displaystyle\quad\operatorname{tr}\left[E_{i}\rho\right]=e_{i}(\lambda)\quad\text{for
all $i=1,\ldots,m$.}\quad$

Moreover, $\sigma(\lambda)$ optimizes:

$\displaystyle\min_{\begin{subarray}{c}\mu\in\mathbbm{R}^{m}\\\
\,\|\mu\|_{\ell_{\infty}}\leq
1\end{subarray}}\log(\mathcal{Z}(\mu))+\beta\sum\limits_{i=1}^{m}\mu_{i}e_{i}(\lambda)\,.$
(16)

###### Proof.

The proof is quite standard, but we include it for completeness. Note that for
any state $\rho\not=\sigma$ that is a feasible point of Eq. (15) we have that:

$\displaystyle S(\sigma(\lambda))-S(\rho)$
$\displaystyle=D(\rho\|\sigma(\lambda))+\operatorname{tr}\left[(\rho-\sigma(\lambda))\log(\sigma(\lambda))\right]$
$\displaystyle=D(\rho\|\sigma(\lambda))-\beta\sum\limits_{i=1}^{m}\lambda_{i}\operatorname{tr}\left[E_{i}\left(\rho-\sigma(\lambda)\right)\right]$
$\displaystyle=D(\rho\|\sigma(\lambda))>0\,,$

where we have used the fact that
$\operatorname{tr}\left[E_{i}\left(\rho-\sigma(\lambda)\right)\right]=0$ for
all feasible points and that the relative entropy between two different states
is strictly positive. This shows that $\sigma(\lambda)$ is the unique solution
of (15). Eq. (16) is nothing but the dual program of Eq. (15). ∎

Eq. (15) above gives a variational principle to find a quantum Gibbs state
corresponding to certain expectation values. As it is well-known, one can use
gradient descent to solve the problem in Eq. (16), as it is a strongly convex
problem. Various recent works have discussed learning of Gibbs states [44, 45,
5] and it is certainly not a new idea to do so through maximum entropy
methods. Nevertheless, we will discuss how to perform the postprocessing in
more detail, as some recent results allow us to give this algorithm stronger
performance guarantees. Finally, it should be said that although we draw
inspiration from [5], our main goal will be to learn a set of parameter
$\mu\in\mathbbm{R}^{m}$ such that the Gibbs states $\sigma(\mu)$ and
$\sigma(\lambda)$ are approximately the same on sufficiently regular
observables and mostly optimizing the sample complexity. This is in contrast
to the goal of [5], which was to learn the vector of parameters $\lambda$.
Learning $\lambda$ corresponds to a stronger requirement, in the sense that if
the vectors of parameters are close, then the underlying states are also
close, as made precise in Prop. A.2.

One of the facts we are going to exploit often is that it is possible to
efficiently estimate the relative entropy between two Gibbs states
$\sigma(\lambda),\sigma(\mu)$ given the parameters $\lambda,\mu$ and the
expectation values of observables in $\mathcal{E}$. This also yields an
efficiently computable bound on the trace distance. Indeed, as observed in
[5], we have that:

######  Proposition A.2.

Let $\sigma(\mu),\sigma(\lambda)\in\mathcal{D}_{d^{n}}$ be Gibbs states with
respect to a set of observables $\mathcal{E}$ at inverse temperature $\beta$.
Denote
$e(\lambda)=(\operatorname{tr}\left[\sigma(\lambda)E_{i}\right])_{i}\in\mathbbm{R}^{m}$.
Then

$\displaystyle\|\sigma(\mu)-\sigma(\lambda)\|_{tr}^{2}\leq
D(\sigma(\mu)\|\sigma(\lambda))+D(\sigma(\lambda)\|\sigma(\mu))=-\beta\,\langle\lambda-\mu|e(\lambda)-e(\mu)\rangle.$
(17)

###### Proof.

The equality in Eq. (17) follows from a simple manipulation. Indeed:

$\displaystyle
D(\sigma(\mu)\|\sigma(\lambda))+D(\sigma(\lambda)\|\sigma(\mu))=-\beta\sum\limits_{i=1}^{m}({\lambda_{i}-\mu_{i}})\operatorname{tr}\left[[\sigma(\mu)-\sigma(\lambda)]E_{i}\right].$

The bound on the trace distance then follows by applying Pinkser’s inequality.
∎

The statement of Proposition A.2 allows us to obtain quantitative estimates on
how well a given Gibbs state approximates another in terms of the expectation
values of known observables. In particular, a simple application of Hölder’s
inequality shows that if two Gibbs states are such that
$|\operatorname{tr}\left[E_{i}[\sigma(\mu)-\sigma(\lambda)]\right]|\leq\epsilon$,
then the sum of their relative entropies is at most

$\displaystyle\beta\,|\langle\lambda-\mu|e(\lambda)-e(\mu)\rangle|\leq\beta\,\|\lambda-\mu\|_{\ell_{\infty}}\,m\epsilon\leq
2m\epsilon\,\beta\,,$ (18)

where the outer bound arises from our assumption that
$\|\lambda\|_{\ell_{\infty}},\,\|\mu\|_{\ell_{\infty}}\leq 1$. Moreover, it is
straightforward to relate the difference of the target function in Eq. (16)
evaluated at two vectors to the difference of relative entropies between the
target state and their corresponding Gibbs states:

######  Lemma A.1.

Let $\sigma(\lambda)\in\mathcal{D}_{d^{n}}$ be a Gibbs state with respect to a
set of observables $\mathcal{E}$ at inverse temperature $\beta>0$ and define

$\displaystyle
f(\mu):=\log(\mathcal{Z}(\mu))+\beta\sum\limits_{i=1}^{m}e_{i}(\lambda)\,\mu_{i}\,.$

Then for any other two vectors $\mu,\xi\in\mathbb{R}^{m}$ with
$\|\mu\|_{\ell_{\infty}}\|\xi\|_{\ell_{\infty}}\leq 1$:

$\displaystyle
f(\mu)-f(\xi)=D(\sigma(\lambda)\|\sigma(\mu))-D(\sigma(\lambda)\|\sigma(\xi)).$

###### Proof.

The proof follows from straightforward manipulations. ∎

Thus, we see that a decrease of the target function $f$ when solving the max
entropy problem is directly connected to the decrease of the relative entropy
between the target state and the current iterate. We will later use this to
show the convergence of gradient descent for solving the max entropy problem
with arbitrary $\mathcal{E}$. However, before that we will discuss how the
convergence of the state $\sigma(\mu)$ to $\sigma(\lambda)$ is related to the
convergence of the parameters $\mu$ to $\lambda$.

###  A.1 Strong convexity and convergence guarantees

The maximum entropy problem (15) being a convex problem, it should come as no
surprise that properties of the Hessian of the optimized function are vital to
understanding its complexity and stability. For the maximum entropy problem,
the Hessian at a point is given by a generalized covariance matrix
corresponding to the underlying Gibbs state. As the results of [5] showcase,
the eigenvalues of such covariance matrices govern both the stability of Eq.
(16) with respect to $\mu$ and the convergence of gradient descent to solve
it. To see why, we recall some basic notions of optimization of convex
functions and refer to [39] for an overview.

######  Definition A.2 (Strong convexity).

Let $C\subset\mathbbm{R}^{m}$ be a convex set. A twice differentiable function
$f:C\to\mathbbm{R}$ is called strongly convex with parameters $U,L>0$ if we
have for all $x\in C$ that:

$\displaystyle UI\geq\nabla^{2}f(x)\geq LI.$

The optimization of strongly convex functions is well-understood. Indeed, we
have:

######  Proposition A.3.

Let $C\subset\mathbbm{R}^{m}$ be a convex set and $f:C\to\mathbb{R}$ be
strongly convex with parameters $L,U$ as in the definition above. Then, for
all $\epsilon>0$, the optimal value $\alpha:=\min_{x\in C}f(x)$ is achieved up
to error $\epsilon$ by the gradient descent algorithm initiated at $x^{0}\in
C$ with step size $U^{-1}$ after $S$ steps for

$\displaystyle
S\leq\frac{U}{L}\log\left(\frac{f(x^{0})-\alpha}{\epsilon}\right).$ (19)

The gradient norm satisfies $\|\nabla f(x^{k})\|_{\ell_{2}}^{2}\leq\delta$
after $S_{\nabla}$ steps with

$\displaystyle
S_{\nabla}\leq\frac{U}{L}\log\left(\frac{2L(f(x^{0})-\alpha)}{\delta}\right)\,.$

Moreover, we have for all $\mu,\lambda\in C$ that:

$\displaystyle\|\mu-\lambda\|_{\ell_{2}}\leq L^{-1}\|\nabla f(\mu)-\nabla
f(\lambda)\|_{\ell_{2}}.$ (20)

###### Proof.

These are all standard results which can be found e.g. in [39, Section 9.1.2
]. ∎

To see the relevance of these results for the maximum entropy problem, we
recall the following Lemma:

######  Lemma A.2.

Let $C=B_{\ell_{\infty}}(0,1)\subset\mathbbm{R}^{m}$, let
$\sigma(\lambda)\in\mathcal{D}_{d^{n}}$ be a Gibbs state with respect to a set
of operators $\mathcal{E}$ at inverse temperature $\beta$ and define
$f:C\to\mathbb{R}$ as:

$\displaystyle
f(\mu):=\log(\mathcal{Z}(\mu))+\beta\sum\limits_{i=1}^{m}e_{i}(\lambda)\,\mu_{i}.$
(21)

Then:

$\displaystyle\left(\nabla
f(\mu)\right)_{i}=\beta\operatorname{tr}\left[(\sigma({\lambda})-\sigma({\mu}))E_{i}\right]$
(22)

and

$\displaystyle\left(\nabla^{2}f(\mu)\right)_{ij}=\frac{\beta^{2}}{2}\operatorname{tr}\left[\\{E_{j},\Phi_{H(\mu)}(E_{i})\\}\sigma(\mu)\right]-\beta^{2}\,e_{i}(\mu)e_{j}(\mu),$
(23)

with

$\displaystyle\Phi_{H(\mu)}(E_{i})=\int\limits_{-\infty}^{+\infty}\nu_{\beta}(t)e^{-iH(\mu)t}E_{i}e^{iH(\mu)t}dt$

where $\nu_{\beta}(t)$ is a probability density function whose Fourier
transform is given by:

$\displaystyle\hat{\nu}_{\beta}(\omega)=\frac{2\tanh\left(\frac{\beta\omega}{2}\right)}{\beta\omega}\,.$

###### Proof.

The quantum belief propagation theorem [46] states that:

$\displaystyle\frac{\partial}{\partial\lambda_{i}}e^{-\beta
H(\lambda)}=-\frac{\beta}{2}\\{e^{-\beta
H(\lambda)},\Phi_{H(\lambda)}(E_{i})\\}.$

The claim then follows from a simple computation. ∎

Thus, we see that in order to compute the gradient of the target function $f$
for the maximum entropy problem, we simply need to compute the expectation
values of observables $\mathcal{E}$ on the current state and on the target
state. Moreover, the Hessian is given by a generalized covariance matrix of
the quantum Gibbs state. That this should indeed be interpreted as a
covariance matrix is most easily seen by considering commuting Hamiltonians.
Then indeed we have:

$\displaystyle\left(\nabla^{2}f(\mu)\right)_{ij}=\beta^{2}\left[\operatorname{tr}\left[\sigma(\mu)E_{i}E_{j}\right]-e_{i}(\mu)e_{j}(\mu)\right].$

For any Gibbs state it holds that:

######  Proposition A.4.

For all $\mu\in B_{\ell_{\infty}}(0,1)\subset\mathbbm{R}^{m}$, inverse
temperature $\beta>0$ and set of operators $\mathcal{E}$ of cardinality $m$,
we have:

$\displaystyle\nabla^{2}f(\mu)\leq 2\beta^{2}m\,I.$

###### Proof.

Note that

$\displaystyle\left|\operatorname{tr}\left[E_{i}\sigma(\mu)e^{iH(\mu)t}E_{j}e^{-iH(\mu)t}\right]\right|\leq
1,$

by Hölder’s inequality, the submultiplicativity and unitary invariance of the
operator norm and the fact that $\|E_{i}\|_{\infty}\leq 1$. Similarly, we have
that $|e_{i}(\mu)|,|e_{j}(\mu)|\leq 1$. Thus, by Lemma A.2,
$|\left(\nabla^{2}f(\mu)\right)_{ij}|\leq 2\beta^{2}$. As $\nabla^{2}f(\mu)$
is an $m\times m$ matrix, it follows from Gershgorin’s circle theorem that
$\nabla^{2}f(\mu)\leq 2\beta^{2}mI$. ∎

The proof above also showcases how exponential decay of correlations can be
used to sharpen estimates on the maximal eigenvalue of $\nabla^{2}f$, since in
that case $\left(\nabla^{2}f(\mu)\right)_{ij}$ will have exponentially
decaying entries. We discuss this more in detail when we focus on many-body
states, for which we also consider the more challenging question of lower
bounds.

###  A.2 Convergence with approximate gradient computation

Proposition A.3 already establishes the convergence of gradient descent
whenever we can compute the gradient exactly and have a bound on $L$.
Moreover, we see from Proposition A.2 that, in order to compute the gradient
of the function $f$ above, it suffices to estimate local expectation values.
Moreover, it is a standard result that gradient descent is strictly decreasing
for strictly convex problems [39].

However, in many settings it is only possible or desirable to compute the
expectation values of quantum Gibbs states approximately. Moreover, the
expectation values of the target state are only known up to statistical
fluctuations. It is then not too difficult to see that gradient descent still
offers the convergence guarantees if we only approximately compute the
gradient. We now state the exact convergence guarantees and precision
requirement for completeness.

######  Proposition A.5 (Convergence of gradient descent with constant
relative precision).

Let $\sigma(\lambda)\in\mathcal{D}_{d^{n}}$ be a quantum Gibbs state at
inverse temperature $\beta$ with respect to a set of operators $\mathcal{E}$,
and for a Gibbs state $\sigma(\mu)$ let $z(\mu)\in\mathbbm{R}^{m}$ be a vector
such that

$\displaystyle\|z(\mu)-\beta(e(\lambda)-e(\mu))\|_{\ell_{2}}\leq\,\frac{\beta}{4c}\|e(\mu)-e(\lambda)\|_{\ell_{\ell_{2}}}\,$
(24)

for some $c>10$. Then we have that:

$\displaystyle
D(\sigma(\lambda)\|\sigma(\mu-\tfrac{z(\mu)}{cU}))-D(\sigma(\lambda)\|\sigma(\mu))\leq-\frac{9\beta^{2}\|e(\mu)-e(\lambda)\|_{\ell_{2}}^{2}}{10\,c\,U}\,,$
(25)

where $U$ is a uniform bound on the operator norm of the Hessian of the
function $f$ defined in Equation (21).

###### Proof.

From a Taylor expansion and strong convexity we have for any two points
$\mu,\xi$ that:

$\displaystyle f(\xi)\leq f(\mu)+\langle\nabla
f(\mu)|(\xi-\mu)\rangle+\frac{U}{2}\|\xi-\mu\|_{\ell_{2}}^{2}.$

Note that $\nabla f(\mu)=\beta\left(e(\lambda)-e(\mu)\right)$. Setting
$\xi=\mu-\frac{z}{cU}=\mu+\frac{1}{cU}\left(-\nabla f(\mu)+\nabla
f(\mu)-z\right)$ we obtain:

$\displaystyle f\left(\mu-\frac{z}{cU}\right)$ $\displaystyle~{}~{}~{}\leq
f(\mu)-\frac{1}{cU}\|\nabla f\|_{\ell_{2}}^{2}+\frac{1}{cU}\langle\nabla
f(\mu)|\nabla f(\mu)-z\rangle+\frac{1}{2c^{2}U}\|-\nabla f(\mu)+\nabla
f(\mu)-z\|_{\ell_{2}}^{2}$ $\displaystyle~{}~{}~{}\leq
f(\mu)-\frac{1}{cU}\|\nabla f\|_{\ell_{2}}^{2}+\frac{1}{cU}\|\nabla
f(\mu)\|_{\ell_{2}}\,\|\nabla
f(\mu)-z\|_{\ell_{2}}+\frac{1}{2c^{2}U}\left(\|\nabla
f(\mu)\|_{\ell_{2}}+\|\nabla f(\mu)-z\|_{\ell_{2}}\right)^{2},$

where in the last step we used the Cauchy-Schwarz inequality. By our
assumption in Eq. (24) we have:

$\displaystyle f(\mu)$ $\displaystyle-\frac{1}{cU}\|\nabla
f\|_{\ell_{2}}^{2}+\frac{1}{cU}\|\nabla f(\mu)\|_{\ell_{2}}\|\nabla
f(\mu)-z\|_{\ell_{2}}+\frac{1}{2c^{2}U}\left(\|\nabla
f(\mu)\|_{\ell_{2}}+\|\nabla f(\mu)-z\|_{\ell_{2}}\right)^{2}$
$\displaystyle\leq f(\mu)-\frac{1}{cU}\|\nabla
f\|_{\ell_{2}}^{2}+\frac{1}{4c^{2}U}\|\nabla
f(\mu)\|_{\ell_{2}}^{2}+\frac{(1+(4c)^{-1})^{2}}{2c^{2}U}\|\nabla
f(\mu)\|_{\ell_{2}}^{2}$
$\displaystyle=f(\mu)-\frac{1}{cU}\left(1-\frac{1}{4c}-\frac{(1+(4c)^{-1})^{2}}{2c}\right)\|\nabla
f(\mu)\|_{\ell_{2}}^{2}$

and it can be readily checked that
$\frac{1}{4c}+\frac{(1+(4c)^{-1})^{2}}{2c}\leq\frac{1}{10}$ for $c\geq 10$. To
conclude the proof, note that by Lemma A.1:

$\displaystyle
D(\sigma(\lambda)\|\sigma(\mu-\tfrac{z}{cU}))-D(\sigma(\lambda)\|\sigma(\mu))=f(\mu-\tfrac{z}{cU}))-f(\mu)\,,$

and insert $\|\nabla
f(\mu)\|_{\ell_{2}}^{2}=\beta^{2}\|e(\mu)-e(\lambda)\|_{\ell_{2}}^{2}$. ∎

Thus, we see that we make constant progress for the gradient descent algorithm
if we only compute the derivative up to constant relative precision. We show
now how to pick our stopping criterium based on approximate computations of
the gradient which ensure convergence in polynomial time.

######  Proposition A.6.

Let $\sigma(\lambda)\in\mathcal{D}_{d^{n}}$ be a quantum Gibbs state at
inverse temperature $\beta$ with respect to a set of operators $\mathcal{E}$.
Suppose that at each time step $t$ of gradient descent we compute an estimate
$e^{\prime}(\mu_{t})$ of $e(\mu_{t})$ that satisfies

$\displaystyle\|e^{\prime}(\mu_{t})-e(\mu_{t})\|_{\ell_{\ell_{2}}}\leq\delta_{\mu}$

and set the stopping criterium to be:

$\displaystyle\|e(\lambda)-e^{\prime}(\mu_{t})\|_{\ell_{2}}<(4c+1)\delta_{\mu}$

for some constant $c>10$. Then gradient descent starting at $\mu=0$ with step
size $\frac{1}{cU}$ will converge to a state $\sigma(\mu_{*})$ satisfying:

$\displaystyle\|\sigma(\lambda)-\sigma(\mu_{*})\|_{tr}^{2}\leq
D(\sigma(\lambda)\|\sigma(\mu_{*}))+D(\sigma(\mu_{*})\|\sigma(\lambda))\leq
2(4c+1)\beta\delta_{\mu}\sqrt{m}$

after at most $\mathcal{O}(U\beta^{-2}n\log(d)\delta_{\mu}^{-2})$ iterations.

###### Proof.

First, we show that the relative precision bound required for Proposition A.5
holds under these assumptions. By our choice of the stopping criterium, at
each time step we have the property that, while we did not stop,

$\displaystyle\|e(\lambda)-e(\mu_{t})\|_{\ell_{2}}$
$\displaystyle=\|e(\mu)-e^{\prime}(\mu_{t})+e^{\prime}(\mu_{t})-e(\lambda)\|_{\ell_{2}}$
$\displaystyle\geq\|e^{\prime}(\mu_{t})-e(\lambda)\|_{\ell_{2}}-\|e(\mu_{t})-e^{\prime}(\mu_{t})\|_{\ell_{2}}$
$\displaystyle\geq(4c+1)\delta_{\mu}-\delta_{\mu}$
$\displaystyle=4c\delta_{\mu}$

by the reverse triangle inequality. As we assumed we have that
$\|e^{\prime}(\mu_{t})-e(\mu_{t})\|_{\ell_{2}}\leq\delta_{\mu}$, it follows
that
$(4c)^{-1}\|e(\lambda)-e(\mu_{t})\|_{\ell_{2}}\geq\|(e^{\prime}(\mu_{t})-e(\lambda))-(e(\mu_{t})-e(\lambda))\|_{\ell_{2}}$.
Multiplying the inequality by $\beta$, we see that the conditions of
Proposition A.5 are satisfied. Let us now show the convergence. By our choice
of initial point, we have that:

$\displaystyle D(\sigma(\lambda)\|\sigma(0))\leq n\log(d)\,.$

Now, suppose that we did not stop before $T$ iterations. It follows from a
telescopic sum and Proposition A.5 that:

$\displaystyle D(\sigma(\lambda)\|\sigma(\mu_{T}))\leq
n\log(d)-T\frac{9\beta^{2}(4c+1)^{2}\delta_{\mu}^{2}}{10\,c\,U}\,,$

since $\|e^{\prime}(\mu_{t})-e(\lambda)\|\geq(4c+1)\delta_{\mu}$ at all
iterations because we did not halt. As the relative entropy is positive, it
follows that $T=\mathcal{O}(\beta^{-2}Uc^{-1}n\log(d)\delta_{\mu}^{-2})$
before the stopping criterium is met. The recovery guarantee whenever the
stopping criterium is met follows from Proposition A.2 and the bound (18). ∎

Since we proved in Proposition A.4 that $U=\mathcal{O}(\beta^{2}m)$, it
follows that the number of iterations is $\mathcal{O}(nm)$. Thus, we see that
having a lower bound on the Hessian is not necessary to ensure convergence,
but it can speed it up exponentially:

######  Proposition A.7 (Exponential convergence of gradient descent with
approximate gradients).

In the same setting as Proposition A.5 we have:

$\displaystyle
f\left(\mu-\frac{z(\mu)}{cU}\right)-f(\lambda)\leq\left(1-\frac{18L}{10cU}\right)(f(\mu)-f(\lambda)).$
(26)

In particular, gradient descent with approximate gradient computations
starting at $\mu_{0}=0$ converges after
$\mathcal{O}\left(\frac{U}{L}\log(m\epsilon^{-1})\right)$ iterations to $\mu$
such that $f(\mu)-f(\lambda)\leq\epsilon$.

###### Proof.

For any strongly convex function $f$ and points $\mu,\xi\in C$ we have that:

$\displaystyle f(\xi)\geq f(\mu)+\langle\nabla
f(\mu)|(\xi-\mu)\rangle+\frac{L}{2}\|\mu-\xi\|_{\ell_{2}}^{2}.$

As explained in [39, chapter 9], the R.H.S. of the equation above is a convex
quadratic function of $\xi$ for $\mu$ fixed. One can then easily show that its
minimum is achieved at $\tilde{\xi}=\mu-\frac{1}{m}\nabla f(\mu)$. From this
we obtain:

$\displaystyle f\left(\mu\right)-f(\lambda)\geq-\frac{1}{2L}\|\nabla
f(\mu)\|_{\ell_{2}}^{2}=-\frac{\beta^{2}\|e(\mu)-e(\lambda)\|_{\ell_{2}}^{2}}{2L}\,,$

where the last identity follows from Eq. (22). By subtracting $f(\lambda)$
from both sides of the inequality Eq. (25) of Proposition A.5 and rearranging
the terms we have that:

$\displaystyle f\left(\mu-\frac{z(\mu)}{cU}\right)-f(\lambda)$
$\displaystyle\leq
f(\mu)-f(\lambda)-\frac{9\beta^{2}\|e(\mu)-e(\lambda)\|_{\ell_{2}}^{2}}{10cU}$
$\displaystyle\leq
f\left(\mu\right)-f(\lambda)-\frac{18L}{10cU}(f\left(\mu\right)-f(\lambda))$
$\displaystyle=\left(1-\frac{18L}{10cU}\right)(f\left(\mu\right)-f(\lambda))$

This yields the claim in Eq. (26). To obtain the second claim, note that
applying Eq. (26) iteratively yields that after $k$ iterations we have that:

$\displaystyle\left(f\left(\mu_{k}\right)-f(\lambda)\right)\leq\left(1-\frac{18L}{10cU}\right)^{k}\left(f(0)-f(\lambda)\right).$

By our choice of initial point we have that $f(0)-f(\lambda)=\mathcal{O}(m)$,
which yields the claim solving for $k$ and noting that
$-\log\left(1-\frac{18L}{10cU}\right)=\mathcal{O}(\frac{U}{L})$. ∎

######  Remark A.1 (Comparison to mirror descent).

It is also worth noting that the convergence guarantees of Proposition A.6 and
update rules of gradient descent are very similar to the ones of mirror
descent with the von Neumann entropy as potential, another algorithm used for
learning of quantum states [47, 48, 49, 50]. In this context, mirror descent
would use a similar update rule. However, instead of computing the whole
gradient, i.e. all expectation values of the basis, for one iteration, mirror
descent just requires us to find one $i$ such that
$|e_{i}(\lambda)-e_{i}(\mu)|\geq\delta$ and updates the Hamiltonian in the
direction $i$. This implies that the algorithm can be run online while we
still estimate some other $e_{i}$, but we will not analyse this variation more
in detail here.

Finally, we assumed so far that we knew the expectation values of the target
state, $e(\lambda)$, exactly. However, it follows straightforwardly from
Proposition A.6 that knowing each expectation value up to an error $\epsilon$
is sufficient to ensure that the additional error due to statistical
fluctuations is at most $\epsilon m$. More precisely, if we have that
$\|e(\lambda)-e(\lambda^{\prime})\|_{\infty}\leq\epsilon$ for some
$\epsilon>0$, then any Gibbs state satisfying
$\|e(\mu_{*})-e(\lambda^{\prime})\|_{\ell_{2}}\leq\delta$ satisfies:

$\displaystyle
D(\sigma(\lambda)\|\sigma(\mu_{*}))+D(\sigma(\mu_{*})\|\sigma(\lambda))\leq
2\beta\delta\sqrt{m}+2\beta\epsilon m$

by Proposition A.2 and a Cauchy-Schwarz inequality. In short, we obtain:

######  Theorem A.1 (Computational complexity).

Let $\sigma(\lambda)\in\mathcal{D}_{d^{n}}$ be a quantum Gibbs state at
inverse temperature $\beta$ with respect to a set of operators $\mathcal{E}$
and $C_{e}$ be the computational cost of computing $e^{\prime}(\mu)$
satisfying

$\displaystyle\|e^{\prime}(\mu_{t})-e(\mu_{t})\|_{\ell_{2}}\leq\delta_{\mu}$

for $\mu\in B_{\ell_{\infty}}(0,1)$ and $\delta_{\mu}>0$. Then gradient
descent starting at $\mu=0$ with step size $\frac{1}{cU}$ converges to a state
$\sigma(\mu_{*})$ satisfying:

$\displaystyle\|\sigma(\lambda)-\sigma(\mu_{*})\|_{tr}^{2}\leq
D(\sigma(\lambda)\|\sigma(\mu_{*}))+D(\sigma(\mu_{*})\|\sigma(\lambda))\leq
2(4c+1)\beta\delta_{\mu}\sqrt{m}$

in

$\displaystyle\mathcal{O}\left(\min\left\\{UC_{e}\beta^{-2}n\log(d)\delta_{\mu}^{-2},\frac{UC_{e}}{L}\log(m\epsilon^{-1})\right\\}\right)$

time.

###### Proof.

Clear. ∎

In short, we see that we can perform the recovery by simply computing the
gradient approximately. In particular, as already hinted at in [5], this
implies that recent methods developed to approximately compute the partition
function of high-temperature quantum Gibbs states can be used to perform the
postprocessing in polynomial time [20, 21, 22, 23].

##  Appendix B Lipschitz constants and transportation cost inequalities

In this section we show one of the main results of this paper. In words, it
identifies conditions under which it is possible to estimate all expectation
values of $k$-local observables up to an error $\epsilon$ by measuring
$\mathcal{O}(\textrm{poly}(k,\log(n),\epsilon^{-1}))$ copies of it, where $n$
is the system size, an exponential speedup in some regimes. To obtain this
result, we will combine the maximum entropy method with techniques from
quantum optimal transport. In order to formalize and prove the result claimed,
we resort to transportation cost inequalities and the notion of Lipschitz
constants of observables, which we now introduce.

###  B.1 Lipschitz constants and Wasserstein metrics

Transportation cost inequalities, introduced by Talagrand in the seminal paper
[26], constitute one of the strongest tools available to show concentration of
measure inequalities. In the quantum setting, their study was initiated in
[51, 52, 14, 16]. Here we are going to show how they can also be used in the
context of quantum tomography. On a high level, a transportation cost
inequality for a state $\sigma$ quantifies by how much the relative entropy
with respect to another state $\rho$ is a good proxy to estimate to what
extent the expectation values of sufficiently regular observables differ on
the states. As maximum entropy methods allow for a straightforward control of
the convergence of the learning in relative entropy (cf. Section A), they can
be combined to derive strong recovery guarantees. But first we need to define
what we mean by a regular observable.

We start by a short discussion of Lipschitz constants and the Wasserstein-1
distance. To obtain an intuitive grasp of these concepts, one way is to first
recall the variational formulation of the trace distance of two quantum states
$\sigma,\rho$:

$\displaystyle\|\rho-\sigma\|_{tr}=\sup\limits_{P=P^{\dagger},\|P\|_{\infty}\leq
1}\operatorname{tr}\left[P(\rho-\sigma)\right].$

Seeing probability distributions as diagonal quantum states, we recover the
variational formulation of the total variation distance by noting that we may
restrict to diagonal operators $P$. Thus, the total variation distance
quantifies by how much the expectation values of arbitrary bounded functions
can differ under the two distributions. However, in many situations we are not
interested in expectation values of arbitrary bounded observables, but rather
observables that are sufficiently regular. E.g., most observables of physical
interest are (quasi)-local. Thus, it is natural to look for distance measures
between quantum states that capture the notion that two states do not differ
by much when restricting to expectation values of sufficiently regular
observables. These concerns are particularly relevant in the context of
tomography protocols, as they should be designed to efficiently obtain a state
that reflects the expectation values of physically motivated observables.

As we will see, one of the ways of ensuring that the sample complexity of the
tomography algorithm reflects the regularity of the observables we wish to
recover is through demanding a good recovery in the Wasserstein distance of
order $1$ [14, 16].

In the classical setting [28, Chapter 3], one way to define a Wasserstein-1
distance between two probability distributions is by replacing the
optimization over all bounded diagonal observables by that over those that are
sufficiently regular: given a metric $d$ on a sample space $\Omega$, we define
the Lipschitz constant of a function $f:\Omega\to\mathbbm{R}$ to be:

$\displaystyle\|f\|_{\operatorname{Lip}}:=\sup\limits_{x,y\in\Omega}\frac{|f(x)-f(y)|}{d(x,y)}\,.$

Denoting the Wasserstein-1 distance by $W_{1}$, it is given for two
probability measures $p,q$ on $\Omega$ by

$\displaystyle W_{1}(p,q):=\sup_{f:\|f\|_{\operatorname{Lip}}\leq
1}|\mathbb{E}_{p}(f)-\mathbb{E}_{q}(f)|\,.$ (27)

That is, this metric quantifies by how much the expectation values of
sufficiently regular functions can vary under $p$ and $q$, in clear analogy to
the variational formulation of the trace distance. We refer to [28, 27] for
other interpretations and formulations of this metric.

####  B.1.1 Quantum Hamming Wasserstein distance

It is not immediately clear how to generalize these concepts to noncommutative
spaces. There are by now several definitions of transport metrics for quantum
states which do not seem to be equivalent [14, 16, 17]. As already noted in
the main text, de Palma et al. defined the Lipschitz constant of an observable
$O\in\mathcal{M}_{d^{n}}$ as :

$\displaystyle\|O\|_{\operatorname{Lip},\square}=\sqrt{n}\max\limits_{1\leq
i\leq n}\underset{\begin{subarray}{c}\rho,\sigma\in\mathcal{D}_{d^{n}}\\\
\operatorname{tr}_{i}[\rho]=\operatorname{tr}_{i}[\sigma]\end{subarray}}{\max}\,\operatorname{tr}\left[O(\rho-\sigma)\right]\,.$
(28)

That is, the Lipschitz constant quantifies how much the expectation value of
an observable can change when evaluated on two states that only differ on one
site. This is in analogy with the Lipschitz constants induced by the Hamming
distance on the hypercube, so we denote it with $\square$. Note that in our
definition we added the $\sqrt{n}$ factor, which will turn out to be
convenient later. Armed with this definition, we can immediately obtain an
analogous definition of the Wasserstein distance in Eq. (27) for two states
$\rho,\sigma$:

$\displaystyle{W}_{1,\square}(\rho,,\sigma):=\sup_{O=O^{\dagger},\|O\|_{\operatorname{Lip},\square}\leq
1}\operatorname{tr}\left[O(\rho-\sigma)\right]\,.$ (29)

The authors also put forth the following equivalent expression for the norm:

$\displaystyle{W}_{1,\square}(\rho,\sigma)=\frac{\min\left\\{\sum_{i=1}^{n}\|X^{(i)}\|_{1}:\rho-\sigma=\sum_{i=1}^{n}X^{(i)},X^{(i)}\in\mathcal{M}_{d^{n}}^{\operatorname{sa}},\operatorname{tr}_{i}[X^{(i)}]=0\right\\}}{2\sqrt{n}}\,.$
(30)

It follows from an application of Hölder’s inequality combined with the
variational formulation of Eq. (28) that
$\|O\|_{\operatorname{Lip},\square}\leq 2\sqrt{n}\|O\|_{\infty}$. However, it
can be the case that $\|O\|_{\operatorname{Lip}}\ll\sqrt{n}\,\|O\|_{\infty}$,
which are exactly those observables that should be thought of as regular. This
is because this signals that changing the state locally leads to significantly
smaller changes to the expectation value of the observable than global ones.
Two examples of this behaviour are given by the observables:

$\displaystyle
O_{1}=\sum\limits_{i=1}^{n}I_{i}\bigotimes\limits_{j\not=i}Z_{j},\quad~{}~{}~{}~{}~{}\text{
and }~{}~{}~{}~{}~{}\quad O_{2}=\sum\limits_{i=1}^{n}Z_{i}\,.$

Clearly, $\|O_{1}\|_{\infty}=\|O_{2}\|_{\infty}=n$. On the other hand, by
considering the states $\rho=|0\rangle\langle 0|^{\otimes n}$ and
$\sigma=|1\rangle\langle 1|\otimes|0\rangle\langle 0|^{\otimes n-1}$, we see
that $\|O_{1}\|_{\operatorname{Lip},\square}\geq\sqrt{n}\,(2n-2)$ whereas
$\|O_{2}\|_{\operatorname{Lip},\square}=2\sqrt{n}$. More generally, it is not
difficult to see that if $O=\sum\limits_{i=1}^{n}O_{i}$ with $\|O_{i}\|$ and
we denote by $\operatorname{supp}(O_{i})$ the qudits on which $O_{i}$ acts
nontrivially, then:

$\displaystyle\|O\|_{\operatorname{Lip},\square}\leq 2\sqrt{n}\max_{1\leq
j\leq n}\sum_{i}|\operatorname{supp}(O_{i})\cap\\{j\\}|.$

That is, the maximal number of intersections of the support of $O_{i}$ on one
qubit. From these examples we see that for local observables, the ratio of the
Wasserstein and Lipschitz constant reflects the locality of the observable.

####  B.1.2 Quantum differential Wasserstein distance

The Wasserstein distance $W_{1,\square}$ generalizes the classical Orstein
distance, that is the Wasserstein distance corresponding to the Hamming
distance on bit strings. Another definition of a Lipschitz constant and its
attached Wasserstein distance was put forth in [14], where the construction is
based on a differential structure that bears more resemblance to that of the
Lipschitz constant of a differentiable function on a continuous sample space,
e.g. a smooth manifold [53]. Let us now define the notion of a noncommutative
differential structure (see [13]):

######  Definition B.1 (Differential structure).

A set of operators $L_{k}\in\mathcal{M}_{d^{n}}$ and constants
$\omega_{k}\in\mathbbm{R}$ defines a differential structure
$\\{L_{k},\omega_{k}\\}_{k\in\mathcal{K}}$ for a full rank state
$\sigma\in\mathcal{D}_{d^{n}}$ if

  * 1

$\\{L_{k}\\}_{k\in\mathcal{K}}=\\{L_{k}^{\dagger}\\}_{k\in\mathcal{K}}$;

  * 2

$\\{L_{k}\\}_{k\in\mathcal{K}}$ consists of eigenvectors of the modular
operator $\Delta_{\sigma}(X):=\sigma X\sigma^{-1}$ with

$\displaystyle\Delta_{\sigma}(L_{k})=e^{-\omega_{k}}L_{k}\,.$ (31)

  * 3

$\|L_{k}\|_{\infty}\leq 1$.

Such a differential structure can be used to provide the set of matrices with
a Lipschitz constant that is tailored to $\sigma$, see e.g. [14, 13] for more
on this. In order to distinguish that constant from the one defined in (28),
we will refer to it as the differential Lipschitz constant and denote it by
$\|X\|_{\operatorname{Lip},\nabla}$. It is given by:

$\displaystyle\|X\|_{\operatorname{Lip},\nabla}:=\left(\sum_{k\in\mathcal{K}}(e^{-\omega_{k}/2}+e^{\omega_{k}/2})\|[L_{k},X]\|_{\infty}^{2}\right)^{1/2}\,.$
(32)

The quantity $[L_{k},X]$ should be interpreted as a partial derivative and is
sometimes denoted by $\partial_{k}X$ for that reason. Then, the gradient of a
matrix $A$, denoted by $\nabla A$ with a slight abuse of notations, refers to
the vector of operator valued coordinates $(\nabla A)_{i}=\partial_{i}A$. For
ease of notations, we will denote the differential structure by the couple
$(\nabla,\sigma)$. The notion of a differential structure is also intimately
connected to that of the generator of a quantum dynamical semigroup converging
to $\sigma$ [13], and properties of that semigroup immediately translate to
properties of the metric. This is because the differential structure can be
used to define an operator that behaves in an analogous way to the Laplacian
on a smooth manifold, which in turn induces a heat semigroup. We refer to [13,
14] for more details.

To illustrate the differential structure version of the Lipschitz constant, it
is instructive to think of the maximally mixed state. In this case, one
possible choice would consist of picking the $L_{k}$ to be all $1-$local Pauli
strings and $\omega_{j}=0$. Then the Lipschitz constant turns out to be given
by:

$\displaystyle\|X\|_{\operatorname{Lip},\nabla}=\left(\sum_{k\in\mathcal{K}}\|P_{k}XP_{k}-X\|_{\infty}^{2}\right)^{1/2},$
(33)

where $P_{k}$ are all $1-$local Pauli matrices. Thus, we see that this
measures how much the operator changes if we act locally with a Pauli unitary
on it. If we think of an operator as a function and conjugating with a Pauli
as moving in a direction, the formula above indeed looks like a derivative. In
fact, it is possible to make this connection rigorous, see [13].

As before, the definition in Eq. (32) yields a metric on states by duality:

$\displaystyle
W_{1,\nabla}(\rho,\sigma):=\sup\limits_{X=X^{\dagger},\,\|X\|_{\operatorname{Lip},\nabla}\leq
1}\left|\operatorname{Tr}\left(X(\rho-\sigma)\right)\right|.$

It immediately follows from the definitions that for any observable $X$:

$\displaystyle\left|\operatorname{tr}\left[X(\rho-\sigma)\right]\right|\leq\|X\|_{\operatorname{Lip},\nabla}\,W_{1,\nabla}(\rho,\sigma)\,.$
(34)

Although this geomemtric connection also opens up other interesting
mathematical connections for this definition, it has the downside of being
state dependent.

####  B.1.3 Wasserstein-like constant generated by quasi-local observables

As we will see in more detail in Section B.2, quasi-local observables are
Lipschitz with respect to $\|.\|_{\operatorname{Lip},\nabla}$ as well as
$\|.\|_{\operatorname{Lip},\square}$. This raises the question of defining the
notion of a Wasserstein-like quantity for this specific class of observables.
The reason for this will become clear in Section B.3 where we will be able to
easily leverage some recent Gaussian concentration bounds for Gibbs states at
high enough temperature [23] in order to prove transportation cost
inequalities for noncommutative Gibbs states, even in the case of long-range
interactions. Here, we consider the set $V:=\\{1\ldots n\\}$ of spins of local
dimension $d$. Then, we define a class of quasi-local observables:

######  Definition B.2 ($(k,g)$ quasi-local observables).

Given $k\in\mathbb{N}$ and $g\geq 0$, an observable $X\in\mathcal{M}_{d^{n}}$
is said to be $(k,g)$ quasi-local if it can be decomposed as

$\displaystyle X=\sum_{|A|\leq k}\,X_{A}~{}~{}~{}~{}\text{ with
}~{}~{}~{}~{}\sum_{A\ni v}\|X_{A}\|_{\infty}\leq g~{}~{}~{}\text{ for
}~{}~{}~{}v\in V\,,$

where each $X_{A}$ is strictly supported on a set $A$ of cardinality $|A|\leq
k$. Notice that this decomposition is unique. We denote by
$\mathcal{M}_{\operatorname{loc}}^{(k,g)}$ the space of $(k,g)$ quasi-local
observables. Moreover, a Gibbs state $\sigma$ on $V$ is called $(k,g)$ quasi-
local if it is the Gibbs state of a $(k,g)$ quasi-local Hamiltonian:

$\displaystyle\sigma:=\frac{e^{-\beta H}}{\mathcal{Z}}\,,$

where $\beta\geq 0$, $\mathcal{Z}:=\operatorname{tr}\left[e^{-\beta H}\right]$
for $H$ $(k,g)$ quasi-local.

Next, we introduce the concept of a quasi-local Wasserstein-like pseudometric:

######  Definition B.3 (Quasi-local Wasserstein pseudometric).

The quasi-local Wasserstein pseudometric on quantum states in
$\mathcal{D}_{d^{n}}$ is defined as:

$\displaystyle
W_{1,\operatorname{loc}}(\rho,\sigma):=\frac{1}{\sqrt{n}}\,\sup_{X\in\mathcal{M}_{\operatorname{loc}}^{(k,g)}}\left|\operatorname{tr}\left[X(\rho-\sigma)\right]\right|\,.$

Note that this is indeed a pseudo-metric, as it cannot distinguish states such
that all $k$-body marginals coincide. Thus, this last version of the
Wasserstein distance is the weaker of the three discussed, as it only
optimizes over a smaller set of observables.

###  B.2 Local evolution of Lipschitz observables

As already discussed in Subsections B.1.1 and B.1.2 when we defined
$\|.\|_{\operatorname{Lip},\nabla}$ and $\|.\|_{\operatorname{Lip},\square}$,
Lipschitz constants can be easily controlled by making assumptions on the
locality of the operators. Thus, if we apply a local circuit to a local
observable, it is straightforward to control the growth of the Lipschitz
constant in terms of the growth of the support of the observable under the
circuit. Indeed, in [16, Proposition 13] the authors show such a statement for
discrete time evolutions with exact lightcones: if we denote by $|L|$ the size
of the largest lightcone of one qubit under a channel $\Phi$, then for any
observable $O\in\mathcal{D}_{d^{n}}$,
$\|\Phi^{*}(O)\|_{\operatorname{Lip},\square}\leq|L|\|O\|_{\operatorname{Lip},\square}$.
Here we will extend such arguments to the evolution under local Hamiltonians
or Lindbladians. By resorting to Lieb-Robinson bounds, we show that the
Lipschitz constants $\|.\|_{\operatorname{Lip},\nabla}$ and
$\|.\|_{\operatorname{Lip},\square}$ of initially local observables evolving
according to a quasi-local dynamics increase at most with the effective
lightcone of the evolution. Thus, short-time dynamics or shallow-depth quantum
channels can only mildly increase the Lipschitz constant. This further
justifies the intuition that observables with small Lipschitz constant reflect
physical observables.

Lieb-Robinson (LR) bounds in essence assert that the time evolution of local
observables under (quasi)-local short-time dynamics have an effective
lightcone. There are many different formulations of Lieb-Robinson bounds.
Reviewing Lieb-Robinson bounds in detail is beyond the scope of this work and
we refer to [54, 55, 56, 57] and references therein for a review on them.
Lieb-Robinson bounds come in various versions. For studying the behaviour of
$\|.\|_{\operatorname{Lip},\nabla}$ under local evolutions, the most natural
one is the _commutator version_. A generator $\mathcal{L}$ of quasi-local
dynamics on $n$ qudits arranged on a graph $G=(V,E)$ is said to satisfy a LR
bound with velocity $v$ if for an observable $O_{A}$ supported on a region $A$
and another observable $B$ supported on a region $B$, we have:

$\displaystyle\|[\Phi_{t}(O_{A}),O_{B}]\|_{\infty}\leq
c\,(e^{vt}-1)\,g(d(A,B))\,\|O_{A}\|_{\infty}\|O_{B}\|_{\infty}\,,$ (LR1)

for $g:\mathbb{N}\to\mathbb{R}_{+}$ a function such that
$\lim_{x\to\infty}g(x)=0$. We again refer to [54, 55, 56, 57] for more
detailed statements depending on the structure of the generator $\mathcal{L}$.
We then have:

######  Proposition B.1 (Growth of differential Lipschitz constant for local
evolutions).

Let $\\{L_{k},\omega_{k}\\}_{k\in\mathcal{K}}$ be a differential structure on
$\mathcal{M}_{d^{n}}$ and let $O=\sum\limits_{i}O_{i}$ be an observable with
$\|O_{i}\|_{\infty}\leq 1$. Let $A_{i}$ denote the support of each $O_{i}$ and
$B_{j}$ of each $L_{j}$. Moreover, let $\Phi_{t}$ be an evolution satisfying
Eq. (LR1) and set $o(i,j)(t)=2$ if $A_{i}\cap B_{j}\not=\varnothing$ and
$c\,(e^{\mu vt}-1)\,g(d(A_{i},B_{j}))$ else. Then:

$\displaystyle\|\Phi_{t}(O)\|_{\operatorname{Lip},\nabla}^{2}\leq\sum_{k\in\mathcal{K}}(e^{\omega_{j}}+e^{-\omega_{j}})\left(\sum_{i}o_{i,j}(t)\right)^{2}.$

###### Proof.

The proof follows almost by definition. We have:

$\displaystyle\|\Phi_{t}(O)\|_{\operatorname{Lip},\nabla}^{2}=\sum_{j\in\mathcal{J}}(e^{\omega_{j}}+e^{-\omega_{j}})\|[\Phi_{t}(O),L_{j}]\|_{\infty}^{2}.$

By a triangle inequality we have that:

$\displaystyle\|[\Phi_{t}(O),L_{i}]\|_{\infty}^{2}\leq\left(\sum_{i}\|[\Phi_{t}(O_{i}),L_{j}]\|_{\infty}\right)^{2}$

For any term in the sum we have $\|[\Phi_{t}(O_{i}),L_{j}]\|\leq 2$ by the
submultiplicativity of the operator norm, a triangle inequality and
$\|L_{j}\|_{\infty}$. In case $A_{i}$ and $B_{j}$ do not overlap, the stronger
bound in Eq. (LR1) holds. ∎

To illustrate this bound more concretely, let us take $O=\sum_{i=1}^{n}Z_{i}$,
$L_{j}$ acting on $[j,j+k]$ for $j=1,\ldots,n-k$, and
$g(d(i,j))=e^{-\mu|i-j|}$, for some constant $\mu$. I.e. we have a $1D$
differential structure and a local time evolution on a $1D$ lattice. Then for
any $j$:

$\displaystyle\sum_{i}o_{i,j}(t)=k+(e^{vt}-1)\left(\sum_{i=1}^{j-1}e^{-\mu|i-j|}+\sum_{i=j+k+1}^{n}e^{-\mu|i-j|}\right)\leq
k+\frac{(e^{vt}-1)e^{-\mu}}{1-e^{-\mu}}.$

Thus,

$\displaystyle\|\Phi_{t}(O)\|_{\operatorname{Lip},\nabla}\leq\sqrt{n-k}\left(k+\frac{(e^{vt}-1)e^{-\mu}}{1-e^{-\mu}}\right).$

We see that for constant times the Lipschitz constant is still of order
$\sqrt{n}k$.

Let us now derive a similar, yet stronger, version of Prop. B.1 for
$\|\Phi_{t}(O)\|_{\operatorname{Lip},\square}$. In some situations, bounds
like (LR1) can be further exploited in order to prove the quasi-locality of
Markovian quantum dynamics [56]: for any region $C\subset D\subset V$, there
exists a local Markovian evolution $\Phi_{t}^{(D)}$ that acts non-trivially
only on region $D$, and such that for any observable $O_{C}$ supported on
region $C$,

$\displaystyle\|\left(\Phi^{*}_{t}-\left(\Phi_{t}^{(D)}\right)^{*}\right)(O_{C})\|_{\infty}\leq
c^{\prime}\,(e^{\mu vt}-1)\,h(d(C,V\backslash D))\|O_{C}\|_{\infty}\,,$ (LR2)

for some other function $h:\mathbb{N}\to\mathbb{R}_{+}$ such that
$\lim_{x\to\infty}h(x)=0$ and constant $c^{\prime}>0$. In other words, at
small times, the channels $\Phi_{t}^{*}$ can be well-approximated by local
Markovian dynamics when acting on local observables. Let us now estimate the
growth of Lipschitz constants for the definition of [16] given a Lieb-Robinson
bond:

######  Proposition B.2.

Assume that $\Phi_{t}$ satisfies the bound (LR2). Then, for any two quantum
states $\rho,\sigma\in\mathcal{D}_{d^{n}}$ and any ordering $\\{1,\cdots,n\\}$
of the graph:

$\displaystyle
W_{1,\square}(\Phi_{t}(\rho),\Phi_{t}(\sigma))\leq(2+2\,c^{\prime}\,(e^{\mu
vt}-1)\,\sum_{i=3}^{n}h(d(\\{i\cdots
n\\},\\{1\\})))\,W_{1,\square}(\rho,\sigma)\,.$ (35)

Moreover, for any observable $H\in\mathcal{M}_{d^{n}}$,

$\displaystyle\|\Phi_{t}(H)\|_{\operatorname{Lip},\square}\leq\left(2+2\,c^{\prime}\,(e^{\mu
vt}-1)\,\sum_{i=3}^{n}h(d(\\{i\cdots
n\\},\\{1\\}))\,\right)\|H\|_{\operatorname{Lip},\square}\,.$ (36)

###### Proof.

From [16], the Wasserstein distance $W_{1,\square}$ arises from a norm
$\|.\|_{W_{1}}$, i.e. $W_{1,\square}(\rho,\sigma)=\|\rho-\sigma\|_{W_{1}}$.
Moreover, the norm $\|.\|_{W_{1}}$ is uniquely determined by its unit ball
$\mathcal{B}_{n}$, which in turn is the convex hull of the set of the
differences between couples of neighboring quantum states:

$\displaystyle\mathcal{N}_{n}=\bigcup_{i\in
V}\,\mathcal{N}_{n}^{(i)}\,,~{}~{}~{}~{}~{}\mathcal{N}_{n}^{(i)}=\\{\rho-\sigma:\,\rho,\sigma\in\mathcal{D}_{d^{n}},\,\operatorname{tr}_{i}(\rho)=\operatorname{tr}_{i}(\sigma)\\}\,.$

Now by convexity, the contraction coefficient of this norm is equal to

$\displaystyle\|\Phi_{t}\|_{W_{1}\to
W_{1}}=\max\big{(}\|\Phi_{t}(X)\|_{W_{1}}:\,X\in\mathcal{M}_{d^{n}}^{\operatorname{sa},0},\,\|X\|_{W_{1}}\leq
1\big{)}=\max_{X\in\mathcal{N}_{n}}\|\Phi_{t}(X)\|_{W_{1}}\,,$

where $\mathcal{M}_{d^{n}}^{\operatorname{sa},0}$ denotes the set of self-
adjoint, traceless observables. Let then $X\in\mathcal{N}_{n}$. By the
expression (30), and choosing without loss of generality an ordering of the
vertices such that $\operatorname{tr}_{1}(X)=0$, we have

$\displaystyle\|\Phi_{t}(X)\|_{W_{1}}$
$\displaystyle\leq\frac{1}{2\sqrt{n}}\sum_{i=1}^{n}\Big{\|}\frac{I}{d^{i-1}}\otimes\operatorname{tr}_{1\cdots
i-1}\circ\Phi_{t}(X)-\frac{I}{d^{i}}\otimes\operatorname{tr}_{1\cdots
i}\circ\Phi_{t}(X)\Big{\|}_{1}$
$\displaystyle=\frac{1}{2\sqrt{n}}\,\sum_{i=1}^{n}\,\Big{\|}\,\int\,d\mu(U_{i})\,\operatorname{tr}_{1\cdots
i-1}\circ(\Phi_{t}(X)-U_{i}\Phi_{t}(X)U_{i}^{\dagger})\Big{\|}_{1}$
$\displaystyle\leq\frac{1}{2\sqrt{n}}\sum_{i=1}^{n}\int\,d\mu(U_{i})\,\|[U_{i},\operatorname{tr}_{1\cdots
i-1}\circ\Phi_{t}(X)]\|_{1}$
$\displaystyle\leq\frac{1}{\sqrt{n}}\,\sum_{i=1}^{n}\|\operatorname{tr}_{1\cdots
i-1}\circ\Phi_{t}(X)\|_{1}$
$\displaystyle\overset{(1)}{=}\frac{1}{\sqrt{n}}\,\sum_{i=1}^{n}\|\operatorname{tr}_{1\cdots
i-1}\circ(\Phi_{t}-\Phi^{(i-k\cdots n)}_{t})(X)\|_{1}$ (37)

where $\mu$ denotes the Haar measure on one qudit, and where $(1)$ follows
from the fact that $\operatorname{tr}_{1}(X)=0$, with $\Phi_{t}^{(i-k\cdots
n)}\equiv\Phi_{t}^{(\\{i-k,\cdots,n\\})}$ defined as in Eq. (LR2) with
$k<i-1$. Next, by the variational formulation of the trace distance and Eq.
(LR2), we have for $k\geq 3$ that

$\displaystyle\|\operatorname{tr}_{1\cdots i-1}\circ(\Phi_{t}-\Phi^{(i-k\cdots
n)}_{t})(X)\|_{1}$ $\displaystyle=\max_{\|O_{i\cdots n}\|_{\infty}\leq
1}\,\big{|}\operatorname{tr}\big{[}X(\Phi_{t}-\Phi^{(i-k\cdots
n)}_{t})(O_{i\cdots n})\big{]}\big{|}$ $\displaystyle\leq\,\max_{\|O_{i\cdots
n}\|_{\infty}\leq 1}\,\|(\Phi_{t}-\Phi^{(i-k\cdots n)}_{t})(O_{i\cdots
n})\|_{\infty}\|X\|_{1}$ $\displaystyle\leq\,c^{\prime}\,(e^{\mu
vt}-1)\,h(d(\\{i\cdots n\\},\\{1\cdots i-k-1\\}))\,\|X\|_{1}$
$\displaystyle\overset{(2)}{\leq}2\,c^{\prime}\,(e^{\mu vt}-1)\,h(d(\\{i\cdots
n\\},\\{1\cdots i-k-1\\}))\,\sqrt{n}\,\|X\|_{W_{1}}\,,$

where $(2)$ follows from [16, Proposition 6]. By picking $k=i-2$ and inserting
this estimate into Eq. (B.2) for $i\geq 3$ and the trivial estimate
$\|\operatorname{tr}_{1\cdots i-1}\circ(\Phi_{t}-\Phi^{(i-k\cdots
n)}_{t})(X)\|_{1}\leq 2\|X\|_{1}$ for $i=1,2$ we obtain Eq. (35). Eq. (36)
follows by duality. ∎

###  B.3 Transportation cost inequalities

Although interesting on their own, the relevance of the quantities introduced
above becomes clearer in our context when we also have a _transportation cost_
inequality [58, 53]. A quantum state $\sigma$ satisfies a transportation cost
inequality with constant $\alpha$ if for any other state $\rho$ it holds that:

$\displaystyle
W_{1}(\rho,\sigma)\leq\sqrt{\frac{1}{2\alpha}\,D(\rho\|\sigma)}\,,$ (38)

where $W_{1}\in\\{{W}_{1,\square},W_{1,\nabla},W_{1,\operatorname{loc}}\\}\,.$
In what follows, we simply write $\|.\|_{\operatorname{Lip}}$ and $W_{1}$ to
denote either of the Lipschitz constants, and their corresponding Wasserstein
distances, defined above. This inequality should be thought of as a stronger
version of Pinsker’s inequality that is tailored to a state $\sigma$ and the
underlying Wasserstein distance.

One of the well-established techniques to establish a transportation cost
inequality for $W_{1,\nabla}$ is by exploiting the fact that it is implied by
a logarithmic Sobolev inequality, another functional inequality. It is beyond
the scope of this paper to explain this connection and we refer to e.g. [14]
and references therein for a discussion on these topics. But for our purposes
it is important to note that in [25] the authors and Capél show modified
logarithmic Sobolev inequalities for several classes of Gibbs states. This
connection with other functional inequalities is one of the main advantages of
$W_{1,\nabla}$ when compared with $W_{1,\square}$, as it endows the metric
with a rich structure. In the Theorem below we summarize the regimes for which
the results of [25] give us a transportation cost inequality:

######  Theorem B.1 (transportation cost for commuting Hamiltonians [25, 18]).

Let $E_{1},\ldots,E_{m}\subset\mathcal{M}_{d^{n}}$ be a set of $k-local$
linearly independent _commuting_ operators with $\|E_{i}\|_{\infty}\leq 1$.
Then $\sigma(\lambda)$ satisfies a transportation cost inequality for all
$\lambda\in B_{\ell_{\infty}}(0,1)$ and $W_{1,\nabla},W_{1,\square}$ with
$\alpha=\Omega(1)$ in case:

  * $\operatorname{(i)}$

The $E_{i}$ are classical and all $\sigma(\lambda)$ satisfy exponential decay
of correlations [25].

  * $\operatorname{(ii)}$

The operators $E_{i}$ are local with respect to a hypergraph $G=(V,E)$ and the
inverse temperature satisfies $\beta<\beta_{c}$, where $\beta_{c}$ only
depends on $k$ and properties of the interaction hypergraph [18, Theorem 3,
Proposition 4].

  * $\operatorname{(iii)}$

The $E_{i}$ are one-dimensional.

Moreover, the underlying differential structure $(\nabla,\sigma)$ consists of
$L_{k}$ acting on at most $\mathcal{O}(k)$ qudits.

Theorem B.1 establishes that transportation cost inequalities are satisfied
for most classes of commuting Hamiltanians known to have exponential decay of
correlations.

In [18, Proposition 5], the authors show that $W_{1,\square}\leq
c(k)W_{1,\nabla}$ holds up to some constant $c(k)$ depending on the locality
of the differential structure. This implies that a transportation cost
inequality for $W_{1,\nabla}$ implies one for $W_{1,\square}$ with the same
constant. Thus, although the authors of [25] only obtain the result for
$W_{1,\nabla}$, we can use it to translate it to $W_{1,\square}$. We conclude
that for commuting Hamiltonians TC are available for essentially all classes
of local Hamiltonians for which they are expected to hold.

In order to consider a larger class of states that satisfy a transportation
cost inequality, it is fruitful to restrict the analysis to quasi-local
observables and consider $W_{1,\operatorname{loc}}$. In this case, the
following result holds:

######  Theorem B.2 (Transportation cost for quasi-local Wasserstein
distances).

For $(k,g)\in\mathbb{N}\times\mathbb{R}_{+}$, we let $\sigma$ be a $(k,g)$
quasi-local Gibbs state at inverse temperature
$\beta<\beta_{c}:=\frac{1}{8e^{3}\,g\,k}$. Then, the following holds for all
$\rho\in\mathcal{D}_{d^{n}}$:

$\displaystyle
W_{1,\operatorname{loc}}(\rho,\sigma)\leq\sqrt{\frac{2g}{\beta_{c}-\beta}\,D(\rho\|\sigma)}\,.$

###### Proof.

We resort to [23, Theorem 1], where it is showed that for $\beta<\beta_{c}$,
any $X\in\mathcal{M}_{\operatorname{loc}}^{(k,g)}$ and
$|\tau|\leq\beta_{c}-\beta$ such that
$\operatorname{tr}\left[X\sigma\right]=0$,

$\displaystyle\ln\big{[}\operatorname{tr}\left[e^{-\tau
X}\sigma\right]\big{]}\leq\frac{\tau^{2}\sum_{|A|\leq
k}\|X_{A}\|_{\infty}}{\beta_{c}-\beta-|\tau|}\leq\frac{\tau^{2}gn}{\beta_{c}-\beta-|\tau|}\,.$

Next, by resorting to Golden-Thompson’s inequality as well as the variational
formulation of the relative entropy, we have

$\displaystyle-\tau\operatorname{tr}\left[\rho
X\right]-D(\rho\|\sigma)\leq\ln\operatorname{tr}\left[e^{-\tau
X+\ln\sigma}\right]\leq\ln\big{[}\operatorname{tr}\left[e^{-\tau
X}\sigma\right]\big{]}\leq\frac{\tau^{2}gn}{\beta_{c}-\beta-|\tau|}\,.$

Next, we can choose $\tau=-\frac{\alpha\operatorname{tr}\left[\rho
X\right]}{gn}$ for some parameter $\alpha$ to be chosen later, so that

$\displaystyle\frac{\alpha\,(\beta_{c}-\beta)}{\beta_{c}-\beta+\alpha}\,\frac{\operatorname{tr}\left[\rho
X\right]^{2}}{gn}\leq D(\rho\|\sigma)\,.$

The result then follows after taking $\alpha=\beta_{c}-\beta$, so that
$\tau=\frac{\beta_{c}-\beta}{2}$, then optimizing over observables
$X\in\mathcal{M}_{\operatorname{loc}}^{(k,g)}$. ∎

Theorem B.2 has the advantage that it holds for esentially all high-
temperature Gibbs states, even with long-range interactions. However, its
weakness is also apparent: it only allows us to obtain recovery guarantees for
observables with support on at most $k$ qubits, while the other transportation
inequalities allowed for arbitrarily large supports. By exploiting the trade-
off between the locality and the temperature, it is possible to extend the
result for observables with higher locality by suitably rescaling the inverse
temperature range.

Finally, let us be a bit more explicit on how to leverage our results to also
cover the outputs of short depth circuits. To this end, let $G=(V,E)$ be a
graph which models the interactions of a unitary circuit and suppose we
implement $L$ layers of a unknown unitary circuit consisting of $1$ and
$2$-qubit unitaries layed out according to $G$. That is, we have an unknown
shallow circuit $\mathcal{U}$ of depth $L$ with respect to $G$. More
precisely,

$\displaystyle\mathcal{U}=\prod_{\ell\in[L]}\,\bigotimes_{e\in
E_{\ell}}\,\mathcal{U}_{\ell,e},$

where each $E_{\ell}\subset E$ are a subset of the edges.

Our goal will be to show how to approximately learn the state

$\displaystyle\ket{\psi}=\mathcal{U}\ket{0}^{\otimes n}.$

We will first find an appropriate Gibbs state to approximate it in Wasserstein
distance.

First, note that for $\beta_{\epsilon}=\log(\epsilon^{-1})$ and
$H_{0}=-\sum_{i}^{n}Z_{i}$ we have that:

$\displaystyle D\left(\mathcal{U}|0\rangle\langle 0|^{\otimes
n}\mathcal{U}^{\dagger}\|\mathcal{U}\frac{e^{-\beta_{\epsilon}H_{0}}}{\operatorname{tr}\left[e^{-\beta_{\epsilon}H_{0}}\right]}\mathcal{U}^{\dagger}\right)=D\left(|0\rangle\langle
0|^{\otimes
n}\|\frac{e^{-\beta_{\epsilon}H_{0}}}{\operatorname{tr}\left[e^{-\beta_{\epsilon}H_{0}}\right]}\right)=nD\left(|0\rangle\langle
0|\|\frac{e^{-\beta_{\epsilon}Z}}{\operatorname{tr}\left[e^{-\beta_{\epsilon}Z}\right]}\right)=$
$\displaystyle ne^{\beta_{\epsilon}}(1-\log(1-e^{-2\beta_{\epsilon}}))\leq
ne^{-\beta_{\epsilon}}=\epsilon n.$

Thus, if the states
$\mathcal{U}\frac{e^{-\beta_{\epsilon}H_{0}}}{\operatorname{tr}\left[e^{-\beta_{\epsilon}H_{0}}\right]}\mathcal{U}^{\dagger}$
satisfy a transportation cost inequality with some constant $\alpha$, this
would allow us to conclude that

$\displaystyle W_{1}\left(\mathcal{U}|0\rangle\langle 0|^{\otimes
n}\mathcal{U}^{\dagger},\mathcal{U}\frac{e^{-\beta_{\epsilon}H_{0}}}{\operatorname{tr}\left[e^{-\beta_{\epsilon}H_{0}}\right]}\mathcal{U}^{\dagger}\right)\leq\sqrt{\frac{\epsilon
n}{2\alpha}}.$

Moreover, defining
$H_{\mathcal{U}}=-\sum_{i}^{n}\mathcal{U}Z_{i}\mathcal{U}^{\dagger}$, we see
that
$\mathcal{U}\frac{e^{-\beta_{\epsilon}H_{0}}}{\operatorname{tr}\left[e^{-\beta_{\epsilon}H_{0}}\right]}\mathcal{U}^{\dagger}=\frac{e^{-\beta_{\epsilon}H_{\mathcal{U}}}}{\operatorname{tr}\left[e^{-\beta_{\epsilon}H_{0}}\right]}$.
As we know the geometry of $\mathcal{U}$, we can bound the support of
$\mathcal{U}Z_{i}\mathcal{U}^{\dagger}$. Thus, we only need to find a suitable
transportation cost inequality to see that this approximation fits into our
framework.

Let us now find a suitable differential structure for the state
$\frac{e^{-\beta_{\epsilon}H_{0}}}{\operatorname{tr}\left[e^{-\beta_{\epsilon}H_{0}}\right]}$.
For simplicity, denote by $\tau_{p}=p|0\rangle\langle 0|+(1-p)|1\rangle\langle
1|$ and note that
$\frac{e^{-\beta_{\epsilon}H_{0}}}{\operatorname{tr}\left[e^{-\beta_{\epsilon}H_{0}}\right]}=\tau_{p_{\epsilon}}^{\otimes
n}$ with
$p_{\epsilon}=\frac{e^{\beta_{\epsilon}}}{e^{\beta_{\epsilon}}+e^{-\beta_{\epsilon}}}$.
Moreover, let

$\displaystyle a_{i}=I^{\otimes i-1}\otimes\ket{0}\bra{1}\otimes I^{\otimes
n-i-1}$

be the anihilation operator acting on qubit $i$. Definining
$L_{i,0}=\left(p(1-p)\right)^{\frac{1}{4}}a_{i}$ and
$L_{i,1}=L_{i,0}^{\dagger}$, we can define a differential structure for
$\tau_{p}^{\otimes n}$ with $\\{L_{i,k},\omega_{i,k}\\}$ for $i=1,\ldots,n$
and $k=0,1$ and $\omega_{i,0}=\frac{p}{1-p}$ and $\omega_{i,1}=\frac{1-p}{p}$.
That this is indeed a differential structure follows from a simple
computation. One can readily check that the induced Lipschitz constant is
given by:

$\displaystyle\|O\|_{\operatorname{Lip},\nabla}^{2}=\sum_{i=1}^{n}\left(\sqrt{\frac{p}{1-p}}+\sqrt{\frac{1-p}{p}}\right)(\|[O,L_{i,0}]\|^{2}+\|[O,L_{i,1}]\|^{2})$
$\displaystyle=\sum_{i=1}^{n}\|[O,a_{i}]\|^{2}+\|[O,a_{i}^{\dagger}]\|^{2}.$

Thus, we see that the Lipschitz constant takes a particularly simple form for
this differential structure. Moreover, it is not difficult to see that
$\\{\mathcal{U}L_{i,k}\mathcal{U}^{\dagger},\omega_{i,k}\\}$ provides a
differential structure for the state $\mathcal{U}\tau_{p}^{\otimes
n}\mathcal{U}^{\dagger}$. Indeed, it is easily checked that this new
differential structure still gives eigenvectors of the modular operator.
Importantly, the result of [59, Theorem 19] establish that the state
$\tau_{p}^{\otimes n}$ satisfies a transportation cost inequality with
constant $\frac{1}{2}$. Putting all these elements together we have:

######  Theorem B.3 (transportation cost for shallow circuits).

Let $\mathcal{U}$ be an unknown depth $L$ quantum circuit on $n$ qubits
defined on a graph $G=(V,E)$ and $\ket{\psi}=\mathcal{U}\ket{0}^{\otimes n}$.
Define for $\epsilon>0$
$H_{\mathcal{U}}=-\sum_{i}\mathcal{U}Z_{i}\mathcal{U}^{\dagger}$ and
$\sigma(U)=\frac{e^{-\beta_{\epsilon}H_{\mathcal{U}}}}{\operatorname{tr}\left[e^{-\beta_{\epsilon}H_{\mathcal{U}}}\right]}$
with $\beta_{\epsilon}=\log(\epsilon^{-1})$. Then for any state $\rho$ and all
observables $O$ we have:

$\displaystyle|\operatorname{tr}\left[O(|\psi\rangle\langle\psi|-\rho)\right]|\leq\|O\|_{\operatorname{Lip},\nabla}\left(\sqrt{\epsilon
n}+\sqrt{D(\rho\|\sigma(\mathcal{U}))}\right)$

with

$\displaystyle\|O\|_{\operatorname{Lip},\nabla}^{2}=\sum_{i=1}^{n}\|[O,\mathcal{U}a_{i}\mathcal{U}^{\dagger}]\|^{2}+\|[O,\mathcal{U}a_{i}^{\dagger}\mathcal{U}^{\dagger}]\|^{2}.$

###### Proof.

We have:

$\displaystyle|\operatorname{tr}\left[O(|\psi\rangle\langle\psi|-\rho)\right]|\leq|\operatorname{tr}\left[O(|\psi\rangle\langle\psi|-\sigma(\mathcal{U}))\right]|+|\operatorname{tr}\left[O(\sigma(\mathcal{U})-\rho)\right]|.$

The claim then follows from the discussion above, as $\sigma(\mathcal{U})$
satisfies a transportation cost inequality with constant $\frac{1}{2}$. ∎

Of course the result above has the downside that the Lipschitz constant
depends on the unknown circuit $\mathcal{U}$. However, as we can estimate the
locality of each term $\mathcal{U}a_{i}\mathcal{U}^{\dagger}$, it is also
possible to estimate the Lipschitz constant by estimating the overlap of the
observable $O$ with each $\mathcal{U}a_{i}\mathcal{U}^{\dagger}$.

##  Appendix C Combining the maximum entropy method with transportation cost
inequalities

With these tools at hand, we are now ready to show that resorting to
transportation cost inequalities it is possible to obtain exponential
improvements in the sample complexity required to learn a Gibbs state. First,
let us briefly review shadow tomography or Pauli regrouping techniques [9, 10,
11, 12]. Although these methods all work under slightly different assumptions
and performance guarantees, they have in common that they allow us to learn
the expectation value of $M$ $k$-local observables
$O_{1},\ldots,O_{M}\in\mathcal{M}_{2^{n}}$ such that $\|O_{i}\|_{\infty}\leq
1$ up to an error $\epsilon$ and failure probability at most $1-\delta$ by
measuring $\mathcal{O}(e^{\mathcal{O}(k)}\log(M\delta^{-1})\epsilon^{-2})$
copies of the state. From now on we will set the probability of failure to
$\delta=\textrm{poly}(n)^{-1}$ and suppress it.

For instance, for the shadow methods of [9], we obtain a
$\mathcal{O}(4^{k}\log(M)\epsilon^{-2})$ scaling by measuring $1$-qubit random
bases. The estimate is then obtained by an appropriate median of means
procedure of the expectation value of each output string. Thus, the complexity
of the procedure is comparable to evaluating the expectation value of product
states on the $E_{i}$ $\mathcal{O}(4^{k}\log(M)\epsilon^{-2})$ times. For
$k$-local observables, this is translated to a postprocessing scaling like
$\mathcal{O}(e^{ck}\log(M)\epsilon^{-2})$ for some $c>0$.

The application of such results to maximum entropy methods is then clear :
given $E_{1},\ldots,E_{m}$ that are at most $k-$local, with probability at
least $1-\delta$ we can obtain an estimate $e^{\prime}(\lambda)$ of
$e(\lambda)$ satisfying:

$\displaystyle\|e^{\prime}(\lambda)-e(\lambda)\|_{\ell_{p}}\leq
m^{\frac{1}{p}}\epsilon$ (39)

using $\mathcal{O}(4^{k}\log(M)\epsilon^{-2}\log(\delta^{-1}))$ copies of
$\sigma(\lambda)$.

We then finally obtain:

######  Theorem C.1 (Fast learning of Gibbs states).

Let $\sigma(\lambda)\in\mathcal{D}_{2^{n}}$ be an $n$-qubit Gibbs state at
inverse temperature $\beta$ with respect to a set of $k$-local operators
$\mathcal{E}=\\{E_{i}\\}_{i=1}^{m}$ that satisfies a transportation cost
inequality with constant $\alpha$. Moreover, assume that
$\log(\mathcal{Z}(\mu))$ is $L,U$ strongly convex for $\mu\in
B_{\ell_{\infty}}(0,1)$. Then

$\displaystyle\mathcal{O}\left(4^{k}\alpha^{-1}\epsilon^{-2}\beta\log(m\delta^{-1})\min\left\\{L^{-1},\beta\epsilon^{-2}\right\\}\right)$

samples of $\sigma(\lambda)$ are sufficient to obtain a state $\sigma(\mu)$
that satisfies:

$\displaystyle|\operatorname{tr}\left[O(\sigma(\lambda)-\sigma(\mu)\right]|\leq\epsilon\,n^{\frac{1}{2}}\|O\|_{\operatorname{Lip}}$
(40)

for all Lipschitz observables $O$ with probability at least $1-\delta$.

###### Proof.

Using the aforementioned methods of [9] we can obtain an estimate
$e^{\prime}(\lambda)$ satisfying the guarantee of Eq. (39) with probability at
least $1-\delta$. Now we solve the maximum entropy problem with
$e^{\prime}(\lambda)$ and set the stopping criterium for the output $\mu_{*}$
as

$\displaystyle\|e^{\prime}(\lambda)-e(\mu_{*})\|_{\ell_{2}}<(4c+1)\epsilon\sqrt{m}.$

This implies that:

$\displaystyle\left|\beta\,\langle\lambda-\mu|e(\lambda)-e(\mu)\rangle\right|\leq\beta\|e(\lambda)-e(\mu_{*})\|_{\ell_{2}}\|\lambda-\mu_{*}\|_{\ell_{2}}=\mathcal{O}\left(\beta\min\\{L^{-1}\epsilon^{2}m,\epsilon
m\\}\right).$

The second term in the minimum follows from the estimate
$\|\lambda-\mu_{*}\|_{\ell_{2}}\leq 2\sqrt{m}$, while the first from the
strong convexity of the maximum entropy problem discussed in Eq. 20, which
yields:

$\displaystyle\|\lambda-\mu_{*}\|_{\ell_{2}}\leq
L^{-1}\|e(\lambda)-e(\mu_{*})\|_{\ell_{2}}.$

Then it follows from Prop. A.6 that we have:

$\displaystyle D(\sigma(\mu_{*})\|\sigma(\lambda))\leq
D(\sigma(\mu_{*})\|\sigma(\lambda))+D(\sigma(\lambda)\|\sigma(\mu_{*}))=\mathcal{O}\left(\beta\min\\{L^{-1}\epsilon^{2}m,\epsilon
m\\}\right).$ (41)

This can then be combined with transportation cost inequalities. Indeed, we
have:

$\displaystyle|\operatorname{tr}\left[O(\sigma(\lambda)-\sigma(\mu_{*})\right]|\leq\|O\|_{\operatorname{Lip}}W_{1}(\sigma(\mu_{*}),\sigma(\lambda))\leq\|O\|_{\operatorname{Lip}}\sqrt{\frac{D(\sigma(\mu_{*}),\sigma(\lambda))}{2\alpha}}.$

Inserting our bound on the relative entropy in Eq. (41) we obtain:

$\displaystyle|\operatorname{tr}\left[O(\sigma(\lambda)-\sigma(\mu_{*})\right]|=\mathcal{O}\left(\|O\|_{\operatorname{Lip}}\sqrt{\beta}\min\\{(\alpha
L)^{-\frac{1}{2}}\epsilon\sqrt{m},\sqrt{\epsilon m}\\}\right).$

We conclude by suitably rescaling $\epsilon$. ∎

The theorem above yields exponential improvements in the sample complexity to
learn the expectation value of certain observables for classes of states that
satisfy a transportation cost inequality with $\alpha=\Omega(\log(n)^{-1})$.
As discussed in Sec. B.2, extensive observables that can be written as a sum
of $l$-local observables have a constant norm that satisfies
$\|O\|_{\operatorname{Lip}}=\mathcal{O}(l\sqrt{n})$. Shadow-like methods would
require $\mathcal{O}(e^{\mathcal{O}(l)}\log(m)\epsilon^{-2}\log(\delta^{-1}))$
samples to learn such observables up to a relative error $n\epsilon$. Our
methods, however, require
$\mathcal{O}(e^{\mathcal{O}(k)}\operatorname{poly}(l,\epsilon^{-1})\log(m)\log(\delta^{-1}))$,
which yields exponential speedups in the regime
$l=\operatorname{poly}(\log(n),\epsilon^{-1})$. Of course it should also be
mentioned that classical shadows do not require any assumptions on the
underlying state.

Furthermore, considering the exponential dependency of the sample complexity
in the locality for shadow-like methods, we believe that our methods yield
practically significant savings already in the regime in which we wish to
obtain expectation values of observables with relatively small support. For
instance, for high-temperature Gibbs state of nearest neighbor Hamiltonians
and observables supported on $15$ qubits, shadows require a factor of $\sim
10^{7}$ more samples on observables supported on $15$ qubits than solving the
maximum entropy problem for obtaining the same precision.

On the other hand, previous guarantees on learning quantum many-body states
[3, 4, 37, 7] required a polynomial in system size precision to obtain a
nontrivial approximation, which implies a polynomial-time complexity. Thus,
our recovery guarantess are also exponentially better compared to standard
many-body tomography results.

##  Appendix D Regimes of efficient postprocessing

The only question we have still to answer is how to perform the postprocessing
efficiently. There have been many recent breakthroughs on approximating
quantum Gibbs states efficiently on classical computers [20, 21, 22, 23, 24].
The gradient descent method only requires us to estimate the gradient of the
partition function for a Gibbs state at each iteration. Thus, any algorithm
that can approximate the log-partition function $\mathcal{Z}(\mu)$ efficiently
or approximate $e(\lambda)$ suffices for our purposes.

For Gibbs states on a lattice, the methods of [21] yield that we can perform
such approximations in time polynomial in $n$ for temperatures
$\beta<\beta_{c}=1/(k8e^{3})$ , where $k$ is the locality of the Hamiltonian,
and $\epsilon$ inverse polynomial in system size. Thus, we conclude that for
this temperature range, which coincides with the range for which the results
of B.2 hold, we can also obtain efficient classical postprocessing.

For the case of $1D$ systems, to the best of our knowledge the best results
available are those of [24]. They show how to efficiently obtain efficient
tensor network approximations for $1D$ Gibbs states for $\beta=o(\log(n))$. As
such tensor networks can be contracted efficiently as well, this gives an
efficient classical algorithm to compute local expectation values of such
states, which suffices for our purposes. Thus, the results of [21, 24] ensure
that for the systems considered in Thm. B.1 we can also perform the
postprocessing efficiently on a classical computer.

It is also interesting to consider what happens if we have access to a quantum
computer to estimate the gradient. We will discuss the implications of this in
the next section for the case of commuiting Hamiltonians.

##  Appendix E A complete picture: commuting Gibbs states