Hamiltonian Monte Carlo

Remarks

The proof works on the joint density proportional to $e^{-H(q,p)}$ on $\R^{2d}$, whose $q$-marginal is $\pi$ and whose $p$-marginal is standard normal. Let $R(q,p) = (q,-p)$; then $F = R \circ \Phi$ is an involution, since each leapfrog step is reversed by flipping the momentum, and $F$ preserves volume, since each of the three substeps is a shear. A Metropolis step with a deterministic, volume-preserving, involutive proposal $F$ and acceptance probability $\min(1, e^{H - H \circ F})$ satisfies detailed balance for $e^{-H}$, and $H \circ R = H$. Resampling $p$ in step 1 also leaves $e^{-H}$ invariant, so the composition does, and hence its $q$-marginal $\pi$ is invariant. The leapfrog map is a symplectic integrator for the Hamiltonian system $\dot q = p$, $\dot p = -\nabla U(q)$. The exact flow would keep $H$ constant and always accept, and the accept step corrects the discretisation error. Example: for $U(q) = \tfrac{1}{2} q^2$ on $\R$, the target is standard normal, and with small $\varepsilon$ the trajectories follow circles in the $(q,p)$ plane, so the chain moves far in each step. Invariance alone does not give convergence. For this $U$, if $\varepsilon L$ is half a period of the exact flow, the chain nearly alternates between $q$ and $-q$. For this reason $\varepsilon$ or $L$ is often randomised in practice. More generally $p$ may be drawn from a normal law with a positive definite mass matrix $M$, with $\tfrac{1}{2} p^T M^{-1} p$ in $H$ and $q' = q + \varepsilon M^{-1} p_{1/2}$.

Going deeper

For a student

HMC is a Markov chain Monte Carlo method that uses the gradient of the log-density to make long, directed moves. It is a special case of the Metropolis–Hastings algorithm on the extended space of positions and momenta. The two properties of the leapfrog map that make the acceptance rule simple are reversibility and volume preservation; energy conservation is only needed for a high acceptance rate. Compared with random-walk proposals, HMC scales much better with dimension for smooth targets, which is why it is the default sampler in several statistics packages. Choosing $\varepsilon$ and $L$ is the main practical difficulty, and adaptive schemes such as the no-U-turn sampler address it.