Hamiltonian Monte Carlo (HMC) Mechanics

While standard Metropolis-Hastings and Gibbs sampling provide foundational approaches for MCMC, they can struggle in high-dimensional or highly correlated posterior distributions. Their proposals often resemble random walks, leading to slow progression and high autocorrelation between samples. Hamiltonian Monte Carlo (HMC) offers a significantly more efficient way to sample the parameter space by borrowing ideas from classical physics, specifically Hamiltonian dynamics.

Think of the negative log posterior density as a potential energy surface. We want to explore this surface efficiently. HMC introduces a fictitious particle moving on this surface. The position of this particle corresponds to the parameters $\theta$ we want to sample, and we introduce an auxiliary momentum variable $p$, typically with the same dimension as $\theta$.

The Hamiltonian Framework

In physics, the total energy of a system, the Hamiltonian $H$, is the sum of potential energy $U$ and kinetic energy $K$. We adapt this for our statistical purpose:

1. Potential Energy $U(\theta)$: This is defined as the negative log of the target density. Since MCMC only requires the target density up to a constant of proportionality, we can use the negative log of the joint probability (likelihood times prior), ignoring the evidence term:

   $$U(\theta) = -\log[p(D|\theta)p(\theta)] = -\log p(\theta, D)$$

   Minimizing potential energy corresponds to maximizing the posterior probability.

2. Kinetic Energy $K(p)$: This depends on the auxiliary momentum variables $p$ and a "mass matrix" $M$, which is typically assumed to be the identity matrix ($M=I$) for simplicity, unless we have reasons to scale dimensions differently.

   $$K(p) = \frac{1}{2} p^T M^{-1} p$$

   If $M=I$, this simplifies to $K(p) = \frac{1}{2} \sum_{i} p_i^2$. This corresponds to sampling the momentum from a multivariate Gaussian distribution, usually $p \sim \mathcal{N}(0, M)$.

3. Hamiltonian $H(\theta, p)$: The total "energy" of our fictitious system is:

   $$H(\theta, p) = U(\theta) + K(p)$$

The core idea is to simulate the particle's movement according to Hamiltonian dynamics. In an ideal physical system, the total energy $H(\theta, p)$ remains constant over time. The particle's trajectory is governed by Hamilton's equations:

$$\frac{d\theta_i}{dt} = \frac{\partial H}{\partial p_i} = (M^{-1}p)_i$$ $$\frac{dp_i}{dt} = -\frac{\partial H}{\partial \theta_i} = -\frac{\partial U(\theta)}{\partial \theta_i}$$

Intuitively, the first equation states that momentum causes the position (parameters) to change. The second equation states that the force acting on the particle (the negative gradient of the potential energy, which is the gradient of the log posterior) causes the momentum to change. Notice that calculating the change in momentum requires the gradient of the log posterior, $\nabla_{\theta} U(\theta)$. This is an important requirement for HMC.

Simulating Dynamics with the Leapfrog Integrator

Solving Hamilton's equations analytically is generally impossible for complex posteriors. Instead, we use numerical integration. A popular and effective method is the leapfrog integrator. It updates position and momentum in interleaved steps, providing better stability and preserving geometric properties compared to simpler integrators like Euler's method.

The leapfrog algorithm takes discrete time steps of size $\epsilon$. To move from time $t$ to $t + \epsilon$, it performs three updates:

1. Half-step update for momentum: $$p(t + \epsilon/2) = p(t) - \frac{\epsilon}{2} \nabla_{\theta} U(\theta(t))$$
2. Full-step update for position: $$\theta(t + \epsilon) = \theta(t) + \epsilon M^{-1} p(t + \epsilon/2)$$
3. Half-step update for momentum: $$p(t + \epsilon) = p(t + \epsilon/2) - \frac{\epsilon}{2} \nabla_{\theta} U(\theta(t + \epsilon))$$

We repeat these steps $L$ times to simulate the trajectory for a total duration of $L\epsilon$.

The HMC Algorithm

An iteration of HMC proceeds as follows:

1. Start with the current parameter state $\theta_{current}$.
2. Sample Momentum: Draw a random momentum vector $p$ from its Gaussian distribution, e.g., $p \sim \mathcal{N}(0, M)$.
3. Simulate Trajectory: Starting from $(\theta_{current}, p)$, run the leapfrog integrator for $L$ steps with step size $\epsilon$. Let the state at the end of the trajectory be $(\theta', p')$.
4. Metropolis Acceptance Step: Because the leapfrog integrator introduces small numerical errors, the Hamiltonian $H$ is not perfectly conserved. We correct for this using a Metropolis acceptance step. Calculate the acceptance probability: $$\alpha = \min\left(1, \frac{\exp(-H(\theta', p'))}{\exp(-H(\theta_{current}, p))}\right) = \min\left(1, \frac{p(\theta', D)\exp(-K(p'))}{p(\theta_{current}, D)\exp(-K(p))}\right)$$ Note that the target density $p(\theta, D)$ appears directly here.
5. Accept or Reject: Draw $u \sim \text{Uniform}(0, 1)$. If $u < \alpha$, accept the proposed state: $\theta_{next} = \theta'$. Otherwise, reject the proposal: $\theta_{next} = \theta_{current}$.
6. Return: The next sample in the Markov chain is $\theta_{next}$. The momentum variable $p$ is discarded; it's only used to generate the proposal. A fresh momentum is drawn at the start of the next iteration.

An illustration of an HMC trajectory moving across contours of a potential energy surface (representing the negative log posterior). The trajectory moves through different parameter values while tending to stay within regions of similar energy.

Why HMC Works Well

HMC uses the gradient information $\nabla_{\theta} U(\theta)$ to guide its proposals intelligently. Instead of a random walk, it proposes moves informed by the shape of the posterior. By simulating dynamics for multiple steps $L$, HMC can generate proposals far from the current state that still have a high probability of acceptance, especially if $\epsilon$ is chosen appropriately. This drastically reduces the autocorrelation between samples and allows for much faster exploration of the target distribution compared to random-walk Metropolis or Gibbs sampling, particularly when dealing with many parameters or strong correlations between them.

Practical Notes

• Gradient Calculation: HMC requires the gradient of the log posterior. Modern probabilistic programming languages (PPLs) like Stan, PyMC, NumPyro, and TensorFlow Probability use automatic differentiation to compute these gradients efficiently and accurately, making HMC practical for complex models.
• Tuning Parameters: HMC introduces two main tuning parameters: the step size $\epsilon$ and the number of steps $L$.
  • $\epsilon$: Controls the discretization error. Too large $\epsilon$ leads to low acceptance rates due to inaccurate Hamiltonian simulation. Too small $\epsilon$ makes the simulation slow and requires many steps $L$ to move significantly.
  • $L$: Controls the length of the simulated trajectory. Too small $L$ results in behavior similar to a random walk (high correlation). Too large $L$ can lead to the trajectory doubling back on itself and wasted computation. Finding optimal $\epsilon$ and $L$ often requires some experimentation or adaptive methods.

HMC represents a significant step forward in MCMC technology. Its ability to use gradients for efficient exploration makes it a workhorse algorithm for many modern Bayesian inference problems. Understanding its mechanics provides insight into how these advanced samplers operate and sets the stage for understanding adaptive extensions like the No-U-Turn Sampler (NUTS), which automates the tuning of $L$.