Calibration & adaptive steps¶
A probabilistic ODE solver returns a Gaussian posterior over the solution: a mean trajectory and an uncertainty band. This notebook is about the quality of that band, which is set by two things.
- Diffusion calibration controls the size of the band. The filter is correct in the mean even with an arbitrary scale, but the band is only correctly sized once a diffusion parameter $\sigma^2$ is estimated from the ODE residuals the filter sees as it integrates.
- Adaptive step-size control decides where effort is spent: small steps through fast transitions, long strides across the smooth plateaus in between.
We walk through every calibration mode in turn on a simple logistic ODE --
the uncalibrated band ("none"), post-hoc MLE, online "dynamic", and
per-component "diagonal"/"diagonal_ekf0" -- explaining what each does and
when to reach for it. Then we meet a multi-scale problem that genuinely needs
per-component calibration, and finish with adaptive stepping and a consistency
check on the band.
Acronyms used below: ODE (ordinary differential equation), IVP (initial value problem), EKF (extended Kalman filter), MLE (maximum-likelihood estimate), IWP (integrated Wiener process prior).
For the prose companions to this notebook, see Diffusion Calibration and Adaptive Step-Size Control.
Setup¶
We import the public API and fix a consistent plotting style. Throughout we use a small colour palette keyed by calibration mode, so the same mode always looks the same across every figure.
import jax
import jax.numpy as np
import matplotlib.pyplot as plt
import numpy as onp
from ode_filters import (
IWP,
ODEInformation,
gaussian_filter,
gaussian_filter_adaptive,
posthoc_mle_sigma_sqr,
quasi_mle_sigma_sqr,
rescale_sqr_seq,
taylor_mode_initialization,
)
jax.config.update("jax_enable_x64", True)
plt.style.use("ode_filters.mplstyle")
# Standard figure sizes.
FIG_WIDE = (9.0, 3.4) # single-axis time plots
FIG_STACK2 = (9.0, 5.4) # two stacked axes
# One colour per calibration mode, so legends line up across figures.
MODE = {
"analytic": dict(color="black", ls="-", alpha=0.55),
"none": dict(color="grey"), # uncalibrated
"post-hoc": dict(color="#1f77b4"), # blue
"dynamic": dict(color="#d62728"), # red
"diagonal": dict(color="#2ca02c"), # green
"diagonal_ekf0": dict(color="#ff7f0e"), # orange
}
Why calibrate? The uncalibrated band¶
We start on the logistic ODE,
$$\dot x = x(1-x), \qquad x(0)=0.1, \qquad t \in [0, 8],$$
a smooth, well-behaved problem with a closed-form solution $x(t) = 1 / (1 + (1/x_0 - 1)e^{-t})$ to check against.
The filter propagates a covariance whose overall magnitude is governed by a
single scalar diffusion $\sigma^2$. With calibration="none" we leave
$\sigma = 1$, the raw scale of the IWP$(2)$ prior. Let's see what the band looks
like before any calibration: solve with calibration="none" and plot the
filtered mean and its $2\sigma$ band against the truth.
def vf_log(x, *, t):
return x * (1 - x)
def logistic_truth(t, x0):
return 1.0 / (1.0 + (1.0 / x0 - 1.0) * onp.exp(-t))
x0_log = np.array([0.1])
tspan_log = (0.0, 8.0)
prior_log = IWP(q=2, d=1)
mu0_log, S0_log = taylor_mode_initialization(vf_log, x0_log, q=2)
measure_log = ODEInformation(vf_log, prior_log.E0, prior_log.E1)
N_log = 80
# Uncalibrated (sigma = 1) fixed-step filter. The mean is already correct; the
# band is just the prior's raw scale. The residuals it records (result.mz /
# result.Pz_sqr) are what the calibration estimators will consume.
res_log = gaussian_filter(
mu0_log, S0_log, prior_log, measure_log, tspan_log, N_log, calibration="none"
)
ts_log = onp.asarray(res_log.t)
x_log = onp.asarray(res_log.m[:, 0])
P_sqr_log = res_log.P_sqr
P_log = np.einsum("nij,nik->njk", P_sqr_log, P_sqr_log)
std_uncal = onp.sqrt(onp.asarray(P_log[:, 0, 0]))
x_exact = logistic_truth(ts_log, float(x0_log[0]))
fig, ax = plt.subplots(figsize=FIG_WIDE)
ax.fill_between(
ts_log,
x_log - 2 * std_uncal,
x_log + 2 * std_uncal,
color=MODE["none"]["color"],
alpha=0.3,
label=r"uncalibrated $2\sigma$ band",
)
ax.plot(ts_log, x_log, color=MODE["dynamic"]["color"], lw=1.4, label="filtered mean")
ax.plot(ts_log, x_exact, **MODE["analytic"], label="analytic truth")
ax.set_xlabel("t")
ax.set_ylabel("x(t)")
ax.set_title(r"Logistic ODE with calibration='none': correct mean, meaningless band")
ax.legend(loc="lower right")
plt.tight_layout()
plt.show()
The mean sits exactly on the analytic solution -- the filter is right -- but
the $2\sigma$ band is absurdly wide, dwarfing the solution itself. That width is
nothing but the IWP$(2)$ prior's arbitrary raw scale; it carries no information
about the actual error. To get a sensibly-sized band we must estimate the diffusion
$\sigma^2$ from the ODE residuals the filter recorded, available as
result.mz and result.Pz_sqr. The next sections do exactly that, one
calibration mode at a time.
Option 1 -- Post-hoc calibration¶
The simplest fix is to fit one global scalar $\widehat\sigma^2$ after the solve and rescale every covariance by it. Calibration asks: what value of $\sigma^2$ would have made the residuals the filter actually saw look like draws from a well-calibrated Gaussian? On a smooth problem there are two equivalent answers.
- Per-step quasi-MLE. Each step produces a predicted-residual marginal $(m_z, S)$. The per-step estimate is the whitened residual size $\widehat\sigma^2_n = m_z^{\top} S^{-1} m_z / d$.
- Post-hoc joint MLE. Under a single constant $\sigma^2$ explaining every step, the joint MLE has a closed form -- the mean of the per-step estimates, $\widehat\sigma^2 = \tfrac1N \sum_n \widehat\sigma^2_n$.
We feed the residual marginals from the uncalibrated run above into both estimators.
mz_log = res_log.mz
Pz_sqr_log = res_log.Pz_sqr
# Per-step quasi-MLE (one number per step) and post-hoc joint MLE.
sigma_sqr_per_step = onp.asarray(
[float(quasi_mle_sigma_sqr(mz_log[i], Pz_sqr_log[i])) for i in range(N_log)]
)
sigma_sqr_post = float(posthoc_mle_sigma_sqr(mz_log, Pz_sqr_log))
print(f"per-step quasi-MLE mean : {sigma_sqr_per_step.mean():.6g}")
print(f"post-hoc joint MLE : {sigma_sqr_post:.6g}")
print(f"absolute difference : {abs(sigma_sqr_per_step.mean() - sigma_sqr_post):.2e}")
per-step quasi-MLE mean : 0.000273259 post-hoc joint MLE : 0.000273259 absolute difference : 1.08e-19
The two numbers agree to machine precision -- exactly as the closed form predicts. Now we apply the calibration: rescale every stored covariance by $\widehat\sigma^2$ and compare the band before and after. (Recover a covariance from its square root via $P = P_{\mathrm{sqr}}^{\top} P_{\mathrm{sqr}}$.)
P_sqr_log_cal = rescale_sqr_seq(P_sqr_log, sigma_sqr_post)
P_log_cal = np.einsum("nij,nik->njk", P_sqr_log_cal, P_sqr_log_cal)
std_cal = onp.sqrt(onp.asarray(P_log_cal[:, 0, 0]))
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=FIG_STACK2, sharex=True)
ax1.plot(ts_log, x_log, color="black", lw=1.2, label="filtered mean")
ax1.fill_between(
ts_log,
x_log - 2 * std_cal,
x_log + 2 * std_cal,
color=MODE["post-hoc"]["color"],
alpha=0.4,
label=r"calibrated $2\sigma$ band",
)
ax1.set_ylabel("x(t)")
ax1.set_ylim(0.0, 1.1)
ax1.set_title("Logistic ODE: mean and post-hoc calibrated uncertainty")
ax1.legend(loc="lower right")
ax2.semilogy(ts_log, std_uncal, color=MODE["none"]["color"], lw=1.2, label=r"uncalibrated ($\sigma=1$)")
ax2.semilogy(
ts_log, std_cal, color=MODE["post-hoc"]["color"], lw=1.4, label="post-hoc calibrated"
)
ax2.set_xlabel("t")
ax2.set_ylabel(r"$\sigma_x(t)$ (log)")
ax2.set_title(
rf"Calibration shrinks the band uniformly by "
rf"$\sqrt{{\widehat\sigma^2}} \approx {sigma_sqr_post**0.5:.2g}$"
)
ax2.legend(loc="lower right")
plt.tight_layout()
plt.show()
Post-hoc calibration multiplies every posterior standard deviation by the same $\sqrt{\widehat\sigma^2}$, bringing the band down to a sensible width. Here that is a large shrinkage, because the IWP$(2)$ prior is loose for such a smooth problem. The shape of the $\sigma_x(t)$ curve is unchanged -- a single scalar is applied uniformly.
Takeaway. This is the cheapest mode: one number for the whole trajectory. You run uncalibrated, then rescale. It is ideal when you only want an honest global error bar after the fact.
Option 2 -- Online dynamic calibration (the default)¶
Post-hoc calibration waits until the end and applies one global scalar. The
online alternative -- calibration="dynamic", which is the library
default -- estimates a per-step $\widehat\sigma^2_n$ at each step and bakes it
into that step's process noise before propagating. The band then carries each
step's local scale forward as the filter integrates, rather than being
rescaled by a single global number at the end.
On a smooth single-component problem the per-step estimates are nearly constant,
so "dynamic" essentially matches post-hoc -- which is why you rarely call
post-hoc by hand. It is the right default for everyday use; you get a calibrated
band in one pass, no second step required. The interesting question is what
happens when the per-step estimates are not nearly constant. That is the next
section.
Option 3 -- Per-component calibration for multi-scale problems¶
A single scalar $\widehat\sigma^2$ assumes every component of the system needs the same amount of slack. When components live on wildly different scales that assumption breaks. The fix is per-component calibration: a separate $\widehat\sigma^2_i$ for each dimension.
Consider two logistic components with staggered transitions that live five orders of magnitude apart:
$$\dot x_i = r\, x_i (1 - x_i), \quad r = 2, \quad x_1(0) = 10^{-5}, \quad x_2(0) = 10^{-10}, \quad t \in [0, 15].$$
Component $x_1$ switches on around $t \approx 5.8$ and $x_2$ around $t \approx 11.5$: two well-separated fast regions, with the components on vastly different scales until each fires. The closed-form solution $x_i(t) = 1 / (1 + (1/x_i(0) - 1)e^{-rt})$ gives us a reference.
We solve with the public adaptive solver gaussian_filter_adaptive under
three calibration modes:
"dynamic"-- one scalar $\widehat\sigma^2$ shared across all components (the default);"diagonal"-- a separate $\widehat\sigma^2_i$ per component;"diagonal_ekf0"-- the same per-component idea with a different, exactly-diagonal denominator.
gaussian_filter_adaptive adapts the step internally and returns the solution
at the times in save_at. We compare the final values against the analytic
solution.
R_RATE = 2.0
X0_VEC = onp.array([1e-5, 1e-10])
def vf_dl(x, *, t):
return np.array([R_RATE * x[0] * (1 - x[0]), R_RATE * x[1] * (1 - x[1])])
def analytic(t, x0, r):
return 1.0 / (1.0 + (1.0 / x0 - 1.0) * onp.exp(-r * t))
tspan_dl = (0.0, 15.0)
prior_dl = IWP(q=3, d=2)
mu0_dl, S0_dl = taylor_mode_initialization(vf_dl, np.asarray(X0_VEC), q=3)
measure_dl = ODEInformation(vf_dl, prior_dl.E0, prior_dl.E1)
save_at = np.linspace(tspan_dl[0], tspan_dl[1], 200)
runs_pub = {
mode: gaussian_filter_adaptive(
mu0_dl, S0_dl, prior_dl, measure_dl, save_at, atol=1e-5, rtol=1e-3, calibration=mode
)
for mode in ("dynamic", "diagonal", "diagonal_ekf0")
}
x_true_end = onp.array(
[analytic(tspan_dl[1], float(X0_VEC[0]), R_RATE), analytic(tspan_dl[1], float(X0_VEC[1]), R_RATE)]
)
print(f"{'mode':14s} {'x1 error':>12s} {'x2 error':>12s}")
print("-" * 40)
for mode, r in runs_pub.items():
x_final = onp.asarray(r.m[-1]) @ onp.asarray(prior_dl.E0.T)
err = onp.abs(x_final - x_true_end)
print(f"{mode:14s} {err[0]:12.2e} {err[1]:12.2e}")
mode x1 error x2 error ---------------------------------------- dynamic 1.26e-13 5.72e-03 diagonal 1.51e-12 1.10e-07 diagonal_ekf0 1.51e-12 1.10e-07
All three modes nail $x_1$. On $x_2$ the scalar "dynamic" mode is markedly less
accurate -- its final-time $x_2$ error is about $6\times10^{-3}$, some five orders
of magnitude worse than the per-component modes ($\sim\!10^{-7}$). As the plot
below shows, "dynamic" switches $x_2$ on late: it lags badly through the
transition and only catches up near the end. Let us see why.
ts_dense = onp.linspace(tspan_dl[0], tspan_dl[1], 600)
x_true_0 = analytic(ts_dense, float(X0_VEC[0]), R_RATE)
x_true_1 = analytic(ts_dense, float(X0_VEC[1]), R_RATE)
def public_components(r):
"""(times, x_1(t), x_2(t)) sampled at save_at from a FilterResult."""
ts = onp.asarray(r.t)
x = onp.asarray(r.m) @ onp.asarray(prior_dl.E0.T)
return ts, x[:, 0], x[:, 1]
ts_dyn, x1_dyn, x2_dyn = public_components(runs_pub["dynamic"])
ts_dia, x1_dia, x2_dia = public_components(runs_pub["diagonal"])
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=FIG_STACK2, sharex=True)
for ax, comp, x_true, x_dyn_c, x_dia_c in [
(ax1, "x_1", x_true_0, x1_dyn, x1_dia),
(ax2, "x_2", x_true_1, x2_dyn, x2_dia),
]:
ax.plot(ts_dense, x_true, **MODE["analytic"], label="analytic")
ax.plot(ts_dyn, x_dyn_c, color=MODE["dynamic"]["color"], lw=1.4, label="dynamic (scalar)")
ax.plot(ts_dia, x_dia_c, color=MODE["diagonal"]["color"], ls="--", lw=1.4, label="diagonal (per-component)")
ax.set_ylabel(f"${comp}(t)$")
ax.set_ylim(-0.1, 1.2)
ax2.set_xlabel("t")
ax1.set_title("Multi-scale staggered logistic: scalar dynamic lags $x_2$, diagonal resolves it")
ax1.legend(loc="center right")
plt.tight_layout()
plt.show()
On $x_1$ both modes agree with the analytic solution. On $x_2$ they part ways: dynamic lags the transition badly -- near the midpoint ($t\approx11.5$) it has barely begun to rise ($\approx0.13$) while the truth is already at $\approx0.5$ -- and only catches up close to the end. Diagonal tracks the transition on time.
The mechanism is the heart of the multi-scale story. The scalar
$\widehat\sigma^2$ is a single residual size averaged across components. Once
$x_1$ is well-tracked its residual is tiny, and it drags the shared scalar down
with it. The process noise -- scaled by that same $\widehat\sigma^2$ -- shrinks
too, starving $x_2$ of the slack it needs to switch on as sharply as it
should, so its (much later) transition is delayed. A per-component
$\widehat\sigma^2_i$ gives each component its own noise budget, so $x_2$ is no
longer held back by $x_1$'s small residual. Both runs used identical tolerances
and the same controller; only calibration= changed.
We can see this in the per-component diffusion trace. In "diagonal" mode the
diffusion estimate has one entry per component; we expect $\widehat\sigma^2_1$
to peak during $x_1$'s transition and $\widehat\sigma^2_2$ during $x_2$'s --
each component flagging its own hard region. To read this trace we use the
lower-level adaptive driver introduced in the next section; for now, here is the
result.
from ode_filters.filters.ode_filter_adaptive import sqr_adaptive_loop
run_dia = sqr_adaptive_loop(
mu0_dl, S0_dl, prior_dl, measure_dl, tspan_dl,
calibration="diagonal", atol=1e-5, rtol=1e-3, h_min=1e-9,
)
sigma_dia = onp.stack([onp.asarray(s) for s in run_dia.sigma_sqr_seq])
ts_dia_steps = onp.asarray(run_dia.t_seq)[1:]
fig, ax = plt.subplots(figsize=FIG_WIDE)
ax.semilogy(
ts_dia_steps, sigma_dia[:, 0], color=MODE["diagonal"]["color"], marker="o", ms=4, lw=1.0,
label=r"$\widehat\sigma_1^2$",
)
ax.semilogy(
ts_dia_steps, sigma_dia[:, 1], color=MODE["diagonal_ekf0"]["color"], marker="s", ms=4, lw=1.0,
label=r"$\widehat\sigma_2^2$",
)
ax.axvspan(5.0, 7.0, color="lightgrey", alpha=0.25, label="$x_1$ transition")
ax.axvspan(10.5, 12.5, color="lightblue", alpha=0.25, label="$x_2$ transition")
ax.set_xlabel("t")
ax.set_ylabel(r"$\widehat\sigma_i^2$")
ax.set_title("Per-component diffusion estimates peak at their own transitions")
ax.legend(loc="lower right", ncol=2)
plt.tight_layout()
plt.show()
Each component's $\widehat\sigma^2_i$ spikes during its own transition and sits orders of magnitude lower in between. That is exactly the per-component slack the scalar mode cannot provide.
What about "diagonal_ekf0"? It is the same per-component idea -- one
$\widehat\sigma^2_i$ per dimension -- and produces nearly identical traces. The
only difference is the denominator used to whiten the residual: diagonal_ekf0
uses an exactly-diagonal form, a clean alternative when you want the
per-component covariance to be diagonal by construction.
Takeaway. Reach for "diagonal" (or "diagonal_ekf0") when components live
on genuinely different scales and a shared scalar would starve the smaller ones.
Choosing a calibration mode¶
| Mode | granularity | use when |
|---|---|---|
"none" |
-- | diagnostics only; calibrate afterwards with posthoc_mle_sigma_sqr / rescale_sqr_seq |
"dynamic" (default) |
one scalar | single-component or single-scale ODEs -- the right default |
"diagonal" |
per-component (EKF Jacobian denominator) | components on genuinely different scales |
"diagonal_ekf0" |
per-component (exactly-diagonal denominator) | clean alternative to "diagonal" |
A guideline: scalar "dynamic" is the right default and is often better on
single-scale systems, where per-component estimation only adds noise. Reach for
the diagonal modes precisely when there is a genuine gap between component
scales, as in the staggered problem above. Use "none" only when you want the
raw uncalibrated residuals and intend to calibrate post-hoc.
Adaptive step-size control¶
So far the multi-scale runs already used gaussian_filter_adaptive, the public
adaptive solver -- it chooses step sizes from the same residuals, taking small
steps through transitions and long strides across plateaus. But it reports only
the solution at save_at, not the per-step internals.
To see the controller at work we need those internals: the accepted step sizes
$h(t)$ and the reject count. For that we drop to the lower-level driver
sqr_adaptive_loop (the same one we used for the diffusion trace above),
which exposes the full accepted-step trajectory (h_seq, n_rejected,
t_seq, ...). The public entry point remains gaussian_filter_adaptive; we
drop down only to read out numbers it does not surface.
runs = {
mode: sqr_adaptive_loop(
mu0_dl,
S0_dl,
prior_dl,
measure_dl,
tspan_dl,
calibration=mode,
atol=1e-5,
rtol=1e-3,
h_min=1e-9,
)
for mode in ("dynamic", "diagonal", "diagonal_ekf0")
}
print(f"{'mode':14s} {'accepted':>8s} {'rejected':>8s} {'h range':>22s}")
print("-" * 56)
for mode, r in runs.items():
h_lo, h_hi = min(r.h_seq), max(r.h_seq)
print(f"{mode:14s} {len(r.h_seq):8d} {int(r.n_rejected):8d} [{h_lo:.3g}, {h_hi:.3g}]")
mode accepted rejected h range -------------------------------------------------------- dynamic 120 6 [0.0614, 0.75] diagonal 49 5 [0.15, 0.75] diagonal_ekf0 54 5 [0.0155, 0.75]
Now the step-size trajectory for the recommended "diagonal" run. Small steps
should cluster at the two transitions, with long strides on the plateaus.
r_rec = runs["diagonal"]
t_step = onp.asarray(r_rec.t_seq)[1:]
h_seq = onp.asarray(r_rec.h_seq)
fig, ax = plt.subplots(figsize=FIG_WIDE)
ax.semilogy(t_step, h_seq, marker="o", color=MODE["diagonal"]["color"], ms=3.5, lw=0.6)
ax.set_xlabel("t")
ax.set_ylabel("$h$ (log)")
ax.set_title(
f"Adaptive step size (diagonal): {len(h_seq)} accepted, "
rf"{int(r_rec.n_rejected)} rejected, $h \in$ [{h_seq.min():.2g}, {h_seq.max():.2g}]"
)
plt.tight_layout()
plt.show()
The controller takes a few short steps early while it calibrates, then settles into long strides, with the reject count staying in single digits. The step sizes respond to the solution's difficulty, not to noise in the diffusion estimate.
Is the band consistent? The whitened-residual check¶
Calibration scales the band to match the residuals -- but how do you verify the result? Look at the whitened residuals. For each step the standardized innovation is $z_n = S_n^{-1/2} m_z^{(n)}$, and on a well-specified problem the squared, dimension-normalized version $\lVert z_n \rVert^2 / d$ should hover around $1$: consistently $\gg 1$ means the band is too tight, $\ll 1$ means too loose, and $\approx 1$ means calibrated. We compute it directly from the per-step residual marginals on the accepted steps.
def whitened_norm_sq(mz, Pz_sqr):
"""||S^{-1/2} m_z||^2 / d for one step, with S = Pz_sqr.T @ Pz_sqr."""
mz = onp.asarray(mz)
Pz_sqr = onp.asarray(Pz_sqr)
z = onp.linalg.solve(Pz_sqr.T, mz) # S^{-1/2} m_z via the square root
return float(z @ z) / mz.shape[0]
for mode in ("dynamic", "diagonal"):
r = runs[mode]
vals = onp.array(
[whitened_norm_sq(mz, Pz) for mz, Pz in zip(r.mz_seq, r.Pz_seq_sqr, strict=True)]
)
print(f"{mode:10s}: mean ||z||^2/d = {vals.mean():.3g} median = {onp.median(vals):.3g}")
dynamic : mean ||z||^2/d = 0.238 median = 0.137 diagonal : mean ||z||^2/d = 0.301 median = 0.172
Both modes keep the mean whitened residual below $1$ (dynamic $\approx 0.24$, diagonal $\approx 0.30$; the medians lower still), so by the rule above their recorded bands are, if anything, a touch loose -- but crucially neither blows up. That is the subtlety: the whitened residual cannot, on its own, reveal the $x_2$ failure. Dynamic mode still looks self-consistent because its scalar $\widehat\sigma^2$ matches the residuals it actually sees, which are dominated by $x_1$; the $x_2$ mistracking stays hidden under the aggregate. A per-component failure can hide inside a scalar summary, so always cross-check against a reference (here, the analytic solution and the per-component diffusion trace) when you can.
Summary¶
We walked through every calibration mode in turn:
"none"leaves $\sigma = 1$: the mean is correct but the band is the prior's arbitrary raw scale -- meaningless until calibrated.- Post-hoc (
posthoc_mle_sigma_sqr+rescale_sqr_seq) fits one global scalar after the solve. Cheapest mode; run uncalibrated, then rescale. "dynamic"(the default) bakes a per-step estimate into the process noise online. On smooth single-scale problems it matches post-hoc, in one pass."diagonal"/"diagonal_ekf0"estimate a separate $\widehat\sigma^2_i$ per component. On multi-scale problems a shared scalar gets dragged down by the well-tracked component and starves the others; per-component calibration gives each its own noise budget and resolves the smaller-scale dynamics.- Adaptive stepping (
gaussian_filter_adaptive) chooses step sizes from the same residuals; the lower-levelsqr_adaptive_loopexposes the per-step diagnostics ($h(t)$, reject counts, per-component $\widehat\sigma^2_i$). - Always sanity-check with the whitened residuals -- but remember a scalar summary can hide a per-component failure, so compare against a reference where you can.
If your ODE carries an exact invariant (energy, mass, a sum-to-one constraint),
you can attach it as a Conservation law to keep the posterior on the
constraint manifold; see the Advanced features notebook for that.
See also: Diffusion Calibration and Adaptive Step-Size Control.