- Email: [email protected]
Higher order hybrid Monte Carlo
Nuclear Physics B (Proc. Suppl.) 9 (1989) 457-462 North-Holland, Amsterdam
457
HIGHER ORDER HYBRID MONTE CARLO A. D. KENNEDY Supercomputer Computations Research Institute, Florida State University, Tallahassee, FL 32306
The Hybrid Monte Carlo algorithm requires a mapping on phase space which is area-preserving, reversible, and at least approximately conserves energy. We are thus naturally lead to investigate schemes for the numerical integration of Hamilton's equations. First we generalize the "leapfrog" scheme to higher order in the time step (~r, and show why we cannot find an initial "half-step" which satisfies the constraints required for Hybrid Monte Carlo. We then introduce a new symmetric integration scheme which can be used to any order in 6~', and we prove that it too is exactly reversible and area-preserving order-by-order in 67". Finally, we show that the symmetric integration scheme conserves energy exactly for the simple harmonic oscillator, and we give numerical results for the anharmonic oscillator.
1. HYBRID MONTE CARLO
and the system may be moved through phase space by some
The Hybrid Monte Carlo algorithm 1,2 generates a
mapping f : (¢,7r) ~-~ (¢',7r') followed by an test which
Markov process which converges to some Bohzmann distri-
accepts the new state (¢,~r) with probability min(l~ e-$H),
bution, such as those needed for field theory or statistical
where 6H ~ H(¢l,Td) - H(¢,~-), and keeps the old state otherwise.
mechanics Monte Carlo computations. This distribution has as its partition function Z = fide] -S(¢)
For this to be a valid procedure we require that the (I)
where [d@l is a suitable functional measure (e.g., Wiener measurefor a scalar field, a product of Haar measuresfor a
lattice gauge theory). The Hybrid Monte Carlo algorithm has been shown to have an asymptotic volume dependence of
V 5/4 when implemented using the lowest-order leapfrog
integration scheme. 3,4,5 Our objective is to find higherorder integration schemes which get closer to a linear vol-
mapping f be reversibJe, f : (~',-or') ~-~ (~,-Tr), to ensure that detailed balance is satisfied; and area-preserving, c9(¢',~rt)/c~(¢,~r) = 1, to guarantee preservation of the phase space measure. Furthermore, it is clearly desirable to choose f so as to keep 6H as small as possible in order to get a large acceptance rate. As
(e-6H) -- 1 it follows
that 6H = ~6H 2 + O(6H 3) >_ 0 (the last inequality holding for small enough 6~'), and thus keeping H constant solving the classical equations of motion - - is the best we can hope to do.
ume dependence. The Hybrid Monte Carlo algorithm introduces a set of momenta conjugate to the fields ¢, and generates the
2. CLASSICAL MECHANICS We now turn to the problem of formulating an exac~/y reversible and
distribution whose partition function is
exac~ly area-preserving finite step size ap-
proximation to classical mechanics. For simplicity we shall
ZHMc = /[d¢][d~] e-H(¢'~),
(2)
restrict our considerations to a simple unconstrained system described by a dynamical variable q and its conjugate
where the Hamiltonian H(¢,~) --- ~,~2 + S(¢). The system is updated by alternating two type of steps in phase
momentum p. The generalization to systems with many
space: the momenta may be "refreshed" by selecting them at random from a Gaussian distribution P(Tr) oc e-~r2/2;
strained systems, such as gauge theories, is discussed in Ref. [0].
0920-5632/89/$03.50 © Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division)
degrees of freedom is trivial, and the generalization to con-
A.D. Kennedy~Higher order hybrid Monte Carlo
458
For our simple system the classical trajectories are solutions of Hamilton's equations,
'~-
OH
OH
:~- c')q
o~'
tions tell us how A changes along a classical trajectory. If
where by S" we mean S"(~(t)). It is to be understood that
th0se equations de ne the point ({(t +
--
(4)
0
OH
P(q'P)- Op Oq
(q(t +
+
/ (with errors of
To era-
phasize the difference the points generated by the discrete
DA = fi = {A, H} = P(q,p)A, OH
in
This new point is only approximately equal
~/dt then we fincl
where
+
%
OA OB Op Oq
to ]9
(11)
terms of the previous points (c~(t - 5r),/3(t - (~r)) and
we introduce the Poisson bracket
and write
}
q(t +6r)-qCt -St) =25r/~(t )- ½5rsS"i6Ct)
(3)
For any function A(q,p) defined on phase space these equa-
OA cgB { A , B } - Oq Op
the leapfrog equations to O(Sr 5) are
(5)
(q,p).
The leapfrog equations have some remarkable proper-
0
Oq Op"
numerical scheme are denoted by (~,/3) while the points on the true classical trajectory are denoted by
(6)
ties. They are manifestly
exact]y reversible order-by-order
in 57-; this follows from the explicit time reversal symmeEquation (5) tells us how the time derivative of quantity
try in the definition of equation (9). This symmetry also
along a classical trajectory may be expressed solely in terms
ensures that only odd powers of 5~- occur in the leapfrog
of the phase space derivatives P(q,p)A. We shall use this
equations.
formalism to construct discrete approximations to Hamilton's equations.
The equations are also
exactly
area-preserving order-
by-order in 5r: this version of Liouville's theorem states that the Jacobian for the change of variable from ( ~ ( t - 5~')
3. LEAPFROG INTEGRATION
to (~(t ÷ 5r) is unity.
3.1. Leapfrog steps Let us denote a classical path through phase space by the function Q(t) ---- (q(t),p(t)), then using Taylor's theorem (or equivalently noting that D is the generator of
There are two simple ways of seeing this. If we explicitly compute the Jacobian
,~ ----det
time translations)
Q(t + 5"r) = e±~rDQ(t),
(7)
0[O(t+~-),O(t)]
(12)
we see that it is the determinant of a triangular matrix for each separate ~ and ~ leap. This is because the left
so we find
side of equation (9) contains all the (~(t-I- 5v) dependence,
Q(t + 5v) - QCt - 5v) -- 2 sinhCSrD)Q(t)
-- 2sinhCSrPQ)QCl).
(8)
whereas the right side (truncated to any order in E~-) is
(9)
only a function of the "spectator" variables Q(t) which are unchanged by a leap.
restatement of Hamilton's
The second way of getting this result is to consider the
equations, but in this form we may expand the phase space
transformation f : Q(t - ~-) F-~ (~(t + 5~-) followed by a
function as a power series in the small time-step 5v, and truncate it to any desired order. For example, choosing the
momentum reversal "P : (~,t3) ~-~ (~,-~), followed by f and "P again. The composite operation "P o f o "P o f ---- 1 by reversibility, and thus the Jacobians satisfy
This equation is just an
exact
Hamiltonian I 2 +S(q), H(p,q)= ~p
,7(7 ~ o f o'P o jr) =if(1) = 1 (10)
= j(;)~j(f)~
_- ~ ( ~ ) ~
(13)
A.D. Kennedy~Higher order hybrid Monte Carlo
459
is a function only of (~(t) for
ing with a final ~ half step as specified by equation (15).
both the forward and backward trajectories. By continuity
Correspondingly, another trajectory is defined by starting
in 6r we thus have i f ( f ) = 1 as required.
and ending with ~ half steps: reversibility requires that the
since ff('P) = ±1 and i f ( f )
It is worth understanding what this proof is telling us. It says that the explicit reversibility of the leapfrog scheme
initial and final half steps be of the same type• We are now in a position to understand what goes
combined with the fact that the right side of the leapfrog
wrong when we try to generalize equations (14) and (15)
equations 9 depends only on a symmetric phase space point
to higher orders in the step size 6~-. Beyond O(6r) we need
- - one which is the same whether we traverse the trajectory
to know both ~(t) and ~(t) in order to leap either ~ or
forwards or backwards - - is enough to guarantee the area-
from time t - 6~- to t + 6r. In order for this higher order
preservation property we need.
leapfrog scheme to be reversible, the endpoint of the final
3.2. Initial and final half steps
leap must coincide with the endpoint of the final half step,
The leapfrog scheme as we have so far presented it is fundamentally incomplete, as we know only the starting point Q(0) of our putative classical trajectory, and the leapfrog equations also require us to specify the values of (~(6r). We therefore must also give a recipe for initial (and
and it is easily verified that this is not true for any obvious generalizations of the half step equations (14) and (15). Is it perhaps possible to invent a cleverer half step which will circumvent this problem? The initial half step must be uniquely determined by the value of the initial point (~(0).
final) "half steps." In the lowest order leapfrog scheme, where we only include the terms linear in ~r in equation (9), it is easy to construct a suitable half step. All we need do is to use the leading term from the Taylor expansion ? about the initial and final points:
By reversibility, the final half step depends
only on the final point (~(r0). The final leap depends only on the value of the potential energy S and a finite number of its derivatives at (~(r 0 - 6r). If we modify the potential in a small enough neighbourhood of the final point, then the last leap is unaffected, so the final half step must be unchanged too. Thus we have shown that a half step cannot depend on the potential in a neighbourhood of its start-
= ,0(o) +
(14)
a Taylor expansion, however we must require that the ap-
= ,0(,o)
: 00o)
ing point. [In deriving the leapfrog equations we implicitly assumed that the potential was smooth enough to admit
Po0( o)
(is)
proximate integration scheme is exactly reversible even if this is not true: it will just give a poor approximation to
The initial and final half steps map into each other under
the solutions of Hamilton's equations in such a case.] Any
time reversal, so the reversibility of the entire leapfrog tra-
such half step scheme would certainly be very peculiar.
jectory is maintained, and it is easily verified that Liouville's theorem is still satisfied.
4. SYMMETRIC INTEGRATION SCHEME
What is a little more subtle is that equations (14) and (15) specify more than is necessary. If the Hamiltonian H(q,p) may be written as the sum of a kinetic
The results of the previous section have shown that although a leapfrog integration scheme correct to any desired order in ~r exists and has all the desired properties, there
energy term which depends only upon p and a potential
appears to be no generalization of the initial and final half
energy term which depends only upon q then, as we are only working to lowest non-trivial order in 5~-, a ~-Ieap only
steps. We therefore turn instead to a new scheme which will allow us to carry out area-preserving reversible approx-
requires the value of c~ at its midpoint, and a c~-Ieap only
imations to Hamilton's equations accurate to any specified
requires the value of 2. A complete trajectory is specified
order in the step size parameter.
by using equation (14) to give only an initial ¢~ half step, following this with a sequence of p-leaps and j-leaps, end-
end of section 3.1. In order to guarantee exact area preser-
The clue to doing this is the observation made at the
460
A.D. Kennedy~Higher order hybrid Monte Carlo
(22)
vation and reversibility it suffices to use a scheme in which Q(t + 6T) - Q(t + ~/") is a function only of quantities
(23)
--_1 [p~(,+,.)Q(~+5.)+p~(oQ(t)]
evaluated at a symmetric point: that is one which does
(241
not depend on the direction in which the trajectory is tra-
={[eQ-'"~.-Q-"IP~+Q+;
versed. In the leapfrog this point is, at least approximately, the midpoint of the classical trajectory; this need not be so, however. In the "symmetric integration scheme" introduced here we shall instead choose the symmetric point to be a symmetric linear combination of the starting point and ending point of an integration step.
The fact that
this point is not on the classical trajectory leads to a few differences in the analysis, but there are no fundamental differences.
(25)
where we have introduced the operator 0 - ( O/Oq, a/ap). Using this result we can write down the general formula for a symmetric integration step:
o
O)
+51 cosh(Q-.
O+
(26)
O)siRh(E/"PQ+)Q+.
This is still an implicit equation for Q_, as Q_ occurs on
.Just as before we start by Taylor expansion, but now
the right hand side as well as on the left; but if we notice
we Taylor expand the final point Q(~ + ~/") about the initial
that the right hand side of equation (26) is always 0(5/")
point Q(t) and vice versa:
then we see that iteratively substituting the right hand side
Q(t +
Er) = e~DQ(t)
QCt) = e-6~D QCt +
}
(16)
6r).
into itself will generate an explicit equation for Q_ in terms of Q+ to any desired order in E/". For example, choosing the Hamiltonian of equation (10)
This allows us to express their difference as
we obtain the following equations truncated at O(5T 7)
Q(t + 6/.)-Q(t) = e&~'DQ(t) - e-~rD Q(t + 6/.)(17) = cosh(grD)[O(t) - O(t + E/.)]
+ sinh(6vD)[O(t)+ O(t + 6r)]. (18)
~_= ½~+ + ~ + s " ( ~ + ) +~67" ' '[~+s'"'+~+s's'"+4~+s":] ^ -1 I t 8 ^~ p_--~,~rS-~,~"r [25't SI I +5'II1 p+]
Let us now define our symmetric variables explicitly,
Q±
1
[QCt+ ~/.) ± QCt)],
(27)
a--~5,'r [p+S +4p+S S +2o~+s,,s,,, +8s,~s,,,+16s, s,,~]. 1
~
^4
IIIII
^~
I
IIII
(19) Unlike the leapfrog equations these equations do not
so [The second line below follows from the first by adding Q_ to both sides of the equation and dividing by two.] Q_ = - cosh(g/.D)Q- + sinh(grD)Q+ =111-cosh(g/.Dl]O-
2
(20)
+ 2sinh(5/.DlO+.(21)
We now want to apply Hamilton's equations in the same way as we did in equation (9). We must be a little careful, as the equations of motion (5) can only be applied to dynamical variables and Q+, being the midpoint of two
tell us the new point in phase space directly. They just give a relationship between the variables Q_ and Q+ which must be satisfied. For small enough 5r the equations can be solved by iteration: starting with ~ + = (~(t) one may find Q_ to the next order in 6T, then set Q(t + 6v) = Q(t) + 2Q_, and then Q+ = ½[(~(t 4- 6r) 4- Q(t)]. This may be iterated until the solution converges to a new exact value for Q(t + 61") - - exact here meaning to machine precision, presumably.
dynamical variable, is not a dynamical variable itself. To
It may be helpful to use the solution of a lower or-
circumvent this we shall take care to apply the equations
der symmetric integration as the starting point for a higher
of motion only to the true dynamical variables Q(t) and
order one, or indeed one could even try using a leapfrog
Q(t + 6"r), and then to extract the behaviour of
method to provide a good initial approximation for the it-
Taylor expansion in phase space:
Q+ by
erative process. In general there will exist multiple solutions
A.D. Kennedy~Higher order hybrid Monte Carlo
[~-e - lJ- - [
to the truncated versions of equation (26), but we are al-
-~
461
6~+ 1 L-6'~+-:~'~ J
ways interested in the one which has (~_ = 0 as ~v --~ 0.
(321
It is easy to verify that the symmetric scheme equa-
F36~++;'(54~+'~-+9,~-)+12;'~+'~]-
tions, truncated to any order in Ev but solved exactly there~r 5
after, are both reversible and area preserving. The proof
.
+ 68T4°I
.
.
.
-36q+-A(60q~+54p2+q+)
L -;~(21~++24~+~+)-2~+
is exactly as given for the leapfrog equations around equation (13).
"
As we cannot find a simple closed form solution to these equations we shall solve them numerically by the iteration scheme described before. The results for ,~ = 15 are shown
5. SIMPLE HARMONIC OSCILLATOR Let us now turn to the simplest example we can find, the one-dimensional simple harmonic oscillator, described by the Hamiltonian H = lp2 + ½q2. The sixth-order symmetric equations (26) give
in Figure 1, in which we plot the logarithm of the magnitude of the energy violation 16HI against the logarithm of the step size 6v. The points to note are that the errors scale with 6T just as expected, and the errors for an en-
1/67" steps,
tire trajectory of length r 0 = 1, which involves
are one power of 6v bigger than those for a single step.
=
+~+
240.)
-,i+ '
The most interesting phenomenon is that for a given step size the errors of the lowest order symmetric scheme are
Defining c~ to be expression in parentheses, we can
about 50 times smaller than those of the leapfrog scheme
solve these equations to obtain _
h=15
:o r 1 Ll~(Oj~ (~+--~:) L-4(o)J'
t
(2!))
.......i:/:iii:iii -1
and letting 0 _= 2tan -1 c~
c
,.'""
-20
,."" ."" ,'"
We see that for the simple harmonic oscillator the symmetric integration scheme follows the exact trajectory
-30
~L_L~_,
through phase space, the only error being the phase ~(6~-). In particular, the Hamiltonian is exactly conserved by the
-8
-7
-6
-5
-4 In 6~
3
2
-1
symmetric equations of motion truncated to any order
in 6"r.
FIGURE
6. ANHARMONIC OSCILLATOR In an attempt to find a model in which 6H behaves non-trivially in the symmetric integration scheme we consider the anharmonic oscillator described by the Hamiltonian
p2
q2
)lq4
H=y+~-+ 4!
The symmetric equations of motion truncated at are
(31) O(gr 7)
1
The magnitude of energy violation ~6H I is shown as a function of the step size ov on a logarithmic plot. The solid symbols are the results after a single step of size 6v, whereas the open symbols are after a sequence of steps of total duration v0 = 1. The top curves are for the lowest order leapfrog algorithm with a /3 half step (the c~half step results look similar). The lower curves are for the symmetric scheme truncated at O(6~-3), O(6r5), and O(6~-7) respectively. The dashed lines are just the curves 6H = 6r n for n = 1, 2 , . . . , 8. The few strange points at the bottom of the graph are just the effects of finite machine precision.
~
A.D. Kennedy~Higher order hybrid Monte Carlo
462
which is nominally correct to the same order. Whether this
even comparing the lowest order symmetric scheme results
desirable behaviour applies to a more complicated system
with those of leapfrog. It remains to be seen how effec-
remains to be seen.
tive the symmetric scheme will be for systems of greater
7. CONCLUSIONS
physical interest.
We have studied higher order generalizations of the leapfrog algorithm, and we have discovered that we cannot find an initial "half step" which will lead to a reversible and area preserving integration scheme beyond lowest order. We have introduced a new "symmetric integration scheme" which is reversible and area preserving to all orders in the step size 5r. This scheme is described implicitly by equation (20), which must be iterated and truncated to some given order in 5v to obtain the equations determining the approximate classical trajectory. This process may easily be done analytically; however, in general the resulting
REFERENCES i. Simon Duane, A. D. Kennedy, Brian J. Pendleton, and Duncan Roweth, Hybrid Monte Carlo, Phys. Left., 195(2), 1987. 2. A. D. Kennedy, Hybrid Monte Carlo, In Field Theory on the Lattice, A. Billoire, R. Lacaze, A. Morel, O. Napoly, and J. Zinn-Justin, editors, pages 576-579, 1988. 1987 International Conference on Field Theory on a Lattice, Seillac, France. 3. Raian Gupta, Gregory W. Kilcup, and Stephen R. Sharpe, Tuning the hybrid Monte Carlo algorithm, Phys. Rev., D38(4):1278-1287, August 1988.
relate the variables Q± rather than directly giving Q(t+Sr) in terms of (~(t).
4. Michael Creutz, Algorithmic alternatives, In Field Theory on the Lattice, A. Billoire, R. Lacaze, A. Morel, O. Napoly, and J. Zinn-Justin, editors, pages 547-556, 1988. 1987 International Conference on Field Theory on a Lattice, Seillac, France.
Our studies in toy systems show that the cost involved in having to solve the equations of motion to machine pre-
5. Michael Creutz, Global Monte Carlo algorithms for many fermion systems, Phys. Rev., D38(4):12281238, August 1988.
cision for every step of the trajectory may be compensated
6. A. D. Kennedy and Pietro Rossi, Classical mechanics on group manifolds, 1988. (In preparation.)
equations of motion have to be solved numerically, as they
for by the much smaller energy violations that are produced,
457
HIGHER ORDER HYBRID MONTE CARLO A. D. KENNEDY Supercomputer Computations Research Institute, Florida State University, Tallahassee, FL 32306
The Hybrid Monte Carlo algorithm requires a mapping on phase space which is area-preserving, reversible, and at least approximately conserves energy. We are thus naturally lead to investigate schemes for the numerical integration of Hamilton's equations. First we generalize the "leapfrog" scheme to higher order in the time step (~r, and show why we cannot find an initial "half-step" which satisfies the constraints required for Hybrid Monte Carlo. We then introduce a new symmetric integration scheme which can be used to any order in 6~', and we prove that it too is exactly reversible and area-preserving order-by-order in 67". Finally, we show that the symmetric integration scheme conserves energy exactly for the simple harmonic oscillator, and we give numerical results for the anharmonic oscillator.
1. HYBRID MONTE CARLO
and the system may be moved through phase space by some
The Hybrid Monte Carlo algorithm 1,2 generates a
mapping f : (¢,7r) ~-~ (¢',7r') followed by an test which
Markov process which converges to some Bohzmann distri-
accepts the new state (¢,~r) with probability min(l~ e-$H),
bution, such as those needed for field theory or statistical
where 6H ~ H(¢l,Td) - H(¢,~-), and keeps the old state otherwise.
mechanics Monte Carlo computations. This distribution has as its partition function Z = fide] -S(¢)
For this to be a valid procedure we require that the (I)
where [d@l is a suitable functional measure (e.g., Wiener measurefor a scalar field, a product of Haar measuresfor a
lattice gauge theory). The Hybrid Monte Carlo algorithm has been shown to have an asymptotic volume dependence of
V 5/4 when implemented using the lowest-order leapfrog
integration scheme. 3,4,5 Our objective is to find higherorder integration schemes which get closer to a linear vol-
mapping f be reversibJe, f : (~',-or') ~-~ (~,-Tr), to ensure that detailed balance is satisfied; and area-preserving, c9(¢',~rt)/c~(¢,~r) = 1, to guarantee preservation of the phase space measure. Furthermore, it is clearly desirable to choose f so as to keep 6H as small as possible in order to get a large acceptance rate. As
(e-6H) -- 1 it follows
that 6H = ~6H 2 + O(6H 3) >_ 0 (the last inequality holding for small enough 6~'), and thus keeping H constant solving the classical equations of motion - - is the best we can hope to do.
ume dependence. The Hybrid Monte Carlo algorithm introduces a set of momenta conjugate to the fields ¢, and generates the
2. CLASSICAL MECHANICS We now turn to the problem of formulating an exac~/y reversible and
distribution whose partition function is
exac~ly area-preserving finite step size ap-
proximation to classical mechanics. For simplicity we shall
ZHMc = /[d¢][d~] e-H(¢'~),
(2)
restrict our considerations to a simple unconstrained system described by a dynamical variable q and its conjugate
where the Hamiltonian H(¢,~) --- ~,~2 + S(¢). The system is updated by alternating two type of steps in phase
momentum p. The generalization to systems with many
space: the momenta may be "refreshed" by selecting them at random from a Gaussian distribution P(Tr) oc e-~r2/2;
strained systems, such as gauge theories, is discussed in Ref. [0].
0920-5632/89/$03.50 © Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division)
degrees of freedom is trivial, and the generalization to con-
A.D. Kennedy~Higher order hybrid Monte Carlo
458
For our simple system the classical trajectories are solutions of Hamilton's equations,
'~-
OH
OH
:~- c')q
o~'
tions tell us how A changes along a classical trajectory. If
where by S" we mean S"(~(t)). It is to be understood that
th0se equations de ne the point ({(t +
--
(4)
0
OH
P(q'P)- Op Oq
(q(t +
+
/ (with errors of
To era-
phasize the difference the points generated by the discrete
DA = fi = {A, H} = P(q,p)A, OH
in
This new point is only approximately equal
~/dt then we fincl
where
+
%
OA OB Op Oq
to ]9
(11)
terms of the previous points (c~(t - 5r),/3(t - (~r)) and
we introduce the Poisson bracket
and write
}
q(t +6r)-qCt -St) =25r/~(t )- ½5rsS"i6Ct)
(3)
For any function A(q,p) defined on phase space these equa-
OA cgB { A , B } - Oq Op
the leapfrog equations to O(Sr 5) are
(5)
(q,p).
The leapfrog equations have some remarkable proper-
0
Oq Op"
numerical scheme are denoted by (~,/3) while the points on the true classical trajectory are denoted by
(6)
ties. They are manifestly
exact]y reversible order-by-order
in 57-; this follows from the explicit time reversal symmeEquation (5) tells us how the time derivative of quantity
try in the definition of equation (9). This symmetry also
along a classical trajectory may be expressed solely in terms
ensures that only odd powers of 5~- occur in the leapfrog
of the phase space derivatives P(q,p)A. We shall use this
equations.
formalism to construct discrete approximations to Hamilton's equations.
The equations are also
exactly
area-preserving order-
by-order in 5r: this version of Liouville's theorem states that the Jacobian for the change of variable from ( ~ ( t - 5~')
3. LEAPFROG INTEGRATION
to (~(t ÷ 5r) is unity.
3.1. Leapfrog steps Let us denote a classical path through phase space by the function Q(t) ---- (q(t),p(t)), then using Taylor's theorem (or equivalently noting that D is the generator of
There are two simple ways of seeing this. If we explicitly compute the Jacobian
,~ ----det
time translations)
Q(t + 5"r) = e±~rDQ(t),
(7)
0[O(t+~-),O(t)]
(12)
we see that it is the determinant of a triangular matrix for each separate ~ and ~ leap. This is because the left
so we find
side of equation (9) contains all the (~(t-I- 5v) dependence,
Q(t + 5v) - QCt - 5v) -- 2 sinhCSrD)Q(t)
-- 2sinhCSrPQ)QCl).
(8)
whereas the right side (truncated to any order in E~-) is
(9)
only a function of the "spectator" variables Q(t) which are unchanged by a leap.
restatement of Hamilton's
The second way of getting this result is to consider the
equations, but in this form we may expand the phase space
transformation f : Q(t - ~-) F-~ (~(t + 5~-) followed by a
function as a power series in the small time-step 5v, and truncate it to any desired order. For example, choosing the
momentum reversal "P : (~,t3) ~-~ (~,-~), followed by f and "P again. The composite operation "P o f o "P o f ---- 1 by reversibility, and thus the Jacobians satisfy
This equation is just an
exact
Hamiltonian I 2 +S(q), H(p,q)= ~p
,7(7 ~ o f o'P o jr) =if(1) = 1 (10)
= j(;)~j(f)~
_- ~ ( ~ ) ~
(13)
A.D. Kennedy~Higher order hybrid Monte Carlo
459
is a function only of (~(t) for
ing with a final ~ half step as specified by equation (15).
both the forward and backward trajectories. By continuity
Correspondingly, another trajectory is defined by starting
in 6r we thus have i f ( f ) = 1 as required.
and ending with ~ half steps: reversibility requires that the
since ff('P) = ±1 and i f ( f )
It is worth understanding what this proof is telling us. It says that the explicit reversibility of the leapfrog scheme
initial and final half steps be of the same type• We are now in a position to understand what goes
combined with the fact that the right side of the leapfrog
wrong when we try to generalize equations (14) and (15)
equations 9 depends only on a symmetric phase space point
to higher orders in the step size 6~-. Beyond O(6r) we need
- - one which is the same whether we traverse the trajectory
to know both ~(t) and ~(t) in order to leap either ~ or
forwards or backwards - - is enough to guarantee the area-
from time t - 6~- to t + 6r. In order for this higher order
preservation property we need.
leapfrog scheme to be reversible, the endpoint of the final
3.2. Initial and final half steps
leap must coincide with the endpoint of the final half step,
The leapfrog scheme as we have so far presented it is fundamentally incomplete, as we know only the starting point Q(0) of our putative classical trajectory, and the leapfrog equations also require us to specify the values of (~(6r). We therefore must also give a recipe for initial (and
and it is easily verified that this is not true for any obvious generalizations of the half step equations (14) and (15). Is it perhaps possible to invent a cleverer half step which will circumvent this problem? The initial half step must be uniquely determined by the value of the initial point (~(0).
final) "half steps." In the lowest order leapfrog scheme, where we only include the terms linear in ~r in equation (9), it is easy to construct a suitable half step. All we need do is to use the leading term from the Taylor expansion ? about the initial and final points:
By reversibility, the final half step depends
only on the final point (~(r0). The final leap depends only on the value of the potential energy S and a finite number of its derivatives at (~(r 0 - 6r). If we modify the potential in a small enough neighbourhood of the final point, then the last leap is unaffected, so the final half step must be unchanged too. Thus we have shown that a half step cannot depend on the potential in a neighbourhood of its start-
= ,0(o) +
(14)
a Taylor expansion, however we must require that the ap-
= ,0(,o)
: 00o)
ing point. [In deriving the leapfrog equations we implicitly assumed that the potential was smooth enough to admit
Po0( o)
(is)
proximate integration scheme is exactly reversible even if this is not true: it will just give a poor approximation to
The initial and final half steps map into each other under
the solutions of Hamilton's equations in such a case.] Any
time reversal, so the reversibility of the entire leapfrog tra-
such half step scheme would certainly be very peculiar.
jectory is maintained, and it is easily verified that Liouville's theorem is still satisfied.
4. SYMMETRIC INTEGRATION SCHEME
What is a little more subtle is that equations (14) and (15) specify more than is necessary. If the Hamiltonian H(q,p) may be written as the sum of a kinetic
The results of the previous section have shown that although a leapfrog integration scheme correct to any desired order in ~r exists and has all the desired properties, there
energy term which depends only upon p and a potential
appears to be no generalization of the initial and final half
energy term which depends only upon q then, as we are only working to lowest non-trivial order in 5~-, a ~-Ieap only
steps. We therefore turn instead to a new scheme which will allow us to carry out area-preserving reversible approx-
requires the value of c~ at its midpoint, and a c~-Ieap only
imations to Hamilton's equations accurate to any specified
requires the value of 2. A complete trajectory is specified
order in the step size parameter.
by using equation (14) to give only an initial ¢~ half step, following this with a sequence of p-leaps and j-leaps, end-
end of section 3.1. In order to guarantee exact area preser-
The clue to doing this is the observation made at the
460
A.D. Kennedy~Higher order hybrid Monte Carlo
(22)
vation and reversibility it suffices to use a scheme in which Q(t + 6T) - Q(t + ~/") is a function only of quantities
(23)
--_1 [p~(,+,.)Q(~+5.)+p~(oQ(t)]
evaluated at a symmetric point: that is one which does
(241
not depend on the direction in which the trajectory is tra-
={[eQ-'"~.-Q-"IP~+Q+;
versed. In the leapfrog this point is, at least approximately, the midpoint of the classical trajectory; this need not be so, however. In the "symmetric integration scheme" introduced here we shall instead choose the symmetric point to be a symmetric linear combination of the starting point and ending point of an integration step.
The fact that
this point is not on the classical trajectory leads to a few differences in the analysis, but there are no fundamental differences.
(25)
where we have introduced the operator 0 - ( O/Oq, a/ap). Using this result we can write down the general formula for a symmetric integration step:
o
O)
+51 cosh(Q-.
O+
(26)
O)siRh(E/"PQ+)Q+.
This is still an implicit equation for Q_, as Q_ occurs on
.Just as before we start by Taylor expansion, but now
the right hand side as well as on the left; but if we notice
we Taylor expand the final point Q(~ + ~/") about the initial
that the right hand side of equation (26) is always 0(5/")
point Q(t) and vice versa:
then we see that iteratively substituting the right hand side
Q(t +
Er) = e~DQ(t)
QCt) = e-6~D QCt +
}
(16)
6r).
into itself will generate an explicit equation for Q_ in terms of Q+ to any desired order in E/". For example, choosing the Hamiltonian of equation (10)
This allows us to express their difference as
we obtain the following equations truncated at O(5T 7)
Q(t + 6/.)-Q(t) = e&~'DQ(t) - e-~rD Q(t + 6/.)(17) = cosh(grD)[O(t) - O(t + E/.)]
+ sinh(6vD)[O(t)+ O(t + 6r)]. (18)
~_= ½~+ + ~ + s " ( ~ + ) +~67" ' '[~+s'"'+~+s's'"+4~+s":] ^ -1 I t 8 ^~ p_--~,~rS-~,~"r [25't SI I +5'II1 p+]
Let us now define our symmetric variables explicitly,
Q±
1
[QCt+ ~/.) ± QCt)],
(27)
a--~5,'r [p+S +4p+S S +2o~+s,,s,,, +8s,~s,,,+16s, s,,~]. 1
~
^4
IIIII
^~
I
IIII
(19) Unlike the leapfrog equations these equations do not
so [The second line below follows from the first by adding Q_ to both sides of the equation and dividing by two.] Q_ = - cosh(g/.D)Q- + sinh(grD)Q+ =111-cosh(g/.Dl]O-
2
(20)
+ 2sinh(5/.DlO+.(21)
We now want to apply Hamilton's equations in the same way as we did in equation (9). We must be a little careful, as the equations of motion (5) can only be applied to dynamical variables and Q+, being the midpoint of two
tell us the new point in phase space directly. They just give a relationship between the variables Q_ and Q+ which must be satisfied. For small enough 5r the equations can be solved by iteration: starting with ~ + = (~(t) one may find Q_ to the next order in 6T, then set Q(t + 6v) = Q(t) + 2Q_, and then Q+ = ½[(~(t 4- 6r) 4- Q(t)]. This may be iterated until the solution converges to a new exact value for Q(t + 61") - - exact here meaning to machine precision, presumably.
dynamical variable, is not a dynamical variable itself. To
It may be helpful to use the solution of a lower or-
circumvent this we shall take care to apply the equations
der symmetric integration as the starting point for a higher
of motion only to the true dynamical variables Q(t) and
order one, or indeed one could even try using a leapfrog
Q(t + 6"r), and then to extract the behaviour of
method to provide a good initial approximation for the it-
Taylor expansion in phase space:
Q+ by
erative process. In general there will exist multiple solutions
A.D. Kennedy~Higher order hybrid Monte Carlo
[~-e - lJ- - [
to the truncated versions of equation (26), but we are al-
-~
461
6~+ 1 L-6'~+-:~'~ J
ways interested in the one which has (~_ = 0 as ~v --~ 0.
(321
It is easy to verify that the symmetric scheme equa-
F36~++;'(54~+'~-+9,~-)+12;'~+'~]-
tions, truncated to any order in Ev but solved exactly there~r 5
after, are both reversible and area preserving. The proof
.
+ 68T4°I
.
.
.
-36q+-A(60q~+54p2+q+)
L -;~(21~++24~+~+)-2~+
is exactly as given for the leapfrog equations around equation (13).
"
As we cannot find a simple closed form solution to these equations we shall solve them numerically by the iteration scheme described before. The results for ,~ = 15 are shown
5. SIMPLE HARMONIC OSCILLATOR Let us now turn to the simplest example we can find, the one-dimensional simple harmonic oscillator, described by the Hamiltonian H = lp2 + ½q2. The sixth-order symmetric equations (26) give
in Figure 1, in which we plot the logarithm of the magnitude of the energy violation 16HI against the logarithm of the step size 6v. The points to note are that the errors scale with 6T just as expected, and the errors for an en-
1/67" steps,
tire trajectory of length r 0 = 1, which involves
are one power of 6v bigger than those for a single step.
=
+~+
240.)
-,i+ '
The most interesting phenomenon is that for a given step size the errors of the lowest order symmetric scheme are
Defining c~ to be expression in parentheses, we can
about 50 times smaller than those of the leapfrog scheme
solve these equations to obtain _
h=15
:o r 1 Ll~(Oj~ (~+--~:) L-4(o)J'
t
(2!))
.......i:/:iii:iii -1
and letting 0 _= 2tan -1 c~
c
,.'""
-20
,."" ."" ,'"
We see that for the simple harmonic oscillator the symmetric integration scheme follows the exact trajectory
-30
~L_L~_,
through phase space, the only error being the phase ~(6~-). In particular, the Hamiltonian is exactly conserved by the
-8
-7
-6
-5
-4 In 6~
3
2
-1
symmetric equations of motion truncated to any order
in 6"r.
FIGURE
6. ANHARMONIC OSCILLATOR In an attempt to find a model in which 6H behaves non-trivially in the symmetric integration scheme we consider the anharmonic oscillator described by the Hamiltonian
p2
q2
)lq4
H=y+~-+ 4!
The symmetric equations of motion truncated at are
(31) O(gr 7)
1
The magnitude of energy violation ~6H I is shown as a function of the step size ov on a logarithmic plot. The solid symbols are the results after a single step of size 6v, whereas the open symbols are after a sequence of steps of total duration v0 = 1. The top curves are for the lowest order leapfrog algorithm with a /3 half step (the c~half step results look similar). The lower curves are for the symmetric scheme truncated at O(6~-3), O(6r5), and O(6~-7) respectively. The dashed lines are just the curves 6H = 6r n for n = 1, 2 , . . . , 8. The few strange points at the bottom of the graph are just the effects of finite machine precision.
~
A.D. Kennedy~Higher order hybrid Monte Carlo
462
which is nominally correct to the same order. Whether this
even comparing the lowest order symmetric scheme results
desirable behaviour applies to a more complicated system
with those of leapfrog. It remains to be seen how effec-
remains to be seen.
tive the symmetric scheme will be for systems of greater
7. CONCLUSIONS
physical interest.
We have studied higher order generalizations of the leapfrog algorithm, and we have discovered that we cannot find an initial "half step" which will lead to a reversible and area preserving integration scheme beyond lowest order. We have introduced a new "symmetric integration scheme" which is reversible and area preserving to all orders in the step size 5r. This scheme is described implicitly by equation (20), which must be iterated and truncated to some given order in 5v to obtain the equations determining the approximate classical trajectory. This process may easily be done analytically; however, in general the resulting
REFERENCES i. Simon Duane, A. D. Kennedy, Brian J. Pendleton, and Duncan Roweth, Hybrid Monte Carlo, Phys. Left., 195(2), 1987. 2. A. D. Kennedy, Hybrid Monte Carlo, In Field Theory on the Lattice, A. Billoire, R. Lacaze, A. Morel, O. Napoly, and J. Zinn-Justin, editors, pages 576-579, 1988. 1987 International Conference on Field Theory on a Lattice, Seillac, France. 3. Raian Gupta, Gregory W. Kilcup, and Stephen R. Sharpe, Tuning the hybrid Monte Carlo algorithm, Phys. Rev., D38(4):1278-1287, August 1988.
relate the variables Q± rather than directly giving Q(t+Sr) in terms of (~(t).
4. Michael Creutz, Algorithmic alternatives, In Field Theory on the Lattice, A. Billoire, R. Lacaze, A. Morel, O. Napoly, and J. Zinn-Justin, editors, pages 547-556, 1988. 1987 International Conference on Field Theory on a Lattice, Seillac, France.
Our studies in toy systems show that the cost involved in having to solve the equations of motion to machine pre-
5. Michael Creutz, Global Monte Carlo algorithms for many fermion systems, Phys. Rev., D38(4):12281238, August 1988.
cision for every step of the trajectory may be compensated
6. A. D. Kennedy and Pietro Rossi, Classical mechanics on group manifolds, 1988. (In preparation.)
equations of motion have to be solved numerically, as they
for by the much smaller energy violations that are produced,