# Learning mixed quantum states in large-scale experiments

Matteo Votto, Marko Ljubotina, Cécilia Lancien, J. Ignacio Cirac, Peter

Zoller, Maksym Serbyn, Lorenzo Piroli, Benoît Vermersch

► To cite this version:

Matteo Votto, Marko Ljubotina, Cécilia Lancien, J. Ignacio Cirac, Peter Zoller, et al.. Learning mixed quantum states in large-scale experiments. 2025. hal-05350401

HAL Id: hal-05350401

https://hal.science/hal-05350401v1

Preprint submitted on 6 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
Learning mixed quantum states in large-scale experiments

Matteo Votto,¹, * Marko Ljubotina,²,³,⁴ Cécilia Lancien,⁵ J. Ignacio Cirac,⁶,³
Peter Zoller,⁷,⁸ Maksym Serbyn,⁴ Lorenzo Piroli,⁹ and Benoît Vermersch¹,¹⁰

¹Université Grenoble Alpes, CNRS, LPMMC, 38000 Grenoble, France
²Technical University of Munich, TUM School of Natural Sciences,
Physics Department, James-Franck-Straße 1, 85748 Garching, Germany

³Munich Center for Quantum Science and Technology (MCQST), Schellingstr. 4, 80799 München, Germany

⁴Institute of Science and Technology Austria (ISTA), Am Campus 1, 3400 Klosterneuburg, Austria

⁵Université Grenoble Alpes, CNRS, Institut Fourier, 38610 Gières, France

⁶Max Planck Institute of Quantum Optics, Hans-Kopfermann-Str. 1, Garching 85748, Germany

⁷Institute for Theoretical Physics, University of Innsbruck, Innsbruck A-6020, Austria

⁸Institute for Quantum Optics and Quantum Information of the Austrian Academy of Sciences, Innsbruck A-6020, Austria

⁹Dipartimento di Fisica e Astronomia, Università di Bologna and INFN,
Sezione di Bologna, via Irnerio 46, I-40126 Bologna, Italy

¹⁰Quobly, Grenoble, France

We present and test a protocol to learn the matrix-product operator (MPO) representation of an experimentally prepared quantum state. The protocol takes as an input classical shadows corresponding to local randomized measurements, and outputs the tensors of a MPO which maximizes a suitably-defined fidelity with the experimental state. The tensor optimization is carried out sequentially, similarly to the well-known density matrix renormalization group algorithm. Our approach is provably efficient under certain technical conditions which are expected to be met in short-range correlated states and in typical noisy experimental settings. Under the same conditions, we also provide an efficient scheme to estimate fidelities between the learned and the experimental states. We experimentally demonstrate our protocol by learning entangled quantum states of up to $N = 96$ qubits in a superconducting quantum processor. Our method upgrades classical shadows to large-scale quantum computation and simulation experiments.

Probing the quantum state of an $N$-qubit system is a crucial yet non-trivial step in the successful implementation of many quantum simulation and quantum computation protocols [1–11]. In this context, randomized measurements provide an effective strategy to probe several physical properties of experimental quantum systems [12–14], such as entanglement [15–21] and fidelities [22–26]. In particular, the framework of classical shadows [27] has become a routine tool to post-process randomized measurement datasets. This is due to (i) provable performance guarantees for statistical errors, that allow us to reach considerable system sizes [21–23, 27], (ii) minimal hardware requirements (applying one layer of single-qubit unitaries before measurements suffices) [13], and (iii) robustness to measurement errors [28–30]. While classical shadows are state-agnostic, typical quantum states generated in noisy quantum devices admit simple descriptions. For instance, matrix-product operators (MPOs) can accurately describe output states of one-dimensional noisy quantum circuits [31, 32], as well as one-dimensional thermal states relevant to quantum simulation [33, 34], with only $O(N)$ parameters. However, it remains unclear how to efficiently reconstruct such quantum states from classical shadows at large system sizes.

In this Letter, we introduce a protocol to learn the MPO representation of quantum states in large-scale experiments, that inherits the key properties (i, ii, iii) of classical shadows. This allows us to probe global properties of quantum states with extensive entropies in experimental settings. We show this by learning many-body states generated by a superconducting quantum processor up to $N = 96$ qubits, whereas randomized measurements have been previously implemented only up to $N = 13$ [30]. The MPO representation provides direct access to several physical properties, without the need of tailored estimation formulas, nor to re-process the dataset as in previous approaches. Leveraging this full mixed-state description, we also investigate and quantify the effect of experimental noise and demonstrate large-scale error mitigation.

The MPO representation of a quantum state $\sigma$ is specified by a set of $N$ tensors $M^{(j)}$ as

$$ \sigma = \sum_{\{s_j\}, \{s'_j\}} M_{s_1, s'_1}^{(1)} M_{s_2, s'_2}^{(2)} \dots M_{s_N, s'_N}^{(N)} |\{s_j\}\rangle \langle\{s'_j\}|, \quad (1) $$

where $|\{s_j\}\rangle = \otimes_{j=1}^N |s_j\rangle$, with $s_j = 0, 1$, the computational basis corresponding to the Hilbert space of $N$ qubits, while $M_{s_j, s'_j}^{(j)}$ are $\chi_j \times \chi_{j+1}$ matrices. The bond dimension $\chi = \max(\chi_j)$ controls the expressivity of this ansatz [35]: while any quantum state can be described by a MPO with $\chi = 2^N$, typical noisy [31–33] or thermal [34] quantum states in one dimension only require $\chi = O(1)$. Other approaches to quantum state tomography based on MPO have been proposed in recent years. They consist of either numerical methods without analytical performance guarantees [36–40], or provably efficient proposals [(i)] [41–46], but incompatible with local
<figure><img src="image_2.png" /><figcaption>FIG. 1. Learning protocol – (a) We present a protocol to learn the MPO representation $\sigma$ of an experimental quantum state $\rho$ from two randomized measurement datasets. We use the first dataset to optimize $\sigma$, and the second one to benchmark the output $\sigma$ by estimating the fidelity $F_{\text{max}}(\rho, \sigma)$. (b) We optimize each tensor $M^{(j)}$ approximately by solving the linear system Eq. (4), represented diagrammatically. Note that it requires experimental data depending only on the reduced density matrix $\rho_{I_j}$, where $I_j = [j - \ell, j + \ell]$. (c) We experimentally demonstrate the protocol by learning the state of the first $N$ qubits out of $N_{\text{tot}} = 96$ qubits on a superconducting quantum processor, where $\rho$ is obtained by implementing Eq. (8) at depth $D = 1$. We set $(\ell, \chi') = (2, 4)$ in the learning algorithm, and use a total of $N_u \times N_M = 2048 \times 1024$ measurements. We compare the AFC fidelity [Eq. (7)] of the experimental state $\rho$ with the learned MPO $\sigma$ (marked lines) and with the target state $|\psi\rangle$ (dotted lines). (d) We estimate the Rényi entropy $S_2$ using both our protocol (dashed line) and the testing set using AFC (markers).</figcaption></figure>

randomized measurements [(ii)] or the robust classical shadow framework [(iii)]. Similar considerations apply to protocols for the tomography of Gibbs states [47–54], another candidate ansatz to learn mixed states in experiments [55]. Here we obtain the properties (i, ii, iii) by combining bounds on statistical errors from the classical shadow framework, with the analytical properties and computational efficiency of MPOs.

Our protocol is illustrated in Fig. 1: we perform randomized measurements on the state $\rho$ in $N_u = N_u^L + N_u^T$ local random bases with $N_M$ shots per basis, and assign $N_u^L$ bases to the learning set, and $N_u^T$ bases to the testing set. The learning set serves as input for a tensor network learning algorithm where each tensor of a MPO $\sigma$ is optimized sequentially, similarly to density-matrix renormalization group (DMRG) and related approaches [56–59]. Then, we use the testing set to benchmark the result by estimating the max fidelity [24, 60]

$$F_{\text{max}}(\rho, \sigma) = \frac{\text{tr}[\rho\sigma]}{\max(\text{tr}[\rho^2], \text{tr}[\sigma^2])} \quad (2)$$

Importantly, we prove that the optimization of the tensors $M^{(j)}$ and the estimation of $F_{\text{max}}(\rho, \sigma)$ can be performed with $O(\text{poly}(N))$ measurements under mild technical assumptions. In the remainder of the manuscript we illustrate our protocol and its experimental implementation, together with a demonstration and discussions of possible applications.

**Learning algorithm**— In order to illustrate and justify all the steps of our algorithm, we make some assumptions on the target state $\rho$ to be learned: (a) $\rho$ is a MPO with a bond dimension $\chi$, and (b) both the density matrices $\rho$

and $\rho' \propto \rho^2$ have finite correlation lengths. The latter assumption allows for an efficient estimation of fidelities, as we comment later, and implies the following approximate factorization conditions (AFC) for the purity [21, 54]

$$\left| \text{tr}[\rho_{ABC}^2]^{-1} \frac{\text{tr}[\rho_{AB}^2] \text{tr}[\rho_{BC}^2]}{\text{tr}[\rho_B^2]} - 1 \right| \le \alpha e^{-|B|/\xi_{\rho}^{(2)}}, \quad (3)$$

where $\alpha > 0$ is a constant, $\xi_{\rho}^{(2)} > 0$ is the correlation length of $\rho' \propto \rho^2$, and $\rho_I = \text{tr}_{\bar{I}}[\rho]$, where A, B, and C form a tripartition for a connected subset of qubits. The correlation length $\xi_{\rho}^{(2)}$ is also referred to as the Rényi-Markov length of $\rho$ [61, 62]. Note that $\xi_{\rho}^{(2)}$ is not related to the injectivity length of matrix-product states (MPSs) [63], which is the relevant length scale in MPS tomography [64, 65]. Later, we will comment on the applicability of our method when these assumptions are only approximately met.

The learning algorithm starts from an initial MPO ansatz $\sigma$ with a given bond dimension $\chi'$. Then, we update sequentially the tensors $M^{(j)}$ in Eq. (1) by maximizing the geometric-mean (GM) fidelity $F_{\text{GM}}(\rho, \sigma) = \frac{\text{tr}[\rho\sigma]}{\sqrt{\text{tr}[\rho^2]\text{tr}[\sigma^2]}}$ [24, 60]. Unlike the max fidelity, the GM fidelity is differentiable, but at the same time it is insensitive to certain types of decoherence [24, 66]. The key result of this work is that, if $\rho$ is a translation-invariant MPO satisfying the assumptions above, and given a parameter $\ell > O(\xi_{\rho}^{(2)})$, the tensor $M^{(j)}$ that maximizes the GM fidelity approximately satisfies [66, 67]

$$\text{tr}_{I_j \setminus \{j\}}[\sigma_{I_j} \partial_{M^{(j)}} \sigma_{I_j}] = \text{tr}_{I_j \setminus \{j\}}[\rho_{I_j} \partial_{M^{(j)}} \sigma_{I_j}], \quad (4)$$

where the trace is taken over all the qubits in $I_j$ but $j$,
and $I_j = [j - \ell, j + \ell]$, see Fig. 1(b) for a diagrammatic representation. More specifically, when the algorithm is close to convergence, we can prove that $||M^{(j)} - M_0^{(j)}||_2 = O(\text{tr}[\rho_{I_j} \sigma_{I_j}]) \chi e^{-\ell/2\xi_{\rho\sigma}}$, where $M_0^{(j)}$ is the tensor solving Eq. (4), and $\xi_{\rho\sigma}$ is the correlation length of the operator $\rho\sigma$. This approximate update is convenient, since the right-hand side of Eq. (4) can be estimated efficiently via classical shadows, as we show later.

For each tensor, Eq. (4) can be solved numerically as a linear system in $M^{(j)}$ if $\chi' \le 4^\ell$ [66]. The algorithm then proceeds by sweeping through the tensors $M^{(j)}$ from $j=1$ to $j=N$ and back, performing a total of $N_S$ sweeps. After each sweep, we estimate $\mathcal{F}_{\text{max}}(\rho, \sigma)$, which is the control parameter of the algorithm. Importantly, it is not necessary to know the values of $\xi_\rho^{(2)}$ and $\chi$ exactly: in order for the algorithm to accurately reconstruct $\rho$, it is only necessary that $\ell > O(\xi_\rho^{(2)})$ and $\chi' \ge \chi$ (however, similarly to standard MPS tomography [64, 65, 68], the number of measurements and post-processing operations increase with them). When $\rho$ is only approximately described by a MPO with bond dimension $\chi$, we could use this approach until we reach the desired accuracy, that will be constrained by truncation errors.

**Randomized measurement protocol**— At this stage, we are left with the task of estimating (a) the right-hand side of Eq. (4) and (b) $\mathcal{F}_{\text{max}}(\rho, \sigma)$ using respectively the learning and the testing set of randomized measurements obtained from $\rho$. A key feature of randomized measurements is that we can re-evaluate both these functions for any $\sigma$ using the same dataset. Using the classical shadow formalism [27], we can estimate $\text{tr}_{I_j \setminus \{j\}}[\rho_{I_j} \partial_{M^{(j)}} \sigma_{I_j}]$ as $\frac{1}{N_u^L N_M} \sum_{r=1}^{N_u^L N_M} \text{tr}_{I_j \setminus \{j\}}[\hat{\rho}_{I_j}^{(r)} \partial_{M^{(j)}} \sigma_{I_j}]$, where each

$$ \hat{\rho}_{I_j}^{(r)} = \bigotimes_{j \in I_j} \left( 3u_j^{(r)\dagger} |s_j^{(r)\rangle \langle s_j^{(r)}| u_j^{(r)} - \mathbb{I}_j \right) \quad (5) $$

is a classical shadow of $\rho_{I_j}$. Here, $\{s_j^{(r)} = 0, 1\}$ are the measurement outcomes in the computational basis, $u_j^{(r)}$ are the random unitary transformations applied to each qubit $j$, and $\mathbb{I}_j$ is the identity. Crucially, this estimation can be done in norm-2 error $\varepsilon$ with $N_u^L = O(2^\ell \chi' \varepsilon^{-2})$ measurement bases [66, 69]. This implies that statistical errors in learning each tensor $M^{(j)}$ can be made arbitrarily small with a number of measurements independent from the total number of qubits $N$.

We now propose a scheme to estimate $\mathcal{F}_{\text{max}}(\rho, \sigma)$ [Eq. (2)] with a controlled measurement budget. While its estimation generally requires $N_u^T N_M = O(4^N)$ randomized measurements [18, 24, 25], it is possible to drastically reduce this overhead under the previous assumption that the correlation lengths $\xi_\rho^{(2)}$ and $\xi_{\rho\sigma}$ are finite.

<figure><img src="image_3.png" /><figcaption>FIG. 2. Experimental results — Results for a kicked Ising model with depth $D=2$, $N_u^L = N_u^T = 4096$ and $N_M = 1024$, and $(\ell, \chi') = (3, 8)$. (a) $\mathcal{F}_{\text{max}}$ as a function of $N$, with respect to the experimental quantum state (markers) and the ideal target state (dotted line). (b) Global purity as a function of $N$, obtained from the learned MPO $\sigma$ (dashed line) and from the testing set using Eq. (6) (markers). (c, d) Expectation values of the magnetization and 2-body correlation functions (averaged over $j \in [3, 88]$), comparing values obtained from the learning protocol (blue lines), classical shadows from the testing set (red lines), and the pure state $|\psi\rangle$ (dashed lines).</figcaption></figure>

In this case, we prove that [66]

$$ \text{tr}[\rho\sigma] \approx \mathcal{O}_{\text{AFC}}^{(k)}(\rho, \sigma) = \frac{\prod_{j=1}^{[N/k]-1} \text{tr}[\rho_{A_j} A_{j+1} \sigma_{A_j} A_{j+1}]}{\prod_{j=2}^{[N/k]-1} \text{tr}[\rho_{A_j} \sigma_{A_j}]}, \quad (6) $$

where $A_j$ are neighboring regions with $|A_j| = k > O(\xi_{\rho\sigma})$. The same assumptions yield a similar result for the purity, i.e. $\text{tr}[\rho^2] \approx \mathcal{P}_{2,\text{AFC}}^{(k)}(\rho) = \mathcal{O}_{\text{AFC}}^{(k)}(\rho, \rho)$ [21]. Defining now the AFC max fidelity

$$ \mathcal{F}_{\text{AFC}}^{(k)}(\rho, \sigma) = \frac{\mathcal{O}_{\text{AFC}}^{(k)}(\rho, \sigma)}{\max(\mathcal{P}_{2,\text{AFC}}^{(k)}(\rho), \mathcal{P}_{2,\text{AFC}}^{(k)}(\sigma))}, \quad (7) $$

Eq. (6) yields $\mathcal{F}_{\text{max}}(\rho, \sigma) \approx \mathcal{F}_{\text{AFC}}^{(k)}(\rho, \sigma)$. This approximation finally allows us to estimate $\mathcal{F}_{\text{max}}(\rho, \sigma)$ efficiently. Indeed, $N_u^T N_M = O(4^{2k} N^3)$ measurements are sufficient to estimate $\mathcal{F}_{\text{AFC}}^{(k)}(\rho, \sigma)$ with finite statistical errors, while the systematic error $|\mathcal{F}_{\text{max}}(\rho, \sigma) - \mathcal{F}_{\text{AFC}}^{(k)}(\rho, \sigma)|$ decreases exponentially with $k$ [21, 66]. Now, we can test quantitatively the output $\sigma$ of the learning algorithm by estimating $\mathcal{F}_{\text{AFC}}^{(k)}(\rho, \sigma)$ for various $k \ge \ell$. Note that this method relies on the assumption that $\rho$ satisfies AFC for some finite $k$, and, in contrast to MPS verification [64, 65], the efficient verification of AFC is an open problem.
Experimental results— We now demonstrate the practical feasibility of our protocol by learning entangled quantum states from a partition of $N_{\text{tot}} = 96$ qubits of the IBM Brisbane superconducting quantum processing unit. We prepare entangled states by running the kicked Ising model quantum circuit at low depth $D$, defined as [70]

$$|\psi\rangle = \left( \prod_{j=1}^{N_{\text{tot}}-1} e^{i \frac{\pi}{4} Z_j Z_{j+1}} \prod_{j=1}^{N_{\text{tot}}} e^{-i \frac{\pi}{8} X_j} \right)^D |0\rangle^{\otimes N}, \quad (8)$$

where $X_j, Y_j, Z_j$ are Pauli operators. The preparation and the detection of the target state $|\psi\rangle$ are however affected by experimental noise, resulting in the observation of a mixed state $\rho$ which we analyze with our protocol.

For $D=1$ and $D=2$, we perform local randomized measurements with $N_M = 1024$ shots per basis and $N_u = N_u^L + N_u^T = 2048 \cdot D$ bases, and post-process the experimental data with the previously introduced learning algorithm. As a slight modification, we perform two-site updates (instead of one-site updates) of the tensors $M^{(j)}$, as this allows us to dynamically increase the bond dimension up to a value $\chi'$ [66]. In addition, we implement common randomized measurements in the estimation of Eq. (4) to reduce statistical errors [71] (although this step does not drastically change the results [66]). We always initialize the algorithm in the ansatz $\sigma = (\mathbb{I}/2)^{\otimes N}$ and perform $N_S = 20$ sweeps.

We report the results in Figs. 1 and 2. First, we benchmark the results of the algorithm by comparing the fidelities $F_{\text{max}}(\rho, \sigma)$ and $\langle\psi|\rho|\psi\rangle$. We observe that $F_{\text{max}}(\rho, \sigma) \gg \langle\psi|\rho|\psi\rangle$, meaning that the learned MPO $\sigma$ correctly captures the effect of experimental noise. For $N=96$ and $D=1,2$ we obtain fidelities of $F_{\text{max}}(\rho, \sigma) \sim 75\%$ and $\sim 50\%$ respectively, mainly due to statistical errors (as we similarly observe in numerical examples [66]). In a product state, the observed values would correspond to a fidelity per qubit ($F_{\text{max}}(\rho, \sigma))^{1/N} > 99\%$. To get a more meaningful insight, we study the magnetization $\langle Z_j \rangle$ and the two-body correlation functions $\langle Y_j Y_{j+d} \rangle$, see Fig. 2(c,d). We compute them for the states $\rho, \sigma$ and $|\psi\rangle$, using standard shadow estimation [27] from the testing set for the former. We observe that the expectation values obtained from $\rho$ and $\sigma$ are in good agreement, and both highlight non-ideal features of the experimental quantum state.

These observation corroborate that $\rho$ is well-described by a mixed state. To quantify this feature, we now analyze the second Rényi entropy $S_2(\rho) = -\log(\text{tr}[\rho^2])$ [72] of the full state, which can also be estimated using randomized measurements together with AFC [21]. In Figs. 1(d),2(b) we compare its AFC estimate $-\log(\mathcal{P}_{2,\text{AFC}}^{(k)})$ from the testing set with $S_2(\sigma)$, observing that they are in perfect agreement at large $k$. While previous experimental approaches were limited to $N =$

<figure><img src="image_4.png" /><figcaption>FIG. 3. Experimental quantum principal component analysis – (a) Running the DMRG algorithm on $H = -\sigma$ allows us to approximately find the principal component $|\psi_0^\sigma\rangle$, which approximates the principal component of the experimental state $\rho$. (b) Fidelity between $|\psi_0^\sigma\rangle$ and the target state $|\psi\rangle$ for the experiment at depth $D=1$, as a function of $\ell$ and $N_u$. (c) Expectation values of a local observable for depth $D=1$, compared between the learned MPO $\sigma$ and the error-mitigated state to the theory prediction. (d) Bipartite Von Neumann entropy as a function of the cut for $D=1,2$, compared between $|\psi_0^\sigma\rangle$ and $|\psi\rangle$.</figcaption></figure>

12, 13 qubits [30, 73], here we measure large-scale entropies to a high degree of accuracy, that correspond to exponentially small purities. By doing so we observe that the entropy of these states obeys a weak-volume law $S_2 \approx \alpha N$, where the coefficient $\alpha \approx 0.14-0.15$, compared to $\alpha=1$ for a fully mixed state, and the modest increase of the global entropy from $D=1$ to $D=2$ suggests that readout errors are dominant in the low-depth regime. Overall, these results demonstrate experimental tomography of $N=96$ qubits entangled quantum states, going beyond the state-of-the-art of $N=20$ qubits obtained for MPS [68] and Gibbs states [50].

**Application: quantum error mitigation**— We now leverage the access to a MPO description of $\rho$ to perform quantum error mitigation. Quantum error mitigation consists in partially removing the effect of errors in the estimation of properties from a quantum experiment via classical post-processing [74]. Conventional approaches are usually based on assuming specific noise models [75], or require additional experiments [76–79]. Our approach consists instead in quantum principal component analysis (QPCA) [80, 81], where we mitigate experimental noise by reconstructing the pure state $|\psi_0^\rho\rangle$ with the largest eigenvalue in the spectral decomposition $\rho = \sum_a \Lambda_a |\psi_a^\rho\rangle \langle\psi_a^\rho|$, where $\Lambda_a > \Lambda_{a+1}$ [82]. While this
idea is powerful, it has been implemented only approxi-
mately by virtual distillation [83–86].

Here we perform QPCA by running the DMRG algo-
rithm [56] on the MPO $H = -\sigma$, obtaining the state
$|\psi_0^\sigma\rangle$ as a MPS. Note that this approach does not re-
quire additional experiments nor data analysis. In Fig. 3
we illustrate the results of DMRG-QPCA on the experi-
mental data using $N = N_{tot}$ [87]. We obtain in this way
$|\psi_0^\sigma\rangle$ for $D = 1, 2$, and their corresponding eigenvalues
$\Lambda_0 = 9.11 \cdot 10^{-3}, 3.13 \cdot 10^{-3}$. Remarkably, we obtain high
fidelities between the MPSs $|\psi_0^\sigma\rangle$ and the target states
$|\psi\rangle$, even larger that 90% for $D = 1$, corroborating the
power of this error mitigation method. This is also re-
flected by expectation values of local observables, that
show a dramatic improvement from the MPO prediction
of $\sigma$. Finally, we leverage the MPS representation to
observe the growth of the bipartite Von Neumann en-
tropy $S_1 = \text{tr}_A[\rho' \log \rho']$ with the circuit depth, where
we take $\rho' = \text{tr}_A[|\psi_0^\sigma\rangle\langle\psi_0^\sigma|]$ and $A = [1, |A|]$, see Fig. 3(d).
These results experimentally demonstrate QPCA at large
scales, and provide a strong case for the application of our
learning protocol to quantum error mitigation.

Conclusions and outlook— Our protocol provides a quantum-to-classical converter [88], that compresses $N_u N_M$ randomized measurement outcomes from an experimental quantum state $\rho$ in a single classical object, the MPO $\sigma$. Having access to an MPO approximation of $\rho$ not only allows us to estimate several physical properties efficiently, i.e. without re-processing experimental data each time, but also to apply powerful error-mitigation based on tensor-network algorithms.

Our ability to combine quantum experiments with clas-
sical tensor-network methods opens up further possibil-
ities. For instance, MPO tomography can be associated
with MPS preparation algorithms which take as input
an MPO/MPS $\sigma$ represented classically and devise the
quantum circuit that efficiently forms the corresponding
state $\rho$ in an universal quantum computer [89, 90]. This
gives us the possibility to interact multiple times with a
quantum system, or even connect different experiments
to perform quantum computing beyond the possibilities
of each individual device.

An interesting use-case consists in quantum circuit cut-
ting: by saving an intermediate state of a quantum al-
gorithm as a MPO and performing QPCA on it, one
could resume the quantum algorithm by reloading the
 error-mitigated) state with a MPS preparation circuit.
Another concrete application consists in modular exper-
 iments where a fault-tolerant quantum computer is used
to analyze the state of a noisy quantum simulator. For
 instance, one could prepare a quantum state of interest
 in a quantum simulation experiment, and perform MPO
tomography; next, one could use a quantum circuit to
 prepare its MPS purification on a fault-tolerant quan-
tum computer, and finally apply a quantum algorithm
to measure the von Neumann entropy [91–93].

ACKNOWLEDGEMENTS

We acknowledge insightful discussions with Antoine
Browaeys, Mari Carmen Bañuls, Soonwon Choi, Thierry
Lahaye, Daniel Stilck-França, Georgios Styliaris, and
Xavier Waintal. The experimental data have been col-
lected using the Qiskit library [94], and post-processed
using the RandomMeas [95] and ITensor [96] libraries.
The work of M.V. and B.V. was funded by the
French National Research Agency via the JCJC project
QRand (ANR-20-CE47-0005), and via the research pro-
grams Plan France 2030 EPIQ (ANR-22-PETQ-0007),
QUBITAF (ANR-22-PETQ-0004) and HQI (ANR-22-
PNCQ-0002). We acknowledge the use of IBM Quantum
Credits for this work. The views expressed are those
of the authors, and do not reflect the official policy or
position of IBM or the IBM Quantum team. M.L. ac-
knowledges support by the Deutsche Forschungsgemein-
schaft (DFG, German Research Foundation) under Ger-
many’s Excellence Strategy – EXC-2111 – 390814868.
The work of C.L. was funded by the French National
Research Agency via the PRC project ESQuisses (ANR-
20-CE47-0014-01). J.I.C. acknowledges funding from the
Federal Ministry of Education and Research Germany
(BMBF) via the project FermiQP (13N15889). Work at
MPQ is part of the Munich Quantum Valley, which is
supported by the Bavarian state government with funds
from the Hightech Agenda Bayern Plus. P.Z. acknowl-
edges support by the European Union’s Horizon Europe
research and innovation programme under grant agree-
ment No 101113690 (PASQANS2). The work of L.P.
was funded by the European Union (ERC, QUANTHEM,
101114881). Views and opinions expressed are however
those of the author(s) only and do not necessarily re-
flect those of the European Union or the European Re-
search Council Executive Agency. Neither the European
Union nor the granting authority can be held respon-
sible for them. We acknowledge support by the Erwin
Schrödinger International Institute for Mathematics and
Physics (ESI).

<footnote>* matteo.votto@lpmmc.cnrs.fr</footnote>

<footnote>[1] R. Blatt and C. F. Roos, Nat. Phys. 8, 277 (2012).</footnote>

<footnote>[2] C. Gross and I. Bloch, Science 357, 995 (2017).</footnote>

<footnote>[3] J. Preskill, Quantum 2, 79 (2018).</footnote>

<footnote>[4] F. Schäfer, T. Fukuhara, S. Sugawa, Y. Takasu, and Y. Takahashi, Nat. Rev. Phys. 2, 411 (2020).</footnote>

<footnote>[5] M. Kjaergaard, M. E. Schwartz, J. Braumüller, P. Krantz, J. I.-J. Wang, S. Gustavsson, and W. D. Oliver, Annu. Rev. Condens. Matter Phys. 11, 369 (2020).</footnote>

<footnote>[6] E. Altman, K. R. Brown, G. Carleo, L. D. Carr, E. Demler, C. Chin, B. DeMarco, S. E. Economou, M. A. Eriksen, K.-M. C. Fu, M. Greiner, K. R. Hazzard, R. G.</footnote>
Hulet, A. J. Kollár, B. L. Lev, M. D. Lukin, R. Ma, X. Mi, S. Misra, C. Monroe, K. Murch, Z. Nazario, K.-K. Ni, A. C. Potter, P. Roushan, M. Saffman, M. Schleier-Smith, I. Siddiqi, R. Simmonds, M. Singh, I. Spielman, K. Temme, D. S. Weiss, J. Vučković, V. Vuletić, J. Ye, and M. Zwierlein, PRX Quantum 2, 017003 (2021).

[7] M. Morgado and S. Whitlock, AVS Quantum Sci. 3, 023501 (2021).

[8] C. Monroe, W. C. Campbell, L.-M. Duan, Z.-X. Gong, A. V. Gorshkov, P. W. Hess, R. Islam, K. Kim, N. M. Linke, G. Pagano, P. Richerme, C. Senko, and N. Y. Yao, Rev. Mod. Phys. 93, 025001 (2021).

[9] Y. Alexeev, D. Bacon, K. R. Brown, R. Calderbank, L. D. Carr, F. T. Chong, B. DeMarco, D. Englund, E. Farhi, B. Fefferman, A. V. Gorshkov, A. Houck, J. Kim, S. Kimmel, M. Lange, S. Lloyd, M. D. Lukin, D. Maslov, P. Maunz, C. Monroe, J. Preskill, M. Roetterler, M. J. Savage, and J. Thompson, PRX Quantum 2, 017001 (2021).

[10] E. Pelucchi, G. Fagas, I. Aharonovich, D. Englund, E. Figueroa, Q. Gong, H. Hannes, J. Liu, C.-Y. Lu, N. Matsuda, J.-W. Pan, F. Schreck, F. Sciarrino, C. Silberhorn, J. Wang, and K. D. Jöns, Nature Reviews Physics 4, 194 (2022).

[11] G. Burkard, T. D. Ladd, A. Pan, J. M. Nichol, and J. R. Petta, Rev. Mod. Phys. 95, 025003 (2023).

[12] A. Elben, B. Vermersch, C. F. Roos, and P. Zoller, Phys. Rev. A 99, 052323 (2019).

[13] A. Elben, S. T. Flammia, H.-Y. Huang, R. Kueng, J. Preskill, B. Vermersch, and P. Zoller, Nature Reviews Physics 5, 9-24 (2022).

[14] P. Cieśliński, S. Imai, J. Dziewior, O. Gühne, L. Knips, W. Laskowski, J. Meinecke, T. Paterek, and T. Vértesi, Physics Reports 1095, 1-48 (2024).

[15] T. Brydges, A. Elben, P. Jurcevic, B. Vermersch, C. Maier, B. P. Lanyon, P. Zoller, R. Blatt, and C. F. Roos, Science 364, 260-263 (2019).

[16] A. Elben, R. Kueng, H.-Y. R. Huang, R. van Bijnen, C. Kokail, M. Dalmonte, P. Calabrese, B. Kraus, J. Preskill, P. Zoller, and B. Vermersch, Phys. Rev. Lett. 125, 200501 (2020).

[17] A. Neven, J. Carrasco, V. Vitale, C. Kokail, A. Elben, M. Dalmonte, P. Calabrese, P. Zoller, B. Vermersch, R. Kueng, and B. Kraus, npj Quantum Information 7, 152 (2021).

[18] A. Rath, C. Branciard, A. Minguzzi, and B. Vermersch, Phys. Rev. Lett. 127, 260501 (2021).

[19] A. Rath, V. Vitale, S. Murciano, M. Votto, J. Dubail, R. Kueng, C. Branciard, P. Calabrese, and B. Vermersch, PRX Quantum 4, 010318 (2023).

[20] J. Carrasco, M. Votto, V. Vitale, C. Kokail, A. Neven, P. Zoller, B. Vermersch, and B. Kraus, Phys. Rev. A 109, 012422 (2024).

[21] B. Vermersch, M. Ljubotina, J. I. Cirac, P. Zoller, M. Serbyn, and L. Piroli, Phys. Rev. X 14, 031035 (2024).

[22] M. P. da Silva, O. Landon-Cardinal, and D. Poulin, Phys. Rev. Lett. 107, 210404 (2011).

[23] S. T. Flammia and Y.-K. Liu, Phys. Rev. Lett. 106, 230501 (2011).

[24] A. Elben, B. Vermersch, R. van Bijnen, C. Kokail, T. Brydges, C. Maier, M. K. Joshi, R. Blatt, C. F. Roos, and P. Zoller, Phys. Rev. Lett. 124, 010504 (2020).

[25] D. Zhu, Z. P. Cian, C. Noel, A. Risinger, D. Biswas, L. Egan, Y. Zhu, A. M. Green, C. H. Alderete, N. H. Nguyen, Q. Wang, A. Maksymov, Y. Nam, M. Cetina, N. M. Linke, M. Hafezi, and C. Monroe, Nature Comm. 13, 2315 (2022).

[26] H.-Y. Huang, J. Preskill, and M. Soleimanifar, (2024), arXiv:2404.07281 [quant-ph].

[27] H.-Y. Huang, R. Kueng, and J. Preskill, Nature Physics 16, 1050-1057 (2020).

[28] D. E. Koh and S. Grewal, Quantum 6, 776 (2022).

[29] S. Chen, W. Yu, P. Zeng, and S. T. Flammia, PRX Quantum 2, 030348 (2021).

[30] V. Vitale, A. Rath, P. Jurcevic, A. Elben, C. Branciard, and B. Vermersch, PRX Quantum 5, 030338 (2024).

[31] K. Noh, L. Jiang, and B. Fefferman, Quantum 4, 318 (2020).

[32] S. Cheng, C. Cao, C. Zhang, Y. Liu, S.-Y. Hou, P. Xu, and B. Zeng, Phys. Rev. Res. 3, 023005 (2021).

[33] F. Verstraete, J. J. García-Ripoll, and J. I. Cirac, Phys. Rev. Lett. 93, 207204 (2004).

[34] T. Kuwahara, A. M. Alhambra, and A. Anshu, Phys. Rev. X 11, 011047 (2021).

[35] J. Guth Jarkovský, A. Molnár, N. Schuch, and J. I. Cirac, PRX Quantum 1, 010304 (2020).

[36] T. Baumgratz, A. Nüßeler, M. Cramer, and M. B. Plenio, New Journal of Physics 15, 125004 (2013).

[37] A. Lidiak, C. Jameson, Z. Qin, G. Tang, M. B. Wakin, Z. Zhu, and Z. Gong, (2022), arXiv:2207.06397 [quant-ph].

[38] G. Torlai, C. J. Wood, A. Acharya, G. Carleo, J. Carrasquilla, and L. Aolita, Nature Communications 14, 2858 (2023).

[39] Y. Guo and S. Yang, Comm. Phys. 7, 322 (2024).

[40] C. Jameson, Z. Qin, A. Goldar, M. B. Wakin, Z. Zhu, and Z. Gong, (2024), arXiv:2408.07115 [quant-ph].

[41] T. Baumgratz, D. Gross, M. Cramer, and M. B. Plenio, Phys. Rev. Lett. 111, 020401 (2013).

[42] M. Holzäpfel, M. Cramer, N. Datta, and M. B. Plenio, J. Math. Phys. 59, 042201 (2018).

[43] Z. Qin, C. Jameson, Z. Gong, M. B. Wakin, and Z. Zhu, IEEE Transactions on Information Theory 70, 5030-5056 (2024).

[44] Z. Qin, C. Jameson, A. Goldar, M. B. Wakin, Z. Gong, and Z. Zhu, (2024), arXiv:2410.02583 [quant-ph].

[45] M. Fanizza, N. Galke, J. Lumbreras, C. Rouzé, and A. Winter, (2024), arXiv:2312.07516 [quant-ph].

[46] Z. Qin, J. M. Lukens, B. T. Kirby, and Z. Zhu, (2025), arXiv:2501.03144 [quant-ph].

[47] B. Swingle and I. H. Kim, Phys. Rev. Lett. 113, 260501 (2014).

[48] A. Anshu, S. Arunachalam, T. Kuwahara, and M. Soleimanifar, Nature Physics 17, 931-935 (2021).

[49] C. Kokail, R. van Bijnen, A. Elben, B. Vermersch, and P. Zoller, Nature Physics 17, 936-942 (2021).

[50] M. K. Joshi, C. Kokail, R. van Bijnen, F. Kranzl, T. V. Zache, R. Blatt, C. F. Roos, and P. Zoller, Nature 624, 539-544 (2023).

[51] C. Rouzé and D. Stilck França, Quantum 8, 1319 (2024).

[52] J. Haah, R. Kothari, and E. Tang, Nature Physics 20, 1027-1031 (2024).

[53] A. Bakshi, A. Liu, A. Moitra, and E. Tang, Proceedings of the 56th Annual ACM Symposium on Theory of Computing STOC 2024, 1470-1477 (2024).

[54] Ángela Capel, Álvaro M. Alhambra, P. Gondolf, A. R. de Alarcón, and S. O. Scalet, (2024), arXiv:2402.18500
[quant-ph].

[55] A. Anshu and S. Arunachalam, Nature Reviews Physics 6, 59 (2024).

[56] U. Schollwöck, Ann. Phys. 326, 96–192 (2011).

[57] E. Stoudenmire and D. J. Schwab, Advances in Neural Information Processing Systems 29 (2016).

[58] Z.-Y. Han, J. Wang, H. Fan, L. Wang, and P. Zhang, Phys. Rev. X 8, 031012 (2018).

[59] T. Ayral, T. Louvet, Y. Zhou, C. Lambert, E. M. Stoudenmire, and X. Waintal, PRX Quantum 4, 020304 (2023).

[60] Y.-C. Liang, Y.-H. Yeh, P. E. M. F. Mendonça, R. Y. Teh, M. D. Reid, and P. D. Drummond, Reports on Progress in Physics 82, 076001 (2019).

[61] S. Sang and T. H. Hsieh, Phys. Rev. Lett. 134, 070403 (2025).

[62] C. Zhang, Y. Xu, J.-H. Zhang, C. Xu, Z. Bi, and Z.-X. Luo, Phys. Rev. B 111, 115137 (2025).

[63] J. I. Cirac, D. Pérez-García, N. Schuch, and F. Verstraete, Reviews of Modern Physics 93, 045003 (2021).

[64] M. Cramer, M. B. Plenio, S. T. Flammia, R. Somma, D. Gross, S. D. Bartlett, O. Landon-Cardinal, D. Poulin, and Y.-K. Liu, Nature Comm. 1, 149 (2010).

[65] B. P. Lanyon, C. Maier, M. Holzäpfel, T. Baumgratz, C. Hempel, P. Jurcevic, I. Dhand, A. S. Buyskikh, A. J. Daley, M. Cramer, M. B. Plenio, R. Blatt, and C. F. Roos, Nature Physics 13, 1158–1162 (2017).

[66] See Supplemental Material for further details.

[67] K. Koor, Y. Qiu, L. C. Kwek, and P. Rebentrost, (2023), arXiv:2312.04858 [quant-ph].

[68] M. K. Kurmapu, V. Tiunova, E. Tiunov, M. Ringbauer, C. Maier, R. Blatt, T. Monz, A. K. Fedorov, and A. Lvovsky, PRX Quantum 4, 040345 (2023).

[69] F. Altomare and M. Campiti, Korovkin-type Approximation Theory and Its Applications (De Gruyter, Berlin, New York, 1994).

[70] Y. Kim, A. Eddins, S. Anand, K. X. Wei, E. van den Berg, S. Rosenblatt, H. Nayfeh, Y. Wu, M. Zaletel, K. Temme, and A. Kandala, Nature 618, 500 (2023).

[71] B. Vermersch, A. Rath, B. Sundar, C. Branciard, J. Preskill, and A. Elben, PRX Quantum 5, 010352 (2024).

[72] We consider the logarithm to be in base 2, such that $S_2((\mathbb{I}/2)^{\otimes N}) = N$.

[73] T. I. Andersen, N. Astrakhantsev, A. H. Karamlou, J. Berndtsson, J. Motruk, A. Szasz, J. A. Gross, A. Schuckert, T. Westerhout, Y. Zhang, et al., Nature 638, 79 (2025).

[74] Z. Cai, R. Babbush, S. C. Benjamin, S. Endo, W. J. Huggins, Y. Li, J. R. McClean, and T. E. O'Brien, Reviews of Modern Physics 95, 045005 (2023).

[75] E. van den Berg, Z. K. Minev, A. Kandala, and K. Temme, Nature Physics 19, 1116–1121 (2023).

[76] T. Giurgica-Tiron, Y. Hindy, R. LaRose, A. Mari, and W. J. Zeng, 2020 IEEE International Conference on Quantum Computing and Engineering (QCE), 306 (2020).

[77] Y. Guo and S. Yang, PRX Quantum 3, 040313 (2022).

[78] S. Filippov, M. Leahy, M. A. C. Rossi, and G. García-Pérez, (2023), arXiv:2307.11740 [quant-ph].

[79] S. Mangini, M. Cattaneo, D. Cavalcanti, S. Filippov, M. A. C. Rossi, and G. García-Pérez, Phys. Rev. Research 6, 033217 (2024).

[80] S. Lloyd, M. Mohseni, and P. Rebentrost, Nature Physics 10, 631–633 (2014).

[81] H.-Y. Huang, M. Broughton, J. Cotler, S. Chen, J. Li, M. Mohseni, H. Neven, R. Babbush, R. Kueng, J. Preskill, and J. R. McClean, Science 376, 1182–1186 (2022).

[82] B. Koczor, New Journal of Physics 23, 123047 (2021).

[83] J. Cotler, S. Choi, A. Lukin, H. Gharibyan, T. Grover, M. E. Tai, M. Rispoli, R. Schittko, P. M. Preiss, A. M. Kaufman, M. Greiner, H. Pichler, and P. Hayden, Phys. Rev. X 9, 031013 (2019).

[84] W. J. Huggins, S. McArdle, T. E. O'Brien, J. Lee, N. C. Rubin, S. Boixo, K. B. Whaley, R. Babbush, and J. R. McClean, Phys. Rev. X 11, 041036 (2021).

[85] A. Seif, Z.-P. Cian, S. Zhou, S. Chen, and L. Jiang, PRX Quantum 4, 010303 (2023).

[86] H. Hakushima, S. Endo, K. Yamamoto, Y. Matsuzaki, and N. Yoshioka, Phys. Rev. Lett. 133, 080601 (2024).

[87] We initialize DMRG with $|\psi_0\rangle \propto |\psi_{RMPS}\rangle + \epsilon|0\rangle^{\otimes N}$ where $|\psi_{RMPS}\rangle$ is a random MPS with small bond dimension, and $\epsilon \gtrsim 10^{-2}$. For $\epsilon = 0$, we numerically observe that the DMRG algorithm does not update the random MPS.

[88] H.-Y. Huang, R. Kueng, G. Torlai, V. V. Albert, and J. Preskill, Science 377, eabk3333 (2022).

[89] D. Malz, G. Styliaris, Z.-Y. Wei, and J. I. Cirac, Phys. Rev. Lett. 132, 040404 (2024).

[90] S. Fomichev, K. Hejazi, M. S. Zini, M. Kiser, J. Fraxanet, P. A. M. Casares, A. Delgado, J. Huh, A.-C. Voigt, J. E. Mueller, and J. M. Arrazola, PRX Quantum 5, 040339 (2024).

[91] T. Gur, M.-H. Hsieh, and S. Subramanian, (2021), arXiv:2111.11139 [quant-ph].

[92] A. Luongo and C. Shao, (2024), arXiv:2011.06475 [quant-ph].

[93] V. Giovannetti, S. Lloyd, and L. Maccone, (2025), arXiv:2504.11049 [quant-ph].

[94] A. Javadi-Abhari, M. Treinish, K. Krsulich, C. J. Wood, J. Lishman, J. Gacon, S. Martiel, P. D. Nation, L. S. Bishop, A. W. Cross, B. R. Johnson, and J. M. Gambetta, (2024), arXiv:2405.08810 [quant-ph].

[95] A. Elben and B. Vermersch, Randommeas: The randomized measurement toolbox in julia (2025).

[96] M. Fishman, S. R. White, and E. M. Stoudenmire, SciPost Phys. Codebases, 4 (2022).
# Supplemental Materials: Learning mixed quantum states in large-scale experiments

## CONTENTS

<table><tr><td>A. Learning algorithm</td><td>1</td></tr><tr><td>1. Exact tensor optimization</td><td>2</td></tr><tr><td>2. Approximate tensor optimization</td><td>2</td></tr><tr><td>3. Convergence guarantees and two-site updates</td><td>4</td></tr><tr><td>B. Local approximation of gradient terms</td><td>4</td></tr><tr><td>1. Translation-invariant matrix-product operators and transfer operators</td><td>4</td></tr><tr><td>2. Proof for local approximation of gradient terms</td><td>6</td></tr><tr><td>Spectral decompositions</td><td>6</td></tr><tr><td>Upper bounds</td><td>7</td></tr><tr><td>Final result</td><td>9</td></tr><tr><td>C. Statistical errors in the learning algorithm</td><td>9</td></tr><tr><td>1. Local randomized measurements and classical shadows</td><td>9</td></tr><tr><td>2. Statistical errors in the optimization of a single tensor</td><td>10</td></tr><tr><td>D. Scalable estimation of mixed-state fidelities</td><td>11</td></tr><tr><td>1. Mixed-state fidelities from randomized measurements</td><td>12</td></tr><tr><td>2. Approximate factorization conditions for fidelities</td><td>12</td></tr><tr><td>3. Approximate factorization conditions for the overlap: proof for TIMPDOs</td><td>13</td></tr><tr><td>4. Statistical errors in AFC fidelity estimation: numerical results</td><td>15</td></tr><tr><td>E. Numerical results for the learning algorithm</td><td>16</td></tr><tr><td>1. Approximate local updates</td><td>16</td></tr><tr><td>2. Statistical errors</td><td>17</td></tr><tr><td>F. Common randomized measurements</td><td>18</td></tr><tr><td>1. Common randomized measurements and tensor updates</td><td>18</td></tr><tr><td>2. Experimental results</td><td>19</td></tr></table>

## Appendix A: Learning algorithm

In this section we derive the optimization step of the learning algorithm described in the main text. In particular, given an experimental state $\rho$ to be learned and a MPO ansatz $\sigma$, defined as

$$ \sigma = \sum_{\{s_j\}, \{s'_j\}} M_{s_1, s'_1}^{(1)} M_{s_2, s'_2}^{(2)} \dots M_{s_N, s'_N}^{(N)} |\{s_j\}\rangle \langle\{s'_j\}|, \quad (A1) $$

we prove that the optimization step for the tensor $M^{(j)}$ proposed in the main text approximately maximizes the geometric-mean fidelity [24, 60]

$$ \mathcal{F}_{\text{GM}}(\rho, \sigma) = \frac{\text{tr}[\rho\sigma]}{\sqrt{\text{tr}[\rho^2]\text{tr}[\sigma^2]}}, \quad (A2) $$

that is differentiable and can be efficiently computed numerically. This is proven in two steps. First, we derive the gradient of $\mathcal{F}_{\text{GM}}(\rho, \sigma)$ with respect to a tensor $M^{(j)}$ of the MPO $\sigma$, and show that the optimization step can be carried out exactly. This step is similar to the derivation of the optimization step in the density matrix renormalization group (DMRG) algorithm [56] and other sequential tensor network optimization schemes [57–59]. Then, we show how under the assumptions that (a) $\rho$ is a MPO with bond dimension $\chi$ and (b) the states $\rho$ and $\rho' \propto \rho^2$ have finite correlation lengths, the optimization can be carried out approximately. This is due to an approximation of the two terms in the gradient, which we prove in detail in App. B. Finally, we explain how to choose the parameters of the learning algorithm to guarantee that the optimization step leads to an accurate solution.
<figure><img src="image_5.png" /><figcaption>FIG. S1. Diagrammatic representation of $C_{\sigma}^{(j)} M^{(j)}$ and $C_{\rho}^{(j)}$ (defined in Eq. (A5)) in terms of $\rho$ and the local tensor of $\sigma$.</figcaption></figure>

## 1. Exact tensor optimization

We start by deriving the equation satisfied by an individual tensor $M^{(j)}$ of the MPO $\sigma$ that maximizes $\mathcal{F}_{\text{GM}}(\rho, \sigma)$. As a first step, we compute the logarithmic gradient of the GM fidelity with respect to a local tensor

$$G^{(j)} = \frac{\partial_{M^{(j)}} \mathcal{F}_{\text{GM}}(\rho, \sigma)}{\mathcal{F}_{\text{GM}}(\rho, \sigma)} = \frac{\text{tr}[\rho \partial_{M^{(j)}} \sigma]}{\text{tr}[\rho \sigma]} - \frac{\text{tr}[\sigma \partial_{M^{(j)}} \sigma]}{\text{tr}[\sigma^2]}, \quad (A3)$$

where we have used Wirtinger derivative rules [67] to obtain that

$$\partial_{M^{(j)}} \text{tr}[\rho\sigma] = \text{tr}[\rho\partial_{M^{(j)}}\sigma], \quad \partial_{M^{(j)}} \text{tr}[\sigma^2] = 2\text{tr}[\sigma\partial_{M^{(j)}}\sigma]. \quad (A4)$$

Assuming that the GM fidelity has only global maxima as a function of $M^{(j)}$, the optimal $M^{(j)}$ satisfies $G^{(j)} = 0$. We show now that it is possible to obtain the local tensor $M^{(j)}$ that satisfies this condition exactly by solving a linear system. Let us define the terms

$$C_{\rho}^{(j)} = \text{tr}[\rho \partial_{M^{(j)}} \sigma], \quad C_{\sigma}^{(j)} M^{(j)} = \text{tr}[\sigma \partial_{M^{(j)}} \sigma] \quad (A5)$$

such that both $C_{\rho}^{(j)}$ and $C_{\sigma}^{(j)}$ do not depend on $M^{(j)}$ (see Fig. S1 for their diagrammatic notation). In terms of these objects, the solution to $G^{(j)} = 0$ is given by

$$C_{\sigma}^{(j)} M^{(j)} = \frac{\text{tr}[\sigma^2]}{\text{tr}[\rho\sigma]} C_{\rho}^{(j)}. \quad (A6)$$

Using that

$$\text{tr}[\rho\sigma] = C_{\rho}^{(j)} \cdot M^{(j)}, \quad \text{tr}[\sigma^2] = C_{\sigma}^{(j)} M^{(j)} \cdot M^{(j)} \quad (A7)$$

it is possible to see that Eq. (A6) is solved by any tensor $M^{(j)}$ of the form

$$M^{(j)} = \lambda C_{\sigma}^{(j)-1} C_{\rho}^{(j)} \quad (A8)$$

for any coefficient $\lambda \neq 0$. We initially choose $\lambda = 1$, which corresponds to imposing

$$\text{tr}[\sigma \partial_{M^{(j)}} \sigma] = \text{tr}[\rho \partial_{M^{(j)}} \sigma], \quad (A9)$$

and then normalize the tensor $M^{(j)}$ such that $\text{tr}[\sigma] = 1$. We express the resulting update rule for the tensor $M^{(j)}$ in the following way

$$M^{(j)} \leftarrow C_{\sigma}^{(j)-1} C_{\rho}^{(j)}, \quad M^{(j)} \leftarrow \frac{1}{\text{tr}[\sigma]} M^{(j)}. \quad (A10)$$

where the first update rule consists in solving a linear system of $4\chi'^2$ equations, where $\chi'$ is the bond dimension of the MPO $\sigma$, and the second one normalizes $\sigma$. This update rule also imposes the Hermiticity of $\sigma$, but not its positivity. Finally, note that such an optimization scheme would require us to estimate the right-hand side of Eq. (A9) from experimental data, which requires $O(2^N)$ measurements using classical shadow framework [27].

## 2. Approximate tensor optimization

We now show that it is possible to approximately optimize the tensor $M^{(j)}$ as

$$\text{tr}[\sigma_{I_j} \partial_{M^{(j)}} \sigma_{I_j}] = \text{tr}[\rho_{I_j} \partial_{M^{(j)}} \sigma_{I_j}] \quad (A11)$$

up to a normalization factor. This can be proven rigorously under the assumptions that
<figure><img src="image_6.png" /><figcaption>FIG. S2. Graphical representation of the condition Eq. (A11) for the approximate update rule for $N = 10$, $\ell = 2$, and $j = 5$. We are considering a right sweep, where the black tensors of $\sigma$ have been already updated, while the gray tensors not yet.</figcaption></figure>

1. $\rho$ is a translation-invariant matrix-product density operator (MPDO, i.e. a quantum state described by a MPO) with finite bond dimension $\chi$ and Rényi-Markov length $\xi_{\rho}^{(2)}$;

2. the algorithm is close to convergence, so that we can assume that $\sigma$ is a translation-invariant MPDO with bond dimension $\chi$ as well.

Note that the condition on the Rényi-Markov length is satisfied by typical MPDOs [21, 63]. Our starting point is the following lemma: given $\rho$ and $\sigma$ two translation-invariant MPDOs with bond dimension $\chi$, and given that the correlation length $\xi_{\rho\sigma}$ of the operator $\rho\sigma$ is finite, under mild technical assumptions the following approximation is valid (see App. B for a rigorous proof)

$$C_{\rho}^{(j)} = a_{\rho}(C_{\rho I_j}^{(j)} + \varepsilon_{\rho}), \qquad (A12)$$

with

$$a_{\rho} = \frac{\operatorname{tr}[\rho \tilde{\sigma}]}{\operatorname{tr}[\rho_{I_j} \tilde{\sigma}_{I_j}]}, \qquad (A13)$$

where $I_j = [j - \ell, j + \ell]$, and $\tilde{\sigma}$ is the MPO obtained by substituting the tensor $M^{(j)}$ in $\sigma$ with any tensor $\tilde{M}^{(j)}$ such that $\operatorname{tr}[\rho\tilde{\sigma}] \neq 0$, and

$$\|\varepsilon_{\rho}\|_2^2 \le \operatorname{tr}[\rho_{I_j} \tilde{\sigma}_{I_j}]^2 K \chi^2 e^{-\ell/\xi_{\rho\sigma}}, \qquad (A14)$$

where $K = O(1)$ is a constant. Analogously, if the Rényi-Markov length $\xi_{\sigma}^{(2)}$ is also finite, the same is valid upon the substitution $\rho \to \sigma$, yielding

$$C_{\sigma}^{(j)} \tilde{M}^{(j)} = a_{\sigma}(C_{\sigma I_j}^{(j)} \tilde{M}^{(j)} + \varepsilon_{\sigma}), \qquad (A15)$$

for any tensor $\tilde{M}^{(j)}$ such that $\operatorname{tr}[\tilde{\sigma}^2] \neq 0$

Let us now take the equation satisfied by the optimal tensor $M^{(j)}$ [Eq. (A6)]

$$C_{\sigma}^{(j)} M^{(j)} = \frac{M^{(j)} \cdot C_{\sigma I_j}^{(j)} M^{(j)}}{M^{(j)} \cdot C_{\rho I_j}^{(j)}} C_{\rho}^{(j)}, \qquad (A16)$$

it is now possible to expand it by applying the approximations Eqs. (A12), (A15), obtaining

$$C_{\sigma I_j}^{(j)} M^{(j)} + \varepsilon_{\sigma} = \frac{M^{(j)} \cdot C_{\sigma I_j}^{(j)} M^{(j)} + M^{(j)} \cdot \varepsilon_{\sigma}}{M^{(j)} \cdot C_{\rho I_j}^{(j)} + M^{(j)} \cdot \varepsilon_{\rho}} (C_{\rho I_j}^{(j)} + \varepsilon_{\rho}). \qquad (A17)$$

By taking $\|\varepsilon_{\rho}\|_2, \|\varepsilon_{\sigma}\|_2 \to 0$, we can expand perturbatively Eq. (A17) in $\varepsilon_{\rho}, \varepsilon_{\sigma}$, obtaining at the leading order

$$C_{\sigma I_j}^{(j)} M_0^{(j)} = \frac{M_0^{(j)} \cdot C_{\sigma I_j}^{(j)} M_0^{(j)}}{M_0^{(j)} \cdot C_{\rho I_j}^{(j)}} C_{\rho I_j}^{(j)}, \qquad (A18)$$
where $||M^{(j)} - M_0^{(j)}||_2 = O(\epsilon)$ with $\epsilon = \max(||\epsilon_\rho||_2, ||\epsilon_\sigma||_2)$. If we assume that, close to convergence, $\xi_{\rho\sigma}, \xi_\sigma^{(2)} \simeq \xi_\rho^{(2)}$, then for $\ell > O(\xi_\rho^{(2)})$, the optimal tensor $M^{(j)}$ can be approximated by the tensor $M_0^{(j)}$. By repeating the same steps as in the previous section, it is possible to express $M_0^{(j)}$ as

$$ M_0^{(j)} = \lambda C_{\sigma I_j}^{(j)-1} C_{\rho I_j}^{(j)} \qquad (A19) $$

where the constant $\lambda$ is fixed by the normalization of $\sigma$. This implies that we can first get $M_0^{(j)}$ for $\lambda = 1$ by solving Eq. (A11), and then fix $\lambda = \frac{1}{\text{tr}[\sigma]}$. Therefore, the learning algorithm with approximate update is controlled by two parameters $(\ell, \chi')$, where the former is related to the Rényi-Markov length $\xi_\rho^{(2)}$, while the latter to the amount of entanglement encoded in $\rho$ [19, 31, 35].

### 3. Convergence guarantees and two-site updates

We now discuss for what choice of $(\ell, \chi')$ the approximate optimization step described above is guaranteed to converge. As described in the main text, each tensor $M^{(j)}$ is optimized individually by solving Eq. (A11), which amounts to a linear system in $M^{(j)}$. This linear system has $4\chi'^2$ variables, and at most $4(4^\ell)^2$ independent equations. This can be seen by performing a sequential singular-value decomposition (SVD) of $\sigma_{I_j}$, which results in a MPO with maximal bond dimension $4^\ell$. This implies that the linear system solver is guaranteed to converge to a unique solution only if

$$ \chi' \le 4^\ell. \qquad (A20) $$

When this is not the case, some parameters of the local tensor $M^{(j)}$ may not be updated in a way which minimizes the global fidelity, resulting in the propagation of errors.

Another possibility consists in performing two-site updates, where at each step we optimize a tensor $\theta^{(j,j+1)} = M^{(j)}M^{(j+1)}$. In this way, we do not have to fix the bond dimension of $\sigma$ from the beginning, but we can increase it dynamically, similarly to what is done in DMRG [56]. However, even if the bond dimension is initialized to $\chi' = 1$, it can grow until $\chi' = 4^{\ell+1}$, as the left and right bond dimensions of $\theta^{(j,j+1)}$ can grow up to $4^\ell$. This means that we need to introduce a cut-off to the bond dimension $\chi'_{\max} \le 4^\ell$ when performing the SVD of $\theta^{(j,j+1)} = M^{(j)}M^{(j+1)}$ after a two-site update. Otherwise, certain local updates might involve solving an under-determined linear system, leading to a possible propagation of uncontrolled errors.

## Appendix B: Local approximation of gradient terms

In this section we provide a rigorous proof for the lemma Eq. (A14) in the case of translation invariant MP-DOs, establishing the validity of the local approximate updates with sufficiently large $\ell$. Before doing so, we define translation-invariant matrix-product density operators (TIMPDOs), and we relate their properties to their transfer operators. In particular, we rigorously define the correlation length of a MPDO. The same language will be used in Sec. D to prove approximation factorization conditions for the fidelity.

### 1. Translation-invariant matrix-product operators and transfer operators

A quantum state $\rho$ is a translation-invariant MPDO (TIMPDO) if it can be written as in Eq. (A1) and $M^{(j)} = M$, i.e. when all the tensors are equal. To have a well-defined translation invariance, we consider MPDOs with periodic boundary conditions, i.e.

$$ \rho = \sum_{\{s_i\}, \{s'_i\}} \text{tr}[M'_{s_1, s'_1} M'_{s_2, s'_2} \dots M'_{s_N, s'_N}] |s_1, \dots, s_N\rangle \langle s'_1, \dots, s'_N|, \qquad (B1) $$

and similarly for $\sigma$

$$ \sigma = \sum_{\{s_j\}, \{s'_j\}} \text{tr}[M_{s_1, s'_1} M_{s_2, s'_2} \dots M_{s_N, s'_N}] |s_1, \dots, s_N\rangle \langle s'_1, \dots, s'_N|. \qquad (B2) $$
<figure><img src="image_7.png" /><figcaption>FIG. S3. Graphical representation of the transfer operators $T_{\rho}^{(n)}$, $T_{\sigma}^{(n)}$ and $T_{\rho\sigma}$, defined in Eqs. (B4), (B5).</figcaption></figure>

For two TIMPDOs $\rho$ and $\sigma$ we can define the transfer operators (see Fig. S3 for a graphical representation)

$$T_{\rho}^{(n)} = \sum_{\{s_i^{(k)}\}} M'_{s_i^{(1)}, s_i^{(2)}} M'_{s_i^{(2)}, s_i^{(3)}} \dots M'_{s_i^{(n)}, s_i^{(1)}}, \quad T_{\sigma}^{(n)} = \sum_{\{s_i^{(k)}\}} M_{s_i^{(1)}, s_i^{(2)}} M_{s_i^{(2)}, s_i^{(3)}} \dots M_{s_i^{(n)}, s_i^{(1)}}, \quad (B3)$$

from which we can obtain the moments

$$\mathrm{tr} \left[ \left( T_{\rho}^{(n)} \right)^{N} \right] = \mathrm{tr}[\rho^{n}], \quad \mathrm{tr} \left[ \left( T_{\sigma}^{(n)} \right)^{N} \right] = \mathrm{tr}[\sigma^{n}]. \quad (B4)$$

Additionally, we define the mixed transfer operator

$$T_{\rho\sigma} = \sum_{s_i, s'_i} M_{s_i, s'_i}^{(i)} M_{s'_i, s_i}^{(i)}, \quad (B5)$$

which analogously yields

$$\mathrm{tr}[\rho\sigma] = \mathrm{tr}[T_{\rho\sigma}^N]. \quad (B6)$$

In what follows, we will focus on the spectral properties of the transfer operators, as obtained from the decompositions

$$T_{\rho}^{(1)} = \sum_{a} \lambda_a |l_a^{\rho}| \langle r_a^{\rho} |, \quad T_{\sigma}^{(1)} = \sum_{a} \mu_a |l_a^{\sigma}| \langle r_a^{\sigma} |, \quad T_{\rho\sigma} = \sum_{a} \nu_a |l_a^{\rho\sigma}| \langle r_a^{\rho\sigma} |, \quad (B7)$$

where the eigenvalues are ordered in absolute value (e.g. $|\lambda_i| \ge |\lambda_{i+1}|$), and we have the orthonormality conditions for left- and right-eigenvectors $\langle l_a^O | r_b^O \rangle = \langle r_a^O | l_b^O \rangle = \delta_{ab}$. In case the spectrum of the transfer operator $T_{\rho}^{(1)}$ is gapped, i.e. $|\lambda_1/\lambda_0| < 1$, the connected two-point correlation functions of $\rho$ will decay exponentially with the distance [63], inducing the following definition of correlation length

$$|\lambda_1/\lambda_0| = e^{-1/\xi_{\rho}^{(1)}}. \quad (B8)$$

Additionally, when a higher order transfer operator $T_{\rho}^{(n)}$ is gapped, it has been shown that the global entanglement properties (like the Rényi entropies) are approximately determined by subsystems of a finite length $O(\xi_{\rho}^{(n)})$, defined in an analogous way [21]. This is also known as approximate conditional independence in terms of Rényi entropies [54], and $\xi_{\rho}^{(n)}$ is called Rényi-Markov length [61, 62].

In what follows, we assume the correlation lengths $\xi_{\rho}^{(1)}$, $\xi_{\rho}^{(2)}$ to be finite. Additionally, we also assume the mixed transfer $T_{\rho\sigma}$ operator to be gapped, i.e. $\xi_{\rho\sigma} = \frac{1}{\log|\nu_0/\nu_1|} = O(1)$. While this last condition is much less explored, we believe it to be true for a generic pair of TIMPDOs with finite correlation lengths. Practically, towards the convergence of the learning algorithm $\rho \sim \sigma$, we expect $\xi_{\rho\sigma} \sim \xi_{\rho}^{(2)} = O(1)$. We make the further technical assumptions that the complex scalars

$$\Gamma_a^{bc} = \langle r_a^{\rho\sigma} | l_b^{\rho} l_c^{\sigma} \rangle, \quad \tilde{\Gamma}_{bc}^a = \langle r_b^{\rho} r_c^{\sigma} | l_a^{\rho\sigma} \rangle \quad (B9)$$

are both non-zero, i.e. $\Gamma_a^{bc}, \tilde{\Gamma}_{bc}^a \neq 0$. We also believe this assumption to be true in the generic case [21].
## 2. Proof for local approximation of gradient terms

In this section we use the transfer operator formalism to prove Eq. (A14). Let us now take $\tilde{\sigma}$ to be the MPO obtained by substituting the tensor $M^{(j)}$ with an arbitrary tensor $\tilde{M}^{(j)}$ such that $\text{tr}[\rho\tilde{\sigma}] \neq 0$, where $\sigma$ is a TIMPDO with bond dimension $\chi$. Furthermore, let us assume $\rho$ is a TIMPDO with bond dimension $\chi$ as well. Then, under the assumptions that $\xi_{\rho\sigma} = O(1)$ and $\Gamma_{a}^{bc}, \tilde{\Gamma}_{bc}^{a} \neq 0$, we have

$$ \left\| \frac{C_{\rho I_j}^{(j)}}{\text{tr}[\rho I_j \tilde{\sigma}_{I_j}]} - \frac{C_{\rho}^{(j)}}{\text{tr}[\rho \tilde{\sigma}]} \right\|_2^2 \leq K \chi^2 e^{-\ell/\xi_{\rho\sigma}}, \qquad (B10) $$

where $C_{\rho}^{(j)} = \text{tr}[\rho\partial_{M^{(j)}}\sigma]$, $C_{\rho I_j}^{(j)} = \text{tr}[\rho I_j \partial_{M^{(j)}}\sigma_{I_j}]$ and $K = O(1)$ is a constant, and the norm-2 is the Hilbert-Schmidt norm

$$ \|C_{\rho}^{(j)}\|_{2}^{2} = \sum_{a,a',s,s'} [\tilde{C}_{\rho}^{(j)}]_{(a,a'),(s,s')} [C_{\rho}^{(j)}]_{(a,a'),(s,s')}, \qquad (B11) $$

where $a, a'$ ($s, s'$) are the two link (site) indices, and $\tilde{C}_{\rho}^{(j)}$ is the complex conjugate of $C_{\rho}^{(j)}$. A proof similar to the following can be used to show that

$$ \left\| \frac{C_{\sigma I_j}^{(j)} \tilde{M}^{(j)}}{\text{tr}[\tilde{\sigma}_{I_j}^2]} - \frac{C_{\sigma}^{(j)} \tilde{M}^{(j)}}{\text{tr}[\tilde{\sigma}^2]} \right\|_2^2 \leq K' \chi^2 e^{-\ell/\xi_{\sigma}^{(2)}}, \qquad (B12) $$

where $\tilde{M}^{(j)}$ is any tensor such that $\text{tr}[\tilde{\sigma}^2]$, and $K'$ is a constant.

As a first step, let us expand the norm-2 distance as follows

$$ \left\| \frac{C_{\rho}^{(j)}}{\text{tr}[\rho \tilde{\sigma}]} - \frac{C_{\rho I_j}^{(j)}}{\text{tr}[\rho I_j \tilde{\sigma}_{I_j}]} \right\|_2^2 = \frac{\tilde{C}_{\rho}^{(j)} \cdot C_{\rho}^{(j)}}{\text{tr}[\rho \tilde{\sigma}]^2} - \frac{2\text{Re}(\tilde{C}_{\rho}^{(j)} \cdot C_{\rho I_j}^{(j)})}{\text{tr}[\rho \tilde{\sigma}] \text{tr}[\rho I_j \tilde{\sigma}_{I_j}]} + \frac{\tilde{C}_{\rho I_j}^{(j)} \cdot C_{\rho I_j}^{(j)}}{\text{tr}[\rho I_j \tilde{\sigma}_{I_j}]^2}. \qquad (B13) $$

By expressing every factor in Eq. (B13) in terms of the spectral decompositions of the transfer operators, we obtain that each one of the three terms can be further decomposed in a term constant in $\ell$, and a term dependent on it. As we now prove in two steps, the constant contribution is equal among the three terms and gets canceled out in the algebraic sum, while the other contribution admits an upper bound proportional to $O(\chi^2 e^{-\ell/\xi_{\rho\sigma}})$.

### Spectral decompositions

Using the spectral decompositions of the transfer operators, we obtain the following expressions for the overlaps

$$ \text{tr}[\rho\tilde{\sigma}] = \sum_{a=0}^{\chi^2-1} t_{aa}\nu_a^{N-1}, \quad \text{tr}[\rho I_j \tilde{\sigma}_{I_j}] = \sum_{a,b=0}^{\chi^2-1} \sum_{c,d=0}^{\chi-1} t_{ab}\nu_a^{\ell}\nu_b^{\ell}\Gamma_b^{cd}(\lambda_c\mu_d)^{N-(2\ell+1)}\tilde{\Gamma}_{cd}^a, \qquad (B14) $$

where we have used $|I_j| = 2\ell + 1$, and defined

$$ t_{ab} = \langle r_a^{\rho\sigma} | \tilde{T}_{\rho\sigma} | l_b^{\rho\sigma} \rangle, \qquad (B15) $$

where $\tilde{T}_{\rho\sigma}$ is the transfer matrix $T_{\rho\sigma}$ after the substitution $M^{(j)} \rightarrow \tilde{M}^{(j)}$. In what follows, we assume a choice of $\tilde{M}^{(j)}$ such that $t_{00} \neq 0$. To compute the numerators in Eq. (B13) we first define the tensor

$$ \Theta_{cd}^{ab} = \partial_{\tilde{M}^{(j)}} [\tilde{T}_{\rho\sigma}]_{ab} \cdot \partial_{M^{(j)}} [T_{\rho\sigma}]_{cd}. \qquad (B16) $$
By applying again the spectral decomposition, we obtain

$$
\begin{align}
\bar{C}_{\rho}^{(j)} \cdot C_{\rho}^{(j)} &= \sum_{a,a'=0}^{\chi^2-1} (\nu_a \bar{\nu}_{a'})^{N-1} \Theta_{a'a'}^{aa}, \\
\bar{C}_{\rho}^{(j)} \cdot C_{\rho_{I_j}}^{(j)} &= \sum_{a,a'=0}^{\chi^2-1} \sum_{b,c=0}^{\chi-1} \bar{\nu}_{a'}^{N-1} \Theta_{a'a'}^{aa''} \nu_a^{\ell} \Gamma_a^{bc} (\lambda_b \mu_c)^{N-(2\ell+1)} \bar{\Gamma}_{bc}^{a''} \nu_{a''}^{\ell}, \tag{B17} \\
\bar{C}_{\rho_{I_j}}^{(j)} \cdot C_{\rho_{I_j}}^{(j)} &= \sum_{a,a'=0}^{\chi^2-1} \sum_{b,b'=0}^{\chi-1} (\nu_a \bar{\nu}_{a'})^{\ell} \Gamma_a^{bc} \bar{\Gamma}_{a'}^{b'c'} (\lambda_b \bar{\lambda}_{b'} \mu_c \bar{\mu}_{c'})^{N-(2\ell+1)} \bar{\Gamma}_{bc}^{a''} \bar{\Gamma}_{b'c'}^{a'''} (\nu_{a''} \bar{\nu}_{a'''})^{\ell} \Theta_{a'a'''}^{aa''}.
\end{align}
$$

At this point, in both the overlaps and the scalar products, we separate the term in the summation where all the a indices are 0 from the rest, i.e.

$$
\frac{\bar{C}_{\rho}^{(j)} \cdot C_{\rho}^{(j)}}{\mathrm{tr}[\rho\tilde{\sigma}]} = \frac{\sum_{a,a'} (\nu_a \bar{\nu}_{a'})^{N-1} \Theta_{a'a}^{aa}}{\sum_{a,a'} (\nu_a \bar{\nu}_{a'})^{N-1} t_{aa} \bar{t}_{a'a'}} = \frac{(\nu_0 \bar{\nu}_{0'})^{N-1} \Theta_{00}^{00} + \sum_{a \neq 0 \lor a' \neq 0} (\nu_a \bar{\nu}_{a'})^{N-1} \Theta_{a'a}^{aa}}{(\nu_0 \bar{\nu}_{0})^{N-1} |t_{00}|^2 + \sum_{a \neq 0 \lor a' \neq 0} (\nu_a \bar{\nu}_{a'})^{N-1} t_{aa} \bar{t}_{a'a'}} = \frac{\Theta_{00}^{00}}{|t_{00}|^2} (1 + \varepsilon_1), \quad (\text{B18})
$$

where we have used that overlaps are real, hence tr[ρ̃σ] = ∑_a t_{aa}ν_a^{N-1} = ∑_a t̄_{aa}ν̄_a^{N-1}, and ε₁ is

$$
\epsilon_1 = \frac{\bar{C}_{\rho}^{(j)} \cdot C_{\rho}^{(j)} / (\mathrm{tr}[\rho\tilde{\sigma}])^2}{(\bar{C}_{\rho}^{(j)} \cdot C_{\rho}^{(j)})_0 / (\mathrm{tr}[\rho\tilde{\sigma}]_0)^2} - 1, \quad (\bar{C}_{\rho}^{(j)} \cdot C_{\rho}^{(j)})_0 = (\nu_0 \bar{\nu}_0)^{N-1} \Theta_{00}^{00}, \quad \mathrm{tr}[\rho\tilde{\sigma}]_0 = t_{00} \nu_0^{N-1}. \quad (\mathrm{B}19)
$$

By repeating this step with the two other terms in Eq. (B13), we obtain the same form

$$
\frac{\bar{C}_{\rho}^{(j)} \cdot C_{\rho_{I_j}}^{(j)}}{\mathrm{tr}[\rho\tilde{\sigma}] \mathrm{tr}[\rho_{I_j} \tilde{\sigma}_{I_j}]} = \frac{\Theta_{00}^{00}}{|t_{00}|^2} (1 + \varepsilon_2), \quad \frac{\bar{C}_{\rho_{I_j}}^{(j)} \cdot C_{\rho_{I_j}}^{(j)}}{\mathrm{tr}[\rho_{I_j} \tilde{\sigma}_{I_j}]} = \frac{\Theta_{00}^{00}}{|t_{00}|^2} (1 + \varepsilon_3). \quad (B20)
$$

Under the assumption that the spectrum of $T_{\rho\sigma}$ is gapped, this way of rewriting the terms in Eq. (B13) separates a leading contribution (given by the largest eigenvalues in the spectra of the transfer operators) from a sub-leading contribution. Being the leading contribution constant, and equal among the three terms in Eq. (B13), we obtain

$$
\left\| \frac{C_{\rho}^{(j)}}{\mathrm{tr}[\rho\tilde{\sigma}]} - \frac{C_{\rho_{I_j}}^{(j)}}{\mathrm{tr}[\rho_{I_j}\tilde{\sigma}_{I_j}]} \right\|^2_2 = \frac{\Theta_{00}^{00}}{|t_{00}|^2} (\varepsilon_1 - 2\mathrm{Re}(\varepsilon_2) + \varepsilon_3). \quad (B21)
$$

Note that for $\Theta_{00}^{00} = 0$ the inequality Eq. (A14) is trivial. Therefore, in what follows we assume $\Theta_{00}^{00} \neq 0$.

Upper bounds

We now provide upper bounds for the terms {$\epsilon_m$}. To do so, we first define the ratios

$$
r_1 = \frac{\bar{C}_{\rho}^{(j)} \cdot C_{\rho}^{(j)}}{(\bar{C}_{\rho}^{(j)} \cdot C_{\rho}^{(j)})_0} - 1, \quad r'_1 = \frac{\mathrm{tr}[\rho\tilde{\sigma}]}{(\mathrm{tr}[\rho\tilde{\sigma}])_0} - 1. \qquad (B22)
$$

Similar ratios will appear in the computation of ε₂ and ε₃. From the absolute values of these ratios, it is possible to
derive upper bounds for ε₁ in the following way

$$
|\epsilon_1| = \left| \frac{1 + r_1}{(1 + r'_1)^2} - 1 \right| = \left| \frac{r_1 - r'_1(2 + r'_1)}{(1 + r'_1)^2} \right| \le \frac{|r_1| + |r'_1|(2 + |r'_1|)}{(1 - |r'_1|)^2}. \quad (B23)
$$

Now we make use of the spectral decompositions of the transfer operators to derive upper bounds on the absolute
values of $r_1$ and $r'_1$. From the spectral decomposition in Eq. (B14), and from the definition of the correlation length
$\xi_{\rho\sigma}$ we obtain

$$
|r'_1| \le \sum_{a \ne 0} |t_{aa}/t_{00}| |\nu_a/\nu_0|^{N-1} \le \alpha \chi^2 |\nu_1/\nu_0|^{N-1} = \alpha \chi^2 e^{-(N-1)/\xi_{\rho\sigma}}, \quad (B24)
$$
where $\alpha = \max_{ab} (|t_{ab}/t_{00}|)$. To further simplify the upper bound Eq. (B23), we impose $|r'_1| \le 1/2$, which is satisfied if

$$ \alpha\chi^2 e^{-(N-1)/\xi_{\rho\sigma}} \le 1/2 \rightarrow N \ge \xi_{\rho\sigma}\log(2\alpha\chi^2) + 1. \qquad (B25) $$

Using the spectral decomposition Eq. (B17), analogously to $r'_1$, for $r_1$ we have

$$ |r_1| \le \left| \sum_{a \ne 0 \lor a' \ne 0} \left( \frac{\nu_a \bar{\nu}_{a'}}{\nu_0 \bar{\nu}_0} \right)^{N-1} \frac{\Theta_{a'a'}^{aa}}{\Theta_{00}^{00}} \right| \le 2\beta\chi^2 e^{-(N-1)/\xi_{\rho\sigma}} + \beta\chi^4 e^{-2(N-1)/\xi_{\rho\sigma}} \le c_1\chi^2 e^{-(N-1)/\xi_{\rho\sigma}}, \qquad (B26) $$

where $\beta = \max_{abcd} \left( \frac{\Theta^{ab}}{\Theta_{00}^{cd}} \right) = O(1)$, and $c_1$ is a constant that depends on $\alpha$ and $\beta$. Putting these results together, using Eq. (B23) and $|r'_1| < 1/2$, we finally obtain the upper bound for $|\varepsilon_1|$

$$ |\varepsilon_1| \le 4(|r_1| + 5|r'_1|/2) \le K_1 \chi^2 e^{-(N-1)/\xi_{\rho\sigma}}, \qquad (B27) $$

where $K_1$ is a constant.

The upper bounds for $|\varepsilon_2|$ and $|\varepsilon_3|$ follow from similar steps. In particular, for $|\varepsilon_2|$ we have the upper bound

$$ |\varepsilon_2| \le \left| \frac{1 + r_2}{(1 + r'_1)(1 + r'_2)} - 1 \right| \le 2 \frac{|r_2| + |r'_1| + 3|r'_2|/2}{1 - |r'_2|}, \qquad (B28) $$

obtained in terms of the ratios

$$ r_2 = \frac{\bar{C}_{\rho}^{(j)} \cdot C_{\rho_{I_j}}^{(j)}}{(\bar{C}_{\rho}^{(j)} \cdot C_{\rho_{I_j}}^{(j)})_0} - 1, \quad r'_2 = \frac{\text{tr}[\rho_{I_j} \sigma_{I_j}]}{(\text{tr}[\rho_{I_j} \sigma_{I_j}])_0} - 1. \qquad (B29) $$

First, from Eq. (B14) and the definition above, we derive the upper bound for $|r'_2|$

$$ |r'_2| = \left| \frac{\sum_{a \ne 0 \lor b \ne 0} \sum_{c,d} t_{ab} \nu_a^\ell \nu_b^\ell \Gamma_b^{cd} (\lambda_c \mu_d)^{N-(2\ell+1)} \tilde{\Gamma}_{cd}^a}{\sum_{c,d} t_{00} \nu_0^{2\ell} \Gamma_0^{cd} (\lambda_c \mu_d)^{N-(2\ell+1)} \tilde{\Gamma}_{cd}^0} \right| \le \alpha\gamma \left( 2 \sum_{a \ne 0} |\nu_a/\nu_0|^\ell + \left( \sum_{a \ne 0} |\nu_a/\nu_0|^\ell \right)^2 \right) \le c'_2 \chi^2 e^{-\ell/\xi_{\rho\sigma}}, \qquad (B30) $$

where $\gamma = \max_{abc} \left( \left| \frac{\Gamma_a^{bc} \tilde{\Gamma}_{bc}^a}{\Gamma_0^{bc} \Gamma_{bc}^0} \right| \right) = O(1)$, $c'_2$ is a constant that depends on $\alpha$ and $\gamma$, and we have assumed

$$ \ell \ge \xi_{\rho\sigma}\log(2\chi^2). \qquad (B31) $$

Also in this case we assume $|r'_2| \le 1/2$, which is satisfied when

$$ \ell \ge \xi_{\rho\sigma}\log(2c'_2\chi^2). \qquad (B32) $$

From Eq. (B17), proceeding in a similar way than for the ratios $r_1$ and $r'_2$, we obtain the upper bound for $|r_2|$

$$ |r_2| \le \beta\gamma \left( 3\chi^2 e^{-\ell/\xi_{\rho\sigma}} + 3\chi^4 e^{-2\ell/\xi_{\rho\sigma}} + \chi^6 e^{-3\ell/\xi_{\rho\sigma}} \right) \le c_2\chi^2 e^{-\ell/\xi_{\rho\sigma}}, \qquad (B33) $$

where $c_2$ is a constant. Combining the upper bounds for $|r'_1|$, $|r_2|$ and $|r'_2|$ in Eq. (B28), we obtain

$$ |\varepsilon_2| \le 4(|r_2| + |r'_1| + 3|r'_2|/2) \le K_2 \chi^2 e^{-\ell/\xi_{\rho\sigma}}, \qquad (B34) $$

where $K_2$ is a constant. Without further assumptions, and repeating similar steps, we also obtain

$$ |\varepsilon_3| \le K_3 \chi^2 e^{-\ell/\xi_{\rho\sigma}}, \qquad (B35) $$

where $K_3$ is a constant.
Final result

As a final step, we combine the expression for the norm-2 distance in Eq. (B21) with the upper bounds for the {|ε<sub>m</sub>|} to obtain Eq. (A11). At this point, we can plug in this last equation the upper bounds for {|ε<sub>m</sub>|} Eqs. (B27),(B34),(B35). Under the assumptions that

1. ξ<sub>ρσ</sub> ≥ 0,

2. Γ<sup>bc</sup><sub>0</sub>Γ<sup>0</sup><sub>bc</sub> ≠ 0, t<sub>00</sub> ≠ 0,

3. N ≥ ξ<sub>ρσ</sub>log(2αχ<sup>2</sup>) + 1,

4. ℓ ≥ ξ<sub>ρσ</sub>max (log(2c'<sub>2</sub>χ<sup>2</sup>), log(2χ<sup>2</sup>)) ,

we obtain

$$
\left\| \frac{C_{\rho}^{(j)}}{\text{tr}[\rho \tilde{\sigma}]} - \frac{C_{\rho I_j}^{(j)}}{\text{tr}[\rho_{I_j} \tilde{\sigma}_{I_j}]} \right\|_2^2 = \frac{\Theta_{00}^{00}}{|t_{00}|^2} (\varepsilon_1 - 2\text{Re}(\varepsilon_2) + \varepsilon_3) \leq K \chi^2 e^{-\ell/\xi_{\rho\sigma}}, \quad (B36)
$$

where K = Θ<sup>00</sup><sub>00</sub> / |t<sub>00</sub>|<sup>2</sup> max(K<sub>1</sub>, K<sub>2</sub>, K<sub>3</sub>) is a constant that depends on α, β and γ.

Note that the same derivation can be now repeated by substituting ρ → σ, obtaining

Appendix C: Statistical errors in the learning algorithm

In this section we study analytically the effect of statistical errors in the learning algorithm due to a finite number of measurements N<sub>u</sub><sup>L</sup>N<sub>M</sub>. We consider the framework of classical shadows obtained from local randomized measurements, which we briefly review in the first subsection. Then, we show how to apply it to our learning algorithm, and prove an upper bound on the effect of statistical errors in the update of a single tensor M<sup>(j)</sup>. For a choice of the parameters (ℓ, χ'), we prove that the number of measurement bases needed to obtain the optimal tensor M<sup>(j)</sup> up to an error ε scales like N<sub>u</sub><sup>L</sup> = O(2<sup>2ℓ</sup>χ'<sup>-2</sup>), independently of the total system size N.

1. Local randomized measurements and classical shadows

In this section we review the classical shadow framework for local randomized measurements, and highlight results which we use throughout the manuscript. We refer the reader to Ref. [13] for a more comprehensive review.

Let us assume we have prepared a N-qubit quantum state ρ in an experiment. Local randomized measurement consists in applying a layer of single-qubit random unitaries ⊗<sub>j=1</sub><sup>N</sup> u<sub>j</sub><sup>(r)</sup>, before measuring in the computational basis, where r is the index associated with the bases, and each u<sub>i</sub><sup>(r)</sup> is sampled from the circular unitary ensemble (CUE). We consider measurements in N<sub>u</sub> different random bases, with N<sub>M</sub> shots per basis, for a total of M = N<sub>u</sub>N<sub>M</sub> experimental repetitions. After each experimental shot, we obtain as a result a bitstring s<sup>(r,m)</sup>, with r ∈ [1, N<sub>u</sub>], m ∈ [1, N<sub>M</sub>]. Such dataset can be then post-processed in various ways to extract properties of ρ. Classical shadows represent a powerful method to post-process such data [27]. For an arbitrary subset of qubits I, we can define for the r-th measurement basis a classical shadow

$$
\hat{\rho}_I^{(r)} = \frac{1}{N_M} \sum_{m=1}^{N_M} \bigotimes_{j \in I} \left( 3 \left( u_j^{(r)} \right)^{\dagger} |s_j^{(r,m)}\rangle \langle s_j^{(r,m)}| u_j^{(r)} - \mathbb{I}_j \right), \quad (C1)
$$

where I<sub>j</sub> is the identity on the j-th qubit. Note that classical shadows only depend on the local unitaries u<sub>j</sub><sup>(r)</sup>, which are classically sampled, and on the measurement outcomes s<sub>j</sub><sup>(r,m)</sup>. While for N<sub>M</sub> = 1 classical shadows are product operators, and therefore numerically efficient, the sum over N<sub>M</sub> can be computed efficiently only for modest subsystem sizes |I|.

The main property of classical shadows is that they are unbiased estimators of ρ, i.e.

$$
\rho_I = \mathbb{E}[\hat{\rho}_I^{(r)}], \qquad (C2)
$$
where $\mathbb{E}$ is the quantum mechanical average over all the possible bases and outcomes. Using this property, it is possible to obtain an estimator for any linear function of $\rho$, e.g. for the expectation value of an observable $O$ with support $I$

$$ \hat{\rho} = \frac{1}{N_u} \sum_r \text{tr}[\hat{\rho}_I^{(r)} O] = \text{tr}[\hat{\rho}_I O], \qquad (C3) $$

where we have defined the averaged shadow $\hat{\rho}_I = \frac{1}{N_u} \sum_r \hat{\rho}_I^{(r)}$. The statistical errors of an estimator are related to its variance, which in this case, for $N_M = 1$, admits the upper bound [16, 27]

$$ \text{Var}[\hat{\rho}] \le \frac{2^{|I|}}{M} \|O\|_2^2. \qquad (C4) $$

Analogously, we can write estimators for multi-copy observables

$$ \hat{\rho}^{(n)} = \frac{1}{n!} \binom{N_u}{n}^{-1} \sum_{r_1 \neq r_2 \neq \dots \neq r_n} \text{tr}[\hat{\rho}_I^{(r_1)} \hat{\rho}_I^{(r_2)} \dots \hat{\rho}_I^{(r_n)} O^{(n)}] \qquad (C5) $$

and compute associated bounds on their variance [18], which generally scale exponentially with the support of $O^{(n)}$.

Relevant examples of such multi-copy observables include the moments of the reduced density matrix $\mathcal{P}_n(\rho_I) = \text{tr}[\rho_I^n]$. However, it is also possible to estimate these moments in a different way, that does not rely on the knowledge of the unitaries $\{u_j^{(r)}\}$ [12, 15]. For example, the purity $\mathcal{P}_2(\rho_I)$ can be estimated as

$$ \hat{\mathcal{P}}_2(\rho_I) = \frac{2^{|I|}}{N_u N_M (N_M - 1)} \sum_{r=1}^{N_u} \sum_{m \neq m'=1}^{N_M} 2^{-D[s_I^{(r,m)}, s_I^{(r',m')}]}, \qquad (C6) $$

where $D[s, s'] = \sum_j |s_j - s'_j|$ is the Hamming distance between two bitstrings, and $s_I^{(r,m)}$ is the bitstring obtained from the $m$-th shot in the $r$-th basis restricted to the subsystem $I$. This estimator also requires $M = O(2^{|I|})$ measurements with $a = O(1)$ for finite accuracy [12], but is computationally less expensive than using classical shadows. Throughout this work, we use Eq. (C6) to estimate purities, and classical shadows otherwise.

## 2. Statistical errors in the optimization of a single tensor

In this section we study the statistical errors in the update of a local tensor using randomized measurement data. The update rule Eq. (A11) requires experimental data to estimate $C_{\rho_{Ij}}^{(j)} = \text{tr}[\rho_{Ij} \partial_{M^{(j)}} \sigma_{Ij}]$. Being $C_{\rho_{Ij}}^{(j)}$ linear in $\rho$, we can use the estimator Eq. (C3) for each one of its terms, obtaining

$$ \hat{C}_{\rho_{Ij}}^{(j)} = \frac{1}{N_u^L} \sum_{r=1}^{N_u^L} \text{tr}[\hat{\rho}_{Ij}^{(r)} \partial_{M^{(j)}} \sigma_{Ij}], \qquad (C7) $$

where $\hat{\rho}_{Ij}^{(r)}$ are classical shadows on the subsystem $I_j = [j - \ell, j + \ell]$, and $N_u^L$ is the number of random measurement bases collected in the learning set. The statistical errors in this estimator propagate to the tensor $M^{(j)}$, which will be then updated $\hat{M}^{(j)} \propto C_{\sigma_{Ij}}^{(j)-1} \hat{C}_{\rho_{Ij}}^{(j)}$. Let us define $M^{(j)} = C_{\sigma_{Ij}}^{(j)-1} C_{\rho_{Ij}}^{(j)} \propto M^{(j)}$. Using the Cauchy-Schwarz inequality, we obtain the propagation of errors

$$ \|\hat{M}^{(j)} - M^{(j)}\|_2 \le \|C_{\sigma_{Ij}}^{(j)-1}\|_{\infty} \|\hat{C}_{\rho_{Ij}}^{(j)} - C_{\rho_{Ij}}^{(j)}\|_2, \qquad (C8) $$

where $\|\cdot\|_2$ is the Hilbert-Schmidt norm (see Eq. (B11)), and $\|O\|_{\infty} = \max_v \frac{\|O_v\|_2}{\|v\|_2}$ is the spectral norm. We expect that the errors in the normalization of $M^{(j)}$ will amount to a constant factor.

We now study the statistical errors in norm-2 $\|\hat{C}_{\rho_{Ij}}^{(j)} - C_{\rho_{Ij}}^{(j)}\|_2$. The statistical errors of the estimator Eq. (C7) admit an upper bound via the multidimensional Chebyshev inequality [69]

$$ \text{Pr} \left( \|\hat{V} - V\|_2 \ge \varepsilon \right) \le \frac{1}{\varepsilon^2} \mathbb{E} \left[ \|\hat{V} - V\|_2^2 \right] = \frac{1}{\varepsilon^2} \sum_j \text{Var}(\hat{V}_j), \qquad (C9) $$
with $\mathbb{E}[\hat{V}] = V$, from which we obtain

$$ \mathrm{Pr} (||\hat{C}_{\rho_{I_j}}^{(j)} - C_{\rho_{I_j}}^{(j)}||_2 \ge \varepsilon) \le \frac{1}{\varepsilon^2} \sum_{s,s'} \sum_{l,r=1}^{\chi'} \mathrm{Var} ([\hat{C}_{\rho_{I_j}}^{(j)}]_{lrss'}) , \quad (C10) $$

where we have assumed that the bond dimension of $\sigma$ is $\chi'$ for all links. Assuming that there is $N_M = 1$ shot per random basis, we can apply the upper bound on the variance for estimators of expectation values Eq. (C4), obtaining

$$ \mathrm{Var} ([\hat{C}_{\rho_{I_j}}^{(j)}]_{lrss'}) \le \frac{2^{2\ell+1}}{N_u^L} \langle rr | C_{\sigma_{I_j}}^{(j)} | ll \rangle , \quad (C11) $$

where $|ll\rangle = |l\rangle \otimes |l\rangle$, $|rr\rangle = |r\rangle \otimes |r\rangle$. Summing over all the indices results in the upper bound

$$ \sum_{l,r,s,s'} \mathrm{Var} ([\hat{C}_{\rho_{I_j}}^{(j)}]_{lrss'}) \le \frac{2^{2\ell+3}\chi'}{N_u^L} \langle \delta | C_{\sigma_{I_j}}^{(j)} | \delta \rangle , \quad (C12) $$

where $|\delta\rangle = \frac{1}{\sqrt{\chi'}} \sum_l |ll\rangle$. We expect that the term $\langle \delta | C_{\sigma_{I_j}}^{(j)} | \delta \rangle$ contributes with a factor proportional to the purity $\mathrm{tr}[\sigma_{I_j}^2]$. This can be seen for translation-invariant MPOs using the transfer operator formalism introduced in App. B. At the leading order in the spectrum of $T_{\sigma}^{(1)}$ and $T_{\sigma}^{(2)}$, we obtain

$$ \langle \delta | C_{\sigma_{I_j}}^{(j)} | \delta \rangle \approx \frac{\langle r_0^{\sigma^2} | \delta \rangle \langle \delta | l_0^{\sigma^2} \rangle}{\nu_0} \mu_0^{N-2\ell-1} \nu_0^{2\ell+1} \langle r_0^{\sigma^2} | l_0^{\sigma} l_0^{\sigma} \rangle \langle r_0^{\sigma} r_0^{\sigma} | l_0^{\sigma^2} \rangle, \quad (C13) $$

while for the purity

$$ \mathrm{tr}[\sigma_{I_j}^2] \approx \mu_0^{N-2\ell-1} \nu_0^{2\ell+1} \langle r_0^{\sigma^2} | l_0^{\sigma} l_0^{\sigma} \rangle \langle r_0^{\sigma} r_0^{\sigma} | l_0^{\sigma^2} \rangle . \quad (C14) $$

This implies that, up to corrections of order $O(\chi'^2 e^{-\ell/\xi_{\sigma}^{(2)}})$, we have

$$ \langle \delta | C_{\sigma_{I_j}}^{(j)} | \delta \rangle \approx c \mathrm{tr}[\sigma_{I_j}^2] , \quad (C15) $$

where $c = \frac{\langle r_0^{\sigma^2} | \delta \rangle \langle \delta | l_0^{\sigma^2} \rangle}{\nu_0} = O(1)$ is a constant factor. Under the assumption that $\ell > O(\xi_{\sigma}^{(2)} \log(\chi'^2))$ we can absorb the sub-leading terms in a constant $c'$, obtaining the final result

$$ \mathrm{Pr} (||\hat{C}_{\rho_{I_j}}^{(j)} - C_{\rho_{I_j}}^{(j)}||_2 \ge \varepsilon) \le c' \mathrm{tr}[\sigma_{I_j}^2] \frac{2^{2\ell+3}\chi'}{\varepsilon^2 N_u^L} . \quad (C16) $$

This implies that the number of measurements needed to achieve a controlled statistical error in norm $-2 ||\hat{M}^{(j)} - M^{(j)}||_2$ scales like $M^L = N_u^L N_M = O(2^{2\ell} \chi' \varepsilon^{-2})$, as reported in the main text, and is independent of the total system size $N$.

## Appendix D: Scalable estimation of mixed-state fidelities

In this section we discuss in more detail the AFC estimation of mixed-state fidelities introduced in the main text. First, leveraging the formalism of transfer operators, we prove approximate factorization conditions (AFC) for the overlap of two TIMPDOs with finite correlation lengths. Then, we argue that combining them with AFC for purities [21] allows us to define the AFC for the fidelities defined in the main text. This implies that a measurement budget of size polynomial in $N$ is sufficient to estimate fidelities in these quantum states. Lastly, we numerically demonstrate that a modest measurement budget is sufficient to estimate AFC fidelities in thermal states of up to $N = 128$ qubits.
# 1. Mixed-state fidelities from randomized measurements

First, we comment on the relation between the two mixed-state fidelities introduced in the main text: the GM fidelity (Eq. (A2)) and the max fidelity, defined as [24, 60]

$$ \mathcal{F}_{\max} = \frac{\mathrm{tr}[\rho\sigma]}{\max(\mathrm{tr}[\rho^2], \mathrm{tr}[\sigma^2])} \qquad (\mathrm{D}1) $$

We can see that $\mathcal{F}_{\mathrm{GM}}(\rho, \sigma) \ge \mathcal{F}_{\max}(\rho, \sigma)$, and $\mathcal{F}_{\max}(\rho, |\psi\rangle\langle\psi|) = \langle\psi|\rho|\psi\rangle$, while the same does not hold for $\mathcal{F}_{\mathrm{GM}}$. Additionally, the GM fidelity is insensitive to certain noise models such as global depolarization, e.g. $\mathcal{F}_{\mathrm{GM}}(\rho, (1-p)\rho+p\mathbb{I}/2^N) \sim 1 - O(\mathrm{tr}[\rho^2]/2^N)$. This makes the max fidelity more suitable to physically distinguish quantum states, while the GM fidelity, being differentiable, is a better-suited cost function for optimization algorithms.

It is possible to evaluate both fidelities from randomized measurements by estimating sequentially the purities $\mathrm{tr}[\rho^2]$ and $\mathrm{tr}[\sigma^2]$ and their overlap $\mathrm{tr}[\rho\sigma]$ (see Sec. C1). In the context of our protocol, we are especially interested in estimating the fidelity between an experimentally prepared state $\rho$ and a known state $\sigma$; applying Eq. (C4), we obtain

$$ \mathrm{Var}[\mathrm{tr}[\hat{\rho}\sigma]] \le \frac{2^N}{M} \mathrm{tr}[\sigma^2]. \qquad (\mathrm{D}2) $$

This implies that the overlap can be estimated up to controlled statistical errors with $M = O(2^N)$ randomized measurements. Given that estimating purities requires instead $M = O(4^N)$ [18], we generally need at least $M = O(4^N)$ measurements to estimate the previously-defined fidelities.

# 2. Approximate factorization conditions for fidelities

We now introduce approximate factorization conditions (AFC) for the above-defined fidelities, and show how they lift the requirement of exponentially-many randomized measurements for fidelity estimation. First, let us recall AFC for purities, defined in Ref. [21]. A quantum state $\rho$ satisfies AFC for purities if it is possible to approximate its global purity by

$$ \mathrm{tr}[\rho^2] \approx \mathcal{P}_{\mathrm{AFC}}^{(k)}(\rho) = \frac{\mathrm{tr}[\rho_{A_1 A_2}^2]\mathrm{tr}[\rho_{A_2 A_3}^2]\dots\mathrm{tr}[\rho_{A_{R-1} A_R}^2]}{\mathrm{tr}[\rho_{A_2}^2]\mathrm{tr}[\rho_{A_3}^2]\dots\mathrm{tr}[\rho_{A_{R-1}}^2]}, \qquad (\mathrm{D}3) $$

where $A_j$ are connected intervals of $|A_j| = k$ qubits, with relative error

$$ \left| \frac{\mathcal{P}_{2, \mathrm{AFC}}^{(k)}(\rho)}{\mathcal{P}_2(\rho)} - 1 \right| \le c_1 \frac{N}{k} e^{-k/\xi}. \qquad (\mathrm{D}4) $$

For instance, a TIMPDO $\rho$ with bond dimension $\chi$ and correlation lengths $\xi_{\rho}^{(1)}, \xi_{\rho}^{(2)} = O(1)$ satisfies AFC for purities, with $\xi = \max(\xi_{\rho}^{(1)}, \xi_{\rho}^{(2)})$, and $c_1 = O(\chi^2)$. From this approximation, it has been proven that for $k > \xi$ the purity estimation is guaranteed to have a controlled error with $M = O(2^{4k}N^3)$ measurements [21].

In a similar fashion, we now define AFC for the overlap between two quantum states. We say that two quantum states $\rho$ and $\sigma$ satisfy AFC for the overlap if the following approximation holds

$$ \mathrm{tr}[\rho\sigma] \approx \mathcal{O}_{\mathrm{AFC}}^{(k)}(\rho, \sigma) = \frac{\mathrm{tr}[\rho_{A_1 A_2}\sigma_{A_1 A_2}]\mathrm{tr}[\rho_{A_2 A_3}\sigma_{A_2 A_3}]\dots\mathrm{tr}[\rho_{A_{R-1} A_R}\sigma_{A_{R-1} A_R}]}{\mathrm{tr}[\rho_{A_2}\sigma_{A_2}]\mathrm{tr}[\rho_{A_3}\sigma_{A_3}]\dots\mathrm{tr}[\rho_{A_{R-1}}\sigma_{A_{R-1}}]}, \qquad (\mathrm{D}5) $$

with $A_j$ connected intervals of $|A_j| = k$ qubits, where the relative error is bounded by

$$ \left| \frac{\mathcal{O}_{\mathrm{AFC}}^{(k)}(\rho, \sigma)}{\mathrm{tr}[\rho\sigma]} - 1 \right| \le c_2 \frac{N}{k} e^{-k/\xi} \qquad (\mathrm{D}6) $$

for some constants $c_2, \xi$. When two states $\rho$ and $\sigma$ satisfy both the AFC for the overlap and, respectively, for purities, we say that they satisfy AFC for the fidelity. To justify this definition, consider a choice of $k$ for a pair of states $\rho$ and $\sigma$ satisfying AFC such that

$$ \left| \frac{\mathcal{O}_{\mathrm{AFC}}^{(k)}(\rho, \sigma)}{\mathcal{O}(\rho, \sigma)} - 1 \right| \le \varepsilon_{\mathcal{O}}, \quad \left| \frac{\mathcal{P}_{2, \mathrm{AFC}}^{(k)}(\rho)}{\mathcal{P}_2(\rho)} - 1 \right| \le \varepsilon'_{\mathcal{P}}, \quad \left| \frac{\mathcal{P}_{2, \mathrm{AFC}}^{(k)}(\sigma)}{\mathcal{P}_2(\sigma)} - 1 \right| \le \varepsilon''_{\mathcal{P}}, \qquad (\mathrm{D}7) $$
where $\epsilon_P = (\epsilon'_P + \epsilon''_P + \epsilon'_P\epsilon''_P) \le 1/2$. Then, it is possible to see that

$$ \left| \frac{\mathcal{F}_{\text{AFC}}^{(k)}(\rho, \sigma)}{\mathcal{F}(\rho, \sigma)} - 1 \right| \le 2(\epsilon_O + \epsilon_P) = O\left(\text{poly}\left(\frac{N}{k}e^{-k/\xi}\right)\right) \quad (\text{D8}) $$

for both $\mathcal{F} = \mathcal{F}_{\text{max}}$ and $\mathcal{F} = \mathcal{F}_{\text{GM}}$, where

$$ \mathcal{F}_{\text{max,AFC}}^{(k)}(\rho, \sigma) = \frac{\mathcal{O}_{\text{AFC}}^{(k)}(\rho, \sigma)}{\max(\mathcal{P}_{2,\text{AFC}}^{(k)}(\rho), \mathcal{P}_{2,\text{AFC}}^{(k)}(\sigma))}, \quad \mathcal{F}_{\text{GM,AFC}}^{(k)}(\rho, \sigma) = \frac{\mathcal{O}_{\text{AFC}}^{(k)}(\rho, \sigma)}{\sqrt{\mathcal{P}_{2,\text{AFC}}^{(k)}(\rho)\mathcal{P}_{2,\text{AFC}}^{(k)}(\sigma)}}. \quad (\text{D9}) $$

This means that, for a pair of states $\rho$ and $\sigma$ satisfying AFC for the fidelity, there exists a value of $k$ such that $\mathcal{F}_{\text{AFC}}^{(k)}(\rho, \sigma) \simeq \mathcal{F}(\rho, \sigma)$. Note that this is trivially true for $\rho = \sigma$, as $\mathcal{F}_{\text{AFC}}^{(k)}(\rho, \rho) = 1$ for any $k$. In App. D3 we prove that two translation-invariant MPDOs $\rho$ and $\sigma$ with finite correlation lengths satisfy the AFC for the overlap (hence for the fidelity), with $c_2 = O(\chi^2)$, $\chi$ being the largest bond dimension between $\rho$ and $\sigma$, and $\xi = \max(\xi_{\rho}^{(1)}, \xi_{\sigma}^{(1)}, \xi_{\rho\sigma})$ (defined in App. B).

Finally, let us comment on the statistical errors in fidelity estimation using AFC. In order to estimate the AFC purity of $N$ qubits with controlled errors $M = O(4^{2k}N^3)$ randomized measurements are sufficient [21]. Repeating the same steps of the proof, it is possible to show that $M = O(2^{2k}N^3)$ randomized measurements are sufficient to estimate $\mathcal{O}_{\text{AFC}}^{(k)}(\rho, \sigma)$ with controlled errors. This implies that, asymptotically, $O(4^{2k}N^3)$ measurements are sufficient to estimate AFC fidelities accurately, and we can in principle use this strategy to benchmark the results of the learning algorithm. In App. D4 we demonstrate numerically that the number of measurements needed is in fact within the reach of current experiments.

### 3. Approximate factorization conditions for the overlap: proof for TIMPDOs

In this section we prove that, under technical assumptions similar to the ones defined in Sec. B, two TIMPDOs satisfy AFC for the overlap. In the language of Sec. B, we take two TIMPDOs $\rho$ and $\sigma$ such that

$$ \xi = \max(\xi_{\rho}^{(1)}, \xi_{\sigma}^{(1)}, \xi_{\rho\sigma}), \quad \Gamma_{0}^{bc}\tilde{\Gamma}_{bc}^{0} \neq 0. \quad (\text{D10}) $$

From the spectral decompositions of the transfer operators, we obtain the following expression for the overlap

$$ \mathcal{O}(\rho, \sigma) = \frac{\text{tr}[\rho\sigma]}{\text{tr}[\rho]\text{tr}[\sigma]} = \frac{\sum_{a=0}^{\chi^2-1} (\nu_a)^N}{\sum_{b,c=0}^{\chi-1} (\lambda_b \mu_c)^N} = \left(\frac{\nu_0}{\lambda_0 \mu_0}\right)^N (1+\delta), \quad (\text{D11}) $$

where, from the definition

$$ \delta = \frac{(1 + \sum_{a=1}^{\chi^2-1} (\nu_a / \nu_0)^N)}{(1 + \sum_{b=1}^{\chi-1} (\lambda_a / \lambda_0)^N)(1 + \sum_{c=1}^{\chi-1} (\mu_a / \mu_0)^N)} - 1 \quad (\text{D12}) $$

and the ones of the correlation lengths, we obtain the upper bound

$$ |\delta| \le 4(\chi^2 e^{-N/\xi_{\rho\sigma}} + \chi e^{-N/\xi_{\rho}^{(1)}} + \chi e^{-N/\xi_{\sigma}^{(1)}} + \chi^2 e^{-2N/\max(\xi_{\rho}^{(1)}, \xi_{\sigma}^{(1)})}) \le 14\chi^2 e^{-N/\xi}, \quad (\text{D13}) $$

with $\xi = \max(\xi_{\rho}^{(1)}, \xi_{\sigma}^{(1)}, \xi_{\rho\sigma})$, and under the further assumptions that

$$ \left| \sum_{b=1}^{\chi-1} (\lambda_a / \lambda_0)^N \right| \le 1/2, \quad \left| \sum_{b=1}^{\chi-1} (\mu_a / \mu_0)^N \right| \le 1/2, \quad (\text{D14}) $$

which are satisfied for

$$ N \ge \xi \log(2\chi). \quad (\text{D15}) $$

Following similar steps, we can also express the overlap on an interval $A_l$, with $|A_l| = k$, as

$$ \mathcal{O}(\rho_{A_l}, \sigma_{A_l}) = \text{tr}[\rho_{A_l} \sigma_{A_l}] = \sum_{abc} \nu_a^k \Gamma_a^{bc} (\lambda_b \mu_c)^{N-k} \tilde{\Gamma}_{bc}^a = \nu_0^k (\lambda_0 \mu_0)^{N-k} \Gamma_0^{00} \tilde{\Gamma}_{00}^0 (1 + \delta_k). \quad (\text{D16}) $$
In this case, $\delta_k$ admits the upper bound

$$|\delta_k| \le \beta(\chi^2 e^{-k/\xi} + 2\chi e^{-(N-k)/\xi} + 2\chi^3 e^{-N/\xi} + \chi^2 e^{-2(N-k)/\xi} + \chi^4 e^{-(2N-k)/\xi}) \le \frac{7}{2}\beta\chi^2 e^{-k/\xi}, \quad (D17)$$

where $\beta = \max_{abc} \left( \left| \frac{\Gamma_a^{bc} \tilde{\Gamma}_{bc}^a}{\Gamma_0^{bc} \Gamma_{bc}^0} \right| \right)$ is a constant, and we have assumed that

$$\chi e^{-(N-2k)/\xi} \le 1/2 \to N \ge 2k + \xi \log(2\chi). \quad (D18)$$

Note that we can obtain a similar form for the overlap $\mathcal{O}(\rho_{A_l A_{l+1}}, \sigma_{A_l A_{l+1}})$ by simply performing a substitution $k \to 2k$. From these two results, we can write the AFC overlap in periodic boundary conditions as

$$\mathcal{O}_{\text{AFC}}^{(k)}(\rho, \sigma) = \frac{\prod_{l=1}^{R} \mathcal{O}(\rho_{A_l A_{l+1}}, \sigma_{A_l A_{l+1}})}{\prod_{l=1}^{R} \mathcal{O}(\rho_{A_l}, \sigma_{A_l})} = \left(\frac{\nu_0}{\lambda_0 \mu_0}\right)^N \left(\frac{1 + \delta_{2k}}{1 + \delta_k}\right)^R = \left(\frac{\nu_0}{\lambda_0 \mu_0}\right)^N (1 + \delta'), \quad (D19)$$

where $R = N/k$, assuming $N$ is a multiple of $k$, and

$$\delta' = \left( \frac{1 + \delta_{2k}}{1 + \delta_k} \right)^R - 1. \quad (D20)$$

At this point, the errors induced by the AFC form can be expressed as

$$\left| \frac{\mathcal{O}_{\text{AFC}}^{(k)}(\rho, \sigma)}{\mathcal{O}(\rho, \sigma)} - 1 \right| = \left| \frac{1 + \delta'}{1 + \delta} - 1 \right|. \quad (D21)$$

Now, let us assume that $|\delta| \le 1/2$, which is satisfied when

$$N \ge \xi \log(28\chi^2). \quad (D22)$$

Under this assumption, we obtain an upper bound on the error induced by the AFC form of the overlap as a function of $|\delta|$ and $|\delta'|$

$$\left| \frac{\mathcal{O}_{\text{AFC}}^{(k)}(\rho, \sigma)}{\mathcal{O}(\rho, \sigma)} - 1 \right| \le 2|\delta' - \delta| \le 2(|\delta| + |\delta'|). \quad (D23)$$

Since we already have an upper bound on $|\delta|$, the last step is to obtain an upper bound on $|\delta'|$. First of all, let us define $\delta'' \ge \max(|\delta_k|, |\delta_{2k}|)$, such that we can write

$$|\delta'| = \left| \left( \frac{1 + \delta_{2k}}{1 + \delta_k} \right)^R - 1 \right| \le \left| \left( \frac{1 + \delta''}{1 - \delta''} \right)^R - 1 \right| \le |e^{4R\delta''} - 1|. \quad (D24)$$

As a last assumption, we take

$$4R\delta'' \le 1, \quad (D25)$$

which is satisfied for

$$k \ge \xi \log(14\beta N \chi^2). \quad (D26)$$

Under this condition, we can use that $e^z - 1 \le 2z$ for $0 \le z \le 1$ and Eq. (D17), obtaining

$$|\delta'| \le |e^{4R\delta''} - 1| \le 8R\delta'' \le 28\beta\chi^2 \frac{N}{k} e^{-k/\xi}. \quad (D27)$$

At this point we can combine this bound with the one on $|\delta|$ (Eq. (D13)) to obtain a bound on the AFC error for the overlap Eq. (D23). To summarize, assuming

1. Eq. (D10)
<figure><img src="image_8.png" /><figcaption>FIG. S4. Numerical results for AFC fidelity estimation between the Ising-Gibbs states $\rho_0$ and $\rho_1$, defined in Eq. (D29). (a) Systematic relative error in the max fidelity $F_{\max}(\rho_0, \rho_1)$ due to the AFC assumption at finite $k$ for various $N$. (b) Fidelity estimation from randomized measurements, assuming AFC with $k=1,3$ of $F_{\max}(\rho_0, \rho_1)$ for various $N$ using $N_u = 1000$ random bases and $N_M = 512$. It is possible to see that the AFC estimation with a finite number of measurements (markers) faithfully captures the exact behavior of the fidelity (dotted line). (c,d) Relative errors between the AFC estimate and the exact value of fidelities for $k=1,3, N_M = 512$, and various $N_u$ and $N$. In the regime considered, we observe a power-law decay of these errors as a function of $N_u N_M / N$.</figcaption></figure>

$$2. N \geq \max(2k + \xi \log(2\chi), \xi \log(28\chi^2))$$

$$3. k \geq \xi \log(14\beta N\chi^2),$$

we finally obtain the upper bound on AFC errors in the overlap

$$\left| \frac{O_{\text{AFC}}^{(k)}(\rho, \sigma)}{O(\rho, \sigma)} - 1 \right| \leq 28\chi^2 e^{-k/\xi} \left( e^{-(N-k)\xi} + 2\beta N/k \right) \leq 28(2\beta + 1)\chi^2 \frac{N}{k} e^{-k/\xi}, \quad (\text{D28})$$

which corresponds to the approximate factorization conditions for the overlap Eq. (D6), with $c_2 = 28(2\beta + 1)\chi^2$ and $\xi = \max(\xi_{\rho}^{(1)}, \xi_{\sigma}^{(1)}, \xi_{\rho\sigma})$, with $\beta = \max_{abc} \left( \left| \frac{\Gamma_{a}^{bc}\Gamma_{bc}^{a}}{\Gamma_{0}^{bc}\Gamma_{bc}^{0}} \right| \right)$.

### 4. Statistical errors in AFC fidelity estimation: numerical results

In this section, we numerically investigate the experimental feasibility of AFC fidelity estimation. While we have previously argued that $O(4^{2k}\text{poly}(N))$ measurements are sufficient for an accurate estimation, this does not give a practical indication of the number of measurements typically needed. Let us consider Ising-Gibbs states, defined as

$$\rho(\beta, g, h) = \frac{1}{Z} e^{-\beta H_{\text{Ising}}(g,h)}, \quad H_{\text{Ising}}(g,h) = \frac{1}{4} \left( \sum_{j=1}^{N-1} Z_j Z_{j+1} + \sum_{j=1}^{N} (g X_j + h Z_j) \right), \quad (\text{D29})$$

where $\beta$ is the inverse temperature, and the factor $Z = \text{tr}[e^{-\beta H_{\text{Ising}}(g,h)}]$ ensures the normalization of the state. These states are known to satisfy AFC for the purity at finite temperature [21, 54]. We consider the two Ising-Gibbs states $\rho_0 = \rho(\beta = 2, g = 1.01, h = 0.04)$ and $\rho_1 = (1, 1.5, 0)$ at sizes up to $N = 128$. We numerically simulate these states as MPDOs using the ITensor library [96] by evolving the infinite temperature state in imaginary time using time-evolving block decimation (TEBD). We fix their bond dimension by fixing the truncation error to $10^{-20}$. This choice results in bond dimension $\chi = (29, 20)$ for the states $(\rho_0, \rho_1)$ respectively.

As a first step, we study the relative error induced by the AFC approximation of the max fidelity $F_{\max}(\rho_0, \rho_1)$ between these two states at finite $k$, defined as

$$\varepsilon_{\text{AFC}} = \left| \frac{\mathcal{F}_{\text{max},\text{AFC}}^{(k)}(\rho_0, \rho_1)}{\mathcal{F}_{\text{max}}(\rho_0, \rho_1)} - 1 \right|, \quad (\text{D30})$$

and observe that it decays exponentially with $k$ (see Fig. S4). This is strong numerical evidence that $\rho_0$ and $\rho_1$ satisfy the AFC for the fidelity. Then, using the RandomMeas library [95], we simulate local randomized measurements on
<figure><img src="image_9.png" /><figcaption>FIG. S5. Numerical performances of the learning algorithm as a function of (ℓ, χ') for $N_u^L \to \infty$. We consider Ising-Gibbs states as defined in Eq. (D29) for $(g, h) = (1.01, 0.04)$ for various inverse temperatures $\beta$. We study both $\mathcal{F}_{\text{GM}}$ and $\mathcal{F}_{\text{max}}$, plotting both the global fidelities (dashed lines) and the AFC fidelities (markers) with $k = \ell + 1$. (a,b) Results for $\beta = 2$ as a function of $(\ell, \chi')$. We observe signatures of an exponential decay of the fidelity errors $1 - \mathcal{F}$ with $\ell$, while the bond dimension $\chi'$ induces a lower bound on the fidelity error. The AFC fidelities generally match the global fidelities. (c, d) Results for $\chi' = 16$ as a function of $\ell$, for various inverse temperatures $\beta$.</figcaption></figure>

$\rho_0$, and estimate its fidelity with both $\rho_1$ and itself as known target states. We consider $N_M = 512$ shots per basis, and up to $N_u = 9000$ local random bases, and report the results in Fig. S4. First of all, it is possible to see that $N_u = 1000$ is already sufficient to accurately estimate these fidelities up to $N = 128$. To perform a full scaling analysis we plot the relative error between the AFC estimation at finite $k = 1, 3$ from $M = N_u N_M$ measurements and the exact value of the fidelity, defined as

$$ \varepsilon_{s,\text{AFC}} = \left| \frac{\hat{\mathcal{F}}_{\text{max},\text{AFC}}^{(k)}(\rho_0, \rho_1)}{\mathcal{F}_{\text{max}}(\rho_0, \rho_1)} - 1 \right| \quad (\text{D31}) $$

First of all we observe that, in this regime, the statistical errors are orders of magnitude larger than the AFC errors. This makes the AFC shadow estimator a good estimator of the fidelity, since its bias is negligible compared to its variance for an experimentally feasible number of measurements. By considering various choices of $N$ and $N_u$, we observe a power-law decay of $\varepsilon_{s,\text{AFC}}$ with $N_u N_M / N$. We note that the power-law decaying behavior of errors with the number of measurements is typical of randomized measurements, see for instance Ref. [24]. Furthermore, the data collapse hints at a linear scaling of the number of measurements required for an accurate estimation with respect to the total number of qubits $N$.

## Appendix E: Numerical results for the learning algorithm

In this section we report on the performances of the learning algorithm by performing numerical simulations on Ising-Gibbs states, defined in Eq. (D29), with $(g, h) = (1.01, 0.04)$. First, we study the performance of the learning algorithm as a function of $(\ell, \chi')$, with $N_u N_M \to \infty$. In this case, the inputs of the learning algorithm are all the reduced density matrices on connected subsets of size $2(\ell + 1)$ (we consider 2-site updates). This allows us to isolate the effect of the approximate local updates introduced in App. A, which we show to be negligible at finite $\ell$ at various temperatures. Finally, we perform a numerical experiment with finite $M^L = N_u^L N_M$ measurements at $\beta = 2$ to study the typical number of measurements required to learn such states. This also serves to further validate the experimental results presented in the main text.

In both cases, we study the behavior of both the geometric-mean fidelity $\mathcal{F}_{\text{GM}}$ and the max fidelity $\mathcal{F}_{\text{max}}$. While the former has a stronger dependence on the parameters of the learning algorithm, they generally display a similar behavior.

### 1. Approximate local updates

Here we numerically study the performance of the learning algorithm as a function of $(\ell, \chi')$ for $N_u N_M \to \infty$ on Ising-Gibbs states (Eq. (D29)). We generate MPOs approximating such states using the TEBD algorithm in imaginary time, implemented with the ITensor library [96]. By imposing a cut-off of $10^{-20}$ on singular values, we obtain MPOs
<figure><img src="image_10.png" /></figure>

FIG. S6. Numerical performances of the learning algorithm with a finite randomized measurement dataset. We consider Ising-Gibbs states (Eq. (D29)) with $(\beta, g, h) = (2, 1.01, 0.04)$, and consider a dataset with $N_u^L + N_u^T = 10000$ bases, with $N_M = 1024$ shots per basis. (a,b) Fidelities as a function of the total number of qubits $N$ for $(\ell, \chi') = (2, 4)$, for which we know the learning algorithm outputs states with fidelities $F \sim 1 - 10^{-4}$. We plot both the exact fidelities (dotted lines) and the estimated AFC fidelity with the testing set for various $k$ (markers) for $N_u^L = 1000$, where only the latter are experimentally accessible. (c) Global purities of the target state $\rho$ and the output $\sigma$ for various $N$. (d,e) Fidelity errors as a function of $N_u^L$ for $N = 16, 128$, $\chi' = 4$, and $\ell = 1, 2$, both numerically exact (dotted lines) and from AFC estimation with $k = \ell + 1$ (markers). (f) Decay of two-point correlations as a function of the distance, compared between the target state (dashed line) and the learned MPO (straight line), averaged over $j \in [10, 113]$.

with bond dimensions $\chi = (29, 45, 57)$ for $\beta = (2, 4, 8)$, and $N = 128$. To mimic the experimental scenario where there is no notion of bond dimension, we learn these states up to a strictly smaller bond dimension $\chi' < \chi$, where $\chi'$ is the largest bond dimension allowed by the two-site updates.

We report the results in Fig. S5. In every plot, we show fidelity errors both computed exactly (dashed lines) and using AFC with $k = \ell + 1$ (markers). The first one allows us to understand the behavior of the learning algorithm, while the second one is what we can reasonably access in an experimental scenario. We generally observe a very good agreement between the two.

For $\beta = 2$ we observe an exponential decay of the fidelity errors $1 - F_{\text{GM}}$ and $1 - F_{\text{max}}$ as a function of $\ell$. Such decay saturates to a value dependent only on $\chi'$, implying that only truncation errors are present at that point. For $\ell = 1$ and $\chi' = 4$ we already observe $1 - F \sim 10^{-4}$ for both fidelities, meaning that for these parameters the state $\rho$ can be learned to a very high accuracy. We observe a similar behavior for other values of $\beta$. By fixing the bond dimension $\chi' = 16$, we still observe an exponential decay of the fidelity errors with $\ell$, with smaller rates for higher $\beta$.

## 2. Statistical errors

In this section we study the effect of statistical errors on the learning algorithm, induced by a finite number of measurements $N_u N_M$. To do so, we numerically sample randomized measurement outcomes from the MPO $\rho$, using the library RandomMeas [95]. Here, we consider $\rho$ to be the Ising-Gibbs state with $(\beta, g, h) = (2, 1.01, 0.04)$, which has bond dimension $\chi = 29$ for a cut-off $10^{-20}$ on singular values. We collect data in a total of $N_u^L + N_u^T = 10^4$ bases, and consider $N_M = 1024$ shots per basis. For each interval $I_j$ we construct the average shadow $\hat{\rho}_{I_j} = \frac{1}{N_u^L N_M} \sum_r \hat{\rho}_{I_j}^{(r)}$ from the randomized measurement data in the learning set, and use them as estimators of the reduced density matrices $\rho_{I_j}$ in the learning algorithm. After each sweep we compute $\hat{F}_{\text{max}}(\rho, \sigma)$ between the target state $\rho$ and the output $\sigma$ using the testing set, and AFC with $k = \ell + 1$. We then choose as final output the MPO $\sigma$ which maximizes this
estimate of the fidelity. In what follows we present the results of a single numerical experiment. By repeating the numerical experiment, we do not observe significant differences in the results presented.

We run the learning algorithm for $\ell = 1, 2$ and $\chi' = 4$, which we have seen to be sufficient to achieve a fidelity error $1 - \mathcal{F} \sim 10^{-4}$. We report the results in Fig. S6. For the parameters considered, we always observe fidelity errors $1 - \mathcal{F} \gg 10^{-4}$, meaning that statistical errors are larger than the systematic errors due to finite $(\ell, \chi')$. We study both $\mathcal{F}_{\text{GM}}$ and $\mathcal{F}_{\text{max}}$, and compute them exactly (dotted lines) and using the testing set under AFC with $k$ (markers). These two estimates are in good agreement for $k > 1$.

First, we study the fidelity errors for $N_u^L = 1000$ as a function of number of qubits $N$, from $N = 16$ to $N = 128$. Similarly to the experiment reported in the main text, we observe a linear decay of both $\mathcal{F}_{\text{GM}}$ and $\mathcal{F}_{\text{max}}$, which hints at a linear propagation of statistical errors in the MPO (statistical errors in each tensor $M^{(j)}$ are independent of $N$, see App. C). For $(\ell, \chi') = (2, 4)$, $N_u^L = 1000$, we observe $\mathcal{F}_{\text{GM}} \sim 75\%$ and $\mathcal{F}_{\text{max}} \sim 70\%$.

For these parameters, we also compare physical properties between the MPOs $\rho$ and $\sigma$ for $N = 128$. In particular, we consider the purity $\mathcal{P}_2(\rho) = \text{tr}[\rho^2]$, and the two-point correlation function $\langle Y_j Y_{j+d} \rangle$, where $Y_j$ is a local Pauli operator. In both cases, we observe an almost perfect agreement between the predictions of $\rho$ and $\sigma$.

Finally, we study fidelity errors as a function of $N_u^L N_M$, with $N_u^L \in [1000, 5500]$, $N_M = 1024$. We observe a power-law decay $1 - \mathcal{F}_{\text{GM}} \propto (N_u^L N_M)^{-0.9}$ in the limit of large $N_u^L N_M$. While $\ell = 2$ already displays this behavior for $N_u^L = 1000$, the converse is not true for $\ell = 1$, which results in larger errors below a certain number of measurements. We note a qualitatively similar behavior for $\mathcal{F}_{\text{max}}$, although its behavior is not strictly monotonous in this regime.

# Appendix F: Common randomized measurements

In this section we give details on the use of common randomized measurements (CRM) [71] in the experimental demonstration of the learning algorithm. CRM allow one to reduce the overhead in number of bases $N_u$ when postprocessing randomized measurement data from a state $\rho$, given that a good classical approximation $\tau$ (also called prior) of $\rho$ is available. In what follows, we review the key elements of common randomized measurements, and explain how to implement them in the learning algorithm. Then, we describe the prior state used in the experimental demonstration, and compare performances with and without using CRM.

## 1. Common randomized measurements and tensor updates

Let us consider an experimental quantum state $\rho$ on which we have collected randomized measurements, which we represent as a set of classical shadows $\{\hat{\rho}^{(r)}\}$. In parallel, let us assume that we have access to a classical simulation of the experiment, which outputs the quantum state $\tau$. CRM consist in defining a new set of unbiased estimators $\{\hat{\rho}_{\tau}^{(r)}\}$ such that, if $||\rho - \tau||_2 \ll ||\rho||_2$, the statistical errors in this new set of estimators is much smaller than conventional classical shadows. From a classical representation of $\tau$, we can compute its classical shadows on a subsystem A for $N_M \to \infty$ in the same bases in which we have performed randomized measurements on $\rho$, obtaining

$$ \hat{\tau}_{A}^{(r)} = \sum_{s,s'} 2^{-D[s,s']} P_{\tau_A, U}(s) \bigotimes_{j \in A} \left( u_j^{(r)\dagger} |s'_j\rangle \langle s'_j| u_j^{(r)} \right), \quad (F1) $$

where $|s\rangle = \bigotimes_{j \in A} s_j$, $|s'\rangle = \bigotimes_{j \in A} s'_j$ are computational basis states on the subsystem A, $P_{\tau_A, U}(s) = \langle s|U^{(r)\dagger}_{\tau_A} U^{(r)}|s\rangle$ is their statistics for the state $\tau_A$ in a given random basis, with $U^{(r)} = \bigotimes_{j=1}^N u_j^{(r)}$, and $D[s, s']$ is the Hamming distance. At this point, the new set of estimators is defined as

$$ \hat{\rho}_{\tau,A}^{(r)} = \hat{\rho}_{A}^{(r)} - \hat{\tau}_{A}^{(r)} + \tau_A. \quad (F2) $$

Note that these new estimators are unbiased, as $\mathbb{E}[\hat{\rho}_{\tau,A}^{(r)}] = \mathbb{E}[\hat{\rho}_{A}^{(r)}] = \rho_A$, and $\hat{\rho}_{\tau,A}^{(r)} = \rho_A$ in the extreme case where $\rho = \tau$ and $N_M \to \infty$. To reach the latter regime, it is sufficient to take at least $N_M \gtrsim 2^{|A|}$ shots per basis, such that we obtain a good inference on the distribution $P_{\rho_A, U}(s)$. More concretely, when estimating the expectation value of a Pauli observable with support on A using measurements in $N_u$ random Pauli bases taking $N_M$ shots per basis, we obtain the variance [71]

$$ \text{Var}[\text{tr}[\hat{\rho}_{\tau,A} O]] \le \frac{3^{|A|}}{N_u} \left( \text{tr}[(\rho - \tau)O]^2 + \frac{1}{N_M} \right), \quad (F3) $$
<figure><img src="image_11.png" /><figcaption>FIG. S7. Estimated Rényi entropies of subsystems of size k on quantum simulation of the kicked Ising model at depth D versus the resulting local depolarization amplitude on the subsystem. We report results for the cases corresponding to $\ell = 1, 2$ in the learning algorithm for $D = 1$, and for $\ell = 3$ for $D = 2$.</figcaption></figure>

where $\hat{\rho}_{\tau,A} = \frac{1}{N_u} \sum_r (\hat{\rho}_A^{(r)} - \hat{\tau}_A^{(r)} + \tau_A)$, which results in a more precise estimation with respect to conventional classical shadows when $|\text{tr}[\rho_A O]| < |\text{tr}[(\rho_A - \tau_A)O]|$, given that $N_M$ is large enough.

We implement CRM in the learning algorithm in a similar way. First, we perform a MPO simulation of the reduced density matrices $\tau_{I_j}$ based on a prior knowledge of $\rho$, such as the quantum circuit preparing $\rho$ plus some details about the noise, or the Hamiltonian of the system when considering systems at finite temperature. Then, we modify the update rule in Eq. (A11) as

$$M^{(j)} \leftarrow C_{\sigma_{I_j}}^{(j)-1} (\hat{C}_{\rho_{I_j}}^{(j)} - \hat{C}_{\tau_{I_j}}^{(j)} + C_{\tau_{I_j}}^{(j)}), \quad M^{(j)} \leftarrow \frac{1}{\text{tr}[\sigma]} M^{(j)}, \qquad (\text{F4})$$

where

$$\hat{C}_{\rho_{I_j}}^{(j)} - \hat{C}_{\tau_{I_j}}^{(j)} + C_{\tau_{I_j}}^{(j)} = \text{tr}[(\hat{\rho}_{I_j} - \hat{\tau}_{I_j} + \tau)\partial_{M^{(j)}}\sigma] \qquad (\text{F5})$$

is the CRM estimator of $C_{\rho_{I_j}}^{(j)}$. While the analytical bound Eq. (F3) only works for Pauli observables, we expect the CRM tensor update Eq. (F4) to be less sensitive to statistical errors when $\tau \sim \rho$. We study its performances empirically in the next section on experimental data.

## 2. Experimental results

We start by describing the prior states $\{\tau_{I_j}\}$ used for CRM estimation from experimental data in the main text. We recall that the quantum circuit implemented in the experiment aims to realize the pure state

$$|\psi\rangle = \left( \prod_{j=1}^{N_{\text{tot}}-1} e^{i\frac{\pi}{4}Z_j Z_{j+1}} \prod_{j=1}^{N_{\text{tot}}} e^{-i\frac{\pi}{8}X_j} \right)^D |0\rangle^{\otimes N}, \qquad (\text{F6})$$

where $D = 1, 2$ is the depth of the circuit, $X_j$ and $Z_j$ are local Pauli operators. The effect of experimental noise effectively result in the preparation of an unknown mixed state $\rho$. For $D = 1, 2$ such states are exact MPS with bond dimension $\chi = 2, 4$. $\tau = |\psi\rangle\langle\psi|$ is not a useful prior state, since the condition $||\rho - \tau||_2 \ll ||\rho||_2$ cannot be satisfied for the values of purity $P_2(\rho)$ observed. In particular, the norm-2 distance is lower bounded by the difference of the purities, i.e. $||\rho - |\psi\rangle\rangle\langle\psi||_2 \ge |1 - P_2(\rho)| \approx 1 - 10^{-4}$ for $D = 1$.

We take instead reduced density matrices $\tau'_{I_j} = \text{tr}_{\bar{I}_j} [|\psi\rangle\langle\psi|]$, and add space-dependent noise using information collected from the experiment. For each reduced density matrix $\tau'_{[i,i+k]}$, we obtain the prior $\tau_{[i,i+k]}$ as

$$\tau_{[i,i+k]} = \left( \prod_{j \in [i, i+k]} \mathcal{E}_j(p_i^{(k)}) \right) [\tau_{[i,i+k]}], \qquad (\text{F7})$$

where each $\mathcal{E}_j(p)[\rho] = (1-p)\rho + p \text{tr}_j[\rho] \otimes \mathbb{I}_j$ is a local depolarizing channel with strength $p$. For each interval $[i, i+k]$ we consider a homogeneous depolarizing strength $p_i^{(k)}$, which we fit by comparing the resulting purity $\text{tr}[\tau_{[i,i+k]}^2]$ to
<figure><img src="image_12.png" /><figcaption>FIG. S8. Performances of the learning algorithm with the experimental data with and without common randomized measurements. (a-c) Estimated AFC fidelity of the output as a function of the number of random bases in the learning set $N_u^L$ for $N = 48, 96$ with (blue lines) and without (red lines) common randomized measurements. We consider $\ell = 1, 2$ and $\chi' = 4$ for $D = 1$, $\chi' = 8$ and $\ell = 3$ for $D = 2$, considering AFC fidelities with parameter $k = 3$ for $D = 1$ and $k = 4$ for $D = 2$. (d) Difference in the estimation of the expectation value $\langle Z_{j-1}X_jZ_{j+1} \rangle$ between randomized measurements and the MPO learned from the full quantum state ($N = 96$ qubits). We plot the result both using common randomized measurements in the learning algorithm (blue lines) and not using them (red lines). (e,f) Global purities as a function of $N$ for $D = 1, 2$, comparing the AFC estimation (dashed lines) for $k = 3$ ($k = 4$) for $D = 1$ ($D = 2$) with the results of the learning algorithm for $(\ell, \chi') = (2, 4)$ ($(\ell, \chi') = (3, 8)$) with (blue markers) and without common randomized measurements (red markers).</figcaption></figure>

$\text{tr}[\rho_{[i,i+k]}^2]$ as estimated from randomized measurement data. In this way we obtain a coarse-grained noise model for $\tau$. We plot in Fig. S7 the observed purities and the corresponding results for the fitting for $D = 1, 2$ and $k = 2(\ell + 1)$, with the values of $\ell$ considered in the main text.

In Fig. S8 we compare the results of the learning algorithm when analyzing the experimental data with and without CRM. In all the cases we investigate, we observe that CRM equipped with the noise model described above allow to reach higher fidelities, as expected. However, apart from the case $D = 1, \ell = 1$, the increase in fidelity is modest, i.e. comparable with roughly doubling $N_u^L$. A similar consideration can be drawn for purities, for which we observe similar predictions with and without CRM for the cases considered in the main text. Finally, we study the errors in estimating observables between our learning protocol and conventional randomized measurements. We define the error in estimating an observable using the learning algorithm as

$$ \delta\langle O \rangle = |\text{tr}[\hat{\rho}O] - \text{tr}[\sigma O]|, \qquad (F8) $$

where $\hat{\rho} = \frac{1}{N_u^T N_M} \sum_{r=1}^{N_u^T N_M} \hat{\rho}^{(r)}$ is the average classical shadow obtained from the testing set, and $\sigma$ is the MPO output by the learning algorithm. For $N_u^L = 512, N_u^T = 1536$, we observe that for $D = 1, \ell = 2$ the errors in estimating $\langle X_{j-1}Z_jX_{j+1} \rangle$ are generally mitigated, with the largest error as a function of $j$ being reduced by half. These results show how CRM can be useful to improve the predictive power of our learning algorithm. Finally, we remark that a different prior $\tau$ can give drastically different results. For instance, we expect that using a more fine-grained noise model (such as inhomogeneous local depolarization) would result in a better improvement of our results, at the cost of increased computational overhead to fit the model.