Outer Optimizers

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

The outer optimizer is the algorithm that minimises the population objective function (OFV) over the transformed parameter vector [log θ, chol(Ω), log σ].

It only applies to methods that have a single outer-loop minimisation over that vector:

method does optimizer apply?
foce / focei yes — picks the algorithm for the population OFV minimisation
gn (pure Gauss-Newton) no — GN runs its own Levenberg-Marquardt loop
gn_hybrid yes for the FOCEI polish phase only; the preceding GN phase ignores it
saem no — SAEM has no single outer minimisation; the M-step uses a hardcoded NLopt algorithm and is not user-selectable. Setting optimizer under method = saem triggers a warning and is otherwise ignored.
[saem, focei] (chain) applies to the FOCEI polish stage only
imp no — IS is a sampling pass, not an optimisation

Set via [fit_options]:

[fit_options]
  method    = focei
  optimizer = auto

The default optimizer = auto picks the algorithm per model (see auto below); name a specific optimizer to override it.

Available optimizers

Automatic selection

Key Algorithm Notes
auto Per-model selection Default. Resolves to nlopt_lbfgs when the exact analytic FOCE/FOCEI gradient is available, and bobyqa when only finite differences are or the model is a mixture. The fit output records the resolved pick as auto (<resolved>). See The auto default.

NLopt algorithms

Key Algorithm Notes
bobyqa Bounded Optimization BY Quadratic Approximation Derivative-free quadratic trust-region; what auto selects on FD-only fits. Avoids the fixed-EBE gradient bias that stalls gradient methods on ill-conditioned fits; consistently reaches a lower OFV on ODE/PD models, sparse data, and Hill-ridge problems without reconverging the inner loop.
slsqp Sequential Least Squares Programming Gradient-based, handles bounds well. Fast on well-conditioned analytical PK models; can stall above the true minimum on ODE/PD / sparse-data fits unless paired with reconverge_gradient_interval = 1 (5–6× cost).
nlopt_lbfgs L-BFGS via NLopt Limited-memory BFGS. Fast on analytical 1-/2-/3-cpt FOCE/FOCEI fits (exact analytic gradient); also useful for high-parameter-count models.
mma Method of Moving Asymptotes Alternative constrained gradient optimizer. Rarely needed.

Deprecated aliases

Key Algorithm Notes
bfgs → nlopt_lbfgs Deprecated alias. Now selects NLopt L-BFGS.
lbfgs → nlopt_lbfgs Deprecated alias. Now selects NLopt L-BFGS.
Note

bfgs and lbfgs previously selected a hand-rolled built-in BFGS / limited-memory L-BFGS. Benchmarking showed that path was strictly worse than NLopt L-BFGS — 3–5× slower on analytic-gradient models and prone to diverging or hanging on harder problems (#483) — so both keys are now aliases for nlopt_lbfgs. The built-in implementation is slated for removal.

Newton trust-region

Key Algorithm Notes
trust_region Newton trust-region with Steihaug CG Uses the analytic outer gradient and a BHHH approximate Hessian H ≈ 4 Σ gᵢgᵢᵀ (always positive semi-definite). The CG budget defaults to ceil(sqrt(n_params)).clamp(5, n_params); pin it with steihaug_max_iters. Best combined with inits_from_nca since it benefits from good starting values. Stops when it can no longer make progress — see convergence reporting.

Trust-region convergence reporting

trust_region stops once it has gone 20 consecutive iterations without measurable progress — no objective improvement beyond outer_ftol and no parameter movement beyond outer_xtol (default 1e-4, relative). outer_ftol is resolved exactly as it is for bobyqa: 1e-8 for a pure non-Gaussian (TTE / categorical) non-ODE objective and 1e-6 otherwise, with an explicit [fit_options] outer_ftol always winning. A run that instead exhausts maxiter is reported as Converged: NO, with a warning naming the budget it ran out of and the gradient norm it stopped at.

Because progress is an or — the objective improved or the parameters moved — and a rejected trust-region step moves neither, 20 no-progress iterations means 20 consecutive rejections: the radius has shrunk by roughly 0.25²⁰ ≈ 1e-12 from the last accepted step. That is the real content of the count, and it is why it does not need tuning per model.

Two consequences worth knowing:

  • Demonstrating settling costs at least 21 outer iterations, so a fit run with maxiter below that can never report Converged: YES however good it is. The warning says so explicitly rather than blaming the fit.
  • A run that rejects every step from the first iteration has not converged, it failed to start. That is reported as its own non-convergence reason — pointing at the starting values, not the budget — and the fit bails out there instead of spending the rest of maxiter on a frozen state.

Before #1000 every trust_region run reported Converged: YES, including one that had simply run out of iterations: the underlying solver has no convergence criterion of its own, so hitting maxiter was its only way to stop. On examples/one_cpt_iv_pooled.ferx with data/one_cpt_iv.csv:

maxiter OFV reported (before) reported (now)
5 11325.8283 Converged: YES Converged: NO (budget below the 21-iteration floor)
300 8504.4349 Converged: YES Converged: NO
2000 −269.6370 Converged: YES Converged: YES (at iteration 1920)

The true optimum is −269.6370 (NONMEM $OMEGA 0 FIX anchors it at −269.63700440), so the first two rows were 8 000–11 000 OFV units off and reported as converged, with standard errors computed at a non-stationary point.

The stopping rule is deliberately not a gradient test. The outer gradient is evaluated at EBEs held fixed within the iteration, which makes its norm an unreliable stationarity measure on an ODE model: on one_cpt_iv_ode the norm at the optimum the default optimizer reaches (OFV −387.0737) is 43.2, while a worse point 0.73 OFV units away carries 4.0 — not monotone in solution quality, so no threshold separates the two. Objective and step size are the noise-robust measures, and they are what bobyqa already stops on.

Stopping early also saves work: the warfarin fit settles at iteration 59 and now returns there rather than grinding out the remaining budget.

Note that “converged” still means this optimizer stopped making progress, not that the point is a verified global optimum — a trust_region fit can settle slightly above the optimum other optimizers reach (#864). Compare against a second optimizer when the OFV matters.

One more thing changes for multi-start users: n_starts ranks candidates by validity, then by converged, then by OFV. Every trust_region start used to report Converged: YES, so that middle key was inert and selection came down to OFV. Now that budget exhaustion reports honestly, a start that settled can be preferred over one that was still descending when it ran out of iterations, even if the latter reached a lower OFV. Give trust_region multi-starts a maxiter generous enough that the good start is allowed to settle.

When to use which

auto (default) — leave this in place for most fits. It inspects the model and picks nlopt_lbfgs when the exact analytic FOCE/FOCEI gradient is available and bobyqa when only finite differences are — the two choices that came out fastest-and-most-reliable in the #490 benchmarking (see The auto default). Override it only when you have a model-specific reason to pin a particular optimizer.

bobyqa — what auto selects on FD-only problems, and a fine explicit choice there. Derivative-free quadratic trust-region; re-evaluates EBEs at every trial point so it avoids the fixed-EBE gradient bias that stalls gradient-based optimizers on ill-conditioned fits. On the cefepime 2-cpt benchmark below it reaches a lower OFV than slsqp (even reconverged at full cost) in less wall time. Works equally well on smooth analytical PK models — converges more slowly than gradient methods per iteration but each iteration is cheap (no FD gradient sweep).

slsqp — switch to this when bobyqa is too slow on a smooth, well-conditioned model with many parameters (it can need many quadratic-interpolation samples to triangulate a high-dimensional surface). Gradient-based, handles box constraints cleanly; pair with reconverge_gradient_interval = 1 if it stalls above an expected OFV.

trust_region — for models with many parameters (large θ + Ω) or when combined with inits_from_nca. The second-order curvature information helps when starting values are already in the basin.

gn / gn_hybrid — these are Gauss-Newton estimation methods, not outer optimizers in the same sense. They replace the FOCE outer loop entirely rather than selecting an algorithm within it. (gn_hybrid polishes via FOCEI and inherits the optimizer setting for that stage — so the polish runs with the auto-selected optimizer unless overridden.)

lbfgs / nlopt_lbfgs — strong choices on analytical 1-/2-/3-cpt FOCE/FOCEI fits, where they get an exact bias-free gradient (see Analytic FOCE / FOCEI gradient) that reaches the true optimum at wall-time comparable to bobyqa (and several × faster than the now-deprecated built-in BFGS). Since #410 they also get the exact gradient on in-scope user-[odes] models (bolus/infusion dosing, ObsCmt/Form-C/per-CMT readout, F, resets, static covariates, and — since #439 — time-varying covariates on bolus dosing, steady-state dosing #473, IOV #466, IOV combined with an ExpressionScale obs_scale = expr divisor #575, and IOV + time-varying covariates + ExpressionScale #590), so lbfgs is a good choice there too; out-of-scope ODE features (input_rate/SDE, LTBS combined with an ExpressionScale divisor, non-IOV TV-covariates combined with an ExpressionScale divisor, IOV combined with LTBS, and TV-covariates combined with infusion/SS/reset) still fall back to the finite-difference gradient and its fixed-EBE bias — prefer bobyqa for those. mma is rarely needed.

Why these two? bobyqa doesn’t use the gradient, so it never sees the fixed-EBE FD gradient bias that drives gradient-based optimizers to local minima hundreds of OFV units above the true optimum on ODE/PD models and sparse data (the Emax PKPD benchmark in saem.md and the cefepime validation below). On analytical PK models the gradient is exact and bias-free, and nlopt_lbfgs exploits it to reach the optimum faster than the derivative-free search. auto (default since #490) routes each model to the right one of the two automatically; before auto the standalone default was bobyqa, and before that slsqp.

The auto default

optimizer = auto (the default) resolves to a concrete optimizer from the model, with no user input:

  • nlopt_lbfgs — when the exact analytic FOCE/FOCEI gradient is available: an analytical PK model (or an in-scope user-[odes] model, see Analytic FOCE / FOCEI gradient) that is in the sensitivity provider’s scope and is not forced onto finite differences via gradient = fd.
  • bobyqa — otherwise: any model outside the analytic scope (for example an out-of-scope ODE/PD model or an SDE model), or any fit with gradient = fd, where the outer loop runs on finite differences.

Two cases override that rule:

  • More than 64 free packed parameters (BOBYQA_MAX_DIM) → nlopt_lbfgs even on a finite-difference gradient. BOBYQA’s interpolation set grows quadratically with the parameter count, so above that size O(n) FD passes are cheaper (#1064).
  • Mixture models ([mixture]) → bobyqa regardless of gradient availability, as the robust choice against label-switching multimodality (#977).

The choice mirrors the analytic-vs-FD decision the outer loop already makes, so auto lands on the gradient-based optimizer exactly when there is an exact gradient to feed it. Limited benchmarking across ~10 real FOCEI datasets (#490) found nlopt_lbfgs fastest-to-optimum on every analytic-gradient problem (bobyqa was faster but often stopped short of the optimum, slsqp converged but much slower), while bobyqa was both fastest and most reliable when only finite differences were available (nlopt_lbfgs/slsqp/mma were very slow and often missed the basin).

For a FOCE/FOCEI fit the resolved optimizer is recorded in the fit output as auto (<resolved>) — e.g. auto (nlopt_lbfgs), or auto (bobyqa) for a mixture model — so a run records which algorithm actually executed (#1540). An explicit trust_region (or, from the API, the built-in BFGS/L-BFGS) that the mixture override replaces with BOBYQA is reported as bobyqa too, with a warning; an explicit NLopt gradient optimizer (slsqp / nlopt_lbfgs / mma) is honoured. One exception: gn_hybrid reports gn even though its FOCEI polish runs the resolved optimizer. These defaults may be refined as the analytic-gradient scope grows and more datasets are studied; pin optimizer explicitly to opt out of the automatic choice.

gradient = fd moves the optimizer too

Because auto resolves off the availability of the analytic gradient, changing one line in [fit_options] — gradient = auto to gradient = fd — changes two factors: the gradient the outer loop computes, and the optimizer that consumes it. On the model that prompted #1381 that one-line change moved the OFV by 4.98; with optimizer = bobyqa pinned in both arms, the gradient alone accounted for 0.38. The other 4.6 was the optimizer swap.

The default is not wrong — picking the optimizer that suits the available gradient is the whole point of auto, and it rests on the #490 benchmarking above. But a controlled one-variable experiment has to pin the other variable:

[fit_options]
gradient = fd
optimizer = bobyqa      # pin it in BOTH arms, so the difference is the gradient

ferx warns when you don’t. Setting gradient = fd while leaving optimizer at auto, on a model whose analytic gradient is in scope, emits W_AUTO_OPTIMIZER_FOLLOWS_GRADIENT naming both the optimizer that ran and the one the unforced arm would use. ferx check reports it too, against the gradient line itself, so the coupling surfaces before the fit is run.

The warning is silent whenever nothing was actually coupled — optimizer pinned; a model out of analytic scope anyway, where auto picks bobyqa either way; a model above BOBYQA_MAX_DIM, where auto takes nlopt_lbfgs on a finite-difference gradient regardless; a mixture model, whose auto is bobyqa whatever the gradient; an evaluation-only (maxiter = 0) run, which constructs no optimizer at all; and an SDE model, which the engine forces onto finite differences whatever you asked for.

Two related fields record the effect after the fact rather than the coupling: optimizer reports the compound "auto (bobyqa)", and final_gradient_source reports "finite_difference" for an arm that silently became derivative-free.

Leaving the initial estimates

Both gradient-based NLopt optimizers start with their quasi-Newton Hessian set to the identity, so their opening search direction is simply \(-\nabla f\). On a typical PK model the scaled analytic gradient is large at the initial estimates (\(10^2\)–\(10^3\) in log/Cholesky space), so that first step lands in a corner of the parameter box, the objective there explodes, and the line search fails — leaving the fit sitting on (or a hair off) its initial estimates for the rest of the budget.

ferx rescales the gradient by a single scalar on that opening step, so it stays inside a sane per-dimension budget. The direction is unchanged. The rescale is deliberately not held on afterwards: L-BFGS reconstructs its Hessian from successive gradient differences, and rescaling those differences corrupts the curvature it has learned (SLSQP re-solves its subproblem from scratch each step, so it is rescaled throughout).

If the fit nevertheless ends up still on its initial estimates, ferx re-runs it once from the same start with the rescale held on until the fit escapes — a fit that never moved has no curvature worth protecting. The second attempt is reported only if it both left the initial estimates and reached a lower OFV; if it stalled as well, the first attempt stands (together with its “did not converge” warning), so a fit reported as non-convergent is never one that quietly picked up a lower objective at a point it never reached. A fit that was descending normally never triggers any of this and follows exactly the trajectory it always did.

A fit that nonetheless fails to leave its initial estimates through this path is reported with converged = false and a “did not converge” warning — never as a converged fit whose estimates and standard errors happen to be the starting values. If you see this, the usual remedies are better initial estimates, a different optimizer, or parameter_scaling (see fit_options).

That guarantee covers the case where the optimizer itself reports failure. It does not cover a derivative-free optimizer that reports success at the start — BOBYQA can meet its own outer_xtol without having travelled anywhere, and then converged = true is a statement about the step size, not about the objective. Since #997 every fit that ends with no free THETA, OMEGA or SIGMA coordinate moved carries a stalled_at_init warning regardless of what the optimizer concluded. The check is deliberately all-or-nothing — it fires only when nothing moved. A per-parameter version would be worse than useless: an initial estimate close to the truth is one a good optimizer is supposed to leave alone, so a ratio test fires hardest on the arm that is right.

Checking a derivative-free converged

converged reports that an optimizer’s stopping rule was met. For the gradient-based optimizers the result also carries final_gradient, so the claim can be checked independently: at a genuine interior optimum it is ≈ 0. A derivative-free run has no gradient to report — which is the wrong way round, because a derivative-free run is where a premature stop is most likely and where a bad answer is hardest to see. On the stiff binding model in #997 a BOBYQA fit stopped 1.54 OFV above the L-BFGS optimum on the same model and data, and reported converged = true with nothing to check it against; widening the parameter box a thousandfold returned a bit-identical objective to six decimals, because it had never travelled far enough for a bound to matter.

So when the optimizer supplies no gradient, ferx computes one after the fit:

field meaning
final_gradient ∇ of the objective the fit minimised, at the reported estimates, in packed space (log-THETA, Cholesky-OMEGA, log-SIGMA)
final_gradient_source "optimizer" — the optimizer’s own, evaluated during the fit; or "finite_difference" — computed afterwards purely for reporting

The finite-difference version is a central difference of the same (penalized) objective, with the EBEs re-solved at every perturbed point, so it is the gradient of the marginal objective rather than a fixed-EBE approximation to it. It costs 2 × n_free objective evaluations — one gradient’s worth, the same as a single step of the optimizers that report one for free — and it never steers anything: the estimates are bit-identical whether it runs or not. Set report_final_gradient = false in [fit_options] to skip it when that cost is not worth paying.

Both kinds are read the same way, but they do not prove the same thing. An "optimizer" gradient near zero is a stationarity claim the optimizer itself acted on. A "finite_difference" gradient near zero is an independent check on a run that had no gradient at all — and a large one says the fit stopped somewhere the objective was still moving, whatever converged says.

What it reads on a fit that is right, and one that is not

Measured on the bundled examples/warfarin.ferx + data/warfarin.csv under method = focei, maxiter = 300, fitting the same model from the same start with the two optimizers:

optimizer converged OFV final_gradient_source ‖∇‖∞
nlopt_lbfgs true −286.0042 optimizer 2.29e−4
bobyqa true −285.9866 finite_difference 2.47e0

Both report success. BOBYQA stops 0.018 OFV short with a gradient four orders of magnitude larger, and before #997 that second column was the only difference visible in the fit object — the optimizer’s own gradient being absent for exactly the arm that needed it.

The two numbers above come from different stencils, so to check that the comparison is fair the same finite-difference stencil was also run at the L-BFGS solution: it reads 3.04e−4 there, against the optimizer’s own 2.29e−4. So the reporting gradient does go to zero at a genuine optimum on this model; the 2.47 is BOBYQA’s, not the stencil’s. The step is h = 1e-4·(1+|x|) per packed coordinate, which sits in the flat part of the step-size curve — sweeping h over 1e-4 … 1e-2 moves the BOBYQA reading by less than 0.1 %, while h = 1e-1 inflates both arms to ~3.5e1 on truncation error alone.

One caveat on reading the magnitude: it is the gradient in packed space (log-THETA, Cholesky-OMEGA, log-SIGMA) and it is not normalised, so there is no universal threshold. It is a comparative quantity — against a second fit of the same model, or against the optimizer gradient of a gradient-based arm — not a pass/fail test. A weakly-identified direction can carry a sizeable ∂OFV/∂x across a genuinely flat objective, which is the same reason ferx does not use a gradient norm as its convergence criterion (see Leaving the initial estimates).

Fixed-EBE gradient bias and reconverge_gradient_interval

By default, gradient-based optimizers (slsqp, lbfgs, etc.) hold each subject’s EBEs fixed while computing the population gradient. This is cheap but omits the response of the inner solution to θ and Ω. On ill-conditioned models the omitted term causes slsqp to stall well above the bobyqa optimum.

Set reconverge_gradient_interval = 1 to re-solve the inner loop at every gradient evaluation, recovering the full gradient surface at roughly 5–6× cost. A value of N reconverges every Nth evaluation and uses the cheap gradient in between — often enough to close most of the gap:

Configuration OFV Wall time
slsqp (fixed-EBE) 68,252 390 s
slsqp + interval = 10 66,118 633 s
slsqp + interval = 1 65,485 1,871 s
bobyqa (derivative-free) 65,598 315 s

On this cefepime 2-compartment dataset bobyqa reaches a near-optimal solution faster than any reconverged slsqp setting. The reconverged gradient is most useful when derivative-free search is too slow (high parameter count) or when a gradient optimizer is required for other reasons.

IOV models (kappa/block_kappa) always reconverge regardless of this setting. Even so, a pure slsqp cold-start on an IOV model can terminate a few OFV units above the minimum, and the exact stopping point is platform-dependent (the re-converged FD gradient’s summation order differs across architectures — issue #160). Prefer optimizer = bobyqa or a methods = [saem, focei] chain for IOV fits; both reach the minimum platform-independently. See Inter-Occasion Variability.

Analytic FOCE / FOCEI gradient (analytical PK models)

For analytical 1-, 2-, and 3-compartment models (IV bolus/infusion, oral, and steady state) under FOCE or FOCEI, the gradient-based optimizers (bfgs, lbfgs, nlopt_lbfgs, slsqp) automatically switch from the finite-difference gradient to an exact closed-form marginal gradient (Almquist et al. 2015). FOCEI differentiates the Laplace marginal; FOCE differentiates ferx’s Sheiner–Beal linearized marginal. Both are computed analytically through second-order dual numbers and include the full EBE response (the term reconverge_gradient_interval recovers by brute force) — so they do not carry the fixed-EBE bias described above, at no extra cost. There is nothing to configure: the method and optimizer choices are the switch, and any model outside the analytical scope falls back to the finite-difference gradient transparently. A model inside the scope can still have individual subjects the provider declines; those subjects alone take the fixed-EBE gradient, and the fit reports them in a warning (see the outer-loop route).

Because the gradient is both exact and bias-free here, lbfgs / nlopt_lbfgs are strong choices on these models: each gradient evaluation is a closed form rather than an O(n) finite-difference sweep, and — unlike the FD gradient — it carries no fixed-EBE bias, so the optimizer reaches the true optimum. Wall-time versus the derivative-free bobyqa is model-dependent and roughly comparable (e.g. warfarin FOCEI: lbfgs ≈ 0.29 s vs bobyqa ≈ 0.43 s; on a 2-cpt fit bobyqa is faster but can stop short of the optimum), and both are several times faster than the now-deprecated built-in BFGS.

The exact gradient lands on the same optimum as the underlying objective, validated against NONMEM on the warfarin 1-cpt oral model: FOCE reaches OFV −280.36 (NONMEM −280.36; TVCL 0.1330, TVKA 0.7252, ω²: 0.0286 / 0.00958 / 0.349) and FOCEI reaches −286.00 (NONMEM −286.00) — both agreeing to ~4–5 significant figures across θ, Ω, and σ.

The exact gradient also covers log-transform-both-sides (log(DV) ~ additive(...)) and constant output scaling ([scaling] obs_scale = k): the analytic provider applies the g = ln(f) jet transform (value, gradient, and Hessian) and the constant divisor in closed form. Validated against NONMEM on the warfarin LTBS model — the gradient-based lbfgs path reaches OFV −675.302 and recovers NONMEM’s MLE (TVCL 0.1327, TVV 7.738, TVKA 0.811) to ~4 significant figures.

The fallback (finite-difference) gradient is used when the model has: an ODE system, inter-occasion variability (kappa), a dose lagtime, time-varying covariates, system resets, expression/per-compartment output scaling, or an overlapping steady-state infusion (T_inf > II, #379). For those, the guidance above (prefer bobyqa, or reconverge slsqp) still applies.