Why Divergences Matter
Hamiltonian Monte Carlo (HMC) can move efficiently through a posterior distribution by using gradients to construct long, informed proposals. But its trajectories are computed numerically, not exactly. If the numerical path fails to follow the posterior’s geometry, the resulting transition can be unreliable. A divergent transition is an HMC diagnostic warning that this has happened: the numerical trajectory encountered an energy error large enough to meet the sampler’s divergence rule, or a non-finite numerical failure occurred.
In a clinical-trial model, a divergence might arise near a narrow posterior region, a strong parameter dependency, or a boundary. It is not a finding about the treatment effect, and it does not by itself prove that the model is wrong. It does mean that the sampler may not be exploring some part of the posterior reliably. A small number of divergences can matter if they cluster in a region that affects an important posterior quantity.
Energy Error and Numerical Trajectories
Let \(q\) denote the model parameters on an unconstrained scale and \(p\) the auxiliary momentum. The target density for \(q\) is proportional to \(\exp[-U(q)]\), where \(U(q)\) is the negative log posterior up to an additive constant. With a positive-definite mass matrix \(M\), the kinetic energy is \(K(p)=\frac{1}{2}p^\mathsf{T}M^{-1}p\), and the Hamiltonian is \(H(q,p)=U(q)+K(p)\). HMC draws momentum, follows the Hamiltonian equations approximately using a numerical integrator, and uses the resulting endpoint to propose a new state.
For an exact trajectory, the Hamiltonian is conserved. A numerical integrator introduces an energy error. Here define \(\Delta H=H_{\mathrm{end}}-H_{\mathrm{start}}\). The Metropolis correction accounts for ordinary integration error, but a sufficiently large error indicates that the trajectory is not providing a trustworthy exploration step. In Stan, the default divergence threshold is 1000 for a positive energy error; non-finite numerical failures are also flagged. This is an implementation-specific diagnostic convention, not a universal mathematical boundary between safe and unsafe trajectories. In particular, a large negative \(\Delta H\) is not equivalent to a large positive \(\Delta H\) under this one-sided threshold check.
Worked Example: Reading an Energy Error
Suppose a hypothetical HMC transition has starting Hamiltonian \(H_{\mathrm{start}}=18.4\) and ending Hamiltonian \(H_{\mathrm{end}}=1019.1\). Apply the stated Stan threshold rule.
The energy error is
This is a positive energy error above 1000, so it meets the default threshold rule and is flagged as divergent. Now suppose a different transition has \(H_{\mathrm{start}}=18.4\) and \(H_{\mathrm{end}}=16.9\). Then \(\Delta H=16.9-18.4=-1.5\). That error is negative and does not meet the positive-error threshold rule. This comparison does not make every unflagged transition accurate; it only applies the stated diagnostic rule to these values.
Why Strong Curvature Challenges HMC
HMC’s numerical integration is most demanding where posterior geometry changes rapidly. One way to see the role of curvature is to consider a one-dimensional Gaussian-shaped potential \(U(q)=\frac{1}{2}\omega^2q^2\). Its curvature is \(\omega^2\), and the Hamiltonian equations are \(\dot q=p\) and \(\dot p=-\omega^2q\). The position changes rapidly when \(\omega\) is large. A step size that works in a broad direction may therefore be too large in a narrow direction.
For this harmonic example, a leapfrog update with step size \(\epsilon\) has a stable oscillatory update only when \(\epsilon\omega<2\). To see the condition, write one update as a matrix mapping the pair \((q,p)\) to its next value. The matrix has determinant 1 and trace \(2-(\epsilon\omega)^2\). Its eigenvalues lie on the unit circle, as required for stable oscillation, when the trace is between \(-2\) and \(2\), which gives \(0<\epsilon\omega<2\). This is a simplified stability result, not a universal cutoff for posterior models; it illustrates why a high-curvature direction can defeat a step size chosen for a less constrained region.
Worked Example: A Step Size That Is Too Large for a Narrow Direction
Consider a hypothetical local posterior direction approximated by the harmonic potential above, with \(\omega=20\). Compare leapfrog step sizes \(\epsilon=0.12\) and \(\epsilon=0.08\).
For \(\epsilon=0.12\), the dimensionless product is \(\epsilon\omega=0.12(20)=2.4\), which exceeds 2. The simplified harmonic update is outside its stable oscillatory range, so this step size is unsuitable for that direction.
For \(\epsilon=0.08\), the product is \(0.08(20)=1.6\), which is below 2. The update falls inside the example’s stability range. That does not guarantee accurate sampling: stability is weaker than accuracy, and a real posterior can have changing curvature and many dimensions. The calculation does show why reducing the step size can help when narrow curvature is the main problem.
Strong curvature can occur in several ways. Parameters may have different scales, posterior correlations may create a long, narrow ridge, or a hierarchical model may have a funnel-shaped region. A likelihood can also become sharply informative in some directions, while a prior or a boundary constrains other directions. These are clues to investigate, not diagnoses that can be made from the model’s subject matter alone.
Investigating and Responding to Divergences
Start by confirming that the divergences are present and locating where they occur. Inspect the sampler’s divergence indicator alongside parameter draws, particularly pairwise plots of the parameters that define the model’s geometry. Divergences concentrated near a boundary or in a narrow part of a parameter relationship are more informative than a count alone. Review whether parameters and predictors are on sensible scales, whether constrained parameters have been represented appropriately, and whether the model’s prior and likelihood produce a well-defined posterior in the region being sampled.
A lower step size can improve numerical resolution. Warmup adaptation can also select a more suitable step size and mass matrix, which scales movement across parameter directions. These are reasonable first responses when the model is well specified and the problem appears to be integration resolution or scale. If divergences persist or cluster in a distinctive geometric region, reparameterization may be more effective. For example, standardizing a predictor can reduce avoidable scale differences. In some hierarchical models, a non-centered parameterization can reduce funnel geometry. Any such change must preserve the intended probability model, including the correct transformed density when a parameter transformation requires a Jacobian adjustment.
Worked Example: Reparameterizing a Hierarchical Trial Model
Suppose a hypothetical trial has treatment effects \(\theta_j\) for four sites. A model describes these effects as \(\theta_j\sim N(\mu,\tau^2)\), with a proper prior on the between-site standard deviation \(\tau\). The centered representation can have a funnel-like geometry when posterior values of \(\tau\) are near zero: the allowable variation in the \(\theta_j\) values becomes very small there. Assume the model’s likelihood and priors are otherwise valid and differentiable on the sampler’s unconstrained representation.
A non-centered representation writes
This is a re-expression of the same hierarchical relationship: conditional on \(\mu\) and \(\tau\), each \(\theta_j\) has mean \(\mu\) and standard deviation \(\tau\). The transformed density must be implemented consistently; the sampler must target the same model rather than silently changing the prior.
Imagine four chains each produce 1,000 post-warmup draws. In a hypothetical centered run, 96 of 4,000 transitions are flagged, or \(96/4000=0.024=2.4\%\). After fitting the non-centered representation, 8 of 4,000 are flagged, or \(8/4000=0.002=0.2\%\). The reduction is evidence that the new representation has improved this diagnostic for these runs, not proof that the computation is fully reliable. The remaining flagged draws and other sampling diagnostics still require review, and the choice between centered and non-centered forms depends on the posterior geometry and data information.
Distinguishing Tuning, Reparameterization, and Model Problems
A practical response should match the suspected cause. If the posterior is otherwise sensibly expressed and divergences appear to result from coarse integration, a smaller step size may help, at the cost of more computation. If the difficulty comes from scale differences, an appropriate mass matrix or rescaling may help. If divergences persist in a narrow or funnel-shaped region, reparameterization is often worth considering. If investigation reveals a coding error, an unintended prior, or an implausible model structure, tuning is not a substitute for correcting it.
These approaches are not interchangeable. Tuning changes how finely the integrator follows a trajectory; it does not remove the underlying geometry. Reparameterization changes coordinates while preserving the target distribution when done correctly. Revising a prior or likelihood changes the statistical model and should be justified scientifically, not used merely to suppress a warning. If HMC remains unsuitable, another sampling method may be considered, but it has its own assumptions and diagnostics. Gibbs sampling, Metropolis methods, and specialized samplers do not make computational validation unnecessary.
The absence of divergences is reassuring but not conclusive. A sampler may fail to visit an important region, and other problems can remain even when no transition crosses the divergence threshold. Conversely, a divergence does not establish that every posterior summary is unusable. The relevant question is whether the computation has explored the parts of the posterior needed for the intended summaries, and whether the diagnostics support that conclusion. In particular, do not dismiss a small number of divergences without checking their locations and potential influence.
Common Mistakes and Limitations
- Treating the threshold as a universal safety boundary. The default value is a software diagnostic convention. It is not a theorem that transitions below it are reliable or transitions above it identify a particular model defect.
- Using absolute energy error when applying Stan’s default threshold. The stated check is for positive energy error above the threshold, with non-finite numerical failures also flagged. A negative error of large magnitude is not automatically equivalent under that rule.
- Assuming a smaller step size solves every problem. It may reduce integration error, but persistent geometry problems can remain and may call for reparameterization or model review.
- Changing the model to remove a warning. Changing a prior or likelihood changes the target posterior. Such changes need substantive justification and should not be treated as mere computational tuning.
- Interpreting zero divergences as proof of correct inference. This only reports that this particular warning was not triggered; it does not establish that the chains explored all relevant posterior regions.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions.
- A transition starts at Hamiltonian 7.2 and ends at 1008.0. What is \(\Delta H\), and how does it compare with the stated default threshold?
- For a harmonic direction with \(\omega=15\), calculate \(\epsilon\omega\) for \(\epsilon=0.10\). What does the simplified stability condition indicate?
- Why might a non-centered representation help a hierarchical model with a funnel-like posterior region?
- If a model has no flagged divergences, what can and cannot be concluded?
- When is changing a prior an inappropriate response to divergences?
Solutions. (1) \(\Delta H=1008.0-7.2=1000.8\), a positive error above 1000, so it meets the stated Stan threshold rule. (2) \(\epsilon\omega=0.10(15)=1.5<2\); this is inside the simplified harmonic stability range, but does not guarantee accuracy in a real posterior. (3) Writing \(\theta_j=\mu+\tau z_j\) with standard-normal \(z_j\) can separate the latent variation from the scale \(\tau\), reducing the funnel geometry that challenges numerical integration; the transformed model must still target the same distribution. (4) It indicates that this divergence warning was not triggered in the run, but does not prove complete exploration or validate every posterior summary. (5) Changing a prior is inappropriate when the sole aim is to silence a computational warning; it changes the statistical model and requires substantive justification.