# B Lipschitz constants and transportation cost inequalities

In this section, we identify conditions under which it is possible to estimate all expectation values of k-local observables up to an error $\epsilon$ by measuring $O(\text{poly}(k, \log(n), \epsilon^{-1}))$ copies of it, where $n$ is the system size, which constitutes an exponential improvement in some regimes. To obtain this result, we combine the maximum entropy method introduced in Section A 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 [61], constitute one of the strongest tools available to show concentration of measure inequalities. In the quantum setting, their study was initiated in [19, 20, 56, 59]. 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$:

$$ \|\rho - \sigma\|_{tr} = \sup_{P=P^{\dagger}, \|P\|_{\infty} \le 1} \text{tr}[P(\rho - \sigma)]. $$

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 extensive observables of the system. 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 [56, 59].

In the classical setting [58, 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 \mathbb{R}$ to be:

$$ \|f\|_{\text{Lip}} := \sup_{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

$$ W_1(p, q) := \sup_{f: \|f\|_{\text{Lip}} \le 1} |\mathbb{E}_p(f) - \mathbb{E}_q(f)|. \qquad (33) $$

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 [58, 65] 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 [26, 45, 56, 59]. 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 [56]:

$$ \|O\|_{\text{Lip},\square} = \sqrt{n} \max_{\substack{1 \le i \le n \\ \rho, \sigma \in \mathcal{D}_{d^n}}} \max_{\text{tr}_i[\rho]=\text{tr}_i[\sigma]} \text{tr}[O(\rho - \sigma)]. \qquad (34) $$
That is, the Lipschitz constant quantifies the amount by which 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. (33) for two states $\rho, \sigma$:

$$W_{1,\square}(\rho, \sigma) := \sup_{O=O^{\dagger}, \|O\|_{\text{Lip},\square} \le 1} \text{tr}[O(\rho - \sigma)]. \qquad (35)$$

The authors of [56] also put forth the following equivalent expression for the norm:

$$W_{1,\square}(\rho, \sigma) = \frac{\min\{\sum_{i=1}^n \|X^{(i)}\|_1 : \rho - \sigma = \sum_{i=1}^n X^{(i)}, X^{(i)} \in \mathcal{M}_{d^n}^{\text{sa}}, \text{tr}_i[X^{(i)}] = 0\}}{2\sqrt{n}}. \qquad (36)$$

It follows from an application of Hölder's inequality combined with the variational formulation in Eq. (34) that $\|O\|_{\text{Lip},\square} \le 2\sqrt{n}\|O\|_{\infty}$. However, it can be the case that $\|O\|_{\text{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 behavior are given by the observables:

$$O_1 = \sum_{i=1}^{n} Z_{i^c}, \quad \text{and} \quad O_2 = \sum_{i=1}^{n} Z_i,$$

where $Z_{i^c}$ acts as identity on $i$ and $Z_i$ else, i.e. $Z_{1^c} = I_1 \otimes Z_2 \otimes Z_3 \otimes \cdots \otimes Z_n$. 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 n-1}$, we see that $\|O_1\|_{\text{Lip},\square} \ge \sqrt{n}(2n-2)$ whereas $\|O_2\|_{\text{Lip},\square} = 2\sqrt{n}$. More generally, it is not difficult to see that if $O = \sum_{i=1}^n O_i$ with $\|O_i\|_{\infty} \le 1$ and we denote by $\text{supp}(O_i)$ the qudits on which $O_i$ acts nontrivially, then:

$$\|O\|_{\text{Lip},\square} \le 2\sqrt{n} \max_{1 \le j \le n} \sum_i |\text{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 operator norm 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 [59], 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 [58]. Let us now define the notion of a noncommutative differential structure (see [21]):

**Definition B.1 (Differential structure).** A set of operators $L_k \in \mathcal{M}_{d^n}$ and constants $\omega_k \in \mathbb{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}} \text{ consists of eigenvectors of the modular operator } \Delta_{\sigma}(X) := \sigma X \sigma^{-1} \text{ with}$$

$$\Delta_{\sigma}(L_k) = e^{-\omega_k} L_k. \qquad (37)$$

$$3 \|L_k\|_{\infty} \le 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. [21, 59] for more on this. In order to distinguish that constant from the one defined in (34), we will refer to it as the differential Lipschitz constant and denote it by $\|X\|_{\text{Lip},\nabla}$. It is given by:

$$\|X\|_{\text{Lip},\nabla} := \left( \sum_{k \in \mathcal{K}} (e^{-\omega_k/2} + e^{\omega_k/2}) \|[L_k, X]\|_{\infty}^2 \right)^{1/2}. \qquad (38)$$
That is, the Lipschitz constant quantifies the amount by which 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. (33) for two states $\rho, \sigma$:

$$W_{1,\square}(\rho, \sigma) := \sup_{O=O^{\dagger}, \|O\|_{\text{Lip},\square} \le 1} \text{tr}[O(\rho - \sigma)]. \quad (35)$$

The authors of [56] also put forth the following equivalent expression for the norm:

$$W_{1,\square}(\rho, \sigma) = \frac{\min\{\sum_{i=1}^n \|X^{(i)}\|_1 : \rho - \sigma = \sum_{i=1}^n X^{(i)}, X^{(i)} \in \mathcal{M}_{d^n}^{\text{sa}}, \text{tr}_i[X^{(i)}] = 0\}}{2\sqrt{n}}. \quad (36)$$

It follows from an application of Hölder's inequality combined with the variational formulation in Eq. (34) that $\|O\|_{\text{Lip},\square} \le 2\sqrt{n}\|O\|_{\infty}$. However, it can be the case that $\|O\|_{\text{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 behavior are given by the observables:

$$O_1 = \sum_{i=1}^{n} Z_{i^c}, \quad \text{and} \quad O_2 = \sum_{i=1}^{n} Z_i,$$

where $Z_{i^c}$ acts as identity on $i$ and $Z_i$ else, i.e. $Z_{1^c} = I_1 \otimes Z_2 \otimes Z_3 \otimes \cdots \otimes Z_n$. 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 n-1}$, we see that $\|O_1\|_{\text{Lip},\square} \ge \sqrt{n}(2n-2)$ whereas $\|O_2\|_{\text{Lip},\square} = 2\sqrt{n}$. More generally, it is not difficult to see that if $O = \sum_{i=1}^n O_i$ with $\|O_i\|_{\infty} \le 1$ and we denote by $\text{supp}(O_i)$ the qudits on which $O_i$ acts nontrivially, then:

$$\|O\|_{\text{Lip},\square} \le 2\sqrt{n} \max_{1 \le j \le n} \sum_i |\text{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 operator norm 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 [59], 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 [58]. Let us now define the notion of a noncommutative differential structure (see [21]):

**Definition B.1 (Differential structure).** A set of operators $L_k \in \mathcal{M}_{d^n}$ and constants $\omega_k \in \mathbb{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}} \text{ consists of eigenvectors of the modular operator } \Delta_{\sigma}(X) := \sigma X \sigma^{-1} \text{ with}$$

$$\Delta_{\sigma}(L_k) = e^{-\omega_k} L_k. \quad (37)$$

$$3 \|L_k\|_{\infty} \le 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. [21, 59] for more on this. In order to distinguish that constant from the one defined in (34), we will refer to it as the differential Lipschitz constant and denote it by $\|X\|_{\text{Lip},\nabla}$. It is given by:

$$\|X\|_{\text{Lip},\nabla} := \left( \sum_{k \in \mathcal{K}} (e^{-\omega_k/2} + e^{\omega_k/2}) \|[L_k, X]\|_{\infty}^2 \right)^{1/2}. \quad (38)$$
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 notation, 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$ [21], 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 the heat semigroup. We refer to [21, 59] 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:

$$ \|X\|_{\text{Lip},\nabla} = \left( \sum_{k \in \mathcal{K}} \|P_k X P_k - X\|_{\infty}^2 \right)^{1/2}, \qquad (39) $$

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 [21].

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

$$ W_{1,\nabla}(\rho, \sigma) := \sup_{X=X^{\dagger}, \|X\|_{\text{Lip},\nabla} \le 1} |\text{Tr}(X(\rho - \sigma))|. $$

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

$$ |\text{tr}[X(\rho - \sigma)]| \le \|X\|_{\text{Lip},\nabla} W_{1,\nabla}(\rho, \sigma). \qquad (40) $$

Although this geometric interpretation opens up other interesting mathematical connections for this definition, the differential Wasserstein distance has the downside of being state dependent. It however induces a stronger topology than the quantum Hamming Wasserstein distance in some situations (see [26, Proposition 5]). In particular, the results of [26, Proposition 5]) imply that for commutative Gibbs states a TC inequality for $W_{1,\nabla}$ implies the corresponding statement for $W_{1,\square}$.

## B.2 Local evolution of Lipschitz observables

As already discussed in Subsections B.1.1 and B.1.2 when we defined $\|\cdot\|_{\text{Lip},\nabla}$ and $\|\cdot\|_{\text{Lip},\square}$, Lipschitz constants can be easily controlled by making assumptions on the locality of the operators. Indeed, 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. More precisely, in [56, 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)\|_{\text{Lip},\square} \le 2|L|\|O\|_{\text{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 $\|\cdot\|_{\text{Lip},\nabla}$ and $\|\cdot\|_{\text{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 and 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 various formulations of Lieb-Robinson bounds. Reviewing those in detail is beyond the scope of this work and we refer to [9, 38, 47, 57] and references therein for more details. For studying the behaviour of $\|\cdot\|_{\text{Lip},\nabla}$ under local evolutions, the most natural formulation is the commutator version: the generator $\mathcal{L}$ of a quasi-local dynamics $t \mapsto \Phi_t = e^{t\mathcal{L}}$ on $n$ qudits arranged on a graph $G = (V, E)$, with graph distance dist, is said to satisfy a LR bound with LR velocity $v$ if for any observable $O_A$ supported on a region $A$ and any other observable $B$ supported on a region $B$, we have:

$$ \|[\Phi_t^*(O_A), O_B]\|_{\infty} \le c(e^{vt} - 1) g(\text{dist}(A, B)) \|O_A\|_{\infty} \|O_B\|_{\infty}, \qquad (LR1) $$

for $g: \mathbb{N} \to \mathbb{R}_+$ a function such that $\lim_{x \to \infty} g(x) = 0$. We then have:
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 notation, 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$ [21], 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 the heat semigroup. We refer to [21, 59] 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:

$$ \|X\|_{\text{Lip},\nabla} = \left( \sum_{k \in \mathcal{K}} \|P_k X P_k - X\|_{\infty}^2 \right)^{1/2}, \qquad (39) $$

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 [21].

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

$$ W_{1,\nabla}(\rho, \sigma) := \sup_{X=X^{\dagger}, \|X\|_{\text{Lip},\nabla} \le 1} |\text{Tr}(X(\rho - \sigma))|. $$

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

$$ |\text{tr}[X(\rho - \sigma)]| \le \|X\|_{\text{Lip},\nabla} W_{1,\nabla}(\rho, \sigma). \qquad (40) $$

Although this geometric interpretation opens up other interesting mathematical connections for this definition, the differential Wasserstein distance has the downside of being state dependent. It however induces a stronger topology than the quantum Hamming Wasserstein distance in some situations (see [26, Proposition 5]). In particular, the results of [26, Proposition 5]) imply that for commutative Gibbs states a TC inequality for $W_{1,\nabla}$ implies the corresponding statement for $W_{1,\square}$.

## B.2 Local evolution of Lipschitz observables

As already discussed in Subsections B.1.1 and B.1.2 when we defined $\|\cdot\|_{\text{Lip},\nabla}$ and $\|\cdot\|_{\text{Lip},\square}$, Lipschitz constants can be easily controlled by making assumptions on the locality of the operators. Indeed, 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. More precisely, in [56, 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)\|_{\text{Lip},\square} \le 2|L|\|O\|_{\text{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 $\|\cdot\|_{\text{Lip},\nabla}$ and $\|\cdot\|_{\text{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 and 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 various formulations of Lieb-Robinson bounds. Reviewing those in detail is beyond the scope of this work and we refer to [9, 38, 47, 57] and references therein for more details. For studying the behaviour of $\|\cdot\|_{\text{Lip},\nabla}$ under local evolutions, the most natural formulation is the commutator version: the generator $\mathcal{L}$ of a quasi-local dynamics $t \mapsto \Phi_t = e^{t\mathcal{L}}$ on $n$ qudits arranged on a graph $G = (V, E)$, with graph distance dist, is said to satisfy a LR bound with LR velocity $v$ if for any observable $O_A$ supported on a region $A$ and any other observable $B$ supported on a region $B$, we have:

$$ \|[\Phi_t^*(O_A), O_B]\|_{\infty} \le c(e^{vt} - 1) g(\text{dist}(A, B)) \|O_A\|_{\infty} \|O_B\|_{\infty}, \qquad (\text{LR1}) $$

for $g: \mathbb{N} \to \mathbb{R}_+$ a function such that $\lim_{x \to \infty} g(x) = 0$. We then have:
**Proposition B.1** (Growth of differential Lipschitz constant for local evolutions). Let $(\nabla, \sigma)$ be a differential structure on $\mathcal{M}_{d^n}$ and let $O = \sum_i O_i$ be an observable with $\|O_i\|_{\infty} \le 1$. Let $A_i$ denote the support of each $O_i$ and $B_j$ that of each $L_j$. Moreover, let $t \mapsto \Phi_t$ be an evolution satisfying Eq. (LR1) and set $o(i, j)(t) = 2$ if $A_i \cap B_j \ne \emptyset$ and $c(e^{vt} - 1)g(\text{dist}(A_i, B_j))$ else. Then:

$$ \|\Phi_t^*(O)\|_{\text{Lip},\nabla}^2 \le \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:

$$ \|\Phi_t^*(O)\|_{\text{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:

$$ \|[\Phi_t^*(O), L_i]\|_{\infty}^2 \le \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]\|_{\infty} \le 2$ by the submultiplicativity of the operator norm, a triangle inequality and $\|L_j\|_{\infty} \le 1$. In case $O_i$ and $L_j$ do not overlap, the stronger bound in Eq. (LR1) holds. $\square$

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, \dots, n - k$, and $g(\text{dist}(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$:

$$ \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) \le k + \frac{(e^{vt} - 1)e^{-\mu}}{1 - e^{-\mu}}. $$

Thus,

$$ \|\Phi_t^*(O)\|_{\text{Lip},\nabla} \le \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 somehow stronger, version of Prop. B.1 for $\|\cdot\|_{\text{Lip},\square}$. In some situations, bounds like (LR1) can be further exploited in order to prove the quasi-locality of Markovian quantum dynamics [9]: for any region $C \subset D \subset V$, there exists a local Markovian evolution $t \mapsto \Phi_t^{(D)}$ that acts non-trivially only on region $D$, and such that for any observable $O_C$ supported on region $C$,

$$ \|(\Phi_t^* - (\Phi_t^{(D)})^*)(O_C)\|_{\infty} \le c' (e^{vt} - 1) h(d(C, V \setminus D)) \|O_C\|_{\infty}, \quad (\text{LR2}) $$

for some other function $h: \mathbb{N} \to \mathbb{R}_+$ such that $\lim_{x \to \infty} h(x) = 0$ and constant $c' > 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 [56] given a Lieb-Robinson bound:

**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, \dots, n\}$ of the graph:

$$ W_{1,\square}(\Phi_t(\rho), \Phi_t(\sigma)) \le \left(8 + 2c'(e^{vt} - 1) \sum_{i=3}^{n} h(d(\{i \cdots n\}, \{1\}))\right) W_{1,\square}(\rho, \sigma). \quad (41) $$

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

$$ \|\Phi_t^*(H)\|_{\text{Lip},\square} \le \left(8 + 2c'(e^{vt} - 1) \sum_{i=3}^{n} h(d(\{i \cdots n\}, \{1\}))\right) \|H\|_{\text{Lip},\square}. \quad (42) $$
**Proposition B.1** (Growth of differential Lipschitz constant for local evolutions). Let $(\nabla, \sigma)$ be a differential structure on $\mathcal{M}_{d^n}$ and let $O = \sum_i O_i$ be an observable with $\|O_i\|_{\infty} \le 1$. Let $A_i$ denote the support of each $O_i$ and $B_j$ that of each $L_j$. Moreover, let $t \mapsto \Phi_t$ be an evolution satisfying Eq. (LR1) and set $o(i, j)(t) = 2$ if $A_i \cap B_j \ne \emptyset$ and $c(e^{vt} - 1)g(\text{dist}(A_i, B_j))$ else. Then:

$$ \|\Phi_t^*(O)\|_{\text{Lip},\nabla}^2 \le \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:

$$ \|\Phi_t^*(O)\|_{\text{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:

$$ \|[\Phi_t^*(O), L_i]\|_{\infty}^2 \le \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]\|_{\infty} \le 2$ by the submultiplicativity of the operator norm, a triangle inequality and $\|L_j\|_{\infty} \le 1$. In case $O_i$ and $L_j$ do not overlap, the stronger bound in Eq. (LR1) holds. $\square$

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, \dots, n - k$, and $g(\text{dist}(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$:

$$ \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) \le k + \frac{(e^{vt} - 1)e^{-\mu}}{1 - e^{-\mu}}. $$

Thus,

$$ \|\Phi_t^*(O)\|_{\text{Lip},\nabla} \le \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 somehow stronger, version of Prop. B.1 for $\|\cdot\|_{\text{Lip},\square}$. In some situations, bounds like (LR1) can be further exploited in order to prove the quasi-locality of Markovian quantum dynamics [9]: for any region $C \subset D \subset V$, there exists a local Markovian evolution $t \mapsto \Phi_t^{(D)}$ that acts non-trivially only on region $D$, and such that for any observable $O_C$ supported on region $C$,

$$ \|(\Phi_t^* - (\Phi_t^{(D)})^*)(O_C)\|_{\infty} \le c' (e^{vt} - 1) h(d(C, V \setminus D)) \|O_C\|_{\infty}, \quad (\text{LR2}) $$

for some other function $h: \mathbb{N} \to \mathbb{R}_+$ such that $\lim_{x \to \infty} h(x) = 0$ and constant $c' > 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 [56] given a Lieb-Robinson bound:

**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, \dots, n\}$ of the graph:

$$ W_{1,\square}(\Phi_t(\rho), \Phi_t(\sigma)) \le \left(8 + 2c'(e^{vt} - 1) \sum_{i=3}^{n} h(d(\{i \cdots n\}, \{1\}))\right) W_{1,\square}(\rho, \sigma). \quad (41) $$

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

$$ \|\Phi_t^*(H)\|_{\text{Lip},\square} \le \left(8 + 2c'(e^{vt} - 1) \sum_{i=3}^{n} h(d(\{i \cdots n\}, \{1\}))\right) \|H\|_{\text{Lip},\square}. \quad (42) $$
Proof. From [56], 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 $B_n$, which in turn is the convex hull of the set of the differences between couples of neighboring quantum states:

$$ \mathcal{N}_n = \bigcup_{i \in V} \mathcal{N}_n^{(i)}, \quad \mathcal{N}_n^{(i)} = \{\rho - \sigma : \rho, \sigma \in \mathcal{D}_{d^n}, \text{tr}_i(\rho) = \text{tr}_i(\sigma)\}. $$

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

$$ \|\Phi_t\|_{W_1 \to W_1} = \max \left( \|\Phi_t(X)\|_{W_1} : X \in \mathcal{M}_{d^n}^{\text{sa},0}, \|X\|_{W_1} \le 1 \right) = \max_{X \in \mathcal{N}_n} \|\Phi_t(X)\|_{W_1}, $$

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

$$ \begin{aligned} \|\Phi_t(X)\|_{W_1} &\le \frac{1}{2\sqrt{n}} \sum_{i=1}^n \left\| \frac{I}{d^{i-1}} \otimes \text{tr}_{1 \dots i-1} \circ \Phi_t(X) - \frac{I}{d^i} \otimes \text{tr}_{1 \dots i} \circ \Phi_t(X) \right\|_1 \\ &= \frac{1}{2\sqrt{n}} \sum_{i=1}^n \left\| \int d\mu(U_i) \text{tr}_{1 \dots i-1} \circ (\Phi_t(X) - U_i \Phi_t(X) U_i^\dagger) \right\|_1 \\ &\le \frac{1}{2\sqrt{n}} \sum_{i=1}^n \int d\mu(U_i) \|[U_i, \text{tr}_{1 \dots i-1} \circ \Phi_t(X)]\|_1 \\ &\le \frac{1}{\sqrt{n}} \sum_{i=1}^n \| \text{tr}_{1 \dots i-1} \circ \Phi_t(X) \|_1 \\ &\stackrel{(1)}{=} \frac{1}{\sqrt{n}} \sum_{i=1}^n \| \text{tr}_{1 \dots i-1} \circ (\Phi_t - \Phi_t^{(i-k \cdots n)})(X) \|_1 \end{aligned} \quad (43) $$

where $\mu$ denotes the Haar measure on one qudit, and where (1) follows from the fact that $\text{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 $i \ge 3$ that

$$ \begin{aligned} \| \text{tr}_{1 \dots i-1} \circ (\Phi_t - \Phi_t^{(i-k \cdots n)})(X) \|_1 &= \max_{\|O_{i \dots n}\|_{\infty} \le 1} \left| \text{tr} [X(\Phi_t^* - \Phi_t^{(i-k \cdots n)*})(O_{i \dots n})] \right| \\ &\le \max_{\|O_{i \dots n}\|_{\infty} \le 1} \|(\Phi_t^* - \Phi_t^{(i-k \cdots n)*})(O_{i \dots n})\|_{\infty} \|X\|_1 \\ &\le c' (e^{vt} - 1) h(\text{dist}(\{i \dots n\}, \{1 \dots i - k - 1\})) \|X\|_1 \\ &\stackrel{(2)}{\le} 2c' (e^{vt} - 1) h(\text{dist}(\{i \dots n\}, \{1 \dots i - k - 1\})) \sqrt{n} \|X\|_{W_1}, \end{aligned} $$

where (2) follows from [56, Proposition 6]. By picking $k = i - 2$ and inserting this estimate into Eq. (43) for $i \ge 3$ and the trivial estimate $\| \text{tr}_{1 \dots i-1} \circ (\Phi_t - \Phi_t^{(i-k \cdots n)})(X) \|_1 \le 2 \|X\|_1$ for $i = 1, 2$ we obtain Eq. (41). Eq. (42) follows by duality. $\square$

## B.3 Transportation cost inequalities

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

$$ W_1(\rho, \sigma) \le \sqrt{\frac{1}{2\alpha} D(\rho||\sigma)}, \quad (44) $$

where $W_1 \in \{W_{1,\square}, W_{1,\nabla}, W_{1,\text{loc}}\}$. In what follows, we simply write $||.||_{\text{Lip}}$ and $W_1$ to denote either of the Lipschitz constants, and their corresponding Wasserstein metrics, 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 functional inequality called the modified logarithmic Sobolev inequality. It is
Proof. From [56], 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 $B_n$, which in turn is the convex hull of the set of the differences between couples of neighboring quantum states:

$$ \mathcal{N}_n = \bigcup_{i \in V} \mathcal{N}_n^{(i)}, \quad \mathcal{N}_n^{(i)} = \{ \rho - \sigma : \rho, \sigma \in \mathcal{D}_{d^n}, \text{tr}_i(\rho) = \text{tr}_i(\sigma) \}. $$

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

$$ \|\Phi_t\|_{W_1 \to W_1} = \max \left( \|\Phi_t(X)\|_{W_1} : X \in \mathcal{M}_{d^n}^{\text{sa},0}, \|X\|_{W_1} \le 1 \right) = \max_{X \in \mathcal{N}_n} \|\Phi_t(X)\|_{W_1}, $$

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

$$ \begin{align} \|\Phi_t(X)\|_{W_1} &\le \frac{1}{2\sqrt{n}} \sum_{i=1}^n \left\| \frac{I}{d^{i-1}} \otimes \text{tr}_{1\dots i-1} \circ \Phi_t(X) - \frac{I}{d^i} \otimes \text{tr}_{1\dots i} \circ \Phi_t(X) \right\|_1 \nonumber \\ &= \frac{1}{2\sqrt{n}} \sum_{i=1}^n \left\| \int d\mu(U_i) \text{tr}_{1\dots i-1} \circ (\Phi_t(X) - U_i \Phi_t(X) U_i^\dagger) \right\|_1 \nonumber \\ &\le \frac{1}{2\sqrt{n}} \sum_{i=1}^n \int d\mu(U_i) \|[U_i, \text{tr}_{1\dots i-1} \circ \Phi_t(X)]\|_1 \nonumber \\ &\le \frac{1}{\sqrt{n}} \sum_{i=1}^n \| \text{tr}_{1\dots i-1} \circ \Phi_t(X) \|_1 \nonumber \\ &\stackrel{(1)}{=} \frac{1}{\sqrt{n}} \sum_{i=1}^n \| \text{tr}_{1\dots i-1} \circ (\Phi_t - \Phi_t^{(i-k\dots n)})(X) \|_1 \tag{43} \end{align} $$

where $\mu$ denotes the Haar measure on one qudit, and where (1) follows from the fact that $\text{tr}_1(X) = 0$, with $\Phi_t^{(i-k\dots n)} \equiv \Phi_t^{\{i-k,\dots,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 $i \ge 3$ that

$$ \begin{align*} \| \text{tr}_{1 \dots i-1} \circ (\Phi_t - \Phi_t^{(i-k \dots n)})(X) \|_1 &= \max_{\|O_{i \dots n}\|_{\infty} \le 1} \left| \text{tr} \left[ X (\Phi_t^* - \Phi_t^{(i-k \dots n)*}) (O_{i \dots n}) \right] \right| \\ &\le \max_{\|O_{i \dots n}\|_{\infty} \le 1} \left\| (\Phi_t^* - \Phi_t^{(i-k \dots n)*}) (O_{i \dots n}) \right\|_{\infty} \|X\|_1 \\ &\le c' (e^{vt} - 1) h(\text{dist}(\{i \dots n\}, \{1 \dots i - k - 1\})) \|X\|_1 \\ &\stackrel{(2)}{\le} 2c' (e^{vt} - 1) h(\text{dist}(\{i \dots n\}, \{1 \dots i - k - 1\})) \sqrt{n} \|X\|_{W_1}, \end{align*} $$

where (2) follows from [56, Proposition 6]. By picking $k = i - 2$ and inserting this estimate into Eq. (43) for $i \ge 3$ and the trivial estimate $\| \text{tr}_{1\dots i-1} \circ (\Phi_t - \Phi_t^{(i-k\dots n)})(X) \|_1 \le 2\|X\|_1$ for $i = 1, 2$ we obtain Eq. (41). Eq. (42) follows by duality. $\square$

## B.3 Transportation cost inequalities

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

$$ W_1(\rho, \sigma) \le \sqrt{\frac{1}{2\alpha} D(\rho \| \sigma)}, \qquad (44) $$

where $W_1 \in \{W_{1,\square}, W_{1,\nabla}, W_{1,\text{loc}}\}$. In what follows, we simply write $||.||_{\text{Lip}}$ and $W_1$ to denote either of the Lipschitz constants, and their corresponding Wasserstein metrics, 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 functional inequality called the modified logarithmic Sobolev inequality. It is
beyond the scope of this paper to explain this connection and we refer to e.g. [59] and references therein for a discussion on these topics. But for our purposes it is important to note that in [18] the authors and Capel show modified logarithmic Sobolev inequalities for several classes of Gibbs states. More recently, one of the authors and De Palma derived such transportation cost inequalities for $W_{1,□}$ in [26]. In Theorem B.1 below we summarize the regimes for which transportation cost inequalities are known to hold:

**Theorem B.1** (transportation cost for commuting Hamiltonians [8, 18, 26]). Let $E_1, \dots, E_m \subset \mathcal{M}_{d^n}$ be a set of $k$-local linearly independent commuting operators with $\|E_i\|_\infty \le 1$. Then $\sigma(\lambda)$ satisfies a transportation cost inequality with constant $\alpha > 0$ for all $\lambda \in B_{\ell_\infty}(0, 1)$ in the following cases:

(i) The $E_i$ are classical or nearest neighbour (i.e. $k = 2$) on a regular lattice and the inverse temperature $\beta < \beta_c$ where $\beta_c$ only depends on $k$ and the dimension of the lattice, for both $W_{1,\nabla}, W_{1,□}$ and $\alpha = \Omega(1)$ [18].

(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 for $W_{1,□}$ and $\alpha = \Omega(1)$ [26, Theorem 3, Proposition 4].

(iii) The $E_i$ are one-dimensional and $\beta > 0$, for both $W_{1,\nabla}$ and $W_{1,□}$ and $\alpha = \Omega(\log(n)^{-1})$ [8].

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 Hamiltonians known to have exponential decay of correlations.

**Remark B.1.** In [26, Proposition 5], the authors show that $W_{1,□} \le 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,□}$ up to $c(k)$. Thus, although the authors of [8, 18] only obtain the result for $W_{1,\nabla}$, we can use it to translate it to $W_{1,□}$. We conclude that for commuting Hamiltonians TC are available for essentially all classes of local Hamiltonians for which they are expected to hold.

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

With these tools at hand, we are now ready to show that by 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 [22, 24, 40, 43]. 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, \dots, O_M \in \mathcal{M}_{2^n}$ such that $\|O_i\|_\infty \le 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.

For instance, for the shadow methods of [40], we obtain a $\mathcal{O}(4^k \log(M\delta^{-1})\epsilon^{-2})$ scaling by measuring in 1-qubit random bases. The estimate is then obtained by an appropriate median-of-means procedure for the expectation value of each output string. The computation for obtaining the expectation value of $E_i$ through this method then entails evaluating the expectation value of the observables on $\mathcal{O}(4^k \log(M\delta^{-1})\epsilon^{-2})$ product states. For $k$-local observables $E_i$, evaluating the expectation value of $E_i$ on a product state takes time $\mathcal{O}(e^{ck \log(M\delta^{-1})\epsilon^{-2}})$ for some $c > 0$. Thus, we see that for $k = \mathcal{O}(\log(n))$ also the postprocessing can be done efficiently.

The application of such results to maximum entropy methods is then clear: given $E_1, \dots, E_m$ assumed to be at most $k$-local, with probability at least $1 - \delta$ we can obtain an estimate $e'(\lambda)$ of $e(\lambda)$ satisfying:

$$ \|e'(\lambda) - e(\lambda)\|_{\ell_p} \le m^{\frac{1}{p}} \epsilon \quad (45) $$

using $\mathcal{O}(4^k \log(m\delta^{-1})\epsilon^{-2})$ copies of $\sigma(\lambda)$. We then finally obtain our main theorem, Theorem 1.1, restated here for the sake of clarity:

**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 $\mu \mapsto \log(\mathcal{Z}(\mu))$ is $L, U$ strongly convex in $B_{\ell_\infty}(0, 1)$. Then

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

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

$$ |\operatorname{tr}[O(\sigma(\lambda) - \sigma(\mu))]| \le \epsilon n^{\frac{1}{2}} \|O\|_{\mathrm{Lip}} \quad (46) $$

for all Lipschitz observables $O$ with probability at least $1 - \delta$.
beyond the scope of this paper to explain this connection and we refer to e.g. [59] and references therein for a discussion on these topics. But for our purposes it is important to note that in [18] the authors and Capel show modified logarithmic Sobolev inequalities for several classes of Gibbs states. More recently, one of the authors and De Palma derived such transportation cost inequalities for $W_{1,□}$ in [26]. In Theorem B.1 below we summarize the regimes for which transportation cost inequalities are known to hold:

**Theorem B.1** (transportation cost for commuting Hamiltonians [8, 18, 26]). Let $E_1, \dots, E_m \subset \mathcal{M}_{d^n}$ be a set of $k$-local linearly independent commuting operators with $\|E_i\|_{\infty} \le 1$. Then $\sigma(\lambda)$ satisfies a transportation cost inequality with constant $\alpha > 0$ for all $\lambda \in B_{\ell_{\infty}}(0, 1)$ in the following cases:

(i) The $E_i$ are classical or nearest neighbour (i.e. $k = 2$) on a regular lattice and the inverse temperature $\beta < \beta_c$ where $\beta_c$ only depends on $k$ and the dimension of the lattice, for both $W_{1,\nabla}, W_{1,□}$ and $\alpha = \Omega(1)$ [18].

(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 for $W_{1,□}$ and $\alpha = \Omega(1)$ [26, Theorem 3, Proposition 4].

(iii) The $E_i$ are one-dimensional and $\beta > 0$, for both $W_{1,\nabla}$ and $W_{1,□}$ and $\alpha = \Omega(\log(n)^{-1})$ [8].

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 Hamiltonians known to have exponential decay of correlations.

**Remark B.1.** In [26, Proposition 5], the authors show that $W_{1,□} \le 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,□}$ up to $c(k)$. Thus, although the authors of [8, 18] only obtain the result for $W_{1,\nabla}$, we can use it to translate it to $W_{1,□}$. We conclude that for commuting Hamiltonians TC are available for essentially all classes of local Hamiltonians for which they are expected to hold.

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

With these tools at hand, we are now ready to show that by 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 [22, 24, 40, 43]. 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, \dots, O_M \in \mathcal{M}_{2^n}$ such that $\|O_i\|_{\infty} \le 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.

For instance, for the shadow methods of [40], we obtain a $\mathcal{O}(4^k \log(M\delta^{-1})\epsilon^{-2})$ scaling by measuring in 1-qubit random bases. The estimate is then obtained by an appropriate median-of-means procedure for the expectation value of each output string. The computation for obtaining the expectation value of $E_i$ through this method then entails evaluating the expectation value of the observables on $\mathcal{O}(4^k \log(M\delta^{-1})\epsilon^{-2})$ product states. For $k$-local observables $E_i$, evaluating the expectation value of $E_i$ on a product state takes time $\mathcal{O}(e^{ck \log(M\delta^{-1})\epsilon^{-2}})$ for some $c > 0$. Thus, we see that for $k = \mathcal{O}(\log(n))$ also the postprocessing can be done efficiently.

The application of such results to maximum entropy methods is then clear: given $E_1, \dots, E_m$ assumed to be at most $k$-local, with probability at least $1 - \delta$ we can obtain an estimate $e'(\lambda)$ of $e(\lambda)$ satisfying:

$$ \|e'(\lambda) - e(\lambda)\|_{\ell_p} \le m^{\frac{1}{p}} \epsilon \quad (45) $$

using $\mathcal{O}(4^k \log(m\delta^{-1})\epsilon^{-2})$ copies of $\sigma(\lambda)$. We then finally obtain our main theorem, Theorem 1.1, restated here for the sake of clarity:

**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 $\mu \mapsto \log(\mathcal{Z}(\mu))$ is $L, U$ strongly convex in $B_{\ell_{\infty}}(0, 1)$. Then

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

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

$$ |\operatorname{tr}[O(\sigma(\lambda) - \sigma(\mu))]| \le \epsilon n^{\frac{1}{2}} \|O\|_{\mathrm{Lip}} \quad (46) $$

for all Lipschitz observables $O$ with probability at least $1 - \delta$.
Proof. Using the aforementioned methods of [40] we can obtain an estimate $e'(\lambda)$ satisfying the guarantee of Eq. (29) with probability at least $1 - \delta$. We will now resort to the results of Thm. A.1 to obtain guarantees on the output of the maximum entropy algorithm. Now we solve the maximum entropy problem with $e'(\lambda)$ and set the stopping criterion for the output $\mu_*$ as

$$ \|e'(\lambda) - e(\mu_*)\|_{\ell_2} < (4c + 1) \epsilon \sqrt{m}. $$

for $c > 10$. Then it follows from Thm. A.1 that:

$$ D(\sigma(\mu_*), \sigma(\lambda)) \le D(\sigma(\mu_*), \sigma(\lambda)) + D(\sigma(\lambda), \sigma(\mu_*)) = \mathcal{O}(\beta \min\{\beta L^{-1} \epsilon^2 m, \epsilon m\}). \quad (47) $$

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

$$ |\operatorname{tr}[O(\sigma(\lambda) - \sigma(\mu_*))]| \le \|O\|_{\mathrm{Lip}} W_1(\sigma(\mu_*), \sigma(\lambda)) \le \|O\|_{\mathrm{Lip}} \sqrt{\frac{D(\sigma(\mu_*), \sigma(\lambda))}{2\alpha}}. $$

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

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

We conclude by suitably rescaling $\epsilon$.

<figure><img src="image_7.png" /></figure>

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 Lipschitz constant that satisfies $\|O\|_{\mathrm{Lip}} = \mathcal{O}(l\sqrt{n})$. Shadow-like methods would require $\mathcal{O}(e^{\mathcal{O}(l)}) \log(m\delta^{-1}) \epsilon^{-2})$ samples to learn such observables up to a relative error of $n\epsilon$. Our methods, however, require $\mathcal{O}(e^{\mathcal{O}(k)}) \operatorname{poly}(l, \epsilon^{-1}) \log(m\delta^{-1}))$, which yields exponential speedups in the regime $l = \operatorname{poly}(\log(n))$. 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 than solving the maximum entropy problem for obtaining the same precision.

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

## C.1 Results for shallow circuits

Let us be 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 that models the interactions of a unitary circuit and suppose we implement $L$ layers of an unknown unitary circuit consisting of 1 and 2-qubit unitaries laid out according to $G$. That is, we have an unknown shallow circuit $\mathcal{U}$ of depth $L$ with respect to $G$. More precisely,

$$ \mathcal{U} = \prod_{\ell \in [L]} \bigotimes_{e \in \mathcal{E}_{\ell}} \mathcal{U}_{\ell, e}, \quad (48) $$

where each $\mathcal{E}_{\ell} \subset E$ are a subset of the edges such that any $e, e' \in \mathcal{E}_{\ell}$ do not share a vertex. Our goal is to show how to approximately learn the state

$$ |\psi\rangle = \mathcal{U}|0\rangle^{\otimes n}. $$

The overall idea consists in finding a Gibbs state approximating $|\psi\rangle$ in Wasserstein distance. We will then find a differential structure for (approximations) of shallow circuits and then showing that the (approximation of the) output satisfies a TC inequality with respect to it. Thus, it suffices to control the relative entropy with this approximation to ensure a good approximation in Wasserstein distance.
Proof. Using the aforementioned methods of [40] we can obtain an estimate $e'(\lambda)$ satisfying the guarantee of Eq. (29) with probability at least $1 - \delta$. We will now resort to the results of Thm. A.1 to obtain guarantees on the output of the maximum entropy algorithm. Now we solve the maximum entropy problem with $e'(\lambda)$ and set the stopping criterion for the output $\mu_*$ as

$$ \|e'(\lambda) - e(\mu_*)\|_{\ell_2} < (4c + 1) \epsilon \sqrt{m}. $$

for $c > 10$. Then it follows from Thm. A.1 that:

$$ D(\sigma(\mu_*), \sigma(\lambda)) \le D(\sigma(\mu_*), \sigma(\lambda)) + D(\sigma(\lambda), \sigma(\mu_*)) = \mathcal{O}(\beta \min\{\beta L^{-1} \epsilon^2 m, \epsilon m\}). \quad (47) $$

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

$$ |\operatorname{tr}[O(\sigma(\lambda) - \sigma(\mu_*))]| \le \|O\|_{\mathrm{Lip}} W_1(\sigma(\mu_*), \sigma(\lambda)) \le \|O\|_{\mathrm{Lip}} \sqrt{\frac{D(\sigma(\mu_*), \sigma(\lambda))}{2\alpha}}. $$

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

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

We conclude by suitably rescaling $\epsilon$.

<figure><img src="image_8.png" /></figure>

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 Lipschitz constant that satisfies $\|O\|_{\mathrm{Lip}} = \mathcal{O}(l\sqrt{n})$. Shadow-like methods would require $\mathcal{O}(e^{\mathcal{O}(l)}) \log(m\delta^{-1}) \epsilon^{-2})$ samples to learn such observables up to a relative error of $n\epsilon$. Our methods, however, require $\mathcal{O}(e^{\mathcal{O}(k)}) \operatorname{poly}(l, \epsilon^{-1}) \log(m\delta^{-1}))$, which yields exponential speedups in the regime $l = \operatorname{poly}(\log(n))$. 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 than solving the maximum entropy problem for obtaining the same precision.

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

## C.1 Results for shallow circuits

Let us be 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 that models the interactions of a unitary circuit and suppose we implement $L$ layers of an unknown unitary circuit consisting of 1 and 2-qubit unitaries laid out according to $G$. That is, we have an unknown shallow circuit $\mathcal{U}$ of depth $L$ with respect to $G$. More precisely,

$$ \mathcal{U} = \prod_{\ell \in [L]} \bigotimes_{e \in \mathcal{E}_{\ell}} \mathcal{U}_{\ell, e}, \quad (48) $$

where each $\mathcal{E}_{\ell} \subset E$ are a subset of the edges such that any $e, e' \in \mathcal{E}_{\ell}$ do not share a vertex. Our goal is to show how to approximately learn the state

$$ |\psi\rangle = \mathcal{U}|0\rangle^{\otimes n}. $$

The overall idea consists in finding a Gibbs state approximating $|\psi\rangle$ in Wasserstein distance. We will then find a differential structure for (approximations) of shallow circuits and then showing that the (approximation of the) output satisfies a TC inequality with respect to it. Thus, it suffices to control the relative entropy with this approximation to ensure a good approximation in Wasserstein distance.
Let us find the appropriate approximation. First, note that for $\beta_{\epsilon} = \log(\epsilon^{-1})$ and $H_0 = -\sum_{i}^{n} Z_i$ where $Z_i$ denotes the Pauli matrix $Z = |0\rangle\langle 0| - |1\rangle\langle 1|$ acting on site $i$, we have that:

$$
\begin{align}
D(\mathcal{U}|0\rangle\langle 0|^{\otimes n} \mathcal{U}^{\dagger} \| \mathcal{U} \frac{e^{-\beta_{\epsilon} H_0}}{\text{tr}[e^{-\beta_{\epsilon} H_0}]} \mathcal{U}^{\dagger}) &= D(|0\rangle\langle 0|^{\otimes n} \| \frac{e^{-\beta_{\epsilon} H_0}}{\text{tr}[e^{-\beta_{\epsilon} H_0}]} ) \nonumber \\
&= nD(|0\rangle\langle 0| \| \frac{e^{\beta_{\epsilon} Z}}{\text{tr}[e^{\beta_{\epsilon} Z}]} ) \nonumber \\
&= -n\beta_{\epsilon} \langle 0|Z|0\rangle + n \log \text{tr} [e^{\beta_{\epsilon} Z}] \nonumber \\
&= n \log (1 + e^{-2\beta_{\epsilon}}) \nonumber \\
&\le ne^{-2\beta_{\epsilon}} \nonumber \\
&= \epsilon^2 n. \tag{49}
\end{align}
$$

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

$$W_1 \left( \mathcal{U}|0\rangle\langle 0|^{\otimes n} \mathcal{U}^{\dagger}, \mathcal{U} \frac{e^{-\beta_{\epsilon} H_0}}{\text{tr}[e^{-\beta_{\epsilon} H_0}]} \mathcal{U}^{\dagger} \right) \le \epsilon \sqrt{\frac{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}}{\text{tr}[e^{-\beta_{\epsilon} H_0}]} \mathcal{U}^{\dagger} = \frac{e^{-\beta_{\epsilon} H_{\mathcal{U}}}}{\text{tr}[e^{-\beta_{\epsilon} H_0}]}$. 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}}{\text{tr}[e^{-\beta_{\epsilon} H_0}]}$. 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}}{\text{tr}[e^{-\beta_{\epsilon} H_0}]} = \tau_{p^{\otimes n}}$ with $p_{\epsilon} = \frac{e^{\beta_{\epsilon}}}{e^{\beta_{\epsilon}} + e^{-\beta_{\epsilon}}}$. Moreover, let

$$a_i = I^{\otimes i-1} \otimes |0\rangle\langle 1| \otimes I^{\otimes n-i-1}$$

be the annihilation operator acting on qubit $i$. Defining $L_{i,0} = (p(1-p))^{\frac{1}{4}} a_i$ and $L_{i,1} = L_{i,0}^{\dagger}$, we get a differential structure for $\tau_{p}^{\otimes n}$ with $\{L_{i,k}, \omega_{i,k}\}$ for $i=1, \dots, n$, $k=0, 1$, $\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:

$$
\begin{align*}
\|O\|_{\text{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) \\
&= \sum_{i=1}^{n} \|[O, a_i]\|^2 + \|[O, a_i^{\dagger}]\|^2.
\end{align*}
$$

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 [11, 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 C.1** (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 $|\psi\rangle = \mathcal{U}|0\rangle^{\otimes n}$. Define for $\epsilon > 0$ $H_{\mathcal{U}} = -\sum_{i}^{n} \mathcal{U} Z_i \mathcal{U}^{\dagger}$ and $\sigma(\mathcal{U}, \epsilon) = \frac{e^{-\beta_{\epsilon} H_{\mathcal{U}}}}{\text{tr}[e^{-\beta_{\epsilon} H_{\mathcal{U}}}]}$ with $\beta_{\epsilon} = \log(\epsilon^{-1})$. Then for any state $\rho$ and all observables $O$ we have:

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

with

$$\|O\|_{\text{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, \quad a_i = I^{i-1} \otimes |0\rangle\langle 1| \otimes I^{n-i-1}$$
Let us find the appropriate approximation. First, note that for $\beta_{\epsilon} = \log(\epsilon^{-1})$ and $H_0 = -\sum_{i}^{n} Z_i$ where $Z_i$ denotes the Pauli matrix $Z = |0\rangle\langle 0| - |1\rangle\langle 1|$ acting on site $i$, we have that:

$$ 
\begin{aligned}
D(\mathcal{U}|0\rangle\langle 0|^{\otimes n}\mathcal{U}^{\dagger} \| \mathcal{U} \frac{e^{-\beta_{\epsilon} H_0}}{\text{tr}[e^{-\beta_{\epsilon} H_0}]\mathcal{U}^{\dagger}}) &= D(|0\rangle\langle 0|^{\otimes n} \| \frac{e^{-\beta_{\epsilon} H_0}}{\text{tr}[e^{-\beta_{\epsilon} H_0}]\mathcal{U}^{\dagger}}) \\
&= nD(|0\rangle\langle 0| \| \frac{e^{\beta_{\epsilon} Z}}{\text{tr}[e^{\beta_{\epsilon} Z}]\mathcal{U}^{\dagger}}) \\
&= -n\beta_{\epsilon} \langle 0|Z|0\rangle + n \log \text{tr}[e^{\beta_{\epsilon} Z}] \\
&= n \log (1 + e^{-2\beta_{\epsilon}}) \\
&\le ne^{-2\beta_{\epsilon}} \\
&= \epsilon^2 n.
\end{aligned} 
\quad (49) $$

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

$$ W_1 \left( \mathcal{U}|0\rangle\langle 0|^{\otimes n}\mathcal{U}^{\dagger}, \mathcal{U} \frac{e^{-\beta_{\epsilon} H_0}}{\text{tr}[e^{-\beta_{\epsilon} H_0}]\mathcal{U}^{\dagger}}\mathcal{U}^{\dagger} \right) \le \epsilon \sqrt{\frac{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}}{\text{tr}[e^{-\beta_{\epsilon} H_0}]\mathcal{U}^{\dagger}} \mathcal{U}^{\dagger} = \frac{e^{-\beta_{\epsilon} H_{\mathcal{U}}}}{\text{tr}[e^{-\beta_{\epsilon} H_0}]\mathcal{U}^{\dagger}}$. 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}}{\text{tr}[e^{-\beta_{\epsilon} H_0}]\mathcal{U}^{\dagger}}$. 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}}{\text{tr}[e^{-\beta_{\epsilon} H_0}]\mathcal{U}^{\dagger}} = \tau_{p}^{\otimes n}$ with $p_{\epsilon} = \frac{e^{\beta_{\epsilon}}}{e^{\beta_{\epsilon}} + e^{-\beta_{\epsilon}}}$. Moreover, let

$$ a_i = I^{\otimes i-1} \otimes |0\rangle\langle 1| \otimes I^{\otimes n-i-1} $$

be the annihilation operator acting on qubit $i$. Defining $L_{i,0} = (p(1-p))^{\frac{1}{4}} a_i$ and $L_{i,1} = L_{i,0}^{\dagger}$, we get a differential structure for $\tau_{p}^{\otimes n}$ with $\{L_{i,k}, \omega_{i,k}\}$ for $i = 1, \dots, n$, $k = 0, 1$, $\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:

$$ 
\begin{aligned}
\|O\|_{\text{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) \\
&= \sum_{i=1}^{n} \| [O, a_i] \|^2 + \| [O, a_i^{\dagger}] \|^2.
\end{aligned} 
$$

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 [11, 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 C.1** (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 $|\psi\rangle = \mathcal{U}|0\rangle^{\otimes n}$. Define for $\epsilon > 0$ $H_{\mathcal{U}} = -\sum_{i}^{n} \mathcal{U} Z_i \mathcal{U}^{\dagger}$ and $\sigma(\mathcal{U}, \epsilon) = \frac{e^{-\beta_{\epsilon} H_{\mathcal{U}}}}{\text{tr}[e^{-\beta_{\epsilon} H_{\mathcal{U}}}]\mathcal{U}^{\dagger}}$ with $\beta_{\epsilon} = \log(\epsilon^{-1})$. Then for any state $\rho$ and all observables $O$ we have:

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

with

$$ \|O\|_{\text{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, \quad a_i = I^{i-1} \otimes |0\rangle\langle 1| \otimes I^{n-i-1} $$
Proof. We have:

$$|\text{tr}[O(|\psi\rangle\langle\psi| - \rho)]| \le |\text{tr}[O(|\psi\rangle\langle\psi| - \sigma(U))]| + |\text{tr}[O(\sigma(U) - \rho)]|.$$

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

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

The result of Thm. C.1 and the discussion preceding it also give us a method of efficiently learning the outputs of shallow circuits, as illustrated by the following proposition:

**Proposition C.1** (Learning the outputs of shallow circuits). Let $U$ be an unknown $n$-qubit quantum circuit with known locality structure as in Eq. (48) and $|\psi\rangle = U|0\rangle^{\otimes n}$. Moreover, define $l_0$ as

$$l_0 = \max_{1 \le i \le n} \text{supp} |UZ_iU^\dagger|.$$

For some $\epsilon > 0$ we have that

$$\mathcal{O}(\epsilon^{-2} 4^{3l_0} \log^4(n4^{l_0} \epsilon^{-1}) \log(4^{l_0} n \delta^{-1})) \quad (50)$$

samples of $|\psi\rangle$ suffice to learn a Gibbs state $\sigma(\mu)$ that satisfies

$$W_{1,\nabla}(|\psi\rangle\langle\psi|, \sigma(\mu)) \le \sqrt{\epsilon n} \quad (51)$$

with probability of success at least $1 - \delta$. Along similar lines we have that

$$\mathcal{O}(\epsilon^{-2} 4^{3l_0} n^2 \log^4(n4^{l_0} \epsilon^{-1}) \log(4^{l_0} n \delta^{-1})) \quad (52)$$

samples suffice to learn a state $\sigma(\lambda)$ such that

$$\|\psi\rangle\langle\psi| - \sigma(\lambda)\|_{\text{tr}} \le \sqrt{\epsilon}.$$

Proof. Let $\mathcal{E}_i$ be a basis of Pauli operators for matrices on the support of $UZ_iU^\dagger$. By our assumption on $l_0$, we know that for each $i$ we have that $|\mathcal{E}_i| \le 4^{l_0}$. For simplicity we will assume that there are no Pauli words that are contained in two distinct $\mathcal{E}_i$ and we will enumerate all different Pauli words as $\{E_{i,j}\}$ for $1 \le i \le n$ and $1 \le j \le 4^{l_0}$ indicating the elements of the different $\mathcal{E}_i$. Thus, there is a $\lambda \in \mathbb{R}^m$ with $m \le n4^{l_0}$ and $\|\lambda\|_{\ell_\infty} \le 1$ such that

$$\sigma(U, \epsilon) \propto \exp \left( -\beta_{\epsilon} \sum_{i=1}^{n} \sum_{E_{i,j} \in \mathcal{E}_i} \lambda_{i,j} E_{i,j} \right).$$

Let $\tilde{\epsilon} > 0$ be given. It follows from Eq. (49) and Pinsker's inequality that picking $\epsilon = \frac{\tilde{\epsilon}}{n16}$ is sufficient to ensure that

$$\|\psi\rangle\langle\psi| - \sigma(U, \frac{\tilde{\epsilon}}{n4})\|_{\text{tr}} \le \frac{\tilde{\epsilon}}{2}. \quad (53)$$

Measuring

$$\mathcal{O}(4^{l_0} \tilde{\epsilon}^{-2} \log(4^{l_0} n \delta^{-1})) \quad (54)$$

copies of $|\psi\rangle$ is sufficient to obtain estimates of $\text{tr} [|\psi\rangle\langle\psi|E_{i,j}]$ up to an $\tilde{\epsilon}/4$ error. By Eq. (53) and a triangle inequality they are also $\tilde{\epsilon}/2$ away from $\text{tr}[\sigma(U, \frac{\tilde{\epsilon}}{n4})E_{i,j}]$.

Thus, running the maximum entropy principle with these estimates and the basis of operators given by $\cup_{i=1}^n \mathcal{E}_i$ will yield us an estimate $\sigma(\mu)$ that satisfies:

$$D(\sigma(\mu)\|\sigma(U, \frac{\tilde{\epsilon}}{n})) \le \log\left(\frac{n}{4\tilde{\epsilon}}\right) \tilde{\epsilon}m \le \log\left(\frac{n}{4\tilde{\epsilon}}\right) \tilde{\epsilon}4^{l_0}n.$$

To obtain the estimate in Wasserstein distance in Eq. (50), we can pick $\tilde{\epsilon} = \mathcal{O}(\epsilon/(4^{l_0} \log^2(n4^{l_0} \epsilon^{-1})))$, as in this case

$$\log\left(\frac{n}{4\tilde{\epsilon}}\right) \tilde{\epsilon}m \le \log\left(\frac{n}{4\tilde{\epsilon}}\right) \tilde{\epsilon}4^{l_0}n = \epsilon \frac{\log(n \log^2(n4^{l_0} \epsilon^{-1})\epsilon^{-1})}{\log^2(n4^{l_0} \epsilon^{-1})} n \le \epsilon n.$$

The claim then follows from the results of Thm. C.1 and substituting $\tilde{\epsilon}$ into Eq. (54).

For the statement on the sample complexity for the trace distance, we may pick $\tilde{\epsilon} = \mathcal{O}(\epsilon/(n4^{l_0} \log^2(n4^{l_0} \epsilon^{-1})))$, which yields the statement after applying Pinsker's inequality. $\square$
Proof. We have:

$$|\text{tr}[O(|\psi\rangle\langle\psi| - \rho)]| \le |\text{tr}[O(|\psi\rangle\langle\psi| - \sigma(U))]| + |\text{tr}[O(\sigma(U) - \rho)]|.$$

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

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

The result of Thm. C.1 and the discussion preceding it also give us a method of efficiently learning the outputs of shallow circuits, as illustrated by the following proposition:

**Proposition C.1** (Learning the outputs of shallow circuits). Let $U$ be an unknown $n$-qubit quantum circuit with known locality structure as in Eq. (48) and $|\psi\rangle = U|0\rangle^{\otimes n}$. Moreover, define $l_0$ as

$$l_0 = \max_{1 \le i \le n} \text{supp} |UZ_iU^\dagger|.$$

For some $\epsilon > 0$ we have that

$$\mathcal{O}(\epsilon^{-2} 4^{3l_0} \log^4(n4^{l_0} \epsilon^{-1}) \log(4^{l_0} n\delta^{-1})) \quad (50)$$

samples of $|\psi\rangle$ suffice to learn a Gibbs state $\sigma(\mu)$ that satisfies

$$W_{1,\nabla}(|\psi\rangle\langle\psi|, \sigma(\mu)) \le \sqrt{\epsilon n} \quad (51)$$

with probability of success at least $1 - \delta$. Along similar lines we have that

$$\mathcal{O}(\epsilon^{-2} 4^{3l_0} n^2 \log^4(n4^{l_0} \epsilon^{-1}) \log(4^{l_0} n\delta^{-1})) \quad (52)$$

samples suffice to learn a state $\sigma(\lambda)$ such that

$$\| |\psi\rangle\langle\psi| - \sigma(\lambda) \|_{\text{tr}} \le \sqrt{\epsilon}.$$

Proof. Let $\mathcal{E}_i$ be a basis of Pauli operators for matrices on the support of $UZ_iU^\dagger$. By our assumption on $l_0$, we know that for each $i$ we have that $|\mathcal{E}_i| \le 4^{l_0}$. For simplicity we will assume that there are no Pauli words that are contained in two distinct $\mathcal{E}_i$ and we will enumerate all different Pauli words as $\{E_{i,j}\}$ for $1 \le i \le n$ and $1 \le j \le 4^{l_0}$ indicating the elements of the different $\mathcal{E}_i$. Thus, there is a $\lambda \in \mathbb{R}^m$ with $m \le n4^{l_0}$ and $\|\lambda\|_{\ell_\infty} \le 1$ such that

$$\sigma(U, \epsilon) \propto \exp \left( -\beta_{\epsilon} \sum_{i=1}^{n} \sum_{E_{i,j} \in \mathcal{E}_i} \lambda_{i,j} E_{i,j} \right).$$

Let $\tilde{\epsilon} > 0$ be given. It follows from Eq. (49) and Pinsker's inequality that picking $\epsilon = \frac{\tilde{\epsilon}}{n16}$ is sufficient to ensure that

$$\| |\psi\rangle\langle\psi| - \sigma(U, \frac{\tilde{\epsilon}}{n4}) \|_{\text{tr}} \le \frac{\tilde{\epsilon}}{2}. \quad (53)$$

Measuring

$$\mathcal{O}(4^{l_0} \tilde{\epsilon}^{-2} \log(4^{l_0} n\delta^{-1})) \quad (54)$$

copies of $|\psi\rangle$ is sufficient to obtain estimates of $\text{tr} [|\psi\rangle\langle\psi| E_{i,j}]$ up to an $\tilde{\epsilon}/4$ error. By Eq. (53) and a triangle inequality they are also $\tilde{\epsilon}/2$ away from $\text{tr} [\sigma(U, \frac{\tilde{\epsilon}}{n4}) E_{i,j}]$.

Thus, running the maximum entropy principle with these estimates and the basis of operators given by $\cup_{i=1}^n \mathcal{E}_i$ will yield us an estimate $\sigma(\mu)$ that satisfies:

$$D(\sigma(\mu)\|\sigma(U, \frac{\tilde{\epsilon}}{n})) \le \log\left(\frac{n}{4\tilde{\epsilon}}\right) \tilde{\epsilon}m \le \log\left(\frac{n}{4\tilde{\epsilon}}\right) \tilde{\epsilon}4^{l_0}n.$$

To obtain the estimate in Wasserstein distance in Eq. (50), we can pick $\tilde{\epsilon} = \mathcal{O}(\epsilon/(4^{l_0} \log^2(n4^{l_0} \epsilon^{-1})))$, as in this case

$$\log\left(\frac{n}{4\tilde{\epsilon}}\right) \tilde{\epsilon}m \le \log\left(\frac{n}{4\tilde{\epsilon}}\right) \tilde{\epsilon}4^{l_0}n = \epsilon \frac{\log(n \log^2(n4^{l_0} \epsilon^{-1}) \epsilon^{-1})}{\log^2(n4^{l_0} \epsilon^{-1})} n \le \epsilon n.$$

The claim then follows from the results of Thm. C.1 and substituting $\tilde{\epsilon}$ into Eq. (54).

For the statement on the sample complexity for the trace distance, we may pick $\tilde{\epsilon} = \mathcal{O}(\epsilon/(4^{l_0} \log^2(n4^{l_0} \epsilon^{-1})))$, which yields the statement after applying Pinsker's inequality. $\square$
This shows that shallow circuits can be learned efficiently as long as $l_0 = O(\log(n))$. However, it is not immediately obvious how to also estimate the expectation values of the Gibbs states required to run the maximum entropy algorithm. Thus, at least with the methods presented here, the postprocessing takes exponential time in the number of qubits.

**Ground states of gapped systems:** in light of our results for shallow circuits, it is natural to ask to what extent our framework can be extended to ground states of gapped Hamiltonians, specially in 1D. Thus, let us briefly comment on the technical barriers on the way of such statements. First, notice that the statement of Thm C.1 required inverse temperatures scaling logarithmically in system size for the approximation of the ground state. Most known TC inequalities have an exponential scaling with the inverse temperature and, thus, TC at this inverse temperature the savings compared to Pinsker are not quadratic, hindering a straightforward extension to gapped systems. There are some nontrivial examples of ground states satisfying a TC inequality with the constant depending inverse linearly with the temperature, like graph states [63]. But, as they can also be prepared from a constant depth quantum circuit almost by definition, they fall into the assumptions of the previous statement.

However, the results of [25] assert that $k$-local density matrices of ground states of gapped Hamiltonians in $1D$ can be approximated by constant depth circuits, giving evidence that our framework should also extend to such states. And to get there, a technical obstacle has to be overcome in the proof of Thm. C.1. Essentially we need to show that local reduced density matrices of the ground state are already well-approximated at inverse temperature $\log(\epsilon^{-1})$. With such a statement we could show that we can still approximate the expectation values of $E_i$ at inverse temperature $\log(\epsilon^{-1})$ from samples from the ground state $|\psi\rangle$.

# D Summary of known strong convexity constants

As we see in the statement of Thm. A.1 and C.1, having a bound on the strong convexity constant $L^{-1}$ can give a quadratic improvement in the sample complexity in terms of the error $\epsilon$. Here we will briefly summarize for which cases estimates on $L$ are known in the literature for the classes of states we considered here.

**General many-body quantum:** first, we should mention the results of [5]. There the authors show bounds on $L^{-1}$ for arbitrary many-body Hamiltonians and temperatures that scales linearly in $m$. Thus, although certainly nontrivial, these bounds do not improve the sample complexity for the regimes we are interested in this work, namely that of logarithmic sample complexity in system size. To obtain improvements in this regime, $L^{-1}$ needs to be at most polylogarithmic in system size.

In the case of high-temperature Gibbs states, the recent work of [34] shows that this is indeed the case. I.e., in their Corollary 4.4 they show that for $\beta = O(k^{-8})$, where $k$ is the locality of the Hamiltonian, we indeed have $L^{-1} = O(\beta^{-2})$. It should also be noted that their results do not hold only for geometrically local Hamiltonians, but rather any Hamiltonian such that each term acts on at most $k$ qubits. This implies that for the high temperature regime, for which we also have the TC inequality in Thm. B.1, the improved sample complexity yielded by our methods holds. Note, however, that there is a slight mismatch between the inverse temperature range for which the two results hold: for the strong convexity we need $\beta = O(k^{-8})$, whereas for TC $\beta = O(k^{-1})$ suffices.

**Commuting Hamiltonians:** as we will prove later in Prop. F.1, in the case of commuting Hamiltonians we have that $L^{-1} = O(\epsilon^\beta \beta^{-2})$. Thus, for any constant inverse temperature $\beta > 0$ we have an improved sample complexity. However, in order to analyse ground states in $1D$, our current proof techniques still require an inverse temperature scaling logarithmically in system size, so for such states we do not obtain improvements through strong convexity. We plan to address this gap in future work.

Besides the cases mentioned above, we also considered the case of local circuits in this work. For those there are no nontrivial estimates on $L$ available to the best of our knowledge.

# E Regimes of efficient postprocessing

The only question we have still to answer is how to perform the postprocessing efficiently, namely how the parameter $C_\epsilon$ appearing in Theorem C.1 scales and how we obtain the bounds in Table 2.

There have been many recent breakthroughs on approximating quantum Gibbs states efficiently on classical computers [35, 46, 49, 51, 52]. The gradient descent method only requires us to estimate the gradient of the
This shows that shallow circuits can be learned efficiently as long as $l_0 = O(\log(n))$. However, it is not immediately obvious how to also estimate the expectation values of the Gibbs states required to run the maximum entropy algorithm. Thus, at least with the methods presented here, the postprocessing takes exponential time in the number of qubits.

**Ground states of gapped systems:** in light of our results for shallow circuits, it is natural to ask to what extent our framework can be extended to ground states of gapped Hamiltonians, specially in 1D. Thus, let us briefly comment on the technical barriers on the way of such statements. First, notice that the statement of Thm C.1 required inverse temperatures scaling logarithmically in system size for the approximation of the ground state. Most known TC inequalities have an exponential scaling with the inverse temperature and, thus, TC at this inverse temperature the savings compared to Pinsker are not quadratic, hindering a straightforward extension to gapped systems. There are some nontrivial examples of ground states satisfying a TC inequality with the constant depending inverse linearly with the temperature, like graph states [63]. But, as they can also be prepared from a constant depth quantum circuit almost by definition, they fall into the assumptions of the previous statement.

However, the results of [25] assert that $k$-local density matrices of ground states of gapped Hamiltonians in $1D$ can be approximated by constant depth circuits, giving evidence that our framework should also extend to such states. And to get there, a technical obstacle has to be overcome in the proof of Thm. C.1. Essentially we need to show that local reduced density matrices of the ground state are already well-approximated at inverse temperature $\log(\epsilon^{-1})$. With such a statement we could show that we can still approximate the expectation values of $E_i$ at inverse temperature $\log(\epsilon^{-1})$ from samples from the ground state $|\psi\rangle$.

# D Summary of known strong convexity constants

As we see in the statement of Thm. A.1 and C.1, having a bound on the strong convexity constant $L^{-1}$ can give a quadratic improvement in the sample complexity in terms of the error $\epsilon$. Here we will briefly summarize for which cases estimates on $L$ are known in the literature for the classes of states we considered here.

**General many-body quantum:** first, we should mention the results of [5]. There the authors show bounds on $L^{-1}$ for arbitrary many-body Hamiltonians and temperatures that scales linearly in $m$. Thus, although certainly nontrivial, these bounds do not improve the sample complexity for the regimes we are interested in this work, namely that of logarithmic sample complexity in system size. To obtain improvements in this regime, $L^{-1}$ needs to be at most polylogarithmic in system size.

In the case of high-temperature Gibbs states, the recent work of [34] shows that this is indeed the case. I.e., in their Corollary 4.4 they show that for $\beta = O(k^{-8})$, where $k$ is the locality of the Hamiltonian, we indeed have $L^{-1} = O(\beta^{-2})$. It should also be noted that their results do not hold only for geometrically local Hamiltonians, but rather any Hamiltonian such that each term acts on at most $k$ qubits. This implies that for the high temperature regime, for which we also have the TC inequality in Thm. B.1, the improved sample complexity yielded by our methods holds. Note, however, that there is a slight mismatch between the inverse temperature range for which the two results hold: for the strong convexity we need $\beta = O(k^{-8})$, whereas for TC $\beta = O(k^{-1})$ suffices.

**Commuting Hamiltonians:** as we will prove later in Prop. F.1, in the case of commuting Hamiltonians we have that $L^{-1} = O(\epsilon^\beta \beta^{-2})$. Thus, for any constant inverse temperature $\beta > 0$ we have an improved sample complexity. However, in order to analyse ground states in $1D$, our current proof techniques still require an inverse temperature scaling logarithmically in system size, so for such states we do not obtain improvements through strong convexity. We plan to address this gap in future work.

Besides the cases mentioned above, we also considered the case of local circuits in this work. For those there are no nontrivial estimates on $L$ available to the best of our knowledge.

# E Regimes of efficient postprocessing

The only question we have still to answer is how to perform the postprocessing efficiently, namely how the parameter $C_\epsilon$ appearing in Theorem C.1 scales and how we obtain the bounds in Table 2.

There have been many recent breakthroughs on approximating quantum Gibbs states efficiently on classical computers [35, 46, 49, 51, 52]. 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 $Z(\mu)$ efficiently or approximate $e(\lambda)$ suffices for our purposes.

For Gibbs states on a lattice, the methods of [51] yield that we can perform such approximations on a classical computer 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 Thm. B.1 hold, $C_\epsilon$ is polynomial in system size and we can also obtain efficient classical postprocessing.

For the case of Gibbs states of 1D systems, to the best of our knowledge, the best results available are those of [52]. 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 [51, 52] 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 commuting Hamiltonians.

Finally, for local quantum circuits we are not aware at this stage of any method that could yield a better postprocessing complexity than computing the partition function explicitly. This would yield a postprocessing that is exponential in the system size, as it requires diagonalizing an exponentially large matrix.

## F A complete picture: commuting Gibbs states

In this section, we discuss two classes of states for which the strongest version of our theorems holds, namely that of commuting 1D Gibbs states, and the one of high-temperature commuting Gibbs states on arbitrary hypergraphs. We already discussed that they satisfy transportation cost inequalities in Thm. B.1. Thus, the last missing ingredient to obtain an optimal performance is to show that the partition function is indeed strongly convex. More precisely, we will now establish that, for these classes of states, both the upper and lower bounds on the Hessian of the log partition function are order 1. In addition to that, with access to a quantum computer, it is possible to perform the post-processing in time $\tilde{O}(m)$. As writing down the vector $\lambda$ takes $\Omega(m)$ time, we conclude that the postprocessing can be performed in a time comparable to writing down a solution to the problem. Thus, our procedure is essentially optimal.

Also in the setting of commuting Gibbs states, it is worth noting that after the completion of this work, we became aware of [4], which gives another method to learn the Gibbs state and its Hamiltonian that neither involves the maximum entropy method nor requires strong convexity guarantees. Their algorithm works by learning local reduced density matrices and showing that the parameters $\lambda$ of the Hamiltonian of a commuting Gibbs state can be efficiently estimated from that. In principle, obtaining a bound on $\lambda$ also suffices for our purposes and we could alternatively use their methods to bypass having to solve the maximum entropy problem for such states. In particular, this means that the postprocessing with their methods could be performed even for temperatures at which the partition function cannot be estimated efficiently but we still have access to samples from the state. However, as we ultimately are interested in the regime in which TC inequalities hold, which corresponds to the high-temperature regime, we do not further comment on their results.

In this section, we consider 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$. In what follows, we also denote by $S(v, r)$ the sphere of radius $r$ centered at vertex $v \in V$, and define for all $r \in \mathbb{N}$:

$$B(r) := \max_{v \in V} |B(v, r)|, \qquad S(r) := \max_{v \in V} |S(v, r)|.$$

Next, we consider a Gibbs state $\sigma(\lambda)$ on $\mathcal{H}_V := \bigotimes_{v \in V} \mathcal{H}_v$, where $\dim(\mathcal{H}_v) = d$ for all $v \in V$. More precisely, $\sigma(\lambda) := e^{-\beta H(\lambda)}/\mathcal{Z}(\lambda)$, where with a slight abuse of notations

$$H(\lambda) = \sum_{i=1}^{m} \lambda_i E_i = \sum_{A \in E} \sum_{\alpha_A \in [d^2]^{|A|}} \alpha_A E_{\alpha_A}, \qquad (55)$$

with $\|E_i\|_\infty = 1$ for all $i \in \{1, \dots, m\}$, $[E_i, E_j] = 0$ for all $i, j \in \{1, \dots, m\}$, and where the sets $\{1, \dots, m\}$ and $\{\alpha_A\}_{A \in E, \alpha_A \in [d^2]^A}$ are in bijection. Note that $B(r_0)$ bounds the maximal locality of the Hamiltonian.

We also denote by $\sigma_A(\lambda)$ the Gibbs state corresponding to the restriction of $H$ onto the region $A$, i.e.

$$\sigma_A(\lambda) := \frac{e^{-\beta H_A(\lambda)}}{\text{tr}[e^{-\beta H_A(\lambda)}]}, \quad \text{where} \quad H_A(\lambda) := \sum_{i \mid \text{supp}(E_i) \cap A \neq \emptyset} \lambda_i E_i.$$
partition function for a Gibbs state at each iteration. Thus, any algorithm that can approximate the log-partition function $Z(\mu)$ efficiently or approximate $e(\lambda)$ suffices for our purposes.

For Gibbs states on a lattice, the methods of [51] yield that we can perform such approximations on a classical computer 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 Thm. B.1 hold, $C_\epsilon$ is polynomial in system size and we can also obtain efficient classical postprocessing.

For the case of Gibbs states of 1D systems, to the best of our knowledge, the best results available are those of [52]. 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 [51, 52] 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 commuting Hamiltonians.

Finally, for local quantum circuits we are not aware at this stage of any method that could yield a better postprocessing complexity than computing the partition function explicitly. This would yield a postprocessing that is exponential in the system size, as it requires diagonalizing an exponentially large matrix.

## F A complete picture: commuting Gibbs states

In this section, we discuss two classes of states for which the strongest version of our theorems holds, namely that of commuting 1D Gibbs states, and the one of high-temperature commuting Gibbs states on arbitrary hypergraphs. We already discussed that they satisfy transportation cost inequalities in Thm. B.1. Thus, the last missing ingredient to obtain an optimal performance is to show that the partition function is indeed strongly convex. More precisely, we will now establish that, for these classes of states, both the upper and lower bounds on the Hessian of the log partition function are order 1. In addition to that, with access to a quantum computer, it is possible to perform the post-processing in time $\tilde{O}(m)$. As writing down the vector $\lambda$ takes $\Omega(m)$ time, we conclude that the postprocessing can be performed in a time comparable to writing down a solution to the problem. Thus, our procedure is essentially optimal.

Also in the setting of commuting Gibbs states, it is worth noting that after the completion of this work, we became aware of [4], which gives another method to learn the Gibbs state and its Hamiltonian that neither involves the maximum entropy method nor requires strong convexity guarantees. Their algorithm works by learning local reduced density matrices and showing that the parameters $\lambda$ of the Hamiltonian of a commuting Gibbs state can be efficiently estimated from that. In principle, obtaining a bound on $\lambda$ also suffices for our purposes and we could alternatively use their methods to bypass having to solve the maximum entropy problem for such states. In particular, this means that the postprocessing with their methods could be performed even for temperatures at which the partition function cannot be estimated efficiently but we still have access to samples from the state. However, as we ultimately are interested in the regime in which TC inequalities hold, which corresponds to the high-temperature regime, we do not further comment on their results.

In this section, we consider 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$. In what follows, we also denote by $S(v, r)$ the sphere of radius $r$ centered at vertex $v \in V$, and define for all $r \in \mathbb{N}$:

$$B(r) := \max_{v \in V} |B(v, r)|, \qquad S(r) := \max_{v \in V} |S(v, r)|.$$

Next, we consider a Gibbs state $\sigma(\lambda)$ on $\mathcal{H}_V := \bigotimes_{v \in V} \mathcal{H}_v$, where $\dim(\mathcal{H}_v) = d$ for all $v \in V$. More precisely, $\sigma(\lambda) := e^{-\beta H(\lambda)}/\mathcal{Z}(\lambda)$, where with a slight abuse of notations

$$H(\lambda) = \sum_{i=1}^{m} \lambda_i E_i = \sum_{A \in E} \sum_{\alpha_A \in [d^2]^{|A|}} \alpha_A E_{\alpha_A}, \qquad (55)$$

with $\|E_i\|_\infty = 1$ for all $i \in \{1, \dots, m\}$, $[E_i, E_j] = 0$ for all $i, j \in \{1, \dots, m\}$, and where the sets $\{1, \dots, m\}$ and $\{\alpha_A\}_{A \in E, \alpha_A \in [d^2]^A}$ are in bijection. Note that $B(r_0)$ bounds the maximal locality of the Hamiltonian.

We also denote by $\sigma_A(\lambda)$ the Gibbs state corresponding to the restriction of $H$ onto the region $A$, i.e.

$$\sigma_A(\lambda) := \frac{e^{-\beta H_A(\lambda)}}{\text{tr}[e^{-\beta H_A(\lambda)}]}, \quad \text{where} \quad H_A(\lambda) := \sum_{i \mid \text{supp}(E_i) \cap A \neq \emptyset} \lambda_i E_i.$$
Note that, in general, $\sigma_A(\lambda) \neq \text{tr}_{A^c}(\sigma(\lambda))$.

## F.1 Upper bound on Hessian for commuting Gibbs states with decay of correlations

In this section, we prove tightened strong convexity constants for the log partition function in the case when the Gibbs state arises from a local commuting Hamiltonian at high-temperature. In fact, the upper constant found in Proposition A.4 can be tightened into a system size-independent one under the condition of exponential decay of the correlations in the Gibbs state. Several notions of exponential decay of correlations exist in the literature [18]. Here, we will say that a Gibbs state $\sigma$ has correlation length $\xi$ if for all observables $O_A, O_B$ supported on non-overlapping regions A and B respectively, we have that:

$$|\text{tr}[\sigma O_A \otimes O_B] - \text{tr}[\sigma O_A] \text{tr}[\sigma O_B]| \le c \|O_A\|_{\infty} \|O_B\|_{\infty} e^{-\xi \text{dist}(A,B)}$$

for some constant $c > 0$. In the classical setting, this condition is known to hold at any inverse temperature for 1D systems, and below a critical inverse temperature $\beta_c$ that depends on the locality of the Hamiltonian when $D \ge 2$ [28]. The same bound also holds in the quantum case for 1D systems [6], or above a critical temperature on regular lattices when $D \ge 2$ [35, 46]. Using these bounds, we obtain the following improvement of Proposition A.4 which shows that for this class of states $U = \mathcal{O}(1)$.

**Lemma F.1** (Strengthened upper strong convexity at high-temperature and for 1D systems). For each $\mu \in B_{\ell_{\infty}}(0, 1)$, let $\sigma(\mu)$ be a Gibbs state at inverse temperature $\beta < \beta_c$ corresponding to the Hamiltonian defined on the hypergraph $G = (V, E)$ in (55). Then

$$\nabla^2 \log(\mathcal{Z}(\mu)) \le (c \beta^2 B(r_0) B(2r_0) d^{2B(r_0)} \sum_{r=1}^{\infty} e^{-\xi r} S(r)) I,$$

where $\xi$ is the correlation length of the state. Moreover, this result holds for all $\beta > 0$ in 1D.

*Proof.* Let us first use the assumption of commutativity to simplify the expression for the Hessian. We find for all $\alpha_A \in [d^2]^{|A|}, \alpha'_B \in [d^2]^{|B|}, A, B \in E$,

$$(\nabla^2 \log(\mathcal{Z}(\mu)))_{\alpha_A \alpha'_B} = \beta^2 \text{tr} [\sigma(\mu) (E_{\alpha_A} - \text{tr}[\sigma(\mu) E_{\alpha_A}]) (E_{\alpha'_B} - \text{tr}[\sigma(\mu) E_{\alpha'_B}])] . \quad (56)$$

The rest follows similarly to Proposition A.4 from Gershgorin's circle theorem together with the decay of correlations arising at $\beta < \beta_c$: for all $\alpha_A$,

$$\sum_{\alpha'_B \neq \alpha_A} |(\nabla^2 \log(\mathcal{Z}(\mu)))_{\alpha_A \alpha'_B}| \le c \beta^2 \sum_{\alpha'_B \neq \alpha_A} e^{-\xi \text{dist}(A,B)},$$

where we also used that the basis operators $\{E_i\}$ have operator norm 1. The claim then follows by bounding the number of basis operators whose support is at a distance $r$ of A for each $r \in \mathbb{N}$: the latter is bounded by the product of (i) the number of vertices in A; (ii) the number of vertices at a distance $r$ of a given vertex; (iii) the number of hyperedges containing a given vertex; and (iv) the number of basis operators corresponding to a given hyperedge. A simple estimate of each of these quantities gives the bound $B(r_0) S(r) B(2r_0) d^{2B(r_0)}$. Therefore:

$$\sum_{\alpha'_B \neq \alpha_A} |(\nabla^2 \log(\mathcal{Z}(\mu)))_{\alpha_A \alpha'_B}| \le c \beta^2 B(r_0) B(2r_0) d^{2B(r_0)} \sum_{r=1}^{\infty} e^{-\xi r} S(r).$$

Note that for D-regular graphs and $r_0 = \mathcal{O}(1)$ we have $B(r_0) B(2r_0) d^{2B(r_0)} = \mathcal{O}(1)$ and $S(r) = \mathcal{O}(r^{D-1})$, giving a scaling of $\nabla^2 \log(\mathcal{Z}(\mu)) = \mathcal{O}(\beta^2 \xi^D)$.

## F.2 Lower bound on Hessian for commuting Gibbs states

Whenever the Gibbs state is assumed to be commuting, the lower strong convexity constant $L$ can also be made independent of $m$, by a direct generalization of the classical argument as found in [66][Lemma 7] or [54] (see also [5]).
Note that, in general, $\sigma_A(\lambda) \neq \text{tr}_{A^c}(\sigma(\lambda))$.

## F.1 Upper bound on Hessian for commuting Gibbs states with decay of correlations

In this section, we prove tightened strong convexity constants for the log partition function in the case when the Gibbs state arises from a local commuting Hamiltonian at high-temperature. In fact, the upper constant found in Proposition A.4 can be tightened into a system size-independent one under the condition of exponential decay of the correlations in the Gibbs state. Several notions of exponential decay of correlations exist in the literature [18]. Here, we will say that a Gibbs state $\sigma$ has correlation length $\xi$ if for all observables $O_A, O_B$ supported on non-overlapping regions A and B respectively, we have that:

$$|\text{tr}[\sigma O_A \otimes O_B] - \text{tr}[\sigma O_A] \text{tr}[\sigma O_B]| \le c \|O_A\|_{\infty} \|O_B\|_{\infty} e^{-\xi \text{dist}(A,B)}$$

for some constant $c > 0$. In the classical setting, this condition is known to hold at any inverse temperature for 1D systems, and below a critical inverse temperature $\beta_c$ that depends on the locality of the Hamiltonian when $D \ge 2$ [28]. The same bound also holds in the quantum case for 1D systems [6], or above a critical temperature on regular lattices when $D \ge 2$ [35, 46]. Using these bounds, we obtain the following improvement of Proposition A.4 which shows that for this class of states $U = \mathcal{O}(1)$.

**Lemma F.1** (Strengthened upper strong convexity at high-temperature and for 1D systems). For each $\mu \in B_{\ell_{\infty}}(0, 1)$, let $\sigma(\mu)$ be a Gibbs state at inverse temperature $\beta < \beta_c$ corresponding to the Hamiltonian defined on the hypergraph $G = (V, E)$ in (55). Then

$$\nabla^2 \log(\mathcal{Z}(\mu)) \le (c \beta^2 B(r_0) B(2r_0) d^{2B(r_0)} \sum_{r=1}^{\infty} e^{-\xi r} S(r)) I,$$

where $\xi$ is the correlation length of the state. Moreover, this result holds for all $\beta > 0$ in 1D.

*Proof.* Let us first use the assumption of commutativity to simplify the expression for the Hessian. We find for all $\alpha_A \in [d^2]^{|A|}, \alpha'_B \in [d^2]^{|B|}, A, B \in E$,

$$(\nabla^2 \log(\mathcal{Z}(\mu)))_{\alpha_A \alpha'_B} = \beta^2 \text{tr} [\sigma(\mu) (E_{\alpha_A} - \text{tr}[\sigma(\mu) E_{\alpha_A}]) (E_{\alpha'_B} - \text{tr}[\sigma(\mu) E_{\alpha'_B}])] . \quad (56)$$

The rest follows similarly to Proposition A.4 from Gershgorin's circle theorem together with the decay of correlations arising at $\beta < \beta_c$: for all $\alpha_A$,

$$\sum_{\alpha'_B \neq \alpha_A} |(\nabla^2 \log(\mathcal{Z}(\mu)))_{\alpha_A \alpha'_B}| \le c \beta^2 \sum_{\alpha'_B \neq \alpha_A} e^{-\xi \text{dist}(A,B)},$$

where we also used that the basis operators $\{E_i\}$ have operator norm 1. The claim then follows by bounding the number of basis operators whose support is at a distance $r$ of A for each $r \in \mathbb{N}$: the latter is bounded by the product of (i) the number of vertices in A; (ii) the number of vertices at a distance $r$ of a given vertex; (iii) the number of hyperedges containing a given vertex; and (iv) the number of basis operators corresponding to a given hyperedge. A simple estimate of each of these quantities gives the bound $B(r_0) S(r) B(2r_0) d^{2B(r_0)}$. Therefore:

$$\sum_{\alpha'_B \neq \alpha_A} |(\nabla^2 \log(\mathcal{Z}(\mu)))_{\alpha_A \alpha'_B}| \le c \beta^2 B(r_0) B(2r_0) d^{2B(r_0)} \sum_{r=1}^{\infty} e^{-\xi r} S(r).$$

Note that for D-regular graphs and $r_0 = \mathcal{O}(1)$ we have $B(r_0) B(2r_0) d^{2B(r_0)} = \mathcal{O}(1)$ and $S(r) = \mathcal{O}(r^{D-1})$, giving a scaling of $\nabla^2 \log(\mathcal{Z}(\mu)) = \mathcal{O}(\beta^2 \xi^D)$.

## F.2 Lower bound on Hessian for commuting Gibbs states

Whenever the Gibbs state is assumed to be commuting, the lower strong convexity constant $L$ can also be made independent of $m$, by a direct generalization of the classical argument as found in [66][Lemma 7] or [54] (see also [5]).
Before we state and prove our result, let us introduce a few useful concepts: given a full-rank quantum state $\sigma$, we denote the weighted 2-norm of an observable Y as

$$ \|Y\|_{2,\sigma} := \mathrm{tr} \left[ \left| \sigma^{\frac{1}{4}} Y \sigma^{\frac{1}{4}} \right|^2 \right]^{1/2} . $$

and refer to the corresponding non-commutative weighted $L_2$ space with inner product $\langle X, Y \rangle_{\sigma} := \mathrm{tr} [\sigma^{\frac{1}{2}} X^{\dagger} \sigma^{\frac{1}{2}} Y]$ as $L_2(\sigma)$. The Petz recovery map corresponding to a subregion $A \subset V$ with respect to $\sigma(\mu)$ is defined as

$$ \mathcal{R}_{A,\sigma(\mu)}(X) = \mathrm{tr}_A(\sigma(\mu))^{-1/2} \mathrm{tr}_A(\sigma(\mu)^{1/2} X \sigma(\mu)^{1/2}) \mathrm{tr}_A(\sigma(\mu))^{-1/2}, $$

We will also need the notion of a conditional expectation $\mathbb{E}_A$ with respect to the state $\sigma(\mu)$ into the region $A \subset \Lambda$ (see [7, 44] for more details). For instance, one can choose $\mathbb{E}_A := \lim_{n \to \infty} \mathcal{R}_{A,\sigma(\mu)}^n$, where $\mathcal{R}_{A,\sigma(\mu)}$ is the Petz recovery map of $\sigma(\mu)$. In other words, the map $\mathbb{E}_A$ is a completely positive, unital map that projects onto the algebra $\mathcal{N}_A$ of fixed points of $\mathcal{R}_{A,\sigma(\mu)}$. This algebra is known to be expressed as the commutant [7]

$$ \mathcal{N}_A := \{\sigma(\mu)^{it} \mathcal{B}(\mathcal{H}_A) \sigma(\mu)^{-it}; t \in \mathbb{R}\}'. $$

Moreover, the maps $\mathbb{E}_A$ commute with the modular operator $\Delta_{\sigma(\mu)}(.) := \sigma(\mu)(.)\sigma(\mu)^{-1}$, and for any $X_A, Z_A \in \mathcal{N}_A$ and all $Y \in \mathcal{B}(\mathcal{H}_V)$, $\mathbb{E}_A[X_A Y Z_A] = X_A \mathbb{E}_A[Y] Z_A$. The commutativity condition for $H(\mu)$ implies frustration freeness of the family of conditional expectations $\{\mathbb{E}_A\}_{A \in E}$: for any $X \in L_2(\sigma(\mu))$, $\|\mathbb{E}_A(X)\|_{2,\sigma(\mu)}^2 + \|(\mathrm{id} - \mathbb{E}_A)(X)\|_{2,\sigma(\mu)}^2 = \|X\|_{L_2(\sigma(\mu))}^2$.

The next technical lemma is essential in the derivation of our strong convexity lower bound. With a slight abuse of notations, we use the simplified notations $\sigma_x(\lambda) := \sigma_{\{x\}}(\lambda)$, $H_x(\lambda) := H_{\{x\}}(\lambda)$ and so on.

**Lemma F.2.** Let $H = \sum_j \mu_j E_j$ be a local commuting Hamiltonian on the hypergraph $G = (V, E)$ defined in (55), each local operator $E_j$ is further assumed to be traceless. The following holds for any $x \in V$:

$$ c(x, \beta) := \max_{\mu \in B_{\ell_{\infty}}(0,1)} \frac{\| \mathcal{R}_{x, \sigma_x(\mu)}(H_x) - \mathrm{tr}[\sigma_x(\mu)H_x] I \|_{2, \sigma_x(\mu)}}{\| H_x - \mathrm{tr}[\sigma_x(\mu)H_x] I \|_{2, \sigma_x(\mu)}} < 1. $$

*Proof.* We first prove that $X = \mathcal{R}_{x,\sigma_x(\mu)}(X)$ is equivalent to $\|\mathcal{R}_{x,\sigma_x(\mu)}(X)\|_{2,\sigma_x(\mu)} = \|X\|_{2,\sigma_x(\mu)}$. One direction trivially holds. For the opposite direction, assume that $X \neq \mathcal{R}_{x,\sigma_x(\mu)}(X)$. This means that $X = Y + Z$, with $Y, Z$ two operators that are orthogonal in $L_2(\sigma_x(\mu))$, with $\mathcal{R}_{x,\sigma_x(\mu)}(Y) = Y$ and $Z \neq 0$. Now, since $\mathcal{R}_{x,\sigma_x(\mu)}$ is self-adjoint and unital, it strictly contracts elements in the orthogonal of its fixed points and we have

$$ \begin{aligned} \| \mathcal{R}_{x, \sigma_x(\mu)}(X) \|_{2, \sigma_x(\mu)}^2 &= \|Y\|_{2, \sigma_x(\mu)}^2 + \| \mathcal{R}_{x, \sigma_x(\mu)}(Z) \|_{2, \sigma_x(\mu)}^2 \\ &< \|Y\|_{2, \sigma_x(\mu)}^2 + \|Z\|_{2, \sigma_x(\mu)}^2 \\ &= \|X\|_{2, \sigma_x(\mu)}^2, \end{aligned} $$

which contradicts the condition of equality of the norms. Now, since the map $\mathcal{R}_{x,\sigma_x(\mu)}$ is unital, it suffices to prove that $\mathcal{R}_{x,\sigma_x(\mu)}(H_x) \neq H_x$, or equivalently $\mathbb{E}_x(H_x) \neq H_x$, in order to conclude. Let us assume instead that equality holds. This means that, for all observables $A_x$ supported on site $x$, and all $t \in \mathbb{R}$:

$$ [H_x, \sigma(\mu)^{it} A_x \sigma(\mu)^{-it}] = \sigma(\mu)^{it} [H_x, A_x] \sigma(\mu)^{-it} = 0 \quad \Rightarrow \quad H_x = \frac{I_x}{d} \otimes \mathrm{tr}_x(H_x). $$

However, this contradicts the fact that $H_x$ is traceless on site $x$. Therefore $\mathcal{R}_{x,\sigma_x(\mu)}(H_x) \neq H_x$ and the proof follows. □

We are ready to prove our strong convexity lower bound.

**Proposition F.1** (Strengthened lower 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)$ defined in (55), where $\mathrm{tr}[E_i E_j] = 0$ for all $i \neq j$ and each local operator $E_j$ is traceless on its support. Then the Hessian of the log-partition function is lower bounded by

$$ \nabla^2 \log \mathcal{Z}(\mu) \geq \beta^2 e^{-\beta(B(2r_0)+2B(4r_0))} d^{-B(2r_0)} (1 - c(\beta)^2) I, $$

where $c(\beta) := \max_{v \in V} c(v, \beta)$.
Before we state and prove our result, let us introduce a few useful concepts: given a full-rank quantum state $\sigma$, we denote the weighted 2-norm of an observable Y as

$$ \|Y\|_{2,\sigma} := \mathrm{tr} \left[ \left| \sigma^{\frac{1}{4}} Y \sigma^{\frac{1}{4}} \right|^2 \right]^{1/2} . $$

and refer to the corresponding non-commutative weighted $L_2$ space with inner product $\langle X, Y \rangle_{\sigma} := \mathrm{tr} [\sigma^{\frac{1}{2}} X^{\dagger} \sigma^{\frac{1}{2}} Y]$ as $L_2(\sigma)$. The Petz recovery map corresponding to a subregion $A \subset V$ with respect to $\sigma(\mu)$ is defined as

$$ \mathcal{R}_{A,\sigma(\mu)}(X) = \mathrm{tr}_A(\sigma(\mu))^{-1/2} \mathrm{tr}_A(\sigma(\mu)^{1/2} X \sigma(\mu)^{1/2}) \mathrm{tr}_A(\sigma(\mu))^{-1/2}, $$

We will also need the notion of a conditional expectation $\mathbb{E}_A$ with respect to the state $\sigma(\mu)$ into the region $A \subset \Lambda$ (see [7, 44] for more details). For instance, one can choose $\mathbb{E}_A := \lim_{n \to \infty} \mathcal{R}_{A,\sigma(\mu)}^n$, where $\mathcal{R}_{A,\sigma(\mu)}$ is the Petz recovery map of $\sigma(\mu)$. In other words, the map $\mathbb{E}_A$ is a completely positive, unital map that projects onto the algebra $\mathcal{N}_A$ of fixed points of $\mathcal{R}_{A,\sigma(\mu)}$. This algebra is known to be expressed as the commutant [7]

$$ \mathcal{N}_A := \{\sigma(\mu)^{it} \mathcal{B}(\mathcal{H}_A) \sigma(\mu)^{-it}; t \in \mathbb{R}\}'. $$

Moreover, the maps $\mathbb{E}_A$ commute with the modular operator $\Delta_{\sigma(\mu)}(.) := \sigma(\mu)(.)\sigma(\mu)^{-1}$, and for any $X_A, Z_A \in \mathcal{N}_A$ and all $Y \in \mathcal{B}(\mathcal{H}_V)$, $\mathbb{E}_A[X_A Y Z_A] = X_A \mathbb{E}_A[Y] Z_A$. The commutativity condition for $H(\mu)$ implies frustration freeness of the family of conditional expectations $\{\mathbb{E}_A\}_{A \in E}$: for any $X \in L_2(\sigma(\mu))$, $\|\mathbb{E}_A(X)\|_{2,\sigma(\mu)}^2 + \|(\mathrm{id} - \mathbb{E}_A)(X)\|_{2,\sigma(\mu)}^2 = \|X\|_{L_2(\sigma(\mu))}^2$.

The next technical lemma is essential in the derivation of our strong convexity lower bound. With a slight abuse of notations, we use the simplified notations $\sigma_x(\lambda) := \sigma_{\{x\}}(\lambda)$, $H_x(\lambda) := H_{\{x\}}(\lambda)$ and so on.

**Lemma F.2.** Let $H = \sum_j \mu_j E_j$ be a local commuting Hamiltonian on the hypergraph $G = (V, E)$ defined in (55), each local operator $E_j$ is further assumed to be traceless. The following holds for any $x \in V$:

$$ c(x, \beta) := \max_{\mu \in B_{\ell_{\infty}}(0,1)} \frac{\| \mathcal{R}_{x, \sigma_x(\mu)}(H_x) - \mathrm{tr}[\sigma_x(\mu)H_x] I \|_{2, \sigma_x(\mu)}}{\| H_x - \mathrm{tr}[\sigma_x(\mu)H_x] I \|_{2, \sigma_x(\mu)}} < 1. $$

*Proof.* We first prove that $X = \mathcal{R}_{x,\sigma_x(\mu)}(X)$ is equivalent to $\|\mathcal{R}_{x,\sigma_x(\mu)}(X)\|_{2,\sigma_x(\mu)} = \|X\|_{2,\sigma_x(\mu)}$. One direction trivially holds. For the opposite direction, assume that $X \neq \mathcal{R}_{x,\sigma_x(\mu)}(X)$. This means that $X = Y + Z$, with $Y, Z$ two operators that are orthogonal in $L_2(\sigma_x(\mu))$, with $\mathcal{R}_{x,\sigma_x(\mu)}(Y) = Y$ and $Z \neq 0$. Now, since $\mathcal{R}_{x,\sigma_x(\mu)}$ is self-adjoint and unital, it strictly contracts elements in the orthogonal of its fixed points and we have

$$ \begin{aligned} \| \mathcal{R}_{x, \sigma_x(\mu)}(X) \|_{2, \sigma_x(\mu)}^2 &= \|Y\|_{2, \sigma_x(\mu)}^2 + \| \mathcal{R}_{x, \sigma_x(\mu)}(Z) \|_{2, \sigma_x(\mu)}^2 \\ &< \|Y\|_{2, \sigma_x(\mu)}^2 + \|Z\|_{2, \sigma_x(\mu)}^2 \\ &= \|X\|_{2, \sigma_x(\mu)}^2, \end{aligned} $$

which contradicts the condition of equality of the norms. Now, since the map $\mathcal{R}_{x,\sigma_x(\mu)}$ is unital, it suffices to prove that $\mathcal{R}_{x,\sigma_x(\mu)}(H_x) \neq H_x$, or equivalently $\mathbb{E}_x(H_x) \neq H_x$, in order to conclude. Let us assume instead that equality holds. This means that, for all observables $A_x$ supported on site $x$, and all $t \in \mathbb{R}$:

$$ [H_x, \sigma(\mu)^{it} A_x \sigma(\mu)^{-it}] = \sigma(\mu)^{it} [H_x, A_x] \sigma(\mu)^{-it} = 0 \quad \Rightarrow \quad H_x = \frac{I_x}{d} \otimes \mathrm{tr}_x(H_x). $$

However, this contradicts the fact that $H_x$ is traceless on site $x$. Therefore $\mathcal{R}_{x,\sigma_x(\mu)}(H_x) \neq H_x$ and the proof follows. $\square$

We are ready to prove our strong convexity lower bound.

**Proposition F.1** (Strengthened lower 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)$ defined in (55), where $\mathrm{tr}[E_i E_j] = 0$ for all $i \neq j$ and each local operator $E_j$ is traceless on its support. Then the Hessian of the log-partition function is lower bounded by

$$ \nabla^2 \log \mathcal{Z}(\mu) \geq \beta^2 e^{-\beta(B(2r_0)+2B(4r_0))} d^{-B(2r_0)} (1 - c(\beta)^2) I, $$

where $c(\beta) := \max_{v \in V} c(v, \beta)$.
Proof. We first use the assumption of commutativity in order to simplify the expression for the Hessian. As in (56), we find

$$ (\nabla^2 \log(\mathcal{Z}(\mu)))_{ij} = \beta^2 \operatorname{tr} [\sigma(\mu) (E_i - \operatorname{tr}[\sigma(\mu)E_i]) (E_j - \operatorname{tr}[\sigma(\mu)E_j])] . $$

Therefore, for any linear combination $H \equiv H(\lambda) = \sum_j \lambda_j E_j$ of the basis vectors, we have

$$ \sum_{ij} \lambda_i \lambda_j (\nabla^2 \log(\mathcal{Z}(\mu)))_{ij} = \beta^2 \operatorname{Var}_{\sigma(\mu)}(H), $$

with $\operatorname{Var}_{\sigma(\mu)}(X) := \|X - \operatorname{tr}[\sigma(\mu)X]\|_{2,\sigma(\mu)}^2$. It is thus sufficient to lower-bound the latter. For this, we choose a subregion $A \subseteq V$ such that any basis element $E_i$ has support intersecting a unique vertex in $A$. We lower-bound the variance by

$$ \|H - \operatorname{tr}[\sigma(\mu)H] I\|_{2,\sigma(\mu)}^2 \ge \|(\operatorname{id} - \mathbb{E}_A)(H - \operatorname{tr}[\sigma(\mu)H] I)\|_{2,\sigma(\mu)}^2 = \|H - \mathbb{E}_A[H]\|_{2,\sigma(\mu)}^2, \quad (57) $$

where the first inequality follows by the $L_2(\sigma(\mu))$ contractivity of $\operatorname{id} - \mathbb{E}_A$. Now, the weighted norm can be further simplified into a sum of local weighted norms as follows: first, for any two $E_i, E_j$ whose supports intersect with a different vertex of $A$, we have

$$ \begin{aligned} \langle \mathbb{E}_A[E_i], \mathbb{E}_A[E_j] \rangle_{\sigma(\mu)} &= \operatorname{tr} [\sigma(\mu)^{1/2} \mathbb{E}_A[E_i] \sigma(\mu)^{1/2} \mathbb{E}_A[E_j]] \\ &= \operatorname{tr} [\sigma(\mu) \Delta_{\sigma(\mu)}^{-1/2} \circ \mathbb{E}_A[E_i] \mathbb{E}_A[E_j]] \\ &= \operatorname{tr} [\sigma(\mu) \mathbb{E}_A[E_i] \mathbb{E}_A[E_j]] \end{aligned} \quad (58) $$

where in the third line we used the commutation of the modular operator $\Delta_{\sigma(\mu)}(X) := \sigma(\mu)X\sigma(\mu)^{-1}$ with $\mathbb{E}_A$ together with the commutativity of $\sigma(\mu)$ and $E_i$. Now, denoting $\operatorname{supp}(E_i) \cap A = \{x\}$ and $\operatorname{supp}(E_j) \cap A = \{y\}$, we show that

$$ \mathbb{E}_A[E_i] = \mathbb{E}_x[E_i] \quad \text{and} \quad \mathbb{E}_A[E_j] = \mathbb{E}_y[E_j]. \quad (59) $$

In order to prove these two identities, we simply need to prove for instance that $\mathbb{E}_x[E_i]$ belongs to the image algebra $\mathcal{N}_A$ of $\mathbb{E}_A$ since $\mathcal{N}_A \subseteq \mathcal{N}_x$ by definition. Hence, it is enough to show that $\mathbb{E}_x[E_i]$ commutes with operators of the form $\sigma(\mu)^{it} X_A \sigma(\mu)^{-it}$ for any $t \in \mathbb{R}$ and any $X_A \in \mathcal{B}(\mathcal{H}_A)$. This claim follows from

$$ \begin{aligned} \mathbb{E}_x[E_i]\sigma(\mu)^{it}X_A\sigma(\mu)^{-it} &= \sigma(\mu)^{it}\sigma(\mu)^{-it}\mathbb{E}_x[E_i]\sigma(\mu)^{it}X_A\sigma(\mu)^{-it} \\ &= \sigma(\mu)^{it}\mathbb{E}_x[\sigma(\mu)^{-it}E_i\sigma(\mu)^{it}]X_A\sigma(\mu)^{-it} \\ &= \sigma(\mu)^{it}\mathbb{E}_x[E_i]X_A\sigma(\mu)^{-it} \\ &= \sigma(\mu)^{it}X_A\mathbb{E}_x[E_i]\sigma(\mu)^{-it} \\ &= \sigma(\mu)^{it}X_A\sigma(\mu)^{-it}\mathbb{E}_x[E_i]. \end{aligned} $$

where the fourth line follows from the fact that the support of $\mathbb{E}_x[E_i]$ does not intersect $A \setminus \{i\}$, together with the fact that $\mathbb{E}_x[E_i]$ is locally proportional to $I$ on site $x$, by definition of $\mathcal{N}_x$. Therefore, using (59) into (58), we get

$$ \begin{aligned} \langle \mathbb{E}_A[E_i], \mathbb{E}_A[E_j] \rangle_{\sigma(\mu)} &= \operatorname{tr} [\sigma(\mu)\mathbb{E}_x[E_i]\mathbb{E}_y[E_j]] \\ &= \operatorname{tr} [\sigma(\mu)\mathbb{E}_x\mathbb{E}_y[E_iE_j]] \\ &= \operatorname{tr} [\sigma(\mu)E_iE_j] \\ &= \langle E_i, E_j \rangle_{\sigma(\mu)}, \end{aligned} $$

where in the second line we used that $\mathbb{E}_x[E_i]$ is a fixed point of $\mathbb{E}_y$, and then that $E_j$ is a fixed point of $\mathbb{E}_x$, by the support conditions of $E_i$ and $E_j$. Therefore, the variance on the right-hand side of (57) can be simplified as

$$ \begin{aligned} \|H - \mathbb{E}_A[H]\|_{2,\sigma(\mu)}^2 &= \sum_{x \in A} \left\| \sum_{j | \operatorname{supp}(E_j) \ni x} \lambda_j (E_j - \mathbb{E}_x[E_j]) \right\|_{2,\sigma(\mu)}^2 \\ &= \sum_{x \in A} \left\| H_x - \mathbb{E}_x[H_x] \right\|_{2,\sigma(\mu)}^2 \end{aligned} $$
Proof. We first use the assumption of commutativity in order to simplify the expression for the Hessian. As in (56), we find

$$ (\nabla^2 \log(\mathcal{Z}(\mu)))_{ij} = \beta^2 \operatorname{tr} [\sigma(\mu) (E_i - \operatorname{tr}[\sigma(\mu)E_i]) (E_j - \operatorname{tr}[\sigma(\mu)E_j])] . $$

Therefore, for any linear combination $H \equiv H(\lambda) = \sum_j \lambda_j E_j$ of the basis vectors, we have

$$ \sum_{ij} \lambda_i \lambda_j (\nabla^2 \log(\mathcal{Z}(\mu)))_{ij} = \beta^2 \operatorname{Var}_{\sigma(\mu)}(H), $$

with $\operatorname{Var}_{\sigma(\mu)}(X) := \|X - \operatorname{tr}[\sigma(\mu)X]\|_{2,\sigma(\mu)}^2$. It is thus sufficient to lower-bound the latter. For this, we choose a subregion $A \subseteq V$ such that any basis element $E_i$ has support intersecting a unique vertex in $A$. We lower-bound the variance by

$$ \|H - \operatorname{tr}[\sigma(\mu)H] I\|_{2,\sigma(\mu)}^2 \ge \|(\operatorname{id} - \mathbb{E}_A)(H - \operatorname{tr}[\sigma(\mu)H] I)\|_{2,\sigma(\mu)}^2 = \|H - \mathbb{E}_A[H]\|_{2,\sigma(\mu)}^2, \quad (57) $$

where the first inequality follows by the $L_2(\sigma(\mu))$ contractivity of $\operatorname{id} - \mathbb{E}_A$. Now, the weighted norm can be further simplified into a sum of local weighted norms as follows: first, for any two $E_i, E_j$ whose supports intersect with a different vertex of $A$, we have

$$ \begin{aligned} \langle \mathbb{E}_A[E_i], \mathbb{E}_A[E_j] \rangle_{\sigma(\mu)} &= \operatorname{tr} [\sigma(\mu)^{1/2} \mathbb{E}_A[E_i] \sigma(\mu)^{1/2} \mathbb{E}_A[E_j]] \\ &= \operatorname{tr} [\sigma(\mu) \Delta_{\sigma(\mu)}^{-1/2} \circ \mathbb{E}_A[E_i] \mathbb{E}_A[E_j]] \\ &= \operatorname{tr} [\sigma(\mu) \mathbb{E}_A[E_i] \mathbb{E}_A[E_j]] \end{aligned} \quad (58) $$

where in the third line we used the commutation of the modular operator $\Delta_{\sigma(\mu)}(X) := \sigma(\mu)X\sigma(\mu)^{-1}$ with $\mathbb{E}_A$ together with the commutativity of $\sigma(\mu)$ and $E_i$. Now, denoting $\operatorname{supp}(E_i) \cap A = \{x\}$ and $\operatorname{supp}(E_j) \cap A = \{y\}$, we show that

$$ \mathbb{E}_A[E_i] = \mathbb{E}_x[E_i] \quad \text{and} \quad \mathbb{E}_A[E_j] = \mathbb{E}_y[E_j]. \quad (59) $$

In order to prove these two identities, we simply need to prove for instance that $\mathbb{E}_x[E_i]$ belongs to the image algebra $\mathcal{N}_A$ of $\mathbb{E}_A$ since $\mathcal{N}_A \subseteq \mathcal{N}_x$ by definition. Hence, it is enough to show that $\mathbb{E}_x[E_i]$ commutes with operators of the form $\sigma(\mu)^{it} X_A \sigma(\mu)^{-it}$ for any $t \in \mathbb{R}$ and any $X_A \in \mathcal{B}(\mathcal{H}_A)$. This claim follows from

$$ \begin{aligned} \mathbb{E}_x[E_i]\sigma(\mu)^{it}X_A\sigma(\mu)^{-it} &= \sigma(\mu)^{it}\sigma(\mu)^{-it}\mathbb{E}_x[E_i]\sigma(\mu)^{it}X_A\sigma(\mu)^{-it} \\ &= \sigma(\mu)^{it}\mathbb{E}_x[\sigma(\mu)^{-it}E_i\sigma(\mu)^{it}]X_A\sigma(\mu)^{-it} \\ &= \sigma(\mu)^{it}\mathbb{E}_x[E_i]X_A\sigma(\mu)^{-it} \\ &= \sigma(\mu)^{it}X_A\mathbb{E}_x[E_i]\sigma(\mu)^{-it} \\ &= \sigma(\mu)^{it}X_A\sigma(\mu)^{-it}\mathbb{E}_x[E_i]. \end{aligned} $$

where the fourth line follows from the fact that the support of $\mathbb{E}_x[E_i]$ does not intersect $A \setminus \{i\}$, together with the fact that $\mathbb{E}_x[E_i]$ is locally proportional to $I$ on site $x$, by definition of $\mathcal{N}_x$. Therefore, using (59) into (58), we get

$$ \begin{aligned} \langle \mathbb{E}_A[E_i], \mathbb{E}_A[E_j] \rangle_{\sigma(\mu)} &= \operatorname{tr} [\sigma(\mu)\mathbb{E}_x[E_i]\mathbb{E}_y[E_j]] \\ &= \operatorname{tr} [\sigma(\mu)\mathbb{E}_x\mathbb{E}_y[E_iE_j]] \\ &= \operatorname{tr} [\sigma(\mu)E_iE_j] \\ &= \langle E_i, E_j \rangle_{\sigma(\mu)}, \end{aligned} $$

where in the second line we used that $\mathbb{E}_x[E_i]$ is a fixed point of $\mathbb{E}_y$, and then that $E_j$ is a fixed point of $\mathbb{E}_x$, by the support conditions of $E_i$ and $E_j$. Therefore, the variance on the right-hand side of (57) can be simplified as

$$ \begin{aligned} \|H - \mathbb{E}_A[H]\|_{2,\sigma(\mu)}^2 &= \sum_{x \in A} \left\| \sum_{j | \operatorname{supp}(E_j) \ni x} \lambda_j (E_j - \mathbb{E}_x[E_j]) \right\|_{2,\sigma(\mu)}^2 \\ &= \sum_{x \in A} \left\| H_x - \mathbb{E}_x[H_x] \right\|_{2,\sigma(\mu)}^2 \end{aligned} $$
where we recall that $H_x := \sum_{j \mid \text{supp}(E_j) \ni x} \lambda_j E_j$. Now, for any $x \in V$, we denote $x\partial := \{A \in E | x \in A\}$ and decompose the Hamiltonian $H(\mu)$ as

$$H(\mu) = H_x(\mu) + K_x^0(\mu) + K_x^1(\mu), \quad \text{where} \quad \left\{ \begin{array}{l} \text{supp}(K_x^1(\mu)) \cap x\partial = \emptyset \\ \text{supp}(K_x^0(\mu)) \cap x\partial \neq \emptyset \end{array} \right. . \qquad (60)$$

Clearly,

$$\sigma(\mu) \geq \sigma_x(\mu) e^{-2\beta \|K_x^0(\mu)\|_{\infty}} \frac{e^{-\beta K_x^1(\mu)}}{\text{tr}_{x\partial^c} [e^{-\beta K_x^1(\mu)}]}. \qquad (61)$$

Now,

$$\begin{aligned} \|H_x - \mathbb{E}_x[H_x]\|_{2,\sigma(\mu)}^2 &= \text{tr} [\sigma(\mu)(H_x - \mathbb{E}_x[H_x])^2] \\ &\geq e^{-2\beta} \|K_x^0(\mu)\|_{\infty} \text{tr} \left[ \sigma_x(\mu) \frac{e^{-\beta K_x^1(\mu)}}{\text{tr}_{x\partial^c} [e^{-\beta K_x^1(\mu)}]} (H_x - \mathbb{E}_x[H_x])^2 \right] \\ &= e^{-2\beta} \|K_x^0(\mu)\|_{\infty} \|H_x - \mathbb{E}_x[H_x]\|_{2,\sigma_x(\mu)}^2 \\ &= e^{-2\beta} \|K_x^0(\mu)\|_{\infty} (\|H_x - \text{tr}[\sigma_x(\mu)H_x] I\|_{2,\sigma_x(\mu)}^2 - \|\mathbb{E}_x[H_x] - \text{tr}[\sigma_x(\mu)H_x] I\|_{2,\sigma_x(\mu)}^2) \\ &\geq e^{-2\beta} \|K_x^0(\mu)\|_{\infty} (\|H_x - \text{tr}[\sigma_x(\mu)H_x] I\|_{2,\sigma_x(\mu)}^2 - \|\mathcal{R}_{x,\sigma(\mu)}[H_x] - \text{tr}[\sigma_x(\mu)H_x] I\|_{2,\sigma_x(\mu)}^2) \\ &\geq (1 - c(\beta)^2) e^{-2\beta} \|K_x^0(\mu)\|_{\infty} \|H_x - \text{tr}[\sigma_x(\mu)H_x] I\|_{2,\sigma_x(\mu)}^2. \end{aligned}$$

The first and second identities above follow once again from the commutativity of the Hamiltonian similarly to (58), where for the second one we also use the disjointness of $x\partial$ and $\text{supp}(K_x^1(\mu))$. The first inequality follows from (61). The third identity is a consequence of the fact that $\mathbb{E}_x$ is a projection with respect to $L_2(\sigma_x(\mu))$. The second inequality follows from $\|\mathcal{R}_{x,\sigma(\mu)}[X]\|_{2,\sigma(\mu)} \geq \|\mathbb{E}_x[X]\|_{2,\sigma(\mu)}$ for all $X$ (see Proposition 10 in [44]). The last inequality is a consequence of Lemma F.2. Finally, we further bound the weighted $L_2$ norm on the last line of the above inequality in terms of the Schatten 2 norm to get

$$\text{Var}_{\sigma(\mu)}(H) \geq (1 - c(\beta)^2) \min_{x \in V} e^{-2\beta \|K_x^0(\mu)\|_{\infty}} \lambda_{\min}(\sigma_x(\mu)) \sum_{x \in A} \|H_x - \text{tr}[\sigma_x(\mu) H_x]\|_2^2 \quad (62)$$

$$\geq (1 - c(\beta)^2) \min_{x \in V} e^{-2\beta \|K_x^0(\mu)\|_{\infty}} \lambda_{\min}(\sigma_x(\mu)) \sum_i \lambda_i^2, \qquad (63)$$

where $\lambda_{\min}(\sigma_x(\mu))$ denotes the minimum eigenvalue of $\sigma_x(\mu)$. The result follows from the simple estimates $\|K_x^0(\mu)\|_{\infty} \leq B(4r_0)$ and $\lambda_{\min}(\sigma_x(\mu)) \geq e^{-\beta B(2r_0)} d^{-B(2r_0)}$.

## F.3 Summary for 1D or high-temperature commuting

Now we genuinely have all the elements in place to essentially give the last word on estimating Lipschitz observables for Gibbs states of nearest neighbor 1D or high-temperature commuting Gibbs states.

**Theorem F.1.** 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=1}^m \mu_j E_j$ on the hypergraph $G = (V, E)$ defined in (55), where $\text{tr}[E_i E_j] = 0$ for all $i \neq j$, each local operator $E_j$ is traceless on its support and acts on a constant number of qubits and $m = \mathcal{O}(n)$. Moreover, assume that $\sigma(\lambda)$ satisfies the conditions of Thm. B.1. Then $\mathcal{O}(\log(n)\epsilon^{-2})$ samples of $\sigma(\lambda)$ suffice to obtain a $\mu \in B_{\ell_\infty}(0,1)$ satisfying

$$\text{tr}[\mathcal{O}(\sigma(\lambda) - \sigma(\mu))] = \mathcal{O}(\epsilon\sqrt{n}\|\mathcal{O}\|_{\text{Lip},\nabla}).$$

with probability at least $1-p$. Moreover, we can find $\mu$ in $\mathcal{O}(\text{poly}(n, \epsilon^{-1}))$ time on a classical computer. With access to a quantum computer, the postprocessing can be performed in $\tilde{\mathcal{O}}(n\epsilon^{-2})$ time by only implementing quantum circuits of $\tilde{\mathcal{O}}(1)$ depth.

*Proof.* To obtain the claim on the sample complexity, note that for such systems $L = \Omega(1)$ by Prop. F.1 and they satisfy a transportation cost inequality by Thm. B.1. Moreover, we can learn the expectation of all $E_i$ up to an error $\epsilon$ and failure probability $\delta$ with $\mathcal{O}(\log(n\delta^{-1})\epsilon^{-2})$ samples using shadows, as they all have constant
where we recall that $H_x := \sum_{j \mid \text{supp}(E_j) \ni x} \lambda_j E_j$. Now, for any $x \in V$, we denote $x\partial := \{A \in E | x \in A\}$ and decompose the Hamiltonian $H(\mu)$ as

$$H(\mu) = H_x(\mu) + K_x^0(\mu) + K_x^1(\mu), \quad \text{where} \quad \left\{ \begin{array}{l} \text{supp}(K_x^1(\mu)) \cap x\partial = \emptyset \\ \text{supp}(K_x^0(\mu)) \cap x\partial \neq \emptyset \end{array} \right. . \qquad (60)$$

Clearly,

$$\sigma(\mu) \geq \sigma_x(\mu) e^{-2\beta \|K_x^0(\mu)\|_{\infty}} \frac{e^{-\beta K_x^1(\mu)}}{\text{tr}_{x\partial^c} [e^{-\beta K_x^1(\mu)}]}. \qquad (61)$$

Now,

$$\begin{aligned} \|H_x - \mathbb{E}_x[H_x]\|_{2,\sigma(\mu)}^2 &= \text{tr} [\sigma(\mu)(H_x - \mathbb{E}_x[H_x])^2] \\ &\geq e^{-2\beta \|K_x^0(\mu)\|_{\infty}} \text{tr} \left[ \sigma_x(\mu) \frac{e^{-\beta K_x^1(\mu)}}{\text{tr}_{x\partial^c} [e^{-\beta K_x^1(\mu)}]} (H_x - \mathbb{E}_x[H_x])^2 \right] \\ &= e^{-2\beta \|K_x^0(\mu)\|_{\infty}} \|H_x - \mathbb{E}_x[H_x]\|_{2,\sigma_x(\mu)}^2 \\ &= e^{-2\beta \|K_x^0(\mu)\|_{\infty}} (\|H_x - \text{tr}[\sigma_x(\mu)H_x] I\|_{2,\sigma_x(\mu)}^2 - \|\mathbb{E}_x[H_x] - \text{tr}[\sigma_x(\mu)H_x] I\|_{2,\sigma_x(\mu)}^2) \\ &\geq e^{-2\beta \|K_x^0(\mu)\|_{\infty}} (\|H_x - \text{tr}[\sigma_x(\mu)H_x] I\|_{2,\sigma_x(\mu)}^2 - \|\mathcal{R}_{x,\sigma(\mu)}[H_x] - \text{tr}[\sigma_x(\mu)H_x] I\|_{2,\sigma_x(\mu)}^2) \\ &\geq (1 - c(\beta)^2) e^{-2\beta \|K_x^0(\mu)\|_{\infty}} \|H_x - \text{tr}[\sigma_x(\mu)H_x] I\|_{2,\sigma_x(\mu)}^2. \end{aligned}$$

The first and second identities above follow once again from the commutativity of the Hamiltonian similarly to (58), where for the second one we also use the disjointness of $x\partial$ and $\text{supp}(K_x^1(\mu))$. The first inequality follows from (61). The third identity is a consequence of the fact that $\mathbb{E}_x$ is a projection with respect to $L_2(\sigma_x(\mu))$. The second inequality follows from $\|\mathcal{R}_{x,\sigma(\mu)}[X]\|_{2,\sigma(\mu)} \geq \|\mathbb{E}_x[X]\|_{2,\sigma(\mu)}$ for all $X$ (see Proposition 10 in [44]). The last inequality is a consequence of Lemma F.2. Finally, we further bound the weighted $L_2$ norm on the last line of the above inequality in terms of the Schatten 2 norm to get

$$\text{Var}_{\sigma(\mu)}(H) \geq (1 - c(\beta)^2) \min_{x \in V} e^{-2\beta \|K_x^0(\mu)\|_{\infty}} \lambda_{\min}(\sigma_x(\mu)) \sum_{x \in A} \|H_x - \text{tr}[\sigma_x(\mu) H_x]\|_2^2 \qquad (62)$$

$$\geq (1 - c(\beta)^2) \min_{x \in V} e^{-2\beta \|K_x^0(\mu)\|_{\infty}} \lambda_{\min}(\sigma_x(\mu)) \sum_i \lambda_i^2, \qquad (63)$$

where $\lambda_{\min}(\sigma_x(\mu))$ denotes the minimum eigenvalue of $\sigma_x(\mu)$. The result follows from the simple estimates $\|K_x^0(\mu)\|_{\infty} \leq B(4r_0)$ and $\lambda_{\min}(\sigma_x(\mu)) \geq e^{-\beta B(2r_0)} d^{-B(2r_0)}$.

□

## F.3 Summary for 1D or high-temperature commuting

Now we genuinely have all the elements in place to essentially give the last word on estimating Lipschitz observables for Gibbs states of nearest neighbor 1D or high-temperature commuting Gibbs states.

**Theorem F.1.** 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=1}^m \mu_j E_j$ on the hypergraph $G = (V,E)$ defined in (55), where $\text{tr}[E_i E_j] = 0$ for all $i \neq j$, each local operator $E_j$ is traceless on its support and acts on a constant number of qubits and $m = \mathcal{O}(n)$. Moreover, assume that $\sigma(\lambda)$ satisfies the conditions of Thm. B.1. Then $\mathcal{O}(\log(n)\epsilon^{-2})$ samples of $\sigma(\lambda)$ suffice to obtain a $\mu \in B_{\ell_{\infty}}(0,1)$ satisfying

$$\text{tr}[\mathcal{O}(\sigma(\lambda) - \sigma(\mu))] = \mathcal{O}(\epsilon\sqrt{n}\|\mathcal{O}\|_{\text{Lip},\nabla}).$$

with probability at least $1-p$. Moreover, we can find $\mu$ in $\mathcal{O}(\text{poly}(n, \epsilon^{-1}))$ time on a classical computer. With access to a quantum computer, the postprocessing can be performed in $\tilde{\mathcal{O}}(n\epsilon^{-2})$ time by only implementing quantum circuits of $\tilde{\mathcal{O}}(1)$ depth.

*Proof.* To obtain the claim on the sample complexity, note that for such systems $L = \Omega(1)$ by Prop. F.1 and they satisfy a transportation cost inequality by Thm. B.1. Moreover, we can learn the expectation of all $E_i$ up to an error $\epsilon$ and failure probability $\delta$ with $\mathcal{O}(\log(n\delta^{-1})\epsilon^{-2})$ samples using shadows, as they all have constant
support. The claimed sample complexity then follows from Thm. C.1. The classical postprocessing result follows from [51]. The postprocessing with a quantum computer follows from the results of [44], [18] combined with the fact that $L^{-1}U = \mathcal{O}(1)$ by also invoking Lemma F.1. Indeed, [18, 44] asserts that we can approximately prepare any $\sigma(\mu)$ on a quantum computer with a circuit of depth $\tilde{\mathcal{O}}(1)$. Moreover, once again resorting to shadows we can estimate $e(\mu)$ with $\tilde{\mathcal{O}}(\epsilon^{-2})$ samples. We conclude that we can run each iteration of gradient descent in $\tilde{\mathcal{O}}(n\epsilon^{-2})$ time. As $L^{-1}U = \mathcal{O}(1)$, Thm. A.1 asserts that we converge after $\tilde{\mathcal{O}}(1)$ iterations, which yields the claim $\square$

The theorem above nicely showcases how the joint use of transportation cost inequalities and strong convexity bounds to improve maximum entropy methods come together. Moreover, with access to a quantum computer, up to polylogarithmic overhead in system size factors, the computational complexity of learning the Gibbs state is comparable to reading out the results of measurements on the system. Together with the polylog sample complexity bounds that we obtain, this justifies our claiming that the result above almost gives the last word on learning such states.

# G Comparison of sample complexity of previous methods

Our work arguably introduces two technical innovations to the literature of learning and tomography of quantum states that underly our exponential speedups in sample complexity. The first is the observation that most observables of physical interest have a small Lipschitz constant, and, thus, it might be more motivated to look for good approximations in Wasserstein distance instead of trace distance. The second is that for states that satisfy a TC inequality, it suffices to obtain an estimate of the state that has a small relative entropy density with the target state to recover Lipschitz observables. And that finding such an estimate can be achieved from a few samples through a combination of maximum entropy and classical shadow methods.

We will now argue that these two innovations are indeed crucial to ensure our exponential speedups. First, we will show that the shadow protocol will yield bad estimates for the expectation value of local observables with high probability even for product states if the number of samples is not exponential in the locality of the underlying observables. This shows that exploiting the locality of the underlying states is crucial to obtaining a polynomial sample complexity in the locality of the observables. After that, we will show in Subsec. G.2 that $\Omega(\sqrt{n\epsilon^{-2}})$ samples are necessary for any algorithm that can recover any Gibbs state on a regular lattice at constant temperature up to trace distance $\epsilon$. This will follow from the results of [27] and showcase the need to focus on the Wasserstein distance instead of the trace distance.

## G.1 Lower bounds for sample complexity of classical shadows

One of the main advantages of our results compared with the classical shadows method of [40] is that whenever a TC inequality is available, we can learn all $k$-local observables with a number of samples that grows polynomially in $k$, whereas classical shadows require a number of samples that grows exponentially in $k$. However, the classical shadows framework does not assume any structure for the underlying state. Thus, it is natural to ask if this undesired exponential scaling of the shadows framework is due to this broader applicability.

In this section, we will demonstrate that this is not the case even for one of the simplest imaginable classes of states, namely tensor products of Pauli eigenstates. We will show that if the number of samples is not exponential in the locality of the desired observables, there will always be a $k$-local observable whose estimate will be wrong with constant probability.

But before we show that, let us briefly recall how the shadow method works. To recover local observables on $n$ qubits, the methods of [40] proceed as follows. First, we sample a random unitary $U = \otimes_{i=0}^{n-1} U_i$, where each $U_i$ is an independent rotation into a Pauli basis. Then we proceed to measure the state of interest $\rho$ in the basis defined by $U$, obtaining a $n$-bit classical string $b_0b_1...b_{n-1}$. The shadow corresponding to this measurement $\hat{\rho}$ is then defined to be given by

$$ \hat{\rho} = \otimes_{i=1}^{n} (3U_i^\dagger |b_i\rangle\langle b_i| U_i - I). $$

We then repeat this procedure $SM$ times for $S, M \in \mathbb{N}$, obtaining shadows $\{\hat{\rho}_{s,m}\}_{1 \le m \le M, 1 \le s \le S}$. We then set our estimate of the expectation value of an observable $O$ to be given by

$$ \hat{O}(\{\hat{\rho}_{s,m}\}) = \text{median} \left\{ S^{-1} \sum_{s=1}^{S} \text{tr}[\hat{\rho}_{s,m=1}O], S^{-1} \sum_{s=1}^{S} \text{tr}[\hat{\rho}_{s,m=2}O], \dots, S^{-1} \sum_{s=1}^{S} \text{tr}[\hat{\rho}_{s,m=M}O] \right\}. \quad (64) $$
support. The claimed sample complexity then follows from Thm. C.1. The classical postprocessing result follows from [51]. The postprocessing with a quantum computer follows from the results of [44], [18] combined with the fact that $L^{-1}U = O(1)$ by also invoking Lemma F.1. Indeed, [18, 44] asserts that we can approximately prepare any $\sigma(\mu)$ on a quantum computer with a circuit of depth $\tilde{O}(1)$. Moreover, once again resorting to shadows we can estimate $e(\mu)$ with $\tilde{O}(\epsilon^{-2})$ samples. We conclude that we can run each iteration of gradient descent in $\tilde{O}(n\epsilon^{-2})$ time. As $L^{-1}U = O(1)$, Thm. A.1 asserts that we converge after $\tilde{O}(1)$ iterations, which yields the claim $\square$

The theorem above nicely showcases how the joint use of transportation cost inequalities and strong convexity bounds to improve maximum entropy methods come together. Moreover, with access to a quantum computer, up to polylogarithmic overhead in system size factors, the computational complexity of learning the Gibbs state is comparable to reading out the results of measurements on the system. Together with the polylog sample complexity bounds that we obtain, this justifies our claiming that the result above almost gives the last word on learning such states.

# G Comparison of sample complexity of previous methods

Our work arguably introduces two technical innovations to the literature of learning and tomography of quantum states that underly our exponential speedups in sample complexity. The first is the observation that most observables of physical interest have a small Lipschitz constant, and, thus, it might be more motivated to look for good approximations in Wasserstein distance instead of trace distance. The second is that for states that satisfy a TC inequality, it suffices to obtain an estimate of the state that has a small relative entropy density with the target state to recover Lipschitz observables. And that finding such an estimate can be achieved from a few samples through a combination of maximum entropy and classical shadow methods.

We will now argue that these two innovations are indeed crucial to ensure our exponential speedups. First, we will show that the shadow protocol will yield bad estimates for the expectation value of local observables with high probability even for product states if the number of samples is not exponential in the locality of the underlying observables. This shows that exploiting the locality of the underlying states is crucial to obtaining a polynomial sample complexity in the locality of the observables. After that, we will show in Subsec. G.2 that $\Omega(\sqrt{n\epsilon^{-2}})$ samples are necessary for any algorithm that can recover any Gibbs state on a regular lattice at constant temperature up to trace distance $\epsilon$. This will follow from the results of [27] and showcase the need to focus on the Wasserstein distance instead of the trace distance.

## G.1 Lower bounds for sample complexity of classical shadows

One of the main advantages of our results compared with the classical shadows method of [40] is that whenever a TC inequality is available, we can learn all $k$-local observables with a number of samples that grows polynomially in $k$, whereas classical shadows require a number of samples that grows exponentially in $k$. However, the classical shadows framework does not assume any structure for the underlying state. Thus, it is natural to ask if this undesired exponential scaling of the shadows framework is due to this broader applicability.

In this section, we will demonstrate that this is not the case even for one of the simplest imaginable classes of states, namely tensor products of Pauli eigenstates. We will show that if the number of samples is not exponential in the locality of the desired observables, there will always be a $k$-local observable whose estimate will be wrong with constant probability.

But before we show that, let us briefly recall how the shadow method works. To recover local observables on $n$ qubits, the methods of [40] proceed as follows. First, we sample a random unitary $U = \otimes_{i=0}^{n-1} U_i$, where each $U_i$ is an independent rotation into a Pauli basis. Then we proceed to measure the state of interest $\rho$ in the basis defined by $U$, obtaining a $n$-bit classical string $b_0b_1...b_{n-1}$. The shadow corresponding to this measurement $\hat{\rho}$ is then defined to be given by

$$ \hat{\rho} = \otimes_{i=1}^{n} (3U_i^\dagger |b_i\rangle\langle b_i| U_i - I). $$

We then repeat this procedure $SM$ times for $S, M \in \mathbb{N}$, obtaining shadows $\{\hat{\rho}_{s,m}\}_{1 \le m \le M, 1 \le s \le S}$. We then set our estimate of the expectation value of an observable $O$ to be given by

$$ \hat{O}(\{\hat{\rho}_{s,m}\}) = \text{median} \left\{ S^{-1} \sum_{s=1}^{S} \text{tr}[\hat{\rho}_{s,m=1} O], S^{-1} \sum_{s=1}^{S} \text{tr}[\hat{\rho}_{s,m=2} O], \dots, S^{-1} \sum_{s=1}^{S} \text{tr}[\hat{\rho}_{s,m=M} O] \right\}. \quad (64) $$
That is, we partition the set of shadow samples into M subgroups, take the average on each of them and then take their median. The authors of [40] then proceed to show in Theorem 1 that if we take $SM = O(4^k \log(N\delta^{-1})\epsilon^{-2})$, then with probability at least $1 - \delta$ the expectation value any given N k-local observables of bounded operator norm will deviate by at most $\epsilon$ from the estimate given in Eq. (64). We will now show that this exponential scaling in k is unavoidable if we want to obtain nontrivial estimates with high probability. More precisely:

**Proposition G.1.** Let $\rho = \otimes_{i=1}^n |\phi_i\rangle\langle\phi_i|$ be an unknown n qubit state where each $|\phi_i\rangle$ is given by a Pauli eigenstate. For $k \ge \log_3(n) \log(n)$, if $SM \le n$ then there is an observable O of the form

$$O = n^{-1} \sum_{i=1}^{n} O_i$$

where each $O_i$ is k-local and such that:

$$\mathrm{tr}[O\rho] = 1, \quad \hat{O}(\{\hat{\rho}_{s,m}\}) = 0 \qquad (65)$$

with probability at least $1 - e^{-1}$.

*Proof.* We will let $O_i = \otimes_{j=i}^{i+k} P_j$, where we take addition modulo $n$ and let $P_i$ be the Pauli matrix such that $P_i |\phi_i\rangle = 1$, where we assume w.l.o.g. that each Pauli eigenstate corresponds to a +1 eigenstate. It is then clear that $\mathrm{tr}[O\rho] = 1$.

Let us now analyse the performance of shadows. Note that as the random unitaries used to obtain each sample are random Pauli bases, we have that $\mathrm{tr}[O_i\hat{\rho}_{s,m}] = 3^k$ if the unitary $U^{s,m}$ corresponds to the identity on qubits $i, i+1, \dots, i+k$ and 0 otherwise. This is because then the string we measure will be rotated to a Pauli basis different from that of $P_i$ on at least one of the qubits in that interval. Thus, as we pick the Pauli bases uniformly at random and there are three different possibilities for each one of them, we see that the probability that a shadow is different from 0 on a given $O_i$ is $3^{-k}$. By a union bound, the probability that a shadow is different than 0 on any of the $n$ $O_i$ is at most $n3^{-k}$. As the different shadows are independent, the probability that all of the SM shadows return 0 is at least $(1 - n3^{-k})^{SM}$. As we picked $k \ge \log_3(n) \log(n)$ and $SM \le n$, we have:

$$(1 - n3^{-k})^{SM} \ge (1 - n^{-1})^n \ge \frac{e^{-1}}{2}$$

for $n$ large enough. Thus, as all shadows will have expectation value 0, the median and means procedure will clearly also output 0, which concludes the proof. □

Despite the fact the proof above used quite rough estimates and simple observables and states, it still gives some intuition as to why classical shadows require an exponential number of samples in the locality. We see that the probability that the shadows "look in the right direction" is exponentially small in the locality of the observables, in the sense that the overlap between the state and a random Pauli eigenstate is likely to be exponentially small. And whenever it does look in a direction in which the underlying state has a significant overlap with that basis it has to compensate that direction exponentially. Thus, if the number of samples is not exponential in the locality, it is unlikely that we will measure in a direction that has significant overlap with our state or there will be significant fluctuations due to the exponential rewarding of the "good" directions.

However, by combining shadows with a locality structure and the maximum entropy principle, as we do in this work, we can bypass the need to measure in random directions for observables with a large locality, bypassing this exponential scaling.

We also note that the authors of [40] already proved the optimality of their protocol by only considering product states in Section 8 of their supplemental material. The main difference between their proof and ours is that we focus on a Lipschitz observable, whereas they focus on observables that only depend on k qubits.

## G.2 Lower bounds for recovery in trace distance

In the main text, we claimed that one of the reasons why we obtain an exponential speedup compared to usual many-body methods is that we focus on a good recovery in the Wasserstein distance instead of trace distance. Moreover, by combining the maximum entropy method with a TC inequality we are able to obtain good recovery guarantees from a constant relative entropy density. That is, as long as two Gibbs states $\sigma(\lambda), \sigma(\mu)$ satisfy

$$D(\sigma(\lambda)\|\sigma(\mu)) \le \epsilon n, \qquad (66)$$

for some $\epsilon > 0$ we already obtain some nontrivial guarantees.
That is, we partition the set of shadow samples into M subgroups, take the average on each of them and then take their median. The authors of [40] then proceed to show in Theorem 1 that if we take $SM = O(4^k \log(N\delta^{-1})\epsilon^{-2})$, then with probability at least $1 - \delta$ the expectation value any given N k-local observables of bounded operator norm will deviate by at most $\epsilon$ from the estimate given in Eq. (64). We will now show that this exponential scaling in k is unavoidable if we want to obtain nontrivial estimates with high probability. More precisely:

**Proposition G.1.** Let $\rho = \otimes_{i=1}^n |\phi_i\rangle\langle\phi_i|$ be an unknown n qubit state where each $|\phi_i\rangle$ is given by a Pauli eigenstate. For $k \ge \log_3(n) \log(n)$, if $SM \le n$ then there is an observable O of the form

$$O = n^{-1} \sum_{i=1}^{n} O_i$$

where each $O_i$ is k-local and such that:

$$\mathrm{tr}[O\rho] = 1, \quad \hat{O}(\{\hat{\rho}_{s,m}\}) = 0 \qquad (65)$$

with probability at least $1 - e^{-1}$.

*Proof.* We will let $O_i = \otimes_{j=i}^{i+k} P_j$, where we take addition modulo $n$ and let $P_i$ be the Pauli matrix such that $P_i |\phi_i\rangle = 1$, where we assume w.l.o.g. that each Pauli eigenstate corresponds to a +1 eigenstate. It is then clear that $\mathrm{tr}[O\rho] = 1$.

Let us now analyse the performance of shadows. Note that as the random unitaries used to obtain each sample are random Pauli bases, we have that $\mathrm{tr}[O_i\hat{\rho}_{s,m}] = 3^k$ if the unitary $U^{s,m}$ corresponds to the identity on qubits $i, i+1, \dots, i+k$ and 0 otherwise. This is because then the string we measure will be rotated to a Pauli basis different from that of $P_i$ on at least one of the qubits in that interval. Thus, as we pick the Pauli bases uniformly at random and there are three different possibilities for each one of them, we see that the probability that a shadow is different from 0 on a given $O_i$ is $3^{-k}$. By a union bound, the probability that a shadow is different than 0 on any of the $n$ $O_i$ is at most $n3^{-k}$. As the different shadows are independent, the probability that all of the SM shadows return 0 is at least $(1 - n3^{-k})^{SM}$. As we picked $k \ge \log_3(n) \log(n)$ and $SM \le n$, we have:

$$(1 - n3^{-k})^{SM} \ge (1 - n^{-1})^n \ge \frac{e^{-1}}{2}$$

for $n$ large enough. Thus, as all shadows will have expectation value 0, the median and means procedure will clearly also output 0, which concludes the proof. □

Despite the fact the proof above used quite rough estimates and simple observables and states, it still gives some intuition as to why classical shadows require an exponential number of samples in the locality. We see that the probability that the shadows "look in the right direction" is exponentially small in the locality of the observables, in the sense that the overlap between the state and a random Pauli eigenstate is likely to be exponentially small. And whenever it does look in a direction in which the underlying state has a significant overlap with that basis it has to compensate that direction exponentially. Thus, if the number of samples is not exponential in the locality, it is unlikely that we will measure in a direction that has significant overlap with our state or there will be significant fluctuations due to the exponential rewarding of the "good" directions.

However, by combining shadows with a locality structure and the maximum entropy principle, as we do in this work, we can bypass the need to measure in random directions for observables with a large locality, bypassing this exponential scaling.

We also note that the authors of [40] already proved the optimality of their protocol by only considering product states in Section 8 of their supplemental material. The main difference between their proof and ours is that we focus on a Lipschitz observable, whereas they focus on observables that only depend on k qubits.

## G.2 Lower bounds for recovery in trace distance

In the main text, we claimed that one of the reasons why we obtain an exponential speedup compared to usual many-body methods is that we focus on a good recovery in the Wasserstein distance instead of trace distance. Moreover, by combining the maximum entropy method with a TC inequality we are able to obtain good recovery guarantees from a constant relative entropy density. That is, as long as two Gibbs states $\sigma(\lambda), \sigma(\mu)$ satisfy

$$D(\sigma(\lambda)\|\sigma(\mu)) \le \epsilon n, \qquad (66)$$

for some $\epsilon > 0$ we already obtain some nontrivial guarantees.
In this section, we will argue that focusing on the Wasserstein recovery instead of trace distance is essential to obtain nontrivial recovery guarantees from a number of samples scaling logarithmically with the system size. To achieve this, we will resort to results of [27, Theorem 1.3]:

**Proposition G.2.** [Lower bound of sample complexity in trace distance] Let $G = (V, E)$ be a graph on $n$ vertices and $m$ edges and for $\lambda \in \mathbb{R}^{15m}$, $\|\lambda\|_{\ell_{\infty}} \le 1$ let $H(\lambda)$ be defined as

$$H(\lambda) = \sum_{i \sim j} \sum_{l=1}^{15} \lambda_{i,j}^{l} H_{i,j}^{l},$$

where the $H_{i,j}^{l}$ correspond to some ordering of the nonidentity Pauli strings acting on sites $i, j$. Then for any $\beta = \Omega(m^{-\frac{1}{2}})$, let $\hat{\sigma}(\lambda)$ be the estimate of $\sigma(\lambda)$ outputted by an algorithm with access to $s$ samples from a state $\sigma(\lambda)$. Then:

$$\sup_{\lambda \in B_{\infty}(0,1)} \|\sigma(\lambda) - \hat{\sigma}(\lambda)\|_{\text{tr}} = \Omega\left(\min\left\{1, \sqrt{\frac{m}{s}}\right\}\right). \qquad (67)$$

*Proof.* This statement immediately follows from [27, Theorem 1.3], which shows the analogous statement when restricted to classical Ising models. As our class of Hamiltonians includes those as a subset, any algorithm that could provide an estimate for this more general class also can find one for the classical instances. Moreover, in the proof of [27, Theorem 1.3] the inverse temperature is absorbed into the coefficients of the Hamiltonian, which are assumed to have 2-norm bounded by a constant independent of the system's size. This is easily seen to be satisfied by our conditions since $\beta = \Omega(m^{-\frac{1}{2}})$. $\square$

The statement above implies in particular that any algorithm that finds an estimate that is $\epsilon$ close in trace distance for all 2-local Gibbs state on a lattice and constant inverse temperature requires $\Omega(n\epsilon^{-2})$ samples. In contrast, we see that it is possible to obtain an estimate that is $O(\epsilon\sqrt{n})$ close in Wasserstein distance from $O(\epsilon^{-2}\log(n))$ samples of the Gibbs states, which is sufficient to already give nontrivial recovery guarantees for Lipschitz observables. Thus, we see from Prop. G.2 that resorting to the Wasserstein distance is essential to obtain recovery guarantees in the regime where the number of samples is logarithmic in the system's size.

Furthermore, it is interesting to note that the proof of [27] is based on a set of Gibbs states of the form:

$$\delta \sum_{i \sim j} s_{i,j} Z_i Z_j, \qquad (68)$$

with $\delta = \Theta(m^{-\frac{1}{2}})$ and $s_{i,j} \in \{\pm 1\}$. Their proof then proceeds by finding a large subset of Gibbs states of the form in Eq. (68) which have a trace distance and relative entropy of constant order. The lower bound on the sample complexity then follows from standard information-theoretic arguments. We believe that this class of examples in the proof further illustrates why the trace distance is not necessarily the adequate distance measure when estimating the error on extensive observables. Indeed, for extensive, local observables the class of Gibbs states from the Hamiltonians in Eq. (68) behaves like the maximally mixed state, as each local term converges to 0 as the system size increases.
In this section, we will argue that focusing on the Wasserstein recovery instead of trace distance is essential to obtain nontrivial recovery guarantees from a number of samples scaling logarithmically with the system size. To achieve this, we will resort to results of [27, Theorem 1.3]:

**Proposition G.2.** [Lower bound of sample complexity in trace distance] Let $G = (V, E)$ be a graph on $n$ vertices and $m$ edges and for $\lambda \in \mathbb{R}^{15m}$, $\|\lambda\|_{\ell_{\infty}} \le 1$ let $H(\lambda)$ be defined as

$$H(\lambda) = \sum_{i \sim j} \sum_{l=1}^{15} \lambda_{i,j}^{l} H_{i,j}^{l},$$

where the $H_{i,j}^{l}$ correspond to some ordering of the nonidentity Pauli strings acting on sites $i, j$. Then for any $\beta = \Omega(m^{-\frac{1}{2}})$, let $\hat{\sigma}(\lambda)$ be the estimate of $\sigma(\lambda)$ outputted by an algorithm with access to $s$ samples from a state $\sigma(\lambda)$. Then:

$$\sup_{\lambda \in B_{\infty}(0,1)} \|\sigma(\lambda) - \hat{\sigma}(\lambda)\|_{\text{tr}} = \Omega\left(\min\left\{1, \sqrt{\frac{m}{s}}\right\}\right). \qquad (67)$$

*Proof.* This statement immediately follows from [27, Theorem 1.3], which shows the analogous statement when restricted to classical Ising models. As our class of Hamiltonians includes those as a subset, any algorithm that could provide an estimate for this more general class also can find one for the classical instances. Moreover, in the proof of [27, Theorem 1.3] the inverse temperature is absorbed into the coefficients of the Hamiltonian, which are assumed to have 2-norm bounded by a constant independent of the system's size. This is easily seen to be satisfied by our conditions since $\beta = \Omega(m^{-\frac{1}{2}})$. $\square$

The statement above implies in particular that any algorithm that finds an estimate that is $\epsilon$ close in trace distance for all 2-local Gibbs state on a lattice and constant inverse temperature requires $\Omega(n\epsilon^{-2})$ samples. In contrast, we see that it is possible to obtain an estimate that is $O(\epsilon\sqrt{n})$ close in Wasserstein distance from $O(\epsilon^{-2}\log(n))$ samples of the Gibbs states, which is sufficient to already give nontrivial recovery guarantees for Lipschitz observables. Thus, we see from Prop. G.2 that resorting to the Wasserstein distance is essential to obtain recovery guarantees in the regime where the number of samples is logarithmic in the system's size.

Furthermore, it is interesting to note that the proof of [27] is based on a set of Gibbs states of the form:

$$\delta \sum_{i \sim j} s_{i,j} Z_i Z_j, \qquad (68)$$

with $\delta = \Theta(m^{-\frac{1}{2}})$ and $s_{i,j} \in \{\pm 1\}$. Their proof then proceeds by finding a large subset of Gibbs states of the form in Eq. (68) which have a trace distance and relative entropy of constant order. The lower bound on the sample complexity then follows from standard information-theoretic arguments. We believe that this class of examples in the proof further illustrates why the trace distance is not necessarily the adequate distance measure when estimating the error on extensive observables. Indeed, for extensive, local observables the class of Gibbs states from the Hamiltonians in Eq. (68) behaves like the maximally mixed state, as each local term converges to 0 as the system size increases.