Variational Monte Carlo

Variational Monte Carlo (VMC) evaluates a trial wavefunction’s energy by sampling its squared modulus. It replaces a high-dimensional integral by local wavefunction ratios or derivatives. The variational upper bound applies to the exact expectation for the trial state; a finite noisy estimate can fall below the ground energy by chance.

Required background. Monte Carlo estimation and MCMC supply sampling controls. The variational principle supplies the energy bound; Rayleigh–Ritz explains its finite-basis implementation.

The local energy estimator

For a normalizable trial state $\psi_\theta(x)$ in the Hamiltonian’s domain, define

$$p_\theta(x)=\frac{|\psi_\theta(x)|^2}{Z_\theta}, \qquad E_L(x)=\frac{(H\psi_\theta)(x)}{\psi_\theta(x)}.$$

Then $E(\theta)=\mathbb E_{p_\theta}E_L$. Nodes have zero sampling probability but can produce large nearby local-energy fluctuations; smoothness, cusp conditions, and finite estimator variance must be checked.

For a discrete basis,

$$E_L(x)=\sum_y H_{xy}\frac{\psi_\theta(y)}{\psi_\theta(x)}.$$

A sparse local Hamiltonian can make this sum inexpensive even when the Hilbert space is large. Evaluate logarithmic amplitudes and phases when direct wavefunction values underflow. Test the Hamiltonian row convention: complex conjugating the wrong index can give a plausible but incorrect real energy.

The oscillator Gaussian calculation

Use $H=-\frac12d^2/dx^2+\frac12x^2$ and $\psi_a(x)\propto e^{-ax^2/2}$ with $a>0$. Differentiation gives

$$E_L(x)=\frac a2+\frac{1-a^2}{2}x^2, \qquad E(a)=\frac14\left(a+\frac1a\right).$$

Since $\langle x^2\rangle=1/(2a)$ and $\operatorname{Var}(x^2)=1/(2a^2)$,

$$\operatorname{Var}(E_L)=\frac{(1-a^2)^2}{8a^2}.$$

At $a=1$, every local-energy sample equals $1/2$: the exact eigenstate has zero local-energy variance. At $a=2$, the mean is $5/8$ and the variance is $9/32$. These exact values test the sampler, derivatives, and uncertainty estimate independently.

Zero full energy variance identifies an eigenstate, though not necessarily the lowest one. The sampled local-energy variance equals that full variance only if the sampling support captures the residual, as the discrete example below shows.

Energy gradients and optimization bias

For real parameters, let $O_i=\partial_{\theta_i}\log\psi_\theta$ on the sampling support. When differentiating the expectation is justified and no contribution from outside that support is omitted, a useful gradient estimator is

$$\partial_i E =2\operatorname{Re} \left[ \langle O_i^*E_L\rangle -\langle O_i^*\rangle\langle E_L\rangle \right].$$

It follows by differentiating the normalized expectation and using Hermiticity. Subtracting the mean energy acts as a covariance baseline. A stochastic-reconfiguration update additionally uses the covariance matrix of logarithmic derivatives as a metric, with regularization for nearly redundant directions.

The optimization sample is not a clean validation set. Repeatedly selecting the lowest noisy estimate introduces selection bias. After fixing the parameters, run independent production chains and report their energy and variance. When reusing samples across nearby parameters, monitor reweighting overlap.

Exact zeros in a discrete trial state

Sampling $|\psi|^2$ cannot visit configurations where $\psi$ vanishes exactly. In a finite orthonormal basis, let $\psi$ be normalized and $S=\{x:\psi(x)\ne0\}$. Splitting the squared residual into its components gives

$$\|(H-E)\psi\|^2 =\mathbb E_p|E_L-E|^2 +\sum_{x\notin S}|(H\psi)(x)|^2.$$

The first term is the sampled variance, with the modulus appropriate for complex local energies. The second is absent from every sample. This is the sampling version of the distinction between a projected and full residual.

For a concrete check, take real $g\ne0$, real $\Delta$, and a real dimensionless parameter $\theta$:

$$H=\begin{pmatrix}0&g\\g&\Delta\end{pmatrix}, \qquad \psi_\theta=\frac{(1,\theta)^T}{\sqrt{1+\theta^2}}, \qquad E(\theta)=\frac{2g\theta+\Delta\theta^2}{1+\theta^2}.$$

At $\theta=0$, only configuration zero is sampled, and its local energy is zero. Every local-energy sample therefore has zero variance, although $H\psi_0=(0,g)^T$ and the full energy variance is $g^2$. The energy mean itself remains correct: $E(0)=0$.

The same zero also defeats the support-only gradient estimator. Its covariance vanishes, but direct differentiation gives $E'(0)=2g$. The missing contribution comes from $\partial_\theta\psi_\theta(1)|_{\theta=0}=1$: changing the parameter opens an unsampled configuration. This support issue in variational forces is analyzed by Sinibaldi et al. (2023). Full support, or parameter derivatives that vanish outside a fixed support, avoids this particular missing-gradient term; other variance and convergence conditions still apply.

A symmetry-restricted support is legitimate when it is invariant under $H$ and the intended variations stay within it. Otherwise, test exact sums on a small system and use a sampling distribution and estimator that include the missing contributions. More samples from the same restricted distribution cannot recover them. Isolated measure-zero nodes in continuous configuration space are different: one must examine integrability and the operator domain, rather than assign them the discrete sum above.

A convergence record for a trial family

Separate sampling uncertainty, optimizer convergence, and ansatz limitation. Increase production length at fixed parameters; then vary optimization tolerances and starts; then enlarge the wavefunction family. A tiny standard error only addresses the first limit.

Complex states require phase-consistent ratios. Singular potentials require physically appropriate local behavior. Sampling only one disconnected nodal pocket can bias observables unless symmetry and weights justify the restriction.

Exercises

Minimize the Gaussian energy analytically and evaluate its curvature at the optimum.

Solution

$E'(a)=(1-a^{-2})/4$ vanishes at $a=1$, and $E''(1)=1/2>0$. For $a=1+\delta$, $E=1/2+\delta^2/4+O(\delta^3)$. Energy errors are quadratic near the optimum, so a very accurate energy can coexist with visibly imperfect wavefunction parameters.

For the two-state example with $\theta\ne0$, find both local energies. As $\theta\to0$, explain why zero sampled variance at $\theta=0$ does not follow by taking the limit of the full-support variance. How many independent draws are needed to have a reasonable chance of visiting configuration one?

Solution

The local energies are $E_L(0)=g\theta$ and $E_L(1)=g/\theta+\Delta$, with $p(1)=\theta^2/(1+\theta^2)$. For every nonzero $\theta$, both states are sampled and the local-energy variance equals the full energy variance. Its limit is $g^2$: the rare state’s probability tends to zero while its squared local energy diverges. At exactly zero that state is excluded, so the sampled variance instead vanishes. The probability of no visits in $N$ independent draws is $(1-p(1))^N$; obtaining an order-one chance of a visit requires $N$ of order $\theta^{-2}$. Reliable error estimation generally needs many such visits, and Markov correlations can increase the cost.

References

• Foulkes, W. M. C., Lubos Mitas, R. J. Needs, and G. Rajagopal, “Quantum Monte Carlo simulations of solids,” Reviews of Modern Physics 73, 33–83 (2001), doi:10.1103/RevModPhys.73.33 — variational and diffusion Monte Carlo with error analysis.
• Sinibaldi, Alessandro, Clemens Giuliani, Giuseppe Carleo, and Filippo Vicentini, “Unbiasing time-dependent Variational Monte Carlo by projected quantum evolution,” Quantum 7, 1131 (2023), doi:10.22331/q-2023-10-10-1131, arXiv:2305.14294 — support-dependent bias in stochastic variational forces, especially Section 2, equations (8)–(9).