Why Numerical Stability Matters
An individual Bayesian fit combines a patient’s concentration data with a population-based prior to form a posterior distribution for that patient’s parameters. Numerical methods are usually needed to find a posterior mode, summarize posterior uncertainty, or integrate a function of the parameters to obtain a forecast. These are distinct computations: a reliable optimizer does not guarantee accurate integration, and accurate integration cannot repair a wrongly specified model or an incorrectly coded likelihood.
For example, an optimizer may stop because its step became tiny even though the objective is still changing substantially in a poorly scaled direction. An integration routine may return a plausible prediction while missing a narrow or skewed region of posterior probability. Numerical stability therefore requires checks on the parameterization, the behavior of the objective and its derivatives, and agreement across independent computations.
Parameterization and Objective Checks
PK parameters can differ greatly in scale. Clearance may be several units, a volume may be several tens of units, and a random effect may be near zero. An optimizer that uses one step size across all parameters can take steps that are too large for one coordinate and too small for another. A poorly scaled objective can also make numerical gradients and curvature calculations unreliable.
A common remedy is to optimize over dimensionless or scaled coordinates. Positive parameters such as clearance or volume can be represented on the log scale, for example \(a=\log(CL_i/CL_{\mathrm{ref}})\), where \(CL_{\mathrm{ref}}\) is a positive reference value in the same units as clearance. This prevents invalid negative values and can make relative changes more natural than absolute changes. Random effects can be scaled by their typical variability. Such transformations often help, but they do not guarantee good conditioning, and they do not fix a flat or weakly identified direction in the likelihood and prior.
Be precise about what is being optimized. If the objective is the negative log posterior density expressed in the original coordinates, re-expressing that same scalar objective in new coordinates preserves its optimum as a mapped point. But a density is defined relative to a coordinate measure. If instead one constructs the density with respect to transformed coordinates, the change-of-variables Jacobian must be included; the numerical location of a density mode can then differ. Posterior probabilities and correctly transformed integrals are invariant when the transformation and Jacobian are handled consistently. A reported “MAP estimate” should therefore state the parameterization in which the density mode is defined.
Worked Example: Diagnosing Scale and Gradient Behavior
Suppose a hypothetical individual fit uses \(CL_i\) in L/h and \(V_i\) in L. For illustration, near the optimum the negative log posterior is $$ q(x,z)=\left(\frac{x}{0.10}\right)^2+\left(\frac{z}{0.20}\right)^2, \quad x=\log(CL_i/5),\quad z=\log(V_i/40). $$ This simple local objective has its minimum at \(x=z=0\), corresponding to \(CL_i=5\) L/h and \(V_i=40\) L.
At a candidate \(x=0.03,\ z=-0.04\), the objective is $$ q=(0.03/0.10)^2+(-0.04/0.20)^2 =0.30^2+(-0.20)^2 =0.13. $$ The gradient in \((x,z)\) coordinates is $$ \nabla q= \left(\frac{2x}{0.10^2},\frac{2z}{0.20^2}\right) =(6,-2). $$ The candidate corresponds to \(CL_i=5e^{0.03}=5.152\) L/h and \(V_i=40e^{-0.04}=38.43\) L.
The raw gradient components have different numerical magnitudes because the coordinates have different scales. Define \(u=x/0.10\) and \(v=z/0.20\). Then \(q=u^2+v^2\), and at this candidate \((u,v)=(0.30,-0.20)\), so \(\nabla q=(0.60,-0.40)\) in the scaled coordinates, with Euclidean norm \(\sqrt{0.60^2+(-0.40)^2}=0.7211\). This is not a converged solution: the objective can still decrease toward zero, and the scaled gradient is not close to zero. The example shows why a stopping report based only on a small step is insufficient; gradient size and objective change should also be considered in appropriately scaled coordinates.
Gradient, Curvature, and Boundary Diagnostics
For a differentiable objective \(q(\boldsymbol{\theta})\), an interior local minimum must have zero gradient. At a twice-differentiable interior local minimum, its Hessian, the matrix of second derivatives, must be positive semidefinite. At a twice-differentiable interior stationary point, a positive-definite Hessian is sufficient to establish a strict local minimum, but it is not necessary: a strict minimum can have a singular Hessian, as \(q(x)=x^4\) at \(x=0\) demonstrates. If the Hessian is singular or nearly singular, further analysis is needed; it may indicate a flat direction, weak information, or numerical difficulty. An indefinite Hessian at a claimed interior minimum is a warning of a saddle point, an inaccurate Hessian, or failed optimization.
These conditions apply to an interior point. If a parameter is constrained to a boundary, the gradient need not be zero in every direction; the relevant condition is that feasible small movements do not decrease the objective. Check whether the estimate is pressed against a bound and whether the boundary is scientifically meaningful or merely an algorithmic safeguard. A boundary estimate can also arise from a poorly chosen parameterization or a likelihood that provides little information.
In practice, inspect more than an optimizer’s success flag. Check the gradient in scaled coordinates, the change in objective over the final iterations, whether estimates are stable under tighter tolerances, and whether the final point satisfies the constraints. Compare analytic or automatic derivatives with finite-difference checks at several points. Finite differences themselves can be inaccurate if the step size is too small (rounding error) or too large (curvature makes the local approximation poor).
Worked Example: Interpreting Curvature Without Overclaiming
Consider two hypothetical one-parameter negative log posterior objectives near their reported optima. For \(q_1(x)=x^2\), the derivative at \(x=0\) is zero and the Hessian is \(q_1''(0)=2\), which is positive. This supports a strict local minimum. For \(q_2(x)=x^4\), the derivative at \(x=0\) is also zero, but \(q_2''(0)=0\). Since \(q_2(x)>q_2(0)\) for every \(x\ne0\), \(x=0\) is still a strict minimum.
Thus a zero or nearly zero Hessian eigenvalue does not by itself prove that the optimizer failed. It does mean that a quadratic approximation to the objective may be unreliable in that direction and that additional checks are warranted. In a patient fit, examine whether small movements along the nearly flat direction materially change the posterior objective or the forecast quantity of interest. Do not report a precise uncertainty summary solely because an optimizer returned a point estimate.
Checking Numerical Integration
For a posterior prediction or summary, the target is often an integral, not the posterior mode. For a function \(g\), the posterior expectation is $$ E[g(\boldsymbol{\theta}_i)\mid y_i,\boldsymbol{\phi}] = \int g(\boldsymbol{\theta}_i) \pi(\boldsymbol{\theta}_i\mid y_i,\boldsymbol{\phi}) \,d\boldsymbol{\theta}_i. $$ Quadrature approximates this integral using weighted evaluations; simulation-based integration approximates it using draws. Check numerical accuracy by increasing quadrature resolution or simulation size, changing integration settings, and comparing independent runs. Convergence of a mode is not evidence that this integral has converged.
For simulation-based integration, compare independent runs using the same target distribution and inspect whether the estimates of the quantities actually needed agree within their Monte Carlo error. Diagnostics such as effective sample size and between-run comparisons can reveal inadequate exploration, but no single diagnostic proves that all posterior regions have been found. For deterministic quadrature, agreement after refining the nodes is reassuring, but can be misleading if all settings miss a narrow region. In either case, compare posterior predictions as well as parameter summaries: a small parameter discrepancy can matter greatly for a nonlinear prediction.
Worked Example: Checking a Posterior Integral Against an Analytic Benchmark
Suppose a hypothetical patient has one concentration observation \(y=11\) mg/L. A simplified local model is \(y=10+2\eta+\epsilon\), where \(\eta\sim N(0,0.25)\) is a dimensionless individual effect and \(\epsilon\sim N(0,1)\) mg/L is independent residual error. The posterior for \(\eta\) is normal because both the prior and likelihood are normal.
The prior precision is \(1/0.25=4\), and the likelihood precision for \(\eta\) is \(2^2/1=4\). Therefore the posterior variance is \(1/(4+4)=0.125\). The observed residual from the model at \(\eta=0\) is \(11-10=1\) mg/L, so the posterior mean is $$ \frac{0\cdot4+2\cdot1\cdot1}{4+4}=0.25. $$ Thus \(\eta\mid y\sim N(0.25,0.125)\).
Suppose the forecast quantity is \(8e^\eta\) mg/L. For a normal random variable with mean \(m\) and variance \(v\), \(E(e^\eta)=e^{m+v/2}\). Hence the exact posterior mean forecast is $$ 8e^{0.25+0.125/2} =8e^{0.3125} =10.93\text{ mg/L}. $$ A quadrature calculation can be checked against this benchmark. A one-node rule at the posterior mean gives \(8e^{0.25}=10.27\) mg/L, which misses the contribution of posterior variance. A three-node Gauss-Hermite calculation gives approximately \(8(1.3668)=10.93\) mg/L; more precisely, the three-node quadrature forecast is about 10.9345 mg/L, close to the exact value \(8e^{0.3125}=10.93\) mg/L. The numerical agreement here is expected for this smooth, one-dimensional example; in a nonlinear, higher-dimensional PK model, refinement and independent checks remain necessary.
Independent Runs and Failure Modes
Independent runs help distinguish a reproducible result from one that depends on initialization, numerical settings, or a particular path through the computation. For optimization, use multiple dispersed starting points and compare the final objective, parameter estimates, scaled gradients, and any clinically relevant predictions. Similar objective values with materially different parameter estimates can signal a ridge or multiple modes; similar parameters and predictions across runs are more reassuring. A lower objective from one run is not automatically the correct answer if that run violates constraints or has not converged.
Worked Example: Comparing Independent Optimization Runs
A hypothetical fit to four patient concentrations uses the same likelihood, population prior, and parameter bounds in every run. The objective is negative log posterior; gradients are reported as scaled norms. The results are:
| Run | Starting point | \(CL_i\) (L/h) | \(V_i\) (L) | Objective | Scaled gradient norm |
|---|---|---|---|---|---|
| 1 | Near population typical values | 4.82 | 31.4 | 12.7310 | 0.0008 |
| 2 | Higher clearance, lower volume | 4.83 | 31.5 | 12.7312 | 0.0011 |
| 3 | Lower clearance, higher volume | 4.82 | 31.4 | 12.7311 | 0.0009 |
| 4 | Near lower clearance bound | 0.20 | 54.0 | 15.9000 | 0.0800 |
Runs 1–3 agree closely in parameters and objective, and all have small scaled gradients. Their agreement supports reproducibility of this local solution under the tested starts and settings. Run 4 has a higher objective by \(15.9000-12.7310=3.1690\), a larger gradient norm, and a clearance at the lower bound. It is not evidence for a competing optimum; it is a warning that this run did not reach the same well-converged solution. The appropriate response is to investigate the run, not to average its estimates with the others. Agreement across three runs is useful evidence, not proof that no other mode exists or that the model is scientifically adequate.
| Approach | Useful check | Main limitation |
|---|---|---|
| Multiple optimization starts | Compare objectives, estimates, constraints, and gradients | Can repeatedly find the same mode while missing another |
| Refined deterministic quadrature | Compare estimates as nodes or tolerances are increased | Can miss concentrated or irregular posterior regions |
| Independent simulation runs | Compare target summaries and Monte Carlo uncertainty across runs | Agreement does not establish exploration of every relevant region |
| Local Hessian assessment | Identify flat, ill-conditioned, or locally non-minimum behavior | Describes local curvature and may be inaccurate near boundaries or for non-quadratic shapes |
Common Pitfalls and Practical Limits
- Treating a success flag as proof of convergence. An algorithm can stop because the step size is small, the iteration limit is reached, or a numerical tolerance is met. Review gradients, objective changes, constraints, and sensitivity to settings.
- Assuming a positive-definite Hessian is required for every strict minimum. It is sufficient at a twice-differentiable interior point, but not necessary. A singular Hessian requires further analysis; a positive-semidefinite Hessian is necessary at a twice-differentiable interior local minimum.
- Confusing a stable mode with a stable posterior prediction. A mode summarizes one point. A prediction may depend on posterior spread, skewness, or multiple regions and must be checked by integration.
- Changing coordinates without accounting for the density measure. A transformed density requires the change-of-variables Jacobian. Keep the intended posterior target unchanged and state the coordinates used for a mode.
- Calling repeated agreement a model validation. Independent numerical runs check computational reproducibility, not whether the structural model, data, prior, or residual assumptions are appropriate.
Numerical diagnostics cannot eliminate all uncertainty. Sparse patient data may leave parameters weakly informed, and a well-behaved computation can still yield a broad or prior-dominated posterior. Conversely, a visibly unstable computation should not be treated as a precise individual estimate. If reasonable reparameterizations, tolerances, or independent runs produce material differences, report the instability and investigate before using the result for a clinical forecast.
Check Your Understanding
Answer each question, then compare with the solutions below.
- An optimizer reports a tiny final step but a large scaled gradient. What does this suggest, and what should be checked next?
- At a twice-differentiable interior minimum, must the Hessian be positive definite? State the necessary condition and the role of positive definiteness.
- In the Gaussian example above, what is the posterior variance of \(\eta\), and why does evaluating the forecast only at the posterior mean underestimate its posterior mean?
- Two independent runs have nearly equal objectives but noticeably different parameter estimates. Give one plausible explanation and one useful follow-up check.
- Why is agreement among repeated optimization runs not evidence that the PK model is scientifically correct?
Solutions. (1) A small step alone does not establish convergence; the large scaled gradient suggests the objective may still decrease. Check scaling, derivatives, objective changes, constraints, and sensitivity to starting values and tolerances. (2) No. The Hessian must be positive semidefinite at a twice-differentiable interior local minimum. A positive-definite Hessian is sufficient for a strict local minimum, but a strict minimum can have a singular Hessian, so further analysis is needed in that case. (3) The posterior variance is \(0.125\). Since \(E(e^\eta)=e^{m+v/2}\), the positive variance term \(v/2\) means \(E(e^\eta)>e^m\) when \(v>0\); substituting only the posterior mean omits this contribution. (4) A weakly identified ridge or multiple posterior modes is plausible. Compare scaled gradients and objective values, inspect local curvature, and compare relevant predictions across additional dispersed starts. (5) Runs assess whether the computation reproduces a result under the tested settings. They do not test whether the model structure, data, prior, or residual assumptions describe the clinical process adequately.