FOCE / FOCEI

Maturity: stable — see Feature Maturity for what this means.

First-Order Conditional Estimation (FOCE) and its interaction variant (FOCEI) are the primary estimation methods in ferx-core. They use a two-level nested optimization to find maximum likelihood estimates of population parameters.

Algorithm Overview

Outer Loop (Population Parameters)

The outer loop optimizes the population parameters \(\theta\), \(\Omega\), and \(\sigma\) by minimizing the population objective function value (OFV):

\[ \text{OFV} = -2 \log L = \sum_{i=1}^{N} \text{OFV}_i \]

Parameters are internally transformed for unconstrained optimization: - Theta: Log-transformed (ensures positivity) - Omega: Cholesky factorization (ensures positive-definiteness) - Sigma: Log-transformed (ensures positivity)

Inner Loop (Empirical Bayes Estimates)

For each subject, the inner loop finds the empirical Bayes estimate (EBE) of the random effects \(\hat{\eta}_i\) by minimizing the individual negative log-likelihood:

\[ -\log p(\eta_i | y_i, \theta, \Omega, \sigma) = \frac{1}{2} \left[ \eta_i^T \Omega^{-1} \eta_i + \log|\Omega| + \sum_j \left( \frac{(y_{ij} - f_{ij})^2}{V_{ij}} + \log V_{ij} \right) \right] \]

The inner loop uses BFGS, falling back to Nelder-Mead simplex if BFGS fails.

Gradient route (analytic vs FD)

The BFGS gradient is the exact analytic η-gradient computed from hand-rolled forward sensitivities (the Dual2 type over the closed-form PK solutions) when the model is in scope, and central finite differences (FD) otherwise.

The route is resolved per subject, so individual subjects fall back to FD for features outside that scope.

The analytic route is the default (gradient_method = auto); one provider evaluation replaces FD’s ~2·n_eta+1 predictions per inner step, so it is exact and faster as the parameter count grows. Use gradient_method = fd only when you want to force finite differences.

What follows is the detail behind that one-line answer:

Reading the resolved route

The startup banner reports the route actually resolved across the population — not just the requested setting — so a silent analytic→FD fallback is visible:

  gradient: analytic (Dual2)  [requested: auto]

When the population splits across routes, the banner shows per-route subject counts, e.g. analytic (Dual2) ×118, FD ×3.

This gradient: line appears for gradient-driven estimators (FOCE/FOCEI/GN) and for imp, which reuses the EBE Hessian built via the same route. SAEM is sampling-based and reports its E-step kernel on a sampler: line instead — see SAEM.

The outer-loop route is reported as a warning

The gradient: banner and the gradient_method_inner field describe the inner (per-subject EBE) loop. The outer (population parameter) loop has its own scope, and its reported method — gradient_method_outer in the fit output — is deliberately a model-level statement: it answers “may this model use the exact analytic outer gradient”, not “did every subject get one”.

The two can disagree. A subject the provider declines — a rate-defined infusion under F ≠ 1, an out-of-scope ODE dose event, two modeled infusion windows that happen to coincide at the current estimates, or any subject whose inner Hessian is not positive-definite at the trial point — is salvaged individually, while the model as a whole stays in scope and keeps reporting analytic (Dual2). Since #1154 that mismatch is stated rather than left to be inferred from the run time:

1 of 32 subjects could not be given the exact analytic outer gradient during this fit
(e.g. subject 7) (a data shape outside the sensitivity provider's scope, or a trial point
where it could not be formed) and used fixed-EBE outer gradients, which omit the
EBE-response term the analytic gradient carries. If the fit stalls,
`reconverge_gradient_interval = N` restores the reconverged gradient every N-th
evaluation. The reported outer gradient method is the model-level route, not the
per-subject one.

What a declined subject is charged

The salvage follows the reconverge_gradient_interval schedule like every other fixed-EBE gradient (#1529): a declined subject gets the gradient at its held EBE — 2·n_free evaluations of that one subject, no inner re-solve — and an evaluation the schedule reconverges re-solves every subject’s EBE anyway. Before #1529 a declined subject was always charged a reconverged finite-difference gradient (2·n_free full EBE solves) whatever the setting. Because declines depend on the trial point, a blown-up line-search trial could send most of the population down that path at once: on a 55-subject parent–metabolite ODE model, 46 subjects declined at one such trial, about 1,200 EBE solves for a single gradient and 82% of the fit’s wall time.

The held-EBE gradient differentiates the same per-subject objective the fit reports. Where a closed form would miss a term of it — a θ-dependent residual magnitude, block_sigma correlations, M3, IIV on residual error — it is a central difference of that objective instead. A declined subject that is repelled at the trial point (its objective is the 1e20 sentinel) contributes a zero gradient: the population objective there carries the sentinel, so no line search accepts the point and its gradient is never used to take a step. If the held-EBE gradient comes back non-finite at a usable objective, that subject alone falls back to the reconverged gradient.

One kind of decline is exempt from the schedule: a structural one, which the subject takes at every parameter point rather than at an occasional trial (#1536). Today that means a block_sigma correlation across observation rows, such as total and unbound assays paired at one time. The analytic assembly serves only a diagonal residual covariance, and the pairing depends on the data alone, so such a subject is never analytic. The held-EBE gradient would then be the only gradient it ever contributes, and its missing EBE-response term biases every step rather than a few. On a 31-subject fluconazole model with paired total/unbound rows, L-BFGS stalled at OFV 807.9 on the held-EBE gradient against 738.0 with the reconverged one (NONMEM: 734.6). Structural declines therefore take the reconverged gradient whatever reconverge_gradient_interval says. The one exception is the blown-up-trial guard (#1520): under SLSQP or MMA, a trial that guard skips gets a zero contribution from a structural decline too, at a point the optimizer rejects anyway. The warning names structural declines separately: used reconverged finite-difference outer gradients because their block_sigma residuals are correlated across observation rows.

IOV models (kappa/block_kappa) always reconverge and ignore reconverge_gradient_interval, so a declined IOV subject still takes the reconverged gradient, and the warning says so: used reconverged finite-difference outer gradients; their results are correct but slower.

How declines are counted

The count is recorded as the fit runs, from the gradient evaluations that actually happened — not predicted by probing the provider at some chosen parameter point. That distinction is load-bearing, because a decline is not purely a property of the subject’s records: the separability check on modeled infusion windows reads the resolved durations and lag times, so a subject can decline at η = 0 and be served at its EBE. A probe would announce a fallback that never occurred.

Recording also means the warning is self-gating. A derivative-free BOBYQA fit (including a [mixture] model left on optimizer = auto, which resolves to BOBYQA), a Gauss-Newton or trust-region fit, a reconverge_gradient_interval = 1 fit that forces the reconverged gradient on every evaluation, and an outer_maxiter = 0 evaluation-only run never select the analytic branch for the optimizer, so they have nothing to record and say nothing. The one outer gradient most of them do evaluate — the flat-θ pre-flight check at the start estimates — keeps its own record, which is discarded. Conversely, any recorded decline is by construction a mismatch with the model-level label — including when every subject declines, which is the case that label gets most wrong.

Declines at a non-positive-definite inner Hessian

The provider is not the only thing that can decline. The outer assembly inverts each subject’s exact inner Hessian \(H = \partial^2 l_i/\partial\eta^2\) to form \(\mathrm{d}\hat\eta/\mathrm{d}\zeta = -H^{-1}M\), and it does so with a Cholesky — so a subject whose \(H\) is not positive definite declines there too, and takes the same salvage as a provider decline. This is common away from the mode: \(H\) carries the residual \(\times\) prediction-curvature term that the Gauss-Newton \(\tilde H\) drops, and 4 of the 10 subjects on the bundled examples/warfarin.ferx fail that Cholesky at \(\eta = 0\) while all 10 are positive definite at their EBEs.

That gate looks one degree too strict, because \(-H^{-1}M\) needs \(H\) only nonsingular. It is not, and the distinction is worth stating because it is easy to get backwards. At an unconstrained local minimum a nonsingular Hessian must be positive definite. So a failed Cholesky on a nonsingular \(H\) is a proof that the \(\hat\eta\) handed to the assembly is not the minimizing EBE — and \(-H^{-1}M\) is the derivative of the inner stationarity condition \(\partial l/\partial\eta|_{\hat\eta} = 0\), which therefore does not hold either. Inverting anyway would follow a saddle or non-stationary branch rather than the profiled minimum the outer objective is defined by. The Cholesky is not a conservative matrix-shape check; it is the precondition detector, and it is the only one on this path.

This was measured rather than assumed (#1513). Substituting a full-pivot LU makes a 100-subject transit-absorption fit 12% faster and eliminates all 16 of its declines, but over those 16 events the resulting gradient disagrees with the reconverged finite-difference one by a median relative \(0.378\) — against a \(1.4 \times 10^{-3}\) noise floor measured on the same fit’s positive-definite subjects — and points in a near-orthogonal direction (cosine \(0.032\)) at the worst of them. The \(H\) matrices are well conditioned throughout, so that is the broken precondition, not a numerical artifact. The fit still lands at the same optimum, which is exactly why an OFV comparison cannot settle a question like this one.

The practical reading for a user: a gradient: banner of analytic (Dual2) together with a per-subject fallback warning means some subjects’ inner loops are finishing away from their optimum, not that ferx is being needlessly cautious. Those subjects get the held-EBE gradient (the reconverged one under IOV) at those points.

Why an inner loop finishes away from its optimum

“Finishing away from the optimum” has a specific, fixable cause on closed-form absorption models, found by instrumenting every one of the 19,924 inner solves in the fit above (#1519).

A transit or inverse-Gaussian closed form is valid only while absorption outruns elimination (\(k_e < K_{tr}\) for transit, \(k_e < 1/(2\,\mathrm{MAT}\,\mathrm{CV}^2)\) for IG). Outside that domain — the flip-flop regime — ferx reroutes the subject to the model’s ODE twin, which is valid in both regimes. The reroute is a function of \((\theta, \eta)\), not of the subject’s records: a fit can start in-domain and closed-form, and walk across the abscissa as \(\theta\) moves. From that point on, that subject’s individual objective is an adaptive RK45 integration, and so carries the solver’s local-error noise floor.

The inner loop already has the right response to a noisy objective — the objective-stall stop, which accepts a point whose objective has flattened and whose gradient has stopped improving, instead of chasing a gnorm < inner_tol target that solver noise makes unreachable. Until #1519 that stop was keyed on the subject-static reroutes only (a time-varying covariate, a TIME reference, IOV), so it stayed off for every flip-flop-rerouted subject. The measured consequence, on examples/one_cpt_transit.ferx + data/datsim_oral.csv:

inner solves exit on gnorm < tol analytic vs. FD \(\partial l/\partial\eta\)
closed form, in domain (1,529) 1,529 (100%) \(6.9 \times 10^{-9}\)
flip-flop, on the ODE twin (18,395) 9,626 (52%) \(1.2\!-\!4.0 \times 10^{-4}\)

Every one of the fit’s 8,769 non-converged inner exits — 4,957 that exhausted inner_maxiter, 3,812 whose line search found no representable decrease — was a flip-flop-rerouted solve, and none of the in-domain closed-form solves failed. Each failure then bought a Nelder-Mead recovery of up to \(5 \times\) inner_maxiter iterations, which is where the wall clock went: 142.6 s before, 3.5 s after, OFV 1215.9771 → 1215.9777. The recovery is not free in accuracy either — on the worst single-subject case measured, it returned an \(\hat\eta\) whose gradient norm was \(2.4 \times 10^{-1}\) where the stalled BFGS point read \(1.1 \times 10^{-4}\).

What this does not do is make the declines go away: 16 decline events before, 15 after. A subject whose objective is an RK45 solve at default tolerances cannot resolve its gradient below roughly \(3 \times 10^{-4}\), so its exact \(H\) can still be indefinite at the point it stops. The gradient those subjects are served agrees with the reconverged finite-difference reference at the same \(1.2 \times 10^{-3}\) median relative error as the fit’s own positive-definite subjects, measured the way #1513’s table was — so the change buys time, not accuracy, and does not spend accuracy to get it.

The salvage is not bought at a blown-up trial point

A reconverged salvage — every declined subject before #1529, and still every declined IOV subject — costs 2·n_free warm EBE re-solves for one subject, and much of that spend used to land on trial points the optimizer was about to reject anyway — line-search steps whose objective sits thousands to \(10^{15}\) units above the incumbent (#1520). Since #1520 a declined subject at such a point contributes nothing to the outer gradient there instead of buying the salvage, and the warning above says how many salvages were skipped.

Two gates, both required, both in the same units — excess of a −2LL over its value at the reference point, per observation:

  • the population is worse than the reference by more than 12 per observation, and
  • the subject’s own contribution is worse than its reference value by more than 12 per observation of its own.

The line is measured, not chosen. Every salvage of every gradient-based fit of the bundled examples was instrumented (118 decline events at points worse than the incumbent, default L-BFGS and SLSQP, 13 example–optimizer pairs): at ordinary rejected trials the subject excess tops out at 6.0 and the population excess at 4.5 per observation; at blown-up ones they start at 21.8 and 36.4 and run to \(10^{19}\). Twelve sits in that gap, the widest in either distribution. It is an excess and not a ratio because a −2LL difference is unit-free and a −2LL is not — the incumbent objective is negative on 9 of the 13 pairs — and per observation because a rich subject legitimately moves by tens of units between ordinary trials.

The gates are only a rejected-trial guarantee when the reference is the point the optimizer’s own acceptance test compares against, and that point is not in general the best evaluation seen: a rejected trial can undercut the current iterate without passing sufficient decrease, and SLSQP’s inexact line search accepts its eleventh trial without sufficient decrease at all (slsqp.c, L200: line > 10), then requests the gradient there as the new iterate. Measured against the best evaluation, a guard could drop a subject’s term at such an accepted point and feed the hole into the optimizer’s Hessian update (PR #1525 review). So the reference is chosen per optimizer, from its acceptance test, or the guard is off:

outer optimizer gradient at a rejected trial point guard
NLopt SLSQP requested only at the first trial of a line search (mode = −2) and discarded when that trial fails the L1-merit test; later trials request no gradient; the accepted iterate is re-evaluated with its gradient (mode = −1) at the same point on — fires only at a fresh point (one whose x differs from the previous evaluation’s, i.e. a first trial) and measures against the last point at which a gradient was requested, which is then the current iterate; a first trial worse than the iterate fails the merit test and is rejected. The line > 10 acceptance is a re-evaluation at a non-fresh point, where the guard is off by construction
NLopt MMA requested at every inner candidate, copied into the model only when fcur < minf on — measures against the best evaluation, which with bounds only is exactly MMA’s acceptance reference
NLopt L-BFGS (Luksan) requested at every trial and read: pnint1 interpolates the next step length from the rejected trial’s directional derivative; acceptance is relative to the current iterate, which ferx cannot observe, and a trial that passes sufficient decrease but fails the curvature condition is rejected yet can be the best evaluation seen off — no reference ferx can form is a rejected-trial guarantee, so the salvage is always bought
built-in BFGS / L-BFGS never requested: the backtracking line search evaluates the objective alone; the gradient runs only at the accepted point, measured against the previous accepted iterate unreachable at a rejected trial by construction
Gauss-Newton, trust region, BOBYQA never reach the per-subject assembly, or never form a gradient not applicable

Measured on the bundled examples: under SLSQP, every fit on which the guard fires (10 examples, 1 to 13 skipped salvages each) keeps every estimate, standard error, OFV, iteration count and sdtab byte-identical to the unguarded run, at up to 4× less wall time where the salvage dominated (mm_multistart, 18.2 s → 4.5 s). Under NLopt L-BFGS, the default for analytic models, nothing changes. The bundled reproducer examples/one_cpt_transit.ferx + data/datsim_oral.csv is unchanged to the byte under either optimizer: after #1519 its 15 remaining declines all sit at one trial whose excess over the incumbent is 33 units on 1,200 observations — nothing there is blown up.

Analytic scope: models and readouts

The analytic route serves the analytical 1-/2-/3-cpt PK models and — since #410 — in-scope user-[odes] models (via a lighter Dual1 forward walk over the ODE solver), including multi-endpoint per-CMT Form-C readouts (y[CMT=N] = …, since #439) and Form-C readouts that reference a θ/η directly (e.g. y = central/V1 * (1 + ETA_CL) + TVBASE, since #486 — each bare θ/η in the readout is desugared into a hidden individual parameter so its ∂y/∂θ/∂y/∂η rides the same individual-parameter sensitivity chain; a readout referencing a neural-network output still falls back to FD).

A structural parameter that reads the event-time built-in TIME/time (a NONMEM $PK IF (TIME.GE.45) CL=… time-dependent switch) is analytic since #486 too: the per-event time is threaded into the same event-driven Dual2 (outer) / Dual1 (inner EBE) walk used for time-varying covariates, so its ∂f/∂θ/∂f/∂η are exact rather than finite-differenced — covering closed-form (non-IOV and IOV — each occasion’s stacked [η_bsv, κ] seeding is evaluated at the event time) and ODE (non-IOV and IOV) models alike. TIME composes with the other event-driven-walk features too: a direct pk(...=TIME) mapping (desugared into a synthetic individual parameter, #486), a built-in absorption input-rate forcing (the walk carries R_in alongside the per-event time, #486), and a non-zero ODE init(...) baseline (the walk seeds the init state, #486) — none of these remains a TIME-specific fallback.

Separately, an [odes] RHS that reads the TAD built-in while the model carries an estimated lagtime is analytic on both loops and under IOV since #1070. TAD’s anchor is the dose’s lagged arrival, so it moves with the lagtime: it is threaded through the sensitivity walks as a dual carrying ∂TAD/∂ALAG = −1, rather than lifted into the RHS as a constant. While it was lifted, that term was missing from the analytic chain everywhere the trajectory is integrated — not merely at a dose event. Predictions were unaffected — only the gradient (the closed-form route had a value defect for the same built-in, fixed separately in #1124; such models now always integrate) — but FOCEI builds its h matrix from that gradient, so the error reached the reported OFV: 66%–68% on the lag axis against finite differences of the production predictor, and 0.176 OFV against NONMEM where the fixed route now agrees to 0.0004 (nonmem_anchor/tad_lag_*). The anchor also carries its jet before the first dose arrives, where both ODE predictors anchor the pre-arrival window at the subject’s first arrival (#1073) and an init(...) baseline gives that window live state. Neighbouring cells are unchanged because they never carried the term: TAFD anchors at the unlagged first dose time, a per-route absorption lag does not move the TAD anchor at all, and a literal constant lagtime carries no derivative — all three were exact before and remain so.

TAD combined with a steady-state dose is the one cell that still routes to FD — not under this rule but under the broader steady-state one stated just above, which declines any non-autonomous RHS (TIME/TAFD/TAD) because it breaks the steady-state cycle recurrence.

Analytic scope: likelihood and scaling

LTBS

Log-transform-both-sides (LTBS) on analytical models is now analytic on both loops for the plain closed-form path — the inner EBE gradient applies the same g = ln(f) jet as the outer since #665, with the covariance step reconverging the EBEs tighter (cov_inner_tol, default min(inner_tol, 1e-8)) and, for a closed-form LTBS fit, the fit’s inner loop tightened to at least 1e-6 so the standard errors of weakly-identified variances are reproducible. Both tightenings follow the analytic inner gradient rather than the model shape, so a closed-form LTBS fit that also carries [iov] takes them too (#486) — its inner loop runs the same ln(f) jet and carries the same sensitivity. LTBS combined with time-varying covariates, a TIME-built-in structural parameter, or an η-dependent ExpressionScale obs_scale now also takes the analytic inner gradient — the event-driven / static inner walks apply the same ln(f) jet last, matching the outer (#486). LTBS × IOV is analytic on both loops on closed-form models (#486): each walk applies the ln(f) jet after its own per-occasion scale quotient, reproducing production’s scale-then-log order — the outer θ/Ω/σ gradient and the inner EBE gradient alike, including combined with an ExpressionScale obs_scale. The transform is one per-observation scalar quotient, identical on every stacked axis, so the occasion κ columns carry it exactly as the BSV η columns do. (LTBS × IOV on ODE models remains finite differences on both loops — that walk applies LTBS in PK-parameter space before the η/θ/κ chain, so production’s scale-then-log order cannot be rebuilt post-walk.)

IIV on residual error

On analytical (closed-form 1-/2-/3-cpt) models, IIV-on-residual-error (iiv_on_ruv) is analytic on both loops, including combined with inter-occasion variability (the residual-variance scaling exp(2·η_ruv) and the η_ruv covariance column ride the stacked [η_bsv, κ] assembly), with a custom / time-varying residual-error magnitude (the magnitude’s direct-θ terms thread the residual-eta d/R coupling on the FOCEI outer gradient — non-IOV since #673, under IOV since #486, as the κ-augmented residual-eta assembly is dimension-generic), and with M3 BLOQ (the censored data term contributes its η_ruv first-order column, and the censored second-order cross-terms enter both the inner Hessian and the Laplace H̃ consistently with quantified rows, at FOCEI order — since #486); the non-IOV ODE iiv_on_ruv × M3 combination is analytic too (since #623).

M3 BLOQ under IOV

M3 BLOQ combined with inter-occasion variability is analytic across the full estimator matrix on closed-form models (FOCEI since #580; non-interaction FOCE and the triple M3 + IOV + iiv_on_ruv since #591): the censored −logΦ((LLOQ−f)/√v) data term and its f-derivatives ride the stacked [η_bsv, κ] assembly — entering the data gradient, the true inner Hessian, and (since #486) the Laplace H̃ at FOCEI order, consistently with quantified rows and exactly as on the non-IOV M3 path. Under FOCEI the inner EBE and the outer θ/Ω/σ gradient match reconverged FD of the censored FOCEI marginal; under FOCE the censored rows leave the augmented Sheiner–Beal marginal and re-enter as the marginal tail probability −logΦ((LLOQ−f0)/√R̃ⱼⱼ), R̃ⱼⱼ = HⱼΣ_bHⱼᵀ + R⁰ⱼ — the same linearized-marginal moments the quantified rows use (#646), matching the way Monolix’s linearization likelihood treats censored data — so plain FOCE stays a consistent Sheiner–Beal objective for the whole subject (no silent promotion to FOCEI). FOCEI (a conditional first-order method) instead keeps the censored term at the conditional f(η̂)/R⁰, so FOCE-M3 and FOCEI-M3 are genuinely different optima. FOCEI is the conditional-family method NONMEM’s METHOD=1 LAPLACE M3 belongs to (NONMEM cannot run M3 without LAPLACE) — but ferx’s FOCEI drops the ∂²f/∂η² second-order term that full Laplace keeps, so it matches NONMEM’s Laplace M3 only up to that FOCEI-vs-Laplace term (ferx has no full-Laplace M3 mode yet). Only the ODE M3 + IOV + iiv_on_ruv triple still falls back to FD.

Output scaling

An η-dependent expression obs_scale is fully analytic on both loops (since #486): on analytical models the inner EBE gradient applies the η-block of the outer’s scale quotient rule; on ODE models the same quotient is applied on the static walk and on the event-driven walk (time-varying covariates and/or the TIME built-in) — the divisor is subject-static (production evaluates it at the subject covariate snapshot), so the non-IOV event-driven walk applies one subject-static jet since #486 — and the IOV ODE walk evaluates one scale jet per occasion group so IOV + time-varying covariates (and TIME) remains analytic since #590 (combined with LTBS it still falls back to FD). A constant ScalarScale output divisor (f/k) is analytic under IOV on both engines since #486 — the trivial covariate/η/κ-independent case of the same quotient: on analytical (closed-form) models it divides the final jet uniformly, and on ODE models the in-walk readout already divides p/k over the stacked dual.

Analytic scope: dose events

On analytical (closed-form 1-/2-/3-cpt) models, modeled-duration/rate doses (RATE=-1/-2 → D{cmt}/R{cmt}) are analytic on both loops since #486, closing the largest closed-form-vs-ODE gap (the ODE path had it via #530/#635): the modeled infusion window resolves from the PK parameters, so the infusion end is a moving boundary in D/R that the event-driven walk carries exactly — to second order — via the dual window length (the sign-mirror of the dose-start lagtime handling), for non-IOV and IOV (each occasion resolves its own window from the per-occasion PK jet, including a κ-coupled modeled slot), on both the outer θ/Ω/σ gradient and the inner EBE η-gradient. Modeled-dose × steady state is analytic too (#486): the closed-form dual SS equilibration threads the modeled-window jet (rate, dur) into each cycle’s active/quiet split, so the moving infusion-end flows through the SS trough exactly as it does through the current pulse (matching the ODE path). Only a rate-defined (RATE=-1) infusion under F ≠ 1 stays on FD (as on the ODE path).

On analytical (closed-form) models, an estimated lagtime (ALAG/LAGTIME) is analytic on both loops even when the subject also carries inter-occasion variability, time-varying covariates, or a TIME-dependent structural parameter (#486). Previously any of those routed the subject to the event-driven walk, which declined a lagtime and dropped the whole fit to finite differences — so an everyday oral model with a lag and IOV never got its exact gradient. The walk now threads each dose’s arrival t + ALAG as a moving boundary carried in the dual sub-interval lengths, which is the same mechanism the modeled infusion end already used: because the closed-form flow is analytic in the interval length, the bolus saltation falls out to all orders and no explicit saltation term is injected (the ODE path must inject one only because its RK45 timeline is f64). This covers a lag that varies per occasion (a κ on the lag: each occasion’s doses arrive at their own offset) and a fixed-duration infusion under a lag (whose window end rides the arrival). Two deliberate exceptions remain on FD: a steady-state dose combined with a lagtime (production overlays the pre-arrival SS tail, which the dual walk has no twin of — serving it would disagree with production in value, not just in derivative), and the measure-zero case of an observation, reset, or covariate change landing exactly on a lagged arrival (the walk resolves one dual position per break time, so it declines rather than return a subtly wrong jet). Steady-state dosing combined with an EVID 3/4 reset is analytic on the closed-form engine too (#486): dose superposition genuinely cannot express it — an SS dose folds in an infinite periodic history that the reset truncates — but the event-driven walk can, and production already computes exactly these subjects with the f64 instance of that same walk, so the dual walk simply mirrors it.

On ODE models the event-driven dual walk (since #439) is analytic for bolus dosing — with or without time-varying covariates, the TIME built-in, and IOV — and now also covers an estimated lagtime (bare LAGTIME/ALAG or per-compartment ALAG{n}, injected as an event-time saltation), finite-duration infusions, modeled-duration/rate doses (RATE=-1/-2 → D{cmt}/R{cmt}, whose infusion window is a moving boundary in the modeled parameter carried by a rate-off saltation — analytic since #530, and combined with IOV since #486, where each occasion resolves its own window from the per-occasion PK jet, including a κ-coupled modeled slot), EVID 3/4 resets, and EVID=2 covariate-only breakpoints on both the non-IOV and IOV ODE paths (κ fixed at zero on the IOV path), for both the outer θ/Ω/σ gradient and the inner EBE η-gradient (since #449/#590/#486).

A time-varying-covariate + init(...) ODE model is analytic, with the initial state seeded at the subject’s first-record covariate snapshot (matching production, #486) — including combined with a reset, an estimated lagtime, a finite infusion, an input-rate forcing, steady state, or a modeled dose.

Analytic scope: absorption and steady state

A built-in absorption input-rate forcing (igd/transit/weibull/first_order) combined with an EVID 3/4 reset or time-varying covariates is analytic since #486: the reset case threads the tracked reset_floor into the shared forcing helper, and the TV-cov case builds each dose’s forcing constants once, from its own dose record’s PK jet (#1569). Combined with an estimated lagtime, igd/transit/first_order are also analytic since #486 — the continuous ∂R_in/∂lag flows through the walk’s dual time-after-dose, and the forcing’s onset at the dose’s lagged arrival is injected as an exact rate-on saltation. zero_order(dur) is also analytic combined with time-varying covariates or an estimated lagtime since #486: the event-driven walk now carries its moving-boundary window (delivered as a per-segment constant, with rate-on / rate-off saltations at the moving window start / end — the rate-off using the general g⁻ − g⁺ saltation so a covariate that varies across the window end is exact). Because the IOV walk is the same event-driven walk, every input-rate kind the non-IOV walk carries is analytic under IOV too (#486): the smooth densities igd/transit/weibull/first_order, the moving-window zero_order, and any composition (mixed, parallel) — each forcing’s rate/window is rebuilt from its dose’s own occasion-seeded PK jet, so κ rides through. Compartment-indexed bioavailability F{cmt} and lagtime ALAG{cmt} are analytic under IOV as well (the walk resolves each dose’s own compartment slot).

A steady-state (SS=1) dose into a built-in absorption compartment is analytic on both loops and both engines since #835: the closed-form periodic steady-state trough u_ss = (I − M)⁻¹·b is carried over the dual type (its implicit-function derivative rides the same linear solve), so first_order/transit/igd/weibull SS absorption gets exact FOCEI gradients — including under IOV, where the SS-dosed occasion’s κ enters the trough and later occasions’ κ enter the forward walk. The only input-rate holdouts under IOV are weibull + estimated lagtime (its onset diverges for shape β < 1, on every path) and a steady-state dose into a zero_order window or combined with an absorption lagtime — the zero-order window’s periodic-window equilibration and the SS+lagtime pre-arrival seed are not yet carried through the kernel (both are rejected upstream with a clear error rather than silently mis-served). Steady-state (SS=1) dosing is analytic too (#473), including combined with a modeled-duration/rate dose, a rate-defined infusion under F ≠ 1, or an estimated lagtime (each closed by #486 — see Steady-state dosing); only SS combined with a non-autonomous RHS (reads TIME/TAFD/TAD) stays on FD.

Known FD fallbacks

The route is resolved per subject and per feature; anything below stays on central finite differences (correct, just slower). Each entry is described in full in the scope sections above.

Feature / combination Engine Loops on FD Why
LTBS × IOV ODE both LTBS is applied in PK-parameter space before the η/θ/κ chain, so production’s scale-then-log order cannot be rebuilt post-walk
η-dependent obs_scale × LTBS ODE both as above
M3 BLOQ + IOV + iiv_on_ruv (the triple) ODE both the only surviving M3 × IOV gap
input_rate forcing + SDE [diffusion] both both —
weibull + estimated lagtime both both the onset can diverge for shape β < 1
RATE=-1 infusion under F ≠ 1 both both —
SS + non-autonomous RHS (reads TIME/TAFD/TAD) ODE both a non-autonomous field breaks the steady-state cycle recurrence
observation / reset / covariate change exactly on a lagged arrival closed-form both the walk resolves one dual position per break time
model wider than the size caps (θ + η > 24 ODE / > 32 closed-form event walk) both both too wide to monomorphise

A steady-state dose into a zero_order window and a steady-state dose combined with an absorption lagtime are not on this list because they are not served at all: both are rejected upstream with a clear error rather than silently falling back.

Two further combinations are described inconsistently in the scope text above and are deliberately left off the table until the code is re-checked: ODE iiv_on_ruv combined with M3 BLOQ (recorded both as an inner-loop-only FD fallback and as analytic since #623), and a closed-form steady-state dose combined with an estimated lagtime (recorded both as a deliberate FD exception and as closed by #486). The scope sections above — likelihood and scaling and dose events — carry both statements as written.

Other per-subject FD fallbacks include an input_rate forcing combined with an SDE [diffusion] block, or weibull combined with an estimated lagtime (its onset can diverge for shape β < 1) — plus, for the inner EBE gradient only, IIV-on-residual-error combined with M3 BLOQ on ODE models (the outer θ/Ω/σ gradient still serves it; plain ODE-LTBS and plain ODE iiv_on_ruv are analytic on both loops since #474, as the Dual1 ODE walk shares solve_ode_g with the objective).

Moving boundaries and one-sided derivatives

A moving boundary — a lagged dose arrival, the rate-on of a lagged infusion, an infusion or zero-order window end — is not a data record, so it supplies no parameters of its own: it subdivides the interval that the next record terminates, and every piece of that interval runs on that record’s snapshot (#1073). A boundary falling strictly inside a record interval therefore has the same covariate snapshot on both sides, and the covariate field jump that used to be injected there is simply absent. Any forcing that straddles the instant — another dose’s infusion, a zero-order window, a built-in R_in — is still carried in both velocities, so the rate and mass jumps stay exact.

A record landing exactly on the boundary

A record landing exactly on the boundary is the one place the derivative really is two-sided: the value stays continuous, but the field assignment flips as the boundary sweeps through the record, so the two one-sided derivatives differ by a finite amount and a central finite difference returns their average, which is neither. The walk returns a one-sided derivative rather than declining, and — since #1068 — it returns the branch its own event ordering implements, because both sides of the boundary read one snapshot. Which branch that is follows from where the boundary sorts against a co-timed record: a dose arrival, a lagged infusion’s rate-on, a built-in absorption onset and a per-route lag= onset all sort before the record, so they report the limit from below; an infusion or zero-order window end sorts after it, so it reports the limit from above. The rule is uniform across the four onset sites by construction rather than by coincidence — a first_order(ka=KA, lag=LAG) forcing and an otherwise identical unlagged one take different timeline events, and reporting opposite branches for them would be a difference of bookkeeping rather than of model. Before #1068 the arrival arms mixed two records at that one instant — the pre side from the co-timed record, the post side from the next one — which attributed a stationary field jump to a moving boundary and returned a value outside both one-sided limits (measured 0.89 % on ∂f/∂η_lag, first order, on a multi-dose subject; invisible on a single-dose one, where the state at the arrival is zero and the spurious term is multiplied away). A second moving boundary colliding with the arrival (#1076) and an input-rate onset’s field jump (#1069) are exact for the same reason the interior case is, and were re-measured at 4.6e-9 and 5.9e-9 respectively.

The remaining approximation: a non-autonomous RHS

The one approximation left in this cell is a straddling built-in R_in whose kernel reads the changing covariate, leaving an explicit ½(∂f⁻/∂t − ∂f⁺/∂t) curvature term uncarried (#1075, tracked in the #1078 epic) — and it is not confined to any co-timed geometry: ẍ = ∂f/∂t + J·f and every saltation site carries only J·f, so any non-autonomous field leaves the explicit term uncarried — including an [odes] RHS that simply reads TAD under an η-carrying lagtime, with no straddling input rate and no colliding boundary. Measured on that geometry at 2.9e-1 on ∂²f/∂η_ALAG² (#1070 improved it from 6.4e-1 but does not close it) — against two independent oracles, double-FD of the production predictor and a central difference of the analytic first-order gradient, which agree with each other to seven significant figures. The straddling-R_in geometry #1075 was originally filed on, at 6.4e-2, is exact today (7.4e-7); the non-autonomous RHS is what is left.

Be precise about what that second-order residual reaches, because it is wider than “the Hessian”: ∂f/∂η and ∂f/∂θ are exact, but the FOCE/FOCEI objective’s θ-gradient is assembled from second derivatives of f — d2f_deta2 feeds the inner Hessian and the ∂log|H̃|/∂η curvature term in sens_outer_gradient.rs, so it enters the analytic outer gradient handed to L-BFGS/SLSQP, as well as the covariance step and standard errors. It does not enter the objective value, which is why the nonmem_anchor/tad_lag_* evaluations agree to 4e-4. On a TAD-reading RHS with an η on the lagtime, prefer optimizer = bobyqa (derivative-free) until #1075 lands if you want the pre-#1070 numerical behaviour; gradient_method = auto now resolves such models to L-BFGS, where it previously chose BOBYQA.

Sampling exactly on a moving dose boundary

With a modeled RATE=-1/-2 dose the infusion window end moves with the estimated D{cmt}/R{cmt}, and an observation may sit exactly on it (sampling at end-of-infusion is a normal design). The prediction is kinked in D at that point: just above the coincidence the infusion is still running at the sample, so only the rate amt/D matters; just below it the dose has finished and a decay term appears. The two one-sided slopes therefore differ, and no two-sided derivative exists — a central finite difference returns their average, which corresponds to neither branch. ferx returns the one-sided derivative — specifically, the derivative of the branch its own event ordering already uses to define the value at that instant (an infusion at its end is still contributing; a dose at its arrival has landed). Both engines agree: the ODE path injects the standard jump/saltation term at the event, and since #486 the closed-form path reaches the identical value by stepping the read-out state back across the zero-length window between the sample and the boundary. The same applies to a sample landing on a lagged dose arrival.

The practical consequence is that at exactly such a coincidence an analytic gradient and a finite-difference gradient will not agree: a central difference averages across a branch switch the model does not actually make there (and across a bolus arrival it would return an enormous value, since it straddles a jump). That is a property of the model, not a defect, and it is transient in a fit — the lagtime and D are estimated, so they only sweep past a fixed sample time momentarily.

Model-size and readout scope limits

Two scope limits were also widened (#486). A [scaling] y = <expr> Form-C readout whose individual parameters spill past the eight structural PK slots (e.g. a sigmoid-Emax readout on a 3-compartment oral model, which exhausts them immediately) is now analytic on the time-varying-covariate walk rather than FD. And the model-size caps that decide when a model is too wide to monomorphise the exact gradient were raised: θ + η from 16 to 24 on the ODE and output-scaling paths, and from 24 to 32 on the closed-form event walk. The old limit was tighter than it looked — a mid-sized covariate model (5 structural θ + 6 covariate-effect θ + 5 η) sat at exactly 16, so adding one more covariate silently dropped the entire fit to finite differences. Beyond the raised caps the fallback is still FD (correct, just slower).

Non-finite individual objectives

If the EBE search wanders into a region where the individual NLL evaluates to a non-finite value (for example, an ODE model whose integration blows up at extreme \(\eta\)), that point is treated as the worst possible objective rather than aborting the fit. The subject is reported as non-converged and estimation continues for the remaining subjects.

Runaway EBEs are re-checked (#958)

A subject’s inner search can be led somewhere its own objective does not actually have a minimum, and then stop there: once the steps are noise rather than descent the objective stops improving and the gradient norm plateaus, which is exactly what the convergence test looks for. The usual cause is a finite-difference inner gradient on an objective whose value carries the ODE solver’s local error. That error is normally far below the signal being differenced, but a prediction driven toward zero under a proportional error model has its residual variance clamped to the floor, and the resulting (y − f)² / R term amplifies solver noise by many orders of magnitude — enough for the FD gradient to come back with the wrong sign. No FD step size repairs that: too small and it is noise, too large and it truncates on a very curved objective.

So whenever the returned η̂ sits further than about ten prior standard deviations per random effect from the prior mean (measured as ηᵀΩ⁻¹η, the prior’s own metric), the inner loop re-solves derivative-free from the prior mean and keeps whichever point has the lower objective. The check can only improve the objective — a genuinely extreme EBE that a derivative-free search cannot beat is the objective’s own minimum and is returned unchanged — so fits whose inner loop was already right are unaffected.

The practical guidance is unchanged: leave gradient_method at its default auto. The analytic route differentiates the same integration it predicts from, so it has no differencing step to be swamped and is immune to this failure mode; gradient_method = fd on a model with near-zero predictions is accuracy-limited by ode_abstol, not by the step size.

Gradient accuracy vs cost (reconverged EBEs)

By default the population gradient holds each subject’s EBEs (\(\hat\eta\)) and FOCE Hessian fixed while perturbing the population parameters — the fixed-EBE gradient. This is cheap, but it omits the response of the inner solution to \(\theta\) and \(\Omega\). On well-conditioned problems that omitted term is negligible and a gradient optimizer matches the derivative-free bobyqa (see Outer Optimizers). On ill-conditioned problems it is not: the omitted term is what separates a true descent direction from a flat-looking plateau, and slsqp can report converged at an OFV far above the bobyqa optimum — which is exactly why the default optimizer = auto falls back to bobyqa on the FD-only fits where this bias bites.

Set reconverge_gradient_interval = 1 (in [fit_options]) to re-solve the inner EBE loop at every finite-difference perturbation, recovering the full surface. This costs roughly 5–6× per gradient — reserve it for fits whose OFV looks suspiciously high. IOV models (kappa/block_kappa) always reconverge, so the setting is a no-op there; it only changes non-IOV fits.

To amortize that cost, reconverge_gradient_interval = N reconverges only every N-th gradient evaluation and uses the cheap fixed-EBE gradient in between — the periodic correction is often enough to keep the optimizer off the plateau at a fraction of the always-on cost. The default 0 disables reconverging entirely (cheap fixed-EBE gradient).

Validation — cefepime pediatric population PK, 2-compartment IV infusion, combined error, no IOV (5937 subjects / 17380 observations). All ferx rows use FOCEI and the same likelihood convention, so OFVs are directly comparable:

Configuration OFV Wall time
slsqp (fixed-EBE gradient) 68,252 390 s
slsqp + reconverge every 10 (interval = 10) 66,118 633 s
slsqp + reconverge every 5 (interval = 5) 66,056 1,004 s
slsqp + reconverge every eval (interval = 1) 65,485 1,871 s
bobyqa (derivative-free) 65,598 315 s
ferx OFV evaluated at NONMEM’s final estimates 67,514 —

The fixed-EBE gradient stalls slsqp ~2,650 OFV units above bobyqa. Reconverging the EBEs on every gradient evaluation closes the entire gap and reaches an optimum marginally below bobyqa’s — and below the point NONMEM converged to (NONMEM’s estimates score 67,514 under ferx’s likelihood), confirming the stall was a gradient-bias artifact, not a worse model.

A larger reconverge_gradient_interval trades that accuracy back for speed: reconverging every 5th or 10th gradient still closes ~80% of the stall, but the OFV degrades monotonically as the interval grows (the cheap fixed-EBE gradients in between bias the direction enough that slsqp declares convergence slightly early). On this problem every interval setting is dominated by bobyqa on both OFV and wall time — which is why the default optimizer = auto selects bobyqa on FD-only fits like this one. The reconverged-slsqp path earns its keep when a gradient optimizer is required (e.g. parameter count high enough that derivative-free search scales poorly, or downstream tooling that consumes the optimizer’s gradient).

FOCE vs FOCEI

Standard FOCE

Uses linearized predictions around the EBEs:

\[ f_0 = f(\hat{\eta}) - H \hat{\eta} \]

where \(H\) is the Jacobian matrix \(\partial f / \partial \eta\). The per-subject objective is:

\[ \text{OFV}_i = (y - f_0)^T \tilde{R}^{-1} (y - f_0) + \log|\tilde{R}| \]

where \(\tilde{R} = H \Omega H^T + R(f(\eta{=}0))\). The residual variance \(R\) is evaluated at the population prediction \(f(\eta{=}0)\) — matching NONMEM’s METHOD=1 (no INTER) semantics — not at the linearized point \(f_0\). On a nonlinear model (e.g. oral absorption) the linearization \(f_0 = f(\hat\eta) - H\hat\eta\) can extrapolate to a near-zero or negative concentration where the true typical-individual prediction is healthy; evaluating a proportional variance \((f\sigma)^2\) there collapses the weight and makes the objective multimodal with an indefinite covariance. \(f(\eta{=}0)\) is always physically sensible, so FOCE+proportional converges deterministically and matches NONMEM FOCE. Additive error is unaffected (\(R\) does not depend on \(f\)).

How \(\tilde{R}^{-1}\) is formed

\(\tilde{R}\) is \(n_{obs} \times n_{obs}\) but its random-effect part \(H\Omega H^T\) has rank \(n_{eta}\), so the inverse is normally taken in the low-rank (Woodbury) form, which factorizes an \(n_{eta} \times n_{eta}\) matrix instead. That identity is subtractive, and it loses significant digits in proportion to how far \(R(f(\eta{=}0))\) sits below \(H\Omega H^T\). The two are evaluated at different points — \(R\) at the typical individual, \(H\) at the subject’s own \(\hat\eta\) — so a subject whose typical prediction has decayed onto the residual-variance floor (\(10^{-12}\); a proportional error model at a steady-state trough, for instance) while its own prediction has not can drive that ratio past \(10^{13}\), at which point the low-rank inverse has no accurate digits left and the analytic FOCE gradient it feeds is wrong rather than merely imprecise.

ferx therefore bounds the conditioning before inverting, from a quantity the low-rank form already computes (\(1 + \mathrm{tr}(\Omega H^T D^{-1} H)\), \(D = \mathrm{diag}(R)\)), and falls back to the observation-sized factorization of \(\tilde{R}\) itself — stable in exactly that regime — whenever the bound is exceeded. The choice is per subject and per evaluation, affects only floating-point accuracy and not the model, and well-conditioned subjects keep the fast path. There is no user-facing setting: the switch is automatic.

The low-rank form is written in its \(\Omega\)-normalized spelling (\(\Omega = LL^T\), factorize \(I + L^T H^T D^{-1} H L\)) so that the matrix actually factorized is the one the bound covers; factorizing \(\Omega^{-1} + H^T D^{-1} H\) instead would inherit the conditioning of \(\Omega\), which a near-rail block_omega correlation degrades independently of the bound. Without the check a fit reaching that regime stalls far from the optimum and reports converged = false (#1498).

FOCEI (Interaction)

Uses individual predictions directly without linearization. ferx-core implements the Laplace approximation form from Almquist et al. (2015):

\[ \text{OFV}_i = (y - \hat{f})^T V^{-1} (y - \hat{f}) + \hat{\eta}^T \Omega^{-1} \hat{\eta} + \log|\tilde{R}| \]

FOCEI is more accurate when the residual variance depends on the predicted value (proportional or combined error models), because it captures the interaction between random effects and residual error.

Almquist, J., et al. (2015). Comparison of maximum a posteriori and conditional likelihood estimation in the FOCE method. PAGE 24, Abstr 3516.

Optimizer Options

The outer optimizer is configured independently of the estimation method. See Outer Optimizers for a full description of all available algorithms (SLSQP, BOBYQA, trust-region, BFGS, L-BFGS, MMA) and guidance on when to use each.

Set via [fit_options]:

[fit_options]
  method    = focei
  optimizer = slsqp    # override the default (auto)

Covariance Step

When covariance = true, ferx-core computes the variance-covariance matrix of the parameter estimates (the R-matrix: the inverse observed Fisher information) at the converged solution. This provides:

  • Standard errors (SE) for all parameters
  • Relative standard errors (%RSE) for assessing estimation precision
  • Omega SEs via delta method from the Cholesky parameterization

Which stage runs it, and when it does not run at all

In a method chain the covariance step runs once, on the last stage that estimates — a trailing evaluation-only stage (imp_eval_only, agq_eval_only) reads out a likelihood at fixed parameters and does not take ownership of it. Recomputing the step after every stage would pay for it several times over for one reported result (#615).

Bayesian estimation is the one method that never runs it: method = bayes reports posterior credible intervals from the sampled chain, not Hessian standard errors, so ferx disables the covariance and SIR steps for it. That makes all six covariance keys — covariance, covariance_method, covariance_fallback, analytic_cov_hessian, fd_hessian_step, cov_inner_tol — inert when a fit ends in bayes, and setting one there is reported as such:

fit option `cov_inner_tol` configures the post-fit covariance step, but this fit
ends in Bayesian estimation, which reports posterior credible intervals instead of
Hessian standard errors and runs no covariance step, so it has no effect.

Where the Bayes stage falls in the chain is what decides this, not its presence: methods = [bayes, focei] ends on FOCEI, which runs the step and consumes every one of those keys, so no notice is issued. methods = [focei, bayes] does not.

Every other estimator — FOCE, FOCEI, Laplace/AGQ, gn, gn_hybrid, SAEM, IMP, IMPMAP, VI — runs the covariance step and honours all six keys. They are framework-level for that reason, rather than listed per method (#956).

How the Hessian is built

The covariance is 2·H⁻¹, where H is the Hessian of the objective. There are two routes to H, and both produce the same object — the observed information at the converged point.

Sensitivity assembly

The analytic R-matrix (analytic_cov_hessian = true, the default since #436) assembles H directly from third-order prediction sensitivities, differencing nothing at the population-objective level. The observed information needs exactly one derivative more than the outer gradient: Almquist et al. (2015) give the first derivative of the FOCE/FOCEI marginal, and the second is third-order in f. ferx obtains that third order by central-differencing the analytic second-order Dual2 jet along each (θ, η) axis — 2·(n_theta + n_eta) + 1 sensitivity evaluations per subject, all at the converged point, with no inner re-solve. For closed forms that jet is machine-precision and the step is ε^(1/3); for [odes] it carries solver error and the step is derived from the effective ode_reltol (capped at 1% of parameter scale). In both cases this first difference amplifies the jet’s numerical error as noise/h, rather than the reconverged objective stencil’s noise/h²; fd_hessian_step does not enter.

For IOV, the random-effect axes include every occasion’s κ, while each shared IOV covariance entry remains one population parameter. Diagonal and correlated IOV blocks are supported; the validation results include a NONMEM R-matrix comparison.

Its scope is narrow, and deliberately so: Gaussian closed-form and [odes] models whose second-order sensitivity provider is available, including IOV and M3 censoring, without LTBS, observation scaling, closed-form Form-C readouts, [initial_conditions], iiv_on_ruv, correlated or custom-magnitude residual error, FREM, a covariate-Selected error spec, or a non-Gaussian endpoint — and not under method = laplace, whose marginal needs fourth-order prediction derivatives. Dual-evaluable ODE y = ... and y[CMT=N] = ... readouts are supported because their derivatives are already part of the ODE sensitivity jet. The third-order prediction blocks are central differences of the existing Dual2 sensitivity jet; closed forms use the machine-precision step and ODE models choose the step from ode_reltol. FOCE uses the linearized-marginal M3 tail described above; FOCEI uses its conditional M3 tail. FOCEI-anchored AGQ has its own analytic assembly, including ODE, IOV, and M3. IOV with modeled lag is also excluded. The scope is applied per subject, not to the population as a whole (#1514). The observed information is the sum \(\sum_i R_i\), and every term in it is the second derivative of one subject’s own marginal, so a subject the analytic assembly declines is finite-differenced on its own — from the same objective, at the same converged point, warm-started from the same modes — while every other subject keeps its exact term. That is a different operation from averaging two approximations of one matrix: each term here estimates itself, which is how the fit-side gradient has assembled the population since #466. What licenses the decomposition is that both halves of the covariance objective are per-subject: the EBE reconvergence warm-starts every perturbed point from the fit’s modes, a constant across the stencil, so \(x \mapsto \hat\eta_i(x)\) is a deterministic per-subject map; and the population marginal is a sum over subjects. The one non-separable objective, the mixture marginal, is excluded from the analytic route outright.

Two limits on the salvage. A model-level exclusion — method = laplace, a mixture, gradient = fd, analytic_cov_hessian = false — is not per-subject and still returns the whole fit to the stencil. And when at least half the population declines, the whole-population stencil is taken instead of salvaging: the win comes from an analytic majority carrying the cost, and at parity there is none. On a 56-subject FOCEI ODE fit where all 56 declined, rebuilding the matrix as 56 per-subject stencils cost +4% wall for cells identical to 7 significant figures.

The measured effect where it does apply: on a 55-subject 2-state ODE FOCEI fit with 13 free parameters, one out-of-scope subject cost 23.4 s of covariance under the old all-or-nothing rule and 2.0 s under the salvage — 11.5× on the step, −38% total wall, with estimates and OFV identical to every printed digit. On warfarin, forcing one subject in ten to decline moves the standard errors by 8.4e-5 relative, against a 2.7e-4 gap between the whole-population FD and whole-population analytic routes the code already ships as interchangeable. A fit that salvages any subject reports a W_COV_ANALYTIC_SALVAGE note naming them. Everything outside the scope keeps the stencil below, subject to numerical convergence and step-size limits.

Moving dose events

Where an axis moves a dose event — an estimated lag, a zero_order duration, F on an infusion — the step is bounded further, so that no observation, EVID=2, reset or dose record changes sides of the event between the two points of the pair: the subject’s events are enumerated by the ODE engine’s own break-time builder at both perturbed points and the step is shrunk until the pair is clear. Before that bound, SE(TVLAG) on a lagged-dose fixture moved 27 % with the step. A subject whose mode sits exactly on a moving event — a corner minimum of the inner objective, which a lagged arrival at a dense early sample does produce — has no third derivative there and declines to the per-subject salvage above rather than differencing across the corner (#1505).

Finite-difference fallback

The finite-difference stencil is the fallback, and remains the only route for every model outside that scope. Two details make its SEs match NONMEM’s $COVARIANCE:

  1. The EBEs are reconverged at every perturbed point. Like NONMEM’s covariance step, ferx re-solves the inner conditional-estimation loop (warm-started from the converged EBEs) at each finite-difference point. Holding the EBEs fixed omits the response of the conditional optimum to the population parameters, which gives a Hessian with the wrong curvature — indefinite even on well-conditioned surfaces such as warfarin — and was the source of the spurious eigenvalue clipping in issue #129.
  2. The Hessian is a reconverged second difference of the objective (issue #209): a 3-point diagonal / 4-point off-diagonal difference of the marginal OFV, with the EBEs reconverged at every perturbed point. Because it re-evaluates the full marginal end-to-end (recomputing a = ∂f/∂η and the log|H̃| EBE-response at each point), it captures the complete curvature with no envelope approximation and no held-fixed a, and serves FOCE, FOCEI and IOV — and additive/proportional/combined error — uniformly. An earlier release also offered a covariance_ofv_hessian = false path that took a central difference of the analytical population gradient (2·n_free gradient evaluations instead of ~2·n_free² objective evaluations); that stencil was a-fixed and only exact to first order in the log|H̃| term (it carried a ~9% TVKA bias and could explode SEs at non-converged points), so it was removed in #639 — the reconverged-OFV difference is now the sole stencil.

Scaling and mode convergence

The factor of two is because the objective is −2·logL: its Hessian is twice the observed information, so the covariance needs 2·H⁻¹ to be on the right scale. It applies to the analytic route as well — the per-subject assembly returns ∂²Fᵢ/∂x² with OFV = 2·Σᵢ Fᵢ, so the population sum is scaled by two before inversion. (Omitting that scale is silent: the matrix stays symmetric and positive-definite, and every SE comes out inflated by exactly √2.)

Both routes reconverge the EBEs at cov_inner_tol before doing anything else. The stencil needs it at every perturbed point; the analytic assembly needs it once, at the converged point, because its envelope and EBE-response terms are derived under stationarity (∂lᵢ/∂η|_η̂ = 0) — a mode left at the fit’s looser inner_tol would contaminate the Hessian with the residual inner gradient.

The covariance objective is 2·pop_nll for both FOCE and FOCEI — the Ω penalty is already inside each per-subject marginal and must not be added again. FOCEI’s Almquist–Laplace marginal carries η̂ᵀΩ⁻¹η̂ + log|Ω| explicitly; FOCE’s Sheiner–Beal marginal carries it through R̃ = HΩHᵀ + R (equivalent by the Woodbury identity). An earlier version added the Ω prior a second time on the FOCE path, double-counting Ω and under-stating the FOCE omega SEs by ~31% — fixed in issue #243.

For a mixed block + diagonal Ω, the structural-zero cross-block off-diagonals (which are not estimated) are excluded from the covariance parameter set like FIX parameters; otherwise their flat Hessian diagonal aborts the step. (The same #243 change exposed and fixed this for both FOCE and FOCEI.)

NONMEM cross-check

On data/warfarin.csv (10 subjects, 1-cpt oral, proportional error) against NONMEM 7.5.1 $COVARIANCE MATRIX=R, FOCEI ($EST METHOD=1 INTER):

Parameter ferx SE NONMEM SE rel. diff
TVCL 0.00712 0.00710 +0.3%
TVV 0.2350 0.2401 −2.1%
TVKA 0.1506 0.1486 +1.3%
PROP_ERR (SD) 0.000832 0.000835 −0.4%
ω²(CL) 0.01291 0.01279 +0.9%
ω²(V) 0.00403 0.00431 −6.4%
ω²(KA) 0.1530 0.1504 +1.8%

and FOCE ($EST METHOD=1, no INTER), where ferx’s OFV/estimates already match NONMEM (OFV −280.17 vs −280.36) and, after the #243 covariance fix, the omega SEs do too:

Parameter ferx SE NONMEM SE rel. diff
TVCL 0.00661 0.00663 −0.3%
TVV 0.2375 0.2340 +1.5%
TVKA 0.1209 0.1245 −2.8%
PROP_ERR (SD) 0.000922 0.000941 −2.0%
ω²(CL) 0.01312 0.01280 +2.5%
ω²(V) 0.00454 0.00430 +5.7%
ω²(KA) 0.1510 0.1607 −6.1%

(Analytic-Jacobian route; the FD-Jacobian fallback agrees to within ~10% on the ω block.) Both are guarded by tests/warfarin_covariance_nonmem.rs within a 20% band.

Covariance estimator: R, S, or sandwich

By default ferx reports the R-matrix covariance R⁻¹ (the inverse observed information), which assumes the model is correctly specified. covariance_method selects an alternative, mirroring NONMEM’s $COVARIANCE MATRIX=:

covariance_method Estimator NONMEM Use when
r (default) R⁻¹ MATRIX=R Model is well-specified; you want the model-based SEs.
s S⁻¹ MATRIX=S You want the empirical-information (cross-product) SEs.
rsr R⁻¹ S R⁻¹ MATRIX=RSR You want SEs robust to model mis-specification (the Huber–White “sandwich”).

Here S = Σᵢ gᵢgᵢᵀ is the cross-product of the per-subject score vectors gᵢ = ∂(−logLᵢ)/∂θ — the same per-subject gradients the Gauss–Newton optimizer uses for its BHHH step. At the MLE of a correctly-specified model the information-matrix equality gives R ≈ S, so all three estimators converge to the same SEs; they diverge when the model is mis-specified, where rsr is the conservative choice.

[fit_options]
  method           = focei
  covariance       = true
  covariance_method = rsr

s and rsr are available for FOCEI, FOCE, and IOV fits. Note that s requires a positive-definite covariance Hessian, analytic R-matrix or FD stencil alike (it reuses r_inv; a non-PD R is not yet rerouted to S⁻¹ directly — see TODO in outer_optimizer.rs).

NONMEM cross-check (S / RSR, FOCEI)

Same warfarin FOCEI fit as the MATRIX=R table above, with two additional NONMEM covariance runs ($COVARIANCE MATRIX=S and MATRIX=RSR); SEs are the .ext ITERATION = -1000000001 row (#266):

Parameter ferx s NONMEM MATRIX=S rel. diff ferx rsr NONMEM MATRIX=RSR rel. diff
TVCL 0.01031 0.00930 +10.9% 0.00742 0.00710 +4.6%
TVV 0.4250 0.4602 −7.7% 0.2418 0.2403 +0.6%
TVKA 0.2157 0.2268 −4.9% 0.1555 0.1487 +4.5%
PROP_ERR (SD) 0.001519 0.001545 −1.7% 0.000790 0.000796 −0.7%
ω²(CL) 0.02005 0.01764 +13.7% 0.01170 0.01093 +7.1%
ω²(V) 0.007139 0.008287 −13.9% 0.003699 0.003978 −7.0%
ω²(KA) 0.2265 0.2596 −12.8% 0.1311 0.1396 −6.1%

(ci/FD build.) The s cross-product matches NONMEM within ~14% and rsr within ~7%. Guarded by covariance_se_matches_nonmem_s_rsr in tests/covariance_method_sandwich.rs (20% band for s, 15% for rsr). The per-subject score S is on the −logL (NLL) scale and, unlike R, carries no factor of two; the close agreement confirms that scaling is correct.

NONMEM cross-check (S / RSR, FOCE)

Same warfarin model, $EST METHOD=1 (no INTER), $COVARIANCE MATRIX=S and MATRIX=RSR (#250):

Parameter ferx s NONMEM MATRIX=S rel. diff ferx rsr NONMEM MATRIX=RSR rel. diff
TVCL 0.00824 0.00838 −1.6% 0.00720 0.00757 −4.9%
TVV 0.3822 0.3996 −4.4% 0.2552 0.2457 +3.8%
TVKA 0.1445 0.1495 −3.3% 0.1373 0.1328 +3.4%
PROP_ERR (SD) 0.001593 0.001573 +1.3% 0.000942 0.000947 −0.6%
ω²(CL) 0.01604 0.01771 −9.4% 0.01038 0.01092 −5.0%
ω²(V) 0.008768 0.008387 +4.5% 0.004057 0.003957 +2.5%
ω²(KA) 0.2443 0.2336 +4.5% 0.1460 0.1549 −5.7%

(ci/FD build.) All within ~10%. Guarded by covariance_se_matches_nonmem_foce_s_rsr in tests/covariance_method_sandwich.rs.

Note that ferx defaults to covariance_method = r whereas NONMEM’s $COVARIANCE default is the rsr sandwich; set covariance_method = rsr when reconciling SEs against a default NONMEM run.

Running the covariance step after a fit (run_covariance)

The covariance step can also be run as a standalone step against a FitResult produced earlier — the covariance-step analogue of run_sir. This is useful when the original fit was expensive and you want to add (or re-run, e.g. with a different covariance_method) the covariance step without re-estimating, or when working with a fit loaded from a .fitrx bundle.

use ferx_core::{fit_from_files, run_covariance, FitOptions};

let mut opts = FitOptions::default();
opts.run_covariance_step = false;               // skip cov during the fit…
let fit = fit_from_files("model.ferx", Some("data.csv"), None, Some(opts.clone()))?;

opts.covariance_method = ferx_core::CovarianceMethod::Sandwich;   // …run it now, as RSR
let fit_with_cov = run_covariance(&fit, None, None, &opts)?;

run_covariance reconstructs the fitted parameters from the FitResult, re-runs the final inner loop (cold-started, exactly as the inline covariance path in fit() does) to rebuild the covariance-step inputs, and calls the same FD-Hessian covariance step fit() runs inline — so the result matches fitting with covariance = true. For an in-memory fit from any packed-Cholesky-space optimizer (BOBYQA/SLSQP/MMA, BFGS, the trust region, or Gauss-Newton) the match is bit-for-bit: the fit carries the optimizer’s exact Cholesky factor and run_covariance reuses it. A fit reloaded from a .fitrx bundle, or produced by SAEM / importance-sampling / Bayes (which rebuild omega from the reported matrix), instead reconstructs the factor by re-decomposing fit.omega, so it matches only up to finite-difference noise (~1e-5, platform dependent). The returned FitResult is a clone of fit with covariance_matrix, se_theta / se_omega / se_sigma / se_kappa, covariance_status, and the condition-number diagnostics refreshed. run_covariance_step on the passed options is ignored — calling the function is itself the request.

A covariance step that runs but cannot produce a usable matrix (a non-positive-definite or structurally-unusable covariance Hessian) is not an error: the returned fit reports covariance_status = Failed with a diagnostic in warnings, mirroring the inline fit() behaviour. Err is reserved for input problems — a missing or hash-mismatched model/dataset, a dimension mismatch, or an IOV model supplied without its population.

The second and third arguments (Option<&CompiledModel> / Option<&Population>) behave exactly as for run_sir: Some(...) is used as-is (caller owns verification), None re-reads from fit.model_path / fit.data_path with SHA-256 integrity checks that refuse stale source. See Caller-supplied vs. re-read inputs.

Hessian regularization

Because the EBEs reconverge, the Hessian is positive-definite on well-conditioned surfaces and no regularization is needed. As a safety net for genuinely ill-conditioned or near-singular problems, ferx-core still inverts via a symmetric eigendecomposition and clips eigenvalues below max(λ_max · 1e-10, 1e-12) to that floor before reconstructing H⁻¹, adding a warning to FitResult.warnings:

Covariance step regularized: eigenvalue floor applied to FD Hessian
  (1 of 7 free-block eigenvalues clipped; min eig = -3.21e-08, floor = 4.56e-09).
Standard errors should be interpreted with care.

SEs in the regularised directions are inflated (since 1/λ_clipped is larger than the would-be true 1/λ). With the reconverging covariance step this warning should now be rare; if it fires, the surface is genuinely near-identifiable and the SEs in those directions warrant scrutiny. A Hessian whose spectrum is genuinely indefinite — every eigenvalue ≤ 0 — still fails outright with the standard Covariance step failed message.

SIR fallback for non-PD Hessians

When a model is poorly identified or converges to a saddle point, the covariance Hessian — whether R came from the analytic route or the FD stencil — can have every free-block eigenvalue ≤ 0, making Hessian inversion impossible. Setting covariance_fallback = sir triggers a SIR (Sampling Importance Resampling) run instead of leaving the covariance step as failed:

[fit_options]
  covariance_fallback = sir

ferx-core constructs a proposal covariance by taking the absolute values of the Hessian eigenvalues, inflating by 4× (heavier tails to account for the non-PD correction), and running the same SIR sampler used by sir = true. The result is reported as 95% credible intervals in the fit YAML; covariance_status is set to sir_fallback.

Use this option when: - The Hessian is non-PD at convergence due to model mis-specification or weak identifiability - You want distributional uncertainty estimates without re-fitting with a more constrained model - You need to proceed to a downstream analysis despite a non-PD Hessian

The SIR fallback inherits all sir_* settings (sir_samples, sir_df, sir_resamples, etc.) from [fit_options].

Convergence

The outer loop terminates when any of: - The gradient norm falls below outer_gtol (default 1e-6) - The maximum number of iterations (maxiter) is reached - The optimizer reports convergence (NLopt XtolReached or FtolReached)

The inner loop terminates when the gradient norm falls below inner_tol (default 1e-5, tight enough that residual EBE noise does not perturb the marginal objective the outer optimizer minimises — see fit-options.md for tuning guidance) or inner_maxiter (default 200) iterations are reached.

When the optimizer stops without converging

NLopt has a fourth way to stop: a bare Failure, which its line search returns when it cannot make progress against the curvature estimate it has built up. That says nothing about where the fit is. It comes back both at a finished optimum — the objective is already flat to ~8 significant figures, so no step improves it — and at a point the fit was still descending through.

ferx tells the two apart on the objective trace rather than on the gradient norm, which at these optima still reads O(1): the run must have a flat tail (several evaluations since the last real improvement), must have moved off its initial estimates, and a cold re-solve of the inner loop at the restored best point must reproduce the objective the optimizer saw. All three hold ⇒ the fit converged and is reported as converged (#751).

When they do not, the optimizer is restarted once from the best point it reached (#1277), on whatever is left of the fit’s maxiter budget. It usually costs nothing in the sense that matters: the restart is adopted only if it ends at a strictly lower objective, so it can never make a reported fit worse. On the two-compartment deep-compartment model in tests/fixtures/two_cpt_dcm_regularized.ferx, L-BFGS quit at evaluation 12 with the trace still falling ~2600 OFV per evaluation; the restart ran to convergence 3886 OFV units lower (3209.81 → −676.77).

Why a fresh start helps depends on which of L-BFGS’s two line-search caps fired. NLopt’s L-BFGS (Luksan’s plis) allows ten step reductions and ten extrapolations per line search, and a cap hit on an iteration where the algorithm has just restarted — the first iteration always is one — ends the run with a bare Failure, discarding every point accepted on the way. The evaluation-12 abort above is the extrapolation cap inside the first line search: the opening-step rescale had shrunk the first gradient by \(1.7 \times 10^5\), so every trial point showed a large real decrease next to a tiny predicted slope and the search extrapolated eleven times. No curvature memory exists yet at that point; what the restart does is re-arm that rescale around a gradient that is now O(1). Later in the same fixture the other cap fires — a subject switches EBE mode, the curvature pair goes degenerate, the next direction is ~\(10^{17}\) long and ten reductions cannot find a decrease — and there the restart’s fresh curvature memory is what helps.

A restarted fit says so: FitResult.warnings carries an optimizer_health entry naming the objective and evaluation count it was restarted from, and the run’s own converged says how the restart ended. It is not gated on the optimizer — the test is “stopped with a bare Failure while still descending”, which SLSQP and the derivative-free methods can also produce, and the strictly-lower rule makes a pointless restart cost one optimization and change nothing.

Three cases are deliberately left alone. A run that spent its maxiter budget stops with NLopt’s MaxEvalReached, not a Failure, and a run that stalls with its budget already spent is not restarted either: the restart runs on the remaining budget, so the two legs together are bounded by the maxiter you set the way a single run is (NLopt checks the bound between line searches, so a few evaluations of overshoot are possible either way), and a stall that cannot be restarted for lack of budget says so with its own “increase maxiter” warning (#1428). A cancelled run is not restarted. And a second stall is reported as a stall: you get converged: false and the Outer optimization did not converge warning, which is the honest signal and bounds the worst case at two optimizations.

Warm Starting

The inner loop warm-starts EBE estimation from the previous outer iteration’s EBEs. This significantly reduces computation time, especially in later iterations when parameters change slowly.

The warm start is taken from the best point seen so far, not from whatever the optimizer evaluated last (#1290). That distinction matters because the individual objective can be multimodal: a trial step the line search is about to reject may land somewhere the inner loop solves into a different, worse basin, and if those EBEs were carried forward, every later evaluation would inherit them. The starting point would then re-evaluate worse than it did a moment earlier — an objective that moves under the optimizer’s feet, which no line search can descend. Anchoring the warm start to the incumbent keeps the objective at the incumbent reproducible; evaluations that do not improve on it leave the cache untouched.

The final inner loop

When the outer optimizer stops, ferx restores the best point it saw (#59) and runs one more inner loop there — the EBEs that reach the sdtab, the covariance step and the reported objective come from that solve, not from the trajectory.

That solve runs cold, from η = 0. When its objective reproduces the best value the optimizer saw, it is what the fit reports and nothing else runs — which is the ordinary case, and it costs exactly what it always did.

When it does not reproduce that value, two more candidates are scored on the same objective: the EBEs the optimizer was minimising against, exactly as it left them, and a re-solve seeded from those. The lowest of the three is reported. Before #833 only the cold solve ran, so where it settled at a different η̂ than the trajectory had, the reported OFV came out above the best-seen value the point had just been restored for (+3.5 on a fluconazole 2-cpt binding model), and the covariance step ran at those different EBEs. The held EBEs are scored as a candidate in their own right, not merely used as a seed, because a re-solve is not guaranteed to improve on the point it starts from.

The comparison is also a diagnostic. A cold re-solve that lands materially worse — or returns no usable number at all — raises the ebe_start_dependent warning (see Warnings), which quotes both objectives. The estimates are unaffected; what depends on the starting point is the EBEs at them, and so every diagnostic built on those (IPRED, IWRES, CWRES, shrinkage) plus the covariance step. Two things produce it, and the warning names both rather than guessing: an individual objective with more than one mode at these estimates, and an inner_maxiter a cold start cannot converge within. Raising inner_maxiter, or inner_restarts for a suspected second mode, tells them apart.

The cold number is also what the convergence self-consistency check reads (#751); a fit whose “optimum” exists only under a warm start is still rejected.

Configuration Example

[fit_options]
  method     = focei
  maxiter    = 500
  covariance = true
  # optimizer omitted → uses the default (auto); set explicitly to switch