Optimization of population annealing Monte Carlo for large-scale spin-glass simulations

Population annealing Monte Carlo is an efficient sequential algorithm for simulating $k$-local Boolean Hamiltonians. Because of its structure, the algorithm is inherently parallel and therefore well suited for large-scale simulations of computationally hard problems. Here we present various ways of optimizing population annealing Monte Carlo using 2-local spin-glass Hamiltonians as a case study. We demonstrate how the algorithm can be optimized from an implementation, algorithmic accelerator, as well as scalable parallelization points of view. This makes population annealing Monte Carlo perfectly suited to study other frustrated problems such as pyrochlore lattices, constraint-satisfaction problems, as well as higher-order Hamiltonians commonly found in, e.g., topological color codes.

Figure 1

Diagram outlining the different optimizations we have implemented for population annealing Monte Carlo. These range from optimizations in the implementation, such as efficient spin selection techniques, to algorithmic accelerators (e.g., the inclusion of cluster updates), as well as parallel implementations. See the main text for details.

Figure 2

Comparison of the entropic population size ${\rho }_s$ for different spin selection methods: random, sequential, and checkerboard updates in three space dimensions. Sequential and checkerboard updates have similar efficiency (b), and both are more efficient than random updates (a).

Figure 3

Panel (a) shows the $\beta$ values as a function of the inverse temperature index $k_{\beta }$ for the different schedules experimented with and panel (b) shows the resulting $\beta$ densities $g\left(\beta \right)$ (the data are cut off at $\beta =3$ for clarity). Note that both TSPL and LBLT schedules have more temperatures at high $T$. Panels (c) and (d) show ${\rho }_s$ as a function of $T_N$ for 2D and 3D simulations, respectively. The vertical shaded line marks the optimum. See the main text for details.

Figure 4

Comparison of the systematic errors for various annealing schedules. The studied observables are energy ($E$), free energy ($F$), and the spin glass Binder cumulant ($g_{\mathrm{SG}}$) for the system size $L=10$. Panels (a) and (b) show the systemic errors for two randomly chosen hard instances, whereas panel (c) illustrates the systematic errors averaged over 100 of the hardest instances. Systematic errors of different observable often have magnitudes largely apart. For this reason, the errors in each observable have been normalized relative to the maximum error across all schedules. For instance, in the top panel the $\mathrm{std}\left(E\right)$ schedule which has the greatest systematic error is normalized to 1 while the rest of the schedules lie below 1. It is seen from the plots that the TSPL schedule is the most efficient. The LBLT schedule, although conveniently simple, competes well with the optimal schedule. Note that we also show ${\rho }_s$ (as a dual $y$ axis) in panel (c). We observe that ${\rho }_s$ greatly correlates with the systematic errors justifying the use of it as an effective optimization criterion.

Figure 5

Energy density distribution of the LBLT annealing schedule for $L=8$ in three space dimensions. Thinner curves show the histograms at all temperatures, whereas the thicker ones are drawn at every 10 temperature steps. There are 200 temperature steps in total. The histograms overlap considerably.

Figure 6

Optimization of the number of annealing steps $N_T$ in two space dimensions [2D, panel (a)] and three space dimensions [3D, panel (b)]. To maximize sampling efficiency, one needs to optimize $1/\left({\rho }_s{N}_T\right)$ with respect to $N_T$. In both panels the points and the solid curves show the disorder average while the dashed envelopes display all 100 studied instances. For smaller system sizes the peak (optimum) is sharp, whereas for systems with more than approximately 1000 spins the peak is broadened, especially in two dimensions. The reason for this broadening can be understood by noticing the increase in the density of chaotic samples as the system grows in size (wiggly lines).

Figure 7

Instance-by-instance comparison for a PAMC simulation with fixed and dynamic population sizes. With a dynamic population size, ${\rho }_s$ and ${\rho }_f$ are well correlated, similarly to the case of uniform population. ${\rho }_s$ is greatly reduced, suggesting that the simulation is much better at the level of averaging over all temperatures. The dynamic population size is also more efficient than the uniform one using the worst-case measure. Here $R_f$ is the final population size.

Figure 8

Population annealing with ICM updates in 2D. Note that replica family is not well defined when ICM updates are included. Therefore, we use ${\rho }_f$ to characterize speed-up. Significant speed-up is observed in 2D.

Figure 9

Population annealing with ICM updates in 3D. Note that replica family is not well defined when ICM updates are included. Therefore, we use ${\rho }_f$ to characterize speed-up. Modest speed-up is observed in 3D.

Figure 10

Scaling of the total wall time as a function of the number of processors $N$ for two system sizes $L=8$ and $L=12$. Launching and initialization time are not included. Note that the efficiency becomes better for larger and harder problems. For $L=12$, the scaling remains $1/N$ up to about 1000 processors. The efficiency then decreases when the time for collecting observables becomes dominant. Note that resampling still takes a relatively small time.

Figure 11

Mean normalized cluster size as a function of $\beta$ for the Wolff algorithm [(a) and (c)] as well as the performance of the algorithm in both 2D and 3D. There is marginal speed-up in 2D (b) and no discernible speed-up in 3D (d).