Mitigating spikes in fermion Monte Carlo methods by reshuffling measurements

We propose a method to mitigate heavy-tailed distributions in fermion quantum Monte Carlo simulations originating from zeros of the fermion determinant. In this case, the second moment of the observables might not be well defined, and we show that by merely changing the synchronization between local updates and computation of observables, one can reduce the prefactor of the heavy-tailed distribution, thus substantially suppressing statistical fluctuations of observables. We also show that the average, or the first moment, is well defined and hence is independent on our measuring scheme. The method is especially suitable for local observables similar to, e.g., double occupancy, where the resulting speedup can reach two orders of magnitude. For observables containing spatial correlators, the speedup is more moderate but still ranges between five and ten. Our results are independent on the nature of the auxiliary field, discrete or continuous, and pave the way to improve measurement strategies for hybrid Monte Carlo simulations.

Figure 1

Sublattice charge fluctuations of Eq. (40) for the global, Eq. (2) (upper plot), T-local, Eq. (3) (middle plot), and ST-local Eq. (4) (lower plot) schemes. Calculations were done on a $6\times 6\times 256$ hexagonal lattice at $U=5.0$ and $\beta =20.0$ using the SU(2) symmetrical discrete Hubbard Stratonovich (HS) field within a standard local update BSS-QMC algorithm [8].

Figure 2

Action of the two-dimensional toy model (13). (a) the action when $K=0.0$ and the fields $x_1$ and $x_2$ are uncorrelated; (b) the action when $K=0.5$. $D=0.0$ in both cases.

Figure 3

Results for the toy model. Figure (a) shows the identical probability distributions for the observable when global or local measurements are used in the case of uncorrelated fields $x_1$ and $x_2$ ($K=0$). Figure (b) is the comparison of the probability distributions for local and global measurements at $K=0.9$. Figure (c) shows the ratio of the widths of the distributions for global and local schemes of measurements [see Eq. (27)]. Larger value of this ratio corresponds to bigger advantage from local measurements. $D=0.0$ in all cases.

Figure 4

Histograms for the distribution of observables after one full sweep. Each plot corresponds to one observable: (a) sublattice spin fluctuations; (b) sublattice charge fluctuations; (c) squared staggered moment; (d) double occupancy. In each case several ways of synchronization between updates and calculation of observable are shown. Calculations were done on the $6\times 6\times 256$ hexagonal lattice at $U=5.0$ and $\beta =20.0$ using the SU(2) symmetrical discrete HS field with standard local Blankenbecler, Scalapino, Sugar (BSS)-QMC updates.

Figure 5

Evolution of the distributions for the sublattice charge fluctuations with increased number of sweeps. The plot shows the share of the configurations, where the observable is less than $-25.0$ (large negative spike) depending on the whole number of configurations. The share remains stable, thus signaling that the distributions shown in the Figs. 4 and 6 are not dependent on the size of statistics. Calculations were done on the $6\times 6\times 256$ hexagonal lattice at $U=5.0$ and $\beta =20.0$.

Figure 6

Comparison of the distribution of observables for different lattices and models: left (right) column shows sublattice charge (spin) fluctuations. First (upper) row corresponds to a $12\times 12\times 256$ hexagonal lattice, the second row shows the results for the $6\times 6\times 512$ hexagonal lattice with a twice smaller step in Euclidean time and the third row corresponds to the square lattice Hubbard model on the $8\times 8\times 256$ lattice. The Hubbard interaction and temperature are the same in all cases: $U=5.0$ and $\beta =20.0$. For comparison with the $6\times 6\times 256$ lattice, one can refer to the upper row in Fig. 4.

Figure 7

Fit of the tail of the histograms to the form $\frac{C}{{(x+A)}^{\xi }}$. Left (right) column shows the sublattice charge (spin) fluctuations. First (upper) row corresponds to the global measurement of an observable, the second row shows the T-local computations and the last (lower) row corresponds to the ST-local scheme of measurement. All data were produced on the $6\times 6\times 256$ lattice at $U=5.0$ and $\beta =20.0$.

Figure 8

(a) Histogram for the distribution of the sublattice charge fluctuations in the calculation with the staggered mass, that eliminates the zeros of the determinant. (c) An attempt to fit the tail of the histogram by a power law (at fixed shift parameter $A$). The inset shows the steady increase of the power in the fit upon increasing $A$. (b) The dependence of ${\chi }^2$ on the shift $A$ for two different cases: with and without the staggered mass term. The computations were carried out on a $6\times 6\times 256$ lattice at $\beta =20, U=5$, and mass $m=0.1$.

Figure 9

Comparison of error bars (width of the 95% confidence interval) for different algorithms depending on the number of sweeps per simulation. We have adjusted the vertical axis to accommodate all three algorithms: global, T-local, and ST-local measurements. The top (bottom) plot shows the results for the sublattice charge (spin) fluctuations. The data set was produced on a $6\times 6\times 256$ hexagonal lattice at $U=5.0$ and $\beta =20.0$. Here, the staggered mass is equal to zero. Insets in both plots show the same data, but the vertical scale adjusted to included only the T-local and ST-local algorithms.

Figure 10

Conservative estimate of the speedup: we take the maximal value within the error bars for the power-law coefficient $\mu$ in case of the T-local measurements and the minimal value of $\mu$ for the ST-local measurements. The speedup is estimated at fixed relative error of the corresponding observable. Data for the power-law fits are taken from Figs. 9 and 14, such that the speedup pertains to simulations on the $6\times 6\times 256$ hexagonal lattice at $U=5.0$ and $\beta =20.0$. Here, the staggered mass is equal to zero.

Figure 11

Comparison of the histograms for the sublattice charge fluctuations obtained in QMC calculations with discrete and continuous auxiliary fields. (a) SU(2) symmetrical discrete fields updated using standard local BSS QMC updates. (b) SU(2) symmetrical continuous fields, updated with the hybrid Monte Carlo algorithm constrained in the dominant thimble [21]. In both cases we employ global measurements: observables are computed on the whole lattice after the update of all fields. Computations were done on a $6\times 6\times 256$ at $\beta =20$ and $U=5$.

Figure 12

Histograms for the distribution of averages of the sublattice charge fluctuations over 25 (a) and 100 (b) sweeps. Vertical lines show the borders of 95% credible interval obtained from Eqs. (A3) and (A2). Computations were done on the $6\times 6\times 256$ lattice at $\beta =20$ and $U=5$. T-local measurements were used to compute the observable.

Figure 13

Dependence of various error estimates on the number of sweeps per simulation. We consider the 95% confidence interval, the 95% credible interval, and the average variance (A4). In all cases we added a power-law fit for the confidence interval data set. Left column corresponds to the ST-local measurements and right column corresponds to the T-local measurements. Different rows show the results for different observables: the first (upper) row demonstrates the error bar for the sublattice charge fluctuations, the second row for the sublattice spin fluctuations, and the last (bottom) row for the double occupancy. All data was produced on the $6\times 6\times 256$ hexagonal lattice at $U=5.0$ and $\beta =20.0$. The staggered mass was set to zero.

Figure 14

Comparison of error bars (width of the 95% confidence interval) for different observables depending on the number of sweeps $N$ per simulation and on the employed algorithm: T-local or ST-local measurements. The top left plot shows the sublattice charge fluctuations, the top right corresponds to the sublattice spin fluctuations, bottom left is the double occupancy and bottom right the squared staggered moment. In all cases we perform power-law fits. The values of the powers $\mu$ for all cases are inserted as labels in the plots. The data set was produced on the $6\times 6\times 256$ hexagonal lattice at $U=5.0$ and $\beta =20.0$. The staggered mass is set to zero.

Figure 15

Comparison of error bars (width of the 95% confidence interval) for the sublattice charge fluctuations (left column) and sublattice spin fluctuations (right column) depending on the number of sweeps per simulation and for different models and lattices. First (upper) row corresponds to a $12\times 12\times 256$ hexagonal lattice, the second row shows the results for a $6\times 6\times 512$ hexagonal lattice and the third row corresponds to the square lattice Hubbard model on an $8\times 8\times 256$ lattice. Here we consider $U=5.0$ and $\beta =20.0$. For comparison with the $6\times 6\times 256$ lattice, one can refer to the upper row in Fig. 14.

Figure 16

Error bars (width of the 95% confidence interval) for the sublattice charge fluctuations as a function of the number of sweeps per simulation. Calculations were done on the $6\times 6\times 256$ hexagonal lattice with the mass term $m=0.1$. The mass term removes the zeros from the fermion determinant. Power-law fits with their respective coefficients $\mu$ are shown in the plot. Here we consider $U=5.0$ and $\beta =20.0$.

Figure 17

Comparison of charge-charge correlators computed with the same statistics (1000 sweeps through the lattice) using ST-local measurements (a) and T-local measurements (b). The distance between sites is in units of the separation between nearest neighbours. Here we consider a $6\times 6\times 256$ lattice at $\beta =20$ and $U=5$.