Add adaptive Jacobian reuse to first-order solvers

[ChrisRackauckas-Claude]

Ignore this PR until it has been reviewed by @ChrisRackauckas.

What

Adds a JacobianReuse(; max_age = 10, max_residual_ratio = 0.1) policy to the first-order solvers (NewtonRaphson, TrustRegion, GaussNewton, LevenbergMarquardt, PseudoTransient, and the first-order polyalgorithms) via a jacobian_reuse keyword. While the residual keeps contracting fast enough and the Jacobian is younger than max_age accepted steps, the Jacobian — and, for an unchanged concrete linear system, its factorization — is reused. A rejected trust-region step, or a failed line search or linear solve on a stale Jacobian, requests a fresh one.

The default is on for length(u0) ≥ 16 and off below it. jacobian_reuse = false forces exact Newton steps and reproduces the previous behavior exactly; jacobian_reuse = JacobianReuse() forces reuse on at any size.

Rebased onto current master (38b61b13).

Design: the policy is a value, not a type

jacobian_reuse is a field of the @concrete algorithm struct, so a policy chosen from a runtime property of the problem would give NewtonRaphson() two types. Measured on the earlier revision of this branch, which encoded "off" as nothing:

NewtonRaphson()                                    → concrete alg / cache / solution
NewtonRaphson(jacobian_reuse = true)               → concrete alg / cache / solution
NewtonRaphson(jacobian_reuse = length(u0) > 25)    → Union{...} at all three levels
                                                     @inferred solve(...) FAILS
                                     inferred cache type: 5131 → 10386 characters

The whole policy now fits in one concrete type, and max_age carries the decision:

JACOBIAN_REUSE_AUTO = -1   # resolve against length(u0) when the cache is built
0 (or 1)                   # reuse off, exact Newton steps
n ≥ 2                      # a Jacobian may serve n accepted steps

nothing, false, true and an explicit JacobianReuse all normalize to JacobianReuse{Float64} before reaching the algorithm struct, so NewtonRaphson(), NewtonRaphson(jacobian_reuse = false) and NewtonRaphson(jacobian_reuse = JacobianReuse()) are one type; resolve_jacobian_reuse then varies only the Int. This mirrors FastShortcutNonlinearPolyalg, where u0_len already picks a start_index::Int rather than a type. A matrix-free Jacobian gets the same cache with the policy switched off rather than no cache, so reset_jacobian_reuse!, jacobian_is_stale and prepare_next_jacobian! each have exactly one method.

Two guards make max_age = 0 exactly equal to the old nothing rather than approximately: reuses_jacobian gates the stale-Jacobian tests that a bare age > 0 would now answer differently on the manual step!(cache; recompute_jacobian = false) path, and the reuse cache skips its residual norm when the policy is off. The max_age = 0 reproduces exact Newton testset asserts identical njacs, nfactors, nsteps and u.

The one leaky part of the encoding is that a negative sentinel shares a field with a count, so validation is max_age >= JACOBIAN_REUSE_AUTO rather than max_age >= 0. The alternative — a separate Bool field — reads worse at the use sites.

Where the cutoff of 16 comes from

Speedup of the default policy over jacobian_reuse = false, min-of-N wall time, three problem families × two algorithms:

n	bratu1d/NR	bratu1d/TR	dense_ad/NR	dense_ad/TR	analyticJ/NR	analyticJ/TR	min
2	1.011	0.993	1.080	0.833	0.996	0.901	0.833
4	1.002	0.998	1.047	1.050	0.951	0.968	0.951
8	1.000	0.986	1.271	1.049	0.972	0.895	0.895
10	0.997	0.989	1.249	1.269	0.948	1.028	0.948
12	1.030	1.019	1.363	1.259	1.057	1.037	1.019
14	1.026	1.019	1.356	1.245	0.963	1.014	0.963
16	1.025	1.021	1.368	1.318	1.061	0.997	0.997
24	1.038	1.020	1.430	1.388	1.139	1.071	1.020
32	1.041	1.029	1.457	1.582	1.196	1.071	1.029
64	1.063	1.049	1.921	1.375	1.436	1.225	1.049

dense_ad is a dense system with an AD Jacobian (C_J ≈ n residual evaluations), analyticJ is the same system with f.jac supplied, bratu1d is a sparse tridiagonal Bratu problem — the cheapest Jacobian of the three. 16 is the smallest size at which nothing measured is below baseline while the expensive-Jacobian family is already at 1.32–1.37x. Below it the cheap-Jacobian families run at 0.90–1.06x.

This is a noise-band crossing, not a sharp one: from n = 12 to n = 24 the cheap-Jacobian families move about ±4% run to run.

Why max_residual_ratio dropped from 1 to 0.1

max_residual_ratio = 1 reuses while the residual improves at all, which is almost always true, so it runs to max_age every time. Over 32 problem-and-algorithm cases at or above the cutoff (dense AD, dense analytic, expensive residual, Brusselator 2D sparse at n = 128…8192, Bratu 1D sparse at n = 100…5000), max_age = 10 throughout:

max_residual_ratio	geomean	worst case
1	1.508	0.767
0.5	1.353	0.713
0.25	1.386	0.764
0.1 (default)	1.383	1.016
0.05	1.280	1.005

A permissive ratio has the better mean and the worse tail — 0.767x on the sparse Brusselator at n = 8192, and with max_age = 100 it turns a Success into a Stalled at n = 512. 0.1 gives up almost none of the mean and never fell below 1.016x.

Other behavior to review

• step! no longer takes the stale-Jacobian retry when the caller passed recompute_jacobian explicitly. OrdinaryDiffEqNonlinearSolve always does (recompute_jacobian = nlsolver.iter == 1 && (cache.W === nothing || cache.new_W)), so turning the default on cannot change how an ODE solver's Newton iteration behaves on a failed line search.
• step!(cache; recompute_jacobian = true) now always recomputes. Previously it also required the cache to want a new Jacobian. This is what the step! docstring in NonlinearSolveBase already promises.
• Matrix-free Jacobians ignore the policy. JacobianCache{<:JacobianOperator} rebinds a StatefulJacobianOperator to the current iterate on every call, so nothing is cached and nothing goes stale; the reuse cache is not built for them.
• Deferred residuals (evaluate_residual = false, from Let a driver defer the residual evaluation that ends a step #1167) drive the policy from refresh_residual!, where the residual at the new iterate actually becomes available.
• Versions: NonlinearSolveFirstOrder 2.4.1 → 2.5.0 (new export plus a defaults change). Root 4.28.1 → 4.29.0 with NonlinearSolveFirstOrder = "2.5" compat — the root polyalgorithm passes jacobian_reuse through unconditionally, so a root resolved against FirstOrder 2.4.x would throw from FastShortcutNonlinearPolyalg().
• The line-search stale-Jacobian retry logs under the existing :linsolve_failed_noncurrent verbosity key rather than adding a new one.

Verification (Julia 1.12, Linux x86_64)

Against master at c498eaff, jacobian_reuse = false on this branch reproduces master exactly — identical nf, njacs and nsteps on all 22 (problem × algorithm) pairs (scalar, SVector, dense n = 2…200, sparse Bratu, Brusselator; NewtonRaphson and TrustRegion), with wall times within noise on a loaded machine.

lib/NonlinearSolveFirstOrder  GROUP=Core
  Adaptive Jacobian reuse      |  125    125   15.8s
  Deferred residual evaluation |  144    144    3.4s
  NewtonRaphson                |  309    309  6m57.8s
  TrustRegion                  | 1176   1176  17m23.0s
  General NLLS Solvers         |  480    480  3m40.6s
  PseudoTransient              |   51     51  1m17.2s
  LevenbergMarquardt           |   33     33  1m00.4s
  Brusselator 2D               |    5      5   18.0s
  Structured Jacobians         |   16     16   26.0s
  SciMLOperator Jacobians      |   14     14    7.1s
     Testing NonlinearSolveFirstOrder tests passed     (exit 0, no failures/errors)

lib/NonlinearSolveFirstOrder  GROUP=QA                  28/28, exit 0

root  NONLINEARSOLVE_TEST_GROUP=PolyAlgorithms
  Basic PolyAlgorithms             |  24   24
  23 Test Problems: PolyAlgorithms |  46   46
     Testing NonlinearSolve tests passed                (exit 0)

root  NONLINEARSOLVE_TEST_GROUP=Core
  PolyAlgorithm Type Inference | 16  16
  Allocation-free cache solve  | 22  22
     Testing NonlinearSolve tests passed                (exit 0)

docs  julia --project=docs docs/make.jl                 exit 0, no errors.
  Warnings are the pre-existing size_threshold_warn on native/solvers.md and
  devdocs/internal_interfaces.md, two oversized @example blocks on an unrelated
  tutorial page, and skipped deployment.

Runic 1.5.1 --check: clean on every changed .jl file.  typos: clean.  git diff --check: clean.

Not verified

• GPU/CUDA groups (no GPU here), downstream packages, Julia pre, and the wrapper groups that need external solvers.
• The docs build needed JULIA_CONDAPKG_BACKEND=Null with the system Python: SciMLBasePythonCallExt cannot precompile on this machine because CondaPkg has no PyPI access, and docs/make.jl loads NonlinearSolveSciPy. That failure is environmental and predates this branch.
• The root Core and PolyAlgorithms runs above predate the last two commits; lib/NonlinearSolveFirstOrder GROUP=Core and GROUP=QA were re-run on the pushed tip. lib/NonlinearSolveFirstOrder GROUP=Core and GROUP=QA were re-run on the pushed commit f3c9c38a and pass there too (Adaptive Jacobian reuse 125/125, TrustRegion 1176/1176, General NLLS 480/480, QA 28/28, exit 0). The three commits the rebase brought in touch only the NLsolveJL extension and two sublibrary version bumps.

What a reviewer should push back on

• The default now changes which root is returned on some multi-root problems. Over the 23 NonlinearProblemLibrary cases, a permissive policy landed on a different solution than exact Newton on 4 of 23; at max_residual_ratio = 0.1 that is 1 of 23. All of those problems are below the cutoff so the shipped default does not hit them, but a user above the cutoff on a multi-root problem can get a different (still valid) root than before. This is a defaults change with a visible output change, shipping under a minor bump.
• length(u0) is knowingly a proxy for the wrong thing. What decides the payoff is the cost of a Jacobian relative to the cost of a nonlinear step. A 10,000-unknown sparse tridiagonal system gains ~1.0x; a 25-unknown dense system with an AD Jacobian gains 1.9x. Both sit above the cutoff. Follow-up with the measured cost model and the static predictors that could replace it: JacobianReuse default should key off Jacobian cost, not length(u0) #1216
• Matrix-free configurations get nothing. With linsolve = KrylovJL_GMRES() and no concrete_jac, njacs == 0 and the policy is inert; every setting timed within 0.2% of baseline at n = 512 and n = 2048. The size default enables a policy that cannot act there.
• The cutoff and the ratio were both picked from a handful of problem families on one loaded machine.

🤖 Generated with Claude Code (model: claude-opus-5[1m])
https://claude.ai/code/session_01CKHupdFDQQCFcEgqThh6Mz

[ChrisRackauckas-Claude]

Implementation scratchpad for commit 6dd21578d:

• preserve exact-Newton defaults and OrdinaryDiffEq's explicit J/W refresh ownership
• add a stateless public policy plus per-cache residual/age tracking
• refresh on bounded age, lost residual progress, or stale solve/globalization failures
• preserve current Jacobians across unchanged-state trust-region rejections
• forward the opt-in policy through first-order constructors and polyalgorithms
• validate with Runic, Julia 1.10/1.12 targeted tests, FirstOrder Core/QA, root PolyAlgorithms, and FBDF/TRBDF2 integration
• retain strict docs link checking; the local build reached document checks and then hit transient failures on pre-existing external links (timeout/HTTP 429)

Benchmarking showed wins for expensive dense/sparse factorizations but overhead and policy sensitivity on small/problem-library cases, so the solver default remains exact Newton.

[ChrisRackauckas-Claude]

Published final stacked head 99348eb0d.

Safety and history checks:

• fetched upstream, fork, and origin immediately before push
• upstream remained at 6264446f4
• remote feature branch remained at the original single-author head 6dd21578d; no contributor commits intervened
• git range-diff reports 6dd21578d = 99348eb0d for the feature patch
• updated with an exact --force-with-lease against 6dd21578d
• every stacked commit retains the required Chris Rackauckas co-author trailer
• invalid Brown expectation commit is absent

Final local outputs include:

• NonlinearSolveFirstOrder Core: passed (feature 64/64, Newton 309/309, TrustRegion 1176/1176, remaining groups passed)
• NonlinearSolveFirstOrder QA: 18/18
• Julia 1.10 focused feature suite after clean resolve: 64/64
• root PolyAlgorithms: passed, including 46/46 problem-library assertions
• root QA: 17/17
• NonlinearSolveBase Core: passed
• NonlinearSolveBase QA: 18/18
• continuation inference: 26/26
• Runic and git diff --check: passed

CI is now running on the combined stack, which includes both focused prerequisite fixes before the unchanged feature commit.

[ChrisRackauckas-Claude]

Final-stack update at 63e043babea0cf86da82009c1b4c918a4ceeeb3e:

• Before replacing the old 99348eb0d remote head, I fetched and audited the exact remote branch history, PR commits, comments, and reviews. There were no intervening contributor commits or reviews. The update used an exact --force-with-lease against that audited head.
• The branch was rebuilt from upstream master 60736b5b6 with the focused prerequisites Eliminate dense continuation step allocations on Julia 1.12 #1061, Use public cache construction for ForwardDiff wrappers #1073, Preserve triangular solves for Broyden Jacobian inversion #1057, Restore root LinearSolve 5 floor #1087, Remove duplicate BibTeX entry #1088, and Restore GPU Broyden inverse initialization #1066. Range-diff checks confirmed that the Jacobian-reuse feature patch is unchanged.
• Exact-head local validation passed: NonlinearSolveBase Core and QA (18/18 QA), NonlinearSolveQuasiNewton Core, NonlinearSolveFirstOrder QA (18/18), root PolyAlgorithms (107/107), root QA (17/17), and Runic 1.7.0.
• The integrated FirstOrder Core feature run passed all groups, including the 64/64 adaptive-reuse regressions. The exact FWW 1.9.3 floor passed Base Core on Julia 1.10.11. The focused V100 CUDA prerequisite passed 22/22, and the strict documentation prerequisite passed 11/11.

The PR description now records the actual stack, benchmark/default decision, baseline-failure investigations, and observed validation. Refreshed combined-stack CI is running.

[ChrisRackauckas-Claude]

Published rebuilt combined head ed4d653e16a378505214b1f601028fda91ad90d7 with the Apple ARM allocation follow-up placed directly behind #1061.

Before replacing the old 63e043babea0cf86da82009c1b4c918a4ceeeb3e remote head, I fetched origin, upstream, and the fork, then audited the exact branch history plus PR commits, comments, and reviews. The PR still had the audited head, no reviews, and no comments from other contributors. The update used an exact --force-with-lease against that commit.

git range-diff marks every commit unchanged; only the new follow-up moved beside its prerequisite. git diff --exit-code confirmed that the reordered head is tree-identical to the locally tested tree. All 11 commits retain the required Chris Rackauckas co-author trailer.

New local evidence includes:

• NonlinearSolveBase Core passed on Julia 1.12 and 1.10.11, including allocation fast paths 6/6
• NonlinearSolveBase QA: 18/18
• unchanged root continuation allocation regression: sweep/arclength 0.0/0.0 B per step
• root PolyAlgorithms after the follow-up: 107/107
• Runic 1.7.0 over the repository and git diff --check: passed

The superseded combined CI run finished 90 successful, 3 skipped, and 4 failed. Two Ubuntu failures were the already-audited Brown unexpected-pass behavior. The two macOS failures were the allocation issue addressed by the cached destination. Fresh CI is now running; Apple ARM must pass before the cross-platform fix is claimed.

[ChrisRackauckas-Claude]

Published rebuilt combined head 5b9e07da11ab32181e86824ef0c33051f3989a29 with the cache-only corrector commit placed directly behind #1061.

Safety audit immediately before publication:

• fetched origin, upstream, and the fork
• upstream remained 60736b5b6
• the remote feature branch remained exactly ed4d653e16a378505214b1f601028fda91ad90d7
• PR Add adaptive Jacobian reuse to first-order solvers #1072 remained draft, with no reviews or comments from other contributors
• updated only this branch using --force-with-lease pinned to that exact audited SHA

git range-diff marks all 12 patches equal; only the new allocation commit moved beside its prerequisite. The new head is tree-identical to the fully tested pre-reorder tree, and all commits have the required co-author trailer.

The prior Apple prerelease job is important negative evidence: it still measured exactly 772.0/400.0 B per sweep/arclength step after the cached Jacobian destination. The remaining allocation came from constructing a complete solution for every interior corrector. The new helper runs those correctors to completion in the reused cache; ordinary solve results, anchors, landings, and failure originals are unchanged.

Final local evidence includes cache-only 0 B with same-cache identity, end-to-end 0.0/0.0 B per step, Base Core on Julia 1.12 and 1.10, all 39 continuation files, Base QA 18/18, ForwardDiff cache reinit 1/1, Runic 1.7, and git diff --check. Fresh Apple ARM CI is required before claiming the cross-platform result.

[ChrisRackauckas-Claude]

Final-tree local integration update for 5b9e07da11ab32181e86824ef0c33051f3989a29:

• GROUP=Core timeout 3600 julialauncher +1.12 --project=lib/NonlinearSolveFirstOrder -e 'using Pkg; Pkg.test()' passed.
• Key totals observed in this run: adaptive reuse 64/64, general NLLS 472/472, NewtonRaphson 309/309, PseudoTransient 51/51, TrustRegion 1176/1176, Levenberg–Marquardt 33/33, Eisenstat–Walker 27/27, sparse autodiff 12/12, structured Jacobians 16/16, and SciMLOperator Jacobians 14/14.
• The generated ignored lib/NonlinearSolveFirstOrder/Manifest.toml was removed; both worktrees are clean.

Fresh Apple ARM Core checks remain queued and are still required for the cross-architecture allocation claim.

[ChrisRackauckas-Claude]

Apple ARM prerelease result on feature head 5b9e07da1:

• tests / Core (julia pre, macos-latest) failed at the unchanged strict continuation slopes: sweep 772.0 B/step, arclength 400.0 B/step.
• The new direct cache-only solve/identity check and preceding continuation tests passed. This disproves the proposed attribution to per-corrector NonlinearSolution construction; the remaining bytes are outside _solve_without_solution!.
• I have not changed the < 48 thresholds or silenced the test.
• Test-only diagnostic commit ec549add5 adds Profile.Allocs aggregation by allocation type and first NonlinearSolve stack site only when a slope fails. It passed locally, including profiler smoke output, and is now triggering a fresh Apple run.

[ChrisRackauckas-Claude]

Diagnostic rerun correction at cfb8e2d18:

• The first diagnostic CI attempt did not reach profiling on Julia 1.12 because Pkg.test() isolates undeclared test dependencies and Profile was absent from the root test target.
• Profile is now declared in the existing root [extras] and targets.test; no test/Project.toml was created.
• The exact allocation test passes locally on Julia 1.13.0-rc1 (exit 0), and TOML parsing confirms both declarations.
• The unchanged strict < 48 B/step assertions remain in place. A fresh Apple ARM run is now queued to produce the allocation-type/site profile.

[ChrisRackauckas-Claude]

Apple allocation profile root cause and upstream follow-up:

• Stable macOS measured sweep/arclength slopes of 957.28 / 496 B per step with the unchanged < 48 assertions.
• LU{Float64, Matrix, Vector{Int32}} construction in the Apple Accelerate LinearSolve path contributed 25,984 bytes / 812 allocations for sweep and 11,648 bytes / 364 allocations for arclength. The profile also showed repeated tuple, Ref, range, and keyword-NamedTuple construction.
• This maps to AppleAccelerateLUFactorization rebuilding its LU wrapper and LAPACK workspace representation on each fresh refactorization, not to the NonlinearSolve cache-only result path.
• Draft upstream fix: SciML/LinearSolve.jl#1099. It keeps factors, pivots, and LAPACK status in a private mutable cache and adds macOS LTS/stable/pre allocation tests.
• No threshold was changed or suppressed. The cross-platform downstream claim remains pending LinearSolve macOS CI and a NonlinearSolve rerun against the fixed dependency.

[ChrisRackauckas-Claude]

Upstream Apple-silicon validation update for LinearSolve #1099:

• Julia 1.10: refactorization reuse 55/55; mixed precision 19 passed with 1 pre-existing broken expectation
• Julia 1.12: refactorization reuse 59/59; mixed precision 19 passed with 1 pre-existing broken expectation
• Julia 1.13: refactorization reuse 59/59; mixed precision 19 passed with 1 pre-existing broken expectation

These are direct LinearSolve macOS tests of the fixed Apple Accelerate path. They do not replace the required NonlinearSolve end-to-end continuation rerun. The combined PR therefore still makes no downstream cross-platform allocation claim until #1099 is merged/released and this branch is retested against that dependency.

[ChrisRackauckas-Claude]

Fresh exact-graph validation of the stacked #1057 prerequisite:

• Baseline NonlinearSolve 60736b5b6 against LinearSolve 0496f135 on Julia 1.10.11 reproduced root Broyden at 94 passed / 1 failed / 1 errored / 19 broken: Generalized Rosenbrock alg Make iip a compile-time constant #4 returned residual 4.3999999999999995, and Brown alg Make iip a compile-time constant #4 unexpectedly passed its existing broken expectation.
• Applying the Preserve triangular solves for Broyden Jacobian inversion #1057 patch unchanged made NonlinearSolveBase Core pass (including workspace 53/53), NonlinearSolveQuasiNewton Core pass (Broyden 810/810, Klement 297/297, limited-memory Broyden 135/135, dedicated Rosenbrock regression 2/2), and changed root Broyden to 95 passed / 0 failed / 1 errored / 19 broken.
• Generalized Rosenbrock therefore passes with the exact Preserve triangular solves for Broyden Jacobian inversion #1057 code already included in this combined stack. The only remaining root error is the separately tracked platform-sensitive Brown unexpected pass in #1056.
• Runic and git diff --check passed; temporary validation patches and generated files were removed.

[ChrisRackauckas-Claude]

Allocation regression update: the GenericLU workspace fix has been reviewed and merged upstream in LinearSolve at SciML/LinearSolve.jl#1195. The blocked-LU packing buffer and pivot storage are now allocated during cache initialization/resize and reused by solve!, with AllocCheck.jl plus runtime zero-allocation QA coverage. On the exact local NonlinearSolve sweep, the before value was 61.76 bytes/step and the fixed LinearSolve checkout measured sweep_per_step=0.0. LinearSolve main has been bumped to 5.9.0, but no v5.9.0 tag/release is published yet, so these existing checks cannot pick up the fix until that release is available or the workflow is rerun against the merged commit.

[ChrisRackauckas]

Claude Code review — 2 finding(s)

normal — lib/NonlinearSolveFirstOrder/src/solve.jl:376-380

The new line-search stale-Jacobian retry at lib/NonlinearSolveFirstOrder/src/solve.jl:381 does InternalAPI.step!(cache; recompute_jacobian = true, cache.kwargs...), but the step! method for GeneralizedFirstOrderAlgorithmCache only accepts the single recompute_jacobian keyword — so splatting any non-empty cache.kwargs throws MethodError. Under a realistic combination (jacobian_reuse = true with a linesearch, plus any user-passed extra kwarg such as alias_u0), the retry crashes the solve. Fix is to drop , cache.kwargs... to match the linear-solve retry a few lines above (line 358).

nit — lib/NonlinearSolveBase/test/linsolve_workspace.jl:66-71

The new testset at lib/NonlinearSolveBase/test/linsolve_workspace.jl:66-71 ("arrays without fast scalar indexing use pinv") is a verbatim duplicate of the pre-existing testset with the same name and identical body at lines 42-47. Remove the duplicate — it wastes CI time and clutters the file, though it does not affect correctness.

---

Generated by Claude Code

[ChrisRackauckas-Claude]

Review findings addressed (2)

Ran an automated deep review over this PR; two findings, both fixed on this branch.

1. step! retry splats cache.kwargs into a method that rejects them (bug)

lib/NonlinearSolveFirstOrder/src/solve.jl:381 — the new line-search stale-Jacobian retry called

InternalAPI.step!(cache; recompute_jacobian = true, cache.kwargs...)

but the step! method for GeneralizedFirstOrderAlgorithmCache (line 311) declares only
recompute_jacobian and has no kwargs... catch-all, so any non-empty cache.kwargs throws.
cache.kwargs holds every keyword __init did not consume by name — alias_u0 is read via
haskey and left in place, and user options like callback/saveat/progress land there too.
So solve(prob, NewtonRaphson(jacobian_reuse = true, linesearch = ...); alias_u0 = false) aborts
the solve the first time the line search rejects a step taken on a stale Jacobian.

The pre-existing linear-solve retry at line 358 already gets this right; the fix is to match it:

InternalAPI.step!(cache; recompute_jacobian = true)

The existing misc_tests__item7.jl testset could not catch this: abstol, reltol and verbose
are all consumed by name in __init, leaving cache.kwargs empty so the splat expanded to
nothing. The test now also passes alias_u0 = false, which is retained in cache.kwargs.

Failing before the fix (test change applied, solve.jl:381 left as-is on this branch):

Line search stale-Jacobian retry: Error During Test at lib/NonlinearSolveFirstOrder/test/misc_tests__item7.jl:127
  Got exception outside of a @test
  MethodError: no method matching step!(::NonlinearSolveFirstOrder.GeneralizedFirstOrderAlgorithmCache{...})
  This method does not support all of the given keyword arguments (and may not support any).

  Closest candidates are:
    step!(::NonlinearSolveFirstOrder.GeneralizedFirstOrderAlgorithmCache; recompute_jacobian) got unsupported keyword argument "alias_u0"
     @ NonlinearSolveFirstOrder lib/NonlinearSolveFirstOrder/src/solve.jl:311
  ...
    [2] step!(cache::...GeneralizedFirstOrderAlgorithmCache{...}; recompute_jacobian::Nothing)
      @ NonlinearSolveFirstOrder lib/NonlinearSolveFirstOrder/src/solve.jl:381

Test Summary:                    | Pass  Error  Total  Time
Line search stale-Jacobian retry |    2      1      3  3.7s
ERROR: LoadError: Some tests did not pass: 2 passed, 0 failed, 1 errored, 0 broken.

Passing after the fix (same file, same command):

Test Summary: | Pass  Total  Time
configuration |   18     18  0.8s
Test Summary:         | Pass  Total  Time
Newton refresh policy |   15     15  1.8s
Test Summary:                    | Pass  Total  Time
Line search stale-Jacobian retry |    6      6  0.6s
Test Summary:              | Pass  Total  Time
manual override and reinit |   11     11  0.2s
Test Summary:     | Pass  Total  Time
TrustRegion reuse |   14     14  1.2s

(2 pass + 1 error becomes 6 pass: with the retry no longer crashing, the assertions after it now run.)

2. Duplicate testset (nit)

lib/NonlinearSolveBase/test/linsolve_workspace.jl carried two byte-identical copies of
@testset "arrays without fast scalar indexing use pinv" (lines 42-47 and 66-71). Removed the
second copy.

Verification

$ GROUP=Core julia --project=. -e 'using Pkg; Pkg.test()'   # lib/NonlinearSolveBase
linsolve_identity!! workspace (#1020) |   55     55  18.3s
     Testing NonlinearSolveBase tests passed

lib/NonlinearSolveFirstOrder GROUP=Core result is posted in a follow-up comment.

Formatting/spelling: Runic.jl 1.7.0 --check clean on all three touched files; typos clean.

Not verified: GPU/CUDA groups (no GPU available here), downstream packages, and the docs
build (no docstrings or docs/ content touched).

🤖 Generated with Claude Code (model: claude-opus-5[1m])
Session: https://claude.ai/code/session_013vbihgEfx7VCEa4W3ccccG

[ChrisRackauckas-Claude]

Pushed as a9849a5.

lib/NonlinearSolveFirstOrder GROUP=Core on the branch with both fixes:

$ GROUP=Core julia --project=. -e 'using Pkg; Pkg.test()'
Test Summary:           | Pass  Total  Time
Adaptive Jacobian reuse |   64     64  5.4s
...
Test Summary:            | Pass  Total  Time
SciMLOperator Jacobians  |   14     14  6.1s
     Testing NonlinearSolveFirstOrder tests passed

Zero errors, zero failures across the group. lib/NonlinearSolveBase GROUP=Core also passed
(reported in the previous comment).

Base branch check: git rev-list --count HEAD..origin/codex/jacobian-reuse-heuristics was 0
before pushing, so nothing was clobbered.

🤖 Generated with Claude Code (model: claude-opus-5[1m])
Session: https://claude.ai/code/session_013vbihgEfx7VCEa4W3ccccG

[ChrisRackauckas-Claude]

Force-pushed a rebase onto master c498eaff (previous tip a9849a54, kept locally as pr-1072-old). The PR body has been rewritten to describe the current tree.

What changed in the rebase

• Dropped the "prerequisite stack": Preserve triangular solves for Broyden Jacobian inversion #1057 (closed), Use public cache construction for ForwardDiff wrappers #1073 (closed; its successor Use public APIs for AutoDePSpecialize wrappers #1121 also closed), the Eliminate dense continuation step allocations on Julia 1.12 #1061 allocation-profiling scaffolding and its Profile test dep, and Restore GPU Broyden inverse initialization #1066 (already merged). 22 files / +603 → 17 files / +503, all on topic.
• misc_tests__item7.jl collided with the deferred-residual tests master gained in Let a driver defer the residual evaluation that ends a step #1167; the reuse tests are now misc_tests__item9.jl.
• Three real fixes on top of the feature, each with a test that fails without it:
  • step!(...; evaluate_residual = false) left the policy reading a stale residual and permanently forcing refreshes — the policy now runs from refresh_residual!.
  • Matrix-free Jacobians no longer get a reuse cache: StatefulJacobianOperator is rebound every step, so "reuse" was a no-op that still triggered a pointless line-search retry.
  • Versions/compat: FirstOrder 2.5.0, root 4.29.0 with NonlinearSolveFirstOrder = "2.5". Without the compat floor a released root against FirstOrder 2.4.x throws unsupported keyword argument from the default FastShortcutNonlinearPolyalg().
• step!(cache; recompute_jacobian = true) now always recomputes (matches the documented contract; OrdinaryDiffEqNonlinearSolve unaffected) — flagged in the body as a behavior change.

Verification on the pushed tree (Julia 1.12.4): FirstOrder GROUP=Core exit 0 (Adaptive Jacobian reuse 70/70, Deferred residual 144/144, TrustRegion 1176/1176, NewtonRaphson 309/309), FirstOrder GROUP=QA 20/20, root GROUP=PolyAlgorithms 46/46, docs/make.jl exit 0, Runic/typos/git diff --check clean. Not run: GPU groups, downstream, root GROUP=Core.

🤖 Generated with Claude Code (model: claude-opus-5[1m])
https://claude.ai/code/session_013vbihgEfx7VCEa4W3ccccG

[ChrisRackauckas-Claude]

CI on 3663b743: 89 pass, 6 fail, 1 skipped. The 6 failures are exactly the six tests / Wrappers jobs (lts/1/pre × ubuntu/macOS), which fail identically on master (https://github.com/SciML/NonlinearSolve.jl/actions/runs/33251921454): NLsolve 5.0 changed autodiff from Symbol to an ADTypes backend and ext/NonlinearSolveNLsolveExt.jl still passes :central. Not related to this PR; fix is in #1212.

Everything this PR touches is green: all NonlinearSolveFirstOrder sublibrary groups on lts/1/pre, NonlinearSolveBase, root Core/PolyAlgorithms/Downgrade, downgrade-sublibraries, Documentation, Runic, typos, and the OrdinaryDiffEq/DiffEqCallbacks downstream jobs.

🤖 Generated with Claude Code (model: claude-opus-5[1m])
https://claude.ai/code/session_013vbihgEfx7VCEa4W3ccccG

[ChrisRackauckas-Claude]

Two commits since the last revision, both about the default rather than the mechanism.

9614792f Make Jacobian reuse a value rather than a type. The policy lived in a type parameter of the algorithm struct, so any default that consulted the problem would have split NewtonRaphson() into two types — @inferred solve fails and the inferred cache type goes from 5131 to 10386 characters. "Off" is now max_age = 1, so the whole policy is one concrete type and the size decision moves only an Int. No behavior change in this commit.

f3c9c38a Default Jacobian reuse on above a length(u0) cutoff. Cutoff 16, and max_residual_ratio 1 → 0.1. Both numbers, and the reasons the previous ratio could not be defaulted on, are in the body. Also: the stale-Jacobian retry is now skipped when the caller passed recompute_jacobian explicitly, which is what keeps OrdinaryDiffEqNonlinearSolve unaffected by the new default.

Rebased onto 38b61b13. GROUP=Core and GROUP=QA for NonlinearSolveFirstOrder pass on the pushed commit; root Core and PolyAlgorithms pass on the identical pre-rebase tree.

Follow-up on replacing the size cutoff with a cost-based heuristic: #1216

[ChrisRackauckas-Claude]

a448611d — review catch: the previous commit left two type-level distinctions behind.

jacobian_reuse = nothing still reached the algorithm struct as Nothing, so NewtonRaphson() and NewtonRaphson(jacobian_reuse = false) were different types; and a matrix-free Jacobian got a nothing reuse cache, so reset_jacobian_reuse!, jacobian_is_stale and prepare_next_jacobian! each carried a second method for it.

Both are gone. nothing normalizes to JacobianReuse(JACOBIAN_REUSE_AUTO, ...) before it reaches the struct, and matrix-free gets the same cache with the policy switched off. max_age now carries the whole decision:

JACOBIAN_REUSE_AUTO = -1   # resolve against length(u0) when the cache is built
0 (or 1)                   # reuse off, exact Newton steps
n ≥ 2                      # a Jacobian may serve n accepted steps

so every spelling of the keyword is one concrete type and every hot-path function has one method. The testset grew by three assertions that check exactly this and would have failed before the commit:

@test typeof(alg()) === typeof(alg(jacobian_reuse = false)) ===
    typeof(alg(jacobian_reuse = true)) === typeof(alg(jacobian_reuse = JacobianReuse()))

Verified on the pushed tip: Adaptive Jacobian reuse 125/125, TrustRegion 1176/1176, General NLLS Solvers 480/480, Testing NonlinearSolveFirstOrder tests passed; GROUP=QA 28/28. Runic and typos clean.

Worth a look from a reviewer: the encoding puts a negative sentinel in the same field as a count, so validation is max_age >= JACOBIAN_REUSE_AUTO rather than max_age >= 0. The alternative was a separate Bool field, which reads worse at the use sites.

[ChrisRackauckas-Claude]

d7fa8986 — docs fix plus a consistency pass.

The one-type commit put a [JACOBIAN_REUSE_AUTO](@ref) link into the JacobianReuse docstring, which a @docs block renders. That constant has a docstring but no @docs entry, and with checkdocs = :exports and no warnonly that terminates the build. Reproduced on a two-file Documenter project of the same shape:

┌ Error: Cannot resolve @ref for md"[`INTERNAL_KNOB`](@ref)" in docs/src/index.md.
ERROR: LoadError: `makedocs` encountered an error [:cross_references] -- terminating build before rendering.

JACOBIAN_REUSE_SIZE_CUTOFF and JACOBIAN_REUSE_AUTO are tuning constants whose values move with the benchmarks, so the fix describes the behavior in prose with the numbers interpolated rather than promoting them to public API.

Swept the rest of the feature's docs for the same class of problem while I was there:

• GaussNewton documented jacobian_reuse in a prose paragraph while the other four solvers use a ### Keyword Arguments bullet — now identical to them.
• docs/src/native/solvers.md used [JacobianReuse()](@ref), a parenthesized cross-reference; now a plain reference to the type.
• The manual's example only showed JacobianReuse() and false, and never mentioned max_age or that 0 is off. It now shows all three spellings and notes that they produce the same algorithm type.

All five solver kwarg docstrings interpolate the cutoff from the constant, so they cannot drift from it. The paragraph claiming an explicit step!(cache; recompute_jacobian = ...) takes precedence over the policy is now actually true, which it was not before the policy_driven gate in f3c9c38a.

No functional change since a448611d, which passed GROUP=Core and GROUP=QA in full. The docs build and Core on this commit were still running when it was pushed; I will follow up with both results.