Table Of ContentADP-11-01/T723
A Novel Multiple-Time Scale Integrator for the
Hybrid Monte Carlo Algorithm
Waseem Kamleh
SpecialResearchCentrefortheSubatomicStructureofMatterandDepartmentofPhysics,
1
UniversityofAdelaide5005,Australia.
1
0
2 Abstract. HybridMonteCarlosimulationsthatimplementthefermionactionusingmultipleterms
n arecommonlyused.Bythenatureoftheirformulationtheyinvolvemultipleintegrationtimescales
a intheevolutionofthesystemthroughsimulationtime.Thesedifferentscalesareusuallydealtwith
J by the Sexton-Weingarten nested leapfrog integrator. In this scheme the choice of time scales is
4 somewhat restricted as each time step must be an exact multiple of the nextsmallest scale in the
sequence.Anovelgeneralisationofthenestedleapfrogintegratorisintroducedwhichallowsforfar
] greaterflexibilityinthechoiceoftimescales,aseachscalenowmustonlybeanexactmultipleof
t
a thesmalleststepsize.
l
- Keywords: LatticeQCD,algorithms,HMC
p
PACS: 11.15.Ha
e
h
[
INTRODUCTION
1
v
1
Hybrid Monte Carlo(HMC)[1] is the algorithm of choice for generating lattice gauge
5
field configurations that include the effects of fermion loops. The main expense in
6
0 such simulationsis evaluatingthe contributionfrom the fermion determinant. A variety
.
1 of improvements to the basic HMC algorithm have been developed in an effort to
0 ameliorate this expense. Many of these variants involve the introduction of additional
1
terms into the fictitious Hamiltonian. For example, Clark and Kennedy[2] introduce
1
: multiple pseudofermion fields into the Rational HMC algorithm, each contributing a
v
i fractionofthefermiondeterminant.Alternatively,Hasenbuch[3]usesthefermionmatrix
X
(withaheavierquarkmass)asapreconditionertothedesiredfermionaction.Thiscauses
r
a thefermionaction to besplitinto twoparts. Anotherexampleis thepolynomialfiltered
HMCalgorithm[4],whichuses shortpolynomialapproximationstotheinversefermion
matrixin orderto separate thedynamics ofthe Hamiltonian.Multiplelevelsof filtering
can be introduced, and this technique is also applicable to single flavour simulations,
allowing the number of terms to be simulated in the Hamiltonian to grow to five or
more. The Schwarz preconditioning technique employed by Luscher[5] also makes use
ofa hierarchyoftimescales.
The traditional way of dealing with a multiple time scale integration is to use the
nestedleapfrogalgorithmbySextonandWeingarten[6].Thisschemeisquiterestrictive
however,asbeginningwiththefinesttimescale,thestep-sizeforthenext(coarser)time
scale must be an exact multipleof the previous scale. This may prevent one from being
able to choose the most efficient set of integration parameters, particularly if there are
many time scales. Here we present an improvement of this scheme, which allows for
greaterfreedomofchoice,aseachtimescalemustonlybeamultipleofthefinestscale.
Hybrid Monte Carlo Overview
In this section we provide a brief overview of the Hybrid Monte Carlo algorithm[1]
in order to provide a framework for the introduction of our novel integrator. We wish
to generate an ensemble {U} of representative gauge fields distributed according to
i
the probability distribution r (U) = e−S[Ui], where the effective action for full QCD
i
S[U] = S [U]+S [U] can be divided into two parts, the gauge action S [U] and the
G F G
fermion action S [U]. Here we have assumed that the fermionic degrees of freedom
F
havebeen integratedoutin theusualway.
In the Hybrid Monte Carlo algorithm, the quantum lattice field theory is embedded
in a higher-dimensional classical system through the introduction of a fictitious (simu-
lation) time. The gauge field U is associated with its (fictitious) conjugate momenta P,
andtheclassicalsystemisdescribed by theHamiltonian,
H[U,P]=T[P]+S[U], (1)
where the fictitious kinetic energy is given by T[P]=(cid:229) x,m 12TrPm (x)2. Given a config-
uration U, a new gauge field U′ is generated by performing a HMC update U →U′,
whichconsistsoftwo steps:
(1) Molecular Dynamics Trajectory: Sample P from a Gaussian ensemble. Integrate
Hamilton’s equations of motion to deterministically evolve (U,P) along a phase
spacetrajectory to(U′,P′).
(2) Metropolis step: Accept or reject the new configuration (U′,P′) with probability
r (U →U′)=min(1,e−D H),D H =H[U′,P′]−H[U,P].
The discretised equations of motionare derived by requiring that the Hamiltonianbe
conserved along the phase space trajectory. We can express the equations of motion in
termsofthetimeevolutionoperatorsinduced bythekineticand potentialenergy terms.
The evolution operators that evolve the gauge field and its conjugate momenta forward
asimulationtimestep hare givenrespectivelyby
V (h):{U,P}→{Uexp ihP),P}, (2)
T
(cid:0) d S
V (h):{U,P}→{U,P−hU }. (3)
S d U
Note that we demand that the evolution must preserve the SU(3) property of the gauge
field.
We must then combine these evolution operators into an overall evolution operator
thatisreversible.Thesimplestsuchintegrationschemeistheleapfrog
h h
V (h)=V ( )V (h)V ( ). (4)
H S T S
2 2
After discretisation, for sufficiently small step sizes h, the integration will conserve
theHamiltonianup toO(h2).
Multiple time-scales in molecular dynamics integrators
Iftheactionand thustheHamiltonianissplitintotwoparts H andH ,
1 2
H =T[P]+S [U]+S [U] (5)
1 2
H H
1 2
| {z } |{z}
thenwedefineintegratorsforH and H asfollows
1 2
V (h)=V (h)V (h)V (h), V (h)=V (h). (6)
H1 S1 2 T S1 2 H2 S2
A compound integrator for the full Hamiltonian can be constructed by using a Sexton-
Weingartenscheme[6]:
m
h h h
V (h)=V ( ) V ( ) V ( ) (7)
H H H H
2 2 (cid:20) 1 m (cid:21) 2 2
wherem∈Z.Thisnestedleapfrogintegratoreffectivelyintroducestwotime-scalesinto
the evolution, h and h/m. Additional time scales may be introduced by repeating the
nestingprocedure.
RESULTS
Webeginbyassumingthat ourHamiltonianconsistsofat leastthreeterms,
H =T +S +S +..., (8)
1 2
where T is the “kinetic energy” due to the conjugate momenta, and the terms S imple-
i
mentthelatticeQCDaction.TypicallywewouldchooseS =S ,thatisS isthegauge
1 G 1
action. Then S ,S ,... will be the terms implementing the fermion action according to
2 3
our algorithm of choice, be it mass-preconditioning or polynomial filtering and so on.
Theonlythingweassumeaboutthefermionactionisthatitisimplementedsuchthatthe
onlyfieldsthatareupdatedduringtheintegrationarethegaugefieldU anditsconjugate
momentaP (e.g. usingpseudofermions).
For each scale we associate a timestep h and a corresponding integer N such that
i i
h = 1/N. We require that i = 1 corresponds to the scale at which the gauge field is
i i
updated.
AssumingthatN >N fori< j a nestedleapfrog algorithmthen requiresthat
i j
N|N ∀i>1. (9)
i i−1
Thatis,N mustbeadivisorofN ,orequivalently,h mustbeanexactmultipleofh .
i i−1 i i−1
This means, for example, that each successive scale must be at least twice as coarse as
the previous scale. It also may lead to being forced to choose a scale that is smaller
than the one desired for a given time scale in order to simultaneously control the finite
step-size errors as well as maintain the required arithmetic relation (9). The restriction
of having to choose successivedivisors for the various N may not be the most efficient
i
orflexibleway ofperformingthemoleculardynamicsintegration.
A generalised leap-frog integrator
We present a generalised integration scheme in which it is only required that the
integrationscalessatisfytherelation
N|N ∀i>1. (10)
i 1
Thisisofcourseequivalenttorequiringthat thestep sizeh isan exact multipleofh .
i 1
Inastandardleapfrogalgorithm,onealternatesbetweenupdatesV tothegaugefield
T
U andupdatesV totheconjugatemomenta.LetV denotetheupdatetoPcorresponding
S i
totheactionS .Now,astheguidebosonsareheldfixedduringanintegrationtheupdates
i
V only depend upon the gauge field. As the updatesV are additive to P, it follows that
i i
thedifferentV commute:
i
h h
i i
V( )V (h )V( )=V (h )V(h)=V(h)V (h ). (11)
i j j i j j i i i i j j
2 2
Definetheintegers
m =N ÷N (12)
i 1 i
tobetheratiosofthescales.Inordertoconstructourreversibleintegratorwefirstdefine
amap
V ifm|k
Q [V;m,k∈N]= (13)
(cid:26)I(theidentity)otherwise.
Letm bethelowestcommonmultipleof{m }, and leth bethesmallesttimestep (in
T i T
ourcase h =h ).Thenourintegratoris
T 1
V(h)=(cid:213) V(hi)×m(cid:213)T−1V (h ) (cid:213) Q [V(h);m,k] V (h )×(cid:213) V(hi), (14)
i T T i i i T T i
2 2
i k=1 n i o i
where h = m h is the total timestep taken by V. The above expression is straightfor-
T T
wardly implemented in code. We demonstrate this with a pseudocode implementation
here. Denote by {a ≡ b modm} the usual notion of congruence modulo m. Then we
can implementthegeneralised integratoras follows.
Pseudo-code forthe generalisedintegrator:
• Foreach term intheactionSi performan initialhalf-stepVi(21hi)updatingP.
• Loopover j=1to N−1
– ApplyV (h)to updateU.
T
– If{0≡ j modm } applyV(h )toupdateP
i i i
• ApplyVT(h)toupdateU.
• Foreach term intheactionSi performa final half-stepVi(21hi)updatingP.
The advantage of the generalised integrator is that it allows finer control over the
differentscales.Ananalysisofthefinite-stepsizeerrorsforthegeneralisedintegratoris
providedin thenextsection.
Integrator Error Analysis
Weperformanerroranalysisofourgeneralisedintegratorforasimplechoiceofstep-
sizes,followingtheprocedurein[6].GivenaHamiltonianH wecanwritetheevolution
operator for our system asV(h)=exp(hHˆ), with step-size h. Here we have defined Hˆ
as the linear operator on the vector space of functions f on phase space (p,q) defined
bythePoissonbracket
Hˆ f =−{H, f}=(cid:229) ¶ H ¶ f −¶ H ¶ f . (15)
(cid:18)¶ p ¶ q ¶ q ¶ p (cid:19)
i i i i i
IfwewritetheHamiltonianas
H =T +S +S +S +S +... (16)
1 2 3 4
then for each term in the Hamiltonian we can correspondingly define a linear operator
using the Poisson bracket relation above. We denote the linear operator associated with
agiventermin theHamiltonianbyaddinga“hat”to theappropriatesymbol.
Proceeding with the error analysis, we make use of the Baker-Campbell-Hausdorff
result,
l 3
el Aˆel Bˆel Aˆ =exp l (2Aˆ+Bˆ)+ ([[Aˆ,Bˆ],Aˆ]+[[Aˆ,Bˆ],Bˆ])+O(l 4) (17)
6
(cid:16) (cid:17)
and apply thisto thegeneralised leapfrog integratorin thesimplecase ofH =T +S +
1
S ,wherethetimescaleforeachterminH correspondstoanumberofintegrationsteps
2
N =6,N =3and N =2respectively.Theintegratorforthissimplestnon-trivialcase
T 1 2
can bewrittenas
VH(h)=eh4Sˆ2e6hSˆ1e3hTˆeh3Sˆ1eh6Tˆe2hSˆ2eh6Tˆe3hSˆ1eh3Tˆeh6Sˆ1eh4Sˆ2. (18)
Repeated application of our BCH result allows us to deduce that the above expression
can bewrittenas
1 1
V (h)=exp hHˆ +h3 [[Sˆ ,Tˆ],Tˆ]+ [[Sˆ ,Tˆ],Sˆ ]+
H 2 2 2
48 96
(cid:16)
(cid:0)
1 1 1
[[Sˆ ,Tˆ],Sˆ ]+ [[Sˆ ,Tˆ],Tˆ]+ [[Sˆ ,Tˆ],Sˆ ] (19)
1 1 1 1 2
216 108 36
(cid:17)
(cid:1)
Fromthisexpressionwecan immediatelyseethattheerrorinthegeneralised integrator
relativeto theleading termisO(h2),justas fortheregularleapfrog.
Ifweexaminetheindividualleapfrogintegratorscorrespondingto
H =T +S , H =T +S (20)
1 1 2 2
weobtain
VH (h)=eh6Sˆ1eh3Tˆe3hSˆ1eh3Tˆeh3Sˆ1e3hTˆeh6Sˆ1 (21)
1
1 1
=exp hHˆ +h3( [[Sˆ ,Tˆ],Tˆ]+ [[Sˆ ,Tˆ],Sˆ ]) , (22)
1 1 1 1
108 216
(cid:16) (cid:17)
and
VH (h)=eh4Sˆ2e2hTˆe2hSˆ2eh2Tˆe4hSˆ2 (23)
2
1 1
=exp hHˆ +h3( [[Sˆ ,Tˆ],Tˆ]+ [[Sˆ ,Tˆ],Sˆ ]) (24)
2 2 2 2
48 96
(cid:16) (cid:17)
Hence we see that the only difference between the individual integrators and our
generalised integrator is the cross term [[Sˆ ,Tˆ],Sˆ ]. The algorithm is identical to a
1 2
standardnestedleapfrog in thecase whereN|N .
i i−1
CONCLUSIONS
We have introduced a novel integration scheme that generalises the nested leapfrog
scheme by Sexton and Weingarten[6]. The new scheme has the advantage that each
terminthefictitiousHamiltonianmaybeassignedatimestepthatonlyneedbeanexact
multiple of the finest time step. This is an improvement on the nested leapfrog, where
each scalemustbean exactmultipleofthenextsmallestscalein thehierarchy.
Each term in the fictitious Hamiltonian will have a corresponding “force term” as-
sociated with it. The typical size of this force term leads one to choose an appropriate
integration time scale for that term. The large the force term, the smaller the time scale
that is required to keep the finite step-size errors under control. With the novel integra-
tionschemeintroducedhereonehasmorefreedomtochooseatime-scalethatisappro-
priate to the force associated with a given term - in the nested leapfrog one may have
been forced to choose a time scale that was smaller than needed due to the restrictions
imposedby thesizeofthenextsmallestscalein thehierarchy.
InaHMCsimulationwithtwodegenerateflavourswherethefermionactionhasbeen
split into two (or more) pieces, we already have at least three different time scales: one
forthegaugeaction,and two(ormore)forthefermionaction.Withthedevelopmentof
multipletimescalealgorithmssuchaspolynomialfilteredHMC[4]thatareapplicableto
bothtwoflavourandsingleflavoursimulations,thenumberoftermsintheHamiltonian
cangrowquitelarge.Itisinthesecasesthatthenovelschemeproposedherewillbecome
particularlyuseful.
Acknowledgments
TheauthorwouldliketothankM.Peardonfordiscussionsrelatingtotheerroranalysis.
REFERENCES
1. S.Duane,A.D.Kennedy,B.J.Pendleton,andD.Roweth,Phys.Lett.B195,216–222(1987).
2. M.A.Clark,andA.D.Kennedy,Phys.Rev.Lett.98,051601(2007),hep-lat/0608015.
3. M.Hasenbusch,Phys.Lett.B519,177–182(2001),hep-lat/0107019.
4. W.Kamleh,andM.J.Peardon,PoSLAT2005,106(2006).
5. M.Lüscher,Comput.Phys.Commun.165,199–220(2005),hep-lat/0409106.
6. J.C.Sexton,andD.H.Weingarten,Nucl.Phys.B380,665–678(1992).
