Microcanonical simulated annealing: Massively parallel Monte Carlo simulations with sporadic random-number generation

TL;DR

MicSA addresses the heavy random-number burden in large-scale Monte Carlo simulations by embedding the spin system in a microcanonical ensemble augmented with daemons/walkers, enabling massively parallel updates with minimal randomness. The method preserves detailed balance through Creutz-like moves and leverages a controlled refresh schedule to steer dynamics toward canonical equilibrium, yielding dynamics that map onto the standard Metropolis evolution via simple time rescaling. Demonstrated on the 3D Edwards–Anderson spin glass, MicSA reproduces equilibrium and out-of-equilibrium results and achieves strong parallel performance on GPUs, with potential applicability to FPGAs and ASICs for extreme-scale simulations. The work suggests broad utility for disordered systems and optimization problems where random-number generation is a bottleneck, offering a practical path to scalable, low-noise Monte Carlo computing.

Abstract

Numerical simulations of models and theories that describe complex experimental systems $\unicode{x2014}$in fields like high-energy and condensed-matter physics$\unicode{x2014}$ are becoming increasingly important. Examples include lattice gauge theories, which can describe, among others, quantum chromodynamics (the Standard Model description of strong interactions between elementary particles), and spin-glass systems. Beyond fundamental research, these computational methods also find practical applications, among many others, in optimization, finance, and complex biological problems. However, Monte Carlo simulations, an important subcategory of these methods, are plagued by a major drawback: they are extremely greedy for (pseudo) random numbers. The total fraction of computer time dedicated to random-number generation increases as the hardware grows more sophisticated, and can get prohibitive for special-purpose computing platforms. We propose here a general-purpose microcanonical simulated annealing (MicSA) formalism that dramatically reduces such a burden. The algorithm is fully adapted to a massively parallel computation, as we show in the particularly demanding benchmark of the three-dimensional Ising spin glass. We carry out very stringent numerical tests of the new algorithm by comparing our results, obtained on GPUs, with high-precision standard (i.e., random-number-greedy) simulations performed on the Janus II custom-built supercomputer. In those cases where thermal equilibrium is reachable (i.e., in the paramagnetic phase), both simulations reach compatible values. More significantly, barring short-time corrections, a simple time rescaling suffices to map the MicSA off-equilibrium dynamics onto the results obtained with standard simulations.

Figures (6)

• Figure 1: After a time rescaling, our microcanonical simulated annealing reproduces canonical dynamics. We study the energy relaxation towards its equilibrium value by plotting the excess energy per spin over its equilibrium value, $e(t)- e_\infty$, versus simulation time $t$, as computed for a spin-glass sample on a $160\times160\times160$ lattice with binary couplings at temperature $T=0.7\approx 0.64 T_\mathrm{c}$ ($T_\mathrm{c}$ is the critical temperature for the phase transition separating the paramagnetic phase at high temperatures from the spin-glass phase at low temperature). The sample is investigated using two dynamics: a standard, random-number greedy Metropolis algorithm (the label can in the plot stands for canonical; these data were obtained from the Janus II supercomputer) and the new microcanonical simulated annealing (MicSA) algorithm introduced in this work. Specifically, we are using the version of the algorithm with 6 walkers per spin, see Sec. \ref{['subsect:our-algorithm']}. We plot both data sets, subtracting the same estimate for $e_\infty$, which was obtained from the canonical simulation as explained in \ref{['ap:time-shift']}. After the MicSA time is rescaled by the factor $k^*=6.034$ (see Table \ref{['tab:k*']}), $t^\text{plot}=k^*\times t^\text{MicSA}$, both energy relaxations are indistinguishable (the naif rescaling would have been $k^\text{naif}=6$, because of the 6 walkers per spin).
• Figure 2: Comparison of Metropolis and MicSA dynamics with 6 walkers in the spin-glass phase at $T\!=\!0.7\approx 0.64 T_\mathrm{c}$. a) Time evolution of the difference of the energy at time $t$ and the equilibrium value obtained through a fit to the canonical data. The solid yellow line is for the best fit described in Table \ref{['tab:fit-info']}. b) The coherence length, $\xi_{12}$, as a function of the simulation time for both models. The solid line is again for the best fit described in Table \ref{['tab:fit-info']}. c) Coherence length, $\xi_{12}$, as a function of the ratio between the coherence lengths $\xi_{23}$ and $\xi_{12}$. d) Difference between the results from the two algorithms for the energy, Eq. \ref{['eq:delta1']}, as a function of the simulation time. e) As in d) but for the coherence length. To compute these differences, we have considered the best fit to our canonical data ---solid lines in a) and b)--- and the rescaling coefficient $k^*=6.034$ (see \ref{['ap:time-shift']} for more information on the estimation of $k^*$). Error bars are smaller than data points. For the sake of clarity points are plotted at regular time intervals in logarithmic scale.
• Figure 3: Six walkers at $T=1.1\simeq T_\mathrm{c}$. As in Fig. \ref{['fig:T0.7']}, but with the time-rescaling factor $k^*=5.876$ (see Table \ref{['tab:k*']}).
• Figure 4: Six walkers in the paramagnetic phase at $T\!=\!1.4\approx1.27T_\mathrm{c}$. As in Fig. \ref{['fig:T0.7']}, but with the time-rescaling factor $k^*=5.904$ (see Table \ref{['tab:k*']}).
• Figure 5: Estimation of the best time shift $k^*$ at $T_\mathrm{c}$. In the main panel, the difference between the logarithm of the original measured time in the MicSA data and the interpolating function $f(k^*t)$ as a function of the coherence length $\xi_{12}$. The horizontal solid line at zero indicates the interval of the fit used to estimate $k^*$ (see Table \ref{['tab:k*']} for details). Inset: $\log(t)$ as a function of $\xi_{12}$ for both algorithms. Solid lines indicate the best fit for our data. In the case of the MicSA data, the only fitting parameter is $k$, while the rest of the parameters $\Theta=\lbrace a_0,\,a_1,\, z\rbrace$ come from the fit to the canonical data described in \ref{['ap:time-shift']}.