Hamiltonian Monte Carlo With ODE Likelihoods

Why Use Gradients in ODE Calibration?

Bayesian calibration of a QSP model can produce a posterior distribution with strong parameter dependence and narrow, curved regions of high probability. A random-walk proposal may move only a short distance before it leaves such a region, or it may propose large numbers of low-probability states. Hamiltonian Monte Carlo (HMC) uses the direction of the log-posterior gradient to construct longer, informed proposals. In an ODE calibration, that gradient must account for how parameters change the model trajectory and therefore the predicted observations.

HMC does not remove the need to define a suitable likelihood, prior, or parameter scale. It is a sampling method for a specified posterior, and its gradient calculations must accurately represent that posterior. The likelihood construction and prior choices are part of the model; here the focus is how their combined log density guides HMC.

From an ODE to a Log-Posterior Gradient

Suppose the mechanistic state \(x(t)\) satisfies \(dx/dt=f(x(t),q,t)\), with a specified initial condition. At observation time \(t_i\), the model predicts \(m_i(q)=h(x(t_i),q)\). For independent Gaussian observations \(y_i\) with known standard deviations \(\sigma_i\), the log-likelihood, omitting constants that do not depend on \(q\), is

For one observation, differentiation by the chain rule gives

Thus, the gradient combines the observation residual, its uncertainty, and the sensitivity of the predicted observation to the parameters. If \(q\) has prior density \(p(q)\), then \(\ell(q)=\log L(q)+\log p(q)\), again up to an additive constant, and \(\nabla_q\ell(q)=\nabla_q\log L(q)+\nabla_q\log p(q)\). If a transformation is used to represent a constrained parameter on an unconstrained scale, the appropriate transformed prior density includes its change-of-variables term.

ODE sensitivities can be computed by differentiating the state equations. For parameter \(q_j\), define \(S_j(t)=\partial x(t)/\partial q_j\). When the model is differentiable, this sensitivity satisfies

The observation sensitivity follows by differentiating \(h(x(t_i),q)\). These equations explain where the likelihood gradient comes from; in practice, automatic differentiation, sensitivity equations, or adjoint calculations may be used. Whatever method is chosen, the resulting gradient must correspond to the same numerical model and log density used to evaluate \(\ell(q)\).

Hamiltonian Dynamics and the Leapfrog Proposal

HMC gives momentum \(p\) a multivariate normal distribution, commonly \(p\sim N(0,M)\), where \(M\) is a positive-definite mass matrix. Define the potential energy \(U(q)=-\ell(q)\) and kinetic energy \(K(p)=p^\mathsf{T}M^{-1}p/2\). The Hamiltonian, or total energy, is \(H(q,p)=U(q)+K(p)\). In exact continuous dynamics, the position and momentum evolve according to

The second equation shows how the log-posterior gradient changes the momentum. HMC cannot usually solve these dynamics exactly, so it uses the leapfrog integrator, which alternates half momentum updates with position updates. For step size \(\epsilon\), one step is

Repeating these steps produces a proposed endpoint. The method then applies a Metropolis correction: given initial and proposed states \((q_0,p_0)\) and \((q_1,p_1)\), the probability of accepting the proposal is

Exact Hamiltonian dynamics conserve energy. Numerical leapfrog steps only approximate that property, so the correction compensates for their energy error. With a reversible, volume-preserving integrator and the momentum reversal used in the proposal construction, this correction preserves the target posterior. Rejected proposals leave the position at its current value. Momentum is refreshed for a subsequent trajectory.

Gradient Checks and Numerical Reliability

A gradient error can systematically steer HMC in the wrong direction, and an ODE solver can make an otherwise correct derivative numerically unreliable. A basic check compares an analytic or automatic gradient with a centered finite difference. For component \(j\), the latter is

where \(e_j\) is the unit vector in coordinate \(j\). The finite-difference step \(h\) should be large enough to avoid domination by numerical noise and small enough to approximate the local derivative. For ODE likelihoods, both evaluations must use sufficiently accurate and consistent solver settings.

Practical Diagnostics, Alternatives, and Limits

HMC performance depends on accurate gradients, a suitable parameter scale, and an integration step size that resolves the posterior geometry. A step size that is too large can produce substantial energy errors and frequent rejections. A step size that is extremely small may yield high acceptance but require many ODE solves to travel a useful distance. Poor scaling can make movement easy in one parameter direction and difficult in another. These symptoms should prompt checks of parameterization, gradient calculations, solver tolerances, and the mass matrix; acceptance alone does not establish that a chain explores the posterior well.

For correlated mechanistic parameters, a joint HMC trajectory can use gradient information to move through their dependence, unlike a simple coordinate-wise random walk. A random-walk Metropolis proposal is easier to implement and does not require derivatives, but may explore a narrow posterior inefficiently. Blocked Metropolis updates can respect parameter groupings, yet still rely on proposal distributions rather than gradients to choose directions. HMC is not automatically preferable: it can be difficult when the posterior is nonsmooth, gradients are unavailable or noisy, or ODE computations are too costly. It also does not guarantee that distinct posterior modes will be found.

A common mistake is to assume that a high acceptance rate proves HMC is working well. A very small step size can produce high acceptance while making only tiny moves; conversely, poor gradients can undermine exploration even when some proposals are accepted. Check gradient agreement and energy errors, and assess the sampled chain with appropriate convergence and efficiency diagnostics. Those checks answer different questions: numerical correctness of the proposal is not the same as adequate posterior exploration.