Hamiltonian Monte Carlo For Bayesian Regression

Why Use Hamiltonian Monte Carlo?

A Bayesian regression posterior can be difficult to explore when its coefficients are correlated or when the posterior mass lies in a narrow region of a continuous parameter space. A random-walk Metropolis proposal moves locally and may spend many iterations taking small steps. Hamiltonian Monte Carlo (HMC) uses the gradient of the log posterior to construct longer, directed proposals. It can make substantial moves while accounting for the shape of the posterior.

HMC is a simulation method, not a different posterior target. Its proposals are generated using a numerical approximation to a dynamical system, then a Metropolis step corrects for integration error. The method requires a differentiable log posterior with gradients that can be evaluated at proposed parameter values. The examples below use Gaussian regression with known error variance, but the same construction applies to other differentiable regression posteriors.

The Posterior Target and Its Gradient

Let \(y\) be the vector of outcomes, \(X\) the design matrix, and \(\beta\) the vector of regression coefficients. For a Gaussian linear regression with known error variance, suppose \(y\mid\beta,\sigma^2 \sim N(X\beta,\sigma^2 I)\), where \(\sigma^2\) is fixed. A prior density \(p(\beta)\) gives the posterior target \(\pi(\beta\mid y)\propto p(y\mid\beta)p(\beta)\). Constants that do not depend on \(\beta\) can be omitted when calculating the gradient or Metropolis ratio.

For the Gaussian regression model, the log posterior, up to an additive constant, is \(\ell(\beta)= -\tfrac{1}{2\sigma^2}(y-X\beta)^\mathsf{T}(y-X\beta)+\log p(\beta)\). Differentiating with respect to \(\beta\) gives \(\nabla_\beta\ell(\beta)=X^\mathsf{T}(y-X\beta)/\sigma^2+\nabla_\beta\log p(\beta)\). The gradient points locally toward increasing log posterior density; the potential-energy gradient has the opposite sign, \(\nabla U(\beta)=-\nabla_\beta\ell(\beta)\). For a Gaussian prior \(\beta\sim N(m,V)\), the prior contribution is \(-V^{-1}(\beta-m)\).

Hamiltonian Dynamics and Leapfrog Integration

In ideal Hamiltonian dynamics, the position \(q\) and momentum \(p\) evolve according to \(\frac{dq}{dt}=M^{-1}p\) and \(\frac{dp}{dt}=-\nabla U(q)\). The momentum is refreshed from \(N(0,M)\) to start a trajectory. If these differential equations could be followed exactly, the Hamiltonian would stay constant, so movement could travel through a posterior region without systematically favoring high- or low-density states. In practice, the equations are integrated numerically.

HMC commonly uses the leapfrog integrator. For step size \(\epsilon>0\), one step consists of a half momentum update, a full position update, and a second half momentum update. Repeating this operation \(L\) times gives a proposed endpoint. The approximate trajectory length is \(\tau=L\epsilon\).

Leapfrog is useful because it is reversible (with momentum reversal) and preserves volume in the joint position-momentum space. Each momentum update is a shear in \(p\), and each position update is a shear in \(q\); each has Jacobian determinant one. The symmetric half-step, full-step, half-step sequence reverses when run backward with momentum negated. These properties make it possible to apply a Metropolis correction to the numerical proposal rather than treating the numerical path as exact.

After \(L\) steps, let the endpoint be \((q^\ast,p^\ast)\). The proposal is accepted with probability \(\min\{1,\exp[-H(q^\ast,p^\ast)+H(q_0,p_0)]\}\). Equivalently, if \(\Delta H=H(q^\ast,p^\ast)-H(q_0,p_0)\), the probability is \(\min(1,\exp(-\Delta H))\). If rejected, the position remains \(q_0\). A momentum sign reversal is often included in the proposal description to make reversibility explicit; because kinetic energy depends quadratically on momentum, the sign does not change the acceptance calculation.

Choosing Step Size and Trajectory Length

The step size \(\epsilon\) controls the numerical resolution of the path. If it is too large, the integrator may make substantial energy errors, causing frequent rejections or unstable trajectories. A smaller step size generally improves numerical accuracy but requires more gradient evaluations for the same trajectory length. The mass matrix \(M\) also matters: it sets the scale of momentum and can help accommodate different scales and correlations among parameters. For example, a well-chosen mass matrix can reduce inefficient movement when coefficients have very different posterior scales.

The number of steps \(L\), together with \(\epsilon\), determines how far the trajectory runs. If \(L\) is too small, proposals may not travel far from their starting values. If it is too large, a trajectory can spend extra computation retracing or circling through regions already visited. The product \(\tau=L\epsilon\) is a useful description of integration length, but the best length depends on the posterior geometry and is not determined by a universal formula.

Applying HMC to Regression Parameters

For a regression with several continuous coefficients, the position vector \(q\) can simply be \(\beta\). If an unknown variance or another continuous parameter is also sampled, it may be included in \(q\) after expressing it on an unconstrained scale. For example, a positive variance can be represented by its logarithm; the target density on the transformed scale must include the appropriate change-of-variables Jacobian. The complete target must be differentiable in the coordinates used for HMC.

In an econometric model, the gradient is valuable because it combines information from all observations into a direction of local posterior increase. That direction can make HMC substantially more efficient than proposals that ignore local posterior geometry, particularly when coefficients are strongly related. But gradients do not remove the need to specify the likelihood and prior correctly, and the acceptance step cannot repair a misspecified model. HMC samples from the posterior density supplied to it, subject to the accuracy of implementation and simulation.

Common Mistakes and Practical Checks

• Reversing the gradient sign. The log-posterior gradient points uphill in log posterior; the potential-energy gradient is its negative. In the leapfrog momentum update, use \(-\nabla U(q)\), equivalently the log-posterior gradient.
• Mixing up parameter order. If \(q=(\beta_0,\beta_1)^\mathsf{T}\), the first gradient entry is for the intercept and the second is for the slope. Keep this order consistent in the design matrix, vectors, and interpretation.
• Assuming leapfrog conserves energy exactly. It approximates the Hamiltonian dynamics. The Metropolis acceptance probability corrects for the resulting energy error.
• Changing step size without considering trajectory length. Since \(\tau=L\epsilon\), changing \(\epsilon\) while holding \(L\) fixed changes both numerical resolution and total path length.
• Treating a high acceptance fraction as proof of good sampling. Acceptance is only one tuning indicator. Short, frequently accepted proposals may still move little; convergence and mixing require separate assessment.