Download Optimizing Temperature Sets in Parallel Tempering Monte Carlo Simulations and more Papers Physics in PDF only on Docsity!
J. Stat. Mech. (2006) P
ournal of Statistical Mechanics:
J^ An IOP and SISSA journal
Theory and Experiment
Feedback-optimized parallel tempering
Monte Carlo
Helmut G Katzgraber^1 , Simon Trebst^1 ,^2 ,^3 , David A Huse^4
and Matthias Troyer^1
(^1) Theoretische Physik, ETH Z¨urich, CH-8093 Z¨urich, Switzerland (^2) Computational Laboratory, ETH Zentrum, CH-8092 Z¨urich, Switzerland (^3) Microsoft Research and Kavli Institute for Theoretical Physics, University of California, Santa Barbara, CA 93106, USA (^4) Department of Physics, Princeton University, Princeton, NJ 08544, USA E-mail: [email protected], [email protected], [email protected] and [email protected]
Received 3 February 2006 Accepted 9 March 2006 Published 29 March 2006
Online at stacks.iop.org/JSTAT/2006/P doi:10.1088/1742-5468/2006/03/P Abstract. We introduce an algorithm for systematically improving the efficiency of parallel tempering Monte Carlo simulations by optimizing the simulated temperature set. Our approach is closely related to a recently introduced adaptive algorithm that optimizes the simulated statistical ensemble in generalized broad-histogram Monte Carlo simulations. Conventionally, a temperature set is chosen in such a way that the acceptance rates for replica swaps between adjacent temperatures are independent of the temperature and large enough to ensure frequent swaps. In this paper, we show that by choosing the temperatures with a modified version of the optimized ensemble feedback method we can minimize the round-trip times between the lowest and highest temperatures which effectively increases the efficiency of the parallel tempering algorithm. In particular, the density of temperatures in the optimized temperature set increases at the ‘bottlenecks’ of the simulation, such as phase transitions. In turn, the acceptance rates are now temperature dependent in the optimized temperature ensemble. We illustrate the feedback-optimized parallel tempering algorithm by studying the two-dimensional Ising ferromagnet and the two-dimensional fully frustrated Ising model, and briefly discuss possible feedback schemes for systems that require configurational averages, such as spin glasses. Keywords: classical Monte Carlo simulations, other numerical approaches, analysis of algorithms
ArXiv ePrint: cond-mat/
©^ c2006 IOP Publishing Ltd and SISSA 1742-5468/06/P03018+22$30.
J. Stat. Mech. (2006) P
Contents
- Introduction 2
- Parallel tempering Monte Carlo 3
- Temperature set optimization 4
- Results 7 4.1. Ferromagnetic Ising model........................... 7 4.2. Fully frustrated Ising model.......................... 12 4.3. Spin glasses................................... 16
- Conclusions 20
Acknowledgments 21 References 21
- Introduction
The free energy landscapes of complex systems are characterized by many local minima that are separated by entropic barriers. The simulation of such systems with conventional Monte Carlo methods is slowed down by long relaxation times due to the suppressed tunnelling through these barriers. Extended ensemble simulations address this problem by broadening the range of phase space which is sampled in the respective reaction coordinate. Recently, an adaptive algorithm [1] has been introduced that explores entropic barriers by sampling a broad histogram in a reaction coordinate and iteratively optimizes the simulated statistical ensemble defined in the reaction coordinate to speed up equilibration. The key idea of the approach is to measure the local diffusivity along the reaction coordinate, thereby identifying the bottlenecks in the simulations and then using this information to systematically shift statistical weights towards these bottlenecks in a feedback loop. The optimized histogram converges to a characteristic form exhibiting peaks at the bottlenecks of the simulation, e.g., in the vicinity of the entropic barriers. The simulation of an optimized ensemble results in equilibration times that can be substantially lower compared to those for other extended ensemble simulations that aim at sampling a flat histogram in the respective reaction coordinate [1, 2]. Flat-histogram algorithms include the multicanonical method [3, 4], broad histograms [5] and transition matrix Monte Carlo [6] when combined with entropic sampling, as well as the adaptive algorithm of Wang and Landau [7, 8]. Parallel tempering Monte Carlo [9] has proven to be a strong ‘workhorse’ in fields as diverse as chemistry, physics, biology, engineering and materials science [10]. Like replica Monte Carlo [11], simulated tempering [12] and extended ensemble methods [13], the algorithm aims to overcome entropic barriers in the free energy landscape by simulating a broad range of temperatures. This allows the system to escape metastable states when wandering to higher temperatures and subsequently relaxing at lower temperatures again on timescales that are substantially smaller than for conventional simulations at a fixed
J. Stat. Mech. (2006) P
depends on the simulated statistical ensemble, that is the choice of temperature points {T 1 , T 2 ,... , TM } in the parallel tempering simulation. In this paper, we present an algorithm that systematically optimizes the simulated temperature set to maximize the number of round-trips between the two extremal temperatures T 1 and TM for each replica and thereby substantially improve equilibration of the system at all temperatures. Conventional approaches assume that to achieve this goal, the simulated temperature set should be chosen in such a way that the probability for replica swap moves at neighbouring temperatures should be ‘flat’, that is approximately independent of temperature. If the specific heat of the system is assumed to be constant, then a good approximation for such a temperature set can be attained with a geometric progression [17]. Given a temperature range [T 1 , TM ], the intermediate M −2 temperatures can be computed via
Tk = T 1
k∏− 1
i=
Ri Ri = M^ −^1
TM
T 1
The geometric progression peaks the number of temperatures around the minimum temperature T 1 where a slower relaxation is assumed. This is not optimal when the system has a diverging specific heat (at an intermediate temperature), because in order to ensure enough overlap between the energy distributions of neighbouring temperatures ∆Ti,i+1 ∼ CV Ti/
N, where CV is the specific heat per spin and N the number of spins, the acceptance probabilities are inversely correlated with the functional behaviour of CV via the inverse beta function law [17]. Thus, for example at a phase transition where the specific heat diverges, the acceptance probabilities for a geometric temperature set will show a pronounced dip (cf section 4.1). More complex methods such as the approach by Kofke [14, 15], its improvement by Rathore et al [16], as well as the method suggested by Predescu [17, 18] aim to obtain acceptance probabilities for the parallel tempering moves that are independent of temperature by compensating for the effects of the specific heat. In particular, Kone and Kofke [19] suggest that an acceptance probability of 23% is optimal. In this work we show that this is not necessarily the optimal case.
- Temperature set optimization
Our goal is to vary the temperature set {Ti} of a parallel tempering simulation in such a way that for a given system we speed up equilibration at all temperatures. To accomplish this, we maximize the rate of round-trips that each replica performs between the two extremal temperatures Tmin = T 1 and Tmax = TM following a similar approach to the ensemble optimization technique presented in [1]. For a given temperature set, we can measure the diffusion of a replica through temperature space by adding a label ‘up’ or ‘down’ to the replica that indicates which of the two extremal temperatures, Tmin and Tmax respectively, the replica has visited most recently. The label of a replica changes only when the replica visits the opposite extreme. For instance, the label of a replica with label ‘up’ remains unchanged if the replica returns to the lowest temperature Tmin, but changes to ‘down’ upon its first visit to Tmax. This is illustrated in figure 1. For each temperature point in the temperature set {Ti}, we record two histograms nup(Ti) and ndown(Ti). Before attempting a sequence of swap moves, we increment that histogram at temperature Ti which has the label of the respective replica currently at temperature
J. Stat. Mech. (2006) P
τrt
T min
T max
τ
τup τdown
Figure 1. Sketch of the random walk that a given replica performs in temperature space in the course of the simulation. Ideally, the replica will wander up (τup) and down (τdown) in the simulated temperature range [Tmin, Tmax]. The goal of the feedback method is to maximize the number of round-trips each replica performs in this temperature range, and thereby minimize the average round-trip time τrt = τup + τdown.
Ti. If a replica has not yet visited either of the two extremal temperatures, we increment neither of the histograms. This allows us to evaluate for each temperature point the fraction of replicas which have visited one of the two extremal temperatures most recently (e.g., Tmin) as
f (Ti) =
nup(Ti) nup(Ti) + ndown(Ti)
The labelled replicas define a steady-state current j from Tmin to Tmax that is independent of temperature, e.g., the rate at which ‘up’ walkers are created at Tmin and— in equilibrium—absorbed at Tmax. In the following we assume that T is a continuous variable, independent of the temperature points in the current temperature set. We can then determine the current j to first order in the derivative as
j = D(T )η(T )
df dT
where D(T ) is the local diffusivity at temperature T and the derivative df /dT is estimated by a linear regression based on the measurements in equation (3); η(T ) is a density distribution indicating the probability for a replica to reside at temperature T. We approximate η(T ) with a step function η(T ) = C/∆T , where ∆T = Ti+1 − Ti is the length of the temperature interval around temperature Ti < T < Ti+1 for the current temperature set. The normalization constant C is chosen such that ∫ (^) TM
T 1
η(T ) dT = C
∫ TM
T 1
dT ∆T
Rearranging equation (4) gives a simple measure of the local diffusivity D(T ) of a replica at temperature T
D(T ) ∼
∆T
df /dT
J. Stat. Mech. (2006) P
- For the given temperature set estimate an optimized probability distribution η′(T ) via
η′(T ) = C′
∆T
df dT
- Obtain the optimized temperatures {T (^) i′ } via ∫ (^) T (^) k′
T 1
η′(T ) dT =
k M
- Increase the number of swaps Nsw ← 2 Nsw.
- Stop once the temperature set {Ti} has converged.
The initial number of swaps Nsw should be chosen large enough that a few of round- trips are recorded, thereby ensuring that steady-state data for nup(T ) and ndown(T ) are measured. The derivative df /dT can be determined by a linear regression, where the number of regression points is flexible. Initial batches with the limited statistics of only a few round-trips may require a larger number of regression points than subsequent batches with smaller round-trip times and better statistics.
- Results
4.1. Ferromagnetic Ising model
In order to illustrate the feedback method, we start with a simple test model, the two- dimensional ferromagnetic Ising model (FM). The Hamiltonian for the model is given by
HFM = −J
〈i,j〉
SiSj , (12)
where J = 1 and Si = ±1 represent Ising spins on a square lattice with N = L^2 spins. In our simulations we apply periodic boundary conditions and consider nearest-neighbour interactions only. The simple Ising model with uniform couplings has no frustration or disorder, and exhibits a second-order phase transition at Tc = 2/ ln(1 +
- ≈ 2 .269 from a magnetically ordered to a paramagnetic phase. For simplicity, we define an initial temperature set {Ti} with M = 21 temperature points performing a geometric progression, equation (2), for a temperature interval [T 1 = 0. 1 , TM = 10.0]. The minimum temperature T 1 is chosen low enough that the system can approach the zero-temperature ground state of the model, and the maximum temperature TM is chosen well above the critical region of the phase transition. For a short parallel tempering simulation (Nsw = 1. 6 × 107 , one parallel tempering swap after each lattice sweep) using this initial temperature set, we measure the diffusive current of replicas wandering from the lowest to the highest temperature using an additional label for each replica as described above. In figure 2 we show the fraction of replicas diffusing ‘up’ in temperature space. For the geometric temperature set a sharp drop between two temperature points is observed close to the critical region of the phase transition at Tc ≈ 2 .269. Similar results are also found for a ‘flat’ temperature set where the acceptance rates are approximately constant around ∼40% (with fluctuations of ∼10%) and independent of temperature, although the drop is not as pronounced.
J. Stat. Mech. (2006) P
Figure 2. Fraction f (T ) of replicas diffusing from the lowest to the highest temperature as a function of the temperature index for the ferromagnetic Ising model. For the initial temperature set based on a geometric progression (filled squares), the fraction shows a sharp drop between two temperature points. A similar behaviour is found for a temperature set where the acceptance rates are constant ∼40% independent of temperature (temperature set with ‘flat’ acceptance rates, open squares). In contrast, for the optimized temperature set (triangles) the fraction constantly decreases. The inset shows the fraction as a function of temperature T. The dashed line in the inset represents the critical temperature of the two-dimensional Ising model, Tc ≈ 2 .269.
Calculating the local diffusivity D(T ) for the random walk in temperature space using equation (6), we find a strong suppression around this critical region as illustrated in figure 3. When increasing the size of the simulated system, the dip in the local diffusivity further proliferates, an additional indicator that the slowdown of the random walk in temperature space is dominated by the occurrence of a phase transition. Note the logarithmic scale of the ordinate axis in figure 3.
When applying the feedback, equation (11), to define a new temperature set, this suppression in the local diffusivity leads to a concentration of temperature points near the critical temperature as shown in figure 4. The feedback thereby tries to compensate for the reduced mobility of replicas in the critical region by reallocating resources towards this temperature range. In contrast, the density of temperatures close to the lowest temperature is greatly reduced, thereby suppressing constant swapping of replicas at low temperatures where for the initial temperature set multiple replicas converged to the fully polarized ground-state configuration.
This effect becomes even more evident when measuring the acceptance probabilities for swap moves as illustrated in figure 5. The acceptance probabilities for the initial temperature set based on a geometric progression saturate close to unity for temperatures
J. Stat. Mech. (2006) P
Figure 4. Temperature sets for the ferromagnetic Ising model for different feedback steps. Starting from a geometric progression temperature set (step 0), we apply a feedback loop until the temperature set converges. While the geometric progression places many temperatures at low temperatures, the density of temperatures after the feedback optimization is highest at the bottleneck of the simulation around the critical temperature (marked by a vertical dashed line). Rapid convergence of the optimized temperature set is found after 3–4 feedback steps and a total of Nsw ≈ 1. 6 × 107 swap moves in our parallel tempering simulations.
reallocates resources towards the bottleneck of the simulation, e.g., in the vicinity of a phase transition or close to the ground state of the system.
Measuring the diffusion of replicas in a subsequent simulation for the optimized temperature set, we find that the current of replicas wandering from the lowest to the highest temperature is now characterized by a constantly decreasing fraction f (T ) in agreement with equation (9) as shown in figure 2 (triangles). In addition, we find that for the optimized temperature set the replicas wander evenly up and down in temperature space, that is τup ≈ τdown. In figure 6 we show the ratio between the mean round-trip times τ (^) rt before optimization from a geometric and ‘flat’ temperature set divided by the mean round- trip times τ optrt after optimization in order to illustrate the speedup in replica diffusion attained by the feedback procedure. The data show clearly for all system sizes that the round-trip times after the optimization of the temperature set do not increase as fast as for a geometric progression or ‘flat’ temperature set. Note that for these temperature sets with a fixed number of temperature points the average round-trip times before and after the feedback optimization scale ∼ exp(aLb) for the system sizes studied. It is important to note that our algorithm identifies the bottleneck of the parallel tempering simulation that in this model occurs in the form of critical slowing down at the phase transition
J. Stat. Mech. (2006) P
Figure 5. Acceptance probabilities A(T ) as a function of temperature T for the ferromagnetic Ising model. The inset shows the acceptance rates as a function of temperature in the optimized case for varying system sizes L and a fixed number of temperature points. The vertical dashed line marks the critical temperature.
Figure 6. Average round-trip times τ (^) rt before the optimization divided by the average round-trip times after the feedback optimization (τ optrt ) as a function of system size. The data for the filled squares are for a system starting from a geometric progression and represent the speedup obtained by the feedback method. The open symbols correspond to a temperature set which initially has ‘flat’ acceptance probabilities. The dashed lines are guides to the eye.
J. Stat. Mech. (2006) P
Figure 7. Fraction f (T ) of replicas diffusing from the lowest to the highest temperature for the fully frustrated Ising model. Displayed are data for an initial ‘flat’ temperature set with M = 21 temperature points that yields temperature- independent acceptance probabilities for swap moves (open squares). In addition, data for a geometric progression (M = 21) are also shown (filled squares). After the optimization, the change in the fraction is independent of the temperature index (triangles). The inset shows the fractions as a function of temperature T. Data for Nsw = 2 × 106.
Figure 8. Local diffusivity D(T ) of a random walk in temperature space for the fully frustrated Ising model as a function of temperature T after the feedback optimization for several system sizes L. Notice the logarithmic vertical scale.
J. Stat. Mech. (2006) P
Figure 9. Energy per spin e(T ) = (1/N )[〈H〉]av as a function of temperature T for the fully frustrated Ising model for several system sizes. The data show that already for T 0 .5 the energy is independent of temperature, thus signalling that the system has reached the ground state. The inset zooms into the temperature range around T = 0.
the diffusivity in the vicinity of this bottleneck remains unchanged with respect to the feedback. By applying the feedback method, additional temperature points are shifted towards this bottleneck which is illustrated in figure 10 for the geometric progression (full symbols) as well as the ‘flat’ temperature set (open symbols). For moderately large system sizes, we again find rapid convergence of the generated temperature sets within 3– feedback steps and a total of Nsw ≈ 7. 5 × 106 swap moves. For the optimized temperature set, the diffusive current is again characterized by a fraction of replicas drifting from the lowest to the highest temperature that linearly decreases with the temperature index; see figure 7 (triangles).
The acceptance probabilities A(T ) for replica swap moves are shown in figure 11 for simulations with a geometric progression and the optimized temperature set. While the acceptance probabilities peak at unity close to zero and dramatically drop in the geometric progression temperature set, in the optimized set most temperatures are reshuffled to the low temperature region slightly above zero temperature where the system enters the highly degenerate ground-state manifold. The inset to figure 11 shows the acceptance probabilities as a function of temperature for a fixed number of temperatures, as well as Tmin = 0.1 and Tmax = 10.0 fixed. As in the case for the ferromagnet, the ‘mean’ acceptance probability away from the ground-state bottleneck seems almost independent of temperature and settles at a value that is determined by the number of temperatures used for a given system size L.
In order to test the efficiency of the feedback method on the FFIM, in figure 12 we show the ratio between the mean round-trip times τ (^) rt before optimization divided by the
J. Stat. Mech. (2006) P
Figure 11. Acceptance probabilities A(T ) for replica swap moves as a function of temperature T for the fully frustrated Ising model. While the acceptance probabilities for a geometric progression temperature set show a pronounced dip close to T = 0, the optimized ensemble shows a peak close to zero temperature where the system enters the ground-state manifold. The inset shows, for a fixed number of temperatures, the acceptance rates as a function of temperature for different system sizes L. As for the Ising model, the ‘mean’ value away from the bottlenecks can be tuned by increasing the number of temperatures. This illustrates that in order to obtain higher acceptance rates away from the bottlenecks of the simulation, the number of temperatures has to be increased with increasing L.
define as the product of the average round-trip time and the number of temperatures, the minimum is more pronounced (inset to figure 13) and clearly illustrates that while a larger number of temperatures has little effect on the round-trip times, the total CPU time increases drastically with increasing M. Because the data for the average round-trip times versus M have an optimal value, it is conceivable to develop a feedback optimization method that optimizes both the position of the temperatures and the number of temperatures M. Furthermore, in addition to optimizing the number and locations of the temperature points, we have also explored varying the ratio of the number of sweep moves attempted to the number of replica swap moves attempted, since this is yet another parameter one must set in a parallel tempering simulation. This ratio can be adjusted globally (the same ratio at all temperatures) or locally (the ratio optimized independently at each temperature). This will be discussed in more detail in a subsequent communication.
4.3. Spin glasses
Since the optimization of temperature sets improves the sampling for the two paradigmatic spin models discussed above, it is a natural step to ask how this feedback optimization technique can be applied to improve state-of-the-art parallel tempering simulations of
J. Stat. Mech. (2006) P
Figure 12. Average round-trip times τ (^) rt before the optimization divided by the average round-trip times after the feedback optimization (τ optrt ) as a function of system size for the fully frustrated Ising model. The data for the filled squares are for a system starting from a geometric progression temperature set and represent the speedup obtained by the feedback method. In addition, we show data for a temperature set with ‘flat’ temperature-independent acceptance rates (open squares). The dashed lines are guides to the eye.
Ising spin glass models, such as the three-dimensional (3D) Edwards–Anderson Ising spin glass [21, 25]:
HSG = −
〈i,j〉
Jij SiSj. (15)
Here the spins lie on the vertices of a cubic lattice with periodic boundary conditions. The bonds Jij are chosen according to a Gaussian distribution with zero mean and standard deviation unity. The system undergoes a spin glass transition at Tc = 0.951(9) [26]–[28]. For the spin glass there is the additional complexity that different disorder realizations can lead to strong sample-to-sample variations. Thus one could also surmise that strong sample-to-sample variations exist for the time it takes to equilibrate individual samples. This becomes evident when measuring the round-trip times for replicas in a standard parallel tempering simulation with a fixed temperature set, as illustrated in figure 14 for the Edwards–Anderson Ising spin glass. We find that the measured round-trip times follow a fat-tailed Fr´echet distribution [29]^6 (solid line, fit performed with R^7 ). The integrated generalized extreme value distribution is given by:
Hξ;μ;β(τ ) = exp
[
1 + ξ
τ − μ β
) 1 /ξ ]
. (16)
(^6) The deviations of the data from the fitting function for large average round-trip times can be ascribed to the limited statistics used. (^7) R Core Team (R Manuals), http://cran.r-project.org
J. Stat. Mech. (2006) P
Figure 14. Distribution of average round-trip times for 5000 different samples of the 3D Edwards–Anderson Ising spin glass with Gaussian disorder and fixed system size L = 4 in the temperature range from Tmin = 0.10 to Tmax = 2.0. The data follow a fat-tailed Fr´echet distribution (solid line) with a shape parameter ξ = 035(1). The inset shows the shape parameter ξ as a function of system size L. Already for L 5, the shape parameter becomes ξ 1 /2, indicating that the variance of the distribution is no longer properly defined. The simulations have been performed using a fixed temperature set with M = 27 temperature points distributed such that, on average, replica swap moves have a nearly flat acceptance rate.
glass systems. Specifically, the finite-size scaling should be sensitive to such systematic errors, as it has been observed that the number of ‘hard’ samples significantly increases with system size; see the inset in figure 14 and [2, 30], and [31]. On the other hand, we can ask whether we can optimize the simulated temperature set and generate a ‘common’ temperature set for samples from the various parts of the distribution. To do so, we apply the feedback optimization outlined above in such a way that we generate a common probability distribution ¯η(T ) for a set of samples that are each characterized by their own diffusivity Di(T ), steady-state fraction fi(T ) and current ji, which are related by
ji = Di(T )¯η(T )
dfi dT
where the index i indicates the samples in the given set. Rearranging this equation, the local diffusivity of an individual sample can be expressed as
Di(T ) =
ji η ¯(T ) · dfi/dT
J. Stat. Mech. (2006) P
To ensure equilibration of all samples we want to simulate each sample for a fixed number of round-trips, despite the strong sample-to-sample variations. In order to minimize the overall computer time spent on simulating such a set of samples, we then minimize the sum of round-trip times
i τi. This is equivalent to minimizing the sum of the inverses of all currents, i.e., j i− 1 , since the current ji is inversely proportional to the round-trip time τi. Following a line of argument similar to that for the original temperature set optimization, we find that the optimal common temperature distribution ¯ηopt(T ) is proportional to the square root of the sum of inverse local diffusivities
η¯opt(T ) ∝
i
Di(T )
Again we can use a feedback loop to find an optimized common temperature set by feeding back the measured local diffusivities
η¯opt(T ) = C
η ¯(T )
i
τi ·
dfi dT
where C is a normalization constant. The common optimized temperature set is then found using a partial integration as given in equation (11).
- Conclusions
We have introduced an approach for systematically optimizing temperature sets for parallel tempering Monte Carlo simulations using an adaptive feedback method that is motivated by a recently developed ensemble optimization technique for broad-histogram Monte Carlo simulations [1]. We have applied the method to two paradigmatic spin models, the ferromagnetic Ising model and the fully frustrated Ising model in two dimensions. For both models we have shown that the feedback technique improves sampling of the phase space by reducing the average round-trip time of replicas diffusing in temperature space. Probably our most important result is the insight that the common wisdom that temperature sets in parallel tempering Monte Carlo should yield temperature-independent acceptance probabilities for the swap moves is not necessarily an optimal choice. Our feedback algorithm shifts temperature points in the optimized temperature sets towards the bottlenecks of the simulation, typically in the vicinity of a phase transition, where equilibration is suppressed. In particular, this has the effect that the acceptance probabilities for replica swap moves are higher around the bottlenecks and are not temperature independent for the so-optimized temperature sets. We also outline an approach to define ‘common’ temperature sets for systems that require configurational averages, such as spin glasses. In addition, we have briefly discussed the effects of sample-to-sample fluctuations with respect to equilibration times of individual spin glass samples, thus pointing towards a potential source of systematic errors in previous numerical studies of spin glasses. Clearly a deeper analysis of feedback-optimized parallel tempering Monte Carlo is needed, in particular in the context of spin glasses, as well as the questions raised at the end of section 4.2.
