Tutorials › Nonlinear Mixed-Effects Models (NLME) › Step Size Selection in Hamiltonian Monte Carlo

Bayesian NLME with MCMC · Tutorial 246 of 797

Step Size Selection in Hamiltonian Monte Carlo

Understand how to choose and interpret an HMC step size by balancing numerical accuracy, acceptance, and the cost of computing each trajectory.

Advanced 12 min read

What You'll Learn

  • Define the HMC step size and its role in leapfrog integration.
  • Calculate Hamiltonian error and the corresponding Metropolis acceptance probability.
  • Derive a stability limit for a Gaussian target with a single curvature scale.
  • Compare the computational cost of smaller steps at fixed integration time.
  • Recognize why a high acceptance rate does not by itself establish efficient sampling.
  • Identify practical checks for step-size choices in Bayesian NLME models.

Why the Step Size Matters

Hamiltonian Monte Carlo (HMC) uses gradients of the log posterior to propose movement through parameter space. Its numerical trajectory is built from repeated leapfrog updates. The step size, written \(\epsilon>0\), is the amount of simulated time advanced at each update. It controls how finely the trajectory is resolved: a smaller step generally follows the continuous Hamiltonian path more accurately, but requires more gradient evaluations to cover the same simulated time.

In a Bayesian NLME model, the parameter vector can contain population parameters, random effects, and residual-error parameters. Posterior curvature can differ greatly across directions and locations. The step size must therefore work adequately across the geometry visited by the chain, rather than only near one convenient point. The mass matrix and parameter scaling affect that geometry; here the focus is how \(\epsilon\) controls numerical integration once those choices are fixed.

Definition: Let \(q\) denote the parameters and let \(U(q)=-\log \pi(q\mid y)\), up to an additive constant, be the negative log posterior density. HMC augments \(q\) with momentum \(p\), with \(p\sim N(0,M)\) for a positive-definite mass matrix \(M\). The Hamiltonian is $$ H(q,p)=U(q)+\frac{1}{2}p^\mathsf{T}M^{-1}p. $$ A leapfrog step of size \(\epsilon\) numerically advances this Hamiltonian system. For a proposed endpoint with Hamiltonian change \(\Delta H=H(q_{\mathrm{new}},p_{\mathrm{new}})-H(q_{\mathrm{old}},p_{\mathrm{old}})\), the Metropolis acceptance probability is $$ \alpha=\min\{1,\exp(-\Delta H)\}. $$

The target posterior is preserved by the Metropolis correction when the proposal construction is reversible and volume-preserving and the acceptance rule is applied correctly. Leapfrog has these properties under the usual smoothness conditions needed for its gradients and updates. This correction does not make a poor step size harmless: inaccurate trajectories tend to produce larger positive energy errors and more rejections, while unstable updates can produce unusable proposals.

Accuracy, Energy Error, and Cost

The exact Hamiltonian flow conserves \(H\). Leapfrog only approximates that flow, so its endpoint usually has a nonzero \(\Delta H\). A positive \(\Delta H\) means the proposed endpoint has higher Hamiltonian and is accepted with probability below one. A negative \(\Delta H\) is accepted with probability one, although it still signals integration error. Thus acceptance is not a direct measure of endpoint accuracy in every individual trajectory.

For a smooth Hamiltonian, leapfrog has local error of order \(\epsilon^3\) per step and global error of order \(\epsilon^2\) over a fixed integration time, provided the numerical solution remains in a stable regime. These orders describe how error scales as the step is made small; they do not give a universal numerical error bound for an arbitrary posterior. The actual error depends on the trajectory and the posterior's curvature.

Key relationship: If a trajectory spans a fixed simulated time \(T\), it uses approximately \(L=T/\epsilon\) leapfrog steps. Halving \(\epsilon\) therefore roughly doubles the number of gradient evaluations for that trajectory, while reducing the leading integration error by a factor of about four in the small-step regime. This is an approximation, not a guarantee of a fourfold reduction in every realized energy error.

A useful stability calculation comes from a one-dimensional Gaussian target. Suppose \(U(q)=\tfrac12\omega^2q^2\), with unit mass, so that the Hamiltonian equations describe an oscillator with angular frequency \(\omega\). One leapfrog step consists of a half momentum update, a full position update, and another half momentum update. In matrix form, $$ \begin{pmatrix}q_{n+1}\\p_{n+1}\end{pmatrix} = \begin{pmatrix} 1-\epsilon^2\omega^2/2 & \epsilon\\ -\epsilon\omega^2(1-\epsilon^2\omega^2/4) & 1-\epsilon^2\omega^2/2 \end{pmatrix} \begin{pmatrix}q_n\\p_n\end{pmatrix}. $$

The update matrix has determinant one and trace \(2-\epsilon^2\omega^2\). Its eigenvalues remain on the unit circle when the trace lies strictly between \(-2\) and \(2\). Thus this oscillator requires \(0<\epsilon\omega<2\) for stability. This is a stability condition for the simple quadratic example, not a recommended universal step size. In a multivariate or non-Gaussian posterior, high-curvature directions and changing local curvature can require a smaller step.

Worked Calculations

Worked Example: One Leapfrog Step and Its Energy Error

Consider a hypothetical one-dimensional Gaussian target with \(U(q)=2q^2=\tfrac12(2^2)q^2\), unit mass, and starting state \(q_0=0,\ p_0=1\). Take \(\epsilon=0.20\). Compute one leapfrog step and its acceptance probability.

The half-step momentum update is $$ p_{1/2}=p_0-\frac{\epsilon}{2}(4q_0) =1-\frac{0.20}{2}(4)(0)=1. $$ The position update gives \(q_1=q_0+\epsilon p_{1/2}=0+0.20(1)=0.20\). The final half-step momentum is $$ p_1=p_{1/2}-\frac{\epsilon}{2}(4q_1) =1-0.10(4)(0.20)=0.92. $$ Initially, \(H_0=2(0)^2+\tfrac12(1)^2=0.5000\). At the endpoint, $$ H_1=2(0.20)^2+\frac12(0.92)^2=0.0800+0.4232=0.5032, \qquad \Delta H=0.5032-0.5000=0.0032. $$ Therefore \(\alpha=\exp(-0.0032)\approx0.9968\). This particular proposal is accepted with high probability. It does not establish that \(\epsilon=0.20\) is adequate for all starting states or for a more complicated posterior.

Worked Example: Curvature Restricts a Stable Step

Suppose a hypothetical posterior is locally approximated, in one direction, by a quadratic potential with curvature \(\lambda=100\), so \(U(q)\approx\tfrac12\lambda q^2\), and suppose the mass in that direction is one. The oscillator frequency is \(\omega=\sqrt{\lambda}=10\). The stability condition from the quadratic calculation is \(\epsilon\omega<2\), giving $$ \epsilon<\frac{2}{10}=0.20. $$ A candidate step \(\epsilon=0.10\) satisfies the bound because \(0.10(10)=1<2\). A step \(\epsilon=0.22\) does not, because \(0.22(10)=2.2>2\), and is unstable for this quadratic direction. The first value is not automatically optimal: other directions, nonlinear changes in curvature, and the desired trajectory length also matter. The calculation shows why a step that seems modest on one parameter scale can be too large in a sharply curved direction.

Worked Example: Interpreting Energy Errors and Average Acceptance

Suppose a hypothetical set of four HMC proposals, generated with the same candidate step size, has energy changes \(\Delta H=0.02,\ 0.10,\ 0.35,\ -0.05\). Their acceptance probabilities are, respectively, $$ e^{-0.02}\approx0.9802,\quad e^{-0.10}\approx0.9048,\quad e^{-0.35}\approx0.7047,\quad 1, $$ because the last energy change is negative. The mean of these four acceptance probabilities is $$ \frac{0.9802+0.9048+0.7047+1}{4} =\frac{3.5897}{4}\approx0.8974. $$ This is the average acceptance probability across these four hypothetical proposals, not a guarantee that a future proposal will be accepted at that rate. The positive errors show that some endpoints had higher Hamiltonian; the negative error is still an integration discrepancy, even though its proposal is accepted with probability one. A larger collection of trajectories and checks of numerical stability are needed to judge a step size in practice.

Choosing and Checking a Step Size

Step-size selection balances at least three considerations: numerical stability, the size of Hamiltonian errors, and computational cost. A step that is too large may yield inaccurate or unstable trajectories. A very small step usually reduces discretization error, but at fixed simulated trajectory length it requires more leapfrog steps and gradient evaluations. If trajectories are instead kept to a fixed number of steps, reducing \(\epsilon\) also shortens their simulated length; the resulting proposals may move less far through parameter space. Step size should therefore be interpreted alongside how trajectory length is set.

In NLME posteriors, a small step may be needed to resolve a narrow, high-curvature region even if much of the posterior is relatively broad. The practical issue is not simply the average curvature: a step suitable for broad regions can behave poorly in a sharply curved part of the posterior. Conversely, shrinking the step to accommodate every local feature can be expensive. Step size is a numerical tuning choice, not a change to the intended posterior distribution when the Metropolis correction and HMC conditions are satisfied.

Useful checks include examining acceptance probabilities and energy errors across a representative range of trajectories, and checking for numerical failures or extreme energy changes. A consistently low acceptance rate can indicate steps that are too large, but a high acceptance rate alone does not prove that proposals travel far enough or that the chain explores the posterior well. Sampling efficiency also depends on trajectory length, posterior geometry, and the mass matrix. The later material on adaptation and HMC diagnostics develops those assessments further.

Step-size choiceLikely numerical consequencePractical trade-off
Too large for local curvatureLarge energy errors or instabilityMore rejections and potentially unreliable numerical trajectories
Moderate for the geometry visitedControlled integration errorBalances gradient cost and proposal quality
Very smallUsually finer integrationMore gradients for a fixed trajectory time; may be inefficient

Compared with random-walk Metropolis, HMC uses gradient information to construct directed proposals, but it introduces numerical integration error and a step-size choice. Within HMC, reducing \(\epsilon\) is a direct way to refine leapfrog integration, while changing the mass matrix addresses parameter scaling and changing trajectory length alters how far a proposal travels. These are distinct adjustments; one cannot generally substitute for another. HMC's Metropolis correction protects the target distribution under its conditions, but it does not guarantee efficient exploration or make a poorly resolved trajectory computationally worthwhile.

Common Pitfall

Wrong reasoning: “The average acceptance probability is high, so the step size must be good.” This overlooks computational cost and the distance covered by proposals. A tiny step can produce very accurate, highly accepted proposals while using many gradients and moving only a short distance if the number of steps is held fixed.

Correction: Interpret acceptance together with energy errors, numerical stability, gradient cost, and trajectory length. Also remember that a negative \(\Delta H\) has acceptance probability one but is not evidence of zero integration error. Neither acceptance nor energy error alone is a complete measure of sampling performance.

Key takeaway: The step size controls leapfrog resolution. Smaller \(\epsilon\) typically improves integration accuracy but costs more at fixed simulated time; larger \(\epsilon\) saves computation per unit simulated time only until curvature-driven errors or instability become unacceptable. Assess the choice in the posterior geometry actually visited by the chain.

Check Your Understanding

Answer each question, then compare with the brief solutions below.

  1. If a proposal has \(\Delta H=0.20\), calculate its acceptance probability and explain the sign of the energy change.
  2. For a quadratic direction with curvature \(\lambda=25\) and unit mass, find the oscillator frequency and the strict stability bound on \(\epsilon\).
  3. At fixed simulated time \(T\), what happens approximately to leapfrog step count if \(\epsilon\) is halved? What happens to leading global error in the small-step regime?
  4. Why does a negative \(\Delta H\) not mean that leapfrog integrated the trajectory exactly?
  5. Give one reason why a high average acceptance probability may coexist with inefficient HMC proposals.

Solutions. (1) \(\alpha=e^{-0.20}\approx0.8187\); the endpoint Hamiltonian is higher by 0.20. (2) \(\omega=\sqrt{25}=5\), so stability requires \(\epsilon<2/5=0.40\). (3) The step count \(T/\epsilon\) approximately doubles; global error over fixed time is of order \(\epsilon^2\), so its leading scale is approximately quartered. (4) Negative \(\Delta H\) means the endpoint Hamiltonian is lower, not equal to its initial value; acceptance is one because of the Metropolis rule. (5) Very small steps can have high acceptance but require many gradients, and with a fixed number of steps they may produce short trajectories that move little.