Irreversible simulated tempering

0.58 MB·

by  Yuji Sakai

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 in-
verse temperature β,
1
∗ [email protected] P (X)= exp[ βE(X)], (1)
β
† [email protected] 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 con-
figuration, 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)<givRen, abnydqz1e,r2o=othqRer,Rw−is1e=(Fi1g.a1n)d.
estimate the value of the free energy for a statistical- Note that q is independent of the configuration X.
r,l
mechanical model. However, if it could be estimated a Then, the procedure of the simulated tempering method
priori,evenapproximately,auniformsamplingforβfrom is described as follows.
high to low temperatures can be put into practice. By
considering an appropriate Markov chain, the state of 1. Arbitrarily chose an initial state (X(0),β(0)).
β wanders on the temperature axis in a random-walk
r
manner. One may expect that it is relatively easy to 2. Iterate the update trials for an original configura-
sample configurations at sufficiently high temperatures, tion X at a fixed β according to the conventional
which would help an efficient sampling at low tempera- Metropolis-Hastings method:
tures through the random walk of β . This is what we
r (a) Suppose that the current state is (X,β ) and
expect to perform in simulated tempering. r
select a configurationX′ with the probability
q(X′,β X,β ).
r r
|
B. Simulated tempering algorithm with detailed (b) Accept the new state (X′,β ) with the prob-
r
balance condition ability W(X′,β X,β ). If it is rejected, set
r r
|
the current state as the next state.
An explicit update scheme of the simulatedtempering
method consists of the following two steps: an update of 3. Iteratetheupdatetrialsforaninversetemperature
a configurationX for a fixed β and an update of β for β according to the following procedure:
r r
a fixedX. Inorderto generate a Markovchain,the cor-
(a) Suppose that the current state is (X,β ) and
respondingtwotransition-probabilitymatricesareintro- r
duced. One is the transitionmatrixfroma state (X,β ) choose βl with the probability qr,l as a candi-
r
to(X′,β )denotedasT(X′,β X,β ). Theotherisone dateforthenextinversetemperature(Fig.1).
r r r
from (X,β ) to (X,β ) as T(X|,β X,β ). They satisfy (b) Accept the next state (X,β ) with the proba-
r l l r l
DBC for the stationary distribution| of Eq. (2). In prac- bilityW(X,β X,β ). Ifitisrejected,setthe
l r
|
tice, the Metropolis-Hastings type [16] of the transition current state as the next state.
3
FIG. 1. Graphical expression of the proposal probabilities qr,l in conventional simulated tempering.
Although any inverse temperature could be chosen as system as a lifting parameter, the state space is dupli-
a candidate in step 3, in practice we only consider the cated. A state in the duplicated state space is denoted
nearest ones to increase the transition rate. By repeat- (X,β ,ε). Then, the extended equilibrium distribution
r
ing steps 2 and 3, our desired Markov chain of the state P (X,β ,ε) for finding a state (X,β ,ε) is given by
ST r r
(X,β) is obtained.
1
Since the simulated tempering was invented, it has P (X,β ,ε)= exp[ β E(X)+g ]. (9)
ST r r r
been widely applied to various problems [4–10]. In addi- 2ZST −
tion, the improvement of simulated tempering has been
Note that the marginal probability for a given configu-
continuously studied [17–22]. In particular, an efficient ration X and an inverse temperature β is exactly the
r
simulated tempering in which DBC breaks and the re- same as P (X,β ) and that it is uniform for a given ε.
ST r
jectionrateis decreasedwiththe violationofDBC using
Inthisstudy,weapplythemethodologyofSDBConly
the Suwa-Todo algorithm [23] was proposed recently in
to the update scheme of the inverse temperature. The
Ref. [22]. Others authors keep DBC unbroken. The typ-
skew detailed balance condition in this context is ex-
ical relaxation time of the algorithm is reduced by some
pressed as
factor, but its temperature dynamics does not change
quantitatively from standard diffusive dynamics. T(X,β ,+X,β ,+)P (X,β ,+)
l r ST r
|
=T(X,β , X,β , )P (X,β , ), (10)
r l ST l
−| − −
III. SIMULATED TEMPERING WITH SKEW where T(X,β ,εX,β ,ε) denotes the transition prob-
l r
DETAILED BALANCE CONDITIONS ability from a sta|te (X,β ,ε) to (X,β ,ε). Again, we
r l
decomposethetransitionprobabilityT(X,β ,εX,β ,ε)
l r
|
In the conventional simulated tempering method ex- into the product of a proposalprobability and an accep-
plained in the previous section, the transition graph of tance probability expressed as
the inverse temperature β under fixed X in Fig. 1 is the
T(X,β ,εX,β ,ε)
same as that of a one-dimensional simple random walk. l r
|
Thus, β behaves as a random walker on the graph. It =q(X,βl,εX,βr,ε)W(X,βl,εX,βr,ε). (11)
| |
is known that a random walk satisfying DBC is essen-
Byusingtheproposalprobabilityq(X,β ,εX,β ,ε),the
tially a diffusive process and its relaxation time is of or- l r
|
der O(Ω2), where Ω denotes the number of states in the generalformofthe Metropolis-Hastings-typeacceptance
probability that satisfies SDBC is explicitly given by
random walk. Several studies have shown, numerically
andanalytically,thattheΩ-dependenceoftherelaxation
W(X,β ,εX,β ,ε)
timeisimprovedbyintroducingthe“lifting”techniqueto l | r
the randomwalk inone dimension[11–14]. Here DBC is q(X,βr, εX,βl, ε)PST(X,βl)
=min 1, − | − , (12)
brokenbyaddingaliftingparameterwhilepreservingthe (cid:20) q(X,β ,εX,β ,ε) P (X,β )(cid:21)
l r ST r
|
globalbalancecondition. Thisstronglysuggeststhatthe
In this study, we construct the proposal probability as
performanceofsimulatedtemperingisimprovedqualita-
tively by the lifting of the update of β. In this section, q(ε), which is independent of the configuration X as fol-
r,l
weapplythemethodologyofSDBC[14,15]tosimulated lows (Fig. 2):
tempering,especiallytotheupdateschemeoftheinverse
temperature. The proposedmethod is called irreversible q(ε) =q(ε) =1, (13)
1,2 R,R−1
simulated tempering.
1 δε
q(ε) = ± , (14)
r,r±1 2
A. Setup if 1 < r < R, and q(ε) = 0 otherwise. The parameter
r,l
δ controls the violation of DBC and satisfies δ < 1.
| |
Let us reconsider the setup of conventional simulated DBC is restored when the parameter δ is set to zero.
tempering to extend it to an irreversible one. By intro- One expects that a finite positive value of δ enhances a
ducing an auxiliary random variable ε +, to the clockwise flow in the dynamics of β in Fig. 2.
∈ { −}
4
FIG. 2. Graphical expression of the proposal probabilities q(ε) in an irreversible simulated tempering.
r,l
B. Simulated tempering algorithm with skew a benchmark and numerically evaluate the efficiency of
detailed balance conditions the algorithm. Let L denotes the linear size of the Ising
modelandletN =L2. Theenergyfunctionofthemodel
For simplicity, we rewrite the acceptance probability is given by
W(X,β ,εX,β ,ε) as W(ε). Then, the update scheme
of an invler|se temrperaturerβ,l and an auxiliary variable ε E(S)= SiSj, (17)
−
is described as follows: Xhiji
1. Iteratetheupdatetrialsforaninversetemperature with S = S N and S = 1. The sum is taken over
β according to the following procedure: { i}i=1 i ±
the nearest neighbor pairs, imposing a periodic bound-
(a) Suppose that the current state is (X,β ,ε) arycondition. Theenergyunitissettounity. Thesetof
r
and choose β with the probability q(ε) as a inverse temperatures in the simulated tempering is pre-
l r,l pared from β = 0.2 and β = 0.5, with the interme-
candidate of the next inverse temperature. 1 R
diate values equally spaced between them. Note that
(b) Acceptthenextstate(X,βl,ε)withtheprob- the critical inverse temperature of the model is known
ability W(ε). as β = ln(1+√2)/2 0.4407 [24], which is inside the
r,l c ≃
temperature region in our simulations.
(c) If the trial is rejected, flip ε and set
(X,β , ε) as the next state with the prob- In the present work, we focus our attention on the
r
− performance of irreversible simulated tempering for an
ability
idealweightparameterg. Thus,theparameterg ,which
r
Λ(ε) should be given a priori, is evaluated by using an exact
λ(ε) r , (15)
r ≡ 1 q(ε)W(ε) free energy numerical method with a polynomial time of
− r,l r,l N [25]availableforfinite-sizeIsingmodelsintwodimen-
Xl6=r sions. The determination of the parameter, which is an
where importantissue whenactually applying the algorithmto
some statistical-mechanics models, is discussed in a sep-
Λ(ε) =max0, ε ε′q(−ε′)W(−ε′). (16) arate paper [29].
r − r,l r,l
εX′=±Xl6=r
 
If it is also rejected, set the current state as A. Relaxation dynamics of the inverse temperature
the next state.
We study the dynamics of the inverse temperature in
2. Update the configurationX with the conventional
simulated tempering. In this work, a Monte Carlo step
Metropolis-Hastingsmethodattheinversetemper-
(MCS) is defined by the time unit where N spin trials
ature β , as explained in Sec. II.
l and a trial of β are performed. Figure 3 illustrates the
Itisstraightforwardtoverifythattheglobalbalancecon- time evolution of the relaxation function of β defined as
dition is satisfied in the above procedure [14, 15].
β(n) β
φ(n) h i−h ieq, (18)
β ≡ β(0) β
eq
IV. BENCHMARK h i−h i
where β(n) denotes the sample average of the inverse
h i
In this section, we apply the proposed algorithm to temperature after n MCS. The expectation with respect
an Ising model in two dimensions on a square lattice as tothetargetdistributioninEq.(2)isdenotedby .
eq
h···i
5
1 107
6 δ = 0.0 O(R2)
10
0.8 δ = 0.3
105 δ = 0.6
0.6 4 δ = 0.9
) x 10
(nβ 0.4 δ = 0.0 ela 3 O(R)
φ δ = 0.3 τr 10
2
0.2 δ = 0.6 10
δ = 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.3. (Coloronline)Timeevolutionoftherelaxation func- FIG.4. (Coloronline)R-dependenceoftherelaxationtimeof
tionoftheinversetemperatureβinsimulatedtemperingwith βinsimulatedtemperingwithSDBC.TheinitialstateS(0) is
SDBC, applied to the Ising model in two dimensions with prepared by the Metropolis-Hastings algorithm at β =βR =
L=25. The initial state S(0) is prepared by the Metropolis- 0.5. Here, L = 25 and the parameter δ, which characterizes
Hastingsalgorithmatβ=βR =0.5. Thenumberofobserved thedeviationfromDBC,ischosenasδ =0.0(redsquare),δ=
temperatures is R = 29. The parameter δ, which character- 0.3(greencircle),δ=0.6(orangetriangle),andδ=0.9(blue
izesthedeviationfromDBC,ischosenasδ=0.0(redsquare), inverted triangle). The average for each point is taken over
δ =0.3 (green circle), δ =0.6 (orange triangle), and δ =0.9 210 samples andtheerrorbarsareoftheorderofthesymbol
(blue inverted triangle), respectively. The history average is sizes. The dotted and dashed line represent asymptotic lines
taken over210 samples and the error bars are of the order of proportional to R2 and R,respectively.
thesymbol sizes.
B. Empirical transition matrix of inverse
temperature
The initial conditions of β are set as β(0) =β . The re-
R
The dynamics of the inverse temperature is explored
laxation function monotonically decays to zero with re-
also through another quantity. In implementing the ir-
laxationtime. Figure3indicatesthatthe convergenceof
reversible simulated tempering algorithm, one can easily
β withδ =0 is morethan10times fasterthan thatwith
6 measuretheempiricaltransitionprobabilitywithrespect
δ =0 andthat the largerδ is, the more the relaxationof
to the inverse temperature β and the auxiliary random
β is accelerated.
variableε. The empiricaltransitionprobability fromthe
Inordertoevaluatequantitativelytheimprovementof state (βr,ε) to (βl,ε′) is defined as
the relaxation dynamics of the inverse temperature, we
T˜(β ,ε′ β ,ε)
measure the relaxation time defined as l r
|
# of transitions from (β ,ε) to (β ,ε′)
r l
. (20)
≡ # of visit (β ,ε)
r
τ (ǫ) inf n>0; φ(n) <ǫ . (19)
relax ≡ { | β | } In this study, we apply the irreversible simulated tem-
pering 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´ndez, 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.