Unbiasing time-dependent Variational Monte Carlo by projected quantum evolution

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.


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 volumelaw 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].
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 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 |Ψ θ(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 U |Ψ θ(t) ⟩ onto the variational manifold M of the ansatz |Ψθ⟩ is solved at each time-step.This is achieved by minimizing a distance in the Hilbert space, which is the infidelity I (as shown in the right panel).
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.

Numerical challenges in timedependent Variational Monte Carlo
We consider a quantum many-body system whose Hilbert space is spanned by the basis states {|σ⟩}, where σ is a set of quantum numbers that is treated as discrete.The state of the system |Ψ⟩ can be efficiently approximated by a variational ansatz |Ψ θ ⟩ whose wave function Ψ θ (σ) ≡ ⟨σ|Ψ θ ⟩ is completely specified by a set of P parameters θ = (θ 1 , . . ., θ P ), thus we have: We consider computationally tractable ansätze, meaning that P is polynomially large in the system size and Ψ θ (σ) can be sampled and queried efficiently [40].
Within this framework, variational dynamics can be encoded onto time-dependent parameters θ(t) such that |Ψ θ(t) ⟩ approximates the physical dynamics.In what follows, we focus on the unitary evolution of a time-independent Hamiltonian H, on complex θ and holomorphic Ψ θ (σ).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 nonholomorphic case.
The McLachlan's variational principle [43] recasts the Schrödinger's equation d|Ψ θ ⟩ dt = −iH |Ψ θ ⟩ at every time onto the optimization problem: where D is the Fubini-Study metric and δt is a small time-step.By keeping only the leading terms in δt in Eq. ( 2), it is possible to derive the following set of explicit equations of motion for θ(t): In the previous relations, the quantity are the log-derivatives of the variational state.
However, we remark that if the wave function and its derivatives have non-identical support, namely there exist configurations σ for which The previous relations show that F MC k and S MC kk ′ are biased if ⟨σ|Ψ θ ⟩ vanishes on some configurations while ⟨σ|∂ θ k Ψ θ ⟩ do not, leading to a mismatch between the tVMC and ideal variational dynamics.
This condition may arise from the variational encoding of several physically relevant states, such as basis states |σ⟩, 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 MC and S 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 MC and S MC approaching zero.The SNR of a function f of random variable σ with distribution Π is: In practical calculations, to ensure an accurate estimation with a finite number of samples N s we require This inequality guarantees that the effective signal (mean value) is larger than the statistical fluctuations and thus it can be resolved.
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.

Paradigmatic examples
We analyze a toy model where the biases are non-zero, and they break the tVMC dynamics.The system is  8) and ( 9) are finite and the stochastic estimates differ from the exact values as: (see Appendix B.1 for the full calculation).In Fig. 2(a) we show a simulation of an initial state |+⟩ which is rotated by the Hamiltonian H = σ y .At t = π/4 the state becomes |Ψ θ ⟩ = |↓⟩ and the tVMC evolution is stuck, as F MC and S MC vanish.Similar considerations hold for more than 1 particle, and in Appendix 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 |Ψ ε ⟩, which is peaked on a single configuration |σ 0 ⟩ and has a constant small amplitude √ ε for all the other basis states (see Fig. 2(d)), namely: where 0 < ε < 1/(2 N −1).We remark that for small ε this ansatz approximates a Hartree-Fock state, which commonly arises from quantum-chemistry Hamiltonians.For ε ̸ = 0 the biases b F/S are zero, but for ε ≈ 0 the leading terms of the SNRs for F MC and S MC are respectively (see Appendix C for the full calculation): 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 ∝ ε −1 .As normalization of the state imposes ε ∝ 2 −N , the number of samples necessary to correctly compute the quantum geometric tensor and the variational forces will diverge as N s ∝ 2 N , eliminating the 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 |Ψ ε ⟩ for t ∈ [0, t f ] according to the Transverse Field Ising (TFI) Hamiltonian, where ⟨i, j⟩ 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 ε, 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 infidelity 1 of the state obtained with tVMC with respect to the exact solution for different ε and N s .We remark that the accuracy of the variational simulation improves when ε 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 ε, proving that N s ∼ 2 N .

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,56] 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 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.

Projected time-dependent Variational Monte Carlo
We consider the general problem of finding the parameters of a variational state |Ψ θ ⟩ such that it approximates the state U |Ψ θ ⟩, where θ are known and U is an arbitrary transformation, in terms of a given distance.Considering the distance to be the infidelity I, this can be expressed as the following optimization problem: 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 noninfinitesimal 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 [57], ADAM [58], Natural Gradient [59,60,44] or similar methods.Since this approach consists of projecting the exactly evolved state U |Ψ θ ⟩ onto the manifold of the variational ansatz |Ψ θ ⟩, 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 I( θ) = E χ [I loc (σ, η)].Many choices for the sampling distribution χ and the local estimator I loc are possible, but assuming that U is unitary we can sample from the joint Born distribution of the two states 32,33].We remark that estimating the infidelity using Eq. ( 17) is efficient if U is K-local [61], namely it acts non-trivially on at most K degrees of freedom (spins, qubits, particles, . . .), where K is polynomially large in the system size.When this is not the case, it can be factored in several sub-terms U = U 1 . . .U N where each term is K-local, and Eq. ( 16) must be solved for every sub-unitary U i .In particular, the unitary propagator of a general Hamiltonian can be decom-posed with the Trotter-Suzuki decomposition [62,63] or with other expansions that are unitary up to leading order, such as the Taylor series [35].
We now analyze the estimator I loc according to the same approach used in the previous section.We find that, in the limit of I → 0, the SNR of the estimator scales as This means that, as the optimization approaches the optimum of I = 0, the number of samples needed to resolve the infidelity increases as I −1 (see Appendix E for 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 [64]: where c ∈ R. The quantity c can be chosen such that the Var χ [I CV loc ] 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 |Ψ θ ⟩ → U |Ψ θ ⟩, c * is exactly −1/2.Therefore, we avoid the high cost of estimating c * at each iteration of an optimization, and directly use the asymptotically ideal estimator Eq. ( 19) with c = −1/2.
To further support this approach, we show in Fig. 3(a) that the CV estimator I CV loc features a variance that is orders of magnitude smaller than the one of the bare estimator Re I 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 I changes, such that for I → 0 we have that: 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 Re I loc goes to zero when I decreases, while the SNR of the corrected estimator decreases at smaller values of I.When c = −1/2, I CV loc has a rescaled SNR which is constant and larger than 1 over the whole range of I considered, suggesting that our strategy is ideal.
Moreover, the reduced variance of the CV estimator implies that the infidelity gradient computed with it results to be more accurate [65].The lower-variance gradient can improve the accuracy of the solution, since its mean value is affected by smaller statistical fluctuations, and can increase the speed of convergence, as it allows for larger learning rates.
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 Appendix 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.

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 onedimension, integrable or evolving via efficiently simulable [66] 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 δt.At each t, two operations are performed on the system: • unitary evolution with U = e −iHTFIδt , where H TFI is the 2D TFI Hamiltonian of Eq. ( 15); • a projective measurement of each spin in the σ 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 −iH TFI δt = e −iHzz δt 2 e −iHxδt e −iHzz δt 2 + O(δt 3 ), (21) where H zz = −J ⟨i,j⟩ σ z i σ z j and H x = −h i σ x i .In the first stages of this work we employed the forward-backward scheme as done in [36,67], but then we move to the Trotterization as it allowed to use larger δt and was more practical for our calculations.
We use the p-tVMC to apply each unitary obtained from the decomposition.Still, we remark that as exp(−iH zz δt/2) is diagonal in the chosen σ z basis, its p-tVMC optimization problem can be solved analytically, as shown in Appendix G. Instead, the unitary containing 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 H TFI , which is in correspondence of the transverse field h c such that h c /J ≈ 3.044 [68,69,70,71,72].In particular, the initial state is the paramagnetic ground state in the limit h → ∞, given by N i=1 |+⟩ i , and this is evolved into the ferromagnetic phase where h < h c .The measurement operators are Measuring spin i with outcome ↓ translates into setting ϕ ↑ i = 0 and ϕ ↓ i ̸ = 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 |Ψ θ ⟩, we monitor the Rényi-2 entanglement entropy S 2 (ρ A ) = − log 2 Tr ρ 2 A , where ρ 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 [73] as: where σ, σ ′ ∈ A and η, η ′ ∈ B (complementary of A).
See Appendix H for a derivation of Eq. ( 23).Fig. 4 shows the evolution of S 2 (ρ A ) for subsystems of size |A| ∈ [1, ⌊N/2⌋] 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 (|A| ≲ ⌊N/4⌋) 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 [74] in the insets of Fig. 4 suggest that with p = 0.01 the steadystates 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 ⌊N/2⌋ 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 [75] valid for local Hamiltonians.

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.

Data availability
All the simulations have been performed using Netket 3 [76,77] with MPI and MPI4jax [78].The code for the p-tVMC method can be found in [79].
where, like in the main text, E loc (σ) is the local energy and O k (σ) are the log-derivatives of the ansatz.F MC k does not have a covariance form like F MC k , therefore it has, in general, larger statistical fluctuations when estimated using a finite number of samples.Moreover, since O k (σ) cannot be used to compute the first term in Eq. ( 24), its computational cost is generally higher than the standard estimator.

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 |Ψ θ ⟩ = α |↓⟩ + β |↑⟩ with parameters θ = (α, β).The evolution is generated by the Hamiltonian H = σ y .For the choice α = 1 and β = 0, we have In these conditions: However, due to the covariance form of the standard tVMC estimates F MC k and S MC k,k ′ evaluated on samples including only one value of σ, we have that ∀ k, k ′ : We verify that, instead, the stochastic estimate F MC with the alternative estimator proposed in the previous section is unbiased, since: as

B.2 Two-spin system
We consider a system of N = 2 spins 1/2 whose state is represented by the variational ansatz We consider the dynamics generated by the TFI Hamiltonian H TFI with coupling J and transverse field h.In these conditions, F and S differ from F MC and S MC as: One can verify that, like for the one spin case, the alternative estimator F MC gives an unbiased estimate for the variational forces.As shown in Eq. ( 28), S MC has a lower rank than S, while F MC is identical to zero meaning that this dynamics starting from |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.

C Signal to noise ratios in Monte Carlo estimates
We consider a system of N spins 1/2 and normalized variational states |Ψ ε ⟩ whose wave function is parameterized by a single parameter ε as: and for a given σ 0 .In the following, we prove that the signal-to-noise ratios (SNRs) of F MC k and S MC kk ′ scale as O( √ ε), namely they diminish indefinitely as the wave function becomes more peaked around σ 0 .In order to ensure normalization of |Ψ ε ⟩ over different system sizes, ε goes as 1/2 N .Therefore, the two SNRs diminish exponentially as N increases, and so the number of samples needed to resolve F MC k and S MC kk ′ with finite precision is exponentially large in the system size.Instead, the SNR of the unbiased estimate F MC k of Eq. ( 24) is O(1) in ε, enabling it to efficiently estimate the forces in the limit ε → 0, where F MC k cannot be used, and independently of the system size.
Proof.The variational log-derivative of Ψ ε (σ) is: We define (1 − (2 N − 1)ε) ≡ α for brevity and we observe that α → 1 when The variance of stochastic estimates of this form is given by: For the local energy we have: where the notation H σσ ′ ≡ ⟨σ| H |σ ′ ⟩ for any σ, σ ′ is used.Therefore, keeping leading order terms in 1/ε, the variances of the stochastic estimates are: Since F MC k ≈ σ̸ =σ0 H σσ0 α/4ε and S MC kk ′ ≈ (2 N − 1)/4ε for ε ≪ 1, we have that: The two SNRs go to zero when the state becomes more peaked because when ε → 0 we have that |Ψ ε (σ)| 2 → 0 for σ ̸ = σ 0 , so these configurations are rarely sampled, but the estimators of the biases for σ ̸ = σ 0 increase as ε is reduced.Indeed, for σ ̸ = σ 0 we have that: It is possible to verify that using the unbiased force estimate F MC k of Eq. ( 24) the SNR remains constant in the limit ε → 0, making it able to efficiently compute the variational forces when the standard estimate F MC k 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 ⟨Ψ θ |σ⟩ without excluding any point.Instead, to obtain the estimator of the first term in the expression of F MC k it is necessary to divide by |Ψ θ (σ)| 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 F MC k , since this is the problematic one in the standard estimate and the second term is common to F MC k , we obtain that: Therefore, since

D Example of vanishing SNRs in tVMC
We consider a chain of N spins 1/2, initially in the ground state of the Hamiltonian − i σ x i , which evolves according to the time-dependent Hamiltonian: where γ(t) oscillates between 0 and 1 with a triangular profile of period T , namely γ(t) = t/T if 0 < t < T and γ(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 ε-close to N i=1 |↓⟩ i , so an instance of the peaked states |Ψ ε ⟩ of the previous section for some ε 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 ≈ 10 3 ), hinting at the exponentially small SNRs of F MC and S 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.

⟨η| U
where I loc (σ, η) is the term in brackets and χ(σ, η) Similarly, the (conjugate) gradient of the infidelity can be computed stochastically as: where O k (σ) = ∂ θk log Ψ θ (σ).We remark that the first term in the previous expression is free from the problem present in F MC k and S MC k when |Ψ θ ⟩ → U |Ψ θ ⟩.Indeed, in this limit we have that for σ where Ψ θ (σ) = 0 the term ⟨∂ θk Ψ θ |σ⟩ ⟨σ| U |Ψ θ ⟩ vanishes, so the bias as in Eq. ( 8) disappears.Since we consider continuous time dynamics generated by a succession of infinitesimal Us, such that |Ψ θ ⟩ ≈ U |Ψ θ ⟩ 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 I loc can be directly computed from the one of F loc = 1 − I loc as: This proves that I loc has a variance bounded above by 1 (because 0 ≤ I ≤ 1), as already noted in Ref. [80] of the main text, and that it features a zero variance principle, since Var χ [I loc (σ, η)] → 0 when the solution is approached (I → 0), as the energy in VMC for ground state search.However, the SNR of I loc is: which vanishes when I → 0. This entails that when approaching the solution the number of samples N s must diverge as 1/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 I loc .Since I is real, Re I loc can be considered in place of I 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 |F loc | 2 satisfies the required properties, since it is correlated to I loc and we have Therefore, we can define the CV infidelity estimator: where c ∈ R. I CV loc has the same mean as Re I loc for any c, while its variance depends on c.Therefore, c can be chosen such that Var χ [I CV loc (σ, η)] is minimized.This optimal value, say c * , is: where Cov χ [A(σ), B(σ)] is the covariance of functions A and B of random variable σ with distribution χ.The corresponding minimized variance is: For reasons similar to the ones previously discussed for the first term in Eq. ( 41), the estimation of the CV factor |1 − I loc | 2 is not affected by a bias or a vanishing SNR as F MC k and S MC kk ′ when |Ψ θ ⟩ → U |Ψ θ ⟩, and it is precisely in this limit where CV acts to keep the SNR constant.Therefore, I CV loc is a well-defined infidelity estimator.

F Infidelity estimator with importance sampling
The estimator Re I loc (σ, η) is unbounded because it contains ratios of wave function amplitudes.Therefore, it may diverge if ⟨σ| U |Ψ θ ⟩ and ⟨σ|Ψ θ ⟩ or ⟨η| U † |Ψ θ ⟩ and ⟨η|Ψ θ ⟩ 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 |Ψ θ ⟩ is close to zero but U † |Ψ θ ⟩ is not.Therefore, despite these configurations are rarely sampled, the estimator I loc (σ, η) 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. [64] of the main text), namely by rewriting: for a given distribution χ ′ .This latter can be chosen such that it minimizes the variance of the estimator, leading to the expression: where w(σ, η) = E χ [|I loc (σ, η)|]/|I loc (σ, η)|.Using I CV loc in place of I loc , importance sampling can be combined with control variates obtaining the estimator: 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 U acts non-trivially.

G Exact application of a diagonal propagator
The exact solution θ of the general optimization problem in the p-tVMC must satisfy: for all configurations σ and any constant C equal for all σ.For certain U and variational states, Eq. ( 50) can be solved exactly.In particular, if U is diagonal in the basis {|σ⟩} and if it is possible to write ⟨σ|Ψ θ ⟩ = ⟨σ|Ψ θ ⟩ ⟨σ|Ψ δθ ⟩ ∀σ with some update δθ, Eq. ( 50) is reduced to: where U σ is the eigenvalue of U for the eigenstate |σ⟩.In the following, we consider spin systems for simplicity.For the RBM, Eq. ( 51) admits a solution when U = R z i (ϕ z ) = e −iϕzσ z i consisting in an update of the visible bias only.If U = R zz ij (ϕ zz ) = e −iϕzzσ z i σ z j , 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 |Ψ θ ⟩, it is possible to add two variational terms such that both R z i and R zz ij can be exactly simulated with a parameter change.Indeed, |Ψ θ ⟩ can be modified into a new ansatz |Ψ J (1) i ,J (2) ij ,θ ⟩ with two additional parameters {J (1) i , J (2) ij } which is defined as: The two exponential terms factorize as required to have Eq.( 51), such that the application of R z i (ϕ z ) translates into the updates δJ

4 H− 2 εFigure 2 :
Figure 2: (a-b) Dynamics of the z-magnetization ⟨σ z i ⟩ for: (a) the N = 1 spin 1/2 toy-model, where the initial state |+⟩ is rotated with σ y , which is simulated by Exact Diagonalization (ED), tVMC and p-tVMC; (b) a system of N = 4 spins 1/2, initialized in |Ψε⟩ with ε = 1.5 • 10 −4 and evolved via H TFI (J = h = 1), which is simulated by ED and tVMC with increasing number of samples Ns (see colorbar).For the tVMC the time-step δt = 10 −3 has been used, while for the p-tVMC δt = 10 −2 .(c) Infidelity I among states evolved with tVMC and with ED after a time t f starting from |Ψε⟩ with increasing Ns (see colorbar) and different ε.The inset shows how the minimal Ns to reach I = 5•10 −3 (indicated by a dashed line in (c)) scales with ε (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 α = 1.(d) Illustration of the Born distribution for the peaked states |Ψε⟩.

5 Figure 3 :
Figure 3: (a) Learning curves for the infidelity I using the bare estimator Re I loc and the CV estimator I CV loc .In the inset the corresponding variances, generically indicated as Varχ[I], are shown.(b) Rescaled signal to noise ratio of the infidelity SNRχ[I] √ Ns as a function of I for the bare and the CV estimators with different values of c.We remark that the slope of the curves of the SNR with c ̸ = −0.5 in the limit I → 0 is approximately 1, which is better than what is given by Eq. (43).This is because the estimator used in practical calculations is Re I loc instead of simply I loc (further discussion in Appendix E).In both (a) and (b) a system of N = 8 spins 1/2 is considered and the transformation used is U = exp −i δt 2 Hzz , where Hzz = −J ⟨i,j⟩ σ z i σ z j with J = 1 and δt = 5 • 10 −2 .In (b) Ns = 10 4 samples have been employed.

2 Figure 4 :
Figure 4: Time evolution of the Rényi-2 entropy S2 simulated with p-tVMC for H TFI interspersed with random local measurements along z in 2D L × L lattices.For L = 4 and L = 5 we also provide ED benchmark results.Subsystems of increasing size |A| up to the maximum ⌊N/2⌋ are considered, and the corresponding markers are indicated with different colors according to the colorbar.The insets show the scaling of S2 as a function of |A| in the steady-state for the three lattices.S2 for subsystems with equal boundary length is indicated using the same marker.The initial state is N i=1 |+⟩ i and the parameters of H TFI are J = 1/2 and h = hc/4.The measurement rate is p = 0.01 and the time interval is δt = 0.1.The results are averaged over 5 trajectories.The ansatz is an RBM with α = 4, endowed with the variational terms to exactly implement the diagonal part of the propagator and the measurements.The number of samples is Ns = 10 4 for L = 4, and Ns = 2 • 10 4 for L = 5, 6.
ij = 0 and δθ = 0, while for R zz ij (ϕ zz ) the changes δJ ij = ϕ zz and δθ = 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 H TFI exactly.