# Scalable Quantum State Tomography with Locally Purified Density Operators
and Local Measurements

Yuchen Guo  0000-0002-4901-2737 State Key Laboratory of Low Dimensional
Quantum Physics and Department of Physics, Tsinghua University, Beijing
100084, China  Shuo Yang  0000-0001-9733-8566<EMAIL_ADDRESS>State
Key Laboratory of Low Dimensional Quantum Physics and Department of Physics,
Tsinghua University, Beijing 100084, China  Frontier Science Center for
Quantum Information, Beijing 100084, China  Hefei National Laboratory, Hefei
230088, China

###### Abstract

Understanding quantum systems is of significant importance for assessing the
performance of quantum hardware and software, as well as exploring quantum
control and quantum sensing. An efficient representation of quantum states
enables realizing quantum state tomography with minimal measurements. In this
study, we propose a new approach to state tomography that uses tensor network
representations of mixed states through locally purified density operators and
employs a classical data postprocessing algorithm requiring only local
measurements. Through numerical simulations of one-dimensional pure and mixed
states and two-dimensional pure states up to size $8\times 8$, we demonstrate
the efficiency, accuracy, and robustness of our proposed methods. Experiments
on the IBM and Quafu Quantum platforms complement these numerical simulations.
Our study opens new avenues in quantum state tomography for two-dimensional
systems using tensor network formalism.

##  I Introduction

Quantum state tomography plays a fundamental role in characterizing and
evaluating the quality of quantum states produced by quantum devices. It
serves as a crucial element in the advancement of quantum hardware and
software, regardless of the underlying physical implementation and potential
applications [1, 2, 3]. However, reconstructing the full quantum state becomes
prohibitively expensive for large-scale quantum systems that exhibit potential
quantum advantages [4, 5], as the number of measurements required increases
exponentially with system size.

Recent protocols try to solve this challenge through two main steps: efficient
parameterization of quantum states and utilization of carefully designed
measurement schemes and classical data postprocessing algorithms. For one-
dimensional (1D) systems with area law entanglement, the matrix product state
(MPS) [6, 7, 8, 9, 10, 11, 12] provides a compressed representation. It
requires only a polynomial number of parameters that can be determined from
local or global measurement results. Two iterative algorithms using local
measurements, singular value thresholding (SVT) [13] and maximum likelihood
(ML) [14], have been demonstrated in trapped-ion quantum simulators with up to
$14$ qubits [15]. However, SVT is limited to pure states and thus impractical
for noisy intermediate-scale quantum (NISQ) systems. Meanwhile, although ML
can handle mixed states represented as matrix product operators (MPOs) [16,
17], it suffers from inefficient classical data postprocessing. Another scheme
reconstructs the quantum state by inverting local measurements, but the
resulting MPO is not necessarily positive [18].

On the other hand, some approaches also base on tensor networks (TN) but
update parameters from the sampling output of global measurements throughout
the system [19, 20]. It is noted that Ref. [21] provides a process tomography
method using global measurements, which can be applied directly to state
tomography due to the Choi-Jamiołkowski isomorphism [22, 23]. However, these
approaches cannot easily incorporate error mitigation techniques that focus
mainly on estimators rather than samplers [24, 25, 26, 27, 28, 29, 30, 31],
especially for readout errors [32, 33]. To scale to larger systems and
integrate with error mitigation techniques, TN-based methods using only local
measurements are preferred. However, TN state tomography remains unexplored
for higher-dimensional quantum systems.

Besides tensor network state tomography, neural networks emerge as strong
competitors for reconstructing unknown quantum states [34, 35, 36, 37, 38],
which can be expected from the deep connection between these two architectures
[39, 40, 41]. Another intriguing avenue is to estimate the observables without
reconstructing the entire density matrix but directly from its classical
shadows. Various techniques in this direction have been proposed for both
quantum states [42, 43, 44, 45] and quantum processes [46, 47].

In this work, we introduce a new approach to reconstructing mixed quantum
states using only local expectation values, for which readout error mitigation
is feasible. We propose representing mixed states as locally purified density
operators (LPDOs) [48] and optimizing their parameters using a variant of
gradient descent for a local loss function. To validate the effectiveness of
our method, we perform numerical simulations using typical quantum states: 1D
critical Ising ground states, 1D gapless spin-$\frac{1}{2}$ Heisenberg ground
states subjected to various noise models, two-dimensional (2D) random
projected entangled pair states (PEPSs) [49, 50, 51], 2D random isometric
tensor network states (isoTNSs) [52, 53, 54, 55], and 2D cluster states.
Furthermore, we implement experiments on real quantum devices accessible
through the IBM and Quafu Quantum platforms.

##  II Quantum state tomography with LPDO and local loss function

In this section, we introduce our new reconstruction approach for scalable
tomography of mixed states. We begin by efficiently parameterizing mixed
states using LPDOs, which have been shown to be efficient in simulating
thermal or dissipative many-body systems in 1D [56, 57].

In the LPDO form, a 1D mixed quantum state is expressed as follows

$\hat{\rho}=\sum_{\\{\bm{\tau},\bm{\omega}\\}}\sum_{\\{\bm{\mu},\bm{\nu},\bm{\kappa}\\}}\prod_{j=1}^{N}[A_{j}]^{\tau_{j},\kappa_{j}}_{\mu_{j-1},\mu_{j}}[A_{j}^{*}]^{\omega_{j},\kappa_{j}}_{\nu_{j-1},\nu_{j}}\ket{\tau_{1},\cdots,\tau_{N}}\bra{\omega_{1},\cdots,\omega_{N}},$
(1)

where $\bm{\mu}$ and $\bm{\nu}$ denote virtual indices that describe quantum
entanglement, and $\bm{\kappa}$ represents inner indices (also known as Kraus
indices) introduced for open systems [57], as depicted in Fig. 1(a).
Importantly, an LPDO is constructed to be Hermitian and positive semidefinite
by design. In scenarios with weak local noise, the Kraus dimension
$d_{\kappa}$ is typically a small constant, independent of the system size
$N$. This enables a significantly more efficient tomography approach compared
to directly reconstructing an MPO. Note that the quantum process tomography
method proposed in Ref. [21] has adopted the LPDO variant to represent quantum
channels. In addition, the ML method proposed in Ref. [14] could be adapted
and integrated into this framework by replacing MPO iterations with
computationally less expensive LPDO operations. However, we do not explore
this further as the exponential computational cost of iteration (shown later)
is not easily reduced.

Next, we propose a loss function constructed from local measurements and a
variant of gradient descent to update the local LPDO tensors. The loss
function for our problem is chosen as

$\Theta=\sum_{i}{||\hat{\sigma}_{\braket{i}}-\hat{\rho}_{\braket{i}}||_{F}^{2}}\equiv\sum_{i}\Theta_{\braket{i}},$
(2)

where $\hat{\sigma}_{\braket{i}}$ and $\hat{\rho}_{\braket{i}}$ are the
reduced density matrices for the sites $\\{i,\cdots,i+L-1\\}$ of the target
state obtained from experiments (see Methods) and the reconstructed state
respectively, as shown in Fig. 1(a). Expanding each term in the loss function
gives

$\Theta_{\braket{i}}={\rm
Tr}{\left[\hat{\sigma}_{\braket{i}}^{2}-2\hat{\sigma}_{\braket{i}}\hat{\rho}_{\braket{i}}+\hat{\rho}_{\braket{i}}^{2}\right]}.$
(3)

To calculate the gradients, terms such as

$\frac{\partial\Theta_{\braket{i}}}{\partial A_{j}^{*}}=2{\rm
Tr}{\left[\left(\hat{\rho}_{\braket{i}}-\hat{\sigma}_{\braket{i}}\right)\frac{\partial\hat{\rho}_{\braket{i}}}{\partial
A_{j}^{*}}\right]},$ (4)

need to be computed, requiring $O(N^{2})$ computational complexity.

However, instead of directly optimizing the entire loss function, we update
each local tensor $A_{j}$ only considering the adjacent terms $\Theta_{\langle
i\rangle}$ that involve the target site $j$ in the loss function.
Specifically, we update $A_{j}$ according to the following rule

$A_{j}\rightarrow
A_{j}-\eta\sum_{i=j-L+1}^{j}\frac{\partial\Theta_{\braket{i}}}{\partial
A_{j}^{*}},$ (5)

where $\eta$ is the learning rate, automatically adjusted using the Adam
optimizer [58]. Evaluating the gradient in each iterative step has a time
complexity of $O(ND^{3})$, where $D$ is the virtual bond dimension. This
approach converges to a high-quality approximation of the target state in only
$O(\log(N))$ iterative steps without encountering the issue of local minima or
barren plateaus, as demonstrated below. LPDO and the local loss function
together constitute our Grad-LPDO method. Importantly, for pure states, one
can set $d_{\kappa}=1$ and optimize the MPS manifold using the same gradient
method, which is treated as a special version of Grad-LPDO.

Here we briefly outline the relationship between our approach and previous
tensor network state tomography methods, which generally fall into two
categories. Our method follows the category of optimizing a local loss
function, encoding only local properties estimated from the measurements for
local observables [13, 14, 15]. The analytical theorem behind this scheme was
discussed in Ref. [13], i.e., for two pure quantum states $\ket{\phi}$
(target) and $\ket{\psi}$ (reconstructed), if $\ket{\psi}$ is the ground state
of a local Hamiltonian $H=\sum_{i}{h_{\braket{i}}}$ with an energy gap
$\Delta$, then the fidelity between these two states satisfies that

$f\geq 1-\frac{{\rm
Tr}{\left[\sum_{i}{\hat{h}_{\braket{i}}\left(\hat{\sigma}_{\braket{i}}-\hat{\rho}_{\braket{i}}\right)}\right]}}{\Delta}.$
(6)

In other words, once the experimental state can be approximated by an MPS,
only local information is sufficient to reconstruct the entire state. For a
general mixed state, though there is no rigorous proof, we believe that such
an argument also holds if the mixed state is generated from a noisy quantum
circuit, where the mixed state can be regarded as a pure state perturbed by
noise. The alternative approach regards the tensor network as a generative
model, using a global loss function to characterize the probability
distributions of samples generated by simultaneous measurements on all qubits.
These measurements involve random bases, and each sample is represented by an
integer-valued string that describes the measurement basis and spin
configuration [19, 20, 21]. The concrete realization and requirements of these
two formalisms will be further compared and discussed later.

##  III Numerical simulations in different systems

###  III.1 Comparison between Grad-LPDO and ML-MPS

Figure 1: Schematic of Grad-LPDO method and performance on 1D critical Ising
ground state with $D=D_{0}=16$ and varying $N$. (a) Mixed state with $N=4$ in
the LPDO form and the reduced density matrix for the sites $\\{i,i+1\\}$.
(b)-(f) Comparison of the ML-MPS and Grad-LPDO methods. (b)-(c) Number of
steps to achieve $f=0.9$ for (b) ML-MPS with $\log(n)$ inset, and (c) Grad-
LPDO with $\log(N)$ inset. (d) Maximal fidelity during iterations for both
methods. (e)-(f) Convergence stability of two methods for different learning
and update rates. (e) The ML-MPS method, where $\epsilon$ is the update rate
defined under Eq. (13). (f) Our Grad-LPDO method, where $\eta$ is the learning
rate in the Adam optimizer.

In this section, we present numerical demonstrations of our Grad-LPDO method
for both pure and mixed states of special interest. We begin by comparing the
performance of a standard approach also using local measurements, namely ML-
MPS [14], and our proposed Grad-LPDO method for reconstructing the ground
state of the 1D transverse field Ising model at the critical point $g=1$ under
open boundary condition (OBC). Both methods share a common initial step for
calculating local reduced density matrices from experimental results for
subsequent construction. Hence, in terms of experimental complexity, these two
methods are identical. The critical comparison therefore lies in the
complexity of classical data postprocessing algorithms.

The Hamiltonian is given by

$H=-\sum_{\braket{i,j}}Z_{i}Z_{j}+g\sum_{i}X_{i}.$ (7)

The ground state is approximated by an MPS with $D_{0}=16$ using the standard
variational method [59, 10], which serves as the target state for subsequent
reconstruction. Here, we assume that all measurements are ideal and directly
construct the reduced density matrices $\hat{\sigma}_{\braket{i}}$ from the
target state. Thus, any reconstruction errors are solely attributed to
insufficient local measurements and inefficient classical data postprocessing.
We perform both algorithms for different system sizes $N$, with $D=D_{0}=16$
and $L=2$, where the initial states are chosen as the paramagnetic ground
state for $g\rightarrow+\infty$. The hyperparameters in the Adam optimizer are
set as $\xi_{1}=\xi_{2}=0.8$ and $\epsilon=10^{-8}$. Fig. 1(b)-(c) shows the
number of iterative steps required to achieve fidelity $f=0.9$ for different
system sizes, representing the convergence speed of the algorithms.

The inset of Fig. 1(b) clearly shows that the time complexity of the ML-MPS
method scales exponentially with system size for a given reconstruction
accuracy. In contrast, the results shown in Fig. 1(c) are highly promising,
indicating that the number of gradient steps required for convergence scales
only as $O(\log{N})$. In other words, the overall time complexity of our Grad-
LPDO method to reconstruct a pure 1D state is $O(N\log(N)D^{3})$. This time
complexity significantly outperforms ML-MPS (and SVT-MPS with time complexity
$O(N^{4}D^{4})$ [13]), meaning that our method saves considerable time
compared to ML-MPS, especially for large systems. Even for small systems, ML-
MPS holds no efficiency advantage since truncating MPS in each iteration, as
required in the ML method (see Methods), is generally more computationally
expensive than simply contracting the environment in Grad-LPDO.

To provide a comprehensive comparison of these two methods, we calculate the
maximal fidelity achieved during iterations for different $N$ in Fig. 1(d),
which demonstrates the higher accuracy of our method. At the same time, we
also consider convergence stability when comparing different methods. Typical
iteration curves plotted in Fig. 1(e)-(f) reveal that Grad-LPDO consistently
converges to maximal fidelity, while the convergence of ML-MPS is less stable,
regardless of the update or learning rate per step. In conclusion, our Grad-
LPDO method surpasses the previous ML-MPS method in three key aspects:
efficiency, accuracy, and stability.

Figure 2: Comparison of different methods for one-dimensional spin-$1/2$
Heisenberg ground states. (a)-(c) Comparison of Grad-LPDO and Learn-LPDO. (a)
Different schemes for experiments and data postprocessing of two methods. (b)
Number of samples to achieve $f=0.95$ of two methods. (c) Classical runtime on
a single laptop CPU to achieve $f=0.95$ of two methods. (d) ZZ correlation
function for the target state and differences in the correlation between the
target state and the reconstructed states for three methods with $N=20$ and
$D=D_{0}=16$..

###  III.2 Comparison between Grad-LPDO and Learn-LPDO

Proceeding to the comparison with a recent process tomography method using
global measurements, referred to as Learn-LPDO [21], we illustrate the
experimental settings of these two formalisms in Fig. 2 (a). In Learn-LPDO,
each sample is generated from the circuit on a random measurement basis in a
single shot. On the other hand, in our method (as well as other local
information-based approaches [13, 14, 15]), each reduced density matrix of
length $L$ must be estimated from the complete local bases of size $4^{L}-1$,
providing $(N-L+1)(4^{L}-1)$ observables to be measured. Fortunately, not all
these measurements are independent, and the total number of independent bases
can be reduced to only $3^{L}$, regardless of the system size $N$ (see
Methods). If each basis is measured $m$ times to obtain an accurate
expectation value, the total number of experimental samples is $M=3^{L}m$.

In the following, we examine the performance of these two methods on the
ground state of the 1D gapless spin-$\frac{1}{2}$ Heisenberg model

$H=\sum_{\braket{i,j}}\mathbf{S}_{i}\cdot\mathbf{S}_{j}$ (8)

with $N=20$ and $D_{0}=16$. With such a model exhibiting non-trivial long-
range antiferromagnetic correlation, one might expect that the global Learn-
LPDO method would perform better than our Grad-LPDO method, since the former
takes into account global measurements while the latter depends only on local
marginals. However, the following results will show that this is not the case.

Fig. 2 (b) shows the number of total samples $M$ required for a given
reconstruction accuracy of $f=0.95$ for each method. Here, both sampling
processes are simulated, unlike in the previous section, to include
statistical errors from measurements. The result reveals that our Grad-LPDO
method has a similar scaling behavior of experimental complexity compared to
Learn-LPDO, with a constant coefficient difference of approximately $\sim
1/2$. In particular, even with the same number of samples, Learn-LPDO would be
more difficult to realize than Grad-LPDO, since all samples in Learn-LPDO are
generated from random bases that differ from each other, while the measurement
settings are much simpler in Grad-LPDO, involving only $3^{L}$ different
configurations. Therefore, our method proves to be slightly more efficient
than Learn-LPDO in terms of experimental complexity.

Figure 3: Results for gapless spin-$\frac{1}{2}$ Heisenberg ground states
($D=D_{0}=16$) without noise or undergoing various types of noise including
depolarizing (DP), bit flipping (BF), amplitude damping (AD), and phase
damping (PD) with error rate $\varepsilon=0.01$. (a-c) Conventional fidelity
$f_{0}$, our alternative fidelity $f$, and the loss function $\Theta$ during
the iteration process. (a) $N=10$ spins with amplitude damping noise,
reconstructed by LPDOs with $L=2$ and $d_{\kappa}=2$. (b) $N=20$ spins without
noise, reconstructed by LPDOs with $L=2$ and $d_{\kappa}=1$. (c) $N=20$ spins
with bit flipping noise, reconstructed by LPDOs with $L=2$ and $d_{\kappa}=2$.
(d-f) Maximal fidelity $f$ and purity $\mathcal{P}$ for different system sizes
$N$ and various types of noise. (d) $N=20$ spins with four types of noise,
reconstructed by LPDOs with $L=2$ and different $d_{\kappa}$. (e) $N=20$ spins
with four types of noise, reconstructed by LPDOs with $d_{\kappa}=2$ and
different $L$. (f) Depolarizing noise added to systems with different $N$,
reconstructed by LPDOs with $d_{\kappa}=2$ and different $L$.

Regarding the classical data postprocessing complexity of the reconstruction
process, we evaluate the CPU runtime required to converge to a fidelity of
$f=0.95$ serving as a figure of merit. Both methods are implemented on a
laptop CPU (Intel i7-11800H). The results shown in Fig. 2 (c) demonstrate a
substantial advantage of orders of magnitude for our Grad-LPDO method, where
the gap further increases with the system size $N$. The significant difference
in the classical runtime can be attributed to the complexity of evaluating
gradients in each iterative step. As discussed in the previous section, our
Grad-LPDO method uses samples from experiments to construct reduced density
matrices $\hat{\sigma}_{\braket{i}}$, then the loss function $\Theta$ in Eq.
(2) involves only $N-L+1$ terms of $\hat{\sigma}_{\braket{i}}$ instead of the
original samples. In essence, Grad-LPDO requires the contraction of $NL$
tensor networks when estimating all the gradients in Eq. (5). In contrast,
Learn-LPDO directly incorporates all samples into the loss function, requiring
the contraction of $MN$ tensor networks, as each sample contributes one tensor
network to one gradient. In summary, to achieve a comparable reconstruction
accuracy, the Learn-LPDO method proves to be more costly both experimentally
and computationally than our Grad-LPDO method. Meanwhile, in our method, each
complexity increases polynomially with system size $N$ rather than
exponentially, indicating the scalability of our method for larger systems.

To gain an intuitive understanding of the performance of the above three
methods, the ZZ correlation function
$C^{\rm{ZZ}}_{ij}=\braket{Z_{i}Z_{j}}-\braket{Z_{i}}\braket{Z_{j}}$ is
visualized for the target state in Fig. 2(d). We also show the differences in
the correlation between the target state and the reconstructed states with
$D=D_{0}=16$ and $L=2$. The result indicates that Grad-LPDO can capture most
of the long-range correlation and antiferromagnetic order from only local
measurements. On the other hand, ML-MPS cannot faithfully reproduce the
correlation between sites with intervals $L\geq 3$, which limits its
application to systems with non-trivial orders. For the Learn-LPDO method, the
overall accuracy is similar to that of Grad-LPDO, but with a much larger
fluctuation in the correlation function due to the intrinsic limitations of
any learning-based method.

###  III.3 Numerical simulations for 1D mixed states

Figure 4:  Generalization of Grad-LPDO to 2D PEPS. (a) Reconstruction fidelity
for random PEPSs with $D_{0}=3$ and different system sizes, with local
measurements of size $1\times 2$ and $2\times 1$ and simple gradient. (b)
Reconstruction fidelity for noisy random isoTNSs with $D=D_{0}=2$ and varying
system sizes, with local measurements of size $1\times L$ and $L\times 1$ and
simple gradients. (c) Reconstruction fidelity for noisy cluster states with
$D=D_{0}=2$ and varying system sizes, with local measurements of size $1\times
L$ and $L\times 1$ and different strategies to calculate gradients (simple or
full). All results are averaged over the $100$ realizations of the target
states.

Next, we consider Heisenberg ground states with $D_{0}=16$ that undergo
different types of local noise. Specifically, we add four types of single-
qubit noise with an equal error rate $\varepsilon=0.01$ to each qubit of the
ideal state, including depolarizing (DP), bit flipping (BF), amplitude damping
(AD), and phase damping (PD). Fig. 3 shows numerical results, where we
reconstruct target states using LPDOs with $D=16$ and different $d_{\kappa}$,
along with accurate local measurements of length $L=2$ without statistical
errors. The fidelity between two density matrices is defined as their inner
product in the operator space
$f\left(\hat{\rho}_{1},\hat{\rho}_{2}\right)\equiv{\rm
Tr}{\left(\hat{\rho}_{1}\hat{\rho}_{2}\right)}/\sqrt{{\rm
Tr}{\left(\hat{\rho}_{1}^{2}\right)}{\rm
Tr}{\left(\hat{\rho}_{2}^{2}\right)}}$, which can be efficiently calculated
for LPDOs (see Methods). To clarify the relationship between different figures
of merits, we compare the conventional definition of fidelity
$f_{0}\left(\hat{\rho}_{1},\hat{\rho}_{2}\right)=\left({\rm
Tr}{\sqrt{\sqrt{\hat{\rho}_{1}}\hat{\rho}_{2}\sqrt{\hat{\rho}_{1}}}}\right)^{2}$
(which is intractable for large systems), our alternative definition $f$, and
the loss function $\Theta$ throughout the iteration process in Fig. 3(a)-(c)
for different system sizes and with or without noise. The results show the
consistent convergence (blue dashed lines) of these indicators in different
situations. This not only confirms the validity of our alternative definition
$f$ but also highlights the sensitivity of the loss function $\Theta$ as a
reliable convergence criterion, especially considering the impracticality of
evaluating fidelity in real experiments.

Introducing noise to the target state will generally reduce the reconstruction
fidelity $f$, which cannot be improved by only increasing $d_{\kappa}$ as
implied in Fig. 3(d). This decrease arises because local measurements alone
cannot distinguish between mixtures in the reduced density matrix from
entanglement with unmeasured qubits (represented by virtual indices) and from
quantum noise (represented by Kraus indices), as confirmed by our numerical
simulations. For example, reconstructing quantum states with depolarizing or
bit-flipping noise poses greater challenges, since these noise models are
stochastic and directly introduce mixtures across different trajectories.
Furthermore, the purity of the target state $\mathcal{P}={\rm
Tr}[\hat{\rho}_{0}^{2}]$ is plotted in Fig. 3(e) for different noise types and
in Fig. 3(f) for depolarizing noise across varying system sizes $N$. They show
a clear dependence of the reconstruction fidelity $f$ on $\mathcal{P}$,
supporting our earlier arguments.

To overcome this challenge, we increase the measurement length $L$ in Fig.
3(e)-(f), significantly improving the accuracy for all types of noise. With
adjacent $4$-site measurements, the reconstruction fidelity $f$ exceeds 0.985
even for the most challenging depolarizing noise with up to $N=20$ qubits.
However, the experimental and data postprocessing costs scale exponentially
with $L$, requiring a trade-off between accuracy and efficiency. Our results
indicate that even for moderate-sized critical systems with typical noise
levels, a small constant $L$ suffices for high-fidelity reconstruction at an
acceptable cost using modern experimental and numerical techniques.

###  III.4 Generalization to 2D systems

Figure 5: Experiments on IBM Quantum platform. (a) Configuration of qubits and
noise information for device ‘ibm_nairobi’. (b) The iteration process of one
circuit realization ($\textrm{seed}=57$) for $d_{\kappa}=2$ and different $D$.
(c) Averaged residual loss $\Theta$ in all circuit realizations for different
$D$ and $d_{\kappa}$. (d)-(i) Expectation values for $L=2$ local observables
obtained from experiments and reconstruction for $d_{\kappa}=2$ and $D=4$.

We now generalize our Grad-LPDO method to pure states of two-dimensional
systems on square lattices, represented as PEPSs. The key idea and procedure
remain the same, where the local measurements in Eq. (2) are applied to the
$L_{1}\times L_{2}$ subsystems. When updating the corresponding local tensor
$A_{j}$, gradients in Eq. (5) contain terms with local measurements covering
the site $j$, which we denote as simple gradients. In general, the contraction
of a 2D TN is computationally expensive. Here, we adopt the standard
truncation method for finite-size systems [60, 10]. The 2D TN with bond
dimension $D$ is contracted layer by layer as the evolution of a 1D MPS,
truncating the bond dimension to $\chi=D^{2}$ after each contraction.

We first simulate random PEPS (i.e., all tensor elements are random complex
numbers) as target states with fixed $D_{0}=3$ and varying system sizes
$N\times N$, where random PEPSs with different $D$ are chosen as initial
states for the iteration. We directly construct $\hat{\sigma}_{\braket{i}}$ on
all $1\times 2$ and $2\times 1$ subsystems of the target states. Fig. 4 (a)
shows the reconstruction fidelity of the $100$ shots for each pair of $N$ and
$D$, with error bars giving the standard deviation. The average fidelity
reaches $0.995$ for random $8\times 8$ PEPSs with $D_{0}=3$, which are
generally noncritical.

However, a PEPS with random tensor elements is generally slightly entangled
and thus does not represent a good 2D generalization. Therefore, we proceed to
two more specific types of 2D system, including the isometric tensor network
(isoTNS) [52] and the cluster state. The former has found wide applications in
topological states [53, 55], quantum circuits [61, 62], and thermal states
[54], while the latter is a typical state with symmetry protected topological
order [63, 64, 65, 66, 67, 68] and serves as a resource state for measurement-
based quantum computation [69, 70, 71, 72, 73] or measurement-induced phase
transitions [74, 75]. As in real experiments for benchmarking quantum devices
or error mitigation, one knows the ideal quantum state representation and
tries to learn the noisy quantum state, we can choose the ideal state as the
initial ansatz for iteration to speed up the convergence. To be more specific,
we first construct an ideal PEPS $\ket{\Psi_{\rm Ideal}}$ and add perturbation
to each local tensor as an effective noise model, where the output state
$\ket{\Psi_{\rm Target}}=\ket{\Psi_{\rm Noisy}}$ is set as the target state
for reconstruction. Meanwhile, the ideal state is chosen as the initial state
for iteration $\ket{\Psi_{0}}=\ket{\Psi_{\rm Ideal}}$.

In Fig. 4 (b), we simulate our tomography method on random isoTNSs, where each
local tensor corresponds to a random haar unitary. In other words, the entire
state can be mapped to a sequential quantum circuit with random quantum gates
[76, 61], whose depth will scale linearly with the system size given only a
constant bond dimension $D$, and therefore has significant expressive power
for complex quantum systems and serves as a strong candidate for benchmarking
our approach. Reconstruction fidelities with different local measurement sizes
$1\times L$ and $L\times 1$ for $D=D_{0}=2$ and varying system sizes are
compared, where the initial fidelity $f_{0}\equiv f(\ket{\Psi_{\rm
Target}},\ket{\Psi_{\rm 0}})=f(\ket{\Psi_{\rm Noisy}},\ket{\Psi_{\rm Ideal}})$
is plotted for benchmarking the ‘noise strength’. The result demonstrates the
effectiveness of our method for typical 2D systems up to $8\times 8$ of
interest to capture most of the noise effects with only local measurements of
size $1\times 3$ and $3\times 1$. In addition, we also compare the simple
gradients and the full gradients (including all terms in the loss function
(2)) in Fig. 4 (c) for 2D cluster states, where full gradients are slightly
more accurate than simple gradients, indicating the effectiveness of only
focusing on local neighbors when implementing classical data postprocessing.

##  IV Real experiments with cloud quantum computation

###  IV.1 Experiments on IBM Quantum platform

Figure 6: Experiments on Quafu Quantum platform. (a) Configuration of qubits
and noise information for device ‘ScQ-P10’. (b) The iteration process for
$d_{\kappa}=2$ and different $D$. (c) Averaged residual loss $\Theta$ for
different $D$ and $d_{\kappa}$. (d) Circuit configuration to generate an
antiferromagnetic GHZ state with $N=10$. (e) Values of the ZZ correlation
function for the target state and the reconstructed state with $D=4$ and
$d_{\kappa}=2$.

To demonstrate our method on real quantum hardware, we conduct experiments on
the IBM Quantum platform. Specifically, we use the ‘ibm_nairobi’ quantum
computer, a superconducting processor with $N=7$ available qubits, whose qubit
configuration and noise information are shown in Fig. 5(a). The input state is
the trivial product state $\ket{\psi_{0}}=\ket{0}^{\otimes N}$, then random
Haar circuits are implemented to generate the output states to be
reconstructed. $3^{L}$ ($L=2$ here) numbers of Pauli strings are measured to
estimate each reduced density matrix $\hat{\sigma}_{i}$ for the target state,
and each observable is measured using $m=10,000$ shots. To mitigate errors in
measurement results, the twirled readout error extinction (T-REx) [32] is
employed. Other error mitigation techniques for circuit errors are not
considered, as one of the main objectives of tomography is to learn about the
noise information in the system. Since the target state is unknown in
practice, we use the loss function $\Theta$ defined in Eq. (2) to represent
the residual error for tomography, which is highly related to the final
reconstruction fidelity that cannot be directly evaluated in practical
scenarios [13].

In Fig. 5(b), the iteration process of one typical circuit realization (seed
$57$ for the random circuit) is plotted for $d_{\kappa}=2$ and different $D$.
The experimentally measured expectation values ${\rm
Tr}{[\hat{\sigma}_{\braket{i}}\hat{P}_{\braket{i}}]}$ and the corresponding
estimated values from the reconstructed states ${\rm
Tr}{[\hat{\rho}_{\braket{i}}\hat{P}_{\braket{i}}]}$ are calculated for each
reduced density matrix with $d_{\kappa}=2$ and $D=4$, where
$\hat{P}_{\braket{i}}$ refers to Pauli strings. These results are visualized
in Fig. 5(d)-(i) to demonstrate their high consistency. Fig. 5(c) shows the
residual error averaged over $64$ random circuit realizations along with the
standard deviation. For each circuit realization, we choose $100$ random LPDOs
as initial states and record the average residual error. These results confirm
that LPDOs with small $d_{\kappa}$ and $D$ serve as good approximations for
quantum states generated from noisy quantum circuits, validating the
performance of our tomography scheme.

###  IV.2 Experiments on Quafu Quantum platform

Besides random quantum circuits, we also apply our approach to test on a
prominent quantum state — the Greenberger-Horne-Zeilinger (GHZ) state
characterized by long-range entanglement, which is of central importance in
both quantum information and condensed matter physics. For this evaluation, we
choose the ‘ScQ-P10’ superconducting quantum chip with $N=10$ qubits available
on the Quafu Quantum platform [77]. The qubit configuration and noise
information are shown in Fig. 6(a), and the quantum circuit designed to
generate a GHZ state is shown in Fig. 6 (d). We add extra $X$ gates on even
qubits to create an antiferromagnetic order, which improves clarity in the
subsequent visualization of the correlation function. For this scenario, we
set $L=2$ and $m=10000$, i.e., $3^{L}$ Pauli strings, each measured with
10,000 shots.

The iteration process shown in Fig. 6 (b) demonstrates the efficient
convergence of our algorithm on the quantum state with long-range
entanglement, which is in principle more difficult to handle than previously
considered random states and short-range entangled states. The residual errors
for different $D$ and $d_{\kappa}$ are plotted in Fig. 6 (c), indicating a
convergence of accuracy at $D=4$ and $d_{\kappa}=2$, greater than those for an
ideal GHZ state ($D=2$ and $d_{\kappa}=1$), implying the correlated noise
effects on a noisy quantum chip in the generated quantum state. Furthermore,
the ZZ correlation function $C_{ij}^{ZZ}$ measured directly from the
experiments and that estimated from the reconstructed state are compared in
Fig. 6 (d). Their high consistency, especially the ability to reveal the
failure of qubit 0 to correlate with other qubits, not only validates our
reconstruction method but also suggests its potential to identify and mitigate
errors in noisy quantum circuits [31, 29].

##  V Discussion

In this study, we introduce a new state tomography scheme that uses the LPDO
parameterization of mixed states and classical data postprocessing to
determine unknown tensors from only local measurements. Our approach
demonstrates enhanced efficiency, accuracy, and robustness compared to
previous methods relying on MPS representations. In particular, we extend TN
state tomography, originally developed for 1D systems, to higher spatial
dimensions. Our optimization method demonstrates efficacy in alleviating the
curse of dimensionality for high-dimensional systems.

The findings suggest that the recently proposed process tomography method [21]
and the error mitigation approach [29] may generalize to higher-dimensional
circuits, enabling a deeper understanding of noise and its effects in such
systems [78]. Our protocol facilitates the realization of quantum state
tomography for large systems with potential quantum advantages, promoting
advances in precise quantum control and complex algorithm implementation [2,
79].

##  VI Methods

###  VI.1 Local measurements and reduced density matrices

To obtain the reduced density matrix $\hat{\sigma}_{\braket{i}}$ for the sites
$\\{i,\cdots,i+L-1\\}$ in experiments, one needs to implement an
informationally complete set of measurements in this subsystem. Specifically,
we consider all possible products of Pauli operators
$\hat{P}_{\braket{i}}^{\bm{m}}$ acting on these adjacent $L$ sites, where
$\bm{m}=\\{m_{i},\cdots,m_{i+L-1}\\}$ with $m_{i}\in\\{I,X,Y,Z\\}$. Then
$\hat{\sigma}_{\braket{i}}$ can be expanded as [1, 13]

$\hat{\sigma}_{\braket{i}}=\frac{1}{2^{L}}\sum_{\bm{m}}{{\rm
Tr}{\left[\hat{\sigma}_{\braket{i}}\hat{P}_{\braket{i}}^{\bm{m}}\right]}\hat{P}_{\braket{i}}^{\bm{m}}}$
(9)

since $\hat{P}_{\braket{i}}^{\bm{m}}$ constitute a set of complete and
orthogonal bases in the operator space. To reconstruct each
$\hat{\sigma}_{\braket{i}}$, the expectation values of all $4^{L}-1$ Pauli
strings must be measured (as the identity operator is trivial). In other
words, a total of $(4^{L}-1)(N-L+1)$ observables are required to implement the
tomography.

To reduce the number of different measurement settings from $(4^{L}-1)(N-L+1)$
to only $3^{L}$, we use several strategies. First, any Pauli string containing
$I$ can be reduced since its expectation value can be estimated from other
Pauli strings, where $I$ is replaced by one of $\\{X,Y,Z\\}$. In other words,
only $3^{L}$ measurement bases are sufficient to reconstruct the reduced
density matrix for a local subsystem of length $L$. Second, the factor
$(N-L+1)$ is effectively eliminated because a single measurement configuration
with multiple shots allows the simultaneous extraction of several observable
expectation values. For example, measuring $Z$ for all qubits $Z_{1}\otimes
Z_{2}\otimes\cdots\otimes Z_{N}$ enables the determination of all local
observables in the form of $Z_{i}\otimes Z_{i+1}\otimes\cdots\otimes
Z_{i+L-1}$ (product-$Z$ strings of length $L$) with sufficient samples. In
fact, all translation-invariant Pauli strings
$\hat{P}_{\braket{1}}^{\bm{m}}\otimes\hat{P}_{\braket{L+1}}^{\bm{m}}\otimes\cdots\hat{P}_{\braket{N-L+1}}^{\bm{m}}$
with a period of $L$ sites, assuming that $N/L$ is an integer, are adequate to
reduce to any local Pauli string of length $L$ after tracing out the
corresponding environment. Here, $\bm{m}={m_{1},\cdots,m_{L}}$ and
$m_{i}\in\\{X,Y,Z\\}$. In essence, only $3^{L}$ measurement bases, each
represented by a translation-invariant global Pauli string with a period of
$L$, are required to compute all the necessary reduced density matrices for
our tomography method.

###  VI.2 SVT-MPS method and ML-MPS method

We briefly review two QST methods for 1D systems based on MPS and local
measurements [13, 14].

The SVT-MPS method [13] is inspired by the SVT algorithm in computer science
for matrix completion [80]. The target state can be approached iteratively by
solving a local Hamiltonian at each step. Specifically, in the $n$-th
iterative step we construct and find the dominant eigenstate of the following
local Hamiltonian

$\hat{Y}_{n+1}=\hat{Y}_{n}+\delta_{n}\left(\sum_{i}\hat{\sigma}_{\braket{i}}-E_{n}\sum_{i}\hat{\rho}_{n\braket{i}}\right).$
(10)

Here, $\delta_{n}$ is the ‘update rate’ for each step,
$\hat{\sigma}_{\braket{i}}$ are reduced density matrices for adjacent $L$
sites of the target state obtained through direct measurements, while
$\hat{\rho}_{n\braket{i}}$ are reduced density matrices of the dominant
eigenstate of $\hat{Y}_{n}$ with eigenvalue $E_{n}$. The authors have shown
that the number of iterative steps to achieve fixed fidelity typically scales
as $O(N^{2})$, where $N$ is the number of qubits. Additionally, in each step,
a variational method is needed to solve a local Hamiltonian, involving $O(N)$
sweep back and forth and $O(N)$ calculations of the environment per sweep.
Updating local tensors for each site also requires calculating the dominant
eigenvector of a $D^{2}d_{p}\times D^{2}d_{p}$ matrix, with complexity at
least $O(D^{4})$ assuming sparsity. Therefore, the total computational cost of
SVT-MPS is $O(N^{4}D^{4})$, which is scalable but limited to pure states.
However, since real experiments encounter mixed states with decoherence, this
reconstruction scheme has limited practicality for characterizing the noise
effect in quantum devices.

The ML-MPS method directly searches for the target state that maximizes the
log-likelihood function

$\log{\mathcal{L}\left(\hat{\rho}\right)}=\sum_{i,j}n_{\braket{i}}^{j}\log{\left({\rm
Tr}\left[\hat{\Pi}_{\braket{i}}^{j}\hat{\rho}\right]\right)}$ (11)

via a fixed-point iterative algorithm. Here, $\hat{\Pi}_{\braket{i}}^{j}$ are
local projectors labeled with $j$ applied at adjacent $L$ sites
$\\{i,\cdots,i+L-1\\}$, and $n_{\braket{i}}^{j}$ are the corresponding
measurement outcomes. The solution $\hat{\rho}_{\textrm{ML}}$ that maximizes
the above function satisfies

$\hat{\rho}_{\textrm{ML}}=\frac{1}{M}\sum_{i,j}\frac{n_{\braket{i}}^{j}}{{\rm
Tr}{\left[\hat{\Pi}_{\braket{i}}^{j}\hat{\rho}_{\text{ML}}\right]}}\hat{\Pi}_{\braket{i}}^{j}\hat{\rho}_{\text{ML}}\equiv\mathcal{R}(\hat{\rho}_{\text{ML}})\hat{\rho}_{\text{ML}},$
(12)

which corresponds to the fixed-point equation

$\hat{\rho}=\mathcal{R}(\hat{\rho})\hat{\rho}\mathcal{R}(\hat{\rho}).$ (13)

In practice, one can replace $\mathcal{R}$ by
$(\mathcal{I}+\epsilon\mathcal{R})/(1+\epsilon)$ with $\epsilon\ll 1$.
Furthermore, under the assumption of pure state, we only need to iterate on
the pure state manifold $\ket{\psi}=\mathcal{R}\ket{\psi}$. To implement this,
we construct the MPO representation of $\mathcal{R}$ in each iteration and
truncate the resulting MPS $\mathcal{R}\ket{\psi}$. Truncation can be done
variationally by minimizing the error
$e=(\bra{\psi}-\bra{\psi^{\prime}})(\ket{\psi}-\ket{\psi^{\prime}})$ or more
efficiently by using SVD from site to site in canonical form, which requires
$O(ND^{3})$ operations.

In the following table, we summarize the computational complexity and
application scope of the previous two methods and our Grad-LPDO method.

Method | Complexity | Application
---|---|---
SVT-MPS | $O(N^{4}D^{4})$ | Pure states
ML-MPS (MPO) | $O(\exp(N)ND^{3})$ | Pure (Mixed) states
Grad-LPDO | $O(\log(N)ND^{3})$ | Mixed states

###  VI.3 Learn-LPDO method

Here we discuss how to apply the process tomography method in Ref. [21] to the
task of state tomography with LPDO. The measurement operators are constructed
as $\Pi_{\bm{\beta}}=\Pi_{\beta_{1}}\otimes\cdots\otimes\Pi_{\beta_{N}}$,
where $\beta_{i}=1,2,\cdots,6$, corresponding to $6$ projectors of one
measurement randomly chosen from $\\{X,Y,Z\\}$ and possible outcomes. These
operators satisfy the condition $\sum_{\beta_{i}}{\Pi_{\beta_{i}}}=I$. For any
quantum state $\rho$, the probability distribution of such global random
measurements is

$P_{\rho}(\bm{\beta})={\rm Tr}{\left[\rho\Pi_{\bm{\beta}}\right]}.$ (14)

Using the Kullbach-Leibler divergence between the target and reconstructed
states

$D_{\rm{KL}}=\sum_{\\{\bm{\beta}\\}}P_{\rm{Exp}}(\bm{\beta})\log{\frac{P_{\rm{Exp}}(\bm{\beta})}{P_{\rm{LPDO}}(\bm{\beta})}},$
(15)

we establish the loss function to minimize as

$\Theta=-\frac{1}{M}\sum_{k=1}^{M}\log{P_{\rm{LPDO}}(\bm{\beta}_{k})},$ (16)

where $k$ is the label of different samples with a total number of $M$.

The LPDO parameters are updated through gradient descent $A_{j}\rightarrow
A_{j}-\eta\frac{\partial\Theta}{\partial A_{j}^{*}}$. The complexity involved
in calculating all gradients in one iterative step is $O(NM)$ since the loss
function contains all original samples. In our numerical simulations for
Learn-LPDO, we use the Julia package PastaQ, which provides built-in functions
to facilitate the sampling process and the implementation of this method.

###  VI.4 Fidelity between two states

The fidelity between two normalized pure states is defined as

$f\left(\ket{\psi},\ket{\phi}\right)=\left|\braket{\psi}{\phi}\right|^{2}.$
(17)

This is usually generalized for two mixed states as

$f_{0}\left(\hat{\rho}_{1},\hat{\rho}_{2}\right)=\left({\rm
Tr}{\sqrt{\sqrt{\hat{\rho}_{1}}\hat{\rho}_{2}\sqrt{\hat{\rho}_{1}}}}\right)^{2}$
(18)

with normalized ${\rm Tr}{\left[\hat{\rho}_{1}\right]}={\rm
Tr}{\left[\hat{\rho}_{2}\right]}=1$. However, this definition cannot be
directly estimated for two mixed states in their LPDO form, which prevents
direct benchmarking of our tomography method for mixed states in large
systems. Therefore, we adopt an alternative definition

$f\left(\hat{\rho}_{1},\hat{\rho}_{2}\right)\equiv{\rm
Tr}{\left(\hat{\rho}_{1}\hat{\rho}_{2}\right)}/\sqrt{{\rm
Tr}{\left(\hat{\rho}_{1}^{2}\right)}{\rm
Tr}{\left(\hat{\rho}_{2}^{2}\right)}},$ (19)

which is the inner product in the operator space and equals the overlap
between two superoperators $|\rho_{1}\rangle\\!\rangle$ and
$|\rho_{2}\rangle\\!\rangle$. In particular, this alternative definition of
fidelity reduces to Eq. (17) for pure states.

###  VI.5 Noise models

In our numerical simulations for mixed states, the noise added after each
state includes four types. The depolarizing noise is defined as

$\mathcal{E}\left(\hat{\rho}\right)=\left(1-\frac{4}{3}\varepsilon\right)\hat{\rho}+\frac{1}{3}\varepsilon\sum_{i=0}^{3}\sigma_{i}\hat{\rho}\sigma_{i}.$
(20)

The bit flipping noise is defined as

$\mathcal{E}\left(\hat{\rho}\right)=\left(1-\varepsilon\right)\hat{\rho}+\varepsilon\sigma_{x}\hat{\rho}\sigma_{x}.$
(21)

The amplitude damping noise is defined by the Kraus operator
$E_{0}=\ket{0}\bra{0}+\sqrt{1-\varepsilon}\ket{1}\bra{1}$ and
$E_{1}=\sqrt{\varepsilon}\ket{0}\bra{1}$ with the operator-sum representation

$\mathcal{E}\left(\hat{\rho}\right)=E_{0}\hat{\rho}E_{0}^{\dagger}+E_{1}\hat{\rho}E_{1}^{\dagger}.$
(22)

Phase damping noise is defined similarly, with the Kraus operator
$E_{0}=\ket{0}\bra{0}+\sqrt{1-\varepsilon}\ket{1}\bra{1}$ and
$E_{1}=\sqrt{\varepsilon}\ket{1}\bra{1}$.

##  VII Data availability

The datasets generated and analyzed during the current study are available
from the corresponding author upon reasonable request.

## References

  * Nielsen and Chuang [2009] M. A. Nielsen and I. L. Chuang,  _Quantum Computation and Quantum Information_ (Cambridge University Press, 2009).
  * Preskill [2018] J. Preskill, Quantum computing in the NISQ era and beyond, Quantum 2, 79 (2018).
  * Gebhart  _et al._ [2023] V. Gebhart, R. Santagati, A. A. Gentile, E. M. Gauger, D. Craig, N. Ares, L. Banchi, F. Marquardt, L. Pezzè, and C. Bonato, Learning quantum systems, Nat. Rev. Phys. 5, 141 (2023).
  * Arute  _et al._ [2019] F. Arute, K. Arya, R. Babbush, D. Bacon, _et al._ , Quantum supremacy using a programmable superconducting processor, Nature 574, 505 (2019).
  * Kim  _et al._ [2023] 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, Evidence for the utility of quantum computing before fault tolerance, Nature 618, 500 (2023).
  * Verstraete and Cirac [2006] F. Verstraete and J. I. Cirac, Matrix product states represent ground states faithfully, Phys. Rev. B 73, 094423 (2006).
  * Pérez-García  _et al._ [2007] D. Pérez-García, F. Verstraete, M. M. Wolf, and J. I. Cirac, Matrix product state representations, Quantum Info. Comput. 7, 401–430 (2007).
  * Verstraete  _et al._ [2008] F. Verstraete, V. Murg, and J. Cirac, Matrix product states, projected entangled pair states, and variational renormalization group methods for quantum spin systems, Adv. Phys. 57, 143 (2008).
  * Schollwöck [2011] U. Schollwöck, The density-matrix renormalization group in the age of matrix product states, Ann. Phys. 326, 96 (2011).
  * Orús [2014] R. Orús, A practical introduction to tensor networks: Matrix product states and projected entangled pair states, Ann. Phys. 349, 117 (2014).
  * Cirac  _et al._ [2021] J. I. Cirac, D. Pérez-García, N. Schuch, and F. Verstraete, Matrix product states and projected entangled pair states: Concepts, symmetries, theorems, Rev. Mod. Phys. 93, 045003 (2021).
  * Verstraete  _et al._ [2023] F. Verstraete, T. Nishino, U. Schollwöck, M. C. Bañuls, G. K. Chan, and M. E. Stoudenmire, Density matrix renormalization group, 30 years on, Nat. Rev. Phys. 5, 273 (2023).
  * Cramer  _et al._ [2010] M. Cramer, M. B. Plenio, S. T. Flammia, R. Somma, D. Gross, S. D. Bartlett, O. Landon-Cardinal, D. Poulin, and Y.-K. Liu, Efficient quantum state tomography, Nat. Commun. 1, 149 (2010).
  * Baumgratz  _et al._ [2013a] T. Baumgratz, A. Nüßeler, M. Cramer, and M. B. Plenio, A scalable maximum likelihood method for quantum state tomography, New J. Phys. 15, 125004 (2013a).
  * Lanyon  _et al._ [2017] B. P. Lanyon, C. Maier, M. Holzäpfel, T. Baumgratz, C. Hempel, P. Jurcevic, I. Dhand, A. S. Buyskikh, A. J. Daley, M. Cramer, M. B. Plenio, R. Blatt, and C. F. Roos, Efficient tomography of a quantum many-body system, Nat. Phys. 13, 1158 (2017).
  * Pirvu  _et al._ [2010] B. Pirvu, V. Murg, J. I. Cirac, and F. Verstraete, Matrix product operator representations, New J. Phys. 12, 025012 (2010).
  * Guth Jarkovský  _et al._ [2020] J. Guth Jarkovský, A. Molnár, N. Schuch, and J. I. Cirac, Efficient description of many-body systems with matrix product density operators, PRX Quantum 1, 010304 (2020).
  * Baumgratz  _et al._ [2013b] T. Baumgratz, D. Gross, M. Cramer, and M. B. Plenio, Scalable reconstruction of density matrices, Phys. Rev. Lett. 111, 020401 (2013b).
  * Han  _et al._ [2018] Z.-Y. Han, J. Wang, H. Fan, L. Wang, and P. Zhang, Unsupervised generative modeling using matrix product states, Phys. Rev. X 8, 031012 (2018).
  * Wang  _et al._ [2020] J. Wang, Z.-Y. Han, S.-B. Wang, Z. Li, L.-Z. Mu, H. Fan, and L. Wang, Scalable quantum tomography with fidelity estimation, Phys. Rev. A 101, 032321 (2020).
  * Torlai  _et al._ [2023] G. Torlai, C. J. Wood, A. Acharya, G. Carleo, J. Carrasquilla, and L. Aolita, Quantum process tomography with unsupervised learning and tensor networks, Nat. Commun. 14, 2858 (2023).
  * Jamiołkowski [1972] A. Jamiołkowski, Linear transformations which preserve trace and positive semidefiniteness of operators, Rep. Math. Phys. 3, 275 (1972).
  * Choi [1975] M.-D. Choi, Completely positive linear maps on complex matrices, Linear Alg. Appl. 10, 285 (1975).
  * Li and Benjamin [2017] Y. Li and S. C. Benjamin, Efficient variational quantum simulator incorporating active error minimization, Phys. Rev. X 7, 021050 (2017).
  * Temme  _et al._ [2017] K. Temme, S. Bravyi, and J. M. Gambetta, Error mitigation for short-depth quantum circuits, Phys. Rev. Lett. 119, 180509 (2017).
  * Endo  _et al._ [2018] S. Endo, S. C. Benjamin, and Y. Li, Practical quantum error mitigation for near-future applications, Phys. Rev. X 8, 031027 (2018).
  * McArdle  _et al._ [2019] S. McArdle, X. Yuan, and S. Benjamin, Error-mitigated digital quantum simulation, Phys. Rev. Lett. 122, 180501 (2019).
  * Cai [2021] Z. Cai, Multi-exponential error extrapolation and combining error mitigation techniques for NISQ applications, npj Quantum Inform. 7, 80 (2021).
  * Guo and Yang [2022] Y. Guo and S. Yang, Quantum error mitigation via matrix product operators, PRX Quantum 3, 040313 (2022).
  * van den Berg  _et al._ [2023] E. van den Berg, Z. K. Minev, A. Kandala, and K. Temme, Probabilistic error cancellation with sparse pauli–lindblad models on noisy quantum processors, Nat. Phys. 19, 1116 (2023).
  * Cai  _et al._ [2023] Z. Cai, R. Babbush, S. C. Benjamin, S. Endo, W. J. Huggins, Y. Li, J. R. McClean, and T. E. O’Brien, Quantum error mitigation (2023), arXiv:2210.00921 .
  * van den Berg  _et al._ [2022] E. van den Berg, Z. K. Minev, and K. Temme, Model-free readout-error mitigation for quantum expectation values, Phys. Rev. A 105, 032620 (2022).
  * Yang  _et al._ [2022] B. Yang, R. Raymond, and S. Uno, Efficient quantum readout-error mitigation for sparse measurement outcomes of near-term quantum devices, Phys. Rev. A 106, 012423 (2022).
  * Torlai  _et al._ [2018] G. Torlai, G. Mazzola, J. Carrasquilla, M. Troyer, R. Melko, and G. Carleo, Neural-network quantum state tomography, Nat. Phys. 14, 447 (2018).
  * Carrasquilla  _et al._ [2019] J. Carrasquilla, G. Torlai, R. G. Melko, and L. Aolita, Reconstructing quantum states with generative models, Nat. Mach. Intell. 1, 155–161 (2019).
  * Torlai and Melko [2020] G. Torlai and R. G. Melko, Machine-learning quantum states in the nisq era, Annu. Rev. Condens. Matter Phys. 11, 325 (2020).
  * Golubeva and Melko [2022] A. Golubeva and R. G. Melko, Pruning a restricted boltzmann machine for quantum state reconstruction, Phys. Rev. B 105, 125124 (2022).
  * Zhu  _et al._ [2022] Y. Zhu, Y.-D. Wu, G. Bai, D.-S. Wang, Y. Wang, and G. Chiribella, Flexible learning of quantum states with generative query neural networks, Nat. Commun. 13, 6222 (2022).
  * Chen  _et al._ [2018] J. Chen, S. Cheng, H. Xie, L. Wang, and T. Xiang, Equivalence of restricted boltzmann machines and tensor network states, Phys. Rev. B 97, 085104 (2018).
  * Huang and Moore [2021] Y. Huang and J. E. Moore, Neural network representation of tensor network and chiral states, Phys. Rev. Lett. 127, 170601 (2021).
  * Wang  _et al._ [2023] M. Wang, Y. Pan, Z. Xu, X. Yang, G. Li, and A. Cichocki, Tensor networks meet neural networks: A survey and future perspectives (2023), arXiv:2302.09019 .
  * Aaronson [2018] S. Aaronson, Shadow tomography of quantum states, in  _STOC 2018_ (2018) p. 325–338.
  * Huang  _et al._ [2020] H.-Y. Huang, R. Kueng, and J. Preskill, Predicting many properties of a quantum system from very few measurements, Nat. Phys. 16, 1050–1057 (2020).
  * Hu and You [2022] H.-Y. Hu and Y.-Z. You, Hamiltonian-driven shadow tomography of quantum states, Phys. Rev. Res. 4, 013054 (2022).
  * Liu  _et al._ [2023] Z. Liu, Z. Hao, and H.-Y. Hu, Predicting arbitrary state properties from single hamiltonian quench dynamics (2023), arXiv:2311.00695 .
  * Kunjummen  _et al._ [2023] J. Kunjummen, M. C. Tran, D. Carney, and J. M. Taylor, Shadow process tomography of quantum channels, Phys. Rev. A 107, 042403 (2023).
  * Levy  _et al._ [2024] R. Levy, D. Luo, and B. K. Clark, Classical shadows for quantum process tomography on near-term quantum computers, Phys. Rev. Res. 6, 013029 (2024).
  * Verstraete  _et al._ [2004] F. Verstraete, J. J. García-Ripoll, and J. I. Cirac, Matrix product density operators: Simulation of finite-temperature and dissipative systems, Phys. Rev. Lett. 93, 207204 (2004).
  * Verstraete  _et al._ [2006] F. Verstraete, M. M. Wolf, D. Perez-Garcia, and J. I. Cirac, Criticality, the area law, and the computational power of projected entangled pair states, Phys. Rev. Lett. 96, 220601 (2006).
  * Schuch  _et al._ [2007] N. Schuch, M. M. Wolf, F. Verstraete, and J. I. Cirac, Computational complexity of projected entangled pair states, Phys. Rev. Lett. 98, 140506 (2007).
  * Pérez-García  _et al._ [2008] D. Pérez-García, F. Verstraete, M. M. Wolf, and J. I. Cirac, Peps as unique ground states of local hamiltonians, Quantum Info. Comput. 8, 650–663 (2008).
  * Zaletel and Pollmann [2020] M. P. Zaletel and F. Pollmann, Isometric tensor network states in two dimensions, Phys. Rev. Lett. 124, 037201 (2020).
  * Soejima  _et al._ [2020] T. Soejima, K. Siva, N. Bultinck, S. Chatterjee, F. Pollmann, and M. P. Zaletel, Isometric tensor network representation of string-net liquids, Phys. Rev. B 101, 085117 (2020).
  * Kadow  _et al._ [2023] W. Kadow, F. Pollmann, and M. Knap, Isometric tensor network representations of two-dimensional thermal states, Phys. Rev. B 107, 205106 (2023).
  * Liu  _et al._ [2024] Y.-J. Liu, K. Shtengel, and F. Pollmann, Topological quantum phase transitions in 2d isometric tensor networks (2024), arXiv:2312.05079 .
  * Werner  _et al._ [2016] A. H. Werner, D. Jaschke, P. Silvi, M. Kliesch, T. Calarco, J. Eisert, and S. Montangero, Positive tensor network approach for simulating open quantum many-body systems, Phys. Rev. Lett. 116, 237201 (2016).
  * Cheng  _et al._ [2021] S. Cheng, C. Cao, C. Zhang, Y. Liu, S.-Y. Hou, P. Xu, and B. Zeng, Simulating noisy quantum circuits with matrix product density operators, Phys. Rev. Res. 3, 023005 (2021).
  * Kingma and Ba [2017] D. P. Kingma and J. Ba, Adam: A method for stochastic optimization (2017), arXiv:1412.6980 .
  * Crosswhite and Bacon [2008] G. M. Crosswhite and D. Bacon, Finite automata for caching in matrix product algorithms, Phys. Rev. A 78, 012356 (2008).
  * Verstraete and Cirac [2004] F. Verstraete and J. I. Cirac, Renormalization algorithms for quantum-many body systems in two and higher dimensions (2004), arXiv:cond-mat/0407066 .
  * Slattery and Clark [2021] L. Slattery and B. K. Clark, Quantum circuits for two-dimensional isometric tensor networks (2021), arXiv:2108.02792 .
  * Lin  _et al._ [2022] S.-H. Lin, M. P. Zaletel, and F. Pollmann, Efficient simulation of dynamics in two-dimensional quantum spin systems with isometric tensor networks, Phys. Rev. B 106, 245102 (2022).
  * Gu and Wen [2009] Z.-C. Gu and X.-G. Wen, Tensor-entanglement-filtering renormalization approach and symmetry-protected topological order, Phys. Rev. B 80, 155131 (2009).
  * Chen  _et al._ [2010] X. Chen, Z.-C. Gu, and X.-G. Wen, Local unitary transformation, long-range quantum entanglement, wave function renormalization, and topological order, Phys. Rev. B 82, 155138 (2010).
  * Pollmann  _et al._ [2010] F. Pollmann, A. M. Turner, E. Berg, and M. Oshikawa, Entanglement spectrum of a topological phase in one dimension, Phys. Rev. B 81, 064439 (2010).
  * Chen  _et al._ [2011] X. Chen, Z.-C. Gu, and X.-G. Wen, Complete classification of one-dimensional gapped quantum phases in interacting spin systems, Phys. Rev. B 84, 235128 (2011).
  * Pollmann and Turner [2012] F. Pollmann and A. M. Turner, Detection of symmetry-protected topological phases in one dimension, Phys. Rev. B 86, 125441 (2012).
  * Chen  _et al._ [2013] X. Chen, Z.-C. Gu, Z.-X. Liu, and X.-G. Wen, Symmetry protected topological orders and the group cohomology of their symmetry group, Phys. Rev. B 87, 155114 (2013).
  * Briegel and Raussendorf [2001] H. J. Briegel and R. Raussendorf, Persistent entanglement in arrays of interacting particles, Phys. Rev. Lett. 86, 910 (2001).
  * Raussendorf and Briegel [2001] R. Raussendorf and H. J. Briegel, A one-way quantum computer, Phys. Rev. Lett. 86, 5188 (2001).
  * Briegel  _et al._ [2009] H. J. Briegel, D. E. Browne, W. Dür, R. Raussendorf, and M. Van den Nest, Measurement-based quantum computation, Nat. Phys. 5, 19 (2009).
  * Raussendorf  _et al._ [2017] R. Raussendorf, D.-S. Wang, A. Prakash, T.-C. Wei, and D. T. Stephen, Symmetry-protected topological phases with uniform computational power in one dimension, Phys. Rev. A 96, 012302 (2017).
  * Raussendorf  _et al._ [2019] R. Raussendorf, C. Okay, D.-S. Wang, D. T. Stephen, and H. P. Nautrup, Computationally universal phase of quantum matter, Phys. Rev. Lett. 122, 090501 (2019).
  * Liu  _et al._ [2022] H. Liu, T. Zhou, and X. Chen, Measurement-induced entanglement transition in a two-dimensional shallow circuit, Phys. Rev. B 106, 144311 (2022).
  * Guo  _et al._ [2023] Y. Guo, J.-H. Zhang, Z. Bi, and S. Yang, Triggering boundary phase transitions through bulk measurements in two-dimensional cluster states, Phys. Rev. Res. 5, 043069 (2023).
  * Garnerone  _et al._ [2010] S. Garnerone, T. R. de Oliveira, and P. Zanardi, Typicality in random matrix product states, Phys. Rev. A 81, 032336 (2010).
  * Chen  _et al._ [2022] C.-T. Chen, Y.-H. Shi, Z. Xiang, Z.-A. Wang, T.-M. Li, H.-Y. Sun, T.-S. He, X. Song, S. Zhao, D. Zheng, K. Xu, and H. Fan, Scq cloud quantum computation for generating greenberger-horne-zeilinger states of up to 10 qubits, Sci. China-Phys. Mech. Astron. 65, 110362 (2022).
  * Guo and Yang [2023] Y. Guo and S. Yang, Noise effects on purity and quantum entanglement in terms of physical implementability, npj Quantum Inform. 9, 11 (2023).
  * Endo  _et al._ [2021] S. Endo, Z. Cai, S. C. Benjamin, and X. Yuan, Hybrid quantum-classical algorithms and quantum error mitigation, J. Phys. Soc. Jpn. 90, 032001 (2021).
  * Cai  _et al._ [2008] J.-F. Cai, E. J. Candes, and Z. Shen, A singular value thresholding algorithm for matrix completion (2008), arXiv:0810.3286 .

##  VIII Acknowledgements

This work is supported by the National Natural Science Foundation of China
(NSFC) (Grant No. 12174214 and No. 92065205), the National Key R&D Program of
China (Grant No. 2018YFA0306504), the Innovation Program for Quantum Science
and Technology (Grant No. 2021ZD0302100), and the Tsinghua University
Initiative Scientific Research Program. The IBM Quantum device ‘ibm_nairobi’
is accessible at https://quantum-computing.ibm.com/. The Quafu Quantum device
‘ScQ-P10’ is accessible at https://quafu.baqis.ac.cn/.

##  IX Author contributions

Y.G. conceived, designed, and performed the experiments. Y.G. and S.Y.
analyzed the data and wrote the paper. S.Y. contributed analysis tools and
supervised the project.

##  X Competing interests

The authors declare no competing interests.