MCMC Using Hamiltonian Dynamics

1 5 mcmc Using Hamiltonian DynamicsRadford M. IntroductionMarkov chain Monte Carlo ( mcmc ) originated with the classic paper of Metropolis et al.(1953), where it was used to simulate the distribution of states for a system of idealizedmolecules. Not long after, another approach to molecular simulation was introduced (Alderand Wainwright, 1959), in which the motion of the molecules was deterministic, followingNewton s laws of motion, which have an elegant formalization asHamiltonian Dynamics . Forfinding the properties of bulk materials, these approaches are asymptotically equivalent,since even in a deterministic simulation, each local region of the material experienceseffectively random influences from distant regions.

2 Despite the large overlap in their appli-cation areas, the mcmc and molecular Dynamics approaches have continued to coexist inthe following decades (see Frenkel and Smit, 1996).In 1987, a landmark paper by Duane, Kennedy, Pendleton, and Roweth united theMCMC and molecular Dynamics approaches. They called their method hybrid MonteCarlo, which abbreviates to HMC, but the phrase Hamiltonian Monte Carlo, retain-ing the abbreviation, is more specific and descriptive, and I will use it here. Duane et HMC not to molecular simulation, but to lattice field theory simulations of quan-tum chromodynamics.

3 Statistical applications of HMC began with my use of it for neuralnetwork models (Neal, 1996a). I also provided a statistically-oriented tutorial on HMC in areview of mcmc methods (Neal, 1993, Chapter 5). There have been other applicationsof HMC to statistical problems ( Ishwaran, 1999; Schmidt, 2009) and statistically-oriented reviews ( Liu, 2001, Chapter 9), but HMC still seems to be underappreciatedby statisticians, and perhaps also by physicists outside the lattice field theory review begins by describing Hamiltonian Dynamics .

4 Despite terminology that maybe unfamiliar outside physics, the features of Hamiltonian Dynamics that are needed forHMC are elementary. The differential equations of Hamiltonian Dynamics must be dis-cretized for computer implementation. The leapfrog scheme that is typically used isquite this introduction to Hamiltonian Dynamics , I describe how to use it to con-struct an mcmc method. The first step is to define a Hamiltonian function in terms of theprobability distribution we wish to sample from. In addition to the variables we are inter-ested in (the position variables), we must introduce auxiliary momentum variables,which typically have independent Gaussian distributions.

5 The HMC method alternatessimple updates for these momentum variables with Metropolis updates in which a newstate is proposed by computing a trajectory according to Hamiltonian Dynamics , imple-mented with the leapfrog method. A state proposed in this way can be distant from the113114 Handbook of Markov Chain Monte Carlocurrent state but nevertheless have a high probability of acceptance. This bypasses the slowexploration of the state space that occurs when Metropolis updates are done Using a simplerandom-walk proposal distribution.

6 (An alternative way of avoiding random walks is to useshort trajectories but only partially replace the momentum variables between trajectories,so that successive trajectories tend to move in the same direction.)After presenting the basic HMC method, I discuss practical issues of tuning the leapfrogstepsize and number of leapfrog steps, as well as theoretical results on the scaling of HMCwith dimensionality. I then present a number of variations on HMC. The acceptance ratefor HMC can be increased for many problems by looking at windows of states at thebeginning and end of the trajectory.

7 For many statistical problems, approximate computa-tion of trajectories ( Using subsets of the data) may be beneficial. Tuning of HMC canbe made easier Using a short-cut in which trajectories computed with a bad choice ofstepsize take little computation time. Finally, tempering methods may be useful whenmultiple isolated modes Hamiltonian DynamicsHamiltonian Dynamics has a physical interpretation that can provide useful two dimensions, we can visualize the Dynamics as that of a frictionless puck that slidesover a surface of varying height.

8 The state of this system consists of thepositionof the puck,given by a two-dimensional vectorq, and themomentumof the puck (its mass times itsvelocity), given by a two-dimensional vectorp. Thepotential energy,U(q), of the puck isproportional to the height of the surface at its current position, and itskinetic energy,K(p),is equal to|p|2/(2m), wheremis the mass of the puck. On a level part of the surface, thepuck moves at a constant velocity, equal top/m. If it encounters a rising slope, the puck smomentum allows it to continue, with its kinetic energy decreasing and its potential energyincreasing, until the kinetic energy (and hencep) is zero, at which point it will slide backdown (with kinetic energy increasing and potential energy decreasing).

9 In nonphysical mcmc applications of Hamiltonian Dynamics , the position will cor-respond to the variables of interest. The potential energy will be minus the log of theprobability density for these variables. Momentum variables, one for each position variable,will be introduced interpretations may help motivate the exposition below, but if you find otherwise,the Dynamics can also be understood as simply resulting from a certain set of Hamilton s EquationsHamiltonian Dynamics operates on ad-dimensionalpositionvector,q, and ad-dimensionalmomentumvector,p, so that the full state space has 2ddimensions.

10 The system is describedby a function ofqandpknown as theHamiltonian,H(q,p). Equations of MotionThe partial derivatives of the Hamiltonian determine howqandpchange over time,t,according to Hamilton s equations: mcmc Using Hamiltonian Dynamics115dqidt= H pi,( )dpidt= H qi,( )fori=1,..,d. For any time interval of durations, these equations define a mapping,Ts,from the state at any timetto the state at timet+s. (Here,H, and henceTs, are assumed tonot depend ont.)Alternatively, we can combine the vectorsqandpinto the vectorz=(q,p)with 2ddimensions, and write Hamilton s equations asdzdt=J H(z),where His the gradient ofH( [ H]k= H/ zk), andJ=00d dId d Id d0d d1( )isa2d 2dmatrix whose quadrants are defined above in terms of identity and zero Potential and Kinetic EnergyFor HMC we usually use Hamiltonian functions that can be written asH(q,p)=U(q)+K(p).