5. Nonlinear solver: scipy.optimize.least_squares with normalization and restarts#

Status: accepted · 2026-10-22

Context. The v0.1 goal is that all 27 NIST StRD nonlinear problems reach the certified values (≥ 6 significant digits) from both official starting points. The original plan wrapped lmfit.minimize. Measured with tight tolerances:

Setup

Start 1

Start 2

lmfit leastsq, default tolerances

17/27

16/27

lmfit least_squares (TRF), tol 1e-12

21/27

22/27

scipy least_squares LM, 2-point Jacobian

20/27

23/27

scipy TRF, 3-point Jacobian, normalized parameters, restarts

27/27

27/27

Three things made the difference:

  1. Normalized parameters u = p / |p_start|. Without them Hahn1 (coefficients between 1 and 1e-7) stops early at LRE 2.3. scipy’s x_scale option does not have the same effect.

  2. Central (3-point) finite differences for the Jacobian. Kirby2 and Lanczos3 need them.

  3. Restarts from the solution with fresh normalization, until χ² stops improving. They fix MGH09, MGH17, BoxBOD and Bennett5 from start 1.

Wrapping lmfit’s minimizer leaves no room for the first and third measures.

Decision. The default backend calls scipy.optimize.least_squares directly, with Trust Region Reflective, 3-point finite differences, normalized parameters and up to 10 restarts. TRF also handles bounds, so bounded and unbounded fits use the same solver. Levenberg–Marquardt stays available as method="lm". For LM, fitwise computes the central differences itself, because older scipy versions (for example 1.12) silently fall back to 2-point for lm. lmfit stays a dependency for profile likelihood (conf_interval, v0.3) and expression models (asteval, v0.2).

The covariance is (JᵀJ)⁻¹ from an SVD of the final Jacobian in normalized coordinates. The Jacobian counts as rank deficient if s_min / s_max ≤ 1e-9. Redundant parameters reach about 1e-12 there, and the worst well-posed NIST problem (Nelson) 3e-6.

Consequences.

  • NIST StRD 27/27 from both starting points, about 0.25 s for all 27 problems.

  • Results agree with curve_fit to better than 1e-5 standard errors, and the covariance agrees to 1e-4 relative (different finite-difference schemes).

  • FitResult.raw is a scipy.optimize.OptimizeResult, not an lmfit object.

  • Restarts cost extra model evaluations (typically 2–3 runs), which is negligible for curve fitting.