Irreversible simulated tempering Yuji Sakai1,∗ and Koji Hukushima1,2,† 1Graduate School of Arts and Sciences, The University of Tokyo, 3-8-1 Komaba, Meguro-ku, Tokyo 153-8902, Japan 2Center for Materials Research by Information Integration, National Institute for Materials Science, 1-2-1 Sengen, Tsukuba, Ibaraki 305-0047, Japan (Dated: January 27, 2016) 6 1 An extended ensemble Monte Carlo algorithm is proposed by introducing a violation of the de- 0 tailed balance condition to the update scheme of the inverse temperature in simulated tempering. 2 Our method, irreversible simulated tempering, is constructed based on the framework of the skew detailed balance condition. By applying this method to the ferromagnetic Ising model in two di- n mensions on a square lattice as a benchmark, the dynamical behavior of the inverse temperature a and an autocorrelation function of the magnetization are studied numerically. It is found that the J relaxation dynamics of the inverse temperature changes qualitatively from diffusive to ballistic by 6 violating the detailed balance condition. Consequently, the autocorrelation time of magnetization 2 is several times smaller than that for the conventional algorithm satisfying the detailed balance condition. ] h PACSnumbers: 02.50.-r,05.10.Ln,02.70.Tt,05.70.Ln c e m I. INTRODUCTION In this study, we propose to apply the idea of a skew - t detailedbalancecondition(SDBC)[14,15]totheupdate a scheme of the inverse temperature in simulated temper- t SinceMetropolisetal.inventedaMarkov-chainMonte s ing, thus conducting the simulated tempering algorithm Carlo (MCMC) method in 1953 [1], it has been widely . t without DBC. As a benchmark, we examine the effi- a implementedinvariousresearchfieldsto evaluateexpec- ciency of our proposed algorithm in a two-dimensional m tationvaluesforahigh-dimensionalprobabilitydistribu- Ising model and show, numerically, that SDBC changes, tion. Meanwhile,someimprovementanddevelopmentof - d the MCMC method have been made for more efficient qualitatively,therelaxationdynamicsoftheinversetem- n sampling. Among them, simulated tempering [2, 3], de- perature in simulated tempering. Furthermore, we ob- o serve that the autocorrelationtime of the magnetization velopedinthefieldofstatisticalphysics,iscategorizedas c is also reduced by the violation of DBC. an extended ensemble method. In simulated tempering, [ the inverse temperature in the Gibbs-Boltzmann distri- This paper is organized as follows. Sec. II intro- 2 bution is treated as a random variable as well as the duces the simulated tempering method satisfying DBC. v configurations and thus the state space of the system is In Sec. III, a simulated tempering algorithmwith SDBC 6 extendedbyaddingthetemperature. AMarkovchainon is constructed. We apply the proposed algorithm to an 8 2 the extended state space is constructed with a detailed Ising model in two dimensions as a benchmark and con- 4 balance condition (DBC), which is a sufficient condition firm the efficiency of our proposed algorithm in Sec. IV. 0 for MCMC. The simulated tempering has been eagerly Section V summarizes the present work. . used for various problems in statistical physics [4–8] and 1 0 protein-folding problems [9, 10]. 6 Recently, a lifting technique in which the detailed bal- II. SIMULATED TEMPERING 1 ancecondition(DBC)intheMarkovchainisbrokenwith : v theglobalbalanceconditionstillholdinghasbeenexten- Inthissection,weexplaintheoutlineofsimulatedtem- i sively studied. Several studies have shown that diffusive X pering [2, 3] in order to fix our notation. relaxation dynamics in a one-dimensional random walk r with DBC is qualitatively improved by using the lifting a technique [11–14]. The transition graph of the inverse temperatureunderafixedconfigurationinthesimulated A. Setup temperingisthesameasthatintherandomwalk. Thus, we expect that the violation of DBC makes the relax- Let X be a configuration to be sampled from a tar- ation dynamics of the inverse temperature in simulated get distribution function in MCMC simulations. In sta- tempering change qualitatively. tistical physics, the target distribution P (X) is often β given by the Gibbs-Boltzmann distribution with an inverse temperature β, 1 ∗ yuji0920@huku.c.u-tokyo.ac.jp P (X)= exp[ βE(X)], (1) β † hukusima@phys.c.u-tokyo.ac.jp Z(β) − 2 where E(X) is a model Hamiltonian and Z(β) is the probabilities is often used. Here, we assume that the partition function of the model. In simulated temper- transition probabilities are decomposed into ing [2, 3], the (inverse) temperature, as well as the configuration, is a random variable. More specifically, β T(X′,βr X,βr) | takes R different values determined before the simula- =q(X′,β X,β )W(X′,β X,β ), (5) r r r r tion, expressed as β ,...,β . Thus, a state is specified | | 1 R by these variables, denoted by (X,βr). Then, an ex- and tended equilibrium distribution P (X,β ) for finding a ST r state (X,β ) is given as T(X,β X,β ) r l r | =q(X,β X,β )W(X,β X,β ), (6) P (X,β )= 1 exp[ β E(X)+g ], (2) l| r l| r ST r r r Z − ST where q denotes the proposal probability and W is the where the weight factor g is a constant depending only acceptance probability. Then, the explicit forms of the r on the inverse temperature and the extended partition Metropolis-Hastingstype of the acceptance probabilities function Z is for Eq. (5) and (6) are given by ST R R W(X′,β X,β ) r r ZST = exp[−βrE(X)+gr]= Z(βr)egr. (3) | q(X,βr X′,βr)PST(X′,βr) Xr=1XX Xr=1 =min(cid:20)1,q(X′,β|X,β ) P (X,β )(cid:21), (7) r r ST r | Foragivenβ ,theprobabilityforfindingaconfiguration r X inEq.(2)isproportionaltothatinEq.(1). Theaver- and age over the sampled configurations conditioned with β is equivalent to the equilibrium average at that temperr- W(X,βl X,βr) | ature. In contrast, the marginal probability for a given q(X,β X,β ) P (X,β ) r l ST l =min 1, | , (8) βr is obtained by (cid:20) q(X,βl X,βr)PST(X,βr)(cid:21) | Z(β )egr P (β )= P (X,β )= r . (4) respectively. ST r ST r XX ZST For simplicity, the set of inverse temperatures are or- dered such that β < β < < β . In addition, 1 2 R The marginal probability is independent of r when gr = throughout this paper we use ·t·h·e proposal probability a−tlnβrZ(oβfrt)h,ewmhiochdeilsspirmouploartteido.naIlntogetnheerbaul,lkitfriseeheanredrgtyo qqrr,,lr±≡1 =q(X1/,2βilf|X1,<βrr)0; φ(n) <ǫ . (19) relax ≡ { | β | } In this study, we apply the irreversible simulated tempering algorithm to the two-dimensional Ising model for 2 105 R MCSs after equilibration and mea- × × Figure 4 represents the R-dependence of the relaxation sure the empirical transition probabilities. In our algo- time τ (ǫ) with ǫ = 0.2. Although no difference be- rithm, explained in the previous section, the non-zero relax tween δ =0 and δ =0 is observed for small R in Fig. 4, component of the transition probability is T˜(βr,εβr,ε), the asymptotic beh6 avior of the relaxation time is quite T˜(β ,εβ ,ε), and T˜(β , εβ ,ε) with ε = .| Fig- r±1 r r r | − | ± differentineachcase. Inthecaseofδ =0,therelaxation ures5 and6 illustrate the β-dependence ofthe empirical time is asymptotically of order O(R2), which indicates transition probabilities. In Fig. 5, T˜(β ,+β ,+) has r+1 r | that the relaxation dynamics of the inverse temperature a dip aroundthe critical inverse temperature for smaller is diffusive. On the other hand, the relaxation time for values of R. However, the dip vanishes with increasing δ = 0 is asymptotically proportional to R, which indi- number of inverse temperatures and the empirical tran- 6 cates that the relaxation dynamics is ballistic. This dif- sition probability becomes flat with respect to β. Thus, ferenceintheasymptoticbehaviorisqualitativelyconsis- the empirical transition matrix can be approximated by tent with previous works on the one-dimensional simple that of the lifted simple random walk in one dimension random walk [11–14]. This shows that the violation of discussed in Ref. [14]. In addition, Fig. 6 shows that the DBC yields an accelerationof the relaxation of β. empirical transition probability for the ε flip is approx- 6 FIG. 5. (Color online) β-dependence of the empirical transition probability from (βr,+) to (βr+1,+) with (a) R = 32 in the left, (b)R=128 in thecenter,and(c) R=512 intherightpanel. Thelinearsize oftheIsing-spinsystem ischosen asL=25 and the parameter δ, which characterizes the deviation from DBC, is chosen as δ = 0.0 (red square), δ = 0.3 (green circle), δ=0.6 (orange triangle), and δ=0.9 (blueinverted triangle), respectively. FIG.6. (Color online) β-dependenceof theempirical transition probability from (βr,+)to(βr,−). Theleft, center,and right panels represent the case (a) R = 32, (b) R = 128, and (c) R = 512, respectively. The linear size of the Ising-spin system is chosenasL=25 andtheparameterδ,whichcharacterizesthedeviationfromDBC,ischosenasδ=0.3(greencircle), δ=0.6 (orange triangle), and δ=0.9 (blueinverted triangle), respectively. Note that thevertical axis is log-scale. imately proportional to R−1. In Ref. [14] it was shown celerationof the relaxationdynamics of the inverse tem- analyticallythattheliftedsimplerandomwalkinonedi- peraturebytheviolationofDBCisnumericallyandthe- mensionwithO(R−1)ε-flipprobabilityfollowsaballistic oretically confirmed. process. Thus, these results imply the reduction of the convergence rate of the inverse temperature. To evaluate the convergence rate from the empirical transition matrix, we define the relaxation time as fol- C. Autocorrelation function lows. Let λ (k = 1,2,...,2R) denote the eigenval- k ues of the empirical transition matrix defined as T˜ = The acceleration of the relaxation of β is expected to T˜(β ,ε′ β ,ε) R2R×2R. Without loss promote the accelerationof the relaxationof the magne- l r (cid:16) | (cid:17)1≤r,l≤R;ε,ε′=± ∈ tizationintheIsingmodel. Inthissubsection,weobserve of generality, the eigenvalues are aligned as 1 = λ > 1 the time evolution of the autocorrelationfunction of the λ λ . Then, the relaxation time of the | 2| ≥ ··· ≥ | 2R| magnetizationwhose initial state is preparedas an equi- inverse temperature is defined as librium state at β , the lowest temperature in our sim- R 1 ulations. Let m = N S /N denote the averaged spin τ˜ . (21) i=1 i relax ≡−ln|λ2| and let h···iβ be thPe expectation value with respect to the Gibbs-Boltzmann distribution in Eq. (1). Then, we Figure 7 illustrates the R-dependence of the relaxation define the (normalized) autocorrelation function of the timeτ˜ obtainedbynumericallydiagonalizingtheem- relax magnetization in simulated tempering as piricaltransitionmatrix. AsshowninFig.7,theasymp- totic behavior of the relaxation time in both the re- m(0)m(n) m m versible and irreversible cases is compatible with the re- C(n) h i−h iβRh ieq, (22) sults obtained in the previous subsection. Thus, the ac- m ≡ hm2iβR −hmiβRhmieq 7 107 1 6 δ = 0.0 O(R2) 10 δ = 0.3 0.8 105 δ = 0.6 x 104 δ = 0.9 n) 0.6 ~τrela 103 O(R) (Cm 0.4 δδ == 00..03 102 0.2 δ = 0.6 δ = 0.9 1 10 0 0 10 0 1 2 3 4 5 0 1 2 3 4 10 10 10 10 10 10 10 10 10 10 10 n [MCS] R FIG. 8. (Color online) Time evolution of the autocorrela- FIG. 7. (Color online) R-dependence of the relaxation time tion function of the magnetization in the simulated temper- of β determined by the empirical transition matrix. The lin- ing with theSDBC. The initial state S(0) is prepared by the ear size of the system is L= 25 and the parameter δ, which characterizes the deviation from DBC, is chosen as δ = 0.0 Metropolis-Hastings algorithm at β = βR = 0.5. The linear size of the system and the numberof temperature points are (redsquare),δ =0.3(greencircle), δ=0.6(orangetriangle), chosen as L = 25 and R = 29, respectively. The parameter and δ = 0.9 (blue inverted triangle), respectively. The dot- δ, which characterizes the deviation from DBC, is chosen as ted anddashed linesrepresent theexpectedasymptoticform δ = 0.0 (red square), δ = 0.3 (green circle), δ = 0.6 (orange proportional to R2 and R,respectively. triangle), and δ = 0.9 (blue inverted triangle), respectively. Each data point is averaged over 210 samples and the error bars are of the order of thesymbol sizes. where m(0)m(n) denotes the sample averageof the cor- h i relation between the initial averaged spin and that after n MCSs. Figure 8 illustrates the time evolution of Cm(n), 107 which indicates that the violation of DBC reduces the 6 δ = 0.0 10 relaxation rate in the autocorrelation function by a fac- δ = 0.3 toraslargeasten(10). Thisreductionis affectedbythe 105 δ = 0.6 acceleration of the relaxation of β. 4 δ = 0.9 10 In order to evaluate quantitatively the improvement rr o of the relaxation dynamics of the autocorrelation, the τc 103 autocorrelationtime is defined as 2 10 τ (ǫ) inf(n>0; C(n) <ǫ), (23) 1 corr ≡ | m | 10 0 and especially τ = τ (ǫ =0.2). Figure 9 represents 10 corr corr 0 1 2 3 4 the R-dependence of the autocorrelation time τ . As 10 10 10 10 10 corr shown in Fig. 9, while there is no difference observed R betweenδ =0andδ =0forsmallR,theautocorrelation 6 timeisimprovedforrelativelylargeR. Theaboveresults FIG. 9. (Color online) R-dependence of the autocorrelation confirmthattheviolationofDBCimprovestheefficiency time in simulated tempering with SDBC. The initial state ofthesimulatedtemperingalgorithmwithrespecttothe S(0) ispreparedbytheMetropolis-Hastingsalgorithmatβ= sampling of both β and X. βR = 0.5. The value L = 25 is used and the parameter δ, which characterizes the deviation from DBC, is chosen as δ = 0.0 (red square), δ = 0.3 (green circle), δ = 0.6 (orange V. SUMMARY AND DISCUSSION triangle), and δ = 0.9 (blue inverted triangle), respectively. Each data point is averaged over 210 samples and the error Wehaveconstructedanirreversiblesimulatedtemper- bars are of the order of thesymbol sizes. ing algorithm by introducing the lifting technique based onthemethodologyofSDBCtotheupdateschemeofthe inverse temperature. Benchmarks for the Ising model the efficiency of extended-ensemble methods. Further- show that our algorithm accelerates the relaxation dy- more, we consider the empirical transition probability namics of inverse temperature and the autocorrelation with respect to the inverse temperature and the lifting function of the magnetization compared to the tradi- parameter to investigate the relaxation dynamics of the tional simulated tempering algorithm based on DBC. inverse temperature in detail. It is easily measured dur- Theseresultsshowthattheliftingtechniquecanimprove ing numerical simulations in the irreversible simulated 8 tempering algorithm. We found that the empirical tran- by an exact numerical method. The choice of the set of sitionmatrixisapproximatelythesameasthetransition inverse temperatures β and parameters g affect the r r { } matrixoftheliftedsimplerandomwalkinonedimension efficiency of simulated tempering. Several studies have discussedin Ref. [14]. Thus, it is theoretically confirmed proposed their efficient choices [18, 26–28]. A promising that the lifting technique accelerates the relaxation dy- way for estimating the weight factor is to implement the namics of the inverse temperature. irreversible simulated tempering algorithm which is our AlthoughweusedourproposedalgorithmfortheIsing current work in progress [29]. model in two dimensions in this paper, our algorithm is, inprinciple,applicabletoanyothersystemsuchasPotts model, Heisenberg spin glass, and protein systems. It ACKNOWLEDGMENTS is also possible to combine other update schemes of the configuration of target systems, such as the Swendsen- Wang algorithm and the Wolff algorithm, instead of The authors are grateful to S. Todo for useful com- the Metropolis-Hastingsalgorithm. Our algorithmcould ments and for bringing the method of Ref. [25] to our takeovertheseadvantagesfromthetraditionalsimulated notice. Y.S. is supported by a Grant-in-Aid from the tempering method. Itis worthinvestigatingwhetherthe Japan Society for Promotion of Science (JSPS) Fellows irreversible simulated tempering combined with such an (Grant No. 267868). K.H. is supported by Grants-in- · update scheme workseffectively in a system with a first- Aid for Scientific Research from MEXT, Japan (Grant order phase transition and spin glasses. Nos. 25610102 and 25120010), and JSPS Core-to-Core In this study, all inverse temperatures were arranged program “Nonequilibrium dynamics of soft matter and atequaldistancesandtheweightfactorg wasestimated information.” r [1] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, [15] K.S.Turitsyn,M.Chertkov,andM.Vucelja,PhysicaD A. H. Teller, and E. Teller, J. Chem. Phys. 21, 1087 240, 410 (2011). (1953). [16] W. Hastings, Biometrika 57, 97 (1970). [2] A. P. Lyubartsev, A. A. Martinovski, S. V. Shevkunov, [17] C. J. Geyer and E. A. Thompson, J. Am. Stat. Assoc. and P. N. Vorontsov-Velyaminow, J. Chem. Phys. 96, 90, 909 (1995). 1776 (1992). [18] A. Mitsutake and Y. Okamoto, Chem. Phys. Lett. 332, [3] E. Marinari and G. Parisi, Europhys. Lett. 19, 451 131 (2000). (1992). [19] A.MitsutakeandY.Okamoto,J.Chem.Phys.121,2491 [4] E. Vicari, Phys.Lett. B 309, 139 (1993). (2004). [5] W.KerlerandP.Rehberg,Phys.Rev.E50,4220(1994). [20] Y. Li, V. A. Protopopescu, and A. Gorin, Phys. Lett. A [6] B. Coluzzi, J. Phys.A: Math. Gen. 28, 747 (1995). 328, 274 (2004). [7] L. A. Fernańdez, E. Marinari, and J. J. Ruiz-Lorenzo, [21] P. H. Nguyen, Y. Okamoto, and P. Derreumaux, J. Phys. I France 5, 1247 (1995). J. Chem. Phys. 138, 061102 (2013). [8] M. Picco and F. Ritort, Physica A 250, 46 (1998). [22] Y. Mori and H. Okumura, J. Comput. Chem. 36, 2344 [9] A. Irba¨ck and F. Potthast, J. Chem. Phys. 103, 10298 (2015). (1995). [23] H. Suwa and S. Todo, Phys. Rev. Lett. 105, 120603 [10] A. Irba¨ck, C. Peterson, and F. Potthast, Phys. Rev. E (2010). 55, 860 (1997). [24] L. Onsager, Phys. Rev.65, 117 (1944). [11] F. Chen, L. Lovasz, and I. Pak, In: Proc. 17th Annual [25] B. Kastening, Phys. Rev.E 64, 066106 (2001). ACM Symposium on Theory of Computing, 275 (1999). [26] U.H.E.HansmannandY.Okamoto,J.Comput.Chem. [12] P.Diaconis,S.Holmes,andR.M.Neal,TechnicalReport 18, 920 (1997). BU-13850-M,BiometricsUnit,CornellUniversity(1997); [27] S.ParkandV.S.Pande,Phys.Rev.E76,016703(2007). Ann.Appl.Probab. 10, 726 (2000). [28] A. Valentim, M. G. E. da uz, and C. E. Fiore, Com- [13] M. Vucelja, arXiv:1412.8762 (2014). put. Phys.Commun. 185, 2046 (2014). [14] Y.Sakai and K. Hukushima,arXiv:1511.08100 (2015). [29] Y. Sakai and K. Hukushima,in preparation.

Irreversible simulated tempering

Similar Irreversible simulated tempering