# Optimal sampling of tensor networks targeting wave function’s fast decaying
tails

Marco Ballarin  0000-0002-3215-3453 Dipartimento di Fisica e Astronomia "G.
Galilei" & Padua Quantum Technologies Research Center, Università degli Studi
di Padova, Italy I-35131, Padova, Italy  INFN, Sezione di Padova, via Marzolo
8, I-35131, Padova, Italy  Pietro Silvi  0000-0001-5279-7064 Dipartimento di
Fisica e Astronomia "G. Galilei" & Padua Quantum Technologies Research Center,
Università degli Studi di Padova, Italy I-35131, Padova, Italy  INFN, Sezione
di Padova, via Marzolo 8, I-35131, Padova, Italy  Simone Montangero
0000-0002-8882-2169 Dipartimento di Fisica e Astronomia "G. Galilei" & Padua
Quantum Technologies Research Center, Università degli Studi di Padova, Italy
I-35131, Padova, Italy  INFN, Sezione di Padova, via Marzolo 8, I-35131,
Padova, Italy  Institute for Complex Quantum Systems, Ulm University, Albert-
Einstein-Allee 11, 89069 Ulm, Germany  Daniel Jaschke  0000-0001-7658-3546
Dipartimento di Fisica e Astronomia "G. Galilei" & Padua Quantum Technologies
Research Center, Università degli Studi di Padova, Italy I-35131, Padova,
Italy  INFN, Sezione di Padova, via Marzolo 8, I-35131, Padova, Italy
Institute for Complex Quantum Systems, Ulm University, Albert-Einstein-Allee
11, 89069 Ulm, Germany

###### Abstract

We introduce an optimal strategy to sample quantum outcomes of local
measurement strings for isometric tensor network states. Our method generates
samples based on an exact cumulative bounding function, without prior
knowledge, in the minimal amount of tensor network contractions. The algorithm
avoids sample repetition and, thus, is efficient at sampling distribution with
exponentially decaying tails. We illustrate the computational advantage
provided by our optimal sampling method through various numerical examples,
involving condensed matter, optimization problems, and quantum circuit
scenarios. Theory predicts up to an exponential speedup reducing the scaling
for sampling the space up to an accumulated unknown probability $\epsilon$
from $\mathcal{O}(\epsilon^{-1})$ to $\mathcal{O}(\log(\epsilon^{-1}))$ for a
decaying probability distribution. We confirm this in practice with over one
order of magnitude speedup or multiple orders improvement in the error
depending on the application. Our sampling strategy extends beyond local
observables, e.g., to quantum magic.

In the era of the early development of quantum processors and simulators
Monroe (2002); Monroe and Kim (2013); Saffman (2018); Wendin (2017); Bruzewicz
_et al._ (2019); Adams  _et al._ (2019); Kaushal  _et al._ (2020); Kjaergaard
_et al._ (2020); Morgado and Whitlock (2021), a high demand arises for
numerical emulators of programmable quantum devices: classical algorithms
capable of replicating the input, real-time processing, and ouput of a quantum
machine at small and intermediate scales Nagy and Savona (2019); Qiskit
contributors (2023); The cuQuantum development team (2023); Efthymiou  _et
al._ (2022). Such demand is motivated by benchmarking and certification of the
programmable quantum device Jaschke  _et al._ (2022) as well as exploring the
extent of quantum advantage: classical simulation of quantum machines has
shown remarkable progress in recent experiments Arute  _et al._ (2019); Zhong
_et al._ (2020) and new numerical algorithms have been developed to overcome
existing obstacles Zhou  _et al._ (2020); Liu  _et al._ (2021); Bulmer  _et
al._ (2022); Daley  _et al._ (2022); Oh  _et al._ (2023).

Tensor networks are often employed to numerically characterize intermediate-
scale quantum systems and devices. As they store spatial correlations as
compressed information, they are an invaluable tool for describing many-body
quantum systems with limited entanglement content. Various classes of tensor
networks have been developed so far, each adapted to a specific correlation
geometry: from matrix product states (MPS) White (1992); Rommer and Östlund
(1997, 1999); Schollwöck (2005); Perez-Garcia  _et al._ (2007); McCulloch
(2007); Schollwöck (2011) via tree tensor networks (TTN) Shi  _et al._ (2006);
Tagliacozzo  _et al._ (2009); Murg  _et al._ (2010); Gerster  _et al._ (2014)
to projected entangled pair states (PEPS) F.  _et al._ (2008); Orús (2014),
and more Zwolak and Vidal (2004); Verstraete  _et al._ (2004); Rizzi  _et al._
(2008); Cincio  _et al._ (2008); McCulloch (2008); Jordan  _et al._ (2008);
Evenbly and Vidal (2009, 2014); Phien  _et al._ (2015); Werner  _et al._
(2016); Silvi  _et al._ (2019); Montangero (2018); Jaschke  _et al._ (2018);
Reinić (2022); Arceci  _et al._ (2022); Bañuls (2023); Okunishi  _et al._
(2022). Also in the context of sampling, tensor network methods have emerged
as a powerful tool: Indeed, tensor network allow efficient output data
acquisition from quantum states and processes Gomez  _et al._ (2022); Yang
_et al._ (2023); Dymarsky and Pavlenko (2022). Sampling of tensor network
states has also been used to determine quantities that are otherwise difficult
to measure, such quantum state magic Lami and Collura (2023), entanglement
characterization Notarnicola  _et al._ (2023), or to sample from distributions
via quantum states Wild  _et al._ (2021). Despite considerable recent
advancements Han  _et al._ (2018); Bravyi  _et al._ (2022), we identified a
sampling strategy meant to outpower known methods.

In this work, we present an efficient sampling procedure, which we label
Optimal tEnsor network Sampling (OPES), applicable to any isometric tensor
network state Zaletel and Pollmann (2020). The algorithm itself falls into the
category of inverse transform sampling algorithms Niederreiter ; Devroye
(2006); Givens and Hoeting (2012); as additional challenge in comparison to
the inverse transform sampling, we have no initial knowledge of the cumulative
probability distribution or its inverse, which even by the end is known only
partially for the many-body quantum state. Our approach exploits the
capability of isometric tensor networks to efficiently acquire exact
conditional probabilities for partial outcome strings, and carefully avoids
re-sampling of previously explored outcomes. Optimal tensor network sampling
method is thus particularly suited to characterize peaked, non-flat
distributions, where we observe a speedup beyond one order of magnitude over
the standard sampling. The speedup to be expected depends foremost on the
probability distribution, then on the number of samples or coverage to be
achieved. To demonstrate the usefulness of the approach, we consider two
popular classes of isometric tensor networks, namely MPS and TTN. We identify
a substantial speedup both for sampling a fixed number of states as well as
for covering a target portion of the probability space; while the exact
probability of non-sampled space is known at all times. The OPES algorithm
benefits from storing in memory all the partial tensor network contractions,
which characterize the conditional probabilities, evaluated until that
instant, to avoid repeating calculations. We analyze the advantage offered by
the OPES method in five practical examples from different quantum
applications, such as image representation, condensed matter or measurement of
quantum magic. The examples include pure and mixed states. The developed code
is open source and included in the _Quantum TEA_ Ballarin  _et al._ (2024a)
library.

This manuscript is structured as follows. First, we focus on the idea
underlying OPES in Sec. I, then on the numerical implementations in Secs. II
and III. In Sec. IV, we introduce classes of different quantum states, each
targeting different applications; we sample with the two methods from these
states and analyze the speedup that one can expect in each application.
Finally, we summarize our results in Sec. V. The appendices contain the
calculation of the exponential speedup of the scaling in App. .1, the full
algorithm in App. .2, as well as additional implementation strategies for
tensor networks, see App. .3.

##  I Theory and standard sampling from tensor network states

Figure 1:  _Different ways of sampling a given probability distribution
function encoded in an $n=3$ qubit state._ a) The motivation to sample by
shots are experiments, where the probabilities are proportional to the
normalized number of shots for a given state and converge in the limit of
large numbers. Thus, the probability bar in green is proportional to the shots
shown as red circles. b) The actual probability distribution function is the
reference distribution for any sampling method. c) One can rewrite the exact
probabilities into a cumulative probability function which highlights that
each state $\alpha$ can be associated with an interval $\delta_{\alpha}$, see
also Eq. (2), where the intervals for the subsystem-state $\ket{0}_{0}$ and
state $\ket{100}$ are indicated as $\delta_{0}$ and $\delta_{100}$. The
intervals are defined via the left and right boundaries of the intervals,
i.e., $\delta_{\alpha}=(b_{\alpha}^{l},b_{\alpha}^{r})$. The probability
stemming from the previous states in the order is shown in black, the
probability of the state itself is shown in cyan. d) The standard sampling
follows a shot-by-shot approach to sample states $\alpha$, but has access to
the exact probabilities $p_{\alpha}$. Additionally, probabilities of states
differing only in the last qubit of the bitstring are known at no additional
cost, see non-zero bars without shot (red circle). The probability of the
green bars do not sum to one, which accounts for the unknown probability
space. e) Using OPES, we avoid resampling and thus all shots are for different
states; in detail, red circles fall into different bars. Moreover, pairs of
probabilities are discovered together analog to _d)_ and even the substates of
the the first $(n-1)$ qubits of the shots in red are unique for each shot.
Four steps are sufficient to sample a space of eight states, which generalizes
to an upper bound of $\mathcal{O}(2^{n-1})$ samples to scan a full Hilbert
space. The distribution of the shots sampled in this way does not converge to
the true distribution anymore as there is either zero or one shot per bin; one
always uses the exact probabilities. f) The OPES scaling is achieved by
removing intervals from the space where the algorithm is allowed to sample;
only previously not sampled states remain as potential next samples. While
running the iterations top-down in this illustration, we gain partial
knowledge of the cumulative probability distribution of _b)_ , i.e., of the
intervals of the full state when reaching the green checkmark as well as
subsystem intervals as indicated. For example in step $n_{\rm r}=1$, we
identify the intervals $\delta_{0}$ and $\delta_{11}$ for subsystems. The
intervals with green checkmarks are then removed for all the following
iterations, see red crosses.

Although we compare exclusively the performance of the standard algorithm and
the OPES algorithm for tensor networks in the following, the motivation
originates in experimental realizations and how sampling works there. Figure
1a illustrates how experiments sample from the underlying probability
distribution by running multiple shots for three qubits. The ratio of the
number of shots for a state $\alpha$ over the total number of shots converges
to the exact distribution in the limit of large number of shots. Due to the
central limit theorem Fischer (2011), the histogram approximates the actual
probability distribution. The basis states $\alpha$ forming the projection in
$z$-direction are an intuitive choice for the experiment, which for example
can be sampled in an atom-qubit experiment with a fluorescence snapshot Monroe
(2002); Saffman (2018); Morgado and Whitlock (2021). Similarly to the
experimental results, both tensor network sampling methods can generate the
number of shots per state $\alpha$.

For simulations and theory, we encounter more possiblities to work with the
actual underlying probability distribution of the quantum state, i.e.,
depicted in Fig. 1b. While sampling generates the outcome of simultaneous
measured local observables of qudits, we focus on qubits as a building block
for the many-body wave function and a projection onto the Pauli operator
$\sigma^{z}$ without loss of generality. The possible measurement outcomes on
a single qubit are its logical states $\alpha\in\\{\ket{0},\ket{1}\\}$. To
build up the notation for the many-body case with a tensor network state, we
consider first probabilities arising from the wave function $\ket{\psi}$ for
$n$ qubits, i.e.,

$\displaystyle\ket{\psi}$
$\displaystyle=\sum_{q_{0},\dots,q_{n-1}=0}^{1}c_{q_{0},\dots,q_{n-1}}\ket{q_{0}\dots
q_{n-1}}=\sum_{\alpha=0}^{2^{n}-1}c_{\alpha}\ket{\alpha}\,.$ (1)

The sampling algorithms aim to reproduce the probabilities
$p_{\alpha}=\left|c_{\alpha}\right|^{2}$. We label the basis states $\alpha$
either with the binary number $q_{0}\dots q_{n-1}$ or its integer
representation. While we focus on a pure state in Eq. (1), the concept
evidently holds as well for mixed states. The order of the basis states
$\alpha$ on the $x$-axis in Fig. 1b is used throughout the manuscript.

Based on this given order, the cumulative probability distribution is
generated and shown in Fig. 1c. We obtain the ordered probability intervals,
where we define the probability $p_{\alpha}=\left|c_{\alpha}\right|^{2}$ as
the width of an interval $\delta_{\alpha}$ as

$\displaystyle\delta_{\alpha}$
$\displaystyle=\left(b^{l}_{\alpha},b^{l}_{\alpha}+|c_{\alpha}|^{2}\right),$
(2) $\displaystyle b^{l}_{\alpha}$
$\displaystyle=\sum_{\beta<\alpha}|c_{\beta}|^{2}\,.$ (3)

The boundaries of each $\delta_{\alpha}$ are well defined within the complete
probability interval $(0,1)$. Thus, the left boundary $b^{l}_{\alpha}$ of the
probability interval associated with the state $\alpha$ can be defined as
reported in Eq. (3). We stress that, as shown in Eq. (2), the right boundary
is just $b^{r}_{\alpha}=b^{l}_{\alpha}+|c_{\alpha}|^{2}$. Inverse transform
sampling algorithms can be used if these intervals are all known Niederreiter
; Devroye (2006); Givens and Hoeting (2012). Then, one can sample a uniform
random number in the interval $(0,1)$, which serves as value on the $y$-axis
of Fig. 1c, and finds the corresponding state $\alpha$ on the $x$-axis via the
inverse cumulative probability distribution. While the inverse cumulative
probability distribution is easily accessible once all $p_{\alpha}$ are known,
they are initially unknown in the tensor network scenario and the number of
states $\alpha$ grows exponentially in system size for many-body systems. But
we can already anticipate based on Eq. (3) that neither all
$\\{|c_{\beta}|^{2}\\}_{\beta\neq\alpha}$ nor even
$\\{|c_{\beta}|^{2}\\}_{\beta<\alpha}$ need to be known. Partial knowledge can
be sufficient to attribute a random number to the corresponding interval and
its state $\alpha$, e.g., knowing the Eq. (3) for an intermediate state with
less than $n$ qubits, i.e., the intervals $\delta_{\alpha^{\prime}}$ with
$\alpha^{\prime}$ being a subsystem of $\alpha$ containing the first $m<n$
qubits; we explain the details in Sec. III.

The standard approach to sample outcomes
$\left(|c_{\alpha}|^{2},\alpha\right)$ is the generation of bitstrings shot by
shot. For each shot, the standard method consecutively projects the qubits
into their basis state Nielsen and Chuang (2007). As we calculate
probabilities shot-by-shot, we obtain the shot count for each bitstring
$\alpha$ at the end of the standard sampling as shown in Fig. 1d. Instead of
calculating the probability as ratio of shot count over number of shots,
tensor networks have direct access to a pair of probabilities $p_{\alpha}$
with each shot. The pair always differs only in the measurement outcome of the
last qubit. Measuring the system qubit by qubit, we accumulate a product of
conditional probabilities

$\displaystyle|c_{\alpha}|^{2}=p(q_{n-1}|q_{n-2}\dots q_{1}q_{0})\dots
p(q_{1}|q_{0})p(q_{0}),$ (4)

i.e., the easily accessible information about the couple
$\left(|c_{\alpha}|^{2},\alpha\right)$. The probability for only the last
qubit being in the other state depends on the same probabilities of the state
$\ket{q_{n-2}\dots q_{1}q_{0}}$ and is obtained for free. Therefore, Fig. 1d
shows the distribution of shots, the exact probability for all states measured
via shots from Eq. (4), as well as the exact probabilities of some states
without shots. The latter appear due to differing only in the last qubit from
a state known due to a shot, e.g., $\ket{010}$ being obtained together with
$\ket{011}$.

We observe that the standard sampling in Fig. 1d repeats the sampling of the
same bitstring multiple times, while the exact probabilities already known,
see state $\ket{100}$. We summarize the problems in the standard sampling
approach that the OPES method overcomes as follows:

  1. 1.

As the trial bitstrings $\ket{q_{0}\dots q_{n-1}}$ via the tensor tensor
network measurement yield the probability from the exact distribution upon the
first encounter, each of them might be sampled repeatedly without gaining any
additional information;

  2. 2.

Similarly, it is statistically unlikely to sample marginal outcomes, thus the
exponentially small tails of peaked distribution either remain largely
unexplored or need an exponential time to be sampled.

The fundamental idea of OPES is highlighted in Fig. 1e and f; details of OPES
follow in Sec. III. Figure 1e shows the complete knowledge of all
probabilities after four shots. Because we can gain the probability of two
states with each measurement of a tensor network, the optimal number of
necessary measurements scales linearly with the size of the Hilbert space,
i.e., as $\mathcal{O}(2^{n-1})$. To reach this scaling, we permanently remove
the sampled outcome strings from the probability measure space in a way that
they cannot be sampled again, as shown in Fig. 1f. After each step, we thus
remove two intervals. Effectively, we renormalize the sampling probabilities
of the remaining, missing outcomes. Thus, the sampling probability of a not
yet sampled string increases at each shot. Let us suppose we just sampled the
couple $(|c_{\alpha^{*}}|^{2},\alpha^{*})$ with probability $p_{\alpha^{*}}$
as the first sample. Then, the probabilities of sampling any
$\alpha\neq\alpha^{*}$ at the following iteration are

$\displaystyle
p_{\alpha\neq\alpha^{*}}\rightarrow\frac{1}{1-p_{\alpha^{*}}}p_{\alpha\neq\alpha^{*}}>p_{\alpha\neq\alpha^{*}}\,,$
(5)

and so forth. Thus, OPES increases the probability to draw a bitstring at
every iteration which results in a faster converging strategy.

We observe an exponential speedup reducing the scaling for sampling from
$\mathcal{O}(\epsilon^{-1})$ to $\mathcal{O}(\log(\epsilon^{-1}))$ for a
decaying probability distribution, with $\epsilon^{-1}$ being the cumulative
unknown probability. For details, see App. .1.

##  II Standard tensor network sampling

We briefly summarize the algorithm for the standard tensor network sampling,
which generates a single shot for two qubits as shown in Fig. 2a. At the end
of the algorithm, we have obtained the conditional probabilities of Eq. (4).
First, we draw a vector of random
number$\vec{u}=\left(u_{0},u_{1},\dots,u_{n-1}\right)\sim
U([0,1]),\;\vec{u}\in\mathbb{R}^{n}$ uniformly distributed in the interval
$[0,1]$. Then, we compute the probability $p_{0}$ of measuring $\ket{0}$ on
the qubit $0$, e.g., via the reduced density matrix $\rho_{0}$. Next, we
project $q_{0}$ of the state $\ket{\psi}$ in the measured subspace based on
the random number $u_{0}$,

$\displaystyle q_{0}=\begin{cases}\ket{0}&\text{if }u_{0}<p_{0}\,,\\\
\ket{1}&\text{if }u_{0}\geq p_{0}\,.\end{cases}$ (6)

Afterward, we repeat the procedure for the second qubit $q_{1}$ with the
random number $u_{1}$ computing the conditional probability
$p(q_{1}=\ket{0}|q_{0}=\ket{0})=p_{0|0}$ of measuring $q_{1}$ in the state
$\ket{0}$ conditioned on the measurement on $q_{0}$. In the two-qubit example
of Fig. 2a, we have all the information to calculate the probability of the
sampled state $p_{\ket{01}}$; moreover, we know the probability
$p_{\ket{11}}=p_{1|1}p_{0}$ which depends on the same sequence of projections
on the previous qubit, or qubits in the general case. We repeat this procedure
until the end of the qubit chain, obtaining the probability of the state
$\alpha$ as already pointed out in Eq. (4). For a detailed description of how
these quantities are efficiently computed via reduced density matrices in an
isometric tensor network, refer to App. .3.

##  III Optimal tensor network sampling

Figure 2:  _Projective measurements per qubit versus OPES sampling for one
sample and two qubits._ We represent with $|?\rangle$ the unknown state a
qubit, with $p_{i}$ the probability of measuring $\ket{i}$ on qubit $0$ and
with $p_{0|1}$ the probability of measuring $\ket{0}$ on qubit $1$ conditioned
by the measurement on qubit $0$. a) We show the standard sampling, where we
determine the measurement results by sampling a random number for each qubit.
The resulting information is just the probability of the final state. We use
the random numbers top-down, first to decide on the state of the first qubit
following the blue arrow to the left; then, the random number for the second
qubit samples state 0 and we follow the blue arrow to the right to the
complete state $\ket{10}$. The reduced density matrix for the second qubits is
calculated after using the projector $\ket{1}\bra{1}_{1}$ on the first qubit
in this example. b) We report OPES, where a single random number is sampled
for all the qubits. We collect knowledge of the probability intervals starting
from an empty bitstring and the interval $(0,1)$ containing the random number
$u_{0}$, see step _1)_. In step _2)_ , we obtain the intervals $\delta_{0}$
and $\delta_{1}$ based on the reduced density matrix of the first qubit. We
apply the projector $\ket{1}\bra{1}_{1}$ because $u_{0}\in\delta_{1}$ as
highlighted in blue. We continue with the second qubit. By the end, we know
the probability intervals $\delta_{10},\;\delta_{1}$ with
$u_{0}\in\delta_{10}\in\delta_{1}$ while generating the bitstring, highlighted
in blue; further, we know of $\delta_{0}$ and $\delta_{11}$. We enumerate the
steps as in the main text, showing the evolution of the state probability
$p_{\alpha^{*}}$ and its left boundary $b^{l}_{\alpha^{*}}$.

Optimal tensor network sampling optimizes the sampling procedure measuring
each state only once, regardless of its probability. We divide the probability
interval $[0,1]$ in unknown but ordered intervals
$\\{\delta_{\alpha}\\}_{\alpha=0,\dots,2^{n}-1}$ following the idea of the
cumulative probability distribution shown in Fig. 1c. Considering one
measurement as presented in Fig. 2b, we highlight the two intervals
$\delta_{10}$ and $\delta_{11}$ gained for the full state as well as the
interval of the partial state $\delta_{0}$. We recall that the index $\alpha$
is the integer representation of the bitstring number. After generating the
sample $\left(\delta_{\alpha},\alpha\right)$ with OPES, we remove
$\delta_{\alpha}$ from the measurable intervals, and in the next iteration, we
can only measure a new state $\beta\neq\alpha$. This exclusion of previously
sampled states enables us to explore exponentially decaying tails efficiently.

The pedagogical and simplified algorithm for measuring simultaneous
observables for two-level systems, depicted in Fig. 2b, is as follows:

  1. 1.

Set the state $\alpha^{*}$ to the empty bitstring, i.e., meaning no
measurement has taken place so far. Therefore, the interval to draw uniform
random numbers from is $\Delta=[0,1]$; set the initial left boundary of the
probability interval $b^{l}_{\alpha^{*}}$ and the state probability
$p_{\alpha^{*}}$ to

$\displaystyle b^{l}_{\alpha^{*}}=0,\quad p_{\alpha^{*}}=1\,.$ (7)

Sample $u_{0}\sim U(\Delta)$, where $U(\Delta)$ is a uniform distribution over
a given interval.

  2. 2.

Compute the probability $p_{0}=p(q_{0}=\ket{0})$ of measuring $\ket{0}$ on
site $0$ and project the wavefunction $\ket{\psi}$ on a given state according
to its probability and the random number $u_{0}$:

$\displaystyle q_{0}$ $\displaystyle=\begin{cases}\ket{0}&\text{if
}u_{0}<p_{0},\\\ \ket{1}&\text{if }u_{0}\geq p_{0},\end{cases}$ (8)
$\displaystyle\ket{\alpha}$ $\displaystyle=\ket{\alpha^{*},q_{0}}\,.$ (9)

Update the left boundary and the state probability according to the measure:

$\displaystyle b^{l}_{\alpha}$
$\displaystyle=\begin{cases}b^{l}_{\alpha^{*}}&\text{if }u_{0}<p_{0},\\\
b^{l}_{\alpha^{*}}+p_{0}&\text{if }u_{0}\geq p_{0}\end{cases}$ (10)
$\displaystyle p_{\alpha}$ $\displaystyle=p_{\alpha^{*}}\cdot p(q_{0}).$ (11)

  3. 3.

Moving to the next qubit $i$, set $\alpha^{*}=\alpha$; the bitstring
$\alpha^{*}$ now contains the qubits up to the $(i-1)^{\mathrm{th}}$ qubit.
Compute the conditional probability
$p_{0|\alpha^{*}}=p(q_{i}=\ket{0}|\alpha^{*})$. Select the measured state
according to $u$ as in the previous point and set $b^{l}_{\alpha},p_{\alpha}$
accordingly, i.e.:

$\displaystyle q_{i}$ $\displaystyle=\begin{cases}\ket{0}&\text{if
}u_{0}<b^{l}_{\alpha^{*}}+p_{0|\alpha^{*}},\\\ \ket{1}&\text{if }u_{0}\geq
b^{l}_{\alpha^{*}}+p_{0|\alpha^{*}},\end{cases}$ (12)
$\displaystyle\ket{\alpha}$ $\displaystyle=\ket{\alpha^{*},q_{i}}\,,$ (13)
$\displaystyle b^{l}_{\alpha}$
$\displaystyle=\begin{cases}b^{l}_{\alpha^{*}}&\text{if
}u_{0}<b^{l}_{\alpha^{*}}+p_{0|\alpha^{*}}\,,\\\
b^{l}_{\alpha^{*}}+p_{0|\alpha^{*}}p_{\alpha^{*}}&\text{if }u_{0}\geq
b^{l}_{\alpha^{*}}+p_{0|\alpha^{*}}\,,\end{cases}$ (14) $\displaystyle
p_{\alpha}$ $\displaystyle=p_{\alpha^{*}}\cdot p(q_{i}|\alpha^{*}).$ (15)

  4. 4.

Repeat step _3._ until all the sites of the system are measured.
Further details are shown in Figure 2b: we have access to the information of
many nested intervals after measuring all the sites. For example, if the final
result is the interval $\delta_{10}$ we have information about the following
intervals:

$\displaystyle\delta_{01}\in\delta_{1}\in[0,1],\;\delta_{0},\delta_{11}.$ (16)

We denote $\delta_{i}$ the interval where only the first qubit has been
measured to simplify the notation. We stress that knowing $\delta_{0}$ is
enough to compute $b^{l}_{\alpha^{*}}$, and we do not need to know
$\delta_{00},\delta_{01}$. In general, the obtained results are:

$\displaystyle b^{l}_{\alpha}$
$\displaystyle=\sum_{i}p\left(\overline{q_{i}}|\\{q_{j}\\}_{j=0}^{i-1}\right)\theta(q_{i}-\frac{1}{2}),$
(17) $\displaystyle p_{\alpha}$ $\displaystyle=p(q_{n-1}|q_{n-2}\dots
q_{0})\dots p(q_{1}|q_{0})p(q_{0}),$ (18)

where $\theta(x)$ is the Heavyside step function and
$\overline{q_{i}}=|1-q_{i}|.$

Furthermore, we highlight that at each iteration we actually compute two
probability intervals, i.e. $\delta_{q_{0}q_{1}\dots q_{n-2}q_{n-1}}$ and
$\delta_{q_{0}q_{1}\dots q_{n-2}\overline{q_{n-1}}}$.

  5. 5.

Remove the measured probability interval from our sampling interval, i.e.
$\Delta=\Delta\setminus\delta_{\alpha}$. Each new sampled number $u_{j}~{}\sim
U(\Delta)$ explores a different probability interval
$\delta_{\beta\neq\alpha}$.

This presentation of the algorithm is pedagogical for understanding the
procedure and generates one sample. Steps 1) to 5) represent in the complete
algorithm only a single superiteration with a number of random numbers sampled
$n_{\rm r}=1$; we report the complete algorithm in Appendix .2. The main
differences are that for efficient computation of the intervals
$\delta_{\alpha}$, multiple random numbers $\\{u_{j}\\}_{j=1,\dots,n_{r}}$ are
drawn for each superiteration, repeating then $n_{\rm\scriptscriptstyle S}$
superiterations. While we do not obtain anymore one new state $\alpha$ for
each random number sampled from $U(\Delta)$, we are able to cache intermediate
results to optimize the procedure; nonetheless, we are still avoiding
additional computation for the same interval. For more information about the
caching strategy using the information of the intervals, see Appendix .3.

As a rule of thumb, if you expect your state to have a peaked distribution
with small tails you should run the first superiteration with $n_{\rm r}>1$,
and then many iterations with a bigger $n_{\rm r}\gg 1$ to explore the space.
Instead, it is much better to always use $n_{\rm r}\gg 1$ if you expect a flat
distribution, since the overhead due to the cached intermediate results might
reduce the performances.

##  IV Benchmarking sampling on applications

In this section, we consider a variety of quantum states, e.g., quantum
images, quantum Ising ground states, Rydberg atoms, and quantum magic; we
compare the computational time and accuracy of the sampling methods. Our
results showcase the strengths of OPES over the standard sampling for states
with decaying tails in the probability distribution. All computational times
for the benchmarking of the sampling are referring to Cineca’s _Galileo100_
cluster, using a single compute node and 32 cores; the node consists of Intel
CascadeLake 8260 with 2xCPU x86 Intel Xeon Platinum 8276-8276L (2.4Ghz).

###  IV.1 Quantum state with Gaussian probability distribution with matrix
product states

First, we evaluate the performance of OPES sampling on a Gaussian quantum
state. The probability density of the state $\ket{\psi}$ has a Gaussian
profile

$\displaystyle\ket{\psi}=\frac{1}{\mathcal{N}}\sum_{i=0}^{2^{n}-1}\sqrt{e^{-\frac{(x_{i})^{2}}{2\sigma^{2}}}}\ket{i},\;x_{i}=\frac{i-\overline{i}}{2^{n}-1},$
(19)

where the state is a superposition of basis states $\ket{i}$ that one samples.
We generate the state exactly and then convert it to an MPS representation
without any truncation. The parameter $n$ sets the dimension of the Hilbert
space in terms of qubits; the standard deviation $\sigma$ defines how close
the state is to a delta-peak and the value $\bar{i}$ shifts the maximum of the
Gaussian between states. The scalar $\mathcal{N}$ ensures that the state
$\ket{\psi}$ is normalized. Then, we compare the computational time $t$
required by two different methods to achieve a certain level of accuracy in
characterizing the state distribution. Specifically, we aim to know the
state’s distribution up to a given threshold $\epsilon$. Therefore, we define
the coverage $\mathcal{C}$ and threshold $\epsilon$ of a set of states with
$\\{\delta_{1},\ldots,\delta_{m}\\}$ as

$\displaystyle\mathcal{C}$
$\displaystyle=\sum_{\alpha=1}^{m}\left|c_{\alpha}\right|^{2}\,,\qquad\epsilon=1-\mathcal{C}\,.$
(20)

We show the results for ten qubits $n=10$, $\overline{i}=2^{n-1}-1$, and
different values of the variance $\sigma^{2}$ in Fig. 3. The computational
time for characterizing the probability interval up to the threshold
$\epsilon$ increases for both methods as the variance increases. However, we
find that OPES sampling outperforms traditional sampling significantly due to
its efficient use of cached information and the tail’s sampling. Let us define
the speedup $\mathcal{S}$

$\displaystyle\mathcal{S}=\frac{t_{\text{standard}}}{t_{\text{OPES}}},$ (21)

with $t_{\text{OPES}}$ and $t_{\text{standard}}$ the computational time for
the OPES and standard sampling, respectively. In Fig. 3b, we show that the
speedup achieved scales with the threshold $\epsilon$, increasing from one
order of magnitude for $\epsilon=10^{-2}$ to three orders of magnitude for
$\epsilon=10^{-4}$. An in-depth argument on the scaling and speedup can be
found in App. .1, but we have no expected theoretical scaling for a comparison
in case of Eq. (19). Finally, the computational time $t$ for OPES sampling is
almost independent of the chosen coverage threshold. These results demonstrate
the superiority of OPES sampling over traditional sampling when one aims to
cover most of the probability distribution.

\begin{overpic}[width=411.93767pt,unit=1mm]{picts/gaussian_comparison.pdf}
\put(-3.0,60.0){a)}
\end{overpic}\begin{overpic}[width=411.93767pt,unit=1mm]{picts/gaussian_speedup.pdf}
\put(-3.0,60.0){b)} \end{overpic}

Figure 3: _Quantum state with a Gaussian probability distribution._ a) The
computational time $t$ of the standard sampling versus OPES for different
probability coverage thresholds $\epsilon$ with $n=10$ qubits. States with a
higher standard deviation $\sigma$ are more difficult to sample due to the
tails of the distribution. b) Computational speedup $\mathcal{S}$ of the OPES
sampling over the standard one with respect to $\epsilon$ at fixed variance
$\sigma^{2}$.

###  IV.2 Quantum image with matrix product states

We proceed and analyze a state obtained after the application of a quantum
circuit. In particular, we are not interested in exploring final states of
quantum algorithms, since these algorithms are built to be easy to sample to
minimize the number of measurements; examples of such algorithms are the
Grover or Shor algorithms Grover (1996); Shor (1994). We instead investigate a
state obtained "at runtime". To perform operations on images using quantum
hardware, we first have to load the image. There are several possible ways of
doing so: amplitude encoding, angle encoding, or flexible representation of
quantum images (FRQI) Le  _et al._ (2011). We focus on the latter case,
encoding images using the algorithm in Ref. Amankwah  _et al._ (2022). The
question at hand is if we can reconstruct the image with high accuracy after
loading it, using tensor network sampling. This procedure helps in
understanding if the image is correctly reproduced even for sizes that are not
accessible with the exact state vector simulation. In particular, we choose to
encode the quantum tea logo, shown in Fig. 4a. The logo is a $600\times 600$
grey-scale image with a total of $M=360000$ pixels encoded in a $20$-qubit
state. We now proceed to a quantitative study of the sampling. In Fig. 4b, we
show the computational time $t$ needed to reach a given mean square error MSE
with the original image, defined as:

$\displaystyle\text{MSE}(t)=\frac{1}{M}\sum_{i=1}^{M}\Big{(}\text{P}_{sv}^{i}-\text{P}_{sm}^{i}(t)\Big{)}^{2},$
(22)

where $\text{P}_{sv}$ is the picture obtained from the complete state vector
of the state and $\text{P}_{sm}(t)$ is the sampled picture at time $t$. The
index $i$ runs over the pixels of the image. Thus, the MSE is improved by over
two orders of magnitude in the figure of merit while saving one order of
magnitude in computational time. We also show the logos retrieved from the
quantum state at different superiterations of the OPES sampling in Fig. 4a.
The FRQI algorithm benefits in our opinion the most from the OPES sampling out
of our examples. As we see in the inline of Fig. 4b, the speedup $\mathcal{S}$
of OPES over the standard method is exponential in the MSE, which is agreement
with App. .1.

\begin{overpic}[width=411.93767pt,unit=1mm]{picts/frqi_qlogos.pdf}
\put(-5.0,33.0){a)}
\end{overpic}\begin{overpic}[width=411.93767pt,unit=1mm]{picts/frqi_image_mse.pdf}
\put(-3.0,60.0){b)} \end{overpic}

Figure 4: _Flexible representation of quantum images (FRQI) algorithm._ a)
Quantum tea logos used in the flexible representation of quantum images
algorithm. From left to right, we show the original image, the reconstructed
image after three superiterations of the OPES sampling, the image at the end
of the OPES sampling. b) Computational time $t$ of the standard sampling
versus OPES for reaching a given mean square error MSE in the image
reconstruction. In the inset, the speedup of OPES over the standard method is
highlighted.

###  IV.3 Quantum Ising model with tree tensor networks

As the next example, we consider the quantum Ising model Sachdev (2011) as an
example from condensed matter. We use the definition of the quantum Ising
model as

$\displaystyle H$
$\displaystyle=J\left(-\sum_{i=1}^{n-1}\sigma_{i}^{z}\sigma_{i+1}^{z}-g\sum_{i=1}^{n}\sigma_{i}^{x}\right)\,,$
(23)

where the parameter $J$ set the energy scale and $g$ tunes the external field,
respectively. The Pauli matrices acting on site $i$ are $\sigma_{i}^{\alpha}$.
From the two limits, we expect three regions with respect to the sampling: (i)
close to the ferromagnetic limit, the states almost align with our measurement
basis; (ii) close to the paramagnetic limit, the local product states are
orthogonal to the local measurement basis; and (iii) the region around the
quantum critical point with high entanglement. The latter two are thus
potentially difficult to sample, especially in the case of the paramagnetic
limit where the probability of each state scales inversely with the dimension
of the Hilbert space.

The sampling of a fixed number of samples $n_{\rm r}$ within one
superiteration scales better either for the standard or OPES approach
according to Fig. 5a, depending on the quantum phase represented by one point
in the ferromagnetic and paramagnetic phase, see $g=0.5$ and $g=1.5$
respectively, and an additional point close to the quantum critical point with
$g=1$. We consider a system size of $n=64$. The probability distribution for
$g=1$ and $g=1.5$ are too similar to a uniform distribution and OPES does not
show a benefit here, but is even slower due to a potential overhead necessary
for caching intermediate results. If we aim to obtain a certain coverage
$\mathcal{C}$ of the probability space instead of a fixed number of
iterations, the task of reaching the coverage becomes computationally more
expensive, especially for flat-distributions. For example, we have to sample
in an exponentially large space in the paramagnetic limit

\begin{overpic}[width=368.57964pt,unit=1mm]{./picts/ising_sampling_fixedns_cmdarg_e.pdf}
\put(-3.0,60.0){a)}
\end{overpic}\begin{overpic}[width=368.57964pt,unit=1mm]{./picts/ising_gfields_coverage_cmdarg_f.pdf}
\put(-3.0,55.0){b)} \end{overpic}

\begin{overpic}[width=368.57964pt,unit=1mm]{./picts/ising_thermal_sampling_cmdarg_g.pdf}
\put(-3.0,60.0){c)} \end{overpic}

Figure 5: _Comparisons for the quantum Ising model and tree tensor networks._
We calculate the speedup $\mathcal{S}$ of OPES as the ratio of the
computational time of the standard method over OPES. a) We sample $n_{\rm r}$
states within one superiteration from the ground state at the critical point
$g=g_{c}=1.0$ for $64$ qubits. b) The speedup $\mathcal{S}$ to reach a
coverage $\mathcal{C}$ for 16 qubits across the phase diagram. The
paramagnetic phase is difficult as the eigenstate does not align with the
measurement basis. c) The speedup $\mathcal{S}$ to produce $n_{\rm r}$ samples
within one superiteration for the density matrix at finite temperatures
$T\in[1.1,0.49,0.01]$ and two external fields $g$.

because the local eigenstates are orthogonal to our measurement basis. Keeping
the measurement basis fixed throughout helps to understand the scaling for
other models where a change into the "good" measurement basis might not be
obvious. But we reduce the number of qubits to $N=16$ and aim for coverages
$\mathcal{C}=30\%$ and $\mathcal{C}=50\%$ in Fig. 5b; without recucing the
number of qubits $N$, we need to sample $2^{N-1}$ states in the paramagnetic
limit for $N=64$ to reach $\mathcal{C}=50\%$. The better scaling of the OPES
method for all fields $g$ originates from the fact that a new iteration does
not resample states from previous iterations. In contrast, the standard method
has a probability of almost 50% to throw a state away because the generated
state is already known when reaching $50\%$ coverage. Thus, the OPES still
pays off in flat distributions for high coverages $\mathcal{C}$, though this
scenario is restricted to moderate number of qubits due to the exponential
growth of the Hilbert space.

The third example moves from pure states, which we consider throughout all the
other examples, to density matrices. We generate finite temperature density
matrices for the external fields $g=0.5$ and $g=1.0$ via an imaginary time-
dependent variational principle (TDVP) time evolution Haegeman  _et al._
(2016) with a locally purified tensor network (LPTN) Werner  _et al._ (2016),
which is then converted into a tree tensor operator (TTO) Reinić (2022). The
behavior for the external field resembles the ground state, i.e., OPES shows a
greater benefit inside the ferromagnetic phase for the data points with
$g=0.5$. Surprisingly, the data indicates a bigger speedup at higher
temperatures although the high-temperature limit itself is the completely-
mixed state with a flat distribution. In summary, any problem where a lower
limit on the unsampled space is required can profit significantly from OPES
even in scenarios similar to the condensed matter case across all phases
treated here.

###  IV.4 Crosstalk errors in Rydberg atoms with tree tensor networks

An error analysis is useful for the benchmarking of quantum processing units,
e.g., for the state preparation. In such a scenario, we expect to be almost
perfectly in the target state, where the remaining probability represents the
error. Depending on the fidelity of the preparation, the measurement of a
state with error has a small probability. The preparation of the Greenberger-
Horne-Zeilinger state (GHZ) with errors including the projective measurement
is shown in Ref. Jaschke  _et al._ (2022) and we consider here the aspect of
the computational scaling of analyzing such a state. Figure 6a shows the
computational time $t$ for an imperfect 64-qubit GHZ state, where errors have
been introduced during the protocol due to crosstalk; while the computation
time for projective measurements scales with the number of samples for the
standard approach, the OPES approach has an almost constant computational time
$t$. This almost constant time $t$ is related to the two large probabilities
for the states $\ket{00\ldots 00}$ and $\ket{11\ldots 11}$. This scenario is
one of the most favorable for OPES.

\begin{overpic}[width=390.25534pt,unit=1mm]{./picts/rydberg_ttn_a_sampling_cmdarg_a.pdf}
\put(-3.0,60.0){a)}
\end{overpic}\begin{overpic}[width=390.25534pt,unit=1mm]{./picts/rydberg_ttn_a_coverage_cmdarg_b.pdf}
\put(-3.0,55.0){b)} \end{overpic}

Figure 6: _Comparison for sampling an imperfect GHZ prepared under crosstalk._
The state is obtained from a simulation of Rydberg atoms encountering
crosstalk during the preparation, which leads to deviations from the GHZ
state. a) We sample a fixed number of samples $n_{\rm r}$ within one
superiteration, where the major probability in the states $\ket{00\ldots 00}$
and $\ket{11\ldots 11}$ leads to the better scaling of OPES. b) We aim for a
coverage $\mathcal{C}\geq 99.8\%$ and show the coverage $\mathcal{C}$ during
the superiterations of the algorithms.

We further extract the coverage $\mathcal{C}$ as defined in Eq. (20) from the
quantum state for each sampling method, i.e., how much of the probability
space has not been sampled. This data could be used to generate data according
to a quantum state’s probability distribution up to a certain accuracy even
without a tensor network representation. Figure 6b shows how the coverage
improves over time with the superiterations of the algorithm, while we sample
$n_{\rm r}=250$ states per superiteration for the default approach in
comparison to $n_{\rm r}=25$ states for OPES. As OPES can target exclusively
previously unexplored probability intervals, we observe a significant
advantage.

###  IV.5 Magic with matrix product states

To stress the generality of the method mentioned in Sec. I, we now generalize
the sampling procedure to measure the quantum magic of a quantum state, also
known as non-stabilizerness. Quantum magic is an important quantity since it
limits the classical simulation of quantum systems for stabilizer simulators
Howard  _et al._ (2014); Seddon  _et al._ (2021), similarly to how the
entanglement limits the simulations with tensor networks. Quantum magic can be
measured through Stabilizer Rényi Entropies (SREs) Leone  _et al._ (2022).
Reference Lami and Collura (2023) introduces a technique to measure the SREs
through standard sampling of an MPS state. Thus, applying OPES can speed up
the process of measuring quantum magic. The magic of a $n$-qubits state using
the $m$-order SREs is defined as:

$\displaystyle\mathcal{M}_{m}=\frac{1}{1-m}\log\left(\frac{1}{\mathcal{N}}\sum_{\mu=0}^{\mathcal{N}}\Pi(\vec{\sigma})^{m-1}\right)-n\log
2,$ (24)

where $\Pi(\vec{\sigma})$ is the probability of measuring a given Pauli
string, $\sigma_{i}\in\\{\mathbb{I},X,Y,Z\\}$. To apply OPES to measure
$\mathcal{M}$, we target the $\Pi(\vec{\sigma})$, keeping into account that
each measure of probability on the single site $\pi(\sigma_{i})$ branches into
four possibilities ($\sigma_{i}\in\\{\mathbb{I},X,Y,Z\\}$) instead of the two
possibilities ($b_{i}=\\{\ket{0},\ket{1}\\}$) treated up to now. For an
extensive explanation of how to measure $\pi(\sigma_{i})$ efficiently with an
MPS, please refer to Ref. Lami and Collura (2023). As an example, we compute
the order-1 SREs of the state

$\displaystyle\ket{\psi(\phi)}=\left(\frac{\ket{0}+e^{i\phi}\ket{1}}{\sqrt{2}}\right)^{\otimes
n}.$ (25)

Notice that $\ket{\psi(0)}$ for $\phi=0$ is a stabilizer state, and thus has
zero magic. The magic of the state is instead maximum for
$\phi=\frac{\pi}{4}$. We call the exact 1-SREs of this state
$\widetilde{\mathcal{M}}$:

$\displaystyle\widetilde{\mathcal{M}}=$
$\displaystyle-\cos^{2}\phi\log\left(\left|\cos\phi\right|\right)$
$\displaystyle-\sin^{2}\phi\log\left(\left|\sin\phi\right|\right).$ (26)

In Fig. 7, we address the computation of $\mathcal{M}_{1}$ for $n=10$ qubits
using for OPES $n_{S}=10$ superiterations with $n_{r}=100$ samples each. The
results are mediated over $10$ realization of the computations. We use for the
standard method the same number of samples. We notice that, due to the
flatness of the distribution, OPES gives no relevant improvement when
$\phi=\frac{\pi}{64}\sim 0$. This result is to be expected, as already
discussed for the paramagnetic phase of the quantum Ising model in Sec. IV.3.
However, we can observe both an improvement in the estimation of the magic and
a speedup in the other two cases, i.e. $\phi=\frac{\pi}{8}$ and
$\phi=\frac{\pi}{4}$. Optimal tensor network sampling is twice as fast as the
standard method, and the estimation of the quantum magic converges much faster
with respect to the standard method.

Figure 7: Magic estimation through sampling. We compare the standard and OPES
methods in computing the magic of the state known by its theory formula, see
$\widetilde{\mathcal{M}}$ in Eq. (26), for $n=10$ qubits and different
parameters $\phi$ with the same number of samples. The points are obtained as
an average over ten realizations. Optimal tensor network sampling is always at
least twice as fast as the standard method and converges faster to the true
value $\mathcal{M}_{1}$ for $\phi=\frac{\pi}{8}$ and $\phi=\frac{\pi}{4}$.

##  V Conclusions

Overall, we have introduced the optimal tensor network sampling (OPES), a
sampling method for tensor network states that avoids sampling known intervals
and caches information for optimal performance. From the computational
statistics point of view, the algorithm represents an extension to the inverse
transform sampling, which is used to sample pseudo-random numbers from an
arbitrary probability distribution function by knowing its inverse probability
distribution function. The extension consists in the fact that we never fully
know the cumulative distribution function or its inverse of the quantum state
represented as a tensor network. The method is applicable to all isometric
tensor network states as well as exact diagonalization. We have investigated
several quantum states and compared the performance of the standard and OPES
methods in characterizing them; these results are obtained for matrix product
states (MPS), tree tensor networks (TTN), and tree tensor operators (TTO).

In detail, we have shown that OPES offers a significant speedup compared to
the standard method, especially when covering most of the probability
distribution is demanded. On the one hand, this speedup is due to the
efficient caching of previous information employing the MPS or TTN; on the
other hand, we work with the underlying probability intervals of the full
Hilbert space which allow us to exclude previously explored intervals from
future superiterations. For example, OPES leads to a speedup of more than an
order of magnitude while achieving lower errors versus the standard sampling
for applications as the flexible representation of quantum images (FRQI)
algorithm. We have shown how the approach is generalizable to the measure of
magic with matrix product state, and its speedup over standard methods. The
speedup in computation time also leads to a reduced integrated energy
consumption of the CPU as our additional data in Ballarin  _et al._ (2024b)
shows for one example.

Future directions of this approach for sampling can go into the direction of
comparing parallelized versions, optimizing the caching strategy, or exploring
the potential benefit in other algorithms and applications. While the standard
sampling trivially parallelizes via data parallelism, OPES can, for example,
parallelize over the groups of the sorted random numbers. For the application
side, there exists a plethora of candidates with potential speedup where the
examples given within this manuscript give a first intuition which
applications benefit most from OPES, i.e., quantum states with a decaying
probability distribution. Furthermore, since OPES only needs the computation
of expectation values of local operators and marginal probabilities it could
be applied to any isometric tensor network such as isoPEPS Zaletel and
Pollmann (2020). Similarly, the strings of projective measurements can be
naturally replaced with Positive Operator-Valued Measures (POVM).

_Code and data availability $-$_ The code implementing the algorithm for both
MPS and TTNs is distributed through the quantum tea leaves Ballarin  _et al._
(2024a) python package of _Quantum TEA_ , while the code and data to generate
the plots of this work is available via zenodo Ballarin  _et al._ (2024c, d).
All figures are available together with metadata under a CC-BY license at
Ballarin  _et al._ (2024b).

_Acknowledgments $-$_ We thank Nora Reinić for providing finite-temperature
states in the tree tensor operator representation. We thank Giovanni Cataldi,
Marcello Dalmonte, Martina Frau, Marco Rigobello, Ilaria Siloi, Poetri
Tarabunga, Marco Tesoro, and Emanuele Tirrito for discussions and feedback. We
acknowledge financial support from the Italian Ministry of University and
Research (MUR) via PRIN2022 project TANQU, and the Departments of Excellence
grant 2023-2027 Quantum Frontiers; from the European Union via H2020 projects
EuRyQa, and TEXTAROSSA, the QuantERA projects QuantHEP and T-NISQ, and the
Quantum Flagship project Pasquans2, from the German Federal Ministry of
Education and Research (BMBF) via the funding program quantum technologies $-$
from basic research to market $-$ project QRydDemo, and from the World Class
Research Infrastructure $-$ Quantum Computing and Simulation Center (QCSC) of
Padova University. This work was performed in part at the Aspen Center for
Physics, which is supported by National Science Foundation grant PHY-2210452;
the participation of D. J. at the Aspen Center for Physics was supported by
the Simons Foundation. We also acknowledge computation time supported support
by the state of Baden-Württemberg through bwHPC and the German Research
Foundation (DFG) through grant no INST 40/575-1 FUGG (JUSTUS 2 cluster), the
Atos Dibona cluster via TEXTAROSSA, as well as computation time on Cineca’s
_Galileo100_ machine.

## References

  * Monroe (2002) C. Monroe, “Quantum information processing with atoms and photons,” Nature 416, 238–246 (2002).
  * Monroe and Kim (2013) C. Monroe and J. Kim, “Scaling the ion trap quantum processor,” Science 339, 1164–1169 (2013).
  * Saffman (2018) Mark Saffman, “Quantum computing with neutral atoms,” National Science Review 6, 24–25 (2018).
  * Wendin (2017) G. Wendin, “Quantum information processing with superconducting circuits: a review,” Reports on Progress in Physics 80, 106001 (2017).
  * Bruzewicz  _et al._ (2019) Colin D. Bruzewicz, John Chiaverini, Robert McConnell, and Jeremy M. Sage, “Trapped-ion quantum computing: Progress and challenges,” Applied Physics Reviews 6, 021314 (2019).
  * Adams  _et al._ (2019) C. S. Adams, J. D. Pritchard, and J. P. Shaffer, “Rydberg atom quantum technologies,” Journal of Physics B: Atomic, Molecular and Optical Physics 53, 012002 (2019).
  * Kaushal  _et al._ (2020) V. Kaushal, B. Lekitsch, A. Stahl, J. Hilder, D. Pijn, C. Schmiegelow, A. Bermudez, M. Müller, F. Schmidt-Kaler, and U. Poschinger, “Shuttling-based trapped-ion quantum information processing,” AVS Quantum Science 2, 014101 (2020).
  * Kjaergaard  _et al._ (2020) Morten Kjaergaard, Mollie E. Schwartz, Jochen Braumüller, Philip Krantz, Joel I.-J. Wang, Simon Gustavsson, and William D. Oliver, “Superconducting qubits: Current state of play,” Annual Review of Condensed Matter Physics 11, 369–395 (2020).
  * Morgado and Whitlock (2021) M. Morgado and S. Whitlock, “Quantum simulation and computing with Rydberg-interacting qubits,” AVS Quantum Science 3, 023501 (2021).
  * Nagy and Savona (2019) Alexandra Nagy and Vincenzo Savona, “Variational quantum Monte Carlo method with a neural-network ansatz for open quantum systems,” Physical review letters 122, 250501 (2019).
  * Qiskit contributors (2023) Qiskit contributors, “Qiskit: An open-source framework for quantum computing,”  (2023).
  * The cuQuantum development team (2023) The cuQuantum development team, “NVIDIA cuQuantum SDK,”  (2023).
  * Efthymiou  _et al._ (2022) Stavros Efthymiou, Marco Lazzarin, Andrea Pasquale, and Stefano Carrazza, “Quantum simulation with just-in-time compilation,” Quantum 6, 814 (2022).
  * Jaschke  _et al._ (2022) Daniel Jaschke, Alice Pagano, Sebastian Weber, and Simone Montangero, “Ab-initio two-dimensional digital twin for quantum computer benchmarking,” arXiv preprint 2210.03763 (2022).
  * Arute  _et al._ (2019) Frank Arute, Kunal Arya, Ryan Babbush, Dave Bacon, Joseph C. Bardin, Rami Barends, Rupak Biswas, Sergio Boixo, Fernando G. S. L. Brandao, David A. Buell, Brian Burkett, Yu Chen, Zijun Chen, Ben Chiaro, Roberto Collins, William Courtney, Andrew Dunsworth, Edward Farhi, Brooks Foxen, Austin Fowler, Craig Gidney, Marissa Giustina, Rob Graff, Keith Guerin, Steve Habegger, Matthew P. Harrigan, Michael J. Hartmann, Alan Ho, Markus Hoffmann, Trent Huang, Travis S. Humble, Sergei V. Isakov, Evan Jeffrey, Zhang Jiang, Dvir Kafri, Kostyantyn Kechedzhi, Julian Kelly, Paul V. Klimov, Sergey Knysh, Alexander Korotkov, Fedor Kostritsa, David Landhuis, Mike Lindmark, Erik Lucero, Dmitry Lyakh, Salvatore Mandrà, Jarrod R. McClean, Matthew McEwen, Anthony Megrant, Xiao Mi, Kristel Michielsen, Masoud Mohseni, Josh Mutus, Ofer Naaman, Matthew Neeley, Charles Neill, Murphy Yuezhen Niu, Eric Ostby, Andre Petukhov, John C. Platt, Chris Quintana, Eleanor G. Rieffel, Pedram Roushan, Nicholas C. Rubin, Daniel Sank, Kevin J. Satzinger, Vadim Smelyanskiy, Kevin J. Sung, Matthew D. Trevithick, Amit Vainsencher, Benjamin Villalonga, Theodore White, Z. Jamie Yao, Ping Yeh, Adam Zalcman, Hartmut Neven, and John M. Martinis, “Quantum supremacy using a programmable superconducting processor,” Nature 574, 505–510 (2019).
  * Zhong  _et al._ (2020) Han-Sen Zhong, Hui Wang, Yu-Hao Deng, Ming-Cheng Chen, Li-Chao Peng, Yi-Han Luo, Jian Qin, Dian Wu, Xing Ding, Yi Hu, Peng Hu, Xiao-Yan Yang, Wei-Jun Zhang, Hao Li, Yuxuan Li, Xiao Jiang, Lin Gan, Guangwen Yang, Lixing You, Zhen Wang, Li Li, Nai-Le Liu, Chao-Yang Lu, and Jian-Wei Pan, “Quantum computational advantage using photons,” Science 370, 1460–1463 (2020).
  * Zhou  _et al._ (2020) Yiqing Zhou, E. Miles Stoudenmire, and Xavier Waintal, “What limits the simulation of quantum computers?” Phys. Rev. X 10, 041038 (2020).
  * Liu  _et al._ (2021) Yong (Alexander) Liu, Xin (Lucy) Liu, Fang (Nancy) Li, Haohuan Fu, Yuling Yang, Jiawei Song, Pengpeng Zhao, Zhen Wang, Dajia Peng, Huarong Chen, Chu Guo, Heliang Huang, Wenzhao Wu, and Dexun Chen, “Closing the "quantum supremacy" gap: Achieving real-time simulation of a random quantum circuit using a new sunway supercomputer,” in  _Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis_, SC ’21 (Association for Computing Machinery, New York, NY, USA, 2021).
  * Bulmer  _et al._ (2022) Jacob F. F. Bulmer, Bryn A. Bell, Rachel S. Chadwick, Alex E. Jones, Diana Moise, Alessandro Rigazzi, Jan Thorbecke, Utz-Uwe Haus, Thomas Van Vaerenbergh, Raj B. Patel, Ian A. Walmsley, and Anthony Laing, “The boundary for quantum advantage in Gaussian boson sampling,” Science Advances 8, eabl9236 (2022).
  * Daley  _et al._ (2022) Andrew J. Daley, Immanuel Bloch, Christian Kokail, Stuart Flannigan, Natalie Pearson, Matthias Troyer, and Peter Zoller, “Practical quantum advantage in quantum simulation,” Nature 607, 667–676 (2022).
  * Oh  _et al._ (2023) Changhun Oh, Minzhao Liu, Yuri Alexeev, Bill Fefferman, and Liang Jiang, “Tensor network algorithm for simulating experimental Gaussian boson sampling,” arXiv preprint arXiv:2306.03709 (2023).
  * White (1992) Steven R. White, “Density matrix formulation for quantum renormalization groups,” Phys. Rev. Lett. 69, 2863–2866 (1992).
  * Rommer and Östlund (1997) Stefan Rommer and Stellan Östlund, “Class of ansatz wave functions for one-dimensional spin systems and their relation to the density matrix renormalization group,” Phys. Rev. B 55, 2164–2181 (1997).
  * Rommer and Östlund (1999) Stefan Rommer and Stellan Östlund, “Thermodynamic limit and matrix-product states,” in  _Density-Matrix Renormalization_ , edited by Ingo Peschel, Matthias Kaulke, Xiaoqun Wang, and Karen Hallberg (Springer Berlin Heidelberg, Berlin, Heidelberg, 1999) pp. 67–89.
  * Schollwöck (2005) U. Schollwöck, “The density-matrix renormalization group,” Rev. Mod. Phys. 77, 259–315 (2005).
  * Perez-Garcia  _et al._ (2007) D. Perez-Garcia, F. Verstraete, M. M. Wolf, and J. I. Cirac, “Matrix Product State Representations,” Quantum Information and Computation 7, 401–430 (2007).
  * McCulloch (2007) Ian P. McCulloch, “From density-matrix renormalization group to matrix product states,” Journal of Statistical Mechanics: Theory and Experiment 2007, P10014 (2007).
  * Schollwöck (2011) Ulrich Schollwöck, “The density-matrix renormalization group in the age of matrix product states,” Annals of Physics 326, 96 – 192 (2011), January 2011 Special Issue.
  * Shi  _et al._ (2006) Y.-Y. Shi, L.-M. Duan, and G. Vidal, “Classical simulation of quantum many-body systems with a tree tensor network,” Phys. Rev. A 74, 022320 (2006).
  * Tagliacozzo  _et al._ (2009) L. Tagliacozzo, G. Evenbly, and G. Vidal, “Simulation of two-dimensional quantum systems using a tree tensor network that exploits the entropic area law,” Phys. Rev. B 80, 235127 (2009).
  * Murg  _et al._ (2010) V. Murg, F. Verstraete, Ö. Legeza, and R. M. Noack, “Simulating strongly correlated quantum systems with tree tensor networks,” Phys. Rev. B 82, 205105 (2010).
  * Gerster  _et al._ (2014) M. Gerster, P. Silvi, M. Rizzi, R. Fazio, T. Calarco, and S. Montangero, “Unconstrained tree tensor network: An adaptive gauge picture for enhanced performance,” Phys. Rev. B 90, 125154 (2014).
  * F.  _et al._ (2008) Verstraete F., Murg V., and Cirac J.I., “Matrix product states, projected entangled pair states, and variational renormalization group methods for quantum spin systems,” Advances in Physics 57, 143–224 (2008).
  * Orús (2014) Román Orús, “A practical introduction to tensor networks: Matrix product states and projected entangled pair states,” Annals of Physics 349, 117 – 158 (2014).
  * Zwolak and Vidal (2004) Michael Zwolak and Guifré Vidal, “Mixed-state dynamics in one-dimensional quantum lattice systems: A time-dependent superoperator renormalization algorithm,” Phys. Rev. Lett. 93, 207205 (2004).
  * 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).
  * Rizzi  _et al._ (2008) Matteo Rizzi, Simone Montangero, and Guifre Vidal, “Simulation of time evolution with multiscale entanglement renormalization ansatz,” Phys. Rev. A 77, 052328 (2008).
  * Cincio  _et al._ (2008) Lukasz Cincio, Jacek Dziarmaga, and Marek M. Rams, “Multiscale entanglement renormalization ansatz in two dimensions: Quantum Ising model,” Phys. Rev. Lett. 100, 240603 (2008).
  * McCulloch (2008) I. P. McCulloch, “Infinite size density matrix renormalization group, revisited,” arXiv preprint 0804.2509 (2008), arXiv:0804.2509 [cond-mat.str-el] .
  * Jordan  _et al._ (2008) J. Jordan, R. Orús, G. Vidal, F. Verstraete, and J. I. Cirac, “Classical simulation of infinite-size quantum lattice systems in two spatial dimensions,” Phys. Rev. Lett. 101, 250602 (2008).
  * Evenbly and Vidal (2009) G. Evenbly and G. Vidal, “Algorithms for entanglement renormalization,” Phys. Rev. B 79, 144108 (2009).
  * Evenbly and Vidal (2014) G. Evenbly and G. Vidal, “Class of highly entangled many-body states that can be efficiently simulated,” Phys. Rev. Lett. 112, 240502 (2014).
  * Phien  _et al._ (2015) Ho N. Phien, Johann A. Bengua, Hoang D. Tuan, Philippe Corboz, and Román Orús, “Infinite projected entangled pair states algorithm improved: Fast full update and gauge fixing,” Phys. Rev. B 92, 035142 (2015).
  * 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).
  * Silvi  _et al._ (2019) Pietro Silvi, Ferdinand Tschirsich, Matthias Gerster, Johannes Jünemann, Daniel Jaschke, Matteo Rizzi, and Simone Montangero, “The Tensor Networks Anthology: Simulation techniques for many-body quantum lattice systems,” SciPost Phys. Lect. Notes , 8 (2019).
  * Montangero (2018) S. Montangero,  _Introduction to Tensor Network Methods_ (Springer Nature Switzerland AG, Cham, CH, 2018).
  * Jaschke  _et al._ (2018) Daniel Jaschke, Simone Montangero, and Lincoln D. Carr, “One-dimensional many-body entangled open quantum systems with tensor network methods,” Quantum Science and Technology 4, 013001 (2018).
  * Reinić (2022) Nora Reinić,  _Tree tensor networks for quantum many-body systems at finite temperature_ , Master’s thesis, University of Zagreb (2022).
  * Arceci  _et al._ (2022) L. Arceci, P. Silvi, and S. Montangero, “Entanglement of formation of mixed many-body quantum states via tree tensor operators,” Phys. Rev. Lett. 128, 040501 (2022).
  * Bañuls (2023) Mari Carmen Bañuls, “Tensor network algorithms: A route map,” Annual Review of Condensed Matter Physics 14, 173–191 (2023).
  * Okunishi  _et al._ (2022) Kouichi Okunishi, Tomotoshi Nishino, and Hiroshi Ueda, “Developments in the tensor network — from statistical mechanics to quantum entanglement,” Journal of the Physical Society of Japan 91, 062001 (2022).
  * Gomez  _et al._ (2022) Abigail McClain Gomez, Susanne F. Yelin, and Khadijeh Najafi, “Reconstructing quantum states using basis-enhanced Born machines,” arXiv preprint 2206.01273 (2022).
  * Yang  _et al._ (2023) Frances Fengyi Yang, Michele Sasdelli, and Tat-Jun Chin, “Training multilayer perceptrons by sampling with quantum annealers,” arXiv preprint 2303.12352 (2023).
  * Dymarsky and Pavlenko (2022) Anatoly Dymarsky and Kirill Pavlenko, “Tensor network to learn the wave function of data,” Physical Review Research 4, 043111 (2022).
  * Lami and Collura (2023) Guglielmo Lami and Mario Collura, “Quantum magic via perfect sampling of matrix product states,” arXiv preprint arXiv:2303.05536 (2023).
  * Notarnicola  _et al._ (2023) Simone Notarnicola, Andreas Elben, Thierry Lahaye, Antoine Browaeys, Simone Montangero, and Benoît Vermersch, “A randomized measurement toolbox for an interacting Rydberg-atom quantum simulator,” New Journal of Physics 25, 103006 (2023).
  * Wild  _et al._ (2021) Dominik S. Wild, Dries Sels, Hannes Pichler, Cristian Zanoci, and Mikhail D. Lukin, “Quantum sampling algorithms, phase transitions, and computational complexity,” Phys. Rev. A 104, 032602 (2021).
  * Han  _et al._ (2018) Zhao-Yu Han, Jun Wang, Heng Fan, Lei Wang, and Pan Zhang, “Unsupervised generative modeling using matrix product states,” Phys. Rev. X 8, 031012 (2018).
  * Bravyi  _et al._ (2022) Sergey Bravyi, David Gosset, and Yinchen Liu, “How to simulate quantum measurement without computing marginals,” Phys. Rev. Lett. 128, 220503 (2022).
  * Zaletel and Pollmann (2020) Michael P. Zaletel and Frank Pollmann, “Isometric tensor network states in two dimensions,” Phys. Rev. Lett. 124, 037201 (2020).
  * (61) Harald Niederreiter, “7. random numbers and pseudorandom numbers,” in  _Random Number Generation and Quasi-Monte Carlo Methods_ (SIAM) pp. 161–176.
  * Devroye (2006) Luc Devroye, “Chapter 4 nonuniform random variate generation,” in  _Simulation_, Handbooks in Operations Research and Management Science, Vol. 13, edited by Shane G. Henderson and Barry L. Nelson (Elsevier, 2006) pp. 83–121.
  * Givens and Hoeting (2012) Geof H. Givens and Jennifer A. Hoeting, “Simulation and Monte Carlo Integration,” in  _Computational Statistics_ (John Wiley & Sons, Ltd, 2012) Chap. 6, pp. 151–199.
  * Ballarin  _et al._ (2024a) Marco Ballarin, Giovanni Cataldi, Aurora Costantini, Daniel Jaschke, Giuseppe Magnifico, Simone Montangero, Simone Notarnicola, Alice Pagano, Luka Pavesic, Marco Rigobello, Nora Reinić, Simone Scarlatella, and Pietro Silvi, “Quantum tea: qtealeaves,”  (2024a).
  * Fischer (2011) Hans Fischer,  _A history of the central limit theorem: from classical to modern probability theory_, Vol. 4 (Springer, 2011).
  * Nielsen and Chuang (2007) Michael A. Nielsen and Isaac L. Chuang,  _Quantum Computation and Quantum Information_, 9th ed. (Cambridge Univ. Press, Cambridge, United Kingdom, 2007).
  * Grover (1996) Lov K. Grover, “A fast quantum mechanical algorithm for database search,” in  _Proceedings of the twenty-eighth annual ACM symposium on Theory of computing_ (1996) pp. 212–219.
  * Shor (1994) P.W. Shor, “Algorithms for quantum computation: discrete logarithms and factoring,” in  _Proceedings 35th Annual Symposium on Foundations of Computer Science_ (1994) pp. 124–134.
  * Le  _et al._ (2011) Phuc Q. Le, Fangyan Dong, and Kaoru Hirota, “A flexible representation of quantum images for polynomial preparation, image compression, and processing operations,” Quantum Information Processing 10, 63–84 (2011).
  * Amankwah  _et al._ (2022) Mercy G. Amankwah, Daan Camps, E. Wes Bethel, Roel Van Beeumen, and Talita Perciano, “Quantum pixel representations and compression for n-dimensional images,” Scientific reports 12, 7712 (2022).
  * Sachdev (2011) Subir Sachdev,  _Quantum Phase Transitions_, 2nd ed. (Cambridge University Press, Cambridge, United Kingdom, 2011).
  * Haegeman  _et al._ (2016) Jutho Haegeman, Christian Lubich, Ivan Oseledets, Bart Vandereycken, and Frank Verstraete, “Unifying time evolution and optimization with matrix product states,” Phys. Rev. B 94, 165116 (2016).
  * Howard  _et al._ (2014) Mark Howard, Joel Wallman, Victor Veitch, and Joseph Emerson, “Contextuality supplies the ‘magic’ for quantum computation,” Nature 510, 351–355 (2014).
  * Seddon  _et al._ (2021) James R. Seddon, Bartosz Regula, Hakop Pashayan, Yingkai Ouyang, and Earl T. Campbell, “Quantifying quantum speedups: Improved classical simulation from tighter magic monotones,” PRX Quantum 2, 010345 (2021).
  * Leone  _et al._ (2022) Lorenzo Leone, Salvatore F. E. Oliviero, and Alioscia Hamma, “Stabilizer rényi entropy,” Phys. Rev. Lett. 128, 050402 (2022).
  * Ballarin  _et al._ (2024b) Marco Ballarin, Pietro Silvi, Simone Montangero, and Daniel Jaschke, “Figures and supplemental material for "Optimal sampling of tensor networks targeting wave function’s fast decaying tails",”  (2024b), supplemental materials contain the figures of the manuscript, additional figures and animations, data for the figures.
  * Ballarin  _et al._ (2024c) Marco Ballarin, Pietro Silvi, Simone Montangero, and Daniel Jaschke, “Simulation scripts for "Optimal sampling of tensor networks targeting wave function’s fast decaying tail",”  (2024c).
  * Ballarin  _et al._ (2024d) Marco Ballarin, Pietro Silvi, Simone Montangero, and Daniel Jaschke, “Input states for "Optimal sampling of tensor networks targeting wave function’s fast decaying tails",”  (2024d).
  * Stoudenmire and White (2010) E. M. Stoudenmire and Steven R. White, “Minimally entangled typical thermal state algorithms,” New Journal of Physics 12, 055026 (2010).
  * Chandra  _et al._ (2001) Rohit Chandra, Leo Dagum, David Kohr, Ramesh Menon, Dror Maydan, and Jeff McDonald,  _Parallel programming in OpenMP_ (Morgan Kaufmann, 2001).

## Supplemental material

The appendices focus on providing additional information which is necessary to
implement the OPES algorithm efficiently. Section .1 gives a brief example on
the scaling of both algorithms. Section .2 extends the pedagogical, simplified
algorithm to measuring multiple shots in superiterations, while Sec. .3
explains the tensor network aspects and how to cache information for the
tensor network.

###  .1 Scaling of the algorithms

We layout a brief calculation of how one can estimate the scaling of the two
algorithms and focus on the scaling for the coverage $\mathcal{C}$. Recall
that the threshold is $\epsilon=1-\mathcal{C}$. During the derivations, we
have to make assumptions on the probability distribution. As a reference
distribution, we introduce the series of probabilities

$\displaystyle p_{i}=2^{-i},i=1,\ldots 2^{n}-1\,,\qquad
p_{i=2^{n}}=2^{-(2^{n}-1)}\,,$ (27)

where the last probability $p_{i=2^{n}}$ ensures normalization. The number of
qubits is $n$.

For the standard sampling, we assume that we have to sample probability
intervals of the order of $\epsilon$ to reach the corresponding coverage. This
statement is exact for the above example in Eq. (27). Then, the number of
samples $j$ needed to find the outcome with probability $\epsilon$ scales as

$\displaystyle\mathcal{O}_{\textrm{Standard}}$
$\displaystyle=\frac{1}{2}\mathcal{O}\left(\frac{1}{\epsilon}\right)\,,$ (28)

where the factor of one-half stems from the fact of obtaining two
probabilities from each sample.

Turning to the OPES algorithm, we want to avoid dealing with the renormalized
probabilities after removing intervals. Therefore, we assume the samples are
drawn according to their probability, highest probability first. Then, we find
the number of probabilities to reach the given coverage, which is in case of
Eq. (27)

$\displaystyle\mathcal{O}_{\mathrm{OPES}}=$
$\displaystyle\frac{1}{2\log(2)}\mathcal{O}\left(\log\left(\frac{1}{\epsilon}\right)\right)\,.$
(29)

In this specific case of Eq. (27), we can also calculate the new probabilities
after $j$ iterations, still assuming we draw them in order

$\displaystyle p_{i>2j}^{j}$ $\displaystyle=2^{2j}\cdot p_{i}\,,$ (30)

where the superscript indicates the iteration in the sampling. Thus, the
favorable scaling of OPES originates in the log-scaling with $\epsilon^{-1}$
in comparison to the first-order polynomial scaling with $\epsilon^{-1}$ for
the standard sampling. We emphasize that this arguments are independent of the
size of the underlying Hilbert space.

To move forward and set up a numerical example, Eq. (27) is decaying too fast
to be captured by double precision. We generalize the decaying probabilities
to

$\displaystyle p_{i}=\frac{d^{-i}}{\mathcal{N}},i=1,\ldots 2^{n}\,,$ (31)

where normalization is guaranteed by the scalar $\mathcal{N}$ and the decay is
parameterized by $d$. In Fig. 8, we report an experiment with $n=10$ qubits
and $d=1.01$, where we show that the numerical results agree with the expected
scalings. The scalings from Eqs. (28) and (29) are shown as additional lines
for confirmation with a suitable prefactor. If the threshold $\epsilon$ is too
close to the minimal probability $\min(p_{i})$ or below, the scaling breaks
down as expected.

Figure 8: _Computational scaling._ We show the computational time $t$ needed
to obtain a given error $\epsilon$ for both methods. The probabilities are
proportional to $p_{i}\propto d^{-i}$ with $d=1.01$ for a system of $n=10$
qubits. The exponential speedup shows in the polynomial versus log scaling
with $\epsilon$.

The generalization to other discrete probability distributions is non-trivial.
For example, the discrete Gaussian distribution from Eq. (19) can be
considered for large system as continuous. Even then, we obtain a scaling
depending on the Gauss error function $\mathrm{erf}$. Moreover, an estimate if
the distribution has a smaller or larger benefit in comparison to Eq. (27) is
not conclusive: in case of the discrete Gaussian, both methods should have a
better scaling due to the smaller tails in the Gaussian versus Eq. (27).
Therefore, a case-by-case analysis is necessary, either empirically as
presented in Fig. 3, by an analytic prediction of the scaling, or numerically
constructing the scaling.

###  .2 Implementation of the algorithm

Section III contains already one version of the algorithm; now, we extend the
description to cover optional features like multiple shots, superiterations,
or cache access. The algorithm takes as input the number of random numbers to
be sampled $n_{r}$, the tensor network, and the set
$\Delta=[0,1]\setminus\overline{\Delta}=\\{\delta_{1},\ldots,\delta_{m}\\}$ of
$m$ known intervals $\delta_{\alpha}$. A superiteration of the algorithm has
the following steps:

  1. 1.

Set the state $\alpha^{*}$ to the empty bitstring, i.e., meaning no
measurement has taken place so far. Therefore, the interval is $\Delta=[0,1]$
and set the initial left boundary of your probability interval
$b^{l}_{\alpha^{*}}$ and the state probability $p_{\alpha^{*}}$:

$\displaystyle b^{l}_{\alpha^{*}}=0,\quad p_{\alpha^{*}}=1\,.$ (32)

  2. 2.

Sample $n_{r}$ uniform random variables $\\{u_{j}\\}_{j=1,\dots,n_{r}}$ from
the interval $\Delta$. Sort the random numbers $\\{u_{j}\\}_{j=1,\dots,n_{s}}$
in ascending order. This step ensures that the cached memory is optimally
managed.

  3. 3.

Pick the next random number $u_{j}$. If the cache is not empty, find the
longest possible bitstring matching $u_{j}$ and set $\alpha^{*}$ to the
longest matching bitstring; at least for the first qubit, there must be a
match amongst the known intervals. Then, check

     * •

if the cache is empty, continue with the next step.

     * •

if the $\alpha^{*}$ contains already all qubits, continue with this step and
$u_{j+1}$, i.e., $u_{j}\in\overline{\Delta}$.

     * •

else: load the corresponding tensors from the cache and continue directly with
step _5_.

  4. 4.

Compute the probability $p_{0}=p(q_{0}=\ket{0})$ of measuring $\ket{0}$ on
site $0$ and project it on a given state according to its probability and the
random number $u_{j}$:

$\displaystyle q_{0}$ $\displaystyle=\begin{cases}\ket{0}&\text{if
}u_{j}<p_{0},\\\ \ket{1}&\text{if }u_{j}\geq p_{0},\end{cases}$ (33)
$\displaystyle\ket{\alpha}$ $\displaystyle=\ket{\alpha^{*},q_{0}}\,.$ (34)

Update the left boundary and the state probability according to the measure:

$\displaystyle b^{l}_{\alpha}$
$\displaystyle=\begin{cases}b^{l}_{\alpha^{*}}&\text{if }u_{j}<p_{0},\\\
b^{l}_{\alpha^{*}}+p_{0}&\text{if }u_{j}\geq p_{0}\end{cases}$ (35)
$\displaystyle p_{\alpha}$ $\displaystyle=p_{\alpha^{*}}\cdot p(q_{0}).$ (36)

  5. 5.

Moving to the next qubit $i$, set $\alpha^{*}=\alpha$; $\alpha^{*}$ is now the
bitstring for the qubits up to the $(i-1)^{\mathrm{th}}$ qubit. Compute the
conditional probability $p_{0|\alpha^{*}}=p(q_{i}=\ket{0}|\alpha^{*})$. Select
the measured state according to $u_{j}$ as in the previous point and set
$b^{l}_{\alpha},p_{\alpha}$ accordingly, i.e.:

$\displaystyle q_{i}$ $\displaystyle=\begin{cases}\ket{0}&\text{if
}u_{j}<b^{l}_{\alpha^{*}}+p_{0|\alpha^{*}},\\\ \ket{1}&\text{if }u_{j}\geq
b^{l}_{\alpha^{*}}+p_{0|\alpha^{*}},\end{cases}$ (37)
$\displaystyle\ket{\alpha}$ $\displaystyle=\ket{\alpha^{*},q_{i}}\,,$ (38)
$\displaystyle b^{l}_{\alpha}$
$\displaystyle=\begin{cases}b^{l}_{\alpha^{*}}&\text{if
}u_{j}<b^{l}_{\alpha^{*}}+p_{0|\alpha^{*}}\,,\\\
b^{l}_{\alpha^{*}}+p_{0|\alpha^{*}}p_{\alpha^{*}}&\text{if }u_{j}\geq
b^{l}_{\alpha^{*}}+p_{0|\alpha^{*}}\,,\end{cases}$ (39) $\displaystyle
p_{\alpha}$ $\displaystyle=p_{\alpha^{*}}\cdot p(q_{i}|\alpha^{*}).$ (40)

  6. 6.

Repeat step _5_ until all the sites of the system are measured. Further
details for this step are explained in Sec. III and remain valid here.

  7. 7.

Remove the measured probability interval from our sampling interval, i.e.
$\Delta=\Delta\setminus\delta_{\alpha}$. Each new sampled number $u_{j}~{}\sim
U(\Delta)$ explores a different probability interval
$\delta_{\beta\neq\alpha}$. Update the cache with the corresponding
information and tensors.

  8. 8.

Continue with step _3_ if there are random numbers left for this
superiteration.

  9. 9.

If the number of superiterations or the desired coverage has been reached,
exit. If the next superiteration is executed with the current cache, continue
with step _2_. For the next superiteration with an empty cache, re-start with
step _1_.

###  .3 Details of the tensor networks implementation

Figure 9:  _Algorithm to sample from an $n$-sites MPS of local dimension $d=2$
and bond dimension $\chi$ with complexity $O(n\chi^{2}d^{2})$_. Tensors are
represented as circles and squares, connected legs are summed over. a)
Computation of the reduced density matrix $\rho^{(0)}$ using the isometry
center of the tensor network. b) Computation of the probability of measuring
the basis states $b_{0}=\ket{0}$ and $b_{1}=\ket{1}$ using the projector
$\\{P_{m}=\ket{b_{m}}\bra{b_{m}}\\}_{m=0,1}$. c) Projection of $T^{(0)}$ on
$b_{0}$ with probability $p_{0}(b_{0})$ or on $b_{1}$ with probability
$1-p_{0}(b_{0})$. d) Movement of the isometry center from $T^{(0)}$ to
$T^{(1)}$, repeating the same computations of a) and b)

The fundamental requirement for the sampling is the ability to efficiently
compute expectation values of local observables and the conditional
probabilities of a given quantum state. Here, we review how these values are
computed focusing on the standard sampling procedure, going afterward into the
details of what we cache to optimize the process.

Suppose we want to sample from the $n$-body quantum state $\ket{\psi}$
represented as an isometric tensor network Zaletel and Pollmann (2020) with
bond dimension $\chi$. In Fig. 9, we report the procedure for an MPS with
local dimension $d=2$. In general, the steps (a) to (d) following Fig. 9 are
the following:

  1. (a)

We set the isometry center to the left-most physical tensor $T^{(0)}$, in such
a way that the contraction over the auxiliary legs leads to the single-site
reduced density matrix $\rho_{0}$ with:

$\displaystyle\rho^{(0)}_{mn}=\sum_{i_{1},i_{2}\dots
i_{n}}T^{(0)}_{m,i_{1},i_{2}\dots i_{n}}\left(T^{(0)}_{n,i_{1},i_{2}\dots
i_{n}}\right)^{*},$ (41)

where $m,n$ represent the physical indexes of local dimension $d$, while the
auxiliary indexes, i.e. indexes not representing the Hilbert space of site
$0$, are defined as $i_{1},i_{2}\dots i_{n}$. The complex conjugate tensor is
indicated by a superscript $*$.

  2. (b)

Defining the computational basis of a $d$-dimensional qudits as
$\\{b_{m}\\}_{m=0,\dots,d-1}$, we compute the probability of site $0$ being in
state $b_{m}$ as:

$\displaystyle
p_{0}(b_{m})=\bra{\psi}P_{i}\ket{\psi}=\mathrm{Tr}\left[\rho^{(0)}P_{m}\right],$
(42)

where $\left\\{P_{m}=\ket{b_{m}}\bra{b_{m}}\right\\}_{m=0,\dots,d-1}$ are the
projectors over the computational basis.

  3. (c)

We randomly select a basis state $b_{m}$ with probability $p_{0}(b_{m})$ and
project the physical tensor $T^{(0)}$ on that state:

$\displaystyle T^{(0)}_{i_{1},i_{2}\dots
i_{n}}=\sum_{m}b_{m}T^{(0)}_{m,i_{1},i_{2}\dots i_{n}}\,.$ (43)

We call $q_{0}$ the measured basis state on site $0$. The projection is done
with a vector and not a projector for efficiency. Effectively, this step
reduces the local Hilbert space to a one-dimensional Hilbert space and reduces
the computational complexity.

  4. (d)

The isometry center is then moved to the next physical tensor $T_{1}$; by
repeating the same procedure, we obtain the conditional probabilities of the
second site:

$\displaystyle p_{1}(b_{m}|q_{0})=\mathrm{Tr}\left[\rho^{(1)}P_{m}\right]\,,$
(44)

and project into the basis state with these probabilities. Calling $q_{1}$ the
measured basis state on site $1$, we denote the probability as
$p(q_{1}|q_{0})$.

Iterating through the whole tensor network returns a single sample
$\alpha=q_{0},\dots,q_{{n-1}}$ from the quantum state:

$\displaystyle p_{\alpha}$ $\displaystyle=p_{q_{0},\dots,q_{{n-1}}}$ (45)
$\displaystyle=p\left(q_{{n-1}}|q_{{n-2}}\dots q_{{0}}\right)\dots
p(q_{1}|q_{0})p(q_{0}).$

Notice that following this procedure the state $\alpha$ is measured with
probability $p_{\alpha}$, thus this sampling method is unbiased, e.g., with
respect to the order of qubits.

The entire procedure scales as $O(n\chi^{2}d^{2})$ for MPS Stoudenmire and
White (2010), while as $O(n\chi^{3})$ for TTNs. Since each sample is
independent from the others the standard sampling can be trivially
parallelized with OPENMP Chandra  _et al._ (2001).

####  .3.1 Caching the MPS intermediate measurements

Figure 10: _Caching for matrix product states (MPS)._ Example of OPES with
$n=3$ sites and local dimension $d=2$. The intervals represent the conditional
probabilities of the state represented above it, i.e. the probability
$p_{0}(0)=0.5$, $p_{1}(0|p_{0}=0)=0.5$, and so on. First, we compute the
probability associated with sample $s_{0}=0.05$ (blue line), obtaining the
probability interval $\delta_{0}=(0,0.1)$ associated to the state $\ket{000}$.
We cache the couples $(0,\widetilde{T}_{0})$, $(00,\widetilde{T}_{1})$, where
with $\widetilde{T}_{i}$ we indicate the tensor after the projection with the
measured basis state. When we continue with the sample $s_{1}=0.29$, we use
the cached information and do not perform any operation on $T_{0}$. At the
level of $T_{1}$, we go to the right instead of the left; thus, we can
eliminate the cached information for $T_{1},T_{2}$. We can do so because any
further samples go to the right of the current interval, due to the ascending
ordering.

One of the reasons for the speedup of OPES over the standard sampling is the
caching of the intermediate states inside a superiteration. We suppose the
left-to-right sampling ordering discussed in Sec. .3, and that all
$\\{T_{j}\\}_{0\leq j<m}$, with $m>0$, have been measured. To restart the
computation from tensor $T_{m}$, it is sufficient to store the measured state
$|b_{0},b_{1},\dots,b_{m-1}\rangle$ and the tensor $T_{m}$ contracted with
tensor $T_{m-1}$ after the projection on state $|b_{m-1}\rangle$. The caching
in the superiteration is optimal because, due to the ascending ordering of the
samples $n_{s}$, we know that once the basis $b_{i}$ of the tensor $T_{i}$
pass from $0$ to $1$ we can delete all the cached information for $j<i$ since
we completely explored that probability interval. The caching structure is
shown in Fig. 10 for an MPS with $n=3$ sites and $d=2$ basis. Performance can
be further optimized by keeping the cache between superiterations.

This procedure can be parallelized to using OPENMP Chandra  _et al._ (2001),
by assigning to different processors different "branches" of the tree, such as
the blue and orange path from tensor $T_{1}$ in Fig. 10. Ensuring an equal
workload for each process is more challenging for the OPES than for the
standard sampling.

####  .3.2 Caching the TTN intermediate measurements

Figure 11: _Caching for tree tensor networks (TTN)._ The caching for the TTN
distinguishes which sites have been measured so far. Therein, an auxiliary
tensor in the upper layers of the TTN can be traversed up to three times. The
cache recognizes the direction based on the length of the state.

The caching in the MPS is simplified due to the fact that the MPS is a chain
and we respect the order of the MPS and treat one site after another; no
tensor in the MPS is visited twice. The TTN does not fulfill this assumption
and tensors can be visited up to three times during the projective
measurements.

Figure 11 sketches the problem for a 16-site TTN and one path towards the
solution. The tensors that are in the bulk, with respect to the layers and
within a layer itself, are traversed three times with QR decompositions while
running the projective measurements on the children. The tensor highlighted in
red in Fig. 11 is one example. Therefore, the cache stores tensors for the
same position of the TTN multiple times, distinguished by the length of the
bitstring measured so far.