Gradient-based local optimizers (D2D-style: TRF/Levenberg-Marquardt + L-BFGS-B) in the async Algorithm loop

[wshlavacek]

Summary

Add gradient-based local optimizers to PyBNF, following the Data2Dynamics (D2D) methodology
(Raue/Schilling/Timmer et al.) that is standard for ODE systems-biology parameter estimation. The
current optimizer suite is metaheuristic only (GA, PSO, SA, scatter search, …). This issue adds the
gradient/least-squares optimizers that consume the residual Jacobian surfaced by #385.

Depends on #385.

Re-scoped 2026-06-26. The original draft (a new pybnf/optimizers/ package built around
blocking scipy.optimize.minimize(jac=…)) was wrong on two counts, both corrected below:

1. PyBNF already has an optimizer package — pybnf/algorithms/optimizers/ (local_base.py,
   powell.py, cmaes.py, …). New optimizers go there, not in a second top-level package.
2. PyBNF's loop is asynchronous — Algorithm (pybnf/algorithms/base.py) drives an
   inversion-of-control propose/score loop (make_job → scheduler → got_result). A blocking
   scipy call that invokes fun/jac synchronously cannot drive PyBNF's distributed scheduler
   (the same incompatibility powell.py documents). The integration is therefore not "pass
   jac=True to scipy."

Reference implementation

examples/becker_d2d_gradient/ is a self-contained, runnable D2D-style fit of the Becker et al.
(2010) EpoR model using BNGsim forward sensitivities: multi-start → trust-region reflective
least-squares (scipy.optimize.least_squares(method="trf")) → profile likelihood. It runs
standalone (BNGsim directly, no PyBNF scheduler), so it is both the methodological target for this
issue and a working example of the residual + Jacobian path. (Profile likelihood itself is tracked
separately — see #446.)

The method (what to build, in priority order)

1. Primary: trust-region least-squares (TRF / Levenberg–Marquardt). This is the workhorse for
   Gaussian / SSR objectives — the common case in PyBNF. It consumes the residual vector + residual
   Jacobian from [Epic] Gradient plumbing: objective gradients & residual Jacobians from the BNGsim output-sensitivity tensor #385 and approximates the Hessian as JᵀJ, giving far better conditioning and
   convergence than feeding a scalar gradient to a generic quasi-Newton method. Build this first.
2. Fallback: L-BFGS-B on the scalar objective + gradient. For the cases where least-squares
   structure does not apply — heavy qualitative/inequality constraints, non-Gaussian likelihoods.
   Consumes Objective.evaluate_gradient ([Epic] Gradient plumbing: objective gradients & residual Jacobians from the BNGsim output-sensitivity tensor #385).

Current state (verified against main)

• pybnf/algorithms/optimizers/ holds the local/derivative-free optimizers (local_base.py,
  powell.py, simplex.py, scatter_search.py, …), each an Algorithm subclass that proposes
  parameter sets and consumes scored results asynchronously.
• No optimizer currently consumes a gradient or residual Jacobian ([Epic] Gradient plumbing: objective gradients & residual Jacobians from the BNGsim output-sensitivity tensor #385 surfaces them).
• fit_type in the .conf schema selects the algorithm by name.

Scope of this issue

1. Least-squares optimizer (primary). A new Algorithm subclass under
   pybnf/algorithms/optimizers/ implementing a bounded trust-region (TRF-style) /
   Levenberg–Marquardt step that consumes the residual + residual-Jacobian form from [Epic] Gradient plumbing: objective gradients & residual Jacobians from the BNGsim output-sensitivity tensor #385.
2. Gradient quasi-Newton optimizer (fallback). A bound-constrained L-BFGS-B Algorithm subclass
   consuming the scalar gradient from [Epic] Gradient plumbing: objective gradients & residual Jacobians from the BNGsim output-sensitivity tensor #385.
3. Multi-start orchestration. Wire both so PyBNF's existing multi-start / refinement machinery can
   seed them (gradient local methods are local — multi-start is how D2D handles non-convexity),
   consistent with the other local optimizers.
4. Config integration. Extend fit_type to accept the new optimizer names (e.g.
   fit_type = trf, fit_type = lbfgs) with optimizer_options parsed as for the existing local
   optimizers. Refuse to start (clear error, suggest a metaheuristic fit_type) when
   BNGSIM_HAS_OUTPUT_SENS is False, when the model has discrete events (BNGsim raises on
   sensitivity requests there), or when no gradient is available for the chosen method.

On wrapping scipy (decision)

The reference notebook calls scipy.optimize.least_squares(method="trf") directly and it works
perfectly — because it drives BNGsim in-process with no PyBNF scheduler. So for a standalone /
manuscript path, scipy TRF is exactly right and need not be reinvented. The async incompatibility
only appears inside PyBNF's propose/score loop, where a blocking least_squares cannot farm its
evaluations to distributed workers. Two ways to resolve it:

• (a) Hand-roll the TRF/LM and L-BFGS-B steps inside the Algorithm loop — full control of the
  async boundary; recommended for the first cut.
• (b) Bridge scipy onto a worker so its synchronous fun/jac callbacks resolve against the
  scheduler.

Prove one gradient method end-to-end via (a) before generalizing. The earlier
OptimizationResult/registry sketch is fine as an internal shape but is not a reason to add a
parallel pybnf/optimizers/ package.

Edge cases and design considerations

• No gradient available → fail fast with a clear message pointing at a metaheuristic fit_type,
  never a silent finite-difference fallback.
• Native vs sampling space. [Epic] Gradient plumbing: objective gradients & residual Jacobians from the BNGsim output-sensitivity tensor #385 delivers the gradient already transformed once via
  priors/scale.py (ADR-0029); the optimizer works in sampling space and must not re-transform.

Deliverables

Status updated 2026-06-29 against main. All eight deliverables done — the discrete-events pre-flight gate landed in 9021e38 (#461) and the convergence-vs-baseline + becker smoke tests in b3dc654 (#462).

• Trust-region / Levenberg–Marquardt least-squares Algorithm subclass under
  pybnf/algorithms/optimizers/, consuming the residual + residual-Jacobian form ([Epic] Gradient plumbing: objective gradients & residual Jacobians from the BNGsim output-sensitivity tensor #385). — trf.py (2bfca26).
• L-BFGS-B Algorithm subclass consuming the scalar gradient ([Epic] Gradient plumbing: objective gradients & residual Jacobians from the BNGsim output-sensitivity tensor #385). — lbfgs.py (7a170ed; upgraded to full L-BFGS-B, generalized Cauchy point + subspace minimization, in f73e3d0).
• Multi-start seeding wired through existing orchestration. — bde946b: GradientOptimizer runs N box-sampled starts concurrently (start 0 = box center, the rest Latin-hypercube) and keeps the global best.
• fit_type + optimizer_options config integration for the new names. — trf / lbfgs are registered and selectable; options landed as co-located Pydantic schemas (TRFConfig / LBFGSConfig, per-key e.g. trf_grad_tol, lbfgs_history) per ADR-0006, not a literal optimizer_options blob.
• Capability / gradient-availability / events gate with clear errors and a suggested fallback. — the edition + sensitivity-backend + output_sensitivities capability gate, the per-evaluation GradientNotSupported refusal, and the explicit discrete-events pre-flight gate (9021e38, Gradient path: add a discrete-events pre-flight gate (refuse cleanly, don't fail mid-run) #461) are all done.
• Tests: each optimizer recovers known parameters on a small FD-validated model; convergence vs a
  metaheuristic baseline on the same model; a smoke test mirroring examples/becker_d2d_gradient/. — recovery tests (tests/test_gradient_optimizer.py: rate/IC, estimated-σ, bound-active, multimodal multi-start) + offline scipy-oracle runner parity (tests/test_gradient_runner.py); the convergence-vs-metaheuristic-baseline test (each gradient fit recovers the truth in ≥10× fewer evaluations than de, to an at-least-as-good objective) and the becker smoke test (trf multi-start through the scheduler on the Antimony → SBML BIOMD0000000271 EpoR model) both landed in b3dc654 (Gradient optimizers: finish #386 deliverable 6 — convergence-vs-baseline test + becker smoke test #462). Merging the examples/becker-d2d-gradient reference branch remains optional and separate.
• User-guide section for the gradient-based fit_types. — docs/gradient_fitting.rst ("Running a gradient fit" + the assembly / cost / parameter-scale sections).
• No behavioral change for existing metaheuristic fits. — the gradient path is purely additive and inactive unless fit_type = trf / lbfgs; existing suites green.

Out of scope

• Third-party optimizer libraries (fides, cyipopt/IPOPT, nlopt) — each a follow-up.
• A generic scipy wrapper/registry as the primary integration path (see decision above).
• Profile-likelihood analysis — tracked in Profile-likelihood analysis (D2D-style identifiability + confidence intervals) #446 (built on the same residual Jacobian).
• Second-order / adjoint sensitivities (forward-mode, first-order only — see [Epic] Gradient plumbing: objective gradients & residual Jacobians from the BNGsim output-sensitivity tensor #385).