NonlinearSolveAlg(LevenbergMarquardt) returns ReturnCode.Success with wrong answers (successor to #4125)

[ChrisRackauckas-Claude]

Summary

NonlinearSolveAlg(LevenbergMarquardt()) returns ReturnCode.Success with answers that are
wrong by many orders of magnitude. This is the last surviving member of the family tracked in
#4125, which was auto-closed by #4142's Fixes #4125 — #4142 repaired the DFSane and
Broyden + BackTracking cases and explicitly did not repair LevenbergMarquardt, so the
tracking for it was lost. This issue is the successor.

Measured on merged master (5cf083900), reltol = 1e-6, abstol = 1e-9, maxiters = 1e7,
error = max absolute deviation from Rodas5P at reltol = 1e-13 / abstol = 1e-15:

problem	ODE alg	inner alg	retcode	naccept	abs err
ROBER	KenCarp4	NewtonRaphson	Success	7177	9.2830e-10
ROBER	KenCarp4	LevenbergMarquardt	Success	31609	5.0291e+07
ROBER	FBDF	NewtonRaphson	Success	634	8.8870e-10
ROBER	FBDF	LevenbergMarquardt	MaxIters	3846912	1.2830e-02
OREGO	KenCarp4	NewtonRaphson	Success	142	2.2699e-01
OREGO	KenCarp4	LevenbergMarquardt	Success	1082095	1.2156e+02
OREGO	FBDF	NewtonRaphson	Success	740	4.2322e-02
OREGO	FBDF	LevenbergMarquardt	Success	103360	1.4069e+03

The FBDF/ROBER row fails loudly, which is fine. The three bolded rows are silent wrong
answers.

Reproduction

using OrdinaryDiffEqSDIRK, OrdinaryDiffEqRosenbrock
using OrdinaryDiffEqNonlinearSolve: NonlinearSolveAlg
using NonlinearSolve: LevenbergMarquardt
using ADTypes

function rober!(du, u, p, t)
    y1, y2, y3 = u
    k1, k2, k3 = p
    du[1] = -k1 * y1 + k3 * y2 * y3
    du[2] = k1 * y1 - k2 * y2^2 - k3 * y2 * y3
    du[3] = k2 * y2^2
    return nothing
end
prob = ODEProblem(rober!, [1.0, 0.0, 0.0], (0.0, 1.0e5), (0.04, 3.0e7, 1.0e4))

ref = solve(prob, Rodas5P(); reltol = 1.0e-13, abstol = 1.0e-15)
sol = solve(
    prob,
    KenCarp4(nlsolve = NonlinearSolveAlg(LevenbergMarquardt(; autodiff = AutoForwardDiff())));
    reltol = 1.0e-6, abstol = 1.0e-9, maxiters = 10^7
)

sol.retcode          # ReturnCode.Success
sol.stats.naccept    # 31609   (NewtonRaphson needs 7177)
sol.u[end]           # wrong by 5.0e+07
ref.u[end]           # [0.017865921142099082, 7.274751468436061e-8, 0.9821340061103996]

Root cause

Instrumented nlsolve! / compute_step! (printing ndz, θ, η, the relative stage
residual, and the linear-model defect ‖W·δ + fu_prev‖ / ‖fu_prev‖ at every outer iteration).

LevenbergMarquardt in NonlinearSolve is a GeneralizedFirstOrderAlgorithm built from a
DampedNewtonDescent (wrapped in GeodesicAcceleration). Its step therefore solves a damped
Newton system, and — measured directly — covers only a fraction c ∈ (0, 1) of the Newton
step:

W·δ = -c·fu_prev      and      fu_new = (1 - c)·fu_prev

The instrumentation confirms this identity exactly. ROBER + KenCarp4 at t ≈ 1.44e5, per outer
iteration (linmodel is ‖W·δ - fu_prev‖/‖fu_prev‖, so linmodel = 1 + c; fratio is
‖fu_new‖/‖fu_prev‖ = 1 - c):

LM  iter=1 ndz=0.00344 theta=0.3286 eta=0.4894 eta*ndz=0.00168 kappa=0.01
    linmodel=1.4986 fratio=0.5014  br_eta=true      <- accepted, c = 0.4986
LM  iter=1 ndz=0.00142 theta=0.3286 eta=0.4894 eta*ndz=0.00069 kappa=0.01
    linmodel=1.4988 fratio=0.5012  br_eta=true      <- accepted, c = 0.4988

For contrast, NewtonRaphson on the same problem at t ≈ 1.47e5 (linmodel = 2.0000000000000013,
i.e. c = 1 to 15 digits — the displacement is the Newton step):

NR  iter=3 ndz=0.18144 theta=0.01452 eta=0.01473 eta*ndz=0.00267 kappa=0.01
    linmodel=1.9999999999982587 fratio=0.01340  br_eta=true

Now the failure. Deeper into the ROBER integration the damping stops covering any of the
Newton step at all, and the outer test accepts anyway. First accepted step whose state departs
from the reference is #3743 at t = 1.4552e5; the two steps around it:

=== TRACE step 3742 t=144653.26968717773 dt=208.0294931836172
NLSOLVE t=144653.26968717773 dt=329.59585745186956 kappa=0.01 eta0=0.1570944686654572
        relres0=(0.0035475549780305117, 1.1502254702762665e-9)
  iter=1 ndz=0.08715022180015697 theta=0.09 eta=0.0989010989010989
        eta*ndz=0.00861925270551003 kappa=0.01
        relres=(0.0036797750288398907, 1.193004796772091e-9)
        br_iter1=false br_eta=true

Read that line: one outer iteration; the stage residual rose from 1.1502e-9 to
1.1930e-9 (c ≈ -0.037, i.e. the step accomplished nothing); and the stage was nonetheless
declared converged because η·ndz = 0.0086 < κ = 0.01.

η there is 0.0989, from θ = 0.09. That θ is stale. nlsolve! reads θ off the
ratio of consecutive displacements, which only exists from the second iteration of a stage
onwards; on the first iteration it falls back to nlsolver.prev_θ, carried over from an
earlier stage. Once these stages start converging at iter == 1, no second iteration ever
runs, prev_θ is never re-measured, and it stays frozen at 0.09 for the rest of the solve.
The convergence test degenerates into the bare displacement threshold ndz < κ/η = 0.101,
with a contraction factor measured on a different stage at a different step size.

The consequence is that every stage is accepted with a fraction of the correction applied. The
embedded error estimator does not see it (both stages are wrong the same way), so the step is
accepted, dt collapses (31609 accepted steps against NewtonRaphson's 7177), and the
solution drifts to 5e7 while retcode stays Success.

Why #3893 and #4142 do not fire

• NSA: gate iter-1 convergence on inner solver retcode (#3817) #3893 (_uninformative_step) requires the inner step! to have left the iterate
  exactly unmoved (iszero(ndz)) with the cache not terminated — a TrustRegion that
  rejected its trial step. LM's accepted steps move the iterate by ndz ≈ 0.087 in the
  weighted norm, ~13 orders of magnitude above the predicate's trigger.
• NonlinearSolveAlg: let the residual veto a stalled inner step #4142 (stalled_inner_step / stage_unsolved) requires the displacement to be below
  8·eps·‖z‖ and the residual not to have decreased. LM's displacement is far above
  roundoff, and over most of the integration its residual does decrease (by the factor
  1 - c ≈ 0.5 above), so the predicate is false on both counts.

Both predicates ask "did the inner solver produce a step at all?". LM does produce a step. The
question neither of them asks is "is the step it produced the Newton step the η/κ
contraction argument assumes?", and for a damped method the answer is no.

Notes for a fix

A plain "the residual must be small" gate does not work here, for the reason recorded in #4125:
for a stiff stage the residual at a perfectly converged iterate is (M - γΔt·J)·dz, amplified
by the stiffness, so no threshold on ‖γΔt·fu‖ separates a converged NewtonRaphson stage
from a stalled LevenbergMarquardt one. (Measured: at declared convergence NewtonRaphson
sits at relative residual ~8e-6 on this problem, not at any machine-precision floor.)

What does discriminate, costs no linear algebra, and is invariant both to the problem's units
and to the inner solver's deliberately-zeroed tolerances, is the residual contraction the
step actually achieved, ‖fu_new‖ / ‖fu_prev‖. Under the linearisation the residual and the
error are related by the constant W, so they contract by the same factor; unlike the
displacement ratio it needs no predecessor iteration, so it is available exactly where the
current code is forced to inherit a stale θ. In the table above it reads 0.5014 where the
damping is covering half the Newton step (correctly permitting convergence: the remaining
error really is ≈ ndz), and 1.037 at the failing stage (correctly refusing).

/cc @ChrisRackauckas

[ChrisRackauckas-Claude]

Correction: the reproduction in the issue body above does not reproduce. My apologies — the snippet builds the Robertson problem with tspan = (0.0, 1.0e5) but quotes naccept and the error from the 1.0e11 problem, so the two halves came from different runs. Anyone pasting it sees a clean solve.

Verified on master 284487ebf, KenCarp4, rtol 1e-6 / atol 1e-9, error vs Rodas5P at 1e-13/1e-15:

problem	inner	retcode	naccept	abs_err
ROBER tspan=(0,1e5) (as written above)	LevenbergMarquardt	Success	2458	1.4558e-06 — correct
ROBER tspan=(0,1e5)	NewtonRaphson	Success	154	9.0954e-10
ROBER tspan=(0,1e11)	LevenbergMarquardt	Success	31609	5.0291e+07
ROBER tspan=(0,1e11)	NewtonRaphson	Success	7177	9.2830e-10
OREGO tspan=(0,30) FBDF	LevenbergMarquardt	Success	103360	1.4069e+03
OREGO tspan=(0,30) FBDF	NewtonRaphson	Success	740	4.2322e-02

The tspan matters because the solution only departs from the reference at t ≈ 4.6e7, when u₁ crosses zero and the drift along the Robertson conservation direction goes unstable — six orders of magnitude past the end of the 1e5 problem.

Corrected reproduction

using OrdinaryDiffEqBDF, OrdinaryDiffEqRosenbrock
using OrdinaryDiffEqNonlinearSolve: NonlinearSolveAlg
using NonlinearSolve: NewtonRaphson, LevenbergMarquardt
using ODEProblemLibrary, ADTypes

prob = ODEProblemLibrary.prob_ode_orego          # tspan (0, 30)
ref  = solve(prob, Rodas5P(); reltol = 1e-13, abstol = 1e-15).u[end]

nsa(alg) = NonlinearSolveAlg(alg)
s_lm = solve(prob, FBDF(nlsolve = nsa(LevenbergMarquardt(; autodiff = AutoForwardDiff())));
             reltol = 1e-6, abstol = 1e-9)
s_nr = solve(prob, FBDF(nlsolve = nsa(NewtonRaphson(; autodiff = AutoForwardDiff())));
             reltol = 1e-6, abstol = 1e-9)

s_lm.retcode                          # ReturnCode.Success
maximum(abs.(s_lm.u[end] .- ref))     # 1.4069e+03
maximum(abs.(s_nr.u[end] .- ref))     # 4.2322e-02

OREGO is the better reproduction: same silent-wrong outcome, far cheaper than integrating Robertson to 1e11.

Also correcting two claims in the analysis above

• The η·ndz < κ branch is not the whole story. Shrinking κ does not fix it — κ = 1e-2 → 1e-6 moves the error only 5.03e7 → 4.80e7 (naccept 31609 → 418814), because the failure migrates to the κ-independent iter == 1 && ndz < 1e-5 branch. At κ=1e-6, 3317 of 9353 accepted stages go through that branch. Any fix confined to the η branch is provably insufficient, measured.
• θ does not freeze for the rest of the solve. It takes 2069 distinct values over 3749 steps; it stalls only in windows, and 862 stages take ≥2 iterations.

The linear-model defect ‖W·δ + fu‖ / ‖fu‖ measured at accepted stages remains a clean discriminator — LevenbergMarquardt [1.96e-2, 1.04] vs NewtonRaphson [1.44e-14, 1.22e-8] and TrustRegion [8.85e-14, 1.63e-8] — but gating on it only changes the outcome if both acceptance branches are gated, and the result is then a loud MaxIters rather than a working solve.