Why Stiffness Can Stop an ODE Fit
An ODE likelihood calculation can fail even when the model equations and parameter values are scientifically reasonable. One cause is stiffness: the system contains substantially different dynamical time scales, and a numerical solver must take very small steps to remain stable while following the fastest dynamics. That step-size restriction can make integration slow or cause it to fail before reaching an observation time.
The central diagnostic question is whether the integration failed at a particular parameter vector, or whether integration succeeded but the optimizer did not meet its stopping criteria. These are different problems. An optimizer cannot reliably compare likelihoods when some candidate likelihood evaluations fail, and an integration warning alone does not establish that the optimizer has failed.
For a system \(y'=f(t,y)\), the Jacobian matrix \(J=\partial f/\partial y\) describes the local response of the rates to changes in the states. Its eigenvalues can help identify fast and slow modes. A large separation between the magnitudes of relevant negative real parts is a warning sign: one mode may decay quickly while another changes over a much longer interval. Eigenvalues are evidence, not a complete diagnosis; transient behavior, state scales, tolerances, and the solver also matter.
A Stability Restriction in a Simple System
A short derivation shows why an explicit solver can struggle even when the rapid component has nearly disappeared. For the scalar test equation \(y'=\lambda y\), with \(\lambda<0\), the explicit Euler update is \(y_{k+1}=y_k+h\lambda y_k=(1+h\lambda)y_k\). Numerical stability requires the magnitude of the multiplier not to exceed one.
Thus, a fast-decaying mode with a large \(|\lambda|\) can force a small step \(h\), even if the quantity of interest varies slowly at that point. In practice, accuracy requirements can impose a still smaller step, and higher-order methods have their own stability regions. The Euler calculation is a justification for the mechanism, not a recommendation to use Euler for fitting.
Worked Example: A Fast and a Slow Mode
Consider the hypothetical system \(y_1'=-1000y_1\), \(y_2'=-y_2\), with \(y_1(0)=y_2(0)=1\). Its exact solutions are \(y_1(t)=e^{-1000t}\) and \(y_2(t)=e^{-t}\). The first mode has a time scale of \(1/1000=0.001\) time units; the second has a time scale of \(1\) time unit.
For the fast mode, explicit Euler requires \(h\leq 2/1000=0.002\) time units for linear stability. At \(h=0.001\), its update multiplier is \(1-1000(0.001)=0\). Starting from \(y_1(0)=1\), Euler therefore gives \(y_1(0.001)=0\), whereas the exact value is \(e^{-1}\approx 0.3679\). Stability alone does not guarantee adequate accuracy.
At \(t=0.01\), the exact states are \(y_1=e^{-10}\approx 0.00004540\) and \(y_2=e^{-0.01}\approx 0.9900\). The fast mode has nearly vanished, but an explicit method may still be constrained by its stability requirement while integrating the system. This contrast between a short fast time scale and a long slow one is the stiffness signal.
Reading Solver Diagnostics and State Trajectories
A solver’s internal diagnostics provide evidence about how it handled a proposed parameter vector. The most informative quantities are not just the final status message, but also the number of accepted and rejected steps, the smallest attempted step, the number of right-hand-side evaluations, and the time or state at which progress stopped. Interpret these relative to the integration interval, tolerances, and solver settings.
- Repeated step rejection or unusually many evaluations can indicate that the solver is working hard to control local error or maintain stability.
- Steps shrinking to a minimum permitted size, followed by failure to advance, are consistent with stiffness, though they can also arise from other numerical problems.
- Failure at a consistent time or state across repeated evaluations at the same parameters suggests an integration-level issue worth investigating.
- Oscillating, negative, or rapidly diverging states may indicate instability or inadequate resolution. Whether a state may be negative depends on the model and its numerical representation, so inspect the equations and expected scale.
A solver may also report successful completion while a trajectory is inaccurate. Compare trajectories under a suitable tighter tolerance or an alternative solver, and inspect predictions at relevant observation times. The earlier material on Numerical Integration Error in ODE Likelihoods covers the role of integration accuracy in the likelihood; here the focus is on recognizing stiffness as a cause of failed or unreliable integration.
Worked Example: Interpreting a Failed Explicit Integration
Suppose a hypothetical pharmacometric model has two state components with the rates in the preceding example and is integrated from time 0 to time 1. At a fixed parameter vector, an explicit solver reports 1,860 accepted steps, 412 rejected steps, a minimum step of \(10^{-8}\) time units, and failure to reach time 1. A second solver using an implicit method reaches time 1 and reports the states at time 0.01 as \(y_1=0.0000454\) and \(y_2=0.9900\).
The exact values at time 0.01 are \(e^{-1000(0.01)}=e^{-10}\approx0.00004540\) and \(e^{-0.01}\approx0.9900\). The implicit solver’s reported values agree with these to the displayed precision. The explicit solver’s many rejections and very small minimum step, together with the fast and slow time scales, support stiffness as a plausible explanation for its failure.
This comparison does not prove that every implicit result is correct or that stiffness is the only possible cause. It does isolate the evidence: the parameter vector and equations were held fixed, while the integration strategy changed. Check the implicit solution at other relevant observation times and under suitable tolerance changes before relying on its likelihood.
Separating Integration Failure from Optimizer Non-Convergence
The optimizer proposes parameter vectors and uses their objective values and, often, derivatives to choose subsequent proposals. If integration fails at a proposed vector, the objective may be unavailable or unreliable there. In contrast, optimizer non-convergence means that objective evaluations were sufficiently available for the algorithm to proceed, but its convergence requirements were not satisfied before termination. A run can exhibit both problems, so inspect the evaluation-level record rather than relying only on the final fit label.
Worked Example: Classifying an Optimizer Log
In a hypothetical fit, the optimizer evaluates the objective successfully at each recorded iteration. The objective decreases from 120.4 to 108.7 over 80 iterations, then stops with the message “maximum iterations reached.” The reported gradient norm is 0.42, while the configured convergence threshold is 0.01. No integration failures are recorded for those evaluations.
The relative objective decrease is \((120.4-108.7)/120.4=11.7/120.4\approx0.0972\), or 9.72%. The objective improved, but that decrease does not establish convergence. The gradient norm is \(0.42/0.01=42\) times the stated threshold, and the optimizer explicitly stopped at its iteration limit. On these facts, this is optimizer non-convergence, not evidence of failed integration.
A useful follow-up is to rerun from the reported estimates with an appropriate iteration limit and inspect whether the objective and parameters stabilize while integration continues to succeed. If instead a fixed estimate repeatedly produces solver failure, that is a separate integration diagnosis; increasing the optimizer’s iteration limit would not repair it.
A Practical Diagnostic Sequence
A controlled sequence helps avoid changing several possible causes at once. Record the parameter vector, initial conditions, integration interval, solver, tolerances, and failure time so that a comparison can be reproduced.
Determine whether the ODE solver failed at a particular evaluation, returned a suspicious trajectory, or completed successfully while the optimizer stopped without meeting its criteria.
Evaluate the same parameter vector without allowing the optimizer to change it. Record the failure time, step statistics, and state values so the integration behavior can be separated from parameter updates.
Where appropriate, compare an explicit method with a stiff-capable implicit method and check whether the resulting trajectories agree at observation times. A solver change is a diagnostic test, not a substitute for validating the solution.
Use tolerances suitable for the state magnitudes and likelihood sensitivity. If a modest tolerance change causes large prediction changes, the numerical solution is not yet trustworthy.
Once objective evaluations are repeatable and trajectories are checked, assess optimizer stopping criteria, gradients, and parameter movement on their own terms.
The main solver alternatives have different trade-offs. Explicit methods are often efficient for nonstiff systems, but their stability restrictions can make them inefficient on stiff systems. Implicit methods are designed to handle a wider range of stiff stability requirements, but each step can require more computation, such as solving nonlinear equations. Neither class guarantees accuracy if tolerances, model equations, or state scaling are inappropriate.
| Evidence or choice | What it supports | What it does not establish |
|---|---|---|
| Many explicit-solver rejections and very small steps | A stability or error-control difficulty; stiffness is one possible cause | That the model equations are correct or that stiffness is the only cause |
| Successful implicit solve with matching checked trajectories | Evidence that the explicit failure was solver-related | That all predictions or likelihood calculations are accurate |
| Successful objective evaluations plus unmet gradient or stopping criteria | Optimizer non-convergence | That the integration is accurate unless it has been checked |
| Optimizer stopped after an integration error | A fit run affected by an integration failure | Optimizer non-convergence as the sole diagnosis |
Common Pitfalls and Limits
- Calling every small step “stiffness.” Small steps can also result from tight tolerances or other difficulties. Use the model’s time scales, step statistics, and a controlled solver comparison together.
- Increasing the iteration limit after a solver failure. This only gives the optimizer more opportunities to encounter problematic evaluations. First establish that the ODE can be evaluated reliably at the relevant parameter vectors.
- Switching solvers and accepting the result automatically. An implicit solver can complete a calculation that an explicit solver cannot, but its trajectory still needs accuracy checks at relevant times.
- Confusing stability with accuracy. A stable numerical trajectory can still be inaccurate, as the Euler calculation showed. Check both successful completion and agreement under appropriate accuracy checks.
Stiffness is not diagnosed by the optimizer’s final status alone, and a solver failure does not by itself identify its cause. A defensible diagnosis combines the equations’ time scales, solver-level evidence, fixed-parameter tests, and checked state trajectories. Keep model parameters fixed while testing solver behavior; changing parameters at the same time makes it harder to tell whether the numerical method or the model evaluation changed.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions.
- For \(y'=-200y\), what is the explicit Euler stability limit on the step size \(h\)?
- A solver repeatedly rejects steps and stops at its minimum step at the same fixed parameter vector. Name one plausible diagnosis and one additional check that would strengthen it.
- An optimizer completes every objective evaluation but reaches its iteration limit with a gradient norm above its threshold. Is this primarily evidence of integration failure or optimizer non-convergence?
- Why is a successful implicit solve not, by itself, proof that the likelihood is accurate?
Solutions. 1. The Euler condition is \(h\leq2/|\lambda|\), so \(h\leq2/200=0.01\) time units. 2. Stiffness is plausible; compare a stiff-capable implicit solver at the same parameter vector and check trajectory agreement at observation times. 3. It is evidence of optimizer non-convergence because objective evaluations succeeded but the stated convergence criterion was unmet. 4. Successful completion establishes that the solver reached the endpoint, not that errors are small enough for the likelihood; check predictions under suitable tolerances and, where appropriate, an alternative solver.