Understanding population annealing Monte Carlo simulations

Population annealing is a recent addition to the arsenal of the practitioner in computer simulations in statistical physics and it proves to deal well with systems with complex free-energy landscapes. Above all else, it promises to deliver unrivaled parallel scaling qualities, being suitable for parallel machines of the biggest caliber. Here we study population annealing using as the main example the two-dimensional Ising model, which allows for particularly clean comparisons due to the available exact results and the wealth of published simulational studies employing other approaches. We analyze in depth the accuracy and precision of the method, highlighting its relation to older techniques such as simulated annealing and thermodynamic integration. We introduce intrinsic approaches for the analysis of statistical and systematic errors and provide a detailed picture of the dependence of such errors on the simulation parameters. The results are benchmarked against canonical and parallel tempering simulations.

Figure 1

(a) Specific heat as estimated from population annealing with Metropolis update and fixed temperature steps $\mathrm{\Delta }\beta =1/300$ for $L=64$. The pure Metropolis simulation shown for comparison uses $R=1$ with a total of $50000$ sweeps and measurements. The exact reference data are calculated according to Refs. [37, 38]. (b) Same as panel (a), but for the heat-bath update. (c) Results of PA simulations with the same parameters as in panel (a), but with the resampling step turned off. To make the comparison as fair as possible, in contrast to the rest of the paper all runs shown here use $\mu =\theta$ measurements per temperature step. The insets in panels (a) and (b) show detail around the peak close to the critical point of the system.

Figure 2

Inverse of the integrated autocorrelation time ${\tau }_{\mathrm{int}}\left(E\right)$ of the energy estimated from simulations of a $L=128$ system with a single replica and different spin updates.

Figure 3

Effective number of independent replicas as estimated from the family statistics for PA simulations with $L=64$ (Metropolis update, $\mathrm{\Delta }\beta =1/300$). (a) Number $f$ of surviving families. (b) $R/{\rho }_t$ according to Eq. (16). (c) $R/{\rho }_s$ according to Eq. (18).

Figure 4

Distance dependence of the overlap correlation function $C_q(i,j)$ according to Eq. (27) for PA runs of a $L=16$ system with $R=50000$ and $\theta =1$, 5, and 10, respectively. (a) At high temperatures, the population is nearly uncorrelated. (b) Close to criticality, significant correlations develop that asymptotically follow the form $C_q(i,j)\sim exp(-|i-j|/{\tau }_R)$ with ${\tau }_R=32.3$ ($\theta =1$), ${\tau }_R=3.0$ ($\theta =5$), and ${\tau }_R=1.8$ ($\theta =10$), respectively. (c) For the overlap, correlations persist in the ordered phase.

Figure 5

Effective number of statistically independent samples (for the energy measurement) in PA and standard canonical simulations for an $L=64$ system and different population sizes (a) for the Metropolis and (b) for the heat-bath updates [please see panel (a) for the legend of data sets]. Panel (c) shows the effective population sizes $R_{\mathrm{eff}}\left(E\right)$ for the energy and $R_{\mathrm{eff}}\left(M\right)$ for the magnetization for $R=25000$, $\theta =2$ in comparison to the quantity $R/{\rho }_t$ based on family statistics alone, illustrating that the latter is a lower bound for the former, but it is far from being tight. All simulations shown use inverse temperature steps of size $\mathrm{\Delta }\beta =1/300$.

Figure 6

Estimated standard deviations of the mean (i.e., error bars) of the internal energy, $\sigma \left(\overline{E}\right)$, for the $L=32$ Ising model as estimated from PA simulations with the blocking technique using Eq. (26) compared to the unbiased estimate resulting from repeating the full simulation 200 times with different sequences of random numbers ($R=50000$, $\theta =10$, $\mathrm{\Delta }\beta =0.01$). All simulations use the Metropolis update.

Figure 7

Ratio of error bars estimated from the jackknife analysis of a single PA run and from the fluctuation between 200 independent runs for PA simulations of the 2D Ising model with $L=32$ ($\mathrm{\Delta }\beta =0.01$). Deviations from the target value 1.0 occur if the number $R_{\mathrm{eff}}$ of independent measurements becomes too small and hence the blocks in the analysis are no longer independent. As illustrated in panel (c), in this case the blocking or jackknifing approach underestimates the errors.

Figure 8

Effective relative number of independent members of the population, $R_{\mathrm{eff}}/R$, for PA simulations of the $L=64$ system with $\mathrm{\Delta }\beta =0.005$ and varying population sizes in dependence of the $\theta$ values. (a) Runs at temperature $T=2.325$, close to the critical point. (b) Runs at $T=2.5$ in the disordered phase. The lines show fits of the functional form (30) to the data with (a) ${\tau }_{\mathrm{eff}}=51.2$ and (b) ${\tau }_{\mathrm{eff}}=2.2$, respectively.

Figure 9

Relative deviation of the effective population size $R_{\mathrm{eff}}$ from the actual population size $R$ as a function of $\mathrm{\Delta }\beta$ for an $L=64$ system close to criticality for a range of different values of $\theta$ ($\beta =0.43$, $R=10000$). The data are averaged over 200 independent runs. The lines show one-parameter fits of the form $1-R_{\mathrm{eff}}/R=A\mathrm{\Delta }\beta$.

Figure 10

Estimated standard deviation of the mean times $\sqrt{R}$ for the energy of the $L=64$ Ising model for PA simulations with a range of different population sizes $R$ for (a) $\theta =10$ and (b) $\theta =1$, respectively ($\mathrm{\Delta }\beta =0.01$).

Figure 11

Variance of the free-energy estimator $\beta \widehat{F}$ for different choices of $\mathrm{\Delta }\beta$. The purple crosses show the variance of the quantity (33) as estimated from 200 independent PA runs. The green squares and blue circles show the estimates derived via Eqs. (39) and (40), respectively, from the same simulations ($L=64$, $R=10000$, $\theta =50$). For clarity of presentation, only every 10th temperature point is plotted with a symbol.

Figure 12

Relative deviation $\epsilon \left(E\right)=(E-E_{\mathrm{exact}})/E_{\mathrm{exact}}$ of the internal energy at the critical inverse temperature ${\beta }_c=\frac{1}{2}ln(1+\sqrt{2})$ for an $L=32$ system sampled in PA runs without resampling on a population of size $R=10000$. The results are averaged over $m=200$ independent runs to reduce statistical errors. (a) Relative deviation as a function of $\theta$ together with fits of the functional form $\epsilon \left(E\right)=a{e}^{-\theta /{\tau }_{\mathrm{rel}}}(1+b/\theta )$ motivated by (56) to the data. (b) Relative deviation as a function of $\mathrm{\Delta }\beta$ with fits of the functional form $\epsilon \left(E\right)=a\mathrm{\Delta }\beta \left(1-e^{-b/\mathrm{\Delta }\beta }\right)$.

Figure 13

(a) Relative deviation $\epsilon \left(E\right)=(E-E_{\mathrm{exact}})/E_{\mathrm{exact}}$ of the internal energy at the critical inverse temperature ${\beta }_c$ for an $L=32$ system sampled in PA runs with resampling on a population of size $R=10000$ as a function of inverse temperature step $\mathrm{\Delta }\beta$. The interpolating lines are merely guides to the eye. The red dashed line shows the exact histogram overlap $\alpha$ at criticality as a function of step size (right scale). (b) Bias $\left|\mathrm{\Delta }E\right|=|E-E_{\mathrm{exact}}|$ for the run with $\theta =10$ as compared to the covariance $cov(E,\beta \widehat{F})$ of Eq. (57) together with a linear fit for small values of $\mathrm{\Delta }\beta$.

Figure 14

Variance ${\sigma }^2\left(\beta \widehat{F}\right)$ and covariance $cov(\beta \widehat{F},E)$ averaged over the temperature range $0\le \beta \le 1$ for $L=32$ as a function of $\mathrm{\Delta }\beta$. Both scale proportional to $\mathrm{\Delta }\beta$ in the limit of small steps, as the linear fits illustrate.

Figure 15

(a) Relative deviation $\left|\epsilon \left(E\right)\right|=|(E-E_{\mathrm{exact}})/E_{\mathrm{exact}}|$ of the internal energy from PA simulations for $L=64$, $\mathrm{\Delta }\beta =0.01$, and $\theta =1$ (heat-bath update) and population sizes ranging from $R={10}^3$ to $R={10}^6$. (b) The same bias for a population $R=1000$ for different choices of $1\le \theta \le 1000$. (c) Energy deviation relative to the statistical error in a logarithmic representation, illustrating that for $\theta \gtrsim 50$, the bias drops below the statistical fluctuations [symbols of data sets as in panel (b)]. Corresponding data sets in all panels amount to the same total computational effort.

Figure 16

Schematic of energy distributions in $\theta =0$ population annealing. (a) $P_{{\beta }_0=0}$ corresponds to the initial population, and $P_{\beta }$ to the distribution at a lower temperature. (b) In a finite population of size $R$ at ${\beta }_0=0$, there are typically no replicas with energies lower than $-E^{*}$ (and similarly none with energies higher than $E^{*}$), where $P_{{\beta }_0}\left(E^{*}\right)=1/R$.

Figure 17

(a) Bias of the (total, not per site) internal energy estimate from PA simulations for $L=32$ and $\theta =0$ and different population sizes. The lines show the conservative bias estimate resulting from Eq. (74) [see panel (c) for key of line types]. (b) Detail of panel (a). The vertical lines indicate the maximum inverse temperature, ${\beta }^{*}$, from which strong biases are expected according to Eq. (70). (c) Value of the normalized cutoff energy $e_{<}=(-E^{*}-{\langle E\rangle }_{\beta })/{\sigma }_{\beta }$.

Figure 18

Effective population size $R_{\mathrm{eff}}$ at $\beta =0.44$ as estimated from the blocking analysis for PA runs for $L=64$ and the indicated values of $R$ and $\theta$. For $R<R_0\left(\theta \right)$, $R_{\mathrm{eff}}$ is essentially independent of $R$, while for $R\ge R_0\left(\theta \right)$, $R_{\mathrm{eff}}\propto R$.

Figure 19

Run times of sample PA simulations for the 2D Ising model and system sizes $L=16,...$, 128 for $R=5000$, $\theta =10$, $\mathrm{\Delta }\beta =16/75L$. The inset shows the relative excess time as compared to simulations with the resampling step turned off.

Figure 20

Relative squared error bars of observable estimates from Metropolis temperature sweeps against PA runs with the same number of spin flips. The resulting “speedups” relate to (a) the specific heat and (b) the magnetic susceptibility, respectively. PA updates are for $R=50000$, $\theta ={(L/16)}^2$, $\mathrm{\Delta }\beta =16/75L$ and Metropolis runs are for $R=1$, $\theta =50000\times {(L/16)}^2$, and the same temperature sequence. (c) Error bar of the specific heat estimate for $L=32$ and Metropolis and parallel tempering simulations as compared to standard and adaptive PA runs.

Figure 21

Ratio of effective population sizes before and after the resampling step as compared to the inverse of the one-step mean-square family size according to Eq. (A2) for two sizes of the inverse temperature steps ($L=32$, $\theta =20$, $R=50000$).

Figure 22

Effective population size $R_{\mathrm{eff},i}$ at inverse temperature $\beta =0.41$, 0.43, and 0.48, respectively, after applying a number of ${\theta }_i$ MCMC steps. The effective population size relaxes to the total population size $R$ according to Eq. (A4) as ${\theta }_i$ is increased. The lines show fits of the form (A4) with parameter ${\tau }_{\mathrm{rel},i}$ to the data, the best fit being achieved for ${\tau }_{\mathrm{rel},i}=5.3$ ($\beta =0.41$), 33.9 ($\beta =0.43$), and 0.8 ($\beta =0.48$), respectively ($L=32$, $\theta =10$, $R=50000$, $\mathrm{\Delta }\beta =0.01$).

Figure 23

(a) Bias and statistical errors of the average internal energy $\overline{E}$ for the $L=32$ square-lattice Ising model as estimated from independent runs as well as from the jackknife blocking method ($R=50000$, $\theta =1$, $\mathrm{\Delta }\beta =0.01$). (b) The same data in a logarithmic scale shown in comparison to the effective population size $R_{\mathrm{eff}}$ (solid line, right scale).