Lecture VII
An Introduction to Markov Chain Monte Carlo
Lecture VII
Lecture VII
Abstract
The Hybrid Monte Carlo Method: Hamiltonian Dynamics, Liouville's Theorem,
Leapfrog Discretization. The NonReversible Directed Metropolis.
Adding Momentum Variables
Recall that to sample from,
f(x) = 
1 Z

exp{ (logf(x)) } 

we use the Metropolis algorithm with energy function as logf(x). As it
will be seen later, it is useful to think of the x Î R^{n} only as the
position variables q = (q_{1},¼,q_{n}) so that the loglikelihood
plays the role of a potential energy. By introducing an independent extra set
of momentum variables p = (p_{1},¼,p_{n}) with i.i.d. standard
Gaussian distributions we can add a kinetic term to the potential to produce
a full Hamiltonian (energy) function on a ficticious phasespace,
H(q,p) = logf(q) + 
1 2


å
i

p_{i}^{2} 

The Cannonical distribution associated to this Hamiltonian is,
p(q,p) = 
1 Z_{H}

exp{ H(q,p) } 

and we can sample from it by using a combination of Gibbs and Metropolis and
finally discard the momentum variables to show only the position
variables that by construction are now comming from the correct marginal
distribution f(q).
The extra burden of enlarging the dimensions of the sampling space to then
through them away at the end, looks like a complete waste, but there is
ofcourse a benefit associated to this procedure. What is gained by adding the
momentum variables is a way to efficiently explore large regions of
phasespace by simmulating the Hamiltonian dynamics in ficticious time. As it
is shown in the next sections, by following the dynamical path in
phasespace, we can propose cadidate moves that are far away from the current
state but that still have a substantial chance of being accepted.
The validity of the Hybrid Monte Carlo Method rests on three standard
properties of Hamiltonian dynamics:
 Time reversibility, i.e. invariance under
t® t,p® p.
 Conservation of Energy, i.e. H(q,p) is the same for all times.
 Liouville's theorem, i.e. conservation of phasespace volumes.
We look at each of these properties in the next sections.
Hamiltonian Dynamics in Ficticious Time
Let us introduce a ficticious time t along which we assume the system
is evolving according to Hamilton's equations derived in
lecture I. In component form they are,



 

 
¶H ¶q_{i}

= 
¶logf(q) ¶q_{i}



 

where for the rightmost handside we have replaced the specific Hamiltonian
defined in the previous section. It is interesting to notice that in this
ficticious dynamics the field of forces is supplied by the score function
(i.e. the derivatives of the loglikelihood).
Time Reversivility
The above system of equations is invariant under the transformations,
To see it, first notice that the Hamiltonian is the same, i.e.,



logf(q¢) + 
1 2


å
i

(p¢_{i})^{2} 
 

logf(q) + 
1 2


å
i

p_{i}^{2} 
 


 

and the equations are also the same,
This is usually stated by saying that the equations of Classical Mechanics
are blind to the direction of time. The equations can not tell if the
movie of the universe is being run forwards or backwards! I don't know which
universe classical mechanics is trying to describe... certainly not mine!
This is the tip of the iceberg of the so called problem of time. Quantum
theory suffers the same ache. The appearance of the arrow of time has been
related to the second law of thermodynamics (MaxEnt^{TM}), the
expansion of the universe, the Big!, and even conciousness!. I have humbly
(yeh, right..) proposed a connection to
uncertainty in the
measurement of location.
Conservation of Energy
With the help of Hamilton's equations we can compute the rate of change with
respect to time (i.e. the time derivative) of any quantity F(q,p) defined
as a function on phasespace. By the chain rule,

d dt

F(q,p) = 
¶F ¶q


. q

+ 
¶F ¶p


. p



were we have suppressed the summation over the coordinates by denoting the
innerproduct of vectors with a plain product. Replacing Hamilton's equations
we get,

d dt

F(q,p) = 
¶F ¶q


¶H ¶p

 
¶F ¶p


¶H ¶q

º { H,F} 

the last formula is known as the Poisson bracket (it is a Lie
product). The derivative wrt time of any function F is then obtained by
computing the Poisson bracket of F with the Hamiltonian for the system
H. The Hamiltonian tells how things change with time. When F = H we
get,
and therefore things may change but energy doesn't. H i.e. total energy
is a constant of time.
Liouville's Theorem
Take a piece of phasespace of some volume v. Now think of each of the
points inside this volume as different systems with the same Hamiltonian but
with different initial conditions. Let time go and the systems evolve
according to Hamilton's equations. After some time the points initially
inside the volume v would be spread all over. However, the equations of
movement assure that the volume of phasespace covered by these points, is
the same at all times. That is Liouville's theorem and that is what we now
prove.
Take an arbitrary time t_{0} and an arbitrary region of phasespace,
D_{0}. Without loss of generality assume
this time to be zero. Let D(t) be the region of phasespace at time
t occupied by the points that at time zero were in D(0) = D_{0}.
As in the picture.
Figure 1: Illustration of Liouville's Theorem
The volume of D(t) is,
v(t) = 
ó õ

D(t)

dq¢dp¢ = 
ó õ

D(0)


det
 (I+tM) dq dp + O(t) 

where we have used the change of variables,
with the little o's holding as t® 0. The matrix M is the
part of the Jacobian with the derivatives of [q\dot] and [p\dot]. We
have,

det
 (I+tM) = 1 + ttr M + o(t) 

where tr is the trace and we have,
tr M = 
¶q

+ 
¶p

= 
¶ ¶q


æ ç
è


¶H ¶p

ö ÷
ø

 
¶ ¶p


æ ç
è


¶H ¶q

ö ÷
ø

= 0 

where we have replaced Hamilton's equations and assumed continuity for the
second order partial derivatives of H in order to have the equality of
the mixed partials. Thus, we can write,
v(t) = 
ó õ

D(0)

(1+0) dq dp + o(t) = v(0) + o(t) 

from where we obtain,
and since t_{0} = 0 was arbitrary we have shown that the derivative is
zero for all times. Hence, the volume is a constant of time.
Leapfrog Discretization
To simmulate the Hamiltonian dynamics we need to discretize the equations of
motion and in general, unless we are carefull, the error introduced by the
discretization destroys time reversibility and the preservation of volumes
which are both needed for metropolis to work without changes. The so called
leapfrog discretization has the desired properties. Leapfrog updates the
coordinates q_{i}, p_{i} in three steps. First it takes a little half
step for the momentum,

^ p

i

(t+e/2) = 
^ p

i

(t) + 
e 2


¶logf ¶q_{i}

( 
^ q

(t)) 

then it takes a full step for the position,

^ q

i

(t+e) = 
^ q

i

(t) +e 
^ p

i

(t+e/2) 

and finally the other half step for the momentum,

^ p

i

(t+e) = 
^ p

i

(t+e/2) + 
e 2


¶logf ¶q_{i}

( 
^ q

(t+e)) 

leapLeapityleap that's why it's called leapfrog. At the end of the three
steps we obtain an approximation to the values of position and momentum at
time t+e from their corresponding values at time t. As
it can be readily check by simple inspection, the leapfrog discretization has
the following necessary property:
 it allmost preserves H, in fact to order O(e^{2}).
 it preserves volumes since the above are just shear transformations.
 it is time reversible.
The Hybrid Monte Carlo Method
The Hybrid Monte Carlo Method is plain vanilla metropolis in phasespace with
the following proposal distribution:
 Choose random direction of time, i.e., with probability 1/2
simulate the dynamics forward or backward in time.
 Given a current state (q,p) = (q(0),p(0)) of energy H. Perform
a random number L of leapfrog steps with (possibly random) small stepsize
e in the random direction choosen above. This produces a candidate
state (q¢,p¢) with energy H¢.
the candidate state is then accepted with the usual metropolis probability of
acceptance,
The probability of acceptance can be made arbitrarily close to one by
decreasing the stepsize of the discretization. Recall, that the value
of the energy reamins constant along the trajectory of the system in
phasespace. The discretization beeing an inexact simulation of the dynamics
of the system, introduces a small error so that H¢H is not equal to zero
but not large.
Hybrid Monte Carlo Works
The validity of the algorithm follows from the symmetry of the proposal
distribution introduced above. To see that the proposal distribution is
in fact symmetric consider the probability of proposing a small region D¢
given that the current point is known to be in a small region D of volume
dv in phasespace. Since, the simulated Hamiltonian dynamics with
leapfrog discretization is invariant under time reversals and it preserves
volumes it follows that D¢ has volume dv as well and the chance of
proposing a state in D starting from D¢ is the same as the chance of
proposing D¢ from D.
The benefit of following the trajectory of the system in phasespace is the
fact that we can explore quickly regions that are far away from the current
state eliminating the random walk aspect of the chain improving mixing and
producing more accurate estimates.
File translated from T_{E}X by T_{T}H, version 2.32.
On 5 Jul 1999, 22:59.