Fix nested ForwardDiff over NonlinearLeastSquaresProblem default vjp

[ChrisRackauckas-Claude]

Draft — please ignore until reviewed by @ChrisRackauckas.

Summary

Fixes the MethodError: no method matching Float64(::ForwardDiff.Dual{...}) that triggers in CI for the core test group on master, e.g.

• lib/SimpleNonlinearSolve/test/core/forward_diff_tests.jl:90 — ForwardDiff.jl Integration NonlinearLeastSquaresProblem
• test/forward_ad_tests.jl:128 — NLLS Hessian SciML/NonlinearSolve.jl#445

These showed up unrelated to the PR that surfaced them (e.g. https://github.com/SciML/NonlinearSolve.jl/actions/runs/25295011499/job/74152282985, https://github.com/SciML/NonlinearSolve.jl/actions/runs/25295011493/job/74152283034) — they regress on master.

Root cause

NonlinearSolveBase.nlls_generate_vjp_function (no user jac/vjp branch) uses DI.pullback with AutoForwardDiff() for small problems. Under DifferentiationInterface, AutoForwardDiff's pullback is emulated via pushforward, and the emulation writes the cotangent into a buffer keyed off eltype(u).

When the returned closure is itself called under an outer ForwardDiff layer — e.g. ForwardDiff.gradient(p -> sum(abs2, solve(NLLS(...; p)).u), p) or any ForwardDiff.hessian of a solve — the inner solve unwraps p to Float64 and computes uu::Vector{Float64}, but the closure runs under nonlinearsolve_∂f_∂p's ForwardDiff.jacobian(Base.Fix1(vjp_fn, uu), p), where p is reseeded as Vector{Dual}. So inside the closure:

• u::Vector{Float64}
• p::Vector{Dual{tag_outer, ...}}
• Base.Fix2(raw_f, p)(u) outputs Vector{Dual{tag_outer, ...}}

DI's pushforward emulation then tries to splat the Dual cotangent values into a Vector{Float64} output buffer (arroftup_to_tupofarr in DI's linalg.jl), and setindex! rejects them.

Fix

When the inner backend is AutoForwardDiff, materialize the Jacobian via DI.jacobian and form 2 * J' * resid directly. DI.jacobian with AutoForwardDiff dispatches to ForwardDiff.jacobian, which handles nested Duals natively using fresh tags per layer.

For the IIP branch, the raw in-place function is wrapped as an OOP closure so that ForwardDiff.jacobian allocates an output buffer of the correctly-promoted Dual type. Using the raw IIP form would force the Dual nesting order of a pre-allocated output buffer, which doesn't compose under nested tags.

Reverse-mode backends — used for larger NLLS problems via select_reverse_mode_autodiff — keep going through DI.pullback (they don't have this issue, since reverse-mode AD has a real pullback rather than pushforward emulation).

Bumps NonlinearSolveBase to 2.25.1.

Test plan

• Local repro of forward_diff_tests.jl:90's OOP NLLS test passes after fix (gradient computed correctly)
• Local repro of forward_ad_tests.jl:128's NLLS Hessian test passes for both the IIP gradient and the IIP Hessian after fix (matches FiniteDiff.finite_difference_gradient / finite_difference_hessian to atol=1e-3)
• Full SimpleNonlinearSolve core group passes locally on Julia 1.11 (35706 Pass, 4 Broken, 35710 Total | 11m11.3s)
• Runic.jl --check clean

🤖 Generated with Claude Code