An Introduction to Hamiltonian Monte Carlo


[Running Title]Hamiltonian Monte Carlo
A novel proposal function
choosing method


[Yifan: edit 2013-11-22] There are many algorithm which can deals with high dimensional sampling. Some are easy to be implemented such as Gibbs, slice Gibbs, block Gibbs, Metropolis Hasting algorithm (Robert 2004), and particle filtering (Arulampalam 2002). But when the target function is highly correlated, Gibbs and HM with Gaussian proposal usually perform poorly (ref). In this manuscript, we mainly described the newly developed MC method Hamiltonian Monte Carlo (HMC). The origin of this algorithm could trace back to molecular dynamics (MD) simulation (Alder 1959) and (Rahman 1964). MD is a computer simulation technique that allows one to predict the time evolution of a system of interacting particles based on classical mechanical Hamiltonian equations, and becomes a very popular procedure for physical research. It is a deterministic technique, which means given initial positions and velocities, the evolution of the system in time is, in principle, completely determined. Besides, it is a statistical mechanics method generating a set of configurations that are distributed according to statistical distribution functions. By ergodecity theory (ref, time average will cause configuration average in large partial system), MC simulations fits well with MD simulations. Basically, MD has been shown work well by nature as a result of its reliance on basis Newtonian mechanics. But it can not be used directly to determine proposal function in Metropolis Hasting algorithm (ref). To overcome this disadvantage, (Duane 1987) introduced Hybrid Monte Carlo. It take the advantage of MD and MH at the same time. On one side, it derives a new position using Hamilton equation in MD. on the other hand, HMC take the advantage of HM type of accept-rejection principle to draw Monte Carlo samples.

Hamiltonian Monte Carlo, as described previously, is highly reliable in high dimensional highly correlated cases. It offers a reasonable proposal function in Metropolis Hasting algorithm. We will first look at some basics of Newtonian mechanics and Hamiltonian mechanics. Then we will introduce Hamiltonian/Hybrid Monte Carlo with leapfrog approach. Some basic properties would be proved. In addition, some tuning methods on leapfrog method would be discussed based on some simulation results.

Math and Equations

5 pages

Hamiltonian Dynamics

There are mainly three ways to describe class mechanics, i.e. Newtonian, Lagrangian and Hamiltonian mechanics. They are all equivalent, but the equations of motion for a given system may appear simpler in one of the approaches than in the others. In some particular fields, such as quantum physics and accelerator physics, Hamiltonian equations have some advantages. Given a function \(H(x,p;t)\) (called Hamiltonian), the equations of motion for a dynamical system are given by Hamilton’s equations:

\[\begin{aligned} %\label{eq:Hamiltonian} \frac{\partial H}{\partial x_i}=- \frac{d p_i(t)}{d t}\\ \frac{\partial H}{\partial p_i}= \frac{d x_i(t)}{d t}\end{aligned}\]

, for all \(i=1,\ldots,d\), where \(x\in \mathbf{R}^n\) is the coordinates vector (hence \(\dot x\) is the velocity vector) and \(p\in \mathbf{R}^n\) is the momentum vector. In simple cases, \(p_i\) is equivalent to the mechanical momentum \(m \frac{d x_i(t)}{d t}\), but this is not always the case. The Hamiltonian plays as the same role in Hamiltonian mechanics as force plays in Newtonian mechanics. It defines a \(2\times{}n\) first oder partial derivative equations instead of second order partial derivative equations defined in Lagrangian mechanics (ref). One import advantage is that Hamiltonian systems are “conservative systems” (ref, Liouville’s Theorem), which keeps conservation of area in phase space. In simple example, Hamiltonian can be written as

\[%label{eq:HamilTotal} H(x,p)=U(x)+K(p)\]

for kinetic energy \(K(p)\) and potential energy \(U(x)\). From this equation, we could find Hamiltonian is very similar to the “total” energy of the system, expressed in terms of the coordinates and conjugate momentum. In fact, this phenomenon ensures Hamiltonian could be extended to more complicated field (ref), which we will not discuss here. In this manuscript, we mainly concentrated on the form (xxx). Back to statistics, we treat

\[\begin{array}{c@{ as }l} U(x)&-loglik(x\mbox{that we wish to sample})+Constant\\ K(p)&-loglik(N(p;0,M))+Constant \end{array}\]

or more straightforward: \[H(x,p)=-logLik(x)-logLik(p), \mbox{assuming }p\sim N_n(0,M)\]

Here \(M>0\) is a “mass matrix” is often chosen as a diagonal matrix. Hence the Hamiltonian we used in this manuscript is reduced to \[\begin{aligned} -\frac{\partial H}{\partial x_i}=\frac{d p_i(t)}{d t}=-\frac {d U(x_i)}{d x_i}\\ \frac{\partial H}{\partial p_i}=\frac{d x_i(t)}{d t}=\frac{d K(p_i)}{d p_i}=[M^{-1}p]_i\end{aligned}\]

Since it choose the direction adaptively, we must verify the ergodecity very carefully. Following theorem ensures this properties. We described the theorem and proved some of them.



Gibbs sampling


2 pages