Hamiltonian Monte Carlo

The Hamiltonian Monte Carlo algorithm (originally known as hybrid Monte Carlo) is a Markov chain Monte Carlo method for obtaining a sequence of random samples whose distribution converges to a target probability distribution that is difficult to sample directly.

This sequence can be used to estimate integrals of the target distribution, such as expected values and moments.

Hamiltonian Monte Carlo corresponds to an instance of the Metropolis–Hastings algorithm, with a Hamiltonian dynamics evolution simulated using a time-reversible and volume-preserving numerical integrator (typically the leapfrog integrator) to propose a move to a new point in the state space.

Compared to using a Gaussian random walk proposal distribution in the Metropolis–Hastings algorithm, Hamiltonian Monte Carlo reduces the correlation between successive sampled states by proposing moves to distant states which maintain a high probability of acceptance due to the approximate energy conserving properties of the simulated Hamiltonian dynamic when using a symplectic integrator.

The reduced correlation means fewer Markov chain samples are needed to approximate integrals with respect to the target probability distribution for a given Monte Carlo error.

The algorithm was originally proposed by Simon Duane, Anthony Kennedy, Brian Pendleton and Duncan Roweth in 1987 for calculations in lattice quantum chromodynamics.

[1] In 1996, Radford M. Neal showed how the method could be used for a broader class of statistical problems, in particular artificial neural networks.

[2] However, the burden of having to provide gradients of the Bayesian network delayed the wider adoption of the algorithm in statistics and other quantitative disciplines, until in the mid-2010s the developers of Stan implemented HMC in combination with automatic differentiation.

be a mass matrix which is symmetric and positive definite, then the Hamiltonian is where

The algorithm requires a positive integer for number of leapfrog steps

Next, the particle will run under Hamiltonian dynamics for time

, this is done by solving the Hamilton's equations numerically using the leapfrog algorithm.

The leapfrog algorithm is an approximate solution to the motion of non-interacting classical particles.

In order to reach a thermodynamic equilibrium distribution, particles must have some sort of interaction with, for example, a surrounding heat bath, so that the entire system can take on different energies with probabilities according to the Boltzmann distribution.

One way to move the system towards a thermodynamic equilibrium distribution is to change the state of the particles using the Metropolis–Hastings algorithm.

is where A full update consists of first randomly sampling the momenta

(independently of the previous iterations), then integrating the equations of motion (e.g. with leapfrog), and finally obtaining the new configuration from the Metropolis-Hastings accept/reject step.

The No U-Turn Sampler (NUTS)[5] is an extension by controlling the number of steps

which corresponds to the potential of a simple harmonic oscillator.

too large, the particle will oscillate and thus waste computational time.

too small, the particle will behave like a random walk.

Loosely, NUTS runs the Hamiltonian dynamics both forwards and backwards in time randomly until a U-Turn condition is satisfied.

In detail, a binary tree is constructed to trace the path of the leap frog steps.

To produce a MCMC sample, an iterative procedure is conducted.

be the position and momentum of the forward particle respectively.

In each iteration, the binary tree selects at random uniformly to move the forward particle forwards in time or the backward particle backwards in time.

Also for each iteration, the number of leap frog steps increase by a factor of 2.

The iterative procedure continues until the U-Turn condition is met, that is or when the Hamiltonian becomes inaccurate or where, for example,

Once the U-Turn condition is met, the next MCMC sample,

, is obtained by sampling uniformly the leap frog path traced out by the binary tree