[1]Institute of Physics, École Polytechnique Fédérale de Lausanne (EPFL),
CH-1015 Lausanne, Switzerland [2]Center for Quantum Science and Engineering,
École Polytechnique Fédérale de Lausanne (EPFL), CH-1015 Lausanne, Switzerland
[3]CPHT, CNRS, École polytechnique, Institut Polytechnique de Paris, 91120
Palaiseau, France

# Unbiasing time-dependent Variational Monte Carlo by projected quantum
evolution

Alessandro Sinibaldi  Clemens Giuliani  Giuseppe Carleo  Filippo Vicentini

###### Abstract

We analyze the accuracy and sample complexity of variational Monte Carlo
approaches to simulate the dynamics of many-body quantum systems classically.
By systematically studying the relevant stochastic estimators, we are able to:
(i) prove that the most used scheme, the time-dependent Variational Monte
Carlo (tVMC), is affected by a systematic statistical bias or exponential
sample complexity when the wave function contains some (possibly approximate)
zeros, an important case for fermionic systems and quantum information
protocols; (ii) show that a different scheme based on the solution of an
optimization problem at each time step is free from such problems; (iii)
improve the sample complexity of this latter approach by several orders of
magnitude with respect to previous proofs of concept. Finally, we apply our
advancements to study the high-entanglement phase in a protocol of non-
Clifford unitary dynamics with local random measurements in 2D, first
benchmarking on small spin lattices and then extending to large systems.

##  1 Introduction

The boundaries of current computational paradigms shape the scope of questions
that can be investigated in many-body quantum systems. Problems such as the
dynamics of high-dimensional interacting systems [1], dissipative phase
transitions out of equilibrium [2], the simulation of digital quantum circuits
[3], or quantum information protocols [4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14]
on states with volume-law entanglement all suffer from a lack of efficient
computational methods. A class of powerful numerical techniques to treat such
problems are variational methods, which rely on an efficient parametrization
of the quantum state and stochastic optimization of its parameters. Compared
to more established Tensor Network (TN) [15, 16] or Quantum Monte Carlo (QMC)
[17, 18] algorithms, variational approaches coupled with Monte Carlo sampling
can simulate high-dimensional or unstructured systems while not suffering from
the sign-problem [19], making them ideal candidates to target molecules [20,
21] and fermionic [22, 23, 24] or frustrated matter [25]. However, the
optimization problem arising in variational calculations is generally non-
convex, making it hard to give general convergence guarantees. While state-of-
the-art results for the calculation of ground states [25, 26, 27] have been
obtained with Variational Monte Carlo (VMC) [28], variational simulations of
dynamics have yet to improve over existing approaches systematically.

Techniques for the variational time evolution rely on so-called variational
principles to recast Schrödinger’s differential equation for the wave function
onto non-linear differential equations for the parameters [29]. These latter
equations can be integrated with an explicit scheme, the _time-dependent
Variational Monte Carlo_ or tVMC [30, 31], or with an implicit method by
solving an optimization problem at every time-step [32, 33, 34, 35, 36]. The
first approach has been applied to many systems [37, 38, 39], but it struggled
to significantly improve upon benchmark methods due to several poorly-
understood challenges in the numerical integration. The second method is
conceptually more powerful than tVMC, but it has yet to be applied to
realistic systems due to an unexpectedly large computational overhead [35].

Figure 1: Sketch of the failure of the tVMC when the state features zeros
(red) and of the dynamics generated by the p-tVMC algorithm (blue). When a
state with zeros (or near zeros) in the wave function is encountered during
the tVMC evolution, such as $\ket*{\Psi_{\theta(t_{i+1})}}$, the variational
dynamics starts to detach from the exact solution due to a bias (or an
vanishing signal to noise ratio). In p-tVMC, the optimization problem of
projecting the exactly evolved state $\mathcal{U}\ket*{\Psi_{\theta(t)}}$ onto
the variational manifold $\mathcal{M}$ of the ansatz
$\ket{\Psi_{\tilde{\theta}}}$ is solved at each time-step. This is achieved by
minimizing a distance in the Hilbert space, which is the infidelity
$\mathcal{I}$ (as shown in the right panel).

In this manuscript, we systematically analyze the accuracy and efficiency of
stochastic variational methods to tackle dynamical problems. First, we
formalize the origin of the numerical challenges affecting tVMC by proving
that they arise from the Monte Carlo sampling, which may hide a bias or an
exponential cost when the wave function contains zeros, as is the case for
many physically-relevant problems. Then, we prove that the high overhead of
the implicit integration arises from poor scaling of the Monte Carlo sampling
and we derive a new scheme that lowers the computational cost by several
orders of magnitude. We call this scheme _projected_ tVMC (p-tVMC). Finally,
we apply the p-tVMC to simulate the dynamics of a quantum system undergoing
unitary evolution interspersed with random measurements, which is a
paradigmatic model for entanglement phase transitions [4, 5, 6, 7, 8, 9, 10,
11, 12, 13, 14]. Established methods can access 1D systems [4, 5, 6, 8, 12,
14] or higher dimensional Clifford dynamics [7, 9, 10, 13], but questions
remain on the nature of this transition in the case of non-Clifford 2D
dynamics. As a proof of concept, we investigate this regime which is
intractable to TNs due to the rapid entanglement growth and the higher
dimensionality but also to tVMC due to the projective measurements enforcing a
large number of zeros in the wave function.

##  2 Numerical challenges in time-dependent Variational Monte Carlo

We consider a quantum many-body system whose Hilbert space is spanned by the
basis states $\\{\ket{\sigma}\\}$, where $\sigma$ is a set of quantum numbers
that is treated as discrete. The state of the system $\ket{\Psi}$ can be
efficiently approximated by a variational ansatz $\ket{\Psi_{\theta}}$ whose
wave function $\Psi_{\theta}(\sigma)\equiv\bra{\sigma}\ket{\Psi_{\theta}}$ is
completely specified by a set of $P$ parameters
$\theta=(\theta_{1},\ldots,\theta_{P})$, thus we have:

$\ket{\Psi}\approx\ket{\Psi_{\theta}}=\sum_{\sigma}\Psi_{\theta}(\sigma)\ket{\sigma}.$
(1)

We consider computationally tractable ansätze, meaning that $P$ is
polynomially large in the system size and $\Psi_{\theta}(\sigma)$ can be
sampled and queried efficiently [40].

Within this framework, variational dynamics can be encoded onto time-dependent
parameters $\theta(t)$ such that $\ket*{\Psi_{\theta(t)}}$ approximates the
physical dynamics. In what follows, we focus on the unitary evolution of a
time-independent Hamiltonian $\mathcal{H}$, on complex $\theta$ and
holomorphic $\Psi_{\theta}(\sigma)$. However, the discussion is general and
also applies to non-Hermitian PT-symmetric Hamiltonians [41], imaginary time
evolution [18], or open quantum systems obeying the Lindblad Master Equation
[42] and can be extended to the non-holomorphic case.

The McLachlan’s variational principle [43] recasts the Schrödinger’s equation
$\frac{d\ket{\Psi_{\theta}}}{dt}=-i\mathcal{H}\ket{\Psi_{\theta}}$ at every
time onto the optimization problem:

$\min_{\dot{\theta}}\mathcal{D}(|\Psi_{\theta(t)+\delta
t\dot{\theta}(t)}\rangle,e^{-i\mathcal{H}\delta t}|\Psi_{\theta(t)}\rangle),$
(2)

where $\mathcal{D}$ is the Fubini-Study metric and $\delta t$ is a small time-
step. By keeping only the leading terms in $\delta t$ in Eq. 2, it is possible
to derive the following set of explicit equations of motion for $\theta(t)$:

$\dot{\theta}_{k}(t)=-i\sum_{k^{\prime}}(S^{-1})_{kk^{\prime}}F_{k^{\prime}}.$
(3)

$F_{k}$ are the variational forces and $S_{kk^{\prime}}$ is the Quantum
Geometric Tensor [44, 45, 46]. These two quantities are defined as:

$\displaystyle F_{k}$
$\displaystyle\textrm{=}\frac{\bra{{\partial_{\theta_{k}}\Psi_{\theta}}}\mathcal{H}\ket{\Psi_{\theta}}}{\bra{\Psi_{\theta}}\ket{\Psi_{\theta}}}\text{\scalebox{0.95}[1.0]{$-$}}\frac{\bra{\partial_{\theta_{k}}\Psi_{\theta}}\ket{\Psi_{\theta}}}{\bra{\Psi_{\theta}}\ket{\Psi_{\theta}}}\frac{\bra{\Psi_{\theta}}\mathcal{H}\ket{\Psi_{\theta}}}{\bra{\Psi_{\theta}}\ket{\Psi_{\theta}}},$
(4) $\displaystyle S_{kk^{\prime}}$
$\displaystyle\textrm{=}\frac{\bra{\partial_{\theta_{k}}\Psi_{\theta}}\ket{\partial_{\theta_{k^{\prime}}}\Psi_{\theta}}}{\bra{\Psi_{\theta}}\ket{\Psi_{\theta}}}\text{\scalebox{0.95}[1.0]{$-$}}\frac{\bra{\partial_{\theta_{k}}\Psi_{\theta}}\ket{\Psi_{\theta}}}{\bra{\Psi_{\theta}}\ket{\Psi_{\theta}}}\frac{\bra*{\Psi_{\theta}}\ket*{\partial_{\theta_{k^{\prime}}}\Psi_{\theta}}}{\bra{\Psi_{\theta}}\ket{\Psi_{\theta}}},$
(5)

where we use the notation $\ket{\Psi_{\theta}}\equiv|\Psi_{\theta(t)}\rangle$.
The equations of motion Eq. 3 are integrated with numerical schemes like Euler
or higher-order Runge-Kutta.

The quantities in Eqs. 4 and 5 are efficiently computed by estimating the
expectation values as statistical averages over the Born distribution
$\Pi(\sigma)=|\Psi_{\theta}(\sigma)|^{2}/\bra{\Psi_{\theta}}\ket{\Psi_{\theta}}$
with Monte Carlo sampling, following the scheme known as _time-dependent
Variational Monte Carlo_ (tVMC) [30, 31]. The resulting stochastic expressions
are given by:

$\displaystyle F_{k}^{\text{MC}}$
$\displaystyle=\mathbb{E}_{\Pi}[O^{*}_{k}(\sigma)(E_{\mathrm{loc}}(\sigma)-\mathbb{E}_{\Pi}[E_{\mathrm{loc}}(\sigma)])]$
(6) $\displaystyle S_{kk^{\prime}}^{\text{MC}}$
$\displaystyle=\mathbb{E}_{\Pi}[O^{*}_{k}(\sigma)(O_{k^{\prime}}(\sigma)-\mathbb{E}_{\Pi}[O_{k^{\prime}}(\sigma)])].$
(7)

In the previous relations, the quantity
$E_{\mathrm{loc}}(\sigma)=\sum_{\sigma^{\prime}}\mathcal{H}_{\sigma,\sigma^{\prime}}\Psi_{\theta}(\sigma^{\prime})/\Psi_{\theta}(\sigma)$
is the local energy and
$O_{k}(\sigma)=\partial_{\theta_{k}}\log\Psi_{\theta}(\sigma)$ are the log-
derivatives of the variational state.

However, we remark that if the wave function and its derivatives have non-
overlapping support, the stochastic estimate $F_{k}^{\text{MC}}$ and
$S_{kk^{\prime}}^{\text{MC}}$ may differ from the exact quantities in Eqs. 4
and 5 by bias terms as:

$\displaystyle
F_{k}=\underbrace{\sum_{\sigma:|\Psi_{\theta}(\sigma)|=0}\frac{\bra{\partial_{\theta_{k}}\Psi_{\theta}}\ket{\sigma}\langle\sigma|\mathcal{H}\ket{\Psi_{\theta}}}{\bra{\Psi_{\theta}}\ket{\Psi_{\theta}}}}_{\textstyle\text{bias}\,b_{F}}+F_{k}^{\text{MC}},$
(8) $\displaystyle
S_{kk^{\prime}}=\underbrace{\sum_{\sigma:|\Psi_{\theta}(\sigma)|=0}\frac{\bra{\partial_{\theta_{k}}\Psi_{\theta}}\ket{\sigma}\langle\sigma|\partial_{\theta_{k^{\prime}}}\Psi_{\theta}\rangle}{\bra{\Psi_{\theta}}\ket{\Psi_{\theta}}}}_{\textstyle\text{bias}\,b_{S}}+S_{kk^{\prime}}^{\text{MC}}.$
(9)

The previous relations show that $F_{k}^{\text{MC}}$ and
$S_{kk^{\prime}}^{\text{MC}}$ are biased if $\bra{\sigma}\ket{\Psi_{\theta}}$
vanishes on some configurations while
$\bra{\sigma}\ket*{\partial_{\theta_{k}}\Psi_{\theta}}$ do not, leading to a
mismatch between the tVMC and ideal variational dynamics.

Figure 2: (a-b) Dynamics of the $z$-magnetization
$\langle\sigma_{i}^{z}\rangle$ for: (a) the $N=1$ spin $1/2$ toy-model, where
the initial state $\ket{+}$ is rotated with $\sigma^{y}$, which is simulated
by Exact Diagonalization (ED), tVMC and p-tVMC; (b) a system of $N=4$ spins
$1/2$, initialized in $\ket{\Psi_{\varepsilon}}$ with $\varepsilon=1.5\cdot
10^{-4}$ and evolved via $\mathcal{H}_{\text{TFI}}$ ($J=h=1$), which is
simulated by ED and tVMC with increasing number of samples $N_{s}$ (see
colorbar). For the tVMC the time-step $\delta t=10^{-3}$ has been used, while
for the p-tVMC $\delta t=10^{-2}$. (c) Infidelity $\mathcal{I}$ among states
evolved with tVMC and with ED after a time $t_{f}$ starting from
$\ket{\Psi_{\varepsilon}}$ with increasing $N_{s}$ (see colorbar) and
different $\varepsilon$. The inset shows how the minimal $N_{s}$ required to
reach $\mathcal{I}=5\cdot 10^{-3}$ (see the dashed line in the main plot)
scales with $\varepsilon$ (markers) and a power-law fit (red line). In (a) for
tVMC the ansatz parametrizing the wave function amplitudes is used, while in
(b-c) a Restricted Boltzmann Machine (RBM) [47] with $\alpha=1$. (d)
Illustration of the Born distribution for the peaked states
$\ket{\Psi_{\varepsilon}}$.

This condition may arise from the variational encoding of several physically
relevant states, such as basis states $\ket*{\sigma}$, anti-symmetric wave
functions (e.g. Slater, Neural Backflow [48], …) or states generated by
digital quantum circuits. Another relevant class of affected states are those
that underwent projective measurements, which are commonly found in trajectory
unravelings of the Lindblad Master Equation [49, 50, 51] or in quantum
information measurement protocols [4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14]. We
remark that for continuous systems (such as a particle in free space), the
bias may only emerge from zeros in the bulk, since at infinity the wave
function and its derivative must vanish. A pictorial representation of the
breakdown of tVMC is shown in Fig. 1.

In realistic calculations, the variational wave function is often arbitrarily
close to zero without encoding nodes exactly. In such cases, the biases
$b_{F}$ and $b_{S}$ are zero. Still, we find that the variances of
$F^{\text{MC}}$ and $S^{\text{MC}}$ grow such that an exponential number of
samples is required to resolve those quantities with finite accuracy. This
phenomenon can be revealed by the signal-to-noise ratios (SNRs) of
$F^{\text{MC}}$ and $S^{\text{MC}}$ approaching zero. The SNR of a function
$f$ of random variable $\sigma$ with distribution $\Pi$ is:

$\text{SNR}_{\Pi}[f]=\sqrt{\frac{N_{s}|\mathbb{E}_{\Pi}[f(\sigma)]|^{2}}{\text{Var}_{\Pi}[f(\sigma)]}},$
(10)

where $N_{s}$ is the number of samples used to estimate the mean value and
variance empirically.

In the following, we first discuss a minimal example where finite biases
emerge, and we then consider a more realistic case where there are no biases
but for which we show that the SNRs go to zero.

###  2.1 Paradigmatic examples

We analyze a toy model where the biases are non-zero, and they break the tVMC
dynamics. The system is made by $N=1$ spin $1/2$ where the state is
parameterized by the ansatz
$\ket*{\Psi_{(\alpha,\beta)}}=\alpha\ket{\downarrow}+\beta\ket{\uparrow}$. For
$\ket{\Psi_{\theta}}=|\Psi_{(1,0)}\rangle=\ket{\downarrow}$, such that
$\bra*{\uparrow}\ket*{\Psi_{\theta}}=0$, and evolution generated by
Hamiltonian $\mathcal{H}=\sigma^{y}$, we have that $b_{F}$ and $b_{S}$ in Eqs.
8 and 9 are finite and the stochastic estimates differ from the exact values
as:

$\displaystyle F$ $\displaystyle=\begin{pmatrix}0\\\ i\\\ \end{pmatrix},$
$\displaystyle S$ $\displaystyle=\begin{pmatrix}0&0\\\ 0&1\\\ \end{pmatrix},$
(11) $\displaystyle F^{\text{MC}}$ $\displaystyle=\begin{pmatrix}0\\\ 0\\\
\end{pmatrix},$ $\displaystyle S^{\text{MC}}$
$\displaystyle=\begin{pmatrix}0&0\\\ 0&0\\\ \end{pmatrix},$ (12)

(see Section B.1 for the full calculation).

In Fig. 2(a) we show a simulation of an initial state $\ket{+}$ which is
rotated by the Hamiltonian $\mathcal{H}=\sigma^{y}$. At $t=\pi/4$ the state
becomes $\ket{\Psi_{\theta}}=\ket{\downarrow}$ and the tVMC evolution is
stuck, as $F^{\text{MC}}$ and $S^{\text{MC}}$ vanish. Similar considerations
hold for more than 1 particle, and in Section B.2 we discuss an example for
the GHZ state of $N=2$ spins.

We now analyze a more realistic case where the variational state does not
exactly encode zeros. In particular, we consider a system of $N$ spins in the
state $\ket{\Psi_{\varepsilon}}$, which is peaked on a single configuration
$\ket{\sigma_{0}}$ and has a constant small amplitude $\sqrt{\varepsilon}$ for
all the other basis states (see Fig. 2(d)), namely:

$\Psi_{\varepsilon}(\sigma)=\begin{cases}\sqrt{\varepsilon},&\text{if
$\sigma\neq\sigma_{0}$},\\\ \sqrt{1-(2^{N}-1)\varepsilon},&\text{if
$\sigma=\sigma_{0}$},\end{cases}$ (13)

where $0<\varepsilon<1/(2^{N}-1)$. We remark that for small $\varepsilon$ this
ansatz approximates a Hartree-Fock state, which commonly arises from quantum-
chemistry Hamiltonians. For $\varepsilon\neq 0$ the biases $b_{F/S}$ are zero,
but for $\varepsilon\approx 0$ the leading terms of the SNRs for
$F^{\text{MC}}$ and $S^{\text{MC}}$ are respectively (see Appendix C for the
full calculation):

$\text{SNR}_{\Pi}[F^{\text{MC}}]\propto\sqrt{\varepsilon},\quad\text{SNR}_{\Pi}[S^{\text{MC}}]\propto\sqrt{\varepsilon}.$
(14)

Intuitively, this suggests that the more peaked the state is, the more samples
will be needed to accurately estimate those quantities, in particular
$N_{s}\propto\varepsilon^{-1}$. As normalization of the state imposes
$\varepsilon\propto 2^{-N}$, the number of samples necessary to correctly
compute the quantum geometric tensor and the variational forces will diverge
as $N_{s}\propto 2^{N}$, eliminating the computational advantage of stochastic
sampling and rendering tVMC computationally ineffective.

We consolidate this argument with a numerical experiment involving a commonly
adopted setup for quantum dynamics. We evolve with tVMC the state
$\ket{\Psi_{\varepsilon}}$ for $t\in[0,t_{f}]$ according to the Transverse
Field Ising (TFI) Hamiltonian,

$\mathcal{H}_{\mathrm{TFI}}=-J\sum_{\langle
i,j\rangle}\sigma_{i}^{z}\sigma_{j}^{z}-h\sum_{i}\sigma_{i}^{x},$ (15)

where $\langle i,j\rangle$ denotes nearest neighbors in a lattice with
periodic boundary conditions. In Fig. 2(b) we show some evolutions obtained
with an increasing number of Monte Carlo samples $N_{s}$ for a fixed
$\varepsilon$, demonstrating that the dynamics is correctly reconstructed only
at large values of $N_{s}$. The scaling of $N_{s}$ with the system size is
studied in Fig. 2(c), where we report the final infidelity111The infidelity
between two arbitrary states $\ket{\psi}$ and $\ket{\phi}$ is defined as
$\mathcal{I}(\psi,\phi)=1-\absolutevalue{\innerproduct{\psi}{\phi}}^{2}/\bra{\psi}\ket{\psi}\bra{\phi}\ket{\phi}$.
of the state obtained with tVMC with respect to the exact solution for
different $\varepsilon$ and $N_{s}$. We remark that the accuracy of the
variational simulation improves when $\varepsilon$ or $N_{s}$ are increased,
as the statistical fluctuations in the estimated quantities are suppressed.
The inset highlights a power-law relation between the $N_{s}$ necessary to
reconstruct the dynamics accurately and $\varepsilon$, proving that $N_{s}\sim
2^{N}$.

###  2.2 Overview

In this first section, we have shown that the tVMC method can be either biased
or require an exponential number of samples when the wave function is exactly
or approximately zero. This highlights the necessity of an efficient
alternative method to tVMC for variational time evolution.

We stress that while our considerations on stochastic estimators arose in the
context of tVMC, they are also applicable to ground-state calculations using
both plain gradient descent or stochastic reconfiguration [52], because they
rely on the same stochastic estimators. However, we believe that in such
calculations, the additional errors contributed by the biasing or small SNR
are mitigated by the iterative optimization scheme, which may avoid the
accumulation of errors that instead affects dynamics. We also remark that
Monte Carlo variational methods for open quantum systems [53, 54, 55] are also
possibly affected by the same issues.

Going forward, in Appendix A we propose a modified estimator for the forces
$F$ for which the bias and the SNR problems are absent, and therefore it can
efficiently estimate the forces when the standard estimator $F^{\text{MC}}$
fails (see Appendices B and C). From our knowledge, this alternative estimator
has never been discussed in the literature, and preliminary investigations
suggest that it already reduces the computational effort needed to reliably
find the ground state of some frustrated or fermionic Hamiltonians.

Unfortunately, we could not find a similarly straightforward modification of
the estimator for the quantum geometric tensor $S$. For that reason, the
following section presents a completely different scheme that avoids using the
tensor.

##  3 Projected time-dependent Variational Monte Carlo

Figure 3: (a) Learning curves for the infidelity $\mathcal{I}$ using the bare
estimator $\real\mathcal{I}_{\text{loc}}$ and the CV estimator
$\mathcal{I}^{\text{CV}}_{\text{loc}}$. In the inset the corresponding
variances, generically indicated as $\text{Var}_{\chi}[\mathcal{I}]$, are
shown. (b) Signal to noise ratio of the infidelity
$\text{SNR}_{\chi}[\mathcal{I}]=\sqrt{N_{s}\mathcal{I}^{2}/\text{Var}_{\chi}[\mathcal{I}}]$
as a function of $\mathcal{I}$ for the bare and the CV estimators with
different values of $c$. In both (a) and (b) a system of $N=8$ spins $1/2$ is
considered. In (b) $N_{s}=10^{4}$ samples have been used.

We consider the general problem of finding the parameters of a variational
state $\ket*{\Psi_{\tilde{\theta}}}$ such that it approximates the state
$\mathcal{U}\ket{\Psi_{\theta}}$, where $\theta$ are known and $\mathcal{U}$
is an arbitrary transformation, in terms of a given distance. Considering the
distance to be the infidelity $\mathcal{I}$, this can be expressed as the
following optimization problem:

$\min_{\tilde{\theta}}\mathcal{I}(\ket*{\Psi_{\tilde{\theta}}},\,\mathcal{U}\ket{\Psi_{\theta}}).$
(16)

Other distance choices, such as the L2 metric, have also been discussed in the
literature [36]. Eq. 16 is similar to Eq. 2, but it can treat arbitrary
unitaries and therefore can be used to simulate non-infinitesimal gates in
quantum circuits [32, 33] or to perform state preparation.

The solution of Eq. 16 can be found with iterative gradient-based optimizers
such as Stochastic Gradient Descent [56], ADAM [57], Natural Gradient [58, 59,
44] or similar methods. Since this approach consists of projecting the exactly
evolved state $\mathcal{U}\ket{\Psi_{\theta}}$ onto the manifold of the
variational ansatz $\ket{\Psi_{\tilde{\theta}}}$, we name it _projected time-
dependent Variational Monte Carlo_ (p-tVMC), and this is pictorially
represented in Fig. 1.

The infidelity in Eq. 16 can be estimated through Monte Carlo sampling as
$\mathcal{I}(\tilde{\theta})=\mathbb{E}_{\chi}[\mathcal{I}_{\text{loc}}(\sigma,\eta)]$.
Many choices for the sampling distribution $\chi$ and the local estimator
$\mathcal{I}_{\text{loc}}$ are possible, but assuming that $\mathcal{U}$ is
unitary we can sample

$\mathcal{I}_{\text{loc}}(\sigma,\eta)=1-\frac{\bra{\sigma}\mathcal{U}\ket{\Psi_{\theta}}}{\bra{\sigma}\ket*{\Psi_{\tilde{\theta}}}}\frac{\bra{\eta}\mathcal{U}^{\dagger}\ket*{\Psi_{\tilde{\theta}}}}{\bra{\eta}\ket{\Psi_{\theta}}},$
(17)

from the joint Born distribution of the two states
$\chi(\sigma,\eta)=|\Psi_{\tilde{\theta}}(\sigma)|^{2}|\Psi_{\theta}(\eta)|^{2}\bra*{\Psi_{\tilde{\theta}}}\ket*{\Psi_{\tilde{\theta}}}\bra{\Psi_{\theta}}\ket{\Psi_{\theta}}$.
For non-unitary $\mathcal{U}$ it is necessary to sample from
$\mathcal{U}\ket{\Psi_{\theta}}$ instead of $\ket{\Psi_{\theta}}$ [32, 33].

We remark that estimating the infidelity using Eq. 17 is efficient if
$\mathcal{U}$ is $K$-local. When this is not the case, it can be factored in
several sub-terms $\mathcal{U}=\mathcal{U}_{1}\dots\mathcal{U}_{N}$ where each
term is $K$-local, and Eq. 16 must be solved for every sub-unitary
$\mathcal{U}_{i}$. In particular, the unitary propagator of a general
Hamiltonian can be decomposed with the Trotter-Suzuki decomposition [60, 61]
or with other expansions that are unitary up to leading order, such as the
Taylor series [35].

We now analyze the estimator $\mathcal{I}_{\text{\text{loc}}}$ according to
the same approach used in the previous section. We find that, in the limit of
$\mathcal{I}\rightarrow 0$, the SNR of the estimator scales as

$\text{SNR}_{\chi}[\mathcal{I}_{\text{loc}}]\propto\sqrt{\mathcal{I}}.$ (18)

This means that, as the optimization approaches the optimum of
$\mathcal{I}=0$, the number of samples needed to resolve the infidelity
increases like $\mathcal{I}^{-1}$ (see Appendix E for an analytical
calculation). This is systematically different from what happens when
minimizing the energy, where the SNR remains constant when close to the
solution, and a constant number of samples can be used to achieve arbitrarily
high precision.

To recover this behaviour in the case of the infidelity optimization we
propose a new estimator based upon the Control Variates (CV) technique [62]:

$\mathcal{I}_{\text{loc}}^{\text{CV}}=\real\mathcal{I}_{\text{loc}}-c(|1-\mathcal{I}_{\text{loc}}|^{2}-1).$
(19)

where $c\in\mathbb{R}$. The quantity $c$ can be chosen such that the
$\text{Var}_{\chi}[\mathcal{I}_{\text{loc}}^{\text{CV}}]$ attains a minimal
value. This optimal value of $c$, say $c^{*}$, depends on the parameters of
the ansatz, so it changes during an infidelity optimization. However, it is
possible to show (see Appendix E for analytical proof) that in the limit
$\ket{\Psi_{\tilde{\theta}}}\rightarrow\mathcal{U}\ket{\Psi_{\theta}}$,
$c^{*}$ is exactly $-1/2$. Therefore, we avoid the high cost of estimating
$c^{*}$ at each iteration of an optimization, and directly use the
asymptotiacally ideal estimator Eq. 19 with $c=-1/2$.

To further support this approach, we show in Fig. 3(a) that the CV estimator
$\mathcal{I}_{\text{loc}}^{\text{CV}}$ features a variance that is orders of
magnitude smaller than the one of the bare estimator
$\real\mathcal{I}_{\text{loc}}$, in such a way that the number of samples
needed for the optimization is reduced by almost three orders of magnitude.
Additionally, the scaling of the variance with $\mathcal{I}$ changes, such
that for $\mathcal{I}\rightarrow 0$ we have that:

$\text{SNR}_{\chi}[\mathcal{I}_{\text{loc}}^{\text{CV}}]=O(1),$ (20)

meaning that the SNR remains asymptotically constant (see analytical proof in
Appendix E). Indeed, in Fig. 3(b) one can see that the SNR of
$\real\mathcal{I}_{\text{loc}}$ goes to zero when $\mathcal{I}$ decreases,
while the SNR of the corrected estimator decreases at smaller values of
$\mathcal{I}$. When $c=-1/2$, $\mathcal{I}_{\text{loc}}^{\text{CV}}$ has a SNR
which is constant and larger than $1$ over the whole range of $\mathcal{I}$
considered, suggesting that our strategy is ideal.

This substantial improvement of the sampling cost using the CV estimator Eq.
19 makes the p-tVMC an efficiently scalable method for simulating large
systems, such that it can address system sizes that have not been investigated
before within this approach. A detailed analysis on the CV infidelity
estimator is present in Appendix E, with further extensions in Appendix F. As
shown in Fig. 2(a), in Section B.2 for the GHZ state and in Appendix D for
adiabatic evolution, the p-tVMC, since it is not affected by biases or
vanishing SNR, can simulate dynamics in cases where tVMC fails or is
inefficient.

##  4 Unitary dynamics with random measurements

In recent years, considerable interest has been devoted to studying
entanglement in many-body quantum systems subject to evolution and random
local measurements. This is a paradigmatic model of a quantum system coupled
to an external environment acting as a measurement apparatus. Therefore, it is
intimately related to the physics of open quantum systems. The competition
between the unitary evolution’s entangling action and the measurements’
localizing effect gives rise to a phase transition between _volume-law_ and
_area-law_ entanglement in the steady state of the dynamics. The order
parameter is the measurement rate. This phenomenology has been originally
investigated in quantum circuits [4, 5, 7, 9, 10, 11, 13, 14], and more
recently for continuous dynamics [6, 8, 12].

To the extent of our knowledge, numerical investigations have focused so far
on systems in one-dimension, integrable or evolving via efficiently simulable
[63] Clifford gates because of algorithmic limitations. However, several open
questions remain on the nature of such transitions in non-integrable 2D
systems or non-Clifford circuits, requiring novel computational paradigms.

In this concluding section, we leverage the p-tVMC to simulate the time
evolution generated by the (non-Clifford) 2D TFI model subject to random local
measurements. This problem cannot be treated efficiently with Tensor Network
methods because of the rapid entanglement growth and the exponential cost of
exact contractions in 2D. Moreover, as projective measurements insert
exponentially many zeros in the wave function amplitudes, the shortcomings of
tVMC discussed in this article emerge, resulting in an exponential cost. We
consider a 2D spin $1/2$ square lattice with side length $L$, such that the
total number of spins is $N=L^{2}$. We evolve the system to time $t_{f}$,
discretized into time-steps of duration $\delta t$. At each $t$, two
operations are performed on the system:

  * •

unitary evolution with $\mathcal{U}=e^{-i\mathcal{H}_{\text{TFI}}\delta t}$,
where $\mathcal{H}_{\text{TFI}}$ is the 2D TFI Hamiltonian of Eq. 15;

  * •

a projective measurement of each spin in the $\sigma^{z}$ basis independently
with probability $p$ (measurement rate).

The overall evolution is stochastic and non-unitary.

As the unitary propagator is not $K$-local, we decompose it with a second-
order Trotter scheme as:

$e^{-i\mathcal{H}_{\mathrm{TFI}}\delta t}=e^{-i\mathcal{H}_{zz}\frac{\delta
t}{2}}e^{-i\mathcal{H}_{x}\delta t}e^{-i\mathcal{H}_{zz}\frac{\delta
t}{2}}+O(\delta t^{3}),$ (21)

where $\mathcal{H}_{zz}=-J\sum_{\langle
i,j\rangle}\sigma_{i}^{z}\sigma_{j}^{z}$ and
$\mathcal{H}_{x}=-h\sum_{i}\sigma_{i}^{x}$.

Figure 4: Time evolution of the Rényi-2 entropy $S_{2}$ simulated with p-tVMC
for $\mathcal{H}_{\text{TFI}}$ interspersed with random local measurements
along $z$ in 2D $L\times L$ lattices. For $L=4$ and $L=5$ we also provide ED
benchmark results. Subsystems of increasing size $|A|$ up to the maximum
$\lfloor N/2\rfloor$ are considered, and the corresponding markers are
indicated with different colors according to the colorbar. The insets show the
scaling of $S_{2}$ as a function of $|A|$ in the steady-state for the three
lattices. $S_{2}$ for subsystems with equal boundary length is indicated using
the same marker. The initial state is $\bigotimes_{i=1}^{N}\ket{+}_{i}$ and
the parameters of $\mathcal{H}_{\text{TFI}}$ are $J=1/2$ and $h=h_{c}/4$. The
measurement rate is $p=0.01$ and the time interval is $\delta t=0.1$. The
results are averaged over $5$ trajectories. The ansatz is an RBM with
$\alpha=4$, endowed with the variational terms to exactly implement the
diagonal part of the propagator and the measurements. The number of samples is
$N_{s}=10^{4}$ for $L=4$, and $N_{s}=2\cdot 10^{4}$ for $L=5,6$.

We use the p-tVMC to apply each unitary obtained from the decomposition.
Still, we remark that as $\exp(-i\mathcal{H}_{zz}\delta t/2)$ is diagonal in
the chosen $\sigma_{z}$ basis, its p-tVMC optimization problem can be solved
analytically, as shown in Appendix G. Instead, the unitary containing
$\mathcal{H}_{x}$ is applied using the p-tVMC, factorizing the propagator into
a product of terms, each of which acts on a small subset of the spins. In this
way, the number of connected elements to compute is not exponentially large
with $N$. The unitary dynamics that we simulate is a global quench across the
critical point of the 2D $\mathcal{H}_{\text{TFI}}$, which is in
correspondence of the transverse field $h_{c}$ such that $h_{c}/J\approx
3.044$ [64, 65, 66, 67, 68]. In particular, the initial state is the
paramagnetic ground state in the limit $h\rightarrow\infty$, given by
$\bigotimes_{i=1}^{N}\ket{+}_{i}$, and this is evolved into the ferromagnetic
phase where $h<h_{c}$.

The measurement operators are $\mathcal{P}^{\pm}_{i}=(1\pm\sigma_{z}^{i})/2$
for $i\in[1,N]$. When a spin $i$ is measured, one of the two projectors
$\mathcal{P}^{\pm}_{i}$ is applied to the state according to the probabilities
$\bra{\Psi}\mathcal{P}^{\pm}_{i}|\Psi\rangle/\bra{\Psi}\Psi\rangle$.
Measurements can be exactly realized on a variational state
$\ket{\Psi_{\theta}}$ by modifying it into a new ansatz
$\ket{\Psi_{\phi,\theta}}$ with additional parameters
$\phi=(\phi_{1}^{\downarrow},\phi_{1}^{\uparrow},\phi_{2}^{\downarrow},\phi_{2}^{\uparrow},\ldots)$
defined as:

$\bra{\sigma}\ket{\Psi_{\phi,\theta}}=\prod_{i=1}^{N}(\phi_{i}^{\downarrow}\delta_{\sigma_{i},\downarrow}+\phi_{i}^{\uparrow}\delta_{\sigma_{i},\uparrow})\bra{\sigma}\ket{\Psi_{\theta}}.$
(22)

Measuring spin $i$ with outcome $\downarrow$ translates into setting
$\phi_{i}^{\uparrow}=0$ and $\phi_{i}^{\downarrow}\neq 0$ and vice versa for
the other outcome. The measurement probabilities are stochastically computed
using Monte Carlo sampling. The protocol of unitary evolution and random
measurements is repeated several times and the final result is obtained by
averaging over all these trajectories, as in the Monte Carlo wave function
method [49]. To study the entanglement growth of $\ket{\Psi_{\theta}}$, we
monitor the Rényi-2 entanglement entropy
$S_{2}(\rho_{A})=-\log_{2}\Tr\rho_{A}^{2}$, where $\rho_{A}$ is the reduced
density matrix of the state on a subsystem $A$. Indeed, the Rényi-2 entropy is
a lower bound for the Von Neumann entanglement entropy and it can be estimated
via Monte Carlo sampling as:

$S_{2}=-\log_{2}\Bigg{(}\mathbb{E}_{\begin{subarray}{c}\Pi(\sigma,\eta)\\\
\hskip
5.69046pt\Pi(\sigma^{\prime},\eta^{\prime})\end{subarray}}\bigg{[}\frac{\Psi_{\theta}(\sigma^{\prime},\eta)\Psi_{\theta}(\sigma,\eta^{\prime})}{\Psi_{\theta}(\sigma,\eta)\Psi_{\theta}(\sigma^{\prime},\eta^{\prime})}\bigg{]}\Bigg{)},$
(23)

where $\sigma,\sigma^{\prime}\in A$ and $\eta,\eta^{\prime}\in B$
(complementary of $A$). See Appendix H for a derivation of Eq. 23.

Fig. 4 shows the evolution of $S_{2}(\rho_{A})$ for subsystems of size
$\absolutevalue{A}\in[1,\lfloor N/2\rfloor]$ in lattices with $L=\\{4,5,6\\}$
and with measurement rate $p=0.01$. To assess the quality of the variational
simulations, we compare with ED for $L=4,5$ and we select a feature density
for the RBM ansatz employed (which determines its expressivity) that gives a
satisfying level of precision. Given the chosen hyper-parameters, there is an
excellent agreement for small subsystem sizes
($\absolutevalue{A}\lesssim\lfloor N/4\rfloor$) and a good agreement for the
largest partitions. The Rényi-2 entanglement entropy is 0 at $t=0$, as
expected for the initial product state, and it grows linearly over time,
plateauing to a value that is proportional to $|A|$. The Page-like curves [69]
in the insets of Fig. 4 suggest that with $p=0.01$ the steady-states belong to
a high-entanglement phase. However, due to possible finite size effects,
whether the scaling of $S_{2}$ with the subsystem size is linear, witnessing a
volume-law, or logarithmic as in critical phases is not obvious. In any case,
a low-entanglement regime with area-law scaling is excluded since, in that
situation, the steady-state $S_{2}$ would not change for subsystems with the
same boundary length (indicated with equal markers in the insets). Instead,
what is observed is that $S_{2}$ increases with $|A|$ independently of the
boundary length, at least far from $\lfloor N/2\rfloor$ where finite size
effects might play a role. The proportionality of the entanglement growth rate
in the initial times with the boundaries of the partitions is in accordance
with the Lieb-Robinson bound [70] valid for local Hamiltonians.

##  5 Conclusions

In this manuscript, we proved that the standard approach to Monte Carlo
variational dynamics, the tVMC, can be limited by a finite bias or by an
exponentially small signal-to-noise ratio when the wave function contains
nodes or is only approximately zero. This implies that the tVMC cannot
efficiently simulate the time evolution of physically-relevant cases such as
completely polarized wave functions, states arising from digital quantum
circuits or measurement processes, including the open dynamics with quantum
jumps. Subsequently, we have formalized an alternative scheme, which consists
in solving an optimization problem at each time step using the infidelity
distance, and we have introduced a novel stochastic estimator which makes this
approach viable and scalable to large systems. Finally, we showed that our
method can solve the lack of efficient algorithms to investigate the high-
entanglement phase in a protocol of non-Clifford unitary dynamics with local
random measurements in 2D. This enables future investigation into the physics
of several classes of systems, including measurement-induced phase transitions
in non-trivial models above 1D and the physics of dissipative systems, all of
which are currently limited by the available computational methods. In
particular, a direct application of the projected method would be the
variational simulation of quantum trajectories arising from unraveling the
Lindblad Master Equation.

##  6 Data availability

All the simulations have been performed using Netket 3 [71, 72] with MPI and
MPI4jax [73]. The code for the p-tVMC method can be found in [74].

##  7 Acknowledgments

F.V. warmly thanks A. Biella for inspirational, insightful discussions and V.
Savona for arguments motivating this study. We gratefully acknowledge J. Nys
and Z. Denis for fruitful comments and discussions. C.G. is supported by the
Swiss National Science Foundation under Grant No. 200021_200336.

## References

  * [1] I. M. Georgescu, S. Ashhab, and F. Nori.  “Quantum simulation”.  Rev. Mod. Phys. 86, 153–185 (2014).  doi: 10.1103/RevModPhys.86.153.
  * [2] F. Minganti, I. I. Arkhipov, A. Miranowicz, and F. Nori.  “Continuous dissipative phase transitions with or without symmetry breaking”.  New Journal of Physics 23, 122001 (2021).  doi: 10.1088/1367-2630/ac3db8.
  * [3] M. A. Nielsen and I. L. Chuang.  “Quantum computation and quantum information: 10th anniversary edition”.  Cambridge University Press. (2010).  doi: 10.1017/CBO9780511976667.
  * [4] B. Skinner, J. Ruhman, and A. Nahum.  “Measurement-induced phase transitions in the dynamics of entanglement”.  Phys. Rev. X 9, 031009 (2019).  doi: 10.1103/PhysRevX.9.031009.
  * [5] Y. Li, X. Chen, and M. P. A. Fisher.  “Measurement-driven entanglement transition in hybrid quantum circuits”.  Phys. Rev. B 100, 134306 (2019).  doi: 10.1103/PhysRevB.100.134306.
  * [6] Q. Tang and W. Zhu.  “Measurement-induced phase transition: A case study in the nonintegrable model by density-matrix renormalization group calculations”.  Phys. Rev. Res. 2, 013022 (2020).  doi: 10.1103/PhysRevResearch.2.013022.
  * [7] X. Turkeshi, R. Fazio, and M. Dalmonte.  “Measurement-induced criticality in $(2+1)$-dimensional hybrid quantum circuits”.  Phys. Rev. B 102, 014315 (2020).  doi: 10.1103/PhysRevB.102.014315.
  * [8] X. Turkeshi, A. Biella, R. Fazio, M. Dalmonte, and M. Schiró.  “Measurement-induced entanglement transitions in the quantum ising chain: From infinite to zero clicks”.  Phys. Rev. B 103, 224210 (2021).  doi: 10.1103/PhysRevB.103.224210.
  * [9] A. Lavasani, Y. Alavirad, and M. Barkeshli.  “Topological order and criticality in $(2+1)\mathrm{D}$ monitored random quantum circuits”.  Phys. Rev. Lett. 127, 235701 (2021).  doi: 10.1103/PhysRevLett.127.235701.
  * [10] O. Lunt, M. Szyniszewski, and A. Pal.  “Measurement-induced criticality and entanglement clusters: A study of one-dimensional and two-dimensional clifford circuits”.  Phys. Rev. B 104, 155111 (2021).  doi: 10.1103/PhysRevB.104.155111.
  * [11] H. Liu, T. Zhou, and X. Chen.  “Measurement-induced entanglement transition in a two-dimensional shallow circuit”.  Phys. Rev. B 106, 144311 (2022).  doi: 10.1103/PhysRevB.106.144311.
  * [12] X. Turkeshi, M. Dalmonte, R. Fazio, and M. Schirò.  “Entanglement transitions from stochastic resetting of non-hermitian quasiparticles”.  Phys. Rev. B 105, L241114 (2022).  doi: 10.1103/PhysRevB.105.L241114.
  * [13] P. Sierant, M. Schirò, M. Lewenstein, and X. Turkeshi.  “Measurement-induced phase transitions in $(d+1)$-dimensional stabilizer circuits”.  Phys. Rev. B 106, 214316 (2022).  doi: 10.1103/PhysRevB.106.214316.
  * [14] J. C. Hoke et al.  “Quantum information phases in space-time: measurement-induced entanglement and teleportation on a noisy quantum processor” (2023).  arXiv:2303.04792.
  * [15] S. R. White.  “Density matrix formulation for quantum renormalization groups”.  Phys. Rev. Lett. 69, 2863–2866 (1992).  doi: 10.1103/PhysRevLett.69.2863.
  * [16] R. Orús.  “A practical introduction to tensor networks: Matrix product states and projected entangled pair states”.  Annals of physics 349, 117–158 (2014).  doi: 10.1016/j.aop.2014.06.013.
  * [17] D. Ceperley and B. Alder.  “Quantum monte carlo”.  Science 231, 555–560 (1986).  doi: 10.1126/science.231.4738.555.
  * [18] W. M. C. Foulkes, L. Mitas, R. J. Needs, and G. Rajagopal.  “Quantum monte carlo simulations of solids”.  Rev. Mod. Phys. 73, 33–83 (2001).  doi: 10.1103/RevModPhys.73.33.
  * [19] M. Troyer and U.-J. Wiese.  “Computational complexity and fundamental limitations to fermionic quantum monte carlo simulations”.  Phys. Rev. Lett. 94, 170201 (2005).  doi: 10.1103/PhysRevLett.94.170201.
  * [20] K. Choo, A. Mezzacapo, and G. Carleo.  “Fermionic neural-network states for ab-initio electronic structure”.  Nature Communications11 (2020).  doi: 10.1038/s41467-020-15724-9.
  * [21] D. Pfau, J. S. Spencer, A. G. D. G. Matthews, and W. M. C. Foulkes.  “Ab initio solution of the many-electron schrödinger equation with deep neural networks”.  Phys. Rev. Research 2, 033429 (2020).  doi: 10.1103/PhysRevResearch.2.033429.
  * [22] Y. Nomura, A. S. Darmawan, Y. Yamaji, and M. Imada.  “Restricted boltzmann machine learning for solving strongly correlated quantum systems”.  Phys. Rev. B 96, 205152 (2017).  doi: 10.1103/PhysRevB.96.205152.
  * [23] J. Stokes, J. R. Moreno, E. A. Pnevmatikakis, and G. Carleo.  “Phases of two-dimensional spinless lattice fermions with first-quantized deep neural-network quantum states”.  Phys. Rev. B 102, 205122 (2020).  doi: 10.1103/PhysRevB.102.205122.
  * [24] J. Nys and G. Carleo.  “Variational solutions to fermion-to-qubit mappings in two spatial dimensions”.  Quantum 6, 833 (2022).  doi: 10.22331/q-2022-10-13-833.
  * [25] K. Choo, T. Neupert, and G. Carleo.  “Two-dimensional frustrated ${J}_{1}\text{$-$}{J}_{2}$ model studied with neural network quantum states”.  Phys. Rev. B 100, 125124 (2019).  doi: 10.1103/PhysRevB.100.125124.
  * [26] O. Sharir, Y. Levine, N. Wies, G. Carleo, and A. Shashua.  “Deep autoregressive models for the efficient variational simulation of many-body quantum systems”.  Phys. Rev. Lett. 124, 020503 (2020).  doi: 10.1103/PhysRevLett.124.020503.
  * [27] D. Wu et al.  “Variational benchmarks for quantum many-body problems” (2023) arXiv:2302.04919.
  * [28] W. L. McMillan.  “Ground state of liquid ${\mathrm{He}}^{4}$”.  Phys. Rev. 138, A442–A451 (1965).  doi: 10.1103/PhysRev.138.A442.
  * [29] X. Yuan, S. Endo, Q. Zhao, Y. Li, and S. C. Benjamin.  “Theory of variational quantum simulation”.  Quantum 3, 191 (2019).  doi: 10.22331/q-2019-10-07-191.
  * [30] G. Carleo, F. Becca, M. Schiró, and M. Fabrizio.  “Localization and glassy dynamics of many-body quantum systems”.  Scientific Reports2 (2012).  doi: 10.1038/srep00243.
  * [31] G. Carleo, L. Cevolani, L. Sanchez-Palencia, and M. Holzmann.  “Unitary dynamics of strongly interacting bose gases with the time-dependent variational monte carlo method in continuous space”.  Phys. Rev. X 7, 031026 (2017).  doi: 10.1103/PhysRevX.7.031026.
  * [32] B. Jónsson, B. Bauer, and G. Carleo.  “Neural-network states for the classical simulation of quantum computing” (2018).  arXiv:1808.05232.
  * [33] M. Medvidović and G. Carleo.  “Classical variational simulation of the quantum approximate optimization algorithm”.  npj Quantum Information 7, 1–7 (2021).  doi: 10.1038/s41534-021-00440-z.
  * [34] S. Barison, F. Vicentini, and G. Carleo.  “An efficient quantum algorithm for the time evolution of parameterized circuits”.  Quantum 5, 512 (2021).  doi: 10.22331/q-2021-07-28-512.
  * [35] K. Donatella, Z. Denis, A. L. Boité, and C. Ciuti.  “Dynamics with autoregressive neural quantum states: application to critical quench dynamics” (2022).  arXiv:2209.03241.
  * [36] I. L. Gutiérrez and C. B. Mendl.  “Real time evolution with neural-network quantum states”.  Quantum 6, 627 (2022).  doi: 10.22331/q-2022-01-20-627.
  * [37] M. Schmitt and M. Heyl.  “Quantum many-body dynamics in two dimensions with artificial neural networks”.  Phys. Rev. Lett. 125, 100503 (2020).  doi: 10.1103/PhysRevLett.125.100503.
  * [38] R. Verdel, M. Schmitt, Y.-P. Huang, P. Karpov, and M. Heyl.  “Variational classical networks for dynamics in interacting quantum matter”.  Phys. Rev. B 103, 165103 (2021).  doi: 10.1103/PhysRevB.103.165103.
  * [39] M. Schmitt, M. M. Rams, J. Dziarmaga, M. Heyl, and W. H. Zurek.  “Quantum phase transition dynamics in the two-dimensional transverse-field ising model”.  Science Advances 8, eabl6850 (2022).  doi: 10.1126/sciadv.abl6850.
  * [40] M. V. den Nest.  “Simulating quantum computers with probabilistic methods” (2010).  arXiv:0911.1624.
  * [41] N. Moiseyev.  “Non-hermitian quantum mechanics”.  Cambridge University Press. (2011).  doi: 10.1017/CBO9780511976186.
  * [42] S. Haroche and J.-M. Raimond.  “Exploring the Quantum: Atoms, Cavities, and Photons”.  Oxford University Press. (2006).  doi: 10.1093/acprof:oso/9780198509141.001.0001.
  * [43] A. McLachlan.  “A variational solution of the time-dependent schrodinger equation”.  Molecular Physics 8, 39–44 (1964).  doi: 10.1080/00268976400100041.
  * [44] J. Stokes, J. Izaac, N. Killoran, and G. Carleo.  “Quantum Natural Gradient”.  Quantum 4, 269 (2020).  doi: 10.22331/q-2020-05-25-269.
  * [45] C.-Y. Park and M. J. Kastoryano.  “Geometry of learning neural quantum states”.  Phys. Rev. Res. 2, 023232 (2020).  doi: 10.1103/PhysRevResearch.2.023232.
  * [46] L. Hackl, T. Guaita, T. Shi, J. Haegeman, E. Demler, and J. I. Cirac.  “Geometry of variational methods: dynamics of closed quantum systems”.  SciPost Phys. 9, 048 (2020).  doi: 10.21468/SciPostPhys.9.4.048.
  * [47] G. Carleo and M. Troyer.  “Solving the quantum many-body problem with artificial neural networks”.  Science 355, 602–606 (2017).  doi: 10.1126/science.aag2302.
  * [48] D. Luo and B. K. Clark.  “Backflow transformations via neural networks for quantum many-body wave functions”.  Phys. Rev. Lett. 122, 226401 (2019).  doi: 10.1103/PhysRevLett.122.226401.
  * [49] K. Mølmer, Y. Castin, and J. Dalibard.  “Monte carlo wave-function method in quantum optics”.  J. Opt. Soc. Am. B 10, 524–538 (1993).  doi: 10.1364/JOSAB.10.000524.
  * [50] H. M. Wiseman and G. J. Milburn.  “Quantum measurement and control”.  Cambridge University Press. (2009).  doi: 10.1017/CBO9780511813948.
  * [51] F. Minganti, D. Huybrechts, C. Elouard, F. Nori, and I. I. Arkhipov.  “Creating and controlling exceptional points of non-hermitian hamiltonians via homodyne lindbladian invariance”.  Phys. Rev. A 106, 042210 (2022).  doi: 10.1103/PhysRevA.106.042210.
  * [52] S. Sorella, M. Casula, and D. Rocca.  “Weak binding between two aromatic rings: Feeling the van der Waals attraction by quantum Monte Carlo methods”.  The Journal of Chemical Physics127 (2007).  doi: 10.1063/1.2746035.
  * [53] F. Vicentini, A. Biella, N. Regnault, and C. Ciuti.  “Variational neural-network ansatz for steady states in open quantum systems”.  Phys. Rev. Lett. 122, 250503 (2019).  doi: 10.1103/PhysRevLett.122.250503.
  * [54] M. J. Hartmann and G. Carleo.  “Neural-network approach to dissipative quantum many-body dynamics”.  Phys. Rev. Lett. 122, 250502 (2019).  doi: 10.1103/PhysRevLett.122.250502.
  * [55] A. Nagy and V. Savona.  “Variational quantum monte carlo method with a neural-network ansatz for open quantum systems”.  Phys. Rev. Lett. 122, 250501 (2019).  doi: 10.1103/PhysRevLett.122.250501.
  * [56] I. Sutskever, J. Martens, G. Dahl, and G. Hinton.  “On the importance of initialization and momentum in deep learning”.  In S. Dasgupta and D. McAllester, editors, Proceedings of the 30th International Conference on Machine Learning.  Volume 28 of Proceedings of Machine Learning Research, pages 1139–1147.  Atlanta, Georgia, USA (2013). PMLR.  url: https://proceedings.mlr.press/v28/sutskever13.html.
  * [57] D. P. Kingma and J. Ba.  “Adam: A method for stochastic optimization” (2017).  arXiv:1412.6980.
  * [58] S.-i. Amari.  “Natural Gradient Works Efficiently in Learning”.  Neural Computation 10, 251–276 (1998).  doi: 10.1162/089976698300017746.
  * [59] S. Sorella.  “Wave function optimization in the variational monte carlo method”.  Phys. Rev. B 71, 241103 (2005).  doi: 10.1103/PhysRevB.71.241103.
  * [60] H. F. Trotter.  “On the product of semi-groups of operators”.  In Proceedings of the American Mathematical Society.  Volume 151, pages 545–551.  (1959).  doi: 10.1090/S0002-9939-1959-0108732-6.
  * [61] M. Suzuki.  “General theory of fractal path integrals with applications to many‐body theories and statistical physics”.  Journal of Mathematical Physics 32, 400–407 (1991).  doi: 10.1063/1.529425.
  * [62] R. Y. Rubinstein and D. P. Kroese.  “Simulation and the monte carlo method”.  John Wiley & Sons. (2016).  doi: 10.1002/9781118631980.
  * [63] D. Gottesman.  “The heisenberg representation of quantum computers” (1998).  arXiv:quant-ph/9807006.
  * [64] R. J. Elliott, P. Pfeuty, and C. Wood.  “Ising model with a transverse field”.  Phys. Rev. Lett. 25, 443–446 (1970).  doi: 10.1103/PhysRevLett.25.443.
  * [65] M. S. L. du Croo de Jongh and J. M. J. van Leeuwen.  “Critical behavior of the two-dimensional ising model in a transverse field: A density-matrix renormalization calculation”.  Phys. Rev. B 57, 8494–8500 (1998).  doi: 10.1103/PhysRevB.57.8494.
  * [66] H. Rieger and N. Kawashima.  “Application of a continuous time cluster algorithm to the two-dimensional random quantum ising ferromagnet”.  The European Physical Journal B - Condensed Matter and Complex Systems 9, 233–236 (1999).  doi: 10.1007/s100510050761.
  * [67] H. W. J. Blöte and Y. Deng.  “Cluster monte carlo simulation of the transverse ising model”.  Phys. Rev. E 66, 066110 (2002).  doi: 10.1103/PhysRevE.66.066110.
  * [68] A. F. Albuquerque, F. Alet, C. Sire, and S. Capponi.  “Quantum critical scaling of fidelity susceptibility”.  Phys. Rev. B 81, 064418 (2010).  doi: 10.1103/PhysRevB.81.064418.
  * [69] D. N. Page.  “Average entropy of a subsystem”.  Phys. Rev. Lett. 71, 1291–1294 (1993).  doi: 10.1103/PhysRevLett.71.1291.
  * [70] S. Bravyi, M. B. Hastings, and F. Verstraete.  “Lieb-robinson bounds and the generation of correlations and topological quantum order”.  Physical review letters 97, 050401 (2006).  doi: 10.1103/PhysRevLett.97.050401.
  * [71] F. Vicentini et al.  “NetKet 3: Machine Learning Toolbox for Many-Body Quantum Systems”.  SciPost Phys. CodebasesPage 7 (2022).  doi: 10.21468/SciPostPhysCodeb.7.
  * [72] G. Carleo et al.  “NetKet: A machine learning toolkit for many-body quantum systems”.  SoftwareXPage 100311 (2019).  doi: 10.1016/j.softx.2019.100311.
  * [73] D. Häfner and F. Vicentini.  “mpi4jax: Zero-copy MPI communication of JAX arrays”.  Journal of Open Source Software6 (2021).  doi: 10.21105/joss.03419.
  * [74] A. Sinibaldi and F. Vicentini.  “NetKet_fidelity Github package, v0.0.1”.  https://github.com/netket/netket_fidelity (2023).
  * [75] V. Havlicek.  “Amplitude Ratios and Neural Network Quantum States”.  Quantum 7, 938 (2023).  doi: 10.22331/q-2023-03-02-938.

## Appendix

##  Appendix A Unbiased force estimate

An unbiased estimate for the variational forces $F_{k}$ can be obtained by
inserting the completeness relation
$\sum_{\sigma}\ket{\sigma}\bra{\sigma}=\mathds{1}$ between $\mathcal{H}$ and
$\ket{\Psi_{\theta}}$, instead of between
$\bra*{\partial_{\theta_{k}^{*}}\Psi_{\theta}}$ and $\mathcal{H}$ as done for
the standard tVMC estimate $F_{k}^{\text{MC}}$, yielding:

$\displaystyle\tilde{F}_{k}^{\text{MC}}$
$\displaystyle=\sum_{\sigma}\frac{\bra{{\partial_{\theta_{k}}\Psi_{\theta}}}\mathcal{H}\ket{\sigma}\bra{\sigma}\ket{\Psi_{\theta}}}{\bra{\Psi_{\theta}}\ket{\Psi_{\theta}}}-\sum_{\sigma,\sigma^{\prime}}\frac{\bra{\partial_{\theta_{k}}\Psi_{\theta}}\ket{\sigma}\bra{\sigma}\ket{\Psi_{\theta}}}{\bra{\Psi_{\theta}}\ket{\Psi_{\theta}}}\frac{\bra*{\Psi_{\theta}}\ket*{\sigma^{\prime}}\bra*{\sigma^{\prime}}\mathcal{H}\ket{\Psi_{\theta}}}{\bra{\Psi_{\theta}}\ket{\Psi_{\theta}}}=$
(24)
$\displaystyle=\mathbb{E}_{\Pi}\bigg{[}\frac{\bra{\partial_{\theta_{k}}\Psi_{\theta}}\mathcal{H}\ket{\sigma}}{\bra{\Psi_{\theta}}\ket{\sigma}}\bigg{]}-\mathbb{E}_{\Pi}[O^{*}_{k}(\sigma)]\mathbb{E}_{\Pi}[E_{\mathrm{loc}}(\sigma)],$

where, like in the main text, $E_{\mathrm{loc}}(\sigma)$ is the local energy
and $O_{k}(\sigma)$ are the log-derivatives of the ansatz.
$\tilde{F}_{k}^{\text{MC}}$ does not have a covariance form like
$F_{k}^{\text{MC}}$, therefore it has, in general, larger statistical
fluctuations when estimated using a finite number of samples. Moreover, since
$O_{k}(\sigma)$ cannot be used to compute the first term in Eq. 24, its
computational cost is generally higher than the standard estimator.

##  Appendix B Examples of biases in tVMC

###  B.1 One-spin system

We consider a system of $N=1$ spin $1/2$ whose state is represented by the
variational ansatz
$\ket{\Psi_{\theta}}=\alpha\ket{\downarrow}+\beta\ket{\uparrow}$ with
parameters $\theta=(\alpha,\beta)$. The evolution is generated by the
Hamiltonian $\mathcal{H}=\sigma_{y}$. For the choice $\alpha=1$ and $\beta=0$,
we have $\ket{\Psi_{\theta}}=\ket{\downarrow}$,
$\ket{\partial_{\alpha}\Psi_{\theta}}=\ket{\downarrow}$ and
$\ket{\partial_{\beta}\Psi_{\theta}}=\ket{\uparrow}$. In these conditions:

$F=\begin{pmatrix}0\\\ i\end{pmatrix},\quad S=\begin{pmatrix}0&0\\\ 0&1\\\
\end{pmatrix}.\,\,\,$ (25)

However, due to the covariance form of the standard tVMC estimates
$F_{k}^{\text{MC}}$ and $S_{k,k^{\prime}}^{\text{MC}}$ evaluated on samples
including only one value of $\sigma$, we have that $\forall\,k,k^{\prime}$:

$F^{\text{MC}}_{k}=O^{*}_{k}(\downarrow)E_{\text{loc}}(\downarrow)-O^{*}_{k}(\downarrow)E_{\text{loc}}(\downarrow)=0,\quad
S^{\text{MC}}_{kk^{\prime}}=O^{*}_{k}(\downarrow)O_{k^{\prime}}(\downarrow)-O^{*}_{k}(\downarrow)O_{k^{\prime}}(\downarrow)=0.$
(26)

We verify that, instead, the stochastic estimate $\tilde{F}^{\text{MC}}$ with
the alternative estimator proposed in the previous section is unbiased, since:

$\tilde{F}_{\alpha}^{\text{MC}}=\frac{\bra{\partial_{\alpha}\Psi_{\theta}}\mathcal{H}\ket{\downarrow}}{\bra{\Psi_{\theta}}\ket{\downarrow}}-O_{\alpha}^{*}(\downarrow)E_{\text{loc}}(\downarrow)=0,\quad\tilde{F}_{\beta}^{\text{MC}}=\frac{\bra{\partial_{\beta}\Psi_{\theta}}\mathcal{H}\ket{\downarrow}}{\bra{\Psi_{\theta}}\ket{\downarrow}}-O_{\beta}^{*}(\downarrow)E_{\text{loc}}(\downarrow)=i,$
(27)

as
$E_{\text{\text{loc}}}(\downarrow)=\mathcal{H}_{\downarrow\downarrow}\Psi_{\theta}(\downarrow)/\Psi_{\theta}(\downarrow)+\mathcal{H}_{\downarrow\uparrow}\Psi_{\theta}(\uparrow)/\Psi_{\theta}(\downarrow)=0$.

###  B.2 Two-spin system

We consider a system of $N=2$ spins $1/2$ whose state is represented by the
variational ansatz
$\ket{\Psi_{\theta}}=\alpha\ket{\downarrow\downarrow}+\beta\ket{\downarrow\uparrow}+\gamma\ket{\uparrow\downarrow}+\delta\ket{\uparrow\uparrow}$
with parameters $\theta=(\alpha,\beta,\gamma,\delta)$. For
$\theta=(1/\sqrt{2},0,0,1/\sqrt{2})$ the ansatz represents the
Greenberger–Horne–Zeilinger state
$\ket{\Psi_{\theta}}=\ket{\text{GHZ}}=(\ket{\uparrow\uparrow}+\ket{\downarrow\downarrow})/\sqrt{2}$.
The variational derivatives are
$\ket{\partial_{\alpha}\Psi_{\theta}}=\ket{\downarrow\downarrow}$,
$\ket{\partial_{\beta}\Psi_{\theta}}=\ket{\downarrow\uparrow}$,
$\ket{\partial_{\gamma}\Psi_{\theta}}=\ket{\uparrow\downarrow}$ and
$\ket{\partial_{\delta}\Psi_{\theta}}=\ket{\uparrow\uparrow}$. We consider the
dynamics generated by the TFI Hamiltonian $\mathcal{H}_{\text{TFI}}$ with
coupling $J$ and transverse field $h$. In these conditions, $F$ and $S$ differ
from $F^{\text{MC}}$ and $S^{\text{MC}}$ as:

$F=\begin{pmatrix}0\\\ -h\sqrt{2}\\\ -h\sqrt{2}\\\ 0\\\ \end{pmatrix},\quad
S=\begin{pmatrix}1/2&0&0&-1/2\\\ 0&1&0&0\\\ 0&0&1&0\\\ -1/2&0&0&1/2\\\
\end{pmatrix},\quad F^{\text{MC}}=\begin{pmatrix}0\\\ 0\\\ 0\\\
0\end{pmatrix},\quad S^{\text{MC}}=\begin{pmatrix}1/2&0&0&-1/2\\\ 0&0&0&0\\\
0&0&0&0\\\ -1/2&0&0&1/2\end{pmatrix}.$ (28)

One can verify that, like for the one spin case, the alternative estimator
$\tilde{F}^{\text{MC}}$ gives an unbiased estimate for the variational forces.
As shown in Eq. 28, $S^{\text{MC}}$ has a lower rank than $S$, while
$F^{\text{MC}}$ is identical to zero meaning that this dynamics starting from
$\ket{\text{GHZ}}$ cannot be evolved with tVMC. On the contrary, the p-tVMC is
not affected by any problem and can perform the evolution, as shown in Fig. 5.

Figure 5: Dynamics of the $z$-correlator
$\langle\mathcal{\sigma}_{1}^{z}\sigma_{2}^{z}\rangle$ computed with Exact
Diagonalization (ED), tVMC and p-tVMC for $N=2$ spins $1/2$ evolved with
$\mathcal{H}_{\text{TFI}}$ from the initial state $\ket{\text{GHZ}}$. The
parameters of $\mathcal{H}_{\text{TFI}}$ are $J=h=1$. The time-steps used are
$\delta t=10^{-3}$ for tVMC and $t=10^{-2}$ for p-tVMC. $N_{s}=10^{3}$ samples
have been employed.

##  Appendix C Signal to noise ratios in Monte Carlo estimates

We consider a system of $N$ spins $1/2$ and normalized variational states
$\ket{\Psi_{\varepsilon}}$ whose wave function is parameterized by a single
parameter $\varepsilon$ as:

$\Psi_{\varepsilon}(\sigma)=\begin{cases}\sqrt{\varepsilon},&\text{if
$\sigma\neq\sigma_{0}$},\\\ \sqrt{1-(2^{N}-1)\varepsilon},&\text{if
$\sigma=\sigma_{0}$},\end{cases}$ (29)

with $0<\varepsilon<1/(2^{N}-1)$ and for a given $\sigma_{0}$. In the
following, we prove that the signal-to-noise ratios (SNRs) of
$F^{\text{MC}}_{k}$ and $S_{kk^{\prime}}^{\text{MC}}$ scale as
$O(\sqrt{\varepsilon})$, namely they diminish indefinitely as the wave
function becomes more peaked around $\sigma_{0}$. In order to ensure
normalization of $\ket{\Psi_{\varepsilon}}$ over different system sizes,
$\varepsilon$ goes as $1/2^{N}$. Therefore, the two SNRs diminish
exponentially as $N$ increases, and so the number of samples needed to resolve
$F_{k}^{\text{MC}}$ and $S_{kk^{\prime}}^{\text{MC}}$ with finite precision is
exponentially large in the system size. Instead, the SNR of the unbiased
estimate $\tilde{F}_{k}^{\text{MC}}$ of Eq. 24 is $O(1)$ in $\varepsilon$,
enabling it to efficiently estimate the forces in the limit
$\varepsilon\rightarrow 0$, where $F_{k}^{\text{MC}}$ cannot be used, and
independently of the system size.

###### Proof.

The variational log-derivative of $\Psi_{\varepsilon}(\sigma)$ is:

$O_{\varepsilon}(\sigma)=\begin{cases}\dfrac{1}{2\varepsilon},&\text{if
$\sigma\neq\sigma_{0}$},\\\
-\dfrac{2^{N}-1}{2(1-(2^{N}-1)\varepsilon)},&\text{if
$\sigma=\sigma_{0}$}.\end{cases}$ (30)

We define $(1-(2^{N}-1)\varepsilon)\equiv\alpha$ for brevity and we observe
that $\alpha\rightarrow 1$ when $\varepsilon\rightarrow 0$. Both
$F_{k}^{\text{MC}}$ and $S_{kk^{\prime}}^{\text{MC}}$ are of the form
$\mathbb{E}_{\Pi}[\bar{A}(\sigma)\bar{B}(\sigma)]$, where
$\bar{A}(\sigma)=A(\sigma)-\mathbb{E}_{\Pi}[A(\sigma)]$ and
$\bar{B}(\sigma)=B(\sigma)-\mathbb{E}_{\Pi}[B(\sigma)]$ for some functions $A$
and $B$ of random variable $\sigma$ with distribution $\Pi$. In
$F_{k}^{\text{MC}}$ we have $A(\sigma)=O_{k}^{*}(\sigma)$ and
$B(\sigma)=E_{\text{loc}}(\sigma)$, while in $S_{kk^{\prime}}^{\text{MC}}$ we
have $A(\sigma)=O_{k}^{*}(\sigma)$ and $B(\sigma)=O_{k^{\prime}}^{*}(\sigma)$.
The variance of stochastic estimates of this form is given by:

$\text{Var}_{\Pi}[\bar{A}(\sigma)\bar{B}(\sigma)]=\mathbb{E}_{\Pi}[|\bar{A}(\sigma)\bar{B}(\sigma)|^{2}]-|\mathbb{E}_{\Pi}[\bar{A}(\sigma)\bar{B}(\sigma)]|^{2}.$
(31)

For the local energy we have:

$\displaystyle E_{\text{loc}}(\sigma\neq\sigma_{0})$
$\displaystyle=\sum_{\sigma^{\prime}\neq\sigma_{0}}\mathcal{H}_{\sigma\sigma^{\prime}}\frac{\Psi_{\varepsilon}(\sigma^{\prime})}{\Psi_{\varepsilon}(\sigma)}+\mathcal{H}_{\sigma\sigma_{0}}\frac{\Psi_{\varepsilon}(\sigma_{0})}{\Psi_{\varepsilon}(\sigma)}=\sum_{\sigma^{\prime}\neq\sigma_{0}}\mathcal{H}_{\sigma\sigma^{\prime}}+\mathcal{H}_{\sigma\sigma_{0}}\sqrt{\frac{\alpha}{\varepsilon}}\stackrel{{\scriptstyle\varepsilon\ll
1}}{{\approx}}\mathcal{H}_{\sigma\sigma_{0}}\sqrt{\frac{\alpha}{\varepsilon}},$
(32) $\displaystyle E_{\text{loc}}(\sigma_{0})$
$\displaystyle=\sum_{\sigma^{\prime}\neq\sigma_{0}}\mathcal{H}_{\sigma_{0}\sigma^{\prime}}\frac{\Psi_{\varepsilon}(\sigma^{\prime})}{\Psi_{\varepsilon}(\sigma_{0})}+\mathcal{H}_{\sigma_{0}\sigma_{0}}\frac{\Psi_{\varepsilon}(\sigma_{0})}{\Psi_{\varepsilon}(\sigma_{0})}=\sum_{\sigma^{\prime}\neq\sigma_{0}}\mathcal{H}_{\sigma_{0}\sigma^{\prime}}\sqrt{\frac{\varepsilon}{\alpha}}+\mathcal{H}_{\sigma_{0}\sigma_{0}}\stackrel{{\scriptstyle\varepsilon\ll
1}}{{\approx}}\mathcal{H}_{\sigma_{0}\sigma_{0}},$ (33)

where the notation
$\mathcal{H}_{\sigma\sigma^{\prime}}\equiv\bra{\sigma}\mathcal{H}\ket{\sigma^{\prime}}$
for any $\sigma,\sigma^{\prime}$ is used. Therefore, keeping leading order
terms in $1/\varepsilon$, the variances of the stochastic estimates are:

$\text{Var}_{\Pi}[F_{k}^{\text{MC}}]\stackrel{{\scriptstyle\varepsilon\ll
1}}{{\approx}}\frac{\alpha}{4\varepsilon^{2}}\sum_{\sigma\neq\sigma_{0}}|\mathcal{H}_{\sigma\sigma_{0}}|^{2},\quad\text{Var}_{\Pi}[S_{kk^{\prime}}^{\text{MC}}]\stackrel{{\scriptstyle\varepsilon\ll
1}}{{\approx}}\dfrac{(2^{N}-1)}{16\varepsilon^{3}}.$ (34)

Since
$F_{k}^{\text{MC}}\approx\sum_{\sigma\neq\sigma_{0}}\mathcal{H}_{\sigma\sigma_{0}}\sqrt{\alpha/4\varepsilon}$
and $S_{kk^{\prime}}^{\text{MC}}\approx(2^{N}-1)/4\varepsilon$ for
$\varepsilon\ll 1$, we have that:

$\text{SNR}_{\Pi}[F_{k}^{\text{MC}}]\stackrel{{\scriptstyle\varepsilon\ll
1}}{{\propto}}\frac{|\sum_{\sigma\neq\sigma_{0}}\mathcal{H}_{\sigma\sigma_{0}}|}{\sqrt{\sum_{\sigma\neq\sigma_{0}}|\mathcal{H}_{\sigma\sigma_{0}}|^{2}}}\sqrt{\varepsilon}\stackrel{{\scriptstyle\varepsilon\rightarrow
0}}{{\longrightarrow}}0,\quad\text{SNR}_{\Pi}[S_{kk^{\prime}}^{\text{MC}}]\stackrel{{\scriptstyle\varepsilon\ll
1}}{{\propto}}\sqrt{(2^{N}-1)\varepsilon}\stackrel{{\scriptstyle\varepsilon\rightarrow
0}}{{\longrightarrow}}0.$ (35)

The two SNRs go to zero when the state becomes more peaked because when
$\varepsilon\rightarrow 0$ we have that
$|\Psi_{\varepsilon}(\sigma)|^{2}\rightarrow 0$ for $\sigma\neq\sigma_{0}$, so
these configurations are rarely sampled, but the estimators of the biases for
$\sigma\neq\sigma_{0}$ increase as $\varepsilon$ is reduced. Indeed, for
$\sigma\neq\sigma_{0}$ we have that:

$|\partial_{\varepsilon}\Psi_{\varepsilon}(\sigma)|^{2}=\frac{1}{4\varepsilon}\stackrel{{\scriptstyle\varepsilon\rightarrow
0}}{{\longrightarrow}}\infty,\quad\partial_{\varepsilon}\Psi_{\varepsilon}(\sigma)\bra{\sigma}\mathcal{H}\ket{\Psi_{\varepsilon}}=\sum_{\sigma^{\prime}\neq\sigma_{0}}\frac{1}{2}\mathcal{H}_{\sigma\sigma^{\prime}}+\frac{1}{2}\mathcal{H}_{\sigma\sigma_{0}}\sqrt{\frac{\alpha}{\varepsilon}}\stackrel{{\scriptstyle\varepsilon\rightarrow
0}}{{\longrightarrow}}\infty.$ (36)

It is possible to verify that using the unbiased force estimate
$\tilde{F}_{k}^{\text{MC}}$ of Eq. 24 the SNR remains constant in the limit
$\varepsilon\rightarrow 0$, making it able to efficiently compute the
variational forces when the standard estimate $F_{k}^{\text{MC}}$ fails. This
is because the original sum of the first term in Eq. 24 already runs over the
points where the wave function is non-zero, so to obtain the corresponding
estimator it is sufficient to divide by $\bra*{\Psi_{\theta}}\ket*{\sigma}$
without excluding any point. Instead, to obtain the estimator of the first
term in the expression of $F_{k}^{\text{MC}}$ it is necessary to divide by
$|\Psi_{\theta}(\sigma)|^{2}$, implicitly excluding from the Monte Carlo
expression the points for which the bias $b_{F}$ can be non-zero. Considering
only the contribution of the first term in $\tilde{F}_{k}^{\text{MC}}$, since
this is the problematic one in the standard estimate and the second term is
common to $F_{k}^{\text{MC}}$, we obtain that:

$\text{Var}_{\Pi}[\tilde{F}_{k}^{\text{MC}}]\stackrel{{\scriptstyle\varepsilon\ll
1}}{{\propto}}\frac{1}{4\varepsilon}\sum_{\sigma\neq\sigma_{0}}\bigg{|}\sum_{\sigma^{\prime}\neq\sigma_{0}}\mathcal{H}_{\sigma^{\prime}\sigma}\bigg{|}^{2}.$
(37)

Therefore, since
$\tilde{F}_{k}^{\text{MC}}\approx\sum_{\sigma\neq\sigma_{0}}\mathcal{H}_{\sigma\sigma_{0}}\sqrt{\alpha/4\varepsilon}=F_{k}^{\text{MC}}$
for $\varepsilon\ll 1$, we have:

$\text{SNR}_{\Pi}[\tilde{F}_{k}^{\text{MC}}]\stackrel{{\scriptstyle\varepsilon\ll
1}}{{\propto}}\frac{|\sum_{\sigma\neq\sigma_{0}}\mathcal{H}_{\sigma\sigma_{0}}|}{\sqrt{\sum_{\sigma\neq\sigma_{0}}|\sum_{\sigma^{\prime}\neq\sigma_{0}}\mathcal{H}_{\sigma^{\prime}\sigma}|^{2}}}\sqrt{\alpha}\stackrel{{\scriptstyle\varepsilon\rightarrow
0}}{{\longrightarrow}}\text{const.}$ (38)

∎

##  Appendix D Example of vanishing SNRs in tVMC

We consider a chain of $N$ spins $1/2$, initially in the ground state of the
Hamiltonian $-\sum_{i}\sigma_{i}^{x}$, which evolves according to the time-
dependent Hamiltonian:

$\mathcal{H}(t)=\gamma(t)\sum_{i}\sigma_{i}^{z}+(\gamma(t)-1)\sum_{i}\sigma^{x}_{i},$
(39)

where $\gamma(t)$ oscillates between $0$ and $1$ with a triangular profile of
period $T$, namely $\gamma(t)=t/T$ if $0<t<T$ and $\gamma(t)=1-(t-T)/T$ if
$t>T$. It is known that for a sufficiently large $T$, so for an adiabatic
evolution, the state at $t=T$ is going to be $\varepsilon$-close to
$\bigotimes_{i=1}^{N}\ket{\downarrow}_{i}$, so an instance of the peaked
states $\ket{\Psi_{\varepsilon}}$ of the previous section for some
$\varepsilon$ depending on $T$.

While in this case the biases in the tVMC estimates are zero, in Fig. 6 we
show that the dynamics is correctly reconstructed for $t>T$ only when choosing
a sufficiently large number of $N_{s}=10^{4}$ samples (the Hilbert space in
this case has size $2^{10}\approx 10^{3}$), hinting at the exponentially small
SNRs of $F^{\text{MC}}$ and $S^{\text{MC}}$. As before, the p-tVMC can
efficiently simulate this dynamics with a fair number of $N_{s}=10^{3}$
samples, for which instead the tVMC fails.

Figure 6: Dynamics of the $z$-magnetization
$\langle\mathcal{\sigma}_{i}^{z}\rangle$ computed with Exact Diagonalization
(ED), tVMC and p-tVMC for $N=10$ spins $1/2$ adiabatically evolved with
$\mathcal{H}(t)$ of Eq. 39 with period $T=8$ from the ground state of
$-\sum_{i}\sigma_{i}^{x}$. The time-steps used are $\delta t=10^{-3}$ for tVMC
and $t=10^{-2}$ for p-tVMC. The ansatz is an RBM with $\alpha=1$.

##  Appendix E Infidelity estimator with control variates

The stochastic estimate of the infidelity between the states
$\ket*{\Psi_{\tilde{\theta}}}$ and $\mathcal{U}\ket{\Psi_{\theta}}$, where
$\ket*{\Psi_{\tilde{\theta}}}$ and $\ket{\Psi_{\theta}}$ are variational
states and $\mathcal{U}$ is an arbitrary transformation, is obtained by
inserting completeness relations and dividing by wave function amplitudes in
its definition as follows:

$\mathcal{I}=1-\frac{\bra{\Psi_{\tilde{\theta}}}\mathcal{U}\ket{\Psi_{\theta}}\bra{\Psi_{\theta}}\mathcal{U}^{\dagger}\ket*{\Psi_{\tilde{\theta}}}}{\bra{\Psi_{\tilde{\theta}}}\ket*{\Psi_{\tilde{\theta}}}\bra{\Psi_{\theta}}\ket{\Psi_{\theta}}}=\sum_{\sigma,\eta}\frac{|\Psi_{\tilde{\theta}}(\sigma)|^{2}}{\bra*{\Psi_{\tilde{\theta}}}\ket*{\Psi_{\tilde{\theta}}}}\frac{|\Psi_{\theta}(\eta)|^{2}}{\bra{\Psi_{\theta}}\ket{\Psi_{\theta}}}\bigg{[}1-\frac{\bra{\sigma}\mathcal{U}\ket{\Psi_{\theta}}}{\bra{\sigma}\ket*{\Psi_{\tilde{\theta}}}}\frac{\bra{\eta}\mathcal{U}^{\dagger}\ket*{\Psi_{\tilde{\theta}}}}{\bra{\eta}\ket{\Psi_{\theta}}}\bigg{]}=\mathbb{E}_{\chi}[\mathcal{I}_{\text{loc}}(\sigma,\eta)],$
(40)

where $\mathcal{I}_{\text{loc}}(\sigma,\eta)$ is the term in brackets and
$\chi(\sigma,\eta)=|\Psi_{\tilde{\theta}}(\sigma)|^{2}|\Psi_{\theta}(\eta)|^{2}/\bra*{\Psi_{\tilde{\theta}}}\ket*{\Psi_{\tilde{\theta}}}\bra{\Psi_{\theta}}\ket{\Psi_{\theta}}$.

Similarly, the (conjugate) gradient of the infidelity can be computed
stochastically as:

$\partial_{\tilde{\theta}^{*}_{k}}\mathcal{I}=\mathbb{E}_{\chi}[O^{*}_{k}(\sigma)\mathcal{I}_{\text{loc}}(\sigma,\eta)]-\mathbb{E}_{\chi}[O^{*}_{k}(\sigma)]\mathbb{E}_{\chi}[\mathcal{I}_{\text{loc}}(\sigma,\eta)],$
(41)

where
$O_{k}(\sigma)=\partial_{\tilde{\theta}_{k}}\log\Psi_{\tilde{\theta}}(\sigma)$.
We remark that the first term in the previous expression is free from the
problem present in $F^{\text{MC}}_{k}$ and $S^{\text{MC}}_{k}$ when
$\ket*{\Psi_{\tilde{\theta}}}\rightarrow\mathcal{U}\ket*{\Psi_{\theta}}$.
Indeed, in this limit we have that for $\sigma$ where
$\Psi_{\tilde{\theta}}(\sigma)=0$ the term
$\bra*{\partial_{\tilde{\theta}_{k}}\Psi_{\tilde{\theta}}}\ket*{\sigma}\bra*{\sigma}\mathcal{U}\ket{\Psi_{\theta}}$
vanishes, so the bias as in Eq. 8 disappears. Since we consider continuous
time dynamics generated by a succession of infinitesimal $\mathcal{U}$s, such
that $\ket*{\Psi_{\tilde{\theta}}}\approx\mathcal{U}\ket*{\Psi_{\theta}}$
already at the beginning of the optimizations, we can use the gradient Eq. 41
without incurring in the issues discussed for the tVMC.

The variance of $\mathcal{I}_{\text{loc}}$ can be directly computed from the
one of $\mathcal{F}_{\text{loc}}=1-\mathcal{I}_{\text{loc}}$ as:

$\text{Var}_{\chi}[\mathcal{I}_{\text{loc}}(\sigma,\eta)]=\mathbb{E}_{\chi}[|\mathcal{F}_{\text{loc}}(\sigma,\eta)|^{2}]-|\mathbb{E}_{\chi}[\mathcal{F}_{\text{loc}}(\sigma,\eta)]|^{2}=1-(1-\mathcal{I})^{2}=2\mathcal{I}-\mathcal{I}^{2}.$
(42)

This proves that $\mathcal{I}_{\text{loc}}$ has a variance bounded above by
$1$ (because $0\leq\mathcal{I}\leq 1$), as already noted in Ref. [75] of the
main text, and that it features a zero variance principle, since
$\text{Var}_{\chi}[\mathcal{I}_{\text{loc}}(\sigma,\eta)]\rightarrow 0$ when
the solution is approached ($\mathcal{I}\rightarrow 0$), as the energy in VMC
for ground state search. However, the SNR of $\mathcal{I}_{\text{loc}}$ is:

$\text{SNR}_{\chi}[\mathcal{I}_{\text{loc}}(\sigma,\eta)]=\sqrt{\frac{N_{s}\mathcal{I}^{2}}{\text{Var}_{\chi}[\mathcal{I}_{\text{loc}}(\sigma,\eta)]}}=\sqrt{\frac{N_{s}\mathcal{I}^{2}}{2\mathcal{I}-\mathcal{I}^{2}}}\stackrel{{\scriptstyle\mathcal{I}\ll
1}}{{\sim}}\sqrt{N_{s}\mathcal{I}},$ (43)

which vanishes when $\mathcal{I}\rightarrow 0$. This entails that when
approaching the solution the number of samples $N_{s}$ must diverge as
$1/\mathcal{I}$ to resolve the infidelity with finite accuracy.

We find that this problem of diverging sampling overhead can be eliminated by
applying the Control Variates (CV) technique on $\mathcal{I}_{\text{loc}}$.
Since $\mathcal{I}$ is real, $\real\mathcal{I}_{\text{loc}}$ can be considered
in place of $\mathcal{I}_{\text{loc}}$ in all the following. CV consists of
adding to an estimator an additional random variable that is correlated to the
estimator and for which the expectation value is known exactly, such that the
fluctuations of this variable cancel out the ones of the original estimator.
For the infidelity, we discovered that $|\mathcal{F}_{\text{loc}}|^{2}$
satisfies the required properties, since it is correlated to
$\mathcal{I_{\text{loc}}}$ and we have
$\mathbb{E}_{\chi}[|\mathcal{F}_{\text{loc}}|^{2}]=\mathbb{E}_{\chi}[|1-\mathcal{I}_{\text{loc}}|^{2}]=1$.
Therefore, we can define the CV infidelity estimator:

$\mathcal{I}_{\text{loc}}^{\text{CV}}=\real\mathcal{I}_{\text{loc}}-c(|1-\mathcal{I}_{\text{loc}}|^{2}-1),$
(44)

where $c\in\mathbb{R}$. $\mathcal{I}_{\text{loc}}^{\text{CV}}$ has the same
mean as $\real\mathcal{I}_{\text{loc}}$ for any $c$, while its variance
depends on $c$. Therefore, $c$ can be chosen such that
$\text{Var}_{\chi}[\mathcal{I}_{\text{loc}}^{\text{CV}}(\sigma,\eta)]$ is
minimized. This optimal value, say $c^{*}$, is:

$c^{*}=-\frac{\text{Cov}_{\chi}[\real\mathcal{F}_{\text{loc}}(\sigma,\eta),|\mathcal{F}_{\text{loc}}(\sigma,\eta)|^{2}]}{\text{Var}_{\chi}[|\mathcal{F}_{\text{loc}}(\sigma,\eta)|^{2}]},$
(45)

where $\text{Cov}_{\chi}[A(\sigma),B(\sigma)]$ is the covariance of functions
$A$ and $B$ of random variable $\sigma$ with distribution $\chi$. The
corresponding minimized variance is:

$\displaystyle\text{Var}_{\chi}[\mathcal{I}_{\text{loc}}^{\text{CV}}(\sigma,\eta)](c^{*})$
$\displaystyle=\text{Var}_{\chi}[\real\mathcal{I}_{\text{loc}}(\sigma,\eta)]\bigg{(}1-\frac{{c^{*}}^{2}\text{Var}_{\chi}[|\mathcal{F}_{\text{loc}}(\sigma,\eta)|^{2}]}{\text{Var}_{\chi}[\real\mathcal{I}_{\text{loc}}(\sigma,\eta)]}\bigg{)},$
(46)

which is smaller than
$\text{Var}_{\chi}[\real\mathcal{I}_{\text{loc}}(\sigma,\eta)]$. $c^{*}$ from
Eq. 45 depends on the parameters $\tilde{\theta}$ of the variational ansatz,
thus its value varies during the optimization. However, it is possible to
prove that when
$\ket*{\Psi_{\tilde{\theta}}}\rightarrow\mathcal{U}\ket{\Psi_{\theta}}$,
$c^{*}$ tends to $-1/2$, such that we can fix $c=-1/2$ when $\mathcal{I}$ is
smaller than a predetermined value and avoid to re-estimate $c^{*}$ at each
iteration. The limit value of $c^{*}$ can be proven analytically considering
that when
$\ket*{\Psi_{\tilde{\theta}}}\rightarrow\mathcal{U}\ket{\Psi_{\theta}}$ we
have $\real\mathcal{F}_{\text{loc}}\rightarrow 1$ and
$|\mathcal{F}_{\text{loc}}|^{2}\rightarrow 1$ with same order as
$\real^{2}\mathcal{F}_{\text{loc}}$. Therefore, computing
$\lim_{\ket*{\Psi_{\tilde{\theta}}}\rightarrow\mathcal{U}\ket{\Psi_{\theta}}}c^{*}$
reduces to solve $\lim_{x\rightarrow 1}-(x^{3}-x)/(x^{4}-1)=-1/2$.

For reasons similar to the ones previously discussed for the first term in Eq.
41, the estimation of the CV factor $|1-\mathcal{I}_{\text{loc}}|^{2}$ is not
affected by a bias or a vanishing SNR as $F_{k}^{\text{MC}}$ and
$S_{kk^{\prime}}^{\text{MC}}$ when
$\ket{\Psi_{\tilde{\theta}}}\rightarrow\mathcal{U}\ket{\Psi_{\theta}}$, and it
is precisely in this limit where CV acts to keep the SNR constant. Therefore,
$\mathcal{I}_{\text{loc}}^{\text{CV}}$ is a well-defined infidelity estimator.

##  Appendix F Infidelity estimator with importance sampling

The estimator $\real\mathcal{I}_{\text{loc}}(\sigma,\eta)$ is unbounded
because it contains ratios of wave function amplitudes. Therefore, it may
diverge if $\bra{\sigma}\mathcal{U}\ket{\Psi_{\theta}}$ and
$\bra{\sigma}\ket*{\Psi_{\tilde{\theta}}}$ or
$\bra{\eta}\mathcal{U}^{\dagger}\ket*{\Psi_{\tilde{\theta}}}$ and
$\bra{\eta}\ket*{\Psi_{\theta}}$ differ of orders of magnitudes. This is a
severe problem when evolving the variational state after local measurements,
since it may happen that on some configurations the state
$\ket{\Psi_{\theta}}$ is close to zero but
$\mathcal{U}^{\dagger}\ket*{\Psi_{\tilde{\theta}}}$ is not. Therefore, despite
these configurations are rarely sampled, the estimator
$\mathcal{I}_{\text{loc}}(\sigma,\eta)$ on them can be very large,
significantly skewing the statistics. A way to include these outliers in the
sampling is to perform importance sampling (see Ref. [62] of the main text),
namely by rewriting:

$\mathcal{I}=\mathbb{E}_{\chi}[\mathcal{I}_{\text{loc}}(\sigma,\eta)]=\mathbb{E}_{\chi^{\prime}}\bigg{[}\mathcal{I}_{\text{loc}}(\sigma,\eta)\frac{\chi(\sigma,\eta)}{\chi^{\prime}(\sigma,\eta)}\bigg{]},$
(47)

for a given distribution $\chi^{\prime}$. This latter can be chosen such that
it minimizes the variance of the estimator, leading to the expression:

$\mathcal{I}=\mathbb{E}_{\chi^{\prime}}[\mathcal{I}_{\text{loc}}(\sigma,\eta)w(\sigma,\eta)],\quad\chi^{\prime}(\sigma,\eta)=\frac{|\mathcal{I}_{\text{loc}}(\sigma,\eta)|\chi(\sigma,\eta)}{\mathbb{E}_{\chi}[|\mathcal{I}_{\text{loc}}(\sigma,\eta)|]},$
(48)

where
$w(\sigma,\eta)=\mathbb{E}_{\chi}[|\mathcal{I}_{\text{loc}}(\sigma,\eta)|]/|\mathcal{I}_{\text{loc}}(\sigma,\eta)|$.
Using $\mathcal{I}_{\text{loc}}^{\text{CV}}$ in place of
$\mathcal{I}_{\text{loc}}$, importance sampling can be combined with control
variates obtaining the estimator:

$\mathcal{I}=\mathbb{E}_{\chi^{\prime}}[\mathcal{I}_{\text{loc}}^{\text{CV}}(\sigma,\eta)w(\sigma,\eta)].$
(49)

The cost of importance sampling is almost $n$ times larger than the cost of
standard sampling, where $n$ is the number of spins on which $\mathcal{U}$
acts non-trivially.

##  Appendix G Exact application of a diagonal propagator

The exact solution $\tilde{\theta}$ of the general optimization problem in the
p-tVMC must satisfy:

$\bra*{\sigma}\ket*{\Psi_{\tilde{\theta}}}=C\bra{\sigma}\mathcal{U}\ket{\Psi_{\theta}},$
(50)

for all configurations $\sigma$ and any constant $C$ equal for all $\sigma$.
For certain $\mathcal{U}$ and variational states, Eq. 50 can be solved
exactly. In particular, if $\mathcal{U}$ is diagonal in the basis
$\\{\ket{\sigma}\\}$ and if it is possible to write
$\bra*{\sigma}\ket*{\Psi_{\tilde{\theta}}}=\bra{\sigma}\ket{\Psi_{\theta}}\bra{\sigma}\ket{\Psi_{\delta\theta}}$
$\forall\sigma$ with some update $\delta\theta$, Eq. 50 is reduced to:

$\bra{\sigma}\ket{\Psi_{\delta\theta}}=CU_{\sigma}^{*},$ (51)

where $U_{\sigma}$ is the eigenvalue of $\mathcal{U}$ for the eigenstate
$\ket{\sigma}$. In the following, we consider spin systems for simplicity. For
the RBM, Eq. 51 admits a solution when
$\mathcal{U}=R_{i}^{z}(\phi_{z})=e^{-i\phi_{z}\sigma_{i}^{z}}$ consisting in
an update of the visible bias only. If
$\mathcal{U}=R_{ij}^{zz}(\phi_{zz})=e^{-i\phi_{zz}\sigma_{i}^{z}\sigma_{j}^{z}}$,
instead, a simple parameter update is not sufficient. The transformation can
be exactly implemented by adding a hidden unit (see Refs. [32, 33] of the main
text). However, in general, starting with an arbitrary ansatz
$\ket{\Psi_{\theta}}$, it is possible to add two variational terms such that
both $R_{i}^{z}$ and $R_{ij}^{zz}$ can be exactly simulated with a parameter
change. Indeed, $\ket{\Psi_{\theta}}$ can be modified into a new ansatz
$|\Psi_{J^{(1)}_{i},J^{(2)}_{ij},\theta}\rangle$ with two additional
parameters $\\{J^{(1)}_{i},J^{(2)}_{ij}\\}$ which is defined as:

$|\Psi_{J^{(1)}_{i},J^{(2)}_{ij},\theta}\rangle=e^{-iJ^{(1)}_{i}\sigma^{z}_{i}}e^{-iJ^{(2)}_{ij}\sigma^{z}_{i}\sigma^{z}_{j}}\ket{\Psi_{\theta}}.$
(52)

The two exponential terms factorize as required to have Eq. 51, such that the
application of $R_{i}^{z}(\phi_{z})$ translates into the updates $\delta
J^{(1)}_{i}=\phi_{z}$, $\delta J^{(2)}_{ij}=0$ and $\delta\theta=0$, while for
$R_{ij}^{zz}(\phi_{zz})$ the changes $\delta J^{(1)}_{i}=0$, $\delta
J^{(2)}_{ij}=\phi_{zz}$ and $\delta\theta=0$ are required. Adding many two-
site terms in the ansatz, it is possible to simulate the dynamics of the
diagonal part of the TFI Hamiltonian $\mathcal{H}_{\text{TFI}}$ exactly.

##  Appendix H Stochastic estimator for the Rényi-2 entropy

An estimator for the Rényi-2 entropy $S_{2}=-\log_{2}\Tr\rho_{A}^{2}$ can be
obtained by explicitly writing the definition of the purity of the reduced
density matrix as:

$\Tr\rho_{A}^{2}=\Tr[\bigg{(}\Tr_{B}\frac{\ket{\Psi_{\theta}}\bra{\Psi_{\theta}}}{\bra{\Psi_{\theta}}\ket{\Psi_{\theta}}}\bigg{)}^{2}]=\sum_{\sigma^{\prime},\eta,\eta^{\prime}}\frac{\bra{\sigma^{\prime},\eta}\ket{\Psi_{\theta}}\bra{\Psi_{\theta}}\ket{\eta}\bra{\eta^{\prime}}\ket{\Psi_{\theta}}\bra{\Psi_{\theta}}\ket{\eta^{\prime},\sigma^{\prime}}}{\bra{\Psi_{\theta}}\ket{\Psi_{\theta}}\bra{\Psi_{\theta}}\ket{\Psi_{\theta}}},$
(53)

where $\Tr_{B}$ indicates the partial trace over $B$, $\sigma^{\prime}\in A$
and $\eta,\eta^{\prime}\in B$. Now, inserting a completeness relation in $A$
and dividing by wave function amplitudes one obtains:

$\displaystyle\Tr\rho_{A}^{2}$
$\displaystyle=\sum_{\begin{subarray}{c}\sigma,\sigma^{\prime}\\\
\eta,\eta^{\prime}\end{subarray}}\frac{|\Psi_{\theta}(\sigma,\eta)|^{2}}{\bra{\Psi_{\theta}}\ket{\Psi_{\theta}}}\frac{|\Psi_{\theta}(\sigma^{\prime},\eta^{\prime})|^{2}}{\bra{\Psi_{\theta}}\ket{\Psi_{\theta}}}\frac{\Psi_{\theta}(\sigma^{\prime},\eta)\Psi_{\theta}(\sigma,\eta^{\prime})}{\Psi_{\theta}(\sigma,\eta)\Psi_{\theta}(\sigma^{\prime},\eta^{\prime})}=\mathbb{E}_{\begin{subarray}{c}\Pi(\sigma,\eta)\\\
\hskip
5.69046pt\Pi(\sigma^{\prime},\eta^{\prime})\end{subarray}}\bigg{[}\frac{\Psi_{\theta}(\sigma^{\prime},\eta)\Psi_{\theta}(\sigma,\eta^{\prime})}{\Psi_{\theta}(\sigma,\eta)\Psi_{\theta}(\sigma^{\prime},\eta^{\prime})}\bigg{]}.$
(54)