Reuse analytic Jacobians for line-search JVPs

[ChrisRackauckas-Claude]

Important

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

Summary

• Preserve explicitly selected JVP, VJP, and reverse-mode line-search backends.
• When no derivative backend was explicitly selected, let SciMLJacobianOperators form the line-search JVP from a supplied analytic Jacobian.
• Add a regression whose residual rejects differentiated inputs.
• Bump NonlinearSolveFirstOrder to 2.3.2 and raise the SciMLJacobianOperators floor to 0.1.17.

Root cause

This follows the direct ModelingToolkitBase InterfaceII failures reported while working on ModelingToolkit.jl#4742. Both failures entered a boundary residual with ForwardDiff.Dual values from the nonlinear line-search JVP.

Commit c6b26f4a began always passing the automatically selected forward-mode JVP backend to LineSearch, which bypasses the analytic-Jacobian fallback in SciMLJacobianOperators.prepare_jvp (that fallback only fires when autodiff === nothing). The introducing boundary was checked against the minimal reproducer: c6b26f4a^ (d21952c3) exits successfully, c6b26f4a errors with Vector{ForwardDiff.Dual} reaching the residual.

The fix passes nothing only when the problem has an analytic Jacobian and no JVP, VJP, or reverse-mode backend was explicitly selected. That preserves the explicit-backend behavior c6b26f4a added while restoring the analytic derivative path.

Prerequisite (resolved)

This depended on a rectangular-output bug in the analytic JVP implementation, split out as #1063. That is now merged (01891bcd) and SciMLJacobianOperators 0.1.17 is registered in General (JuliaRegistries/General#163962), so the 0.1.17 compat floor here resolves — including on downgrade CI.

Verification

Rebased onto master at 01891bcd. Julia 1.12.4, with lib/ path sources so both local SciMLJacobianOperators 0.1.17 and local NonlinearSolveFirstOrder are exercised (asserted via pathof).

The regression discriminates. With lib/NonlinearSolveFirstOrder/src/solve.jl reverted to master and the test left in place:

item4: Error During Test
  Got exception outside of a @test
  LoadError: residual must not be differentiated
  Stacktrace:
   [1] error(s::String)
     @ Base ./error.jl:44
   [2] (::AnalyticJacResidual)(res::Vector{ForwardDiff.Dual{ForwardDiff.Tag{NonlinearSolveBase.NonlinearSolveTag, Float64}, Float64, 1}}, u::Vector{ForwardDiff.Dual{...}}, p::Vector{Float64})
     @ Main .../test/misc_tests__item4.jl:38
   ...
Test Summary: | Pass  Error  Total  Time
item4         |    1      1      2   8.2s
ERROR: LoadError: Some tests did not pass: 1 passed, 0 failed, 1 errored, 0 broken.

With the fix restored, same file, same environment:

Test Summary: | Pass  Total  Time
item4         |    3      3   7.3s

Full GROUP=Core suite for NonlinearSolveFirstOrder — Pkg.test() exited 0, zero failures and zero errors across every testset:

General NLLS Solvers   |  480 pass |  480 total |  3m06.7s
NewtonRaphson          |  309 pass |  309 total |  7m39.3s
PseudoTransient        |   51 pass |   51 total |  1m18.1s
TrustRegion            | 1176 pass | 1176 total | 17m53.0s
LevenbergMarquardt     |   33 pass |   33 total |    58.0s
Custom JVP             |    4 pass |    4 total |     9.1s
Structured Jacobians   |   16 pass |   16 total |    23.6s
SciMLOperator Jacobians|   14 pass |   14 total |     6.2s
Line search uses forward-mode autodiff: Issue #837 | 3 pass | 3 total
     Testing NonlinearSolveFirstOrder tests passed

Structured Jacobians and SciMLOperator Jacobians matter here specifically: the restored prepare_jvp fallback materializes a dense jac_cache and calls f.jac into it, so those groups are what exercise non-dense jac_prototypes through the new path.

Also run clean: Runic 1.7.0 (--check, exit 0), typos (exit 0), git diff --check (exit 0).

Not verified

• QA group, docs build, and the other subpackages — this touches no docstring, no public API, and no other lib/.
• Julia LTS and pre; only 1.12.4 was run locally.
• The original ModelingToolkitBase bvproblem.jl integration was verified on the pre-rebase branch, not re-run after this rebase. The code change is identical; only the Project.toml conflict resolution differs.

For a reviewer to push back on

• SciMLBase.has_jac(f) is the whole trigger condition. A problem that supplies jac but whose jac is more expensive per call than a single forward-mode pushforward now pays that cost in the line search. The prior behavior (before c6b26f4a) was the same, so this is a restoration, not a new tradeoff — but it is a performance-relevant default.
• The compat floor jumps 0.1.1 → 0.1.17, dropping a wide range of resolvable older versions. That is required for correctness (the rectangular-output fix), not cosmetic.

[ChrisRackauckas-Claude]

Investigation scratchpad for commit 19fd7c2b:

1. Reproduced the downstream failure on a clean ModelingToolkitBase master checkout.
2. Traced both errors through LineSearch JVP preparation rather than changing the downstream tolerance or boundary residual.
3. Checked the behavior at the exact introducing boundary: d21952c3 passes and c6b26f4a fails.
4. Added the smallest owning-package change: retain explicit backend precedence and otherwise reuse an analytic Jacobian.
5. Verified the focused regression on stable/pre, the full NonlinearSolveFirstOrder Core group, and the complete downstream boundary-value file on stable/pre.

Next step is CI monitoring; any follow-up commit and its local verification will be recorded in a new comment.

[ChrisRackauckas-Claude]

Verification correction: I am withdrawing the previously reported full Core result after a process-path audit found a separate active Pkg.test process resolving outside the current workspace. That process has been stopped. A replacement Core run is in progress with a fresh workspace-local depot and environment; it prints and asserts pathof(NonlinearSolveFirstOrder) against this PR checkout before invoking the standard package test runner. I will post the exact path and result when it completes.

[ChrisRackauckas-Claude]

Follow-up commit a5a889a4c fixes the rectangular analytic-JVP failure exposed by the first CI run.

The out-of-place analytic fallback reshaped jac(u, p) * v to size(u). That is incorrect for an overdetermined problem whose residual shape differs from the state shape; it now reshapes to size(fu). A 3×1 residual / 3×2 Jacobian regression covers the output shape and value.

Final local verification:

• SciMLJacobianOperators Core on Julia 1.12.6: 480 scalar, 288 in-place, 482 out-of-place, and 11 copy tests passed; Pkg.test exited 0.
• Exact root test/Core/core_tests__item1.jl: exited 0 with local path assertions.
• NonlinearSolveFirstOrder Core on Julia 1.12.6: full group exited 0 with local pathof and workspace TMPDIR assertions; notable groups include NLLS 472/472, NewtonRaphson 309/309, TrustRegion 1176/1176, and operator Jacobians 14/14.
• Final ModelingToolkitBase/test/bvproblem.jl on Julia 1.12.6: exited 0; Lotka 4/4, cost compilation 6 pass + 2 expected broken, parameter estimation 2/2.
• Julia 1.13.0-rc1 rectangular analytic-JVP smoke: exited 0 with the expected 3×1 result.
• Runic 1.7.0 over the whole repository and git diff --check: exited 0.

[ChrisRackauckas-Claude]

I split the rectangular analytic-JVP fix out of this PR into prerequisite #1063. The branch now contains only the line-search change plus the NonlinearSolveFirstOrder 2.2.1 / SciMLJacobianOperators 0.1.16 release floor.

Post-split verification used the exact #1063 checkout and asserted both version and package path before running the regression: 5/5 tests passed on Julia 1.12.6. The previous combined full-Core and downstream results remain applicable because they ran the same two source changes together.

I fetched the remote branch and confirmed that it contained only the two agent-authored commits before replacing the duplicated SciMLJacobianOperators commit with the dependency-floor commit. The update used --force-with-lease; no contributor commit was overwritten.

[ChrisRackauckas-Claude]

The post-split CI run reached the expected release-order gate: every package-resolving job stopped immediately because SciMLJacobianOperators 0.1.16 is not registered yet. Runic and spell checks passed. There is no runtime test failure to debug in this revision; I will rerun the package jobs after #1063 merges and 0.1.16 is tagged.

[ChrisRackauckas-Claude]

Version-order note: draft #1060 independently uses the same provisional next versions (SciMLJacobianOperators 0.1.16 and NonlinearSolveFirstOrder 2.2.1). If #1060 merges first, this branch must rebase to NonlinearSolveFirstOrder 2.2.2 and require the post-#1063 SciMLJacobianOperators 0.1.17 release. The current floors are valid only if #1063/#1062 land first; they are not evidence that #1060 contains the residual-shape fix.

[ChrisRackauckas-Claude]

Rebased onto current master at 41681c5e4 after #1060 merged and updated the prerequisite floor to SciMLJacobianOperators 0.1.17; the branch head is now ad591a395.

Final post-rebase validation used the exact #1063 source overlay and asserted local paths plus versions (NonlinearSolveFirstOrder 2.2.2, SciMLJacobianOperators 0.1.17). The full Julia 1.12 Core harness exited 0: General NLLS 472/472, NewtonRaphson 309/309, PseudoTransient 51/51, TrustRegion 1176/1176, LevenbergMarquardt 33/33, Structured Jacobians 16/16, and SciMLOperator Jacobians 14/14. Whole-repository Runic and git diff --check also passed.

The temporary #1063 source overlay was removed before committing; this PR contains only its intended FirstOrder changes and declares the released-version dependency.

[ChrisRackauckas-Claude]

The five immediate failures are all prerequisite-resolution failures before tests run:

• Documentation, Adjoint LTS, Bounds, and QA resolve the registry's current SciMLJacobianOperators 0.1.16, which cannot satisfy this PR's required 0.1.17 floor.
• NonlinearSolveFirstOrder downgrade develops the sibling source at its master version 0.1.16 and exits with the same empty intersection against the 0.1.17 floor.

This is the declared #1063 merge/release dependency. All five jobs stop during environment resolution before loading or running branch code. The exact #1063 0.1.17 source was developed locally for the final full Core run, which passed with path and version assertions; no CI failure is being bypassed or hidden here.

[ChrisRackauckas-Claude]

Scheduled CI follow-up at 2026-07-13 18:59 EDT: I fetched and inspected all 37 completed failing job logs. All 37 stop during dependency resolution on the unreleased SciMLJacobianOperators 0.1.17 floor (Unsatisfiable requirements / empty 0.1.16-vs-0.1.17 intersection); none reached this branch's runtime tests. Fourteen checks pass, three matrix placeholders are skipped, and twenty jobs remain queued. The actionable gate remains merging/releasing prerequisite #1063; no #1062 code fix is indicated by this run.

[ChrisRackauckas-Claude]

Local downstream verification on BoundaryValueDiffEq master d31a88a430c is complete (Julia 1.10.11, official package Core groups, no tests skipped or reclassified).

With the released dependency set, MIRK reported:

MIRK Basic Tests | 189 pass | 2 fail | 4 error | 195 total

With this PR at ad591a395 plus prerequisite #1063 at af904522 developed into the environment (NonlinearSolveFirstOrder 2.2.2, SciMLJacobianOperators 0.1.17, LinearSolve 4.3.0), MIRK reported:

MIRK Basic Tests | 194 pass | 2 fail | 1 error | 197 total

All three line-search Dual regressions are gone:

• Swirling Flow III: 1/1 passed.
• Test interpolant evaluation with big defect: 3/3 passed.
• Test maxsol and minsol: 2/2 passed.

The two remaining failures are unchanged missing_second_order_ad warnings. The one remaining error is the unchanged OptimizationDISparse constraint-Jacobian arity failure in the parameter-tuning optimize path. I independently bisected those to OptimizationIpopt 1.3.0 / Optimization.jl #1257; they are unrelated to this PR.

The Shooting Core group also exited 0 with the same combined sources:

Basic Shooting Tests             | 82 pass | 82 total
NLLS BVProblem Tests             | 56 pass | 56 total
Periodic Orbit Shooting Tests    | 40 pass | 40 total

This confirms the analytic-Jacobian line-search path fixes both BoundaryValueDiffEq downstreams while preserving the non-line-search baseline results.