Reuse homotopy caches across implicit ODE stages

[ChrisRackauckas-Claude]

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

Dependencies and scope

SciML/NonlinearSolve.jl#1100 is merged and released in NonlinearSolve 4.24 / NonlinearSolveBase 2.38. It provides reusable HomotopySweep and KantorovichHomotopy init/reinit!/solve! caches through public APIs.

The related ImplicitDiscreteSolve cache work, #4042, and its registered NonlinearSolveBase 2.38 compatibility correction, #4055, are merged.

This PR only changes implicit ODE stage solves through OrdinaryDiffEqNonlinearSolve.HomotopyNonlinearSolveAlg.

What this changes

HomotopyNonlinearSolveAlg previously rebuilt a HomotopyProblem and used one-shot solve for every implicit stage. This PR:

• initializes one continuation cache for HomotopySweep and KantorovichHomotopy;
• refreshes the stage initial guess, parameter tuple (including stage time), and tolerances with public reinit!, then calls public solve!;
• rebuilds the continuation cache after integrator resizing;
• retains one-shot solve for homotopy algorithms without this reusable cache path, including arc-length continuation;
• imports the continuation algorithms from their declaring public module, NonlinearSolveBase;
• documents the cached and fallback behavior.

IDSolve should not be replaced by this wrapper. Its controller is a timestep/path-continuation policy that uses the previous accepted state and every Newton contraction ratio, whereas HomotopyNonlinearSolveAlg solves an individual implicit ODE stage equation. The merged IDS path now has its own reusable init/reinit!/solve! implementation.

Solver comparison

On the exact 51-solve IDS recurrence benchmark, using a fresh environment pinned to registered SciMLBase 3.39.1, NonlinearSolve 4.24, and NonlinearSolveBase 2.38, two independent 20-sample runs measured:

• cached IDS: 0.09914-0.10570 ms, 192 bytes, 561 residual calls;
• cached Kantorovich: 0.27951-0.29414 ms, 32 bytes, 2,040 residual calls;
• cached Sweep: 0.30429-0.30484 ms, 32 bytes, 2,346 residual calls.

IDS was 2.64x-2.97x faster than Kantorovich and 2.88x-3.08x faster than Sweep on its native problem. To reduce sequential-load bias, a separate 31-trial paired benchmark with 200 recurrences per timing measured 250.94 microseconds for Kantorovich versus 268.39 microseconds for Sweep, a 6.55% median paired time reduction. Kantorovich also used 13.04% fewer residual calls.

A rotated 31-trial variant benchmark also measured 250.27 microseconds / 2,040 residual calls for strict Kantorovich's constant predictor, 243.56 microseconds / 1,989 calls with its trust-monitored secant predictor, and 242.25 microseconds / 1,989 calls for secant plus non-strict acceptance. The useful Sweep predictor machinery is therefore already available in Kantorovich; this single smooth recurrence is not evidence to change its defaults.

That does not make strict Kantorovich the safer default. In a controlled slow-corrector check with a 0.99 residual-contraction ratio, strict Kantorovich rejected the converged corrector from u0 = 0, while non-strict Kantorovich and Sweep both reached u = 0.1. At a fold, Kantorovich failed and the Sweep-to-arc-length polyalgorithm succeeded. Strict Kantorovich is therefore more conservative/less robust on these checks, not numerically “unstable”; Sweep remains the default, with arc length needed to round folds.

Local verification

Final rebased head 967902f1b3 on current upstream/master (fd3f95bb99):

• The official homotopy_nlsolve_tests.jl file passed 38/38. This includes cache identity across stages, time-dependent parameter refresh, resize rebuild, exact warmed @allocated == 0 for both cached algorithms, out-of-place/static arrays, arc-length fallback, fold rejection, and DAE rejection.
• GROUP=QA julia --project=lib/OrdinaryDiffEqNonlinearSolve -e 'using Pkg; Pkg.test()' exited 0: JET 13/13 with 33 pre-existing broken cases; Aqua 19/19; package tests passed.
• Runic 1.7 --check on every changed Julia file and git diff --check upstream/master exited 0. Whole-tree Runic on current master is tracked independently by Run Runic on the disco-optimizations files #4064; the homotopy diff is not involved.
• A 20-sample, 1,000-stage ImplicitEuler comparison measured 7.4617735 microseconds/stage and 0 bytes/stage with the reused cache versus 9.316223 microseconds/stage and 5,520 bytes/stage when forcing a rebuild before every stage: 1.2485x faster / 19.91% lower median time, with identical endpoints.
• The current full Core run did not complete: after 1h50 and successful groups through CheckInit, the nested-AD test raised UndefVarError: issquare not defined in LinearSolveForwardDiffExt with registered LinearSolve 5.3, before reaching the homotopy file. This dependency regression is being handled separately; it was not skipped or silenced here.

[ChrisRackauckas-Claude]

CI dependency checkpoint: the docs job failed during package resolution before loading OrdinaryDiffEq code. General currently provides NonlinearSolve through 4.22.0, while this stacked PR requires the new reusable cache API in NonlinearSolve 4.23 from SciML/NonlinearSolve.jl#1100.

This is the expected prerequisite boundary, not a source/test failure. No transient [sources] entry is being committed. Once #1100 merges and NonlinearSolve 4.23 is registered, rerun this PR's CI.

[ChrisRackauckas-Claude]

Clean-rebase update: this branch is now stacked only on #3984 and depends on SciML/NonlinearSolve.jl#1100. The unrelated #3983 commit has been removed.

I reran the official ODEDIFFEQ_TEST_GROUP=Core package path on Julia 1.12 against the local NonlinearSolve 4.23/Base 2.37 implementation and LinearSolve 5.0.1. Pkg.test() exited 0 with OrdinaryDiffEqNonlinearSolve tests passed; the homotopy block passed 38/38, including exact warmed-stage zero-allocation assertions for both HomotopySweep and KantorovichHomotopy. Whole-repository Runic --check and git diff --check exited 0.

CI has restarted on the clean stack. Dependency resolution for this PR remains expected to wait for NonlinearSolve 4.23/#1100 to be released; no transient source override is committed.

[ChrisRackauckas-Claude]

Current-master rebase and exact-head validation update (b1f231d622):

• Remote branch history was audited (only the two prior bot-authored commits), then replaced with the single focused commit using an explicit force-with-lease. The PR is now current with upstream/master and mergeable.
• Official ODEDIFFEQ_TEST_GROUP=Core on Julia 1.12.6 exited 0 with Testing OrdinaryDiffEqNonlinearSolve tests passed; the homotopy blocks passed 38/38 and 16/16.
• Official QA exited 0: JET 13 passed / 33 existing broken; Aqua 19/19.
• Whole-repository Runic 1.7 and git diff --check exited 0.
• Exact local-head 1,000-stage benchmark: reuse 7.460 μs and 0 bytes/stage versus forced rebuild 9.244 μs and 5,520 bytes/stage (1.239×), identical endpoints.

During the requested IDS comparison I also found a separate clean-master ImplicitDiscreteProblem reinitialization bug: initial init performs the discrete initialization map, while generic reinit! currently omits it, making fresh solve and init → reinit! → solve! differ by one update. That is being fixed and tested as a separate focused OrdinaryDiffEq PR; it does not require changing this homotopy-cache implementation or replacing the merged IDS controller.

[ChrisRackauckas-Claude]

Exact-head CI dependency checkpoint: the apparent broad red matrix has one common pre-test cause.

Representative logs from Documentation, root CI, downgrade sublibraries (including OrdinaryDiffEqNonlinearSolve), and downstream integration all stop during Pkg resolution with:

NonlinearSolve possible versions: 0.1.0 - 4.23.0
restricted to versions 4.24.0 - 4 by OrdinaryDiffEqNonlinearSolve — no versions left

No failing job reached package loading or this PR's tests. This is the expected SciML/NonlinearSolve.jl#1100 release boundary; local exact-head Core/QA used that PR's 4.24 implementation and passed. I am leaving the public dependency explicit rather than committing a transient source override. CI should be rerun after #1100 merges and NonlinearSolve 4.24 is registered.

[ChrisRackauckas-Claude]

Rebased onto current upstream/master 938dae56bf and pushed final head 01929d4f5b after auditing the remote branch. No contributor commits were overwritten; the remote contained only the previous PR commit. git range-diff eeb76608b8^! 01929d4f5b^! reports =.

The full official OrdinaryDiffEqNonlinearSolve Core group has now completed locally on that patch-identical implementation with exit code 0 and the final line Testing OrdinaryDiffEqNonlinearSolve tests passed. In particular, the homotopy blocks passed 38/38 and 16/16, including exact warmed zero-allocation cache tests for Sweep and Kantorovich. The PR body now has the complete test, allocation, performance, IDS comparison, and dependency evidence. Fresh CI is starting on the rebased head.

[ChrisRackauckas-Claude]

Early CI note: the immediate documentation, downstream, downgrade, and StochasticDiffEq failures all resolve the current master copy of ImplicitDiscreteSolve and stop at its NonlinearSolveBase = "2.40" bound, while General contains only 2.38. Representative logs show the direct 2.38-vs-2.40 unsatisfiable requirement; the downgrade resolver reports the same missing UUID after filtering the nonexistent version. This is the already-isolated master metadata regression fixed by #4051, not a homotopy-cache failure. I will reclassify/rerun after #4051 lands; the jobs that can resolve without the umbrella package are continuing.

[ChrisRackauckas-Claude]

Rebased onto current fd3f95bb99 and lease-protected pushed as 967902f1b3. Before the push I fetched/audited the remote branch; it still contained only the prior PR commit. The range-diff is implementation-equivalent, with the sole intentional adjustment that the new test imports ArcLengthContinuation and KantorovichHomotopy from their declaring public module (NonlinearSolveBase). Rebased local results: homotopy tests 38/38, JET 13/13 (+33 pre-existing broken), Aqua 19/19, changed-file Runic and diff check clean, and cache benchmark 7.4618 vs 9.3162 μs/stage with 0 vs 5,520 bytes/stage.

[ChrisRackauckas-Claude]

Fresh benchmark audit after checking the environment manifest:

The prior workspace environment had subsequently been developed onto local NonlinearSolve sources, so its current manifest could not substantiate the registered-version label. I discarded that as version evidence and rebuilt a fresh local environment. The new manifest pins registered SciMLBase 3.39.1, NonlinearSolve 4.24.0, and NonlinearSolveBase 2.38.0. Only OrdinaryDiffEq's DiffEqBase, Core, and ImplicitDiscreteSolve sublibraries are developed locally, from the explicit test-only stack 5d0dfe43f8 (this PR plus the SciMLBase-floor prerequisite and #4064); the two prerequisites do not change the benchmarked solver implementations.

Two independent 20-sample runs of the exact 51-solve recurrence measured:

• cached IDS: 0.09914–0.10570 ms, 192 bytes, 561 residual calls;
• cached Kantorovich: 0.27951–0.29414 ms, 32 bytes, 2,040 residual calls;
• cached Sweep: 0.30429–0.30484 ms, 32 bytes, 2,346 residual calls.

A separate interleaved check (31 paired trials, 200 complete recurrences per timing) measured 250.94 μs for Kantorovich versus 268.39 μs for Sweep, a 6.55% median paired time reduction.

A rotated predictor/acceptance check (31 trials × 100 recurrences) measured:

• constant + strict Kantorovich: 250.27 μs, 2,040 residual calls;
• constant + non-strict: 250.14 μs, 2,040 calls;
• secant + strict: 243.56 μs, 1,989 calls;
• secant + non-strict: 242.25 μs, 1,989 calls;
• Sweep + secant: 279.75 μs, 2,346 calls.

Thus Kantorovich already has the useful trust-monitored secant machinery from Sweep. This smooth recurrence is not enough evidence to change defaults: on the S-fold, constant/strict and secant/strict Kantorovich returned ConvergenceFailure, secant/non-strict Kantorovich and Sweep returned MaxIters, while Kantorovich → arc length returned Success at the connected root u = 2.1038034027355366.

The PR body now contains the corrected pinned-stack evidence. No code or branch change was made for this audit.

[ChrisRackauckas-Claude]

Current-head CI triage (967902f1b3):

• Documentation, root test, and OrdinaryDiffEqNonlinearSolve sublibrary jobs all stop
  before package loading because General currently ends at SciMLBase 3.39.1 while the
  current master sublibraries require 3.40. The clean-master failure was reproduced and
  bisected to d78527d4 (Import the SciMLOperators names from their owner, and require SciMLBase 3.40 #4059); focused fix: Restore SciMLBase 3.39 compatibility for owner imports #4066.
• Runic reaches the source tree and reports only the independently merged DISCO files,
  principally test/Integrators_I/disco_tests.jl. Clean master was reproduced and
  bisected to 408a7a30; formatting-only fix: Run Runic on the disco-optimizations files #4064.

On #4066, the full official OrdinaryDiffEqNonlinearSolve Core group passed with
registered SciMLBase 3.39.1. Its SciMLBase 3.40.1 control passed every affected suite
and then stopped at the separately assigned LinearSolveForwardDiffExt.issquare
regression in Nested AD.

The homotopy branch is unchanged. A test-only stack of this head plus the two focused
fixes cherry-picked without conflicts. With local OrdinaryDiffEq sublibraries,
registered SciMLBase 3.39.1, NonlinearSolve 4.24.0, and NonlinearSolveBase 2.38.0:

• homotopy_nlsolve_tests.jl: 38/38 passed, and
  homotopy_default_nlsolve_tests.jl: 16/16 passed;
• whole-tree Runic 1.7 and git diff --check: exit 0;
• 1,000-stage cache benchmark: 7.53971 microseconds and 0 bytes/stage reused versus
  9.334121 microseconds and 5,520 bytes/stage forced-rebuild (1.23799x), identical
  endpoint.

No CI-only workaround or unrelated change has been added to #3985. It should be rerun
after the two clean-master prerequisites land.

[ChrisRackauckas-Claude]

Broader registered-version comparison after the IDS-specific benchmark:

Environment: registered SciMLBase 3.39.1, NonlinearSolve 4.24.0, and NonlinearSolveBase 2.38.0. I compared HomotopySweep, strict/constant Kantorovich, strict/secant Kantorovich, non-strict/secant Kantorovich, and arc length on eight no-fold paths derived from the homotopy regression cases plus one coupled curved path. Every algorithm used the same NewtonRaphson() inner corrector, the same 0.1 initial-step default, min_dλ = 1e-4, maxsteps = 2000, and solve maxiters = 100.

Timing used 31 trials of 50 solves per algorithm with algorithm order rotated each trial. All returned accurate endpoints in every trial.

problem	Sweep μs / residual calls	default strict Kantorovich μs / calls	best strict Kantorovich configuration
smooth cubic	26.046 / 54	23.873 / 40	constant: 23.873 / 40
quadratic blend	31.838 / 87	25.991 / 51	constant: 25.991 / 51
atan turn	28.949 / 65	27.276 / 53	constant: 27.276 / 53
linear cubic path	22.143 / 38	26.388 / 55	secant: 21.848 / 33
sharp tanh turn	210.909 / 728	224.541 / 744	secant: 169.837 / 539
anchored branch	31.934 / 87	25.959 / 51	constant: 25.959 / 51
decreasing span	21.608 / 37	19.559 / 22	constant: 19.559 / 22
coupled curved path	30.147 / 57	27.394 / 43	constant: 27.394 / 43

Default Kantorovich was faster on 6/8 paths; the median of the eight paired default-Kantorovich/Sweep ratios was about 0.91. One of the standard Kantorovich predictor settings was fastest on every path, but neither setting dominated: constant won six, while secant mattered on the linear and sharp-turn paths. Testing secant with predictor_order = 1 did not dominate either; it helped smooth/atan/linear paths but regressed the sharp turn to 226.568 μs and 749 calls. This supports keeping the existing explicit predictor/order controls rather than changing defaults from this small suite.

Strict and non-strict secant had identical residual counts on all eight no-fold paths, so the strict safety check had no observed ordinary-path cost. On the cubic S-fold, strict Kantorovich correctly returned ConvergenceFailure; with the larger corrector budget, Sweep and non-strict Kantorovich reached the target in 86 and 67 calls, respectively, but monotone natural-parameter continuation cannot continuously round that fold, so those successes necessarily switch branches. Arc length reached the connected endpoint in 197 calls. The strict Kantorovich to arc-length polyalgorithm also succeeded, but used 631 calls because the strict first stage shrank repeatedly before handing off.

Conclusion: Kantorovich is often faster than Sweep, but is not categorically better. Strict rejection is path-safety behavior, not random instability; arc length remains the correct fold tool. The results do not justify replacing the default Sweep to arc-length polyalgorithm, and they strengthen the conclusion that generic homotopy should not replace the specialized IDS recurrence, which remained 2.64–2.97x faster than cached Kantorovich in the exact IDS benchmark. No code or branch change was made for this comparison.

[ChrisRackauckas-Claude]

Rebased onto current upstream/master a7080c6f6d and lease-protected pushed as 7c027e4ce8.

Before pushing I fetched and audited the fork branch and PR: it still contained exactly the previous single PR commit (967902f1b3), with no reviews or contributor commits. The push used an explicit SHA lease. The final branch is 1 commit ahead / 0 behind master and GitHub reports it mergeable.

Validation:

• git range-diff 8b86e7dc04^! 7c027e4ce8^! reports =; the final master-only update (Declare fix_dt_at_bounds!, handle_tstop!, initialize_tstops, initialize_saveat, initialize_d_discontinuities public #4081) does not overlap this PR.
• The full official GROUP=OrdinaryDiffEqNonlinearSolve_Core rerun on the patch-identical rebased implementation exited 0 after 3326 seconds. The package and root harnesses both reported tests passed; homotopy sets passed 38/38 and 16/16, and the remaining matrix-free/polyalgorithm sets completed 14/14 and 8/8.
• On the exact final SHA, the two existing homotopy regression files ran in a fresh local environment using final-head local sublibraries (including OrdinaryDiffEqCore 4.13.0), NonlinearSolve 4.24.0, NonlinearSolveBase 2.39.1, and SciMLBase 3.41.0: 38/38 and 16/16 passed. The stage-cache block passed 10/10, including warmed zero-allocation checks.
• Official GROUP=OrdinaryDiffEqNonlinearSolve_QA exited 0: JET 13 passed / 33 pre-existing broken; Aqua 25/25.
• Whole-tree Runic --check --diff, git diff --check, and clean-worktree checks exited 0 on the final SHA.

Fresh CI is running on 7c027e4ce8.

[ChrisRackauckas-Claude]

Post-rebase hosted CI audit for 7c027e4ce84cbb1aa3fa8f5dca8a832de82a77ad:

• Green: OrdinaryDiffEqNonlinearSolve Core on Julia current/LTS/pre, QA, downgrade compat, both benchmark lanes, ImplicitDiscreteSolve downgrade, and SimpleImplicitDiscreteSolve downgrade.
• The root GPU red is an exact clean-master match: Hairer4 + CSC Jacobian + diagonal CuArray mass matrix reports the same 3.62e-3 / 7.02e-3 CPU-GPU mismatch on current master and an earlier master run. This is outside the changed files and is tracked by Hegedüs warm-start default (#3991) amplifies round-off ~1e13x on the Newton-Krylov path (GPU DAE test failure) #4034.
• The DiffEqDevTools downgrade job lost its self-hosted runner connection. Both the CLI token and GitHub App were denied permission to rerun that single job (403), so I did not create a no-op commit to retrigger CI.
• OrdinaryDiffEqBDF [GPU] eventually received a runner and passed; hosted CI is now fully terminal.
• Other red checks were reproduced on clean master and are being handled as separate follow-ups; no unrelated fix was added to this branch.

Final branch audit: current upstream/master is a7080c6f6d992edee79fa40450e0bd6eb0b8f735; the branch is 0 behind / 1 ahead, local and remote heads match, the worktree is clean, and GitHub reports the draft PR mergeable.