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 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. It is resolved per subject, so individual subjects fall back to FD for features outside that scope. 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.) 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). 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 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. 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). 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. 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 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 hoists the forcing’s dose-invariant constants fresh per segment as the PK snapshot changes. 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. A time-varying-covariate + init(...) ODE model is analytic for the plain-bolus subset — the initial state is seeded at the subject’s first-record covariate snapshot (matching production, #486); init(...) combined with a reset / estimated lagtime / finite infusion / input-rate forcing / steady-state / modeled-dose still falls back to FD. It 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.

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.

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.

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).

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.

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, default)	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\)).

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 (bobyqa)

Global Search

Enable gradient-free pre-search to help escape local minima:

[fit_options]
  method        = foce
  global_search = true
  global_maxeval = 2000

The pre-search uses NLopt’s CRS2-LM algorithm (Controlled Random Search) to explore the parameter space before handing off to the local optimizer. The number of evaluations auto-scales with model complexity if global_maxeval is not set.

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

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.

The analytic R-matrix (analytic_cov_hessian = true, the default since #436) assembles H directly from third-order sensitivities of the closed-form prediction, differencing nothing at the population 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 exact 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. Differencing a machine-precision quantity amplifies round-off as ε/h rather than the second-difference stencil’s ε/h², so the ε^(1/3) step is the textbook optimum and fd_hessian_step does not enter.

Its scope is narrow, and deliberately so: plain analytical (closed-form) Gaussian models without IOV, LTBS, [scaling], Form-C readouts, [initial_conditions], iiv_on_ruv, correlated or custom-magnitude residual error, M3 censoring, FREM, a covariate-Selected error spec, or a non-Gaussian endpoint — and not under method = laplace/agq, whose marginal this assembly does not describe. The decision is all-or-nothing across the population: one out-of-scope subject returns the whole fit to the finite-difference stencil, because a Hessian half-assembled from two different approximations would be neither. Everything outside the scope keeps the stencil below, which is correct for all of it.

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.

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 FD Hessian (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 FD 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 FD Hessian 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.

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.

Configuration Example

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