1. Motivation: Why Hamiltonian Monte Carlo?

Standard MCMC methods like:

  • Gibbs sampling
  • Random-walk Metropolis

often behave like a random walk:

  • The chain moves in small, local steps.
  • It “zigzags” slowly through the parameter space.
  • In high dimensions, this can be very inefficient.

Even with:

  • good reparameterization, and
  • good Metropolis proposals,

the local random-walk behavior remains, especially for complicated, high-dimensional posteriors.

Hamiltonian Monte Carlo (HMC) is designed to fix exactly this issue:

It suppresses random-walk behavior and lets the Markov chain move smoothly and quickly through the space, following something like “physical trajectories.”

That’s why HMC is often much more efficient than random-walk methods.

2. Basic Idea of HMC: Positions and Momenta

Let:

  • θ = (θ₁, …, $\theta_d$) be the parameters you care about (position variables).
  • φ = (φ₁, …, $φ_d$) be auxiliary momentum variables added by HMC.

So HMC works in a 2d-dimensional space (θ, φ).

  • You define a joint density:
    • $p(\theta, \phi \mid y) = p(\theta \mid y) \, p(\phi)$
  • The posterior $p(\theta \mid y)$ is your usual target.
  • The momentum distribution $p(\phi)$ is artificial; it is introduced just to help the algorithm move more effectively.

We simulate from this joint distribution, but we only keep θ. φ is an auxiliary variable.

This construction allows HMC to:

  • use φ to “push” θ along long trajectories,
  • move through parameter space in a coherent way rather than taking tiny random steps.

This is why HMC is sometimes called hybrid Monte Carlo: it combines:

  • MCMC (Metropolis accept–reject),
  • with deterministic simulation of Hamiltonian dynamics.

3. Gradient Requirement

HMC needs not just the log-posterior, but also its gradient.

  • You need to be able to compute:
    • $\nabla_{\theta} \log p(\theta \mid y) = \left( \frac{\partial \log p(\theta \mid y)}{\partial \theta_1}, \dots, \frac{\partial \log p(\theta \mid y)}{\partial \theta_d} \right)$

This must be computed analytically in practice:

  • Numerical finite differences would require too many function evaluations and be too slow.

Recommended practice:

  • Implement the analytic gradient.
  • Also implement a numerical gradient (finite differences) just for debugging.
  • Check that they match to several decimal places.
    • If they don’t, there is almost certainly a bug in your gradient.

For most standard Bayesian models, deriving the gradient is not too hard.

4. The Momentum Distribution $p(\phi)$ and Mass Matrix $M$

HMC typically chooses:

  • $\phi \sim N(0, M)$

where:

  • M is the mass matrix (by analogy to physics),
  • M is a d × d positive definite matrix,
  • often taken as diagonal for simplicity.

If M is diagonal:

  • $\phi_j \sim N(0, M_{jj})$ independently for j = 1,…,d.

Good practice:

  • It’s helpful if M roughly matches the inverse posterior covariance:
    • $M \approx \left[\text{Var}(\theta \mid y)\right]^{-1}$
    • because then the dynamics better match the geometry of the posterior.
  • However, HMC still works even if M is just the identity; better scaling just improves efficiency.

5. One HMC Iteration: The Three Main Steps

Each HMC iteration consists of:

  1. Resample momentum φ
  2. Simulate Hamiltonian dynamics via leapfrog steps
  3. Metropolis accept–reject

Let’s go through them.

Step 1: Sample φ

At the start of each iteration:

  1. Sample $\phi \sim N(0, M)$.

Since φ is independent of θ (in the joint density), its conditional and prior are the same.

This “refreshes” the momentum randomly at every iteration, ensuring exploration is not deterministic and preserving the Markov property.

Step 2: Leapfrog updates of (θ, φ)

This is the core of HMC.

We perform L leapfrog steps, each with step size ε.

In each leapfrog step, θ and φ are updated in a way that approximates Hamiltonian dynamics. Each iteration does:

  1. Half-step for φ using the gradient of log p(θ | y):
    • $\phi \leftarrow \phi + \frac{\epsilon}{2} \nabla_{\theta} \log p(\theta \mid y)$
  2. Full step for θ using φ:
    • $\theta \leftarrow \theta + \epsilon M^{-1} \phi$
  3. Another half-step for φ using the new θ:
    • $\phi \leftarrow \phi + \frac{\epsilon}{2} \nabla_{\theta} \log p(\theta \mid y)$

Except at the very beginning and end, step 3 of one leapfrog and step 1 of the next can be merged, but conceptually it’s “half-step–full-step–half-step”.

Why this structure?

  • The leapfrog integrator is symplectic and time-reversible.
  • It approximately conserves the “Hamiltonian” (energy), so the joint density $p(\theta, \phi \mid y)$ stays approximately constant.
  • In the limit ε → 0, the leapfrog dynamics preserve the joint density exactly.

Intuition:

  • If θ is in a flat region of the posterior:
    • The gradient is near zero.
    • φ remains almost unchanged.
    • So θ moves in (almost) a straight line with constant velocity.
  • If θ moves toward a region of lower posterior density:
    • The gradient points back toward higher density, so φ is updated to slow down.
    • The kinetic energy (log p(φ)) decreases to compensate for the drop in log p(θ).
    • Eventually the trajectory slows, stops, and curves back or around the low-density region.
  • If θ moves toward a region of higher posterior density:
    • The gradient pushes φ to increase in that direction.
    • The trajectory accelerates into the high-density region, then eventually curves and slows past the mode.

This behavior allows HMC to take long, coherent, curved trajectories through the posterior, rather than small random wiggles.

Discretization error:

  • For finite ε, energy is not preserved exactly; there is small error each step.
  • Leapfrog has the nice property that the error tends to oscillate rather than accumulate linearly.
  • Keeping ε small keeps this error under control.

Step 3: Metropolis Accept–Reject

Let:

  • $(\theta_{t-1}, \phi_{t-1})$: state at the start of the iteration.
  • $(\theta^\*, \phi^\*)$: state after L leapfrog steps.

Compute:

$r = \frac{p(\theta^\* \mid y) p(\phi^\*)}{p(\theta_{t-1} \mid y) p(\phi_{t-1})}$

Then:

  • With probability $\min(r, 1)$, accept: $\theta_t = \theta^\*$
  • Otherwise, reject and keep the previous state: $\theta_t = \theta_{t-1}$

φ doesn’t matter after this step because it will be resampled in the next iteration, so we don’t bother to store $φ_t$.

We repeat this procedure across iterations until convergence, monitoring:

  • $\hat{R}$ (R-hat) near 1,
  • effective sample size large enough.

6. Constraints and Zero-Density Regions

HMC assumes the target density is positive wherever we move.

If during a leapfrog trajectory, θ enters a region where the posterior density is zero, e.g.:

  • a parameter that must be positive goes negative,
  • violating a constraint like θ ∈ (0, 1),

then:

  1. Simple strategy (reject whole trajectory):
    • Stop the leapfrog steps.
    • Reject this proposal.
    • Stay at the previous θ for that iteration.
    • This preserves detailed balance.
  2. Bouncing strategy:
    • If a step goes outside the allowed region,
    • Flip the sign of the momentum φ in that direction (bounce back).
    • This also preserves detailed balance.
    • Typically more efficient than rejecting the whole trajectory, especially with hard boundaries.
  3. Transformations:

Instead of working with constrained parameters directly, you can map them to an unconstrained space:

  • Positive θ: use log θ.
  • θ ∈ (0, 1): use logit(θ).
  • Ordered parameters (θ₁ < θ₂ < θ₃), or simplex constraints (α₁ + … + α₄ = 1), use more complex bijective transforms.

If you transform:

  • You must include the Jacobian of the transformation in the log-posterior,
  • and compute the gradient in the transformed space correctly.

This is what most modern HMC implementations (e.g. Stan) do internally.

7. Tuning HMC: M, ε, and L

HMC has three main tuning components:

  1. Mass matrix M (defines the momentum distribution $p(\phi)$)
    • Often diagonal.
    • Each diagonal element sets the scale of φ in that dimension.
  2. Step size ε
    • How far you move in one leapfrog micro-step.
  3. Number of leapfrog steps L
    • How long each trajectory is.

7.1 Warm-up (adaptation) vs sampling

Adaptive tuning is tricky:

  • If you adapt ε, L, or M based on past iterations while collecting samples for inference, you can destroy detailed balance and the chain might no longer converge to the correct target distribution.
  • So the usual approach is:
    • Have a warm-up (or burn-in) phase, during which:
      • You tune M, ε, and L based on preliminary runs.
      • You are allowed to adapt here.
    • Discard all samples from this warm-up phase.
    • Then fix these parameters.
    • Run HMC again (or continue) and use only the samples from the fixed-parameter phase for inference.

You can repeat this warm-up and retuning multiple times as long as you only keep samples after the last adaptation.

7.2 Practical initial choices

  1. Mass matrix M
    • Start with a rough estimate of the scale of θ.
    • If no information, use the identity matrix (diagonal with ones).
  2. Set the product εL ≈ 1
    • This means: if you start near the center of the posterior, one trajectory of length ≈1 in “scaled units” should roughly go across the bulk of the distribution.
    • A basic default:
      • ε = 0.1, L = 10 (so εL = 1).
  3. Target acceptance rate ≈ 65%
    • Theory suggests HMC is most efficient when the acceptance probability is around 65%.
    • This is analogous to the 23% optimal rate for high-dimensional random-walk Metropolis.
    • If average acceptance is too low:
      • Leapfrog steps are too aggressive → reduce ε and increase L (keeping εL ≈ 1).
    • If average acceptance is too high:
      • Steps are too conservative → increase ε and decrease L.

These rules are not perfect but provide a reasonable guideline.

8. Varying Tuning Parameters During the Run

There are two different kinds of changes:

8.1 Dangerous: adaptation based on past acceptance during sampling

If we adapt ε using acceptance outcomes during the sampling phase, for example:

  • increase ε after accepted proposals,
  • decrease ε after rejected proposals,

then:

  • we would violate detailed balance,
  • we would no longer be guaranteed that the stationary distribution is the correct posterior.
  • So such adaptation must be restricted to the warm-up phase only.

8.2 Safe: purely random variation within fixed ranges

If we randomly vary ε and L from iteration to iteration (within chosen ranges), independently of acceptance history:

  • This does not break detailed balance.
  • It can help:
    • sometimes taking long trajectories,
    • sometimes short ones,
    • giving the chain different kinds of movement.

The downside:

  • You lose some optimality—sometimes the chain will attempt long trajectories in cramped regions, or short ones when long would be fine.
  • But correctness (convergence to the true posterior) is preserved.

9. Locally Adaptive HMC: NUTS and Riemannian Adaptation

For difficult posteriors, we’d like HMC to adapt locally:

  • M adjusts to local curvature (like a local metric).
  • ε gets smaller where curvature is large (to maintain precision).
  • L is large enough to make meaningful moves but not so large that the trajectory loops around.

Two advanced methods:

9.1 No-U-Turn Sampler (NUTS)

NUTS modifies HMC by choosing L adaptively in each iteration:

  • Instead of fixing L, it simulates forward (and backward) along the trajectory.
  • It stops when the trajectory starts to turn back on itself:
    • Formally, when the dot product between the current momentum φ and the displacement from the starting θ becomes negative.

This “no-U-turn” condition:

  • Prevents wasting computation making the trajectory loop around.
  • Effectively sends the trajectory as far as it can go in a productive direction.

However:

  • If you only apply the naïve stopping rule, it would not preserve detailed balance.
  • The full NUTS algorithm is more complex, building a binary tree of states and carefully choosing samples to maintain detailed balance.

NUTS also typically includes a warm-up procedure that automatically:

  • tunes the mass matrix M,
  • tunes step size ε.

After warm-up, these are held fixed.

9.2 Riemannian adaptation

Here, the mass matrix M is allowed to vary with θ, adapting to the local curvature of the log-posterior, like a Riemannian metric:

  • In regions with strong curvature, M adapts to match that geometry.
  • This can greatly improve efficiency for posteriors with highly varying curvature or non-exponential tails.

However:

  • The algorithms become more complex.
  • Computing and updating a position-dependent metric is expensive.
  • This is often impractical in very high-dimensional problems.

Riemannian adaptation can be combined with NUTS for even more flexibility.

Limitations:

  • NUTS is efficient and self-tuning but can still struggle with extremely short- or long-tailed distributions.
  • Riemannian methods handle curvature and tail behavior better but don’t scale well to very high dimensions.

10. Combining HMC with Gibbs Sampling and Discrete Parameters

There are two main ways Gibbs-style ideas appear in an HMC context.

10.1 Blocking structure in hierarchical models

Consider a hierarchical model with parameters:

$\theta = (\eta^{(1)}, \dots, \eta^{(J)}, \phi)$

where:

  • $\eta^{(j)}$ are group-specific parameter vectors,
  • φ is a vector of hyperparameters,
  • Posterior factorizes as:
    • $p(\theta \mid y) \propto p(\phi) \prod_{j=1}^{J} p(\eta^{(j)} \mid \phi) p(y^{(j)} \mid \eta^{(j)})$

Even if you could update all of θ at once using HMC, it might be more efficient to block:

  • Use HMC to update each $\eta^{(j)}$ given φ and $y^(j)$,
  • Then HMC (or another method) to update φ given all $\eta^{(j)}$,
  • Cycle through these J+1 updates.

Advantages:

  • Each step only needs one likelihood factor $p(y^{(j)} \mid \eta^{(j)})$,
  • Potentially cheaper per update,
  • Better mixing in complicated hierarchical structures,
  • You can also combine with parameter expansion to improve mixing further.

10.2 Discrete parameters

Hamiltonian dynamics work only for continuous variables.

If some parameters are discrete, for example:

  • mixture component indicators,
  • spike-and-slab indicators (parameters that are exactly zero with some probability),

then:

  • HMC cannot update them directly.
  • You can use:
    • Gibbs sampling (if full conditionals are available),
    • Metropolis updates,
    • slice sampling, etc., for the discrete variables.

A common approach:

  1. Partition parameters into:
    • Continuous part: use HMC.
    • Discrete part: use Gibbs / Metropolis / slice.
  2. Alternate:
    • HMC steps on the continuous parameters given current discrete values,
    • Discrete updates given the current continuous parameters.

This hybrid scheme respects the correct joint posterior and often works well in practice.

Final Summary

Hamiltonian Monte Carlo:

  • Adds momentum variables φ and simulates Hamiltonian dynamics in (θ, φ).
  • Uses leapfrog integration plus a Metropolis accept–reject step.
  • Requires gradients of the log-posterior.
  • Avoids random-walk behavior and works especially well for high-dimensional continuous posteriors.
  • Needs careful tuning (M, ε, L), usually done in a warm-up phase.
  • Can be extended and improved via:
    • NUTS (automatic trajectory length),
    • Riemannian metrics (local curvature adaptation),
    • Blocking and combination with Gibbs updates for hierarchical and/or discrete-parameter models.