# Looking elsewhere: improving variational Monte Carlo gradients by importance sampling

Antoine Misery, Luca Gravina, Alessandro Santini, Filippo Vicentini

► To cite this version:

Antoine Misery, Luca Gravina, Alessandro Santini, Filippo Vicentini. Looking elsewhere: improving variational Monte Carlo gradients by importance sampling. Machine Learning: Science and Technology, In press, 10.1088/2632-2153/ae387f hal-05484417

HAL Id: hal-05484417

https://hal.science/hal-05484417v1

Submitted on 29 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Looking elsewhere: improving variational Monte Carlo gradients by importance
sampling

Antoine Misery,¹,²,³,⁴,\* Luca Gravina,⁵,⁶ Alessandro Santini,¹,²,³,⁴ and Filippo Vicentini¹,²,³,⁴

¹CPHT, CNRS, Ecole Polytechnique, Institut Polytechnique de Paris, 91120 Palaiseau, France

²Collège de France, Université PSL, 11 place Marcelin Berthelot, 75005 Paris, France

³Inria Paris-Saclay, Bâtiment Alan Turing, 1, rue Honoré d'Estienne d'Orves - 91120 Palaiseau

⁴LIX, CNRS, Ecole polytechnique, Institut Polytechnique de Paris, 91120 Palaiseau, France

⁵Institute of Physics, Ecole Polytechnique Fédérale de Lausanne (EPFL), CH-1015 Lausanne, Switzerland

⁶Center for Quantum Science and Engineering, Ecole Polytechnique
Fédérale de Lausanne (EPFL), CH-1015 Lausanne, Switzerland

(Dated: January 29, 2026)

Neural-network quantum states (NQS) offer a powerful and expressive ansatz for representing quantum many-body wave functions. However, their training via Variational Monte Carlo (VMC) methods remains challenging. It is well known that some scenarios - such as sharply peaked wave functions emerging in quantum chemistry - lead to high-variance gradient estimators hindering the effectiveness of variational optimizations. In this work we investigate a systematic strategy to tackle those sampling issues by means of adaptively tuned importance sampling. Our approach is explicitly designed to target the gradient estimator instead of the loss function, and be computationally inexpensive. We benchmarked our approach across the ground-state search of a wide variety of hamiltonians, including frustrated spin systems and ab-initio quantum chemistry. Overall, our approach can reduce the computational cost of vanilla VMC considerably, up to a factor of 100x when targeting highly peaked quantum chemistry wavefunctions.

# I. INTRODUCTION

Variational Monte Carlo (VMC) algorithms are a powerful tool for managing the otherwise exponential complexity of quantum many-body systems. These methods rely on two key components: (i) a variational ansatz that effectively compresses the exponentially large wavefunction into a tractable set of polynomially many parameters, and (ii) Monte Carlo sampling, which reduces the computational cost of evaluating observables, loss functions, and their gradients. In modern applications, VMC algorithms are often combined with Neural-Network (NN) ansatze for the wavefunction, commonly referred to as Neural Quantum States (NQS) [1]. The promise of NNs lies in their ability to efficiently represent both area-law and volume-law entangled states [2] in arbitrary geometries. Furthermore, networks can be engineered to encode lattice symmetries [3–5], bosonic and fermionic exchange statistics [6–8], as well as gauge symmetries in field theories [9–11] and non-equilibrium density matrices [12–17]. The expressivity of NNs can be tuned by changing the number of parameters [18, 19]. These parameters are typically optimized by following a stochastic estimate of the natural gradient of the objective function [20–22].

While the expressiveness of variational ansatze has significantly advanced since the introduction of NQS, the Monte Carlo estimators used for estimating gradients have remained largely unchanged since the early days of VMC. The conventional approach draws samples from

the Born probability distribution $|\psi_θ(x)|^2$ and employs these samples to empirically estimate expectation values of objective functions (including energy and fidelity), their gradients, and the metric tensor. Although this sampling strategy might be justified for estimating observables, its use in computing gradients of objectives and metric tensors has been challenged. Reference [23] demonstrated that such an approach can introduce systematic bias and lead to exploding-variance issues that can severely impair VMC performance. More recently, Ref. [24] showed that the stability and asymptotic accuracy of natural gradient descent are highly sensitive to the amount of MC noise of stochastic gradient estimates.

## a. Sampling issues in ground-state calculations

Sampling-related challenges are well-documented in various spin models, such as the $J_1-J_2$ model at high frustration [18], or the XXZ model in the anti-ferromagnetic phase [25]. Peaked wave functions emerging from 2nd quantized molecular Hamiltonians also lead to the failure of VMC and require the use of techniques such as non-stochastic sampling on a preselected subspace or sampling without replacement [26–30]. The problem is often twofold. First, the structure of the wavefunction can hinder the performance of Markov Chain Monte Carlo (MCMC) sampling, leading to poor mixing times and thus a smaller effective number of samples. Second, the intrinsic variance of the gradient estimated with the Born distribution can be prohibitively large. These two issues can be addressed independently. While for the first one, autoregressive sampling strategies [2, 31–33] or improved MCMC algorithms [34, 35] have been developed, the second one remains largely unexplored and will be the focus of this manuscript.

The brute-force approach to increase the signal to noise

<footnote>\* Correspondence to: antoine.misery@polytechnique.edu</footnote>
ratio is to increase the number of samples, but this is far from optimal. In fact, modern formulations of natural gradient descent—such as NTK or minSR—require inverting a matrix whose size scales quadratically with the number of samples, imposing a practical upper limit of roughly $2^{14}$ samples unless further approximations are employed. We demonstrate that the variance of empirical gradient estimates strongly depends on the properties of the distribution from which samples are drawn. Crucially, this variance is an inherent property of the estimator, independent of the specific sampling strategy used. We show that sampling from the Born distribution is typically suboptimal and that the variance can be largely reduced by adaptive importance sampling strategies at no additional cost.

b. **Importance sampling** A commonly used method to reduce the variance of Monte Carlo estimators is importance sampling (IS), where samples are drawn from an alternative distribution $q$, typically chosen as a function of the quantity being estimated. In deep learning, IS has been applied to accelerate training by emphasizing the most informative data points [36–38], thereby effectively reducing gradient variance. In variational inference, it has also been used both to reduce variance [39, 40] and to facilitate the sampling of energy-based models [41]. In the VMC literature, a few authors have proposed sampling from $|\psi(x)|$, a simple form of IS, to enhance training performance [18, 42]. However, this approach is rarely adopted, remains poorly understood, and lacks a clear theoretical justification. Nevertheless, this sampling scheme can be interpreted as a special case of overdispersed importance sampling [39], where samples are drawn from $|\psi|^\alpha$. For $0 < \alpha < 2$ this distribution has heavier tails than $|\psi|^2$, mitigating sampling issues while introducing only a single new hyper-parameter.

In this manuscript, we systematically investigate IS as a variance reduction technique for optimizing neural quantum states. Inspired by [39], we restrict to a one-parameter family of overdispersed sampling distributions and propose an adaptive algorithm to tune this additional parameter for improved gradient estimation. The manuscript is structured as follows: in section II, we review the VMC formalism under arbitrary sampling distributions. In section IID, we present our adaptive strategy for optimizing the sampling distribution. Numerical results are reported for a range of benchmark problems: frustrated spin systems on the square lattice (section IIIA), molecular Hamiltonians (section IIIB).

## II. IMPORTANCE SAMPLING IN VARIATIONAL MONTE CARLO

In the following, we consider systems composed of $N$ particles. We denote the Hilbert space of such systems by $\mathcal{H}$ and a generic variational state by $|\psi_\theta\rangle \in \mathcal{H}$, where $\theta \in \mathbb{R}^{N_p}$ is the vector of variational parameters. We work in the computational basis $\{|x\rangle\}$ of $\mathcal{H}$. The wavefunction

amplitude for a given basis state $x$ is given by $\psi_\theta(x) = \langle x|\psi_\theta\rangle$.

### A. Importance sampling for gradient estimation

The prototypical VMC framework focuses on minimizing objective functions expressed as statistical averages over the Born distribution $p(x) = |\psi_\theta(x)|^2$ with normalization constant $Z_\theta = \langle\psi_\theta|\psi_\theta\rangle = \sum_x p(x)$. The general loss function takes the form $\mathcal{L}(\theta) = \mathbb{E}_p[\ell(x)]$, for some local quantity $\ell(x)$. For example, in ground-state optimization problems $\ell(x) = H_{\text{loc}}(x) = \langle x|H|\psi_\theta\rangle / \langle x|\psi_\theta\rangle$ is the local energy of some Hamiltonian $H$. In state compression tasks, the same form is used with $H = |\phi\rangle\langle\phi| / \langle\phi|\phi\rangle$ and for a fixed target state $|\phi\rangle$.

The optimization of $\mathcal{L}$ is driven to convergence by consecutive estimates of the gradient $\hat{F} = \nabla_\theta \mathcal{L}(\theta) \in \mathbb{R}^{N_p}$, which can itself be written as a statistical average: $F_i = \mathbb{E}_p[f_i(x)] \in \mathbb{R}$. In standard NQS implementations such as those mentioned above, the local gradient components $f_i(x)$ are given by

$$f_i(x) = 2 \operatorname{Re}\{\partial_{\theta_i} \log \psi_{\theta}(x) \Delta \ell(x)^*\}, \quad (1)$$

where $\Delta A(x) \equiv A(x) - \mathbb{E}_p[A]$ denotes a centered observable. Practical finite-sample estimates of the true gradient $\hat{F}$ are obtained by drawing $N_s$ i.i.d. samples from $p(x)$ and forming the Monte Carlo estimator

$$\hat{F}_i = \frac{1}{N_s} \sum_{\mu=1}^{N_s} f_i(x_\mu), \quad x_\mu \sim p \quad (2)$$

such that $\mathbb{V}[\hat{F}_i] = \mathbb{V}[F_i]/N_s$. To ensure reliable convergence, it is crucial to accurately estimate the direction of steepest descent $\hat{F}$ from the background noise introduced by finite Monte Carlo sampling. Our ability to do so can be quantified in terms of the signal-to-noise ratio (SNR) vector, defined as $\text{SNR}_p(\hat{F}) = (\text{SNR}_p(f_1), \dots, \text{SNR}_p(f_{N_p}))$ with

$$\text{SNR}_p(f_i) = \sqrt{\frac{|\mathbb{E}_p[f_i(x)]|^2}{\mathbb{V}_p[f_i(x)]}}. \quad (3)$$

Each parameter direction can thus be resolved with varying degrees of statistical accuracy. When $\text{SNR}_p(\hat{F}_i) = \sqrt{N_s}$, $\text{SNR}_p(f_i) < 1$, noise dominates the signal, rendering accurate estimation of the gradient difficult or impossible, thereby hindering reliable progress in the optimization. As a normalized, dimensionless quantity, the SNR provides more reliable insight than the raw variance, particularly when dealing with vector estimators as it is insensitive to the scale of each component. This notion has already appeared in the NQS literature, mainly in the study of dynamics [23, 24, 43].

To improve the resolution with which the gradient is resolved, we use importance sampling (IS): we sample
from an alternative distribution $q(x)$ that is more correlated with the quantity of interest, thereby reweighting the estimator under this new distribution. Specifically, we employ the identity

$$F_i = \mathbb{E}_p[f_i(x)] = \mathbb{E}_q[W(x)f_i(x)], \quad (4)$$

where $q(x)$ is a generic IS distribution with normalization $Z_q = \sum_x q(x)$, and

$$W(x) = \frac{|\psi_{\theta}(x)|^2 Z_q}{q(x) Z_{\theta}} \quad (5)$$

are the importance weights. When $q(x) = |\psi_{\theta}(x)|^2$, the weights reduce to $W(x) = 1$ and we recover the standard estimator. The finite-sample estimator of $F_i$ under $q(x)$ is

$$\hat{F}_i = \frac{1}{N_s} \sum_{\mu=1}^{N_s} f_i^q(x_\mu), \quad f_i^q(x) = W(x)f_i(x), \quad (6)$$

where $x_\mu \sim q$. When sampling from $q$, any centered quantity is considered centered with respect to the sampling distribution. Specifically,

$$\Delta A(x) = A(x) - \mathbb{E}_p[A(x)] = A(x) - \mathbb{E}_q[W(x)A(x)]. \quad (7)$$

Accordingly, this affects the definition of $f_i$ in eq. (1). Different choices of $q(x)$ produce estimators with different variance properties and thus different SNR vectors $\text{SNR}_q(\mathbf{F}) = (\text{SNR}_q(f_1^q), \dots, \text{SNR}_q(f_{N_p}^q))$. We evaluate the performance of a given distribution $q$ using the average SNR across all parameter directions:

$$\mathcal{L}_{\text{IS}}(q) = \frac{1}{N_p} \sum_{i=1}^{N_p} \text{SNR}_q(f_i^q). \quad (8)$$

Importantly, $\mathcal{L}_{\text{IS}}$ can be expressed as a ratio of expectations under $q$, which in turn can be estimated by Monte Carlo sampling. Specifically, the variance $\mathbb{V}_q[f_i^q(x)]$ appearing at the denominator of Eq. (3) can be expressed as

$$\mathbb{V}_q[f_i^q(x)] = \mathbb{E}_q [W^2(x) |f_i(x) - F_i|^2]. \quad (9)$$

A detailed derivation of eq. (9) is provided in appendix B 2. A practical consequence of this is that $\mathcal{L}_{\text{IS}}$ can be evaluated during the optimization at no additional cost, provided the full Jacobian has already been materialized. This is indeed the case for vanilla SR/NTK implementations as those described in Refs. [18, 44]. Nonetheless, we remark that the square modulus in eq. (9) prevents us from applying the layer-wise contraction trick from Ref. [45] that is known to further speedup the computation of the NTK.

Finally, we remark that in VMC the normalizing constants $Z_\theta$ and $Z_q$ are not known a priori and the ratio $Z_q/Z_\theta$ appearing in eq. (5) must itself be empirically estimated. This can be done by expressing

$W(x) = w(x)/\mathbb{E}_q[w(x)]$ with $w(x) = |\psi_\theta(x)|^2/q(x)$ and $\mathbb{E}_q[w(x)] = Z_\theta/Z_q$. This leads to the so-called self-normalized importance sampling estimator which is in general biased. For a detailed discussion, see appendix B 2.

## B. Optimal importance sampling distribution

Since the gradient $\mathbf{F}$ is a vector of statistical averages, maximal variance reduction across each component in general requires sampling from $N_p$ distinct distributions, each tailored to the local structure of the corresponding gradient direction. According to Hesterberg [46], the optimal sampling distribution that minimizes the variance of the self-normalized importance sampling estimator of a given gradient component $F_i$ is

$$q_{\text{opt}}^i(x) \propto |\psi_{\theta}(x)|^2 |f_i(x) - F_i|. \quad (10)$$

In practice, this strategy is impractical in VMC applications due to its high computational cost. It requires prior knowledge of the exact value of $F_i$, evaluation of derivatives at every Monte Carlo step, and mostly the ability to sample separately from a different distribution for each parameter $i$, which becomes prohibitively expensive in high-dimensional settings. Nevertheless, this distribution remains valuable as a theoretical benchmark, offering insight into how an ideal sampler would allocate probability mass to reduce variance most effectively.

## C. Overdispersed distributions

In this work we consider the family of sampling distributions defined as

$$\left\{ q_{\alpha}(x) = \frac{|\psi_{\theta}(x)|^{\alpha}}{Z_{\alpha}}, \quad \alpha \ge 0 \right\}, \quad (11)$$

where $Z_\alpha = \sum_x q_\alpha(x)$ is the normalization constant. When $\alpha < 2$, these distributions are known as overdispersed [39], as they flatten the tails of the Born distribution and converge to the uniform distribution in the limit $\alpha \to 0$. We will show that, in general, there exists an $\alpha$ for which the variance of the gradient estimator in Eq. (9) is lower than that obtained using standard Born sampling ($\alpha = 2$). While it may seem counterintuitive to underweight the regions where the wave function amplitude is largest, overdispersed sampling often leads to more accurate gradient estimates. There are two key reasons for this improvement. First, the optimal variational state $|\psi_{\theta^*}\rangle$ may be highly localized, causing standard sampling to explore only a narrow region of the Hilbert space, potentially too limited to capture meaningful gradient information. Second, flatter sampling distributions typically lead to improved mixing times and acceptance rates in MCMC methods, enabling more efficient exploration of the configuration space.
<figure><img src="image_2.png" /><figcaption>FIG. 1. Component-wise averaged signal-to-noise ratio (SNR), as defined in eq. (8), measured along the ground-state optimization trajectory of the N₂ molecule using a neural network backflow with 16 hidden units and two layers. The reference trajectory is computed in full summation, i.e., by evaluating all statistical averages exactly by summing over the entire Hilbert space. For each state visited during the optimization, we report the SNR of noisy energy gradient estimates obtained using four sampling strategies. First, the standard Born distribution |ψ<sub>θ</sub>(x)|<sup>2</sup> (green). Second, the per-component optimal importance sampling distribution q<sup>i</sup><sub>opt</sub> from eq. (10) (orange), which requires N<sub>p</sub> distinct sampling processes. Third the uniform mixture of optimal importance sampling distributions q<sup>IS</sup><sub>opt</sub> ∝ Σ<sub>i</sub> |ψ<sub>θ</sub>(x)|<sup>2</sup> |f<sub>i</sub>(x)| (blue). Fourth, the overdispersed parametric distribution |ψ<sub>θ</sub>(x)|<sup>α</sup> with α optimized at each optimization step to maximize the SNR (black). The horizontal axis shows the relative error of the variational energy with respect to the exact ground-state energy, rather than the iteration index, to enhance interpretability.</figcaption></figure>

In fig. 1, we report the performance of various gradient estimators along a prototypical ground state optimization trajectory. The performance is quantified using the component-averaged SNR, L<sub>IS</sub>(q), as defined in eq. (8). We compare four sampling strategies. First, sampling from the standard Born distribution |ψ<sub>θ</sub>(x)|<sup>2</sup>, whose SNR rapidly vanishes as convergence is approached. Second, independent sampling from the N<sub>p</sub> optimal distributions q<sup>i</sup><sub>opt</sub>, each tailored to a single gradient component, as described in section II B. In line with [46], this method achieves the slowest decay in SNR but is prohibitively expensive for large system sizes. Third, we consider sampling from a single distribution defined as a mixture of the optimal IS distributions, q<sup>IS</sup><sub>opt</sub> = 1/N<sub>p</sub> Σ<sub>i</sub> |ψ<sub>θ</sub>(x)|<sup>2</sup> |f<sub>i</sub>(x)| [47]. A detailed discussion of the rationale behind this choice can be found in appendix B 3. Finally, we examine sampling from the overdispersed distribution |ψ<sub>θ</sub>(x)|<sup>α</sup>, where the value of α is optimized at each step to maximize the SNR. We find that this approach performs comparably to q<sup>IS</sup><sub>opt</sub>, with both achieving nearly the same accuracy as the optimal component-wise sampling. This supports the use of |ψ<sub>θ</sub>(x)|<sup>α</sup> as a practical and effective parametric family for IS.

Note that sampling from the overdispersed distribution
introduces no additional computational overhead, but it

adds a new hyperparameter α to set, and possibly tune,
for each problem. The computational cost of VMC opti-
mizations motivates the development of an adaptive tun-
ing strategy for α.

D. Adaptive overdispersion strategy

Inspired by Ref. [39], we restrict our search space to the one-parameter family of overdispersed distributions q<sub>α</sub>(x) ∝ |ψ<sub>θ</sub>(x)|<sup>α</sup> and adaptively tune α to maximize L<sub>IS</sub>. Rather than relying on costly hyperparameter searches, we update α dynamically at each optimization step using a local gradient-based rule:

$$
\alpha' = \alpha + \eta_{\text{IS}} \partial_{\alpha} \mathcal{L}_{\text{IS}}(q_{\alpha}), \quad (12)
$$

where η<sub>IS</sub> > 0 is a fixed, conservative, learning rate. Since the optimal value of α depends on both the Hamil-
tonian and the current variational state, this approach
naturally adjusts to changes in the variational state as
the optimization progresses. The gradient of the IS ob-
jective L<sub>IS</sub>(q<sub>α</sub>) with respect to the overdispersion param-
eter α is found by recursive application of the chain rule:

$$
\partial_{\alpha} \mathcal{L}_{\text{IS}}(q_{\alpha}) = \frac{1}{N_p} \sum_{i=1}^{N_p} \partial_{\alpha} \text{SNR}_{q_{\alpha}}(f_i^{q_{\alpha}}), \quad (13)
$$

$$
\partial_{\alpha} \mathrm{SNR}_{q_{\alpha}}(f_{i}^{q_{\alpha}}) = \frac{1}{2} \mathrm{SNR}_{q_{\alpha}}(f_{i}^{q_{\alpha}}) \partial_{\alpha} \log (\mathbb{V}_{q_{\alpha}}[f_{i}^{q_{\alpha}}]),
$$

where

$$
\partial_{\alpha} \mathbb{V}_{q_{\alpha}} [f_{i}^{q_{\alpha}}(x)] = - \mathrm{Cov}_{q_{\alpha}} (\partial_{\alpha} \log q_{\alpha}(x), g_{i}(x)^{2}), \quad (14)
$$

$$
g_i(x) = W(x) |f_i(x) - F_i|.
$$

To preserve Markov chain thermalization—since the sam-
pling distribution evolves during training—we clip the
per-step update of α typically restricting it to a maxi-
mum increment of 0.01 so as not to excessively diverge
from the previous iterate. In practice, we find that a
learning rate of η<sub>IS</sub> = 0.1 yields stable and effective con-
vergence across all tested cases, enabling robust, fully
automated tuning of the overdispersion factor. We stress
that this parameter can still have an impact on the adap-
tive tuning trajectory, especially if it is set to be too
small. When it is too high, the exponent α will typically
change by as much as the maximum increment allows.

E. Stochastic reconfiguration

Successful ground-state and fidelity optimizations of-
ten require stochastic reconfiguration to effectively navi-
gate the rugged optimization landscape characteristic of
such problems [18, 20, 21, 49]. Stochastic reconfiguration
preconditions the bare gradient update ∇<sub>θ</sub> L(θ) using the
inverse of the Fubini-Study metric tensor, which encodes
<figure><img src="image_3.png" /><figcaption>FIG. 2. Scaling of relative error with sample size for a square $J_1-J_2$ model with $L = 6$, for several importance sampling strategies and couplings. In both cases, the adaptive strategy (black) outperforms the fixed overdispersion coefficient. (Left): $J_2 = 0$ (Heisenberg); (Middle): $J_2 = 0.5$; (Right): $J_2 = 1.0$. The adaptive trajectories for a single optimization per sample size are shown in the bottom panels, as a function of the relative error during training. One can see that in this specific case, $\alpha$ seems to decrease as we approach the ground state to greater precision by using more samples. We use $n_l = 4$ layers, $d_{\text{model}} = 24$ features and $h = 12$ heads. We impose only translational symmetry, disregarding other symmetries, and, consistent with previous ViT-based approaches, we do not impose a Marshall prior on the basis [48]. With this setting the network has $N_p \approx 10^4$ parameters, and is thus optimized using the tangent kernel method. The plotted variational energies are averaged over 10 runs.</figcaption></figure>

the local geometry of the Hilbert space at the current parameter configuration. Specifically, the update rule takes the form

$$ \theta_{k+1} = \theta_k - \eta_k (S + \lambda_k)^{-1} \nabla_{\theta} \mathcal{L}(\theta), \quad (15) $$

where $\eta_k$ is the learning rate and $\lambda_k$ is a form of Tikhonov regularization for $S$ [50]. The components of $S$ are given by

$$ S_{ij} = \text{Re}\{\mathbb{E}_q[W(x) \Delta\partial_{\theta_i}\log\psi_{\theta}(x)^* \Delta\partial_{\theta_j}\log\psi_{\theta}(x)]\}, \quad (16) $$

where the centering (denoted by $\Delta$) is performed with respect to the sampling distribution $q$, as described in eq. (7). For practical implementations, it is important that both the gradient and the metric tensor are computed over the same set of samples [51]. This ensures that the stochastic errors in both estimators are correlated, leading to reduced variance in the natural gradient update. Accordingly, the sampling distribution $q$ should be shared between the two estimators.

When targeting SR, one may opt for a sampling distribution that improves the SNR of (i) the quantum geometric tensor, (ii) the variational forces, (iii) or some quantity that depends on both. In section 3 of Ref. [24] it was shown that optimal estimators can be identified by focusing solely on the gradient estimator, with the geometric tensor playing only a minor role. Furthermore,

stable VMC optimizations with geometric tensors computed over as little as $\approx 100$ samples have been reported by Malyshev et al. [27]. In line with these findings, in this work we focus exclusively on maximizing the SNR of the bare gradient when comparing different sampling distributions, as discussed in the previous sections. This assumption will be validated a posteriori by the performance of our approach.

### III. APPLICATIONS AND NUMERICAL RESULTS

In this section, we present a comprehensive numerical study of the effects of overdispersed sampling and our autotuning strategy on the simulation of various physical systems, from spin systems (section III A) to molecular Hamiltonians (section III B).

#### A. Spin systems: $J_1-J_2$ and Heisenberg model

To evaluate the performance of IS in realistic VMC simulations, we consider the spin-1/2 $J_1-J_2$ model on the square lattice—a widely used benchmark for NQS. The
Hamiltonian is given by

$$H = J_1 \sum_{\langle i,j \rangle} \mathbf{S}_i \cdot \mathbf{S}_j + J_2 \sum_{\langle\langle i,j \rangle\rangle} \mathbf{S}_i \cdot \mathbf{S}_j, \quad (17)$$

where $\mathbf{S}_i = (S_i^x, S_i^y, S_i^z)$ denotes the spin-1/2 operator at site $i = 1, \dots, L^2$, and $\langle \cdot, \cdot \rangle$ and $\langle\langle \cdot, \cdot \rangle\rangle$ denote pairs of nearest-neighbor and next-nearest-neighbor sites, respectively. It has been shown that finding the ground state of this model close to the point of maximum frustration $J_2/J_1 = 0.5$ is a challenging problem for NQS. Bukov et al. [49] conjectured that the optimization landscape in this regime is highly rugged, with many local minima, making training in such a non-convex space especially difficult. Accurate gradient estimation is therefore crucial for ensuring convergence. This can be achieved either by increasing the number of samples or by employing IS techniques. Because of its proven performance on this very system [44], we choose a factored-attention Vision Transformer with the order of $10^4$ parameters as variational ansatz [52]. We impose only translational symmetry and, consistently with previous ViT-based approaches, we do not impose a Marshall prior on the basis [48]. Architectural details are provided in appendix C 1 a.

Figure 2 illustrates how IS improves performance relative to standard Born sampling. We analyze the scaling of the relative error $\Delta\epsilon$ with the number of samples $N_s$, across several values of the overdispersion coefficient $\alpha \in [0, 2]$. While all curves show similar scaling, their prefactors differ. While gains are modest for simpler frustration points ($J_2/J_1 = 0, 1$), the improvement is considerable at the challenging configuration $J_2/J_1 = 0.5$, where the relative error is reduced by nearly half an order of magnitude across all sample sizes. Importantly, our auto-tuned overdispersed IS procedure, displayed as a black curve, consistently yields the best performance across all frustration regimes and sample sizes.

## B. Molecular Hamiltonians in second quantization

The Hamiltonians of molecular systems,

$$H = \sum_{i,j,\sigma} t_{ij} c_{i,\sigma}^{\dagger} c_{j,\sigma} + \frac{1}{2} \sum_{i,j,k,l,\sigma,\sigma'} V_{ijkl} c_{i,\sigma}^{\dagger} c_{j,\sigma'}^{\dagger} c_{l,\sigma'} c_{k,\sigma}, \quad (18)$$

where $t_{ij}$ and $V_{ijkl}$ depend on the exact choice of molecule and basis, in general lead to very peaked wavefunctions, which makes them challenging for VMC-based approaches to tackle [7]. Alternatives like autoregressive models [27, 32] have been proposed to balance for the inefficiency of the MCMC sampling. Sophisticated sampling without replacement techniques like the Gumbel top-k trick [53] has been applied to enlarge the exploration of the Hilbert space, either on the whole space with autoregressive models [32] or on an evolving subset of configurations [28].

<figure><img src="image_4.png" /><figcaption>FIG. 3. (Top): scaling of the best VMC absolute error with respect to the Full Configuration Interaction (FCI) energy of Lithium Oxide in the STO-3G basis as a function of the number of samples. The dashed line denotes chemical accuracy. The VMC absolute error is computed using $2^{15}$ samples drawn from the Born distribution, with a colored area corresponding to the Monte-Carlo error. The ansatz is a neural network backflow with 2 layers and 256 hidden units, for a total of 7440 parameters. (Bottom): acceptance rate for several sampling distributions over an optimization performed with $N_s = 2^{12}$ samples.</figcaption></figure>

In this section, we investigate the impact of IS on the ground-state search for the Li$_2$O molecule whose geometry is retrieved from the PubChem [54] database, using the STO-3G basis and canonical Hartree-Fock (HF) orbitals; in that setting, $\dim(H) = 4.14 \times 10^7$. We use a neural network backflow with parameters identical to the literature [28, 55]. While it has been identified as a particularly challenging molecule for which standard MCMC struggles, we are able to improve the results by at least two orders of magnitude using automatic overdispersed sampling, surpassing chemical accuracy. In the top panel of fig. 3, we plot the scaling of the absolute energy error with respect to the Full Configuration Interaction (FCI) energy of Lithium Oxide in the STO-3G basis as a function of the number of samples. We are able to approximately reach the same accuracy as the one obtained in Ref. [55] with fixed-size selected configuration (FSSC) scheme and $2^{14}$ samples, using as little as $2^{11}$ samples, while the error of traditional Born sampling remains two orders of magnitude bigger. It is worth noting that the optimization remains particularly chal-
<figure><img src="image_5.png" /><figcaption>FIG. 4. Mean gradient Signal-to-Noise ratio, for several importance sampling distribution and systems. (Top): show the evolution of the SNR as a function of the relative error to the exact ground state. The black curve denotes the maximal SNR one can attain using overdispersed sampling distributions. (Bottom): evolution of the overdispersion coefficient which maximizes the SNR. Left: Nitrogen in the STO-3G basis ; Middle: square $J_1-J_2$, $L = 4$, $J_2/J_1 = 0.5$; Right: Heisenberg model, i.e. $J_2 = 0$.</figcaption></figure>

lenging and is also limited by the low acceptance rate of MCMC chains, which for $|\psi|^2$ approaches $O(10^{-4})$ near convergence, while overdispersed distributions keep a relatively high acceptance rate, as depicted in the bottom panel of fig. 3.

The acceptance could be further improved by designing a transition rule more fitted to the Hamiltonian, which introduces long-range correlations. Currently, the sampler proposes new configuration by exchanging two orbitals of the same spin along the edges of the graph defined by the adjacency matrix $A_{ij} = t_{ij} + \sum_k V_{kikj}$.

We found that using the SPRING optimizer [56] with $\mu = 0.9$ improved the overall stability; we refer the reader to appendix C 1 b for more details on the experimental setup.

### C. Signal-to-Noise Ratio analysis

To isolate the intrinsic benefits of IS from Monte Carlo noise, in fig. 4 we perform exact evaluations over the full Hilbert space, summing over all basis states without sampling. This allows us to compute $\mathcal{L}_{\text{IS}}$ exactly and determine the optimal overdispersion parameter at various stages of training.

We find that the effect of IS on the SNR depends strongly on the target system. In the $J_1-J_2$ model, in the case of a 4 × 4 system, significant SNR gains appear only near maximal frustration ($J_2/J_1 = 0.5$), where the Born distribution performs slightly worse than its overdispersed counterpart. By contrast, we examine the

nitrogen molecule $N_2$ in the STO-3G basis ($\text{dim}(\mathcal{H}) = 1.44 \times 10^4$), a classic strongly correlated system with a highly peaked ground state. In this regime, the extremely small values of $|\psi_\theta(x)|^2$ impedes the estimation of the gradient because some configurations with a high gradient magnitude will almost never be sampled. Our analysis confirms this effect: near convergence, $\mathcal{L}_{\text{IS}}(|\psi_\theta|^2)$ vanishes. Moreover, we observe at most a two-orders-of-magnitude difference between the SNR of $|\psi_\theta(x)|^2$ and the optimal *overdispersed* distribution. This difference implies that the number of required samples for a Monte-Carlo estimator under $|\psi_\theta(x)|^2$ to achieve a comparable SNR to the optimal one must be increased by at least a factor of 1000.

A practical stability criterion, proposed in Sinibaldi et al., [23], states that a Monte Carlo estimator of the gradient with $N_s$ samples becomes reliable once $\sqrt{N_s}\mathcal{L}_{\text{IS}}(q) \gtrsim 1$. This could provide a criterion for adaptive sample size; Indeed, during the VMC optimization one can dynamically adjust the sample size to maintain $\sqrt{N_s}\mathcal{L}_{\text{IS}}(q) \gtrsim 1$.

Finally, it is important to mention that VMC can still perform decently even though the SNR of the gradient is low. Indeed, stochastic sampling paired with small parameter updates at every step allow the algorithm to average out noise over many steps, gradually recovering a reliable descent direction. However, a higher SNR systematically yields greater stability, accuracy, and faster convergence in practice.
IV. CONCLUSION

In this work, we have shown that importance sampling
offers a simple yet powerful lever to reduce gradient noise
and accelerate Variational Monte Carlo (VMC) training
of neural quantum states. We have shown that even an
extremely simple one-parameter family of overdispersed
sampling distributions, $q_{\alpha}(x) \propto |\psi_{\theta}(x)|^{\alpha}$ can lead to con-
siderable improvements. Moreover, we have introduced
a technique to optimally tune the single hyperparameter
$\alpha$ at a negligible cost when performing natural-gradient
optimization.

Across several realistic VMC calculations, such as the
ground-state search in the frustrated $J_1-J_2$ Heisenberg
model, the Li₂O and N₂ molecules, we observe consis-
tent, sometimes dramatic, reductions in required sample
counts. We remark that at the beginning of this research
project we also considered more complex and costly im-
portance sampling distributions incorporating informa-
tion about the Hamiltonian or the model's Jacobian, but
found that they ultimately underperformed or were only
marginally better when compared against the much sim-
pler and cheaper overdispersed family.

Our procedure does not increase the computational
complexity of Stochastic Reconfiguration yet yields large
stability and convergence benefits. We believe it can be
readily adopted in a wide range of contexts including
impurity problems, open quantum systems, strongly cor-
related fermionic models, and non-equilibrium dynamics.

As a future research direction, we believe it might be

interesting to investigate the use of a generative model
as an IS distribution, dynamically optimized to keep the
SNR minimal while ensuring perfectly uncorrelated sam-
ples. Finally, we believe that the success of this approach
based on the SNR suggests that more about VMC con-
vergence could be understood by looking at this quantity.

The methods presented in this paper are fairly straight-
forward to add to an existing VMC codebase. The equa-
tions necessary to obtain an implementation are pro-
vided in the recapitulative appendix C 3. An implemen-
tation of the adaptive sampling algorithms discussed in
this manuscript, built on top of NetKet [57], as well as
all accompanying data can be found at the repository
GitHub@NeuralQXLab/importance_sampling_nqs.

ACKNOWLEDGMENTS

Simulations were performed with NetKet [57, 58], and
at times parallelized with mpi4JAX [59]. This software
is built on top of JAX [60] and Flax [61].

We acknowledge Ahmedeo Shokry for providing the
initial importance sampling code, as well as the Neu-
ral Network backflow code. We acknowledge insight-
ful discussions with Marylou Gabrié and Adrien Kahn,
and Markus Schmitt. F.V. acknowledges support by the
French Agence Nationale de la Recherche through the
NDQM project ANR-23-CE30-0018. This project was
provided with computing HPC and storage resources by
GENCI at IDRIS thanks to the grant 2024-A0170515698
on the supercomputer Jean Zay's V100 partition. We
used approximately 5k V100 hours and 1.5k A100 hours.

[1] G. Carleo and M. Troyer, Solving the quantum many-body problem with artificial neural networks, Science **355**, 602 (2017).

[2] O. Sharir, Y. Levine, N. Wies, G. Carleo, and A. Shashua, Deep Autoregressive Models for the Efficient Variational Simulation of Many-Body Quantum Systems, Physical Review Letters **124**, 020503 (2020).

[3] K. Choo, G. Carleo, N. Regnault, and T. Neupert, Symmetries and Many-Body Excitations with Neural-Network Quantum States, Physical Review Letters **121**, 167204 (2018).

[4] M. Reh, M. Schmitt, and M. Gärttner, Optimizing design choices for neural quantum states, Physical Review B **107**, 10.1103/physrevb.107.195115 (2023).

[5] R. P. Nutakki, A. Shokry, and F. Vicentini, Design principles of deep translationally-symmetric neural quantum states for frustrated magnets (2025).

[6] D. Luo and B. K. Clark, Backflow Transformations via Neural Networks for Quantum Many-Body Wave Functions, Physical Review Letters **122**, 226401 (2019).

[7] K. Choo, A. Mezzacapo, and G. Carleo, Fermionic neural-network states for ab-initio electronic structure, Nature Communications **11**, 2368 (2020).

[8] J. Robledo Moreno, G. Carleo, A. Georges, and J. Stokes,

Fermionic wave functions from neural-network con-
strained hidden states, Proceedings of the National
Academy of Sciences **119**, e2122059119 (2022).

[9] D. Luo, G. Carleo, B. K. Clark, and J. Stokes, Gauge equivariant neural networks for quantum lattice gauge theories, Physical Review Letters **127**, 10.1103/phys-revlett.127.276402 (2021).

[10] D. Luo, Z. Chen, K. Hu, Z. Zhao, V. M. Hur, and B. K. Clark, Gauge-invariant and anyonic-symmetric autoregressive neural network for quantum lattice models, Physical Review Research **5**, 10.1103/physrevresearch.5.013216 (2023).

[11] A. Apte, C. Córdova, T.-C. Huang, and A. Ashmore, Deep learning lattice gauge theories, Physical Review B **110**, 10.1103/physrevb.110.165133 (2024).

[12] G. Torlai and R. G. Melko, Latent space purification via neural density operators, Physical Review Letters **120**, 10.1103/physrevlett.120.240503 (2018).

[13] F. Vicentini, R. Rossi, and G. Carleo, Positive-definite parametrization of mixed quantum states with deep neural networks (2022), 2206.13488.

[14] F. Vicentini, A. Biella, N. Regnault, and C. Ciuti, Variational neural-network ansatz for steady states in open quantum systems, Physical Review Letters **122**,
10.1103/physrevlett.122.250503 (2019).

[15] D. Luo, Z. Chen, J. Carrasquilla, and B. K. Clark, Autoregressive neural network for simulating open quantum systems via a probabilistic formulation, Physical Review Letters **128**, 10.1103/physrevlett.128.090501 (2022).

[16] D. Eeltink, F. Vicentini, and V. Savona, Variational dynamics of open quantum systems in phase space (2023).

[17] M. Reh, M. Schmitt, and M. Gärttner, Time-dependent variational principle for open quantum systems with artificial neural networks, Physical Review Letters **127**, 230501 (2021).

[18] A. Chen and M. Heyl, Empowering deep neural quantum states through efficient optimization, Nature Physics **20**, 1476 (2024).

[19] S. Dash, L. Gravina, F. Vicentini, M. Ferrero, and A. Georges, Efficiency of neural quantum states in light of the quantum geometric tensor, Communications Physics **8**, 10.1038/s42005-025-02005-4 (2025).

[20] S. Sorella, Green Function Monte Carlo with Stochastic Reconfiguration, Physical Review Letters **80**, 4558 (1998).

[21] S.-i. Amari, Natural gradient works efficiently in learning, Neural Computation **10**, 251 (1998).

[22] J. Martens, New insights and perspectives on the natural gradient method, Journal of Machine Learning Research **21**, 1 (2020).

[23] A. Sinibaldi, C. Giuliani, G. Carleo, and F. Vicentini, Unbiasing time-dependent Variational Monte Carlo by projected quantum evolution, Quantum **7**, 1131 (2023).

[24] L. Gravina, V. Savona, and F. Vicentini, Neural Projected Quantum Dynamics: A systematic study, Quantum **9**, 1803 (2025).

[25] C.-Y. Park and M. J. Kastoryano, Expressive power of complex-valued restricted Boltzmann machines for solving nonstochastic Hamiltonians, Physical Review B **106**, 134437 (2022).

[26] X. Li, J.-C. Huang, G.-Z. Zhang, H.-E. Li, C.-s. Cao, D. Lv, and H.-S. Hu, A Non-stochastic Optimization Algorithm for Neural-network Quantum States, Journal of Chemical Theory and Computation **19**, 8156 (2023).

[27] A. Malyshev, M. Schmitt, and A. I. Lvovsky, Neural Quantum States and Peaked Molecular Wave Functions: Curse or Blessing? (2024), arXiv:2408.07625 [quant-ph].

[28] A.-J. Liu and B. K. Clark, Efficient optimization of neural network backflow for *ab initio* quantum chemistry, Phys. Rev. B **112**, 155162 (2025).

[29] X. Li, J.-C. Huang, G.-Z. Zhang, H.-E. Li, Z.-P. Shen, C. Zhao, J. Li, and H.-S. Hu, Improved optimization for the neural-network quantum states and tests on the chromium dimer, The Journal of Chemical Physics **160**, 234102 (2024).

[30] X. Cao, Z. Zhong, and Y. Lu, Vision transformer neural quantum states for impurity models, Phys. Rev. B **112**, 235155 (2025).

[31] M. Hibat-Allah, M. Ganahl, L. E. Hayward, R. G. Melko, and J. Carrasquilla, Recurrent neural network wave functions, Physical Review Research **2**, 023358 (2020).

[32] A. Malyshev, J. M. Arrazola, and A. I. Lvovsky, Autoregressive Neural Quantum States with Quantum Number Symmetries (2023), arXiv:2310.04166 [quant-ph].

[33] T. Zhao, J. Stokes, and S. Veerapaneni, Scalable neural quantum states architecture for quantum chemistry, Machine Learning: Science and Technology **4**, 025034 (2023).

[34] A. A. Bagrov, A. A. Iliasov, and T. Westerhout, Kinetic samplers for neural quantum states, Physical Review B **104**, 104407 (2021).

[35] S. Bravyi, G. Carleo, D. Gosset, and Y. Liu, A rapidly mixing Markov chain from any gapped quantum many-body system, Quantum **7**, 1173 (2023).

[36] G. Alain, A. Lamb, C. Sankar, A. Courville, and Y. Ben-gio, Variance Reduction in SGD by Distributed Importance Sampling (2016), arXiv:1511.06481 [stat].

[37] A. Katharopoulos and F. Fleuret, Not All Samples Are Created Equal: Deep Learning with Importance Sampling, in *Proceedings of the 35th International Conference on Machine Learning* (PMLR, 2018) pp. 2525–2534.

[38] R. Liu, T. Wu, and B. Mozafari, Adam with Bandit Sampling for Deep Learning, ArXiv (2020).

[39] F. J. R. Ruiz, M. K. Titsias, and D. M. Blei, Overdispersed Black-Box Variational Inference (2016), arXiv:1603.01140 [stat].

[40] X. Li, C. Li, J. Chi, and J. Ouyang, Variance Reduction in Black-box Variational Inference by Adaptive Importance Sampling, in *Proceedings of the Twenty-Seventh International Joint Conference on Artificial Intelligence* (International Joint Conferences on Artificial Intelligence Organization, Stockholm, Sweden, 2018) pp. 2404–2410.

[41] L. Grenioux, E. Moulines, and M. Gabrié, Balanced Training of Energy-Based Models with Adaptive Flow Sampling (2024), arXiv:2306.00684.

[42] K. Inui, Y. Kato, and Y. Motome, Determinant-free fermionic wave function using feed-forward neural networks, Physical Review Research **3**, 10.1103/physrevresearch.3.043126 (2021).

[43] M. Schmitt and M. Heyl, Quantum many-body dynamics in two dimensions with artificial neural networks, Physical Review Letters **125**, 100503 (2020).

[44] R. Rende, L. L. Viteritti, L. Bardone, F. Becca, and S. Goldt, A simple linear algebra identity to optimize large-scale neural network quantum states, Communications Physics **7**, 1 (2024).

[45] R. Novak, J. Sohl-Dickstein, and S. S. Schoenholz, Fast finite width neural tangent kernel (2022).

[46] T. Hesterberg, Weighted Average Importance Sampling and Defensive Mixture Distributions, Technometrics : a journal of statistics for the physical, chemical, and engineering sciences **37**, 185 (1995), 1269620.

[47] We do not include a mixture of the exact optimal distributions $q_i^{\text{opt}}$ as estimating each $F_i$ at every MCMC step would be computationally prohibitive and largely redundant.

[48] W. Marshall, Antiferromagnetism, Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences **232**, 48 (1955), 99682.

[49] M. Bukov, M. Schmitt, and M. Dupont, Learning the ground state of a non-stochastic quantum Hamiltonian in a rugged neural network landscape, SciPost Physics **10**, 147 (2021).

[50] J. Nocedal and S. J. Wright, *Numerical Optimization*, Springer Series in Operations Research and Financial Engineering, Vol. 2 (Springer New York, 2006).

[51] F. Becca and S. Sorella, *Quantum Monte Carlo Approaches for Correlated Systems* (Cambridge University Press, Cambridge, 2017).

[52] L. L. Viteritti, R. Rende, and F. Becca, Transformer Variational Wave Functions for Frustrated Quantum Spin Systems, Physical Review Letters **130**, 236401 (2023).
[53] W. Kool, H. van Hoof, and M. Welling, Stochastic Beams and Where to Find Them: The Gumbel-Top-k Trick for Sampling Sequences Without Replacement (2019), arXiv:1903.06059 [cs].

[54] PubChem Substance and Compound databases | Nucleic Acids Research | Oxford Academic.

[55] A.-J. Liu and B. K. Clark, Neural network backflow for ab initio quantum chemistry, Physical Review B 110, 115137 (2024).

[56] G. Goldshlager, N. Abrahamsen, and L. Lin, A Kaczmarz-inspired approach to accelerate the optimization of neural network wavefunctions, Journal of Computational Physics 516, 113351 (2024).

[57] F. Vicentini, D. Hofmann, A. Szabó, D. Wu, C. Roth, C. Giuliani, G. Pescia, J. Nys, V. Vargas-Calderon, N. Astrakhantsev, and G. Carleo, NetKet 3: Machine Learning Toolbox for Many-Body Quantum Systems, Scientific Physics Codebases, 7 (2022).

[58] G. Carleo, K. Choo, D. Hofmann, J. E. T. Smith, T. Westerhout, F. Alet, E. J. Davis, S. Efthymiou, I. Glasser, S.-H. Lin, M. Mauri, G. Mazzola, C. B. Mendl, E. van Nieuwenburg, O. O'Reilly, H. Théveniaut, G. Torlai, F. Vicentini, and A. Wietek, NetKet: A machine learning toolkit for many-body quantum systems, SoftwareX 10, 100311 (2019).

[59] D. Häfner and F. Vicentini, Mpi4jax: Zero-copy MPI communication of JAX arrays, Journal of Open Source Software 6, 3419 (2021).

[60] J. Bradbury, R. Frostig, P. Hawkins, M. J. Johnson, C. Leary, D. Maclaurin, G. Necula, A. Paszke, J. VanderPlas, S. Wanderman-Milne, and Q. Zhang, JAX: Composable transformations of Python+NumPy programs (2018).

[61] J. Heek, A. Levskaya, A. Oliver, M. Ritter, B. Ronde-pierre, A. Steiner, and M. van Zee, Flax: A neural network library and ecosystem for JAX (2024).

[62] G. W. Oehlert, A note on the delta method, The American Statistician 46, 27 (1992), 2684406.

<figure><img src="image_6.png" /><figcaption>FIG. 5. Adaptive importance sampling results on the 10 × 10 square lattice for the $J_1-J_2$ model at $J_2/J_1 = 0.5$. (Top): optimization curves sampling from the Born distribution (grey) and using adaptive importance sampling (orange), using 8192 samples and the same seed to initialize the parameters. A moving average over 10 steps is employed to smoothen the curves. The relative error is computed with respect to the exact energy extrapolation reported in [18]. The variational energy obtained by Rende et al. [44] with the same network using 6000 samples is denoted by the dashed line. (Bottom): trajectory of $\alpha$ during the optimization.</figcaption></figure>

# Appendix A: Additional results

## 1. Scaling to larger system sizes

The tests of section III A on the $J_1-J_2$ model at $J_2/J_1 = 0.5$ have been extended to a $10 \times 10$ square lattice, using the same network as Rende et al. [44], namely a Vision Transformer with 8 layers, 72 features and 12 heads, totaling more than $10^5$ parameters. In fig. 5, we compare standard sampling from the Born distribution to our adaptive importance sampling method, for an un-symmetrized network. We show both optimization curves in the top panel, proving that adaptive importance sampling is able to converge to a lower variational energy. The result is consistent with fig. 2 where we also showed a small accuracy improvement.

The trajectory of $\alpha$ is also shown in the bottom panel, which converges to $\alpha \approx 1.6$. This supports the observation that the behavior of the optimal sampling distribution strongly depends on the network, which here is a variant with more parameters of the ViT employed in

section III A.

# Appendix B: Methods

## 1. Variational Monte Carlo

As explained in the main text, Variational Monte Carlo aims at minimizing the expectation value of the hamiltonian under a variational state $|\psi_\theta\rangle$:

$$E_{\theta} = \frac{\langle \psi_{\theta} | H | \psi_{\theta} \rangle}{\langle \psi_{\theta} | \psi_{\theta} \rangle}.$$

By inserting the sum of projector along each element of the orthogonal basis $\{|x\rangle\}$, namely $\mathbb{I} = \sum_x |x\rangle\langle x|$ one can rewrite the previous expression as an expectation value under the Born amplitude $p(x) = \frac{|\psi_\theta(x)|^2}{\langle \psi_\theta | \psi_\theta \rangle}$:

$$E_{\theta} = \mathbb{E}_{p}[H_{loc}(x)],$$
where:

$$H_{loc}(x) = \frac{1}{\psi_{\theta}(x)} \sum_{y, H_{xy} \neq 0} H_{xy} \psi_{\theta}(y)$$

is called the local energy. This quantity is indeed local as it can be evaluated with the knowledge of the NQS $\psi_{\theta}(\cdot)$ and $|x\rangle$.

## 2. The self-normalized Importance Sampling estimator

In this section, we review the mathematical aspects of the self-normalized IS estimator (SNIS). We start by defining the SNIS estimator, then computes its variance and bias. Afterwards, we derive the derivative of the variance of the SNIS estimator for a variational sampling distribution.

### a. Definition

The definition of the Self-normalized IS estimator is given by eq. (B2) below. To define it, and discuss its properties we start by introducing some notation.

* $P(x) = p(x)/Z_p$ is the original (or target) Born distribution we want moments under.

* $Q(x) = q(x)/Z_q$ is the proposal distribution from which we actually intend to sample.

* $x_i \sim Q$ for $i = 1, \dots, n$ are the samples.

* $X$ and $Y$ are arbitrary functions of interest and $X_i = X(x_i)$ and $Y_i = Y(x_i)$.

* $w_i = w(x_i) = p(x_i)/q(x_i)$ are the raw importance weights and $W_1 = \sum_i w_i$, $W_2 = \sum_i w_i^2$ are the first and second weight-moments.

* $\tilde{w}_i = w_i/W_1$ are the self-normalized weights.

We note for future reference that since

$$\mathbb{E}_Q[w(x)] = \sum_x \frac{q(x)}{Z_q} \frac{p(x)}{q(x)} = \frac{\sum_x p(x)}{Z_q} = \frac{Z_p}{Z_q}, \quad (B1)$$

the expectation value under $P$ of some local quantity $f$ can be written as a ratio,

$$\mathbb{E}_P[f(x)] = \sum_x \frac{p(x)}{Z_p} f(x) = \sum_x \frac{p(x)}{Z_p} \frac{q(x)}{q(x)/Z_q} f(x) \\ = \frac{\sum_x \frac{q(x)}{Z_q} w(x) f(x)}{Z_p/Z_q} = \frac{\mathbb{E}_Q[w(x) f(x)]}{\mathbb{E}_Q[w(x)]}.$$

We now define the self-normalized importance sampling estimator as the plug-in estimator of the above ratio, that is

$$\begin{aligned} \mathbb{E}_P[f(x)] &= \frac{\mathbb{E}_Q[w(x)f(x)]}{\mathbb{E}_Q[w(x)]} \approx \frac{\frac{1}{n} \sum_{i=1}^n w_i f(x_i)}{\frac{1}{n} \sum_{i=1}^n w_i} \\ &= \frac{\sum_{i=1}^n w_i f(x_i)}{W_1} = \sum_{i=1}^n \tilde{w}_i f(x_i). \end{aligned} \quad (B2)$$

Let $X$ and $Y$ be two random variables. Let $\mu_X = \mathbb{E}_P[X]$, $\mu_Y = \mathbb{E}_P[Y]$, $\mu_{XY} = \mathbb{E}_P[XY]$, and $C_{XY} = \text{Cov}_P(X, Y) = \mathbb{E}_P[(X - \mu_X)(Y - \mu_Y)]$. Equation (B2) can be applied to the means of $X(x)$ and $Y(x)$ for which we have that

$$\mu_X = \mathbb{E}_P[X] \approx \hat{\mu}_X = \sum_{i=1}^{n} \tilde{w}_i X_i,$$

and similarly for $\mu_Y$. For the expression of the covariance of $X$ and $Y$, we have an analogous replacement:

$$\begin{aligned} C_{XY} &= \mathbb{E}_P[(X(x) - \mu_X)(Y(x) - \mu_Y)] \\ &= \frac{\mathbb{E}_Q[w(x)(X(x) - \mu_X)(Y(x) - \mu_Y)]}{\mathbb{E}_Q[w(x)]}, \end{aligned}$$

with the plug-in estimator being

$$\hat{C}_{XY}^{\text{simple}} = \sum_{i=1}^{n} \tilde{w}_i (X_i - \hat{\mu}_X)(Y_i - \hat{\mu}_Y) = \sum_{i=1}^{n} \tilde{w}_i X_i Y_i - \hat{\mu}_X \hat{\mu}_Y. \quad (B3)$$

Recall that now $x_i$ are sampled from $Q$, not from $P$. We now want to estimate the variance of this estimator. A finite-sample closed form is algebraically complicated because of the normalizing denominator. Following the overwhelming majority of the literature, we rely on the first-order multi-variate delta method [62] to obtain the asymptotic variance of both estimators.

### b. Variance of the SNIS sample covariance

In this appendix we will prove that the variance of the covariance estimator (eq. (B3)) of the self-normalized importance sampling estimator is

$$\mathbb{V}(\hat{F}_i^q) = \frac{1}{N_s} \frac{\mathbb{E}_q[w^2(x)|f_i(x) - F_i|^2]}{(\mathbb{E}_q[w(x)])^2}. \quad (B4)$$

The proof follows. We define the random vectors

$$\mathbf{V}_i = \frac{Z_q}{Z_p} \begin{pmatrix} w_i \\ w_i X_i \\ w_i Y_i \\ w_i X_i Y_i \end{pmatrix}, \quad \mathbf{\Sigma}_n = \frac{1}{n} \sum_{i=1}^n \mathbf{V}_i = \begin{pmatrix} A \\ B \\ C \\ D \end{pmatrix}$$
with

$$A = \frac{1}{n} \sum_{i} w_{i},$$

$$B = \frac{1}{n} \sum_{i} w_{i} X_{i},$$

$$C = \frac{1}{n} \sum_{i} w_{i} Y_{i},$$

$$D = \frac{1}{n} \sum_{i} w_{i} X_{i} Y_{i}.$$

Note that the factor $Z_q/Z_p$ in the definition of $V_i$ is deterministic and only serves the purpose of normalizing the expectation values. Indeed, the expectations under $Q$ are

$$m = \mathbb{E}_Q[V] = (1, \mu_X, \mu_Y, \mu_{XY})^T = \mathbb{E}_Q[\Sigma_n] \quad (B5)$$

where $V$ is a representative of the sample statistics $V_i$. Replacing into Eq. (B3) we obtain

$$\hat{C}_{XY}^{\text{simple}} = g(\Sigma_n) = g(A, B, C, D) = \frac{D}{A} - \frac{BC}{A^2}.$$

To apply the delta method [62] we need to compute the gradient of $g(\Sigma_n)$ and evaluate it at $m$:

$$\frac{\partial g}{\partial A} = -\frac{D}{A^2} + \frac{2BC}{A^3} \xrightarrow{(B5)} -C_{XY} + \mu_X \mu_Y,$$

$$\frac{\partial g}{\partial B} = -\frac{C}{A^2} \xrightarrow{(B5)} -\mu_Y$$

$$\frac{\partial g}{\partial C} = -\frac{B}{A^2} \xrightarrow{(B5)} -\mu_X$$

$$\frac{\partial g}{\partial D} = \frac{1}{A} \xrightarrow{(B5)} 1,$$

so that

$$\nabla g(A, B, C, D) \bigg|_{m} = (-C_{XY} + \mu_X \mu_Y, -\mu_Y, -\mu_X, 1)^T.$$

To conclude, we need the covariance matrix of one draw, that is

$$\Sigma = \mathrm{Cov}_q[V] = \mathbb{E}_q[VV^T] - mm^T$$

From the delta method [62], we have that to leading order in $n$

$$\mathbb{V}[\hat{C}_{XY}^{\text{simple}}] \approx \frac{1}{n} \nabla g(m)^T \Sigma \nabla g(m) + o(n^{-1}). \quad (B6)$$

From here on out we will neglect the terms vanishing faster than $1/n$. Computing the leading term by explicitly carrying out the contraction is straightforward but can be quite tedious. Instead we realize that

$$\nabla g(m)^T \Sigma \nabla g(m) = \nabla g(m)^T \mathrm{Cov}_Q[V] \nabla g(m) = \mathbb{V}_Q[\nabla g(m)^T V],$$

where we used that $\nabla g(m)^T V$ is a scalar quantity to replace the covariance with the variance. Since $\nabla g(m)^T$ is deterministic, we have that $\mathbb{E}_Q(\nabla g(m)^T V) = m$. Now

$$\nabla g(m)^T (V - m) = \frac{Z_q}{Z_p} \cdot \frac{w(XY - \mu_Y X - \mu_X Y - C_{XY})}{(X - \mu_X)(Y - \mu_Y) - C_{XY}} + \underbrace{C_{XY} - \mu_X \mu_Y + 2\mu_X \mu_Y - \mu_{XY}}_{= 0} = \frac{Z_q}{Z_p} w(C(X, Y) - C_{XY}) = \frac{Z_q}{Z_p} w H(X, Y),$$

where

$$H(X, Y) \equiv C(X, Y) - C_{XY} \equiv (X - \mu_X)(Y - \mu_Y) - C_{XY}.$$

Note that $H(X, Y)$ is a functional and should be seen as $H(X, Y)(x) = H(X(x), Y(x))$. Finally we have that

$$\begin{aligned} \mathbb{V}[\hat{C}_{XY}^{\text{simple}}] &= \frac{1}{n} \nabla g(m)^T \Sigma \nabla g(m) \\ &= \frac{1}{n} \mathbb{V}_Q \left[ \frac{Z_q}{Z_p} w H(X, Y) \right] = \frac{1}{n} \frac{\mathbb{V}_Q[w H(X, Y)]}{(Z_p/Z_q)^2} \\ &= \frac{1}{n} \frac{\mathbb{E}_Q[w^2 H^2(X, Y)]}{(\mathbb{E}_Q[w])^2}, \end{aligned}$$

where in the last equality we used that $wH$ is centered in $Q$. The plug-in estimator of this quantity is

$$\begin{aligned} \hat{\mathbb{V}}[\hat{C}_{XY}^{\text{simple}}] &= \frac{1}{n} \frac{\frac{1}{n} \sum_i w_i^2 H^2(X_i, Y_i)}{\left(\frac{1}{n} \sum_i w_i\right)^2} \\ &= \sum_i \left(\frac{w_i}{\sum_i w_i}\right)^2 H^2(X_i, Y_i) = \sum_i \tilde{w}_i^2 H^2(X_i, Y_i) \\ &= \sum_i \tilde{w}_i^2 [(X_i - \hat{\mu}_X)(Y_i - \hat{\mu}_Y) - \hat{C}_{XY}^{\text{simple}}]^2 \end{aligned}$$

Going back to the notation of the paper, we see that we get eq. (B4).

### c. Bias

In this appendix we prove that the bias of the SNIS estimator is

$$\mu_X \left[ \frac{1}{n} \left( \underbrace{\frac{\mathbb{E}_Q[w^2(\mu_X - X)]}{\mu_X \mathbb{E}_Q[w]^2}}_{\rho_0} \right) \right]. \quad (B7)$$

The proof follows. To simplify the calculation, we will make an approximation by neglecting the covariance structure of the estimator, and only the consider the bias related to the ratio structure. This approximation is motivated by the fact that the covariance will contribute higher order terms which we assume are not dominant here (we also remark that without this approximation
it would be very hard to get interpretable analytical results).

We recall that the SNIS estimator is denoted as $F = \frac{1}{n} \sum_{i=1}^{n} \tilde{w}_i X_i$. Similarly to the previous section, we will use the delta method to compute the bias. In that case, we need to go to second order as $\mathbb{E}_Q(\frac{Z_q}{Z_p}(wX) - \mu_X) = \mathbb{E}_Q[\frac{Z_q}{Z_p}w - 1] = 0$.

Defining:

$$V_i = \frac{Z_q}{Z_p} \begin{pmatrix} w_i \\ w_i X_i \end{pmatrix}, \quad \Sigma_n = \frac{1}{n} \sum_{i=1}^{n} V_i = \begin{pmatrix} A \\ B \end{pmatrix},$$

one observes that given $g(A, B) = \frac{B}{A}$ and $\mathbf{m} = (\mu_X, 1)^T$, we have:

$$\nabla^2 g(\mathbf{m}) = \begin{pmatrix} 2\mu_X & -1 \\ \cdot & 0 \end{pmatrix}.$$

Let $r = \frac{Z_p}{Z_q}$. One can thus write:

$$\begin{aligned} \mathbb{E}(F) &\approx \mu_X + \frac{1}{2} \left[ 2\mu_X \frac{1}{n} \mathbb{V}_Q\left(\frac{w}{r}\right) - \frac{2}{n} \mathrm{Cov}_Q\left(\frac{w}{r}, \frac{wX}{r}\right) \right] \\ &= \mu_X \left[ 1 + \frac{1}{n} \left( \underbrace{\frac{\mathbb{E}_Q[w^2(\mu_X - X)]}{\mu_X r^2}}_{\rho_0} \right) \right], \end{aligned}$$

where we denoted as $\rho_0$ the shrinking factor, independent from the sample size $n$. One has especially:

$$\frac{\mathbb{V}_Q(w)}{r^2} = \frac{1}{\text{ESS}} - 1,$$

where the effective sample size (ESS), defined as $\frac{\mathbb{E}_Q[w(x)]^2}{\mathbb{E}_Q[w(x)]^2}$, is a common metric to measure the similarity between two probability distributions.

It is often employed in the usual context of IS where one tries to approximate an intractable distribution as a mean to sample from it. It varies from 0 for distributions with disjoint support, to 1 for identical distributions.

As a consequence, one can infer that the more dissimilar $q$ and $p$ are, the greater the bias will be.

In fig. 6 we plot the signal-to-noise ratio as well as the shrinkage factor $\rho_0$ defined above for a variational state taken at the end of the ground state search for the nitrogen, as described in section III C. One sees that in that case, $\rho_0 \sim O(10^1)$ which is negligible compared to the $\frac{1}{n}$ prefactor as $n$ is usually taken in the range $10^3 - 10^4$. However, $\rho_0$ quickly increases as $\alpha$ approaches 0, which suggests that using a distribution that is too different from the original sampling distribution may lead to problems, except if one find a way to mitigate this bias by employing a refined estimator.

### d. Gradient of the SNR

Let $(q_\alpha)_{\alpha \in \mathbb{R}^m}$ be a variational family in $\mathcal{F}(\{0, 1\}^N, \mathbb{R}^+)$. Let $Q_\alpha(x) = \frac{q_\alpha(x)}{Z_\alpha}$ the associated

<figure><img src="image_7.png" /><figcaption>FIG. 6. Effective sample size and SNR curves over $\alpha$ computed at the end of a ground state search of the $N_2$ molecule, one of the cases where the optimal $\alpha$ drops to almost 0.5. The shrinkage factor is relatively small compared to the usually used sample sizes of $N_s \sim O(10^3 - 10^4)$</figcaption></figure>

normalized probability distribution. Differentiating $Q_\alpha$ with respect to $\alpha$ gives:

$$\forall i \in \{1...m\}, \frac{\partial_{\alpha_i} q_{\alpha}(x)}{Z_{\alpha}} - \sum_{y} \frac{\partial_{\alpha_i} q_{\alpha}(y)}{Z_{\alpha}} \frac{Q_{\alpha}(x)}{Z_{\alpha}},$$

which can in general be simplified as:

$$\begin{aligned} \partial_{\alpha_i} Q_{\alpha}(x) &= Q_{\alpha}(x)(\partial_{\alpha} \log q_{\alpha}(x) - \mathbb{E}_{q_{\alpha}}[\partial_{\alpha_i} \log q_{\alpha}(x)]) \\ &= q_{\alpha}(x) \Delta_{q_{\alpha}} \partial_{\alpha} \log q_{\alpha}(x), \end{aligned}$$

where

$$\Delta_{q_{\alpha}} \partial_{\alpha} \log q_{\alpha}(x) = \partial_{\alpha} \log q_{\alpha}(x) - \mathbb{E}_{q_{\alpha}}[\partial_{\alpha} \log q_{\alpha}(x)].$$

Thus one can write, taking the notations of the main text:

$$\begin{aligned} \partial_{\alpha_i} \mathbb{V}[\hat{F}_i] &= \partial_{\alpha_i} \left( \sum_x \frac{P(x)^2}{Q_\alpha(x)} |f_i(x) - F_i|^2 \right) \\ &= -\sum_x Q_\alpha(x) \Delta_q \partial_{\alpha_i} \log q_\alpha(x) g_i^2(x) \\ &= -\mathbb{E}_{q_\alpha} [\Delta_{q_\alpha} \partial_{\alpha_i} \log q_\alpha(x) g_i^2(x)], \end{aligned}$$

where $g_i(x) = W(x)|f_i(x) - F_i|$. The gradient of $\mathcal{L}_{IS}$ with respect to $\alpha$ can thus be obtained through the chain rule.

### 3. Alternative definitions of the objective function

As mentioned in the main text, one knows analytically that the distribution minimizing the variance of each gradient component is:

$$q_{\text{opt}}^{i}(x) \propto |\psi(x)|^{2}(x) |f_{i}(x) - F_{i}|.$$

We will denote $q_{\text{opt}}^{\text{SNIS}}(x) = \frac{1}{N_p} \sum_{i=1}^{N_p} q_{\text{opt}}^{i}(x)$ the average of the above distributions, and $q_\alpha$ a variational family of
distribution. One can choose a measure of divergence between probability distribution as the objective function, like the Kullback-Leibler divergence [40]:

$$D_{\text{KL}}(q_{\alpha}|q_{\text{opt}}^{\text{SNIS}}) = \mathbb{E}_{q_{\alpha}}\left[\log\left(\frac{q_{\alpha}(x)}{q_{\text{opt}}^{\text{SNIS}}(x)}\right)\right].$$

One could also consider the following distribution, which is the average of the optimal distributions for standard IS:

$$q_{\text{opt}}^{\text{IS}}(x) \propto \sum_{i=1}^{N_p} |\psi(x)|^2 |f_i(x)|.$$

In fig. 7, we show that the minimums of KL divergences almost coincides with the maximum SNR for $q_{\text{opt}}^{\text{IS}}$ but not $q_{\text{opt}}^{\text{SNIS}}$. Indeed, averaging all probability distribution may break optimality in a more dramatic way for $q_{\text{opt}}^{\text{SNIS}}$.

<figure><img src="image_8.png" /><figcaption>FIG. 7. Comparison of several objective for the adaptive tuning. We compare Signal-to-noise ratio (grey) with Kullback Leibler divergence with respect to near optimal distributions for importance sampling (orange and blue). The comparison is carried on the state obtained at the end of the ground state search for the nitrogen molecule.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>FIG. 8. Comparison of SPRING with $\mu = 0.9$ vs SGD SR using two different sample sizes. SPRING improves the optimization at low sample sizes.</figcaption></figure>

upper index runs on the number of heads. These vectors are then passed through a series of multi-head factored attention layers, whose central part is the following operation:

$$\alpha_i^\mu = \sum_{j=1}^{N} A_{i,j}^\mu X_j^\mu,$$

where $N$ denotes the number of patches, and $A_{i,j}^\mu$ is a tensor of variational parameters.

Eq. C1 a essentially mixes the patches’ embeddings, allowing the model to capture long-range correlations between them. The full ViT is a stack of these factored attentions combined with a patchwise Multi-Layer Perceptron to process the information aggregated by the factored attention in each patch vector. This architecture makes the ViT translational invariant on patch level. To impose full translational invariance, we just have to apply it at the level of the patch, by averaging the model over one-site translations.

This same architecture is also used for the SNR analysis on the 4x4 lattice, the only difference being that we employed a patch size of 1, which greatly improved the results for this lattice size.

# Appendix C: Experimental details

## 1. $J_1-J_2$ results

### a. ViT architecture

The Vision Transformer, as introduced in [52], is a deep neural network architecture originally applied to Computer Vision problems. First, one patches the input by grouping together nearby sites on the lattice. It is particularly interesting in the case of the square $J_1-J_2$ model as using patches of lateral size 2 corresponds to the unit cell of the model. Then, a linear embedding layer transforms these patches into $h$ vectors of dimension $d_{\text{model}}/h$, where $h$ denotes the number of heads and $d_{\text{model}}$ the total number of features. We will denote this matrix $X_i^\mu$, where the
### b. Hyperparameters

For training the ViT on the J₁-J₂ model, we used a fixed budget of 6000 iterations using Stochastic Recon- figuration, in which all simulations converged. To regu- larize the Quantum Geometric Tensor (QGT), we add a diagonal shift, initialized at 10⁻² and decreased to 10⁻⁴ during the first 1000 steps using a cosine decay schedule. The learning rate is annealed in a similar manner, from O(10⁻²) to O(10⁻³). To remove bias, every variational energy result is computed as an expectation values under |ψ|².

## 2. Ab-initio quantum chemistry

### a. Neural network backflow

We use an NQS with the same architecture as the one used in Ref. [6] and Ref. [28]. For Li₂O, we use 2 hidden layers, and 256 hidden units. For N₂, the number of hidden units is reduced to 16. To improve the results, we also enforce spin flip symmetry.

### b. Optimization

The Hamiltonians of N₂ and Li₂O are expressed in the STO-3G basis and canonical HF molecular orbitals. We observed that with MCMC, using natural orbitals from CCSD makes the optimization slightly more difficult.

We employ SR or MinSR as an optimization technique depending on the number of variational parameters. We also found that SPRING with a momentum of 0.9 im- proves the training stability, as shown in figure fig. 8.

For lithium oxide, we train for 4000 iterations the learning is decayed linearly from 0.1 to 0.001 in the first 1000 iterations while the diagonal shift is also decreased linearly from 10⁻³ to 10⁻⁴ in 500 steps. These parame- ters are very conservative for many simulations, but could give the overall best results, especially for vanilla VMC using |ψ(x)|².

## 3. Equation summary (implementation aid)

In this section we give all the necessary equations for the reader to implement the method presented in this paper. For details about the notation and derivation, please refer to the main text and later appendices.

### a. VMC

We aim at minimizing:

$$ \mathcal{L}(\theta) = \mathbb{E}_p[\ell(x)] \quad (C1) $$

for some local objective ℓ(x). To estimate the gradient **F** = ∇θ**L**(θ), we resort to:

$$ F_i = \mathbb{E}_p \left[ \underbrace{2 \operatorname{Re}\{\Delta \partial_{\theta_i} \log \psi_{\theta}(x)^* \Delta \ell(x)\}}_{\equiv f_i(x)} \right], $$

where ΔA(x) ≡ A(x) − Eₚ[A]. The self normalized im- portance sampling estimator of the gradient is:

$$ F_i = \mathbb{E}_q \left[ \underbrace{W(x) f_i(x)}_{\equiv f_i^q(x)} \right], $$

where we define W(x) = w(x)/Eₚ[w(x)] and w(x) = p(x)/q(x) = |ψ(x)|²/q(x). In a similar manner, we es- timate the Quantum Geometric Tensor (QGT) as:

$$ S_{ij} = \mathbb{E}_q \left[ \operatorname{Re} \left\{ W(x) \Delta \partial_{\theta_i} \log \psi_{\theta}(x)^* \Delta \partial_{\theta_j} \log \psi_{\theta}(x) \right\} \right], $$

### b. Adaptive tuning

The variance of the gradient estimator can be esti- mated as:

$$ \mathbb{V}_q[f_i^q(x)] = \mathbb{E}_q \left[ W^2(x) |f_i(x) - F_i|^2 \right]. $$

We aim at maximizing the following objective function:

$$ \mathcal{L}_{\text{IS}}(q) = \frac{1}{N_p} \sum_{i=1}^{N_p} \sqrt{\frac{|\mathbb{E}_q[f_i^q(x)]|^2}{\mathbb{V}_q[f_i^q(x)]}} = \frac{1}{N_p} \sum_{i=1}^{N_p} \text{SNR}(f_i^q), $$

over the family of overdispersed sampling distributions:

$$ \left\{ q_{\alpha}(x) = \frac{|\psi_{\theta}(x)|^{\alpha}}{Z_{\alpha}}, \quad \alpha \ge 0 \right\}. $$

We use:

$$ \alpha' = \alpha + \eta_{\text{IS}} \partial_{\alpha} \mathcal{L}_{\text{IS}}(q_{\alpha}). $$

The gradient of the variance with respect to the varia- tional parameters can be expressed as:

$$ \partial_{\alpha} \mathbb{V}_{q_{\alpha}} [f_i^{q_{\alpha}}(x)] = - \mathbb{E}_{q_{\alpha}} [ \Delta_{q_{\alpha}} \partial_{\alpha_i} \log f_{\alpha}(x) W^2(x) |f_i(x) - F_i|^2 ], $$

where

$$ \Delta_{q_{\alpha}} \partial_{\alpha} \log f_{\alpha}(x) = \partial_{\alpha} \log f_{\alpha}(x) - \mathbb{E}_{q_{\alpha}} [\partial_{\alpha} \log f_{\alpha}(x)]. $$

This leads to:

$$ \partial_{\alpha} \text{SNR}(f_i^{q_{\alpha}}(x)) = \frac{1}{2} \partial_{\alpha} \mathbb{V}_{q_{\alpha}} [f_i^{q_{\alpha}}(x)] \frac{|\mathbb{E}_{q_{\alpha}}[f_i^{q_{\alpha}}(x)]|}{\mathbb{V}_{q_{\alpha}}[f_i^{q_{\alpha}}(x)]^{3/2}}, $$

from which we build an estimator of ∂_α**L**IS(q_α).