Infinite variance in fermion quantum Monte Carlo calculations

For important classes of many-fermion problems, quantum Monte Carlo (QMC) methods allow exact calculations of ground-state and finite-temperature properties without the sign problem. The list spans condensed matter, nuclear physics, and high-energy physics, including the half-filled repulsive Hubbard model, the spin-balanced atomic Fermi gas, and lattice quantum chromodynamics calculations at zero density with Wilson Fermions, and is growing rapidly as a number of problems have been discovered recently to be free of the sign problem. In these situations, QMC calculations are relied on to provide definitive answers. Their results are instrumental to our ability to understand and compute properties in fundamental models important to multiple subareas in quantum physics. It is shown, however, that the most commonly employed algorithms in such situations have an infinite variance problem. A diverging variance causes the estimated Monte Carlo statistical error bar to be incorrect, which can render the results of the calculation unreliable or meaningless. We discuss how to identify the infinite variance problem. An approach is then proposed to solve the problem. The solution does not require major modifications to standard algorithms, adding a “bridge link” to the imaginary-time path integral. The general idea is applicable to a variety of situations where the infinite variance problem may be present. Illustrative results are presented for the ground state of the Hubbard model at half-filling.

Figure 1

Distribution of the computed ground-state energies from 100 independent determinantal QMC calculations. In the left panel, each shade band indicates one standard deviation, as computed from the 100 data points, which are plotted vs the (arbitrary) run index. In the right panel, the histogram of the 100 data points are shown with a bin size of the computed standard deviation. For comparison, the theoretical Gaussian distribution from the CLT is also shown. (A shift is applied on the horizontal axis so the exact result is shown at zero, the vertical red line.) The computed mean from the 100 data are shown by the vertical blue line, with its thickness indicating the statistical error bar. Note the logarithmic scale on the vertical axis.

Figure 2

Normalized histograms of MC measurements for $\alpha =0.2$ (top row) and $\alpha =0.0$ (bottom row). The calculations at each $\alpha$ are done with $5\times {10}^6$ MC samples. Each histogram is obtained by grouping a different number (“block”) of samples together to make one entry of $\langle y\left(\alpha \right)\rangle$. In the top row, they converge quickly to Gaussian distributions as “block” is increased and reach agreement with the red (dashed) curves, which indicate the Gaussian as defined according to the CLT. In the bottom row, in contrast, the histograms do not converge to Gaussians. A persistent discrepancy is seen between the CLT prediction and the data. The insets, which display the un-normalized histogram values (semilog scale), highlight the long tail on the right.

Figure 3

Comparison of finite and infinite variance calculations. The variance, statistical error on the expectation value, and the statistical error on the computed variance are shown for $\alpha =0.2$ (top row) and $\alpha =0$ (bottom row). The statistical errors are estimated from 30 independent MC runs. For $\alpha =0.2$, the computed variance remains consistent with analytic results as the sample size $M$ is varied, and the computed statistical errors on the variance and on the expectation values decreases as $1/\sqrt{M}$. For $\alpha =0$ the MC variance shows increasing fluctuations as $M$ is increased and does not converge. (As a guide, the dashed red line plots Eq. (14) with $\alpha$ replaced by $1/M$.) The statistical errors do not decrease following $1/\sqrt{M}$, as is especially evident in the error on the variance.

Figure 4

Detection and further analysis of the infinite variance problem in the Hubbard model calculation in Sec. 2. The top panel shows a reblocking analysis similar to that in Fig. 2. The histograms of the computed ground-state energy do not converge to Gaussians and do not follow the CLT. In the bottom panel, the computed variance and the statistical error on the variance are shown vs sample size, similarly to Fig. 3. The variance does not converge to a finite value. Its error bar grows with sample size in contrast with the expected $1/\sqrt{M}$ decay. (The magenta line is a linear fit.)

Figure 5

The new method applied to the problem in Fig. 4. The top row shows histograms of the expectation values in the (a) denominator and (b) numerator of the new ground-state energy estimator, compared with the CLT analysis. The middle row shows the respective variances, together with the computed error bars on the variances, versus sample size. The purple dashed lines, which plot $s/\sqrt{M}$, indicate the expected behavior of the error bars. The bottom panel plots the size of the computed error bars on the variances vs $1/\sqrt{M}$. The red solid lines show a linear fit whose slopes give the values of $s$ above for the denominator and numerator, respectively.

Figure 6

Comparison between the standard (top panel) and new (bottom panel) methods for spin correlations in the $4\times 4$ Hubbard model with periodic boundary conditions and $U=8t$. Results from exact diagonalization (ED) are shown by the red dashed lines. The insets show the deviations from ED. The three separate QMC results in each set, with increasing sample size, are displayed with small horizontal shifts for clarity. In the new method the statistical error bars decrease as expected, while with the standard algorithm a drastic increase is seen in the largest run. Note that the vertical scale in inset (a) is 5 times that in inset (b).

Figure 7

Comparison of the standard and new algorithms: Spin correlations (staggered) in a larger lattice. Results from the two sets of calculations are shown side by side, with a small horizontal shift for clarity. To aid the eye, those from the new method are connected by a red dashed line. The system is an $8\times 8$ Hubbard model with periodic boundary condition and $U=0.5t$. The inset shows a zoom of the segment indicated by the dotted purple rectangular box.

Figure 8

Accurate and reliable predictions of long-range order. The main figure plots spin correlations at $U=8t$, again for an $8\times 8$ Hubbard model. The inset is an enlargement of the part indicated by the dashed box. The two runs from the standard algorithm with different random number sequences but otherwise identical parameters show drastically different results. The new algorithm, using the same parameters, provides results with small and reliable error estimates to allow determination of the magnetic correlations.