Hamiltonian Monte Carlo

Suppose the target distribution to sample is ${\displaystyle f(\mathbf {x} )}$ and a chain of samples ${\displaystyle \mathbf {X} _{0},\mathbf {X} _{1},\mathbf {X} _{2},\ldots }$ is required. The Hamilton's equations are

${\displaystyle {\frac {{\text{d}}x_{i}}{{\text{d}}t}}={\frac {\partial H}{\partial p_{i}}}}$

and

${\displaystyle {\dfrac {{\text{d}}p_{i}}{{\text{d}}t}}=-{\dfrac {\partial H}{\partial x_{i}}}}$

where ${\displaystyle x_{i}}$ and ${\displaystyle p_{i}}$ are the ${\displaystyle i}$th component of the position and momentum vector respectively and ${\displaystyle H}$ is the Hamiltonian. Let ${\displaystyle M}$ be a mass matrix which is symmetric and positive definite, then the Hamiltonian is

${\displaystyle H(\mathbf {x} ,\mathbf {p} )=U(\mathbf {x} )+{\dfrac {1}{2}}\mathbf {p} ^{\text{T}}M^{-1}\mathbf {p} }$

where ${\displaystyle U(\mathbf {x} )}$ is the potential energy. The potential energy for a target is given as

${\displaystyle U(\mathbf {x} )=-\ln f(\mathbf {x} )}$

which comes from the Boltzmann's factor.

The algorithm requires a positive integer for number of leap frog steps ${\displaystyle L}$ and a positive number for the step size ${\displaystyle \Delta t}$. Suppose the chain is at ${\displaystyle \mathbf {X} _{n}=\mathbf {x} _{n}}$. Let ${\displaystyle \mathbf {x} _{n}(0)=\mathbf {x} _{n}}$. First, a random Gaussian momentum ${\displaystyle \mathbf {p} _{n}(0)}$ is drawn from ${\displaystyle {\text{N}}\left(\mathbf {0} ,M\right)}$. Next, the particle will run under Hamiltonian dynamics for time ${\displaystyle L\Delta t}$, this is done by solving the Hamilton's equations numerically using the leap frog algorithm. The position and momentum vectors after time ${\displaystyle \Delta t}$ using the leap frog algorithm are

${\displaystyle \mathbf {p} _{n}\left(t+{\dfrac {\Delta t}{2}}\right)=\mathbf {p} _{n}(t)-{\dfrac {\Delta t}{2}}\nabla \left.U(\mathbf {x} )\right|_{\mathbf {x} =\mathbf {x} _{n}(t)}}$

${\displaystyle \mathbf {x} _{n}(t+\Delta t)=\mathbf {x} _{n}(t)+\Delta tM^{-1}\mathbf {p} _{n}\left(t+{\dfrac {\Delta t}{2}}\right)}$

${\displaystyle \mathbf {p} _{n}(t+\Delta t)=\mathbf {p} _{n}\left(t+{\dfrac {\Delta t}{2}}\right)-{\dfrac {\Delta t}{2}}\nabla \left.U(\mathbf {x} )\right|_{\mathbf {x} =\mathbf {x} _{n}(t+\Delta t)}}$

These equations are to be applied to ${\displaystyle \mathbf {x} _{n}(0)}$ and ${\displaystyle \mathbf {p} _{n}(0)}$ ${\displaystyle L}$ times to obtain ${\displaystyle \mathbf {x} _{n}(L\Delta t)}$ and ${\displaystyle \mathbf {p} _{n}(L\Delta t)}$.

Because the leap frog algorithm is a numerical method and does not solve the Hamilton's equations exactly, a Metropolis–Hastings step is used. The transition from ${\displaystyle \mathbf {X} _{n}=\mathbf {x} _{n}}$ to ${\displaystyle \mathbf {X} _{n+1}}$ is

${\displaystyle \mathbf {X} _{n+1}|\mathbf {X} _{n}=\mathbf {x} _{n}={\begin{cases}\mathbf {x} _{n}(L\Delta t)&{\text{with probability }}\alpha \left(\mathbf {x} _{n}(0),\mathbf {x} _{n}(L\Delta t)\right)\\\mathbf {x} _{n}(0)&{\text{otherwise}}\end{cases}}}$

where

${\displaystyle \alpha \left(\mathbf {x} _{n}(0)_{n},\mathbf {x} _{n}(L\Delta t)\right)={\text{min}}\left(1,{\dfrac {\exp \left[-H(\mathbf {x} _{n}(L\Delta t),\mathbf {p} _{n}(L\Delta t))\right]}{\exp \left[-H(\mathbf {x} _{n}(0),\mathbf {p} _{n}(0))\right]}}\right)}$

This is repeated to obtain ${\displaystyle \mathbf {X} _{n+1},\mathbf {X} _{n+2},\mathbf {X} _{n+3},\ldots }$.

The No U-Turn Sampler (NUTS)[2] is an extension by controlling ${\displaystyle L}$ automatically. Tuning ${\displaystyle L}$ is critical. For example in the one dimensional ${\displaystyle {\text{N}}(0,1/{\sqrt {k}})}$ case, the potential is ${\displaystyle U(x)=kx^{2}/2}$ which corresponds to the potential of a simple harmonic oscillator. For ${\displaystyle L}$ too large, the particle will oscillate and this waste computational time. For ${\displaystyle L}$ 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. When that happens, a random point from the path is chosen for the MCMC sample and the process is repeated from that new point.

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. A slice variable ${\displaystyle U_{n}\sim {\text{Uniform}}(0,\exp(-H[\mathbf {x} _{n}(0),\mathbf {p} _{n}(0)]))}$ is sampled. Let ${\displaystyle \mathbf {x} _{n}^{+}}$ and ${\displaystyle \mathbf {p} _{n}^{+}}$ be the position and momentum of the forward particle respectively. Similarly, ${\displaystyle \mathbf {x} _{n}^{-}}$ and ${\displaystyle \mathbf {p} _{n}^{-}}$ for the backward particle. 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. For example, in the first iteration, the forward particle moves forwards in time using 1 leap frog step. In the next iteration, the backward particle moves backwards in time using 2 leap frog steps.

The iterative procedure continues until the U-Turn condition is met, that is

${\displaystyle (\mathbf {x} _{n}^{+}-\mathbf {x} _{n}^{-})\cdot \mathbf {p} _{n}^{-}<0\quad {\text{or}}\quad (\mathbf {x} _{n}^{+}-\mathbf {x} _{n}^{-})\cdot \mathbf {p} _{n}^{+}<0}$

or when the Hamiltonian becomes inaccurate

${\displaystyle \exp \left[-H(\mathbf {x} _{n}^{+},\mathbf {p} _{n}^{+})+\delta \right]<U_{n}}$

or

${\displaystyle \exp \left[-H(\mathbf {x} _{n}^{-},\mathbf {p} _{n}^{-})+\delta \right]<U_{n}}$

where, for example, ${\displaystyle \delta =1000}$.

Once the U-Turn condition is met, the next MCMC sample, ${\displaystyle \mathbf {x} _{n+1}}$, is obtained by sampling uniformly the leap frog path traced out by the binary tree ${\displaystyle \{\mathbf {x} _{n}^{-},\ldots ,\mathbf {x} _{n}(-\Delta t),\mathbf {x} _{n}(0),\mathbf {x} _{n}(\Delta t),\ldots ,\mathbf {x} _{n}^{+}\}}$ which satisfies

${\displaystyle U_{n}<\exp \left[-H(\mathbf {x_{n+1}} ,\mathbf {p_{n+1})} \right]}$

This is usually satisfied if the remaining HMC parameters are sensible.

• Dynamic Monte Carlo method
• Software for Monte Carlo molecular modeling
• Stan

• Neal, Radford M (2011). "MCMC Using Hamiltonian Dynamics" (PDF). In Steve Brooks; Andrew Gelman; Galin L. Jones; Xiao-Li Meng (eds.). Handbook of Markov Chain Monte Carlo. Chapman and Hall/CRC. ISBN 9781420079418.
• Betancourt, Michael (2018). "A Conceptual Introduction to Hamiltonian Monte Carlo". arXiv:1701.02434. Bibcode:2017arXiv170102434B. Cite journal requires |journal= (help)
• Liu, Jun S. (2004). Monte Carlo Strategies in Scientific Computing. Springer Series in Statistics, Springer. pp. 189-203. ISBN 978-0-387-76369-9.

• Betancourt, Michael. "Efficient Bayesian inference with Hamiltonian Monte Carlo". MLSS Iceland 2014 – via YouTube.
• Hamiltonian Monte Carlo from scratch
• Optimization and Monte Carlo Methods