# Enabling large-depth simulation of noisy quantum circuits with positive
tensor networks

Ambroise Müller  Eviden Quantum Lab, Les Clayes-sous-Bois, France  Thomas
Ayral  Eviden Quantum Lab, Les Clayes-sous-Bois, France  Corentin Bertrand
<EMAIL_ADDRESS>Eviden Quantum Lab, Les Clayes-sous-Bois, France

###### Abstract

MPDO are tensor network representations of locally purified density matrices
where each physical degree of freedom is associated to an environment degree
of freedom. MPDO have interesting properties for mixed state representations:
guaranteed positivity by construction, efficient conservation of the trace and
computation of local observables. However, they have been challenging to use
for noisy quantum circuit simulation, as the application of noise increases
the dimension of the environment Hilbert space, leading to an exponential
growth of bond dimensions. MPDO also lack a unique canonical form, due to the
freedom in the choice of basis for the environment Hilbert space, which leads
to a vast variation of bond dimensions.

In this work, we present a systematic way to reduce the bond dimensions of
matrix product density operators by disentangling the purified state. We
optimize the basis for the environment Hilbert space by performing density
matrix renormalization group (DMRG)-like sweeps of local 2-qubit basis
optimization. Interestingly, we find that targeting only the disentanglement
of the purified state leads to a reduction of the environment dimension. In
other words, a compact MPDO representation requires a low-entanglement
purified state.

We apply our compression method to the emulation of noisy random quantum
circuits. Our technique allows us to keep bounded bond dimensions, and thus
bounded memory, contrary to previous works on MPDOs, while keeping reasonable
truncation fidelities.

##  I Introduction

With the progress of experimental platforms for quantum computing, it is
becoming increasingly critical to produce accurate emulations of these
machines on classical hardware, taking into account their imperfections.
Indeed, finding algorithms relevant for noisy intermediate scale quantum
(NISQ) hardware and appropriate noise-mitigation techniques [1] requires an
understanding of the effect of noise. The controlled environment of a
classical emulator for large noisy quantum circuits is therefore a precious
tool.

Besides perfect (compressionless) emulators, whose memory and run time scale
exponentially in the number of qubits or gates (or both), approximate
emulators using the fruitful concept of tensor networks have emerged in the
last decade as methods of choice for this task [2, 3]. The simplest form of
TN, a one-dimensional network called the matrix product state (MPS), allows
for good approximations of weakly entangled pure states [4, 5], and for the
emulation of noiseless quantum circuits with hundreds of qubits [6, 7, 8, 2].
Other types of TNs better suited to entanglement patterns in higher than 1D
geometries have been used, such as tree tensor networks [9], multi-scale
entanglement renormalization ansatz [10], projected entangled pair states [5],
or more recently isometric tensor network states [11, 12].

In noisy quantum computing, the state is represented by a density operator
$\rho$, as opposed to a pure state $\ket{\Psi}$, which a priori leads to an
increased memory footprint. Yet, TN techniques can be applied to density
operators as well, with the intuition that more noise will lead to less
entanglement and therefore more representative power of compressed tensor
networks [13, 14]. A natural extension of MPS to density operators are matrix
product operators [15, 16, 17, 18, 13, 14, 19]: a product of matrices, each of
them acting locally on a single (or a few) qubits. Despite being powerful
compression tools for mixed states, MPOs suffer from some important flaws.
First, the positive-semidefinite character of the density matrix is not easily
enforced in MPO, and may be lost in particular during the truncation of bond
dimensions, leading to unphysical states. Second, taking expectation values
and normalizing the state require global operations, whereas MPS can be cast
in a canonical form that allows for local operations only. Finally, the
approximation made with an MPO is controlled by the operator entanglement
entropy (OEE) [20]—instead of the entanglement entropy for pure states— whose
relationship with entanglement is not straightforward.

There exists two alternatives to MPOs, both rooted in the idea of coming back
to pure states. The first is to unravel the density operator into a
statistical population of pure states of the same Hilbert space, leading to
powerful Monte Carlo methods [21]. The second, which is the focus of this
work, is to purify the density operator in a larger Hilbert space, extending
each original qubit with an “environment” or “ancilla” qudit of dimension $r$.
This purified state can be compressed into an MPS, with an internal bond
dimension $\chi$. This gives a locally purified TN representation with a
memory footprint proportional to $r\chi^{2}$, that we call in this work MPDO
111MPDOs are sometimes called locally purified tensor networks, locally
purified density operators or locally purified forms. Also, MPOs are sometimes
called MPDOs. See App. F. [15, 22, 23, 21, 19]. This form is guaranteed to be
positive semi-definite, and it is guaranteed to represent a physical state.
Moreover, in its canonical form, expectation values and normalization can be
obtained with local operations only. MPDO have been mostly used in quantum
many-body physics, e.g. for the unitary evolution of mixed states [24, 25],
imaginary-time evolutions [21, 24], or the evolution of open quantum systems
under a Lindblad equation [23, 21]. In Ref. [26], MPDOs are used as an
approximation for pure states, conserving local correlations only. In quantum
computing MPDOs have also been used to perform tomography [27].

Nevertheless, MPDOs have been challenging to use for the simulation of noisy
quantum circuits. Similarly as with MPOs, entangling gates cause an increase
of the bond dimension $\chi$, which needs to be reduced to avoid an
exponential growth in the number of entangling layers. A first reduction of
$\chi$ can be obtained with routine MPS compression techniques based on the
singular value decomposition (SVD) [28]. In addition, the freedom in the
choice of basis for the environment can be beneficial in reducing the
entanglement of the purified state. Yet, a global optimization is required to
take full benefit of this freedom, and can be done with successions of 2-sites
local optimizations [24, 29]. Other ideas were also studied, such as Uhlmann
parallel transport [25], or the application of backward time evolution to the
environment [30, 31, 32]. For MPOs, different techniques have also been
developed, such as density matrix truncation (DMT) [33].

On top of the increase in $\chi$, however, applying a noisy gate to an MPDO
increases the environment dimension $r$, leading to an exponential growth as a
function of the number of noisy layers, a serious drawback compared to MPOs.
Note that this issue is absent from real or imaginary-time evolution of closed
systems. A straightforward idea to try to reduce $r$ is to perform truncated
SVD on the environment bonds, which boils down to an optimization of single-
site environment bases followed by a truncation of unused dimensions [23, 21,
34]. Nevertheless, this local compression (LC) technique, due to its single-
site nature, does not avoid an exponential increase of $r$, as can be seen
from Ref. [34]. Extending it to a global optimization of the environment basis
seems a natural way to go, but raises some questions. While the optimization
problem to reduce $\chi$ has a clear cost function, i.e. the entanglement
entropy of the purified state, here we face the coupled problem of reducing
both $\chi$ and $r$. From an optimization point of view, there is no obvious
choice of cost function. Note that ideas that do not involve the environment
basis freedom have also been explored, such as modifying the environment sites
distribution along the MPDO, e.g. to gather all the environment dimensions on
a single site [35].

In this work, we show that the global optimization of the environment basis
targeting minimal entanglement of the purification leads to drastic reduction
in $r$, and thus forms a good compression scheme for MPDOs. Such an
optimization is expected to reduce $\chi$ [24, 29], but its impact on $r$ was
never thoroughly studied. We perform a global change of the environment basis
with DMRG-like sweeps of local 2-qubit optimizations, each of which
disentangles the two qubits. We will refer to this procedure as the iterative
purification disentanglement (IPD) technique. We find that aiming at
disentangling the purified state is enough to reduce both $r$ and $\chi$,
giving a simple cost function for the compression procedure. In fact, we
observe that we can systematically reduce $r$ to a value close to the optimal
value of $2$, with a very high truncation fidelity. Thanks to this
compression, we are able to emulate a noisy quantum random circuit with MPDOs
at arbitrary depth with bounded memory, which was impossible with the LC
technique.

This article is structured as follows. We start by a pedagogical introduction
to MPDOs in Sec. II.1, followed to a description of the compression method in
Sec. II.2. Results are presented in Sec. III. First, the compression of a
simple two-qudit state is studied (III.1), before moving to a many-qubit state
(III.2). Application to a full circuit simulation is shown in Sec. III.3.
Discussion and conclusion are regrouped in Sec. IV.

##  II Method

In this section we first give an overview of the representation of mixed
quantum states with MPDOs. We then detail the IPD technique used to obtain
compressed MPDO representations.

###  II.1 MPDO states: general considerations

Figure 1:  (a) Tensor network representation of a MPDO state. (b) Extraction
of singular values on entanglement bond at (bond) index $i=2$. (c) Extraction
of singular values on purification bond at index $i=2$. (d) Normalization
identities for singular value distributions.  Figure 2:  Tensor network
representation of the application to a MPDO state of (a) a one-site noisy
channel with Kraus operators $\\{K_{i}\\}_{i=1}^{k}$, and (b) a two-site noisy
channel. For the two-site case a QR decomposition is first performed on the
two-qubit Kraus operators (allowing for an even splitting of the Kraus index)
and the resulting tensors are then contracted with the site tensors. The
purification bond dimensions are multiplied by the local Kraus ranks $k$,
$k_{1}$, $k_{2}$, and the entanglement bond dimension is multiplied by the QR
rank $q$.

Given a Hilbert space $\mathcal{H}$ and a mixed state $\rho$ on $\mathcal{H}$,
a purification of $\rho$ is a pure state $\ket{\Psi}$ defined on an extended
Hilbert space $\mathcal{H}\otimes\mathcal{H}^{\prime}$, where
$\mathcal{H}^{\prime}$ is an ancillary system, and such that

$\rho=\text{Tr}_{\mathcal{H}^{\prime}}\ket{\Psi}\bra{\Psi}.$ (1)

$\text{Tr}_{\mathcal{H}^{\prime}}$ denotes the partial trace on the ancillary
system. Purifications always exist and are not unique in general. The
ancillary Hilbert space $\mathcal{H}^{\prime}$ need not have the same
dimension as $\mathcal{H}$, though its dimension must be at least equal to the
rank of $\rho$. In particular there always exists a purification of rank
$2^{N}$ for an arbitrary state of $N$ qubits. As any pure state, a
purification may be written in MPS form:

$\ket{\Psi}=\sum_{\\{s_{i}\\}}\sum_{\\{a_{i}\\}}\Tr{\prod_{i=1}^{N}A_{i}^{s_{i}a_{i}}}\ket{s_{1}a_{1},\dots\
,s_{N}a_{N}},$ (2)

where $s_{i}\in\\{1,...,d_{i}\\}$ are the indices of physical degrees of
freedom ($d_{i}\equiv 2$ for a system of qubits) and
$a_{i}\in\\{1,...,r_{i}\\}$ correspond to the ancillary degrees of freedom.
For each $i$, $A_{i}^{s_{i}a_{i}}$ is a matrix of dimension
$\chi_{i}\times\chi_{i+1}$. We will refer to $\chi_{i}$ as _entanglement bond
dimensions_ and to the corresponding bonds as _entanglement bonds_. By
contrast, we refer to $r_{i}$ as _purification bond dimensions_ , and to the
corresponding bonds as _purification bonds_.

We can combine Eqs. (1) and (2) into the TN diagrammatic expression of Fig.
1a. The physical Hilbert space has dimension $d=\prod_{i=1}^{N}d_{i}$ and the
ancillary system has dimension $r=\prod_{i=1}^{N}r_{i}$. The physical state is
obtained by tracing out the purified state over the ancillary degrees of
freedom, resulting in internal bonds. Note that since the tensors in the
bottom row are adjoints to those in the top row, one only needs to store one
of them (i.e. the purified state vector or its dual) in memory.

The singular value decomposition (SVD) is an essential tool in TN-based
methods. It allows to associate a singular value distribution to any internal
bond. SVD can be used to fix gauge freedom intrinsic to (non-looped) tensor
networks and introduce canonical forms [36], yielding left-unitary and right-
unitary tensors (which we denote by $L$ and $R$ respectively) and ‘center’ the
representation on a particular bond, of which singular values are extracted.
Figures 1b and 1c show the extraction of singular values across an
entanglement bond and a purification bond, respectively, in the appropriately-
centered canonical forms. The singular value distributions associated with
entanglement bonds are represented by the white diamonds, and will be denoted
$\quantity{\sigma_{j}^{(i)}}_{j=1}^{\chi_{i}}$ for bond index $i$. On the
other hand, the singular value distributions associated with purification
bonds are shown as a blue diamond, and will be denoted
$\quantity{\lambda_{j}^{(i)}}_{j=1}^{r_{i}}$ for bond index $i$.

These singular values satisfy the normalization conditions shown in
diagrammatic notation in Fig. 1d and which reads

$\Tr\rho=\sum_{j=1}^{r_{i}}\lambda_{j}^{(i)}=\sum_{j=1}^{\chi_{i}}\quantity(\sigma_{j}^{(i)})^{2}$
(3)

for all $i=1,\ldots,N$. For a normalized state ($\Tr\rho=1$) this enables us
to associate a probability distribution to each MPDO bond, to which we can
associate an entropy measure. In this vein, if $X$ is a positive semidefinite
matrix whose singular values $\quantity{x_{j}}_{j=1}^{n}$ form a probability
distribution, we define the Rényi entropies [37]

$E_{\alpha}\left[X\right]=\frac{1}{1-\alpha}\log\left(\sum_{j=1}^{n}x_{j}^{\alpha}\right),$
(4)

where $\alpha\geq 0$ is the entropy order. The (common) choices $\alpha=1$ and
$\alpha=2$ simplify to

$\displaystyle E_{1}[X]$ $\displaystyle=-\Tr{X\log X},$ (5) $\displaystyle
E_{2}[X]$ $\displaystyle=-\log\Tr{X^{2}}.$ (6)

We will refer to the entropy of the distribution associated to a bond as the
_entanglement bond entropy_ or _purification bond entropy_ , denoted
respectively by

$\displaystyle\mathcal{E}_{\alpha}^{i}=E_{\alpha}\left[(\Sigma^{(i)})^{2}\right]\
\text{and}\ \mathcal{P}_{\alpha}^{i}=E_{\alpha}\left[\Lambda^{(i)}\right],$
(7)

where we defined the diagonal singular value matrices
$\Sigma^{(i)}=\text{diag}\left(\quantity{\sigma_{j}^{(i)}}_{j=1}^{\chi_{i}}\right)$
and
$\Lambda^{(i)}=\text{diag}\left(\quantity{\lambda_{j}^{(i)}}_{j=1}^{r_{i}}\right)$.
Note that $\mathcal{E}_{\alpha}^{i}$ is the $\alpha$-entanglement entropy of
the purified state between left and right of index $i$ (both sides including
system and ancilla degrees of freedom). $\mathcal{E}_{2}^{i}$ is therefore
related to the purity of the corresponding reduced density matrix. Moreover,
$\mathcal{P}_{\alpha}^{i}$ is the $\alpha$-entanglement entropy between the
ancilla degree of freedom at index $i$ and the all other degrees of freedom.

The simulation of noisy quantum circuits with MPDOs is rather straightforward
once the circuit has been decomposed into one-qubit and nearest-neighbour two-
qubit channels. Distant two-qubit channels can be made nearest-neighbour by
introducing a series of SWAP gates [38, 39], even though this incurs a large
increase in the number of entangling gates. A generic quantum channel is a
completely positive trace-preserving map [40], and admits the decomposition

$\rho\longrightarrow\sum_{\ell=1}^{k}K_{\ell}\rho K^{\dagger}_{\ell},$ (8)

where $\\{K_{\ell}\\}_{\ell=1}^{k}$ are known as Kraus operators. The sum over
the Kraus index (of dimension $k$) can be construed as an additional bond in
TN diagrammatic notation, as illustrated in Fig. 2. Applying the channel
consists in contracting the Kraus operators with the qubit site tensors and
grouping the Kraus index with the purification bond. This has the result of
increasing the purification bond dimension by a factor $k$. In an exact
circuit simulation, the purification bond dimensions thus grow exponentially
with circuit depth. Two-qubit gates act in a similar way, though also
increasing the dimensions of entanglement bonds, much in the same way as with
MPS or MPO simulations.

###  II.2 IPD

Our proposed method relies on the fact that purifications of a mixed quantum
state are not unique. For any mixed state one is free to construct a
purification requiring an ancilla system much larger than the physical system,
yielding large purification bond dimensions (in MPDO form), or which are
highly entangled even for a separable physical state, yielding large
entanglement bond dimensions. In order to reach a memory-efficient
representation, it is thus necessary to find a purification of small rank with
minimal entanglement, i.e. with minimal quantum correlations between ancillary
degrees of freedom but still retaining the quantum correlations in the
physical system. Given two purifications
$\ket{\Psi}\in\mathcal{H}\otimes\mathcal{H}^{\prime}$ and
$\ket{\Psi}^{\prime}\in\mathcal{H}\otimes\mathcal{H}^{\prime\prime}$ of a
quantum state $\rho$ there exists a partial isometry
$U:\mathcal{H}^{\prime}\rightarrow\mathcal{H}^{\prime\prime}$ (i.e.
$U^{\dagger}U$ is a projector onto $\mathcal{H}^{\prime}$, and, equivalently,
$UU^{\dagger}$ is a projector onto $\mathcal{H}^{\prime\prime}$) such that

$\ket{\Psi}^{\prime}=(\mathds{1}_{\mathcal{H}}\otimes U)\ket{\Psi}.$ (9)

In particular, $U$ is unitary when
$\dim\mathcal{H}^{\prime}=\dim\mathcal{H}^{\prime\prime}$[40] and can be
interpreted as a change of basis for the ancilla system. This translates to a
gauge freedom in the MPDO representation, whereby a unitary transformation may
be applied across any number of purification bonds without altering the
physical state.

Applying this gauge freedom with single-site operators $U$ allows to change
the local basis of the ancilla system. An optimal local basis to reduce the
purification bond dimensions can be obtained with an SVD, truncating the
lowest singular values, as in Refs. [23, 34]. We will refer to this as the
local compression (LC) technique.

Now, consider two adjacent sites, at indices $i$ and $i+1$. An optimal basis
for these two sites cannot be found simply with an SVD, as the size of the
entanglement bond between the two sites must be taken into account. Instead,
we first contract the sites into a tensor $A$ and group the purification bonds
into a single bond of dimension $r=r_{i}r_{i+1}$. We apply an optimized
$r\times r$ unitary $U$ on the ancilla index. As illustrated in Fig. 3, this
unitary can be understood as a bond partitioning operation, whereby the
contracted local purification bond of dimension $r$ is factorized again into
two bonds of dimension $r_{i}$ and $r_{i+1}$. Denoting by $\tilde{A}$ the
gauge-transformed tensor, where left and right indices are grouped to form a
matrix of shape $d_{i}\chi_{i}r_{i}\times d_{i+1}\chi_{i+1}r_{i+1}$, one can
extract the entanglement bond entropy
$\mathcal{E}_{\alpha}^{i}=E_{\alpha}\quantity[\tilde{A}^{\dagger}\tilde{A}]$.

We choose to optimize $U$ so as to minimize this entropy. This means that we
need to solve the problem

$\displaystyle\underset{U}{\text{minimize}}\ \mathcal{E}_{\alpha}^{i}$ (10)
$\displaystyle\text{such that}\ UU^{\dagger}=U^{\dagger}U=\mathds{1}.$

Perhaps surprisingly, we choose a cost function that depends only on the
entanglement bond, and no the purification bond. This is motivated by the
intuition that a purification with a simpler entanglement structure might be
able to be described by fewer degrees of freedom. We numerically check that
both bond dimensions will be reduced with this cost function in Sec. III.

Any minimum of Eq. (LABEL:eq:entropy-opt) is necessarily non-unique, since we
can identify equivalence classes of gauge unitaries. Indeed, two unitaries
which are related through the tensor product of unitaries of shape
$r_{i}\times r_{i}$ and $r_{i+1}\times r_{i+1}$ (acting on the separated
purification bonds), i.e. where only the single-site gauge choice differs,
will yield the same entanglement bond entropy.

The optimization problem in Eq. (LABEL:eq:entropy-opt), with either entropy
measure, can be solved via gradient descent. In our results we used the
implementation of the conjugate gradient (CG) algorithm provided by SciPy.
Details of the gradient computation are shown in App. A. An a priori advantage
of $\mathcal{E}_{2}^{i}$ over $\mathcal{E}_{1}^{i}$, as already observed in
Ref. [24], is the absence of a matrix logarithm, thus simplifying the function
and gradient evaluations.

Figure 3:  Upon contracting two adjacent sites of the purification MPS into a
tensor $A$, one is free to apply a gauge unitary $U$ on the ancillary indices
without affecting the physical state. The choice of unitary will affect the
distribution of singular values $\Sigma^{(i)}$ across the entanglement bond
upon re-separation of the sites via SVD of the transformed tensor
$\tilde{A}=L\Sigma^{(i)}R$.

To generalize to $N$ qubits, we do a sequence of sweeps along the MPDO, in
alternating directions. A sweep is performed in a two-site DMRG fashion, that
is, over pairs of adjacent sites at indices $i$ and $i+1$, for $i=1$ to
$i=N-1$. On each pair of sites, we solve problem (LABEL:eq:entropy-opt) and
apply the corresponding disentangling operator $U$. Upon separating the two
sites via SVD, we usually see an increase in bond dimension, due to the
introduction of small singular values. To get the benefits of the change of
representation, we truncate singular values that fall under a given threshold
$\varepsilon$. After every sweep, we also truncate singular values on local
purification bonds below this threshold. A small $\varepsilon$ guarantees a
small truncation fidelity after a set number of sweeps, but a low compression
rate. A high $\varepsilon$ gives rise to a better compression rate, but at the
risk of a reduced fidelity. Truncation error is guaranteed to be bounded by
the total weight of removed singular values, for both the entanglement and
purification bond, thanks to the canonical form.

Figure 4:  Purification optimization on a random state composed of two qudits
of dimension $d=4$. For each gradient descent step we show (a) the $\alpha=1$
and $\alpha=2$ entanglement bond entropies, (b) the $\alpha=1$ purification
bond entropy on the left site and (c) the $\alpha=1$ purification bond entropy
on the right site. We plot singular value distributions before and after
optimization on (d) the entanglement bond, (e) the left purification bond and
(f) the right purification bond. Optimizing only the entanglement bond
entropy, we are able to reduce the number of significant singular values on
purification bonds down to the physical dimension $d=4$.

We highlight the fact that similar techniques for optimizing a purified state
entanglement have been introduced for finding the purification entanglement
entropy [29] or to compress a MPDO during unitary or imaginary-time evolution
[24]. However, in these works the purification bond dimensions was not
changing and was fixed to $2$, eliminating the fundamental issue we are
contributing to solve here. Ref. [24] also employed a different optimization
technique to find the disentangler operator, solving a fixed point equation,
while we perform a gradient descent. We believe that gradient descent is a
better option, since fixed point iterations are not guaranteed to converge.

The problem of finding a bond factorization that minimizes entanglement was
also identified by [23] but applied only to the decomposition of two-qubit
gates. Our method for disentangling two qubits is fairly similar to the
technique used in Ref. [23] to disentangle two-qubit noisy gates.

##  III Results

In this section, we make a detailed numerical analysis of each step of our
compression method. First we focus on the local optimization in a two-qudit
case, before exploring the optimization of a multi-qubit MPDO with sweeps.
Finally, we integrate the technique within a noisy circuit emulator.

###  III.1 One-shot compression of a two-qudit MPDO

In Fig. 4 we showcase the compression of a two-qudit MPDO with the IPD
technique. We initialize a random MPDO of two qudits of dimension $d=4$, with
(suboptimal) local purification dimensions $r=8$ and entanglement bond
dimension $\chi=32$. The constituting tensors were generated with random
entries uniformly distributed over the region $[-1,1)\times[-1,1)$ of the
complex plane.

In panels (a)-(c) we record the variation of the bond entropies against the
number of gradient descent steps. We then optimize the gauge unitary so as to
minimize the 2-Rényi entropy across the two sites. As expected, we observe in
panel (a) that this procedure decreases the entropy across the entanglement
bond; not only the $\alpha=2$ Rényi entropy, minimized by construction, but
also the $\alpha=1$ entropy. As a result, the distribution of singular values
becomes much steeper (panel (d)).

Interestingly, this procedure has also the effect of reducing the purification
bond entropies (panels (b) and (c)), thereby simplifying the correlations
between the physical and ancilla systems. This translates to a sharp drop in
the corresponding (sorted) singular values at index $4$ (panels (e) and (f)),
suggesting that the optimization has steered the MPDO state to a purification
of optimal rank. Indeed, we are now able to truncate the low singular values
of index larger than 4 with very high fidelity, thereby reaching a total
purification dimension of $4^{2}$, which is also the dimension of the Hilbert
space. A similar behaviour was observed for other starting values of $d$ and
$r$ (not shown). This demonstrates that it is possible to reduce local
purification dimensions through an optimization procedure focused only on the
entanglement bonds.

###  III.2 One-shot compression of a generic MPDO

(a)

(b)

(c)

Figure 5:  Local purification disentanglement (16 sweeps, truncation threshold
$\varepsilon=10^{-3}$) of a $N=6$ qubits random MPDO state with initial
entanglement bond dimensions $\chi=12$ and local purification dimensions
$r=8$. We show in (a) the variation of total fidelity and memory usage, and
the initial and final singular value distributions on (b) entanglement bonds
and (c) local purification bond.

Fig. 5 shows the performance of the IPD technique on a random MPDO state of
$N=6$ qubits. Details on the construction of the random state are given in
App. D. We used a noise level $p=10\%$ over ten circuit layers, truncating
entanglement bonds to a dimension $\chi=12$ and purification bonds to
dimension $r=8$. We set the truncation threshold to $\varepsilon=10^{-3}$,
close to the smallest singular value on any bond in the starting MPDO
representation. Fidelities are computed by contracting the tensor networks and
using the Uhlmann quantum state fidelity [41] defined by

$\mathcal{F}(\rho,\sigma)=\left[\Tr{\sqrt{\sqrt{\rho}\sigma\sqrt{\rho}}}\right]^{2}.$
(11)

For the $s$-th sweep, we plot the total fidelity
$\mathcal{F}(\rho_{s},\rho_{0})$, the memory used to store $\rho_{s}$ in MPDO
form, where $\rho_{s}$ is the state after the $i$-th disentangling sweep, as
well as the average entanglement and purification bond dimensions in the MPDO
representation (panel (a)). We obtain a gain of an order of magnitude in
purification bond dimension and overall memory, while conserving a very good
fidelity $>0.99$.

We observe that, as in the two qudit case, the optimization of local
entanglement entropies has a substantial effect on the singular value
distributions across local purification bonds. Indeed, successive sweeps
gradually decrease the singular values starting from the third index. After 16
sweeps these have decreased by several orders of magnitude and can be
truncated with high fidelity (panel (c)). We also observe that singular value
distributions on entanglement bonds have become steeper as a result of the
local entanglement entropy minimization (panel (b)). More detailed data
including singular value distributions in the intermediate states can be found
in App. B, where we also discuss the convergence of this scheme.

Figure 6:  Compression results of a truncated SVD scheme applied on the state
used in Fig. 5. The truncation threshold corresponds to the ratio of singular
values discarded on each bond. The compression ratio is the ratio of memory
usage after compression to that before compression.

One may ask if this method yields a more precise compression than the LC
technique, as used in e.g. [33] and [34]. Indeed our proposed scheme performs
a significantly higher number of truncations due to its iterative nature. In
Fig. 6 we show the compression fidelity and compression ratio obtained by
naive compression of the same quantum state as in Fig. 5 via truncated SVD for
different truncation thresholds. We see that our method is able, for a certain
memory compression ratio, to find a better approximation of the initial
quantum state; alternatively, for an infidelity budget of $\sim 10^{-2}$, it
finds a MPDO representation an order of magnitude smaller in memory than that
found by the LC technique.

This improvement in the one-shot scenario comes with a drawback, which is that
the memory usage has to be temporarily increased to accommodate growing
entanglement bond dimensions, until further singular values on entanglement
bonds fall below the truncation threshold. In Fig. 5 one can see that the
memory usage increases temporarily by approximately an order of magnitude,
before decaying to a value approximately an order of magnitude smaller than
its initial size. This new compression procedure also experiences an increase
in computation time, primarily due to local optimization routines, which
become the limiting factor in terms of runtime.

The compression was performed by minimizing the $\alpha=2$ Rényi entropy, as
we found that it performed better than minimizing $\alpha=1$. Not only the
gradient is easier to compute (see App. A for the expressions of each
gradient), but the minimization itself converges faster, as described in App.
C.

###  III.3 Application to simulating noisy quantum circuits

Figure 7:  Comparison of noisy random circuit ($N=8$ qubits, noise level
$p=10\%$) simulations with MPO, LC-MPDO and IPD-MPDO emulators. After each
circuit layer we show (a) the pseudo-fidelity (defined in Eq. (13)) with the
exact simulation result, (b) the operator entanglement entropy (OEE), (c) the
entanglement entropy of the purified state, (d) the average entanglement bond
dimension, (e) the average purification bond dimension and (f) the memory
usage after the compression routine. Values are averaged over 8 random
circuits.  Figure 8:  Same as Fig. 7, with noise level $p=2\%$. Results for a
single random circuit.  Figure 9:  Simulation of a noisy random circuit ($N=8$
qubits, noise level $p=5\%$) with fixed bond maximum dimensions and different
number of purification disentangling sweeps. After each circuit layer we show
(a) the pseudo-fidelity (Eq. (13)) with the exact simulation result, (b) the
operator entanglement entropy (OEE) and (c) the entanglement entropy of the
purified state.

We apply the purification disentanglement routine to MPDO-based noisy circuit
simulations to limit the growth in memory usage. We chose to simulate random
circuits sampled as detailed in App. D, with $N=8$ qubits and depolarizing
noise

$\rho\longrightarrow(1-p)\rho+p\frac{\mathds{1}}{2}$ (12)

of various strengths $p$. Such a depolarizing noise leads at high depths to a
maximally mixed state (MMS), i.e. $\rho_{\rm MMS}=\mathds{1}/2^{N}$. Since we
want to be able to compress arbitrary states, and not only the trivial MMS, we
quantified how far we are from it in App. E.

We also compare the results of our simulations with implementations of
emulators employing the LC technique, in which memory growth is prevented only
by performing truncated SVDs on entanglement and purification bonds.

In this section we substitute the usual quantum state fidelity with the
pseudo-fidelity defined in [42] by

$\mathcal{F}_{\text{alt}}(\rho,\sigma)=\frac{\left|\Tr{\rho^{\dagger}\sigma}\right|}{\sqrt{\Tr{\rho^{\dagger}\rho}\Tr{\sigma^{\dagger}\sigma}}}$
(13)

for two density operators $\rho$ and $\sigma$. This is done to accommodate the
comparison with MPO simulators, for which Eq. (11) may be ill-defined, due to
the fact that positive-semidefiniteness is not enforced.

####  III.3.1 Simulations with adaptive bond dimensions

Figure 7 shows data for the simulation of a $N=8$ qubit random circuit with
noise level $p=10\%$. We compare the IPD technique with MPO and MPDO
simulators using the LC technique, where truncation is ‘adaptive’, i.e.
singular values are truncated up to a certain threshold (tolerance).
Truncation tolerances for each simulator were tuned so as to obtain comparable
memory usage growth at short depths. Values plotted were obtained with
$\varepsilon_{MPO}=0.01$, $\varepsilon_{SC}=0.015$ and
$\varepsilon_{DGC}=0.001$. For the IPD technique 8 sweeps were performed after
every circuit layer. A safety maximum bond dimension of 64 (regardless of the
bond type) was set for all three simulators. The plain blue line corresponds
to our method, the plain black line to the LC technique, and the dotted black
line to an MPO simulation. The data is averaged over $8$ random circuits, with
error bars estimated from the standard deviation.

First notice that the IPD technique performs similarly to MPO in terms of
fidelity (panel (a)) and memory (panel (f)), while improving largely compared
to the LC technique. On panel (a), one can see that the maximal infidelity is
half that of the LC technique, and is further reduced at large depths $>13$,
while being similar to MPO simulation.

Importantly, in our case the memory footprint is bounded, while it grows
exponentially with the circuit depth with the LC technique. This is explained
by the evolution of both the entanglement and purification bond dimensions,
which can be seen averaged over sites respectively in panels (d) and (e). In
the absence of disentanglement, they both increase exponentially with depth,
due to the application of entangling gates (for entanglement bonds) and noisy
gates (for purification bonds). With disentanglement, in contrast, the
entanglement bond dimensions reach a maximum before going down to one,
similarly to the bond dimension of MPO; and, strikingly, the purification bond
dimensions are reduced down to two, consistently along the circuit at this
level of noise. This value of two is optimal in the sense that one cannot hope
to do better for a generic density matrix, as its rank is bounded by the
Hilbert space dimension $2^{N}$.

Panel (b) shows the OEE for the three methods. The OEE is obtained by
contracting the MPDO states along purification bonds, yielding MPO
representations, and computing the normalized bond entropy on the middle bond
in mixed canonical form. Since OEE is a physical quantity, discrepancies are
explained by the approximate nature of the simulations, and not by the choice
of purification. Once again, we observe a very good agreement of our method
with the MPO simulation, better than the LC technique.

Finally, panel (c) shows that the entanglement entropy of the purified state,
which coincides with the $\alpha=1$ entanglement bond entropy computed on the
middle entanglement bond. We observe that while the entropy obtained with the
LC technique grows in a linear fashion with depth, independently of the OEE
(understood as a measure of physical entanglement), the IPD technique allows
the simulator to restrain the growth in the entanglement entropy of the
purification. This makes explicit the fact that disentangling sweeps are able
to find an efficient representation of the state at each depths by finding a
purification with low ‘virtual’ entanglement (i.e. low entanglement between
ancilla qubits).

Overall, the IPD technique produces a representation of the state of much
better quality than the LC technique, and of similar quality as in a MPO
simulation, with all the advantages of the purification form. We see in Fig. 8
that this conclusion also holds for a weaker noise level $p=2\%$. In fact, we
argue that we improve a bit compared to MPO, since we produce, after depth
$25$, a representation of significantly smaller memory than MPO (see panel
(f)), caused by a smaller entanglement bond dimension (panel (d)).

In Fig. 8 truncation tolerances for the three emulators were once again tuned
so as to produce comparable memory growth at short depths, yielding
$\varepsilon_{MPO}=0.015$, $\varepsilon_{SC}=0.005$ and
$\varepsilon_{DGC}=0.0005$. More significantly, the number of sweeps performed
after every circuit layer was increased to 32.

Another difference with the stronger noise case is that the purification bond
dimensions (panel (e)) show a maximum ($<5$) before being reduced to the
optimal value of two. A lower noise level indeed leads to less entanglement
destruction, and thus an entanglement pattern that is more complex to
disentangle for our algorithm.

####  III.3.2 Performance at fixed bond dimensions

Thankfully, in our approach one can increase the number of sweeps to improve
the quality of the optimization, and thus of the purification. This is
illustrated in Fig. 9, which compares the pseudo-fidelity, the OEE and the
entanglement entropy of the purified state (the latter two computed on the
middle bond) with an increasing number of disentangling sweeps per circuit
layer. Here, to make a fair comparison, we avoid the somewhat arbitrary nature
of the truncation tolerance, by fixing all bond dimensions of the MPDO and
focus on entropies rather than bond dimensions. Entanglement bonds are
truncated to a maximum dimension of $\chi=32$, while purification bonds are
truncated to the rank-optimal dimension $r=d=2$. This truncation is performed
after every circuit layer (and after the disentangling sweeps in the case of
the IPD technique). The circuit simulated is a single random quantum circuit
with $8$ qubits as used in the results of Fig. 7. We used a depolarizing noise
level of $p=5\%$. We clearly see an improvement in fidelity (panel (a)) and a
systematic decrease in OEE (panel (b)) and purification entropy (panel (c))
with more sweeps. An decrease in entropy translates into a reduction in
memory, as it directly shows that less singular values need to be kept. Note
how the OEE gets closer to the MPO simulation values (panel (b)), which seems
to act as lower bound.

These results indicate that, with a correct disentangling technique such as
the one we propose, the representative power of MPDO is as good as MPO all
along the circuit. In their recent work [43], Guo and Yang argue on the
contrary that MPDOs have a good representative power only at low depths and at
high depths, and lack representative power in between when compared to MPOs.
They also observe that the separation between these two regimes is close to
the maximum of OEE. If we indeed observe a minimum of fidelity and a maximum
of memory with our MPDOs at intermediate depths, the values are still fairly
close to the MPO simulation, and thus reflect an increase of difficulty that
is shared with MPO techniques. In addition, we see no co-location of this
difficulty barrier with the maximum of OEE. For example in Fig. 7, the maximum
of OEE lies around depth $3$, while the minimum of fidelity moves between $10$
and $14$ depending on the number of sweeps, and the maximum of purification
entropy is between depths $6$ and $10$.

This indicates that the barrier observed in the data of Ref. [43] is caused by
the MPDO compression they perform and by the MPO-to-MPDO mapping they use. The
IPD technique, although more involved computationaly, leads to MPDOs of better
quality (higher fidelity, lower memory), and are thus able to overcome the
complexity barrier of Ref. [43]. One can hypothesize that we explore a larger
subspace of equivalent purifications, which allows for a more compact
representation of the state.

##  IV Discussion and conclusion

In conclusion, we presented a compression technique for MPDOs exploiting the
freedom of choosing the basis for the environment qudits in a purification.
Instead of optimizing this choice locally on each qubit, we perform it
globally, by applying optimized changes of basis on successive pairs of
qudits, organized in sweeps. This amounts to disentangling the purified state
in view of finding a representation with low entanglement and purification
bond entropies. We then truncate the smallest singular values to finish the
compression.

We find that a good criterion for this local optimization is to minimize the
entanglement bond entropy. With this choice of cost function, we observe,
surprisingly, a reduction of _both_ the entanglement and purification entropy.
This is a surprisingly simple cost function for the minimization of two
seemingly distinct entanglements sources. We recall that entanglement bond
entropy is produced by entangling gates, while purification bond entropy is
produced by noisy gates. If a relationship exists between these, it is not
obvious with MPDOs. For instance, noisy gates, when they are a source of
decoherence, usually destroy entanglement[13], but this does not automatically
translate, in MPDOs, to the reduction of entanglement bond entropies.

In a way, one can see our compression scheme as a way to recover this
relationship. Another way to interpret this choice of cost function is that
the IPD technique is similar in spirit as trying to find the entanglement of
purification [44], which is in our language the minimal entropy of a given
entanglement bond, over all equivalent purifications. This quantity captures
the overall complexity of a bipartite mixed state, including both quantum and
classical correlations [44]. As classical correlations come from noisy gates,
it seems natural in this context that reducing entanglement bond entropy leads
also to a reduced purification bond entropy.

In fact, for every state we explored, we were able to reduce the purification
bond dimensions to a value close to $2$, with a very good truncation fidelity.
This value of $2$ is optimal for a generic state, in the sense that it
corresponds to a decomposition $\rho=X^{\dagger}X$ with $X$ of dimension
$2^{N}\times 2^{N}$, which is minimal given a generic $\rho$ with rank
$2^{N}$. It is left to future work to establish if our technique can reduce
$r$ below this value when the density matrix is of lower rank, such as when
noise is low.

We also compared two measures of entanglement in App. C and found that locally
minimizing the 2-Rényi entanglement bond entropy yields a better final
purification with faster convergence.

This work also brings a new perspective on the question of the
representability of practical states with MPDOs. The existence of states that
are easily represented as MPOs but difficult as MPDOs has been proved in Ref.
[22]. However, this does not close the debate for practical states, i.e.
states encountered in typical quantum computing programs. Our results seem to
indicate that MPDOs have a similar representation power compared to MPOs, at
least for noisy random circuits. These contradict the conclusions of Ref.
[43], which observed a sharp loss of representability close to the maximum of
OEE, but used simpler compression methods. Using the IPD technique, it seems
that a good MPDO representation of practical states can always be found. It is
possible however, that such a good representation is harder to reach than for
MPOs. This would not be surprising given the advantages of MPDOs over MPOs:
enforcing positive-semidefiniteness, and providing efficient computation of
observables, trace and reduced density matrices thanks to its canonical form.
Further work is needed to sharpen this point.

Despite being efficient at compressing a MPDO, the IPD technique in its
present version is not computationally cheap. Further studies are required to
better understand how the compression ratio scales with the number of qubits,
and how it is impacted by 2D geometry, in order to evaluate if it is a good
option for the simulation of circuits of practical sizes and connectivities.
In addition, the tuning of parameters (number of sweeps, truncation
thresholds, etc) could be optimized depending on noise level. Performances
would also benefit from improvements of the optimization under unitary
constraints. For all these reasons, the precise assessment of computation
times is left for further works.

###### Acknowledgements.

This work is part of HQI initiative (www.hqi.fr) and is supported by France
2030 under the French National Research Agency award number
“ANR-22-PNCQ-0002”. Computations have been done on Eviden Qaptiva.

## References

  * Cai  _et al._ [2022] Z. Cai, R. Babbush, S. C. Benjamin, S. Endo, W. J. Huggins, Y. Li, R. Jarrod, and T. E. O. Brien, (2022), arXiv:2210.00921 .
  * Ayral  _et al._ [2023] T. Ayral, T. Louvet, Y. Zhou, C. Lambert, E. M. Stoudenmire, and X. Waintal, PRX Quantum 4, 020304 (2023), arXiv:2207.05612 [cond-mat, physics:quant-ph] .
  * Tindall  _et al._ [2023] J. Tindall, M. Fishman, M. Stoudenmire, and D. Sels, (2023), arXiv:2306.14887 [quant-ph] .
  * Paeckel  _et al._ [2019] S. Paeckel, T. Köhler, A. Swoboda, S. R. Manmana, U. Schollwöck, and C. Hubig, Annals of Physics 411, 167998 (2019).
  * Cirac  _et al._ [2021] J. I. Cirac, D. Pérez-García, N. Schuch, and F. Verstraete, Reviews of Modern Physics 93, 045003 (2021).
  * Vidal [2003] G. Vidal, Physical Review Letters 91, 147902 (2003).
  * Bañuls  _et al._ [2006] M. C. Bañuls, R. Orús, J. I. Latorre, A. Pérez, and P. Ruiz-Femenía, Physical Review A 73, 022344 (2006).
  * Dang  _et al._ [2019] A. Dang, C. D. Hill, and L. C. L. Hollenberg, Quantum 3, 116 (2019).
  * Kloss  _et al._ [2020] B. Kloss, D. Reichman, and Y. Bar Lev, SciPost Physics 9, 070 (2020).
  * Evenbly and Vidal [2014] G. Evenbly and G. Vidal, Journal of Statistical Physics 157, 931 (2014).
  * Zaletel and Pollmann [2020] M. P. Zaletel and F. Pollmann, Physical Review Letters 124, 037201 (2020).
  * Lin  _et al._ [2022] S.-H. Lin, M. P. Zaletel, and F. Pollmann, Physical Review B 106, 245102 (2022).
  * Noh  _et al._ [2020] K. Noh, L. Jiang, and B. Fefferman, Quantum 4, 318 (2020).
  * Oh  _et al._ [2021] C. Oh, K. Noh, B. Fefferman, and L. Jiang, Physical Review A 104, 022407 (2021), arXiv:2101.11234 .
  * Verstraete  _et al._ [2004] F. Verstraete, J. J. García-Ripoll, and J. I. Cirac, Physical Review Letters 93, 207204 (2004), arXiv:0406426v2 [cond-mat] .
  * Zwolak and Vidal [2004] M. Zwolak and G. Vidal, Physical Review Letters 93, 207205 (2004).
  * Prosen and Žnidarič [2009] T. Prosen and M. Žnidarič, Journal of Statistical Mechanics: Theory and Experiment 2009, P02035 (2009).
  * Guth Jarkovský  _et al._ [2020] J. Guth Jarkovský, A. Molnár, N. Schuch, and J. I. Cirac, PRX Quantum 1, 010304 (2020), arXiv:2003.12418 .
  * De Las Cuevas and Netzer [2020] G. De Las Cuevas and T. Netzer, Journal of Mathematical Physics 61, 041901 (2020), arXiv:1907.03664 .
  * Rath  _et al._ [2023] A. Rath, V. Vitale, S. Murciano, M. Votto, J. Dubail, R. Kueng, C. Branciard, P. Calabrese, and B. Vermersch, PRX Quantum 4, 010318 (2023).
  * Jaschke  _et al._ [2018] D. Jaschke, S. Montangero, and L. D. Carr, Quantum Science and Technology 4, 013001 (2018).
  * Cuevas  _et al._ [2013] G. D. l. Cuevas, N. Schuch, D. Pérez-García, and J. I. Cirac, New Journal of Physics 15, 123021 (2013), arXiv:1308.1914 [cond-mat, physics:math-ph, physics:quant-ph] .
  * Werner  _et al._ [2016] A. H. Werner, D. Jaschke, P. Silvi, M. Kliesch, T. Calarco, J. Eisert, and S. Montangero, Physical Review Letters 116, 237201 (2016), arXiv:1412.5746 [cond-mat, physics:quant-ph] .
  * Hauschild  _et al._ [2018] J. Hauschild, E. Leviatan, J. H. Bardarson, E. Altman, M. P. Zaletel, and F. Pollmann, Physical Review B 98, 235163 (2018), arXiv:1711.01288 [cond-mat] .
  * Pokart  _et al._ [2023] T. Pokart, C. Lehmann, and J. C. Budich, (2023), arXiv:2309.07961 [cond-mat.str-el] .
  * Surace  _et al._ [2019] J. Surace, M. Piani, and L. Tagliacozzo, Physical Review B 99, 235115 (2019).
  * Guo and Yang [2023a] Y. Guo and S. Yang,  (2023a), arXiv:2307.16381 [quant-ph] .
  * Schollwöck [2013] U. Schollwöck, in  _Strongly Correlated Systems: Numerical Methods_, Springer Series in Solid-State Sciences, edited by A. Avella and F. Mancini (Springer, Berlin, Heidelberg, 2013) pp. 67–98.
  * Nguyen  _et al._ [2018] P. Nguyen, T. Devakul, M. G. Halbasch, M. P. Zaletel, and B. Swingle, Journal of High Energy Physics 2018, 98 (2018).
  * Karrasch  _et al._ [2012] C. Karrasch, J. H. Bardarson, and J. E. Moore, Physical Review Letters 108, 227206 (2012).
  * Karrasch  _et al._ [2013] C. Karrasch, J. H. Bardarson, and J. E. Moore, New Journal of Physics 15, 083031 (2013).
  * Barthel [2013] T. Barthel, New Journal of Physics 15, 073010 (2013).
  * White  _et al._ [2017] C. D. White, M. Zaletel, R. S. K. Mong, and G. Refael, Physical Review B 97, 035127 (2017), arXiv:1707.01506 .
  * Cheng  _et al._ [2021] S. Cheng, C. Cao, C. Zhang, Y. Liu, S.-Y. Hou, P. Xu, and B. Zeng, Phys. Rev. Res. 3, 023005 (2021).
  * Guo [2022] C. Guo, Physical Review B 105, 195152 (2022).
  * Bridgeman and Chubb [2017] J. C. Bridgeman and C. T. Chubb, Journal of Physics A: Mathematical and Theoretical 50, 223001 (2017).
  * Rényi [1961] A. Rényi, in  _Proceedings of the Fourth Berkeley Symposium on Mathematical Statistics and Probability, Volume 1: Contributions to the Theory of Statistics_, Vol. 4.1 (University of California Press, 1961) pp. 547–562.
  * Li  _et al._ [2019] G. Li, Y. Ding, and Y. Xie,  (2019), arXiv:1809.02573 .
  * Martiel and de Brugière [2020] S. Martiel and T. G. de Brugière,  (2020), arXiv:2012.09663 .
  * Renes [2022] J. M. Renes,  _Quantum Information Theory_ (De Gruyter Oldenbourg, Berlin, Boston, 2022).
  * Uhlmann [1976] A. Uhlmann, Reports on Mathematical Physics 9, 273 (1976).
  * Wang  _et al._ [2008] X. Wang, C.-S. Yu, and X. Yi, Physics Letters A 373, 58 (2008).
  * Guo and Yang [2023b] Y. Guo and S. Yang,  (2023b), arXiv:2312.02854 [quant-ph] .
  * Terhal  _et al._ [2002] B. M. Terhal, M. Horodecki, D. W. Leung, and D. P. DiVincenzo, Journal of Mathematical Physics 43, 4286 (2002).

##  Appendix A Derivation of the gradient for gauge optimizations

In this section we compute the gradient used to solve the optimization problem
in Eq. (LABEL:eq:entropy-opt). We adopt the notation of Fig. 3 where $A$
denotes the contraction of two adjacent sites in mixed canonical form and
$\tilde{A}$ is the contraction of $A$ with the gauge unitary $U$.

Letting $X=\tilde{A}^{\dagger}\tilde{A}$, the chain rule reads

$\frac{\partial E}{\partial U}=\sum_{ij}\frac{\partial X_{ij}}{\partial
U}\frac{\partial E}{\partial X_{ij}},$ (14)

where $\frac{\partial X_{ij}}{\partial U}$ can be obtained via the tensor
expression in Fig. A.1. Below we compute $\frac{\partial E}{\partial U}$ for
the Rényi entropies in Eq. (4) with $\alpha=1$ and $\alpha=2$.

###  A.1 Rényi entropy with $\alpha=1$

The 1-Rényi entropy is given by

$E(X)=-\Tr{X\ln X}.$ (15)

Differentiating with respect to elements of $X$ yields

$\frac{\partial E(X)}{\partial X_{ij}}=-\Tr{\frac{\partial X}{\partial
X_{ij}}\ln X+X\frac{\partial\ln X}{\partial X_{ij}}}.$ (16)

The first term is straightforward to simplify, since

$\left(\frac{\partial X}{\partial
X_{ij}}\right)_{k\ell}=\delta_{ik}\delta_{j\ell}.$ (17)

One then obtains

$\displaystyle\Tr{\frac{\partial X}{\partial X_{ij}}\ln X}$
$\displaystyle=\sum_{k}\sum_{\ell}\left(\frac{\partial X}{\partial
X_{ij}}\right)_{k\ell}(\ln X)_{\ell k}=(\ln X)_{ji}.$ (18)

We now turn to the second term of (16). For a positive definite matrix $M$ of
unit trace, the logarithm is well defined and the series

$\ln X=\sum_{n=1}^{\infty}\frac{(-1)^{n-1}}{n}(X-\mathds{1})^{n}$ (19)

converges since $\|X-\mathds{1}\|_{\infty}<1$ . By the product rule we have

$\frac{\partial}{\partial
X_{ij}}(X-\mathds{1})^{n}=\sum_{k=0}^{n-1}(X-\mathds{1})^{k}\frac{\partial
X}{\partial X_{ij}}(X-\mathds{1})^{n-1-k},$ (20)

and using the cyclic property of the trace, we obtain

$\Tr{X\frac{\partial\ln X}{\partial X_{ij}}}=\Tr{X\frac{\partial X}{\partial
X_{ij}}\sum_{n=0}^{\infty}(\mathds{1}-X)^{n}}.$ (21)

Moreover, since $X=\mathds{1}-(\mathds{1}-X)$, we obtain a telescopic sum, so
that only the $n=0$ term survives, and

$\Tr{X\frac{\partial\ln X}{\partial X_{ij}}}=\Tr{\frac{\partial X}{\partial
X_{ij}}}=\delta_{ij}.$ (22)

Finally,

$\frac{\partial E(X)}{\partial X}=-\left[(\ln X)^{T}+\mathds{1}\right].$ (23)

Figure A.1: Tensor network expression of the partial derivative of
$X=\tilde{A}^{\dagger}\tilde{A}$ with respect to the unitary $U$.

(a)

(b)

Figure A.2:  Convergence analysis for the one-shot purification
disentanglement of a random MPDO shown in Fig. 5. The convergence of local
bond entropies and the variation of singular value distributions on each bond
over the number of sweeps are shown for entanglement bonds in panel (a) and
purity bonds in panel (b).

###  A.2 Rényi entropy with $\alpha=2$

We now compute the partial derivative of the 2-Rényi entropy

$E_{2}(X)=-\ln\Tr{X^{2}}.$ (24)

with respect to elements of $X$. First we note that

$\frac{\partial}{\partial
X_{ij}}\left(-\ln\Tr{X^{2}}\right)=-\frac{1}{\Tr{X^{2}}}\Tr{\frac{\partial
X^{2}}{\partial X_{ij}}}.$ (25)

Then, by cyclicity

$\Tr{\frac{\partial X^{2}}{\partial X_{ij}}}=2\Tr{\frac{\partial X}{\partial
X_{ij}}X},$ (26)

and

$\Tr{\frac{\partial X}{\partial X_{ij}}X}=\sum_{k\ell}\frac{\partial
X_{k\ell}}{\partial X_{ij}}X_{\ell
k}=\sum_{k\ell}\delta_{ik}\delta_{j\ell}X_{\ell k}=X_{ji},$ (27)

from which

$\frac{\partial E_{2}(X)}{\partial X}=-\frac{2X^{T}}{\Tr{X^{2}}}.$ (28)

##  Appendix B Details of convergence in the one-shot scenario

Here we study in more detail the optimization performed in Fig. 5 and discuss
convergence of the IPD technique.

Figure A.3:  Convergence study of the procedure of Fig. 5. We show the
infidelity of the $s$-th sweep for the state (plain blue line, left-hand-side
scale) and for its purification (dashed orange line, right-hand-side scale).
$\ket{\Psi_{s}}$ and $\rho_{s}$ denote respectively the purification and the
state obtained after the $s$-th sweep.

Fig. A.2 shows a more complete optimization report of the one-shot
purification disentanglement procedure performed in Fig. 5. We show the
variation of bond entropies (of order $\alpha=1$) $\mathcal{E}_{1}^{i}$ and
$\mathcal{P}_{1}^{i}$ (defined in Eq. (7)) with the number of sweeps $s$, and
the bond singular value distributions $\\{(\sigma_{j}^{(i)})^{2}\\}_{j}$ and
$\\{\lambda_{j}^{(i)}\\}_{j}$ after different numbers of sweeps. We observe
that bond entropies, for both entanglement bonds and purification bonds,
decrease with the number of sweeps and seem to converge to non-zero values.

One would be justified in asking whether this iterative disentanglement
technique converges to one final purified state. To answer this we repeat the
computations reported in Fig. 5, this time without truncating purification
bonds between sweeps. This allows us to keep the dimension of the ancilla
system constant and study the variation of the purified state. In Fig. A.3 we
observe that the distance between purifications obtained after successive
sweeps, quantified using the fidelity, reliably decreases to zero (dashed
orange line). After about 15-20 sweeps, the local disentangling operations
converge in one or a few gradient descent steps, and the procedure has thus
converged to a locally optimized purification.

Fig. A.3 also records the relative error of the state due to the necessary
truncation of entanglement bonds. The behavior for the first 10-15 sweeps is
similar to the purification relative distance measure, and can be explained by
the fact that early sweeps tend to increase entanglement bond dimensions
significantly, yielding a non-negligible error upon truncation. Moreover, we
observe a peak in the value of $1-\mathcal{F}[\rho_{s},\rho_{s-1}]$ in the
range of 15-20 sweeps. This can be identified in Fig. A.2 with the region in
which a significant number of singular values on entanglement bonds fall below
the truncation threshold, thus accounting for most of the truncation error.

Figure A.4: Comparison of $\alpha=1$ and $\alpha=2$ cost functions in the one-
shot disentanglement scenario. We plot the exact compression fidelity (top)
and the memory usage (bottom) in both cases.

##  Appendix C Comparison of performance with different entropy metrics

In this section we justify numerically our choice of using the $\alpha=2$
entanglement bond entropy as the objective function in Eq. (LABEL:eq:entropy-
opt), as opposed to $\alpha=1$.

In Fig. A.4 we compare the results of the one-shot disentanglement compression
of a random quantum state, using as an objective function the $\alpha=1$ (blue
line) and $\alpha=2$ (orange line) entanglement bond entropies.

The quantum state is generated by simulating 6 layers of a 6-qubit random
circuit with depolarizing noise level $p=5\%$. At each layer both entanglement
bonds and purifications bond are held fixed by truncated SVD, at dimensions
$\chi=r=8$. In the purification disentanglement sweeps we used a truncation
threshold of $\varepsilon=10^{-3}$ and a maximum allowed bond dimension
$\chi_{max}=64$. Local optimizations were implemented with a conjugate
gradient algorithm with tolerance $10^{-5}$ and a maximum of 500 iterations.

We observe that locally minimizing the $\alpha=2$ entropies allow one to
converge to a low-memory solution faster and with much better fidelity that
minimizing $\alpha=1$ entropies. A closer study of singular value
distributions across sweeps (not shown) reveals that the $\alpha=1$ entropy
optimization results in a faster decay of entanglement bond dimensions, but
does not exhibit the “plateau descent” of singular values on purification
bonds observed in Fig. A.2 with $\alpha=2$. A similar behavior was observed
for other parameters of the quantum state generation or the disentangling
procedure, which also translated to generally better performances for
$\alpha=2$ optimizations in circuit simulations.

##  Appendix D Random noisy circuits

To benchmark our circuit emulators we use random quantum circuits. We consider
similar circuits as in Ref. [34]. A circuit layer consists of a layer of one-
qubit gates followed by a layer of nearest-neighbor (in a 1D topology) two-
qubit gates, alternating as illustrated in Fig. A.5. Such circuits were also
used for generating random MPDOs in Sec. III.2 and App. B and C, with
truncation to fixed bond dimensions after each circuit layer.

Figure A.5: Illustration of three layers of a random quantum circuit as used
in [34]. Figure A.6:  Convergence to a maximally mixed state (MMS) at large
depths for $N=8$ qubits random circuits. This is quantified at various
depolarizing noise levels $p$ with three quantities: purity (a), fidelity to
the MMS (b) and OEE (c). Values are averaged over 16 random circuits and error
bars show the corresponding standard deviation.

The one-qubit gates are parametrized rotations defined by the unitary matrix
$\exp{i\alpha(\sigma_{x}\sin\theta\cos\phi+\sigma_{y}\sin\theta\sin\phi+\sigma_{z}\cos\theta)}$,
where the parameters $\alpha$, $\theta$, $\phi$ are chosen uniformly at random
in $[0,2\pi)$ for each gate, and $\sigma_{x},\sigma_{y},\sigma_{z}$ are the
three Pauli matrices. The two-qubit gates are chosen to be either controlled-
NOT (CNOT) or controlled-Z (CZ) gates with equal probability. A depolarizing
noise channel,

$\rho\longrightarrow(1-p)\rho+p\frac{\mathds{1}}{2},$ (29)

with noise level $p$ is applied qubit-wise after each layer.

##  Appendix E Proximity to the maximally mixed state (MMS)

As depolarizing noise leads ultimately to the maximally mixed state (MMS),
i.e. $\rho_{\rm MMS}=\mathds{1}/2^{N}$, we want to understand if the
performances of our compression is caused by the proximity to this state.
Indeed, the MMS is a trivial state to represent with an MPDO of purification
bond dimensions $r_{i}=2$ and entanglement bond dimensions $\chi_{i}=1$. One
could question if the reduction of purification bonds we obtain is simply a
demonstration that a MMS has been reached.

In Fig. A.6, we quantify the proximity to the MMS along a noisy random circuit
through three quantities: purity, fidelity against the MMS and OEE. For the
MMS, these should be zero, one and zero respectively. Exact quantum states for
each depth are computed and subsequently mapped to a MPO via successive SVDs
to measure the OEE.

For noise level $p=10\%$ (blue line), about ten layers are needed to reach the
MMS. However, we see in Fig 7 that the IPD yields purification bond dimensions
$r_{i}\approx 2$ even for the 10 first layers. This means the compression of
the purification bond dimensions was obtained for states that are far from the
MMS, with substantial purity and OEE.

Performing the same analysis with the $p=2\%$ case (green line), and comparing
to Fig. 8, gives the same conclusion. This time, purification bond dimensions
$r_{i}\approx 2$ is reached after about 15 layers, at which the OEE and
infidelity to the MMS are not negligible. We note that during the first
layers, where purification bond dimensions are a bit higher ($r_{i}\leq 5$),
purity is also larger. This suggests a relation between optimized purification
bond dimension and purity.

Note also that the random state used in App. C has purity $\approx 0.1$ and
fidelity $\approx 0.5$ with respect to the MMS.

##  Appendix F List of naming conventions for MPOs and MPDOs

The literature about MPDOs and MPOs produced a surprising amount of different
and sometimes mutually exclusive naming conventions. To help the new reader
navigate through them, we compiled a list of those we encountered in Table 1.

Reference | name for MPDO | name for MPO
---|---|---
Cheng  _et al._ [34] | MPDO |
De Las Cuevas and Netzer [19] | LPF222Locally purified form, not an acronym in the text. | MPDO / MPO
Guo [35], Surace  _et al._ [26] | MPDO | MPO
Guo and Yang [27], Guo and Yang [43] | LPDO | MPO
Werner  _et al._ [23] | LPTN | MPO
Jaschke  _et al._ [21] | LPTN | MPDO
Hauschild  _et al._ [24] | purification | MPO
Cuevas  _et al._ [22] | local purification | MPDO
Verstraete  _et al._ [15] | MPDO | MPDO
Prosen and Žnidarič [17], Oh  _et al._ [14], Noh  _et al._ [13] |  | MPO
Zwolak and Vidal [16] |  | MPD
White  _et al._ [33], Guth Jarkovský  _et al._ [18] |  | MPDO
Table 1: Guide for the reader of the MPDO and MPO literature.

  *[DMRG]: density matrix renormalization group
  *[MPDO]: matrix product density operator
  *[MPDOs]: matrix product density operator
  *[NISQ]: noisy intermediate scale quantum
  *[TN]: tensor network
  *[MPS]: matrix product state
  *[TNs]: tensor network
  *[MPOs]: matrix product operator
  *[MPO]: matrix product operator
  *[OEE]: operator entanglement entropy
  *[SVD]: singular value decomposition
  *[DMT]: density matrix truncation
  *[LC]: local compression
  *[IPD]: iterative purification disentanglement
  *[MMS]: maximally mixed state
  *[SVDs]: singular value decomposition