NonlinearSolveCorrector: NonlinearSolve.jl algorithms as the PALC/Natural corrector by ParadaCarleton · Pull Request #345 · bifurcationkit/BifurcationKit.jl

ParadaCarleton · 2026-10-02T17:16:01Z

Adds NonlinearSolveCorrector (package extension NonlinearSolveExt, loaded by using NonlinearSolve), which lets any NonlinearSolve.jl algorithm be the corrector of PALC and Natural and the solver of solve(prob, corrector, NewtonPar()), with Jacobian reuse within and across correctors, matrix-free Jacobians through a Krylov linsolve, and [p_min, p_max] kept by algorithms with native bounds; the default corrector (nothing, BifurcationKit's own Newton) is unchanged. test/continuation/test_nonlinearsolve_corrector.jl checks that the branch matches BifurcationKit's Newton point by point past the fold, and that the reuse, trust-region, bounded and matrix-free configurations stay on it. It stays a draft until SciML/NonlinearSolve.jl#1353 (reinit!(...; reuse_jacobian = true)) is released; the projecting test runs only once NonlinearSolve defines BoundsProjection (SciML/NonlinearSolve.jl#1354).

…gh NonlinearSolve.jl A package extension (NonlinearSolveExt) solves the bordered system of PALC, the fixed-parameter Newton system of Natural and solve(prob, corrector, NewtonPar) with any NonlinearSolve algorithm, one cache per corrector so its Jacobian reuse and factorizations survive between solves. corrector = nothing keeps BifurcationKit's Newton method.

…orrectors via reinit!(...; reuse_jacobian = true)

… steps, bounds for algorithms that support them

…arSolve has project_bounds

…(), the keyword NonlinearSolve settled on

rveltz · 2026-10-04T14:03:33Z

That part looks pretty cool, thank you!

rveltz · 2026-10-04T14:05:39Z

+        x, p = palc_state(w), w[end]
+        Fx = BK.residual(ctx.prob, x, BK.set(ctx.par, ctx.lens, p))
+        Nx = BK.arc_length_eq(ctx.dotθ, x, ctx.z0u, p - ctx.z0p, ctx.τu, ctx.τp, ctx.θ, ctx.ds)
+        return vcat(Fx, Nx)


this part is bad for GPU, that is why I developed the BorderedArrays.

Agreed: vcat(Fx, Nx), w[begin:end-1] and w[end] force a flat array with scalar indexing. I'd use BorderedArray directly, but NonlinearSolve needs u0::AbstractArray, and BorderedArray isn't one.

The closest thing I see is RecursiveArrayTools.ArrayPartition(x, [p]), which is already a BifurcationKit dependency. It's an AbstractVector that NonlinearSolve accepts, and it broadcasts part by part, so x stays a CuArray and there's no scalar indexing or concatenation. p would live in a one-element array of the same type as x. Would that work for you, or would you rather make BorderedArray an AbstractVector?

I think it would work.

A few remarks though:

I want to separate u and parameters p: this pattern appears everywhere: periodic orbit (there is a parameter-like constraint), hopf/ fold continuations, continuation of bifurcations of PO (2 constraints), etc.

if you cannot use BorderedArray, at least allow for the user to pass AbstractArray and not just AbstractVector, this is handy for solving PDE

we have the issue BifurcationKit not working with ArrayPartitions ? #112

if you pass an ArrayPartition to a nonlinear solve, does it return a ArrayPartition?

rveltz · 2026-10-04T14:07:01Z

+    bordered_matrix(ctx::PALCContext, J::AbstractMatrix, dFdp) = vcat(hcat(J, dFdp), hcat(LA.adjoint(ctx.row), ctx.c))
+
+    # products of `[J ∂ₚF; row c]` with `v = (vx, vp)`
+    function palc_jvp(v, w, ctx::PALCContext)


Maybe use BLS from Linsolve.jl ?

Do you mean running BifurcationKit's own BLS as the linear solver inside NonlinearSolve? That works, I think: linsolve takes a LinearSolve.jl algorithm, and LinearSolve.LinearSolveFunction wraps any function as one. The function would read J, ∂ₚF and the arclength row from the problem's parameter (the PALCContext) and call BK.solve_bls with the continuation's bls. The bordered system would then be solved with BorderingBLS/MatrixFreeBLS, as in BifurcationKit's Newton, instead of building the dense [J ∂ₚF; row c] (bordered_matrix) or forming its products here. If that's what you had in mind, I'll make that the default and keep the dense matrix only for a user-given dense linsolve.

I meant https://github.com/SciML/LinearSolve.jl and https://github.com/SciML/BorderedLinearSolve.jl

…rrector

rveltz · 2026-10-06T07:11:38Z

+
+    #━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
+    # The problem's Jacobian (and `∂ₚF` for PALC) at the last point asked for, built once for the products and the dense matrix of one iterate
+    struct Linearization{Tx, TJ, Tdp}


Maybe call it BorderedJacobian

rveltz · 2026-10-06T07:25:11Z

+    #━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
+    # Solves of the bordered system of PALC in `(x, p)`. The parameter of the NonlinearProblem is a `PALCContext`.
+    # `row` and `c` are the last row `[θ A τu' (1-θ) τp]` of the bordered matrix of `MatrixBLS`.
+    struct PALCContext{Tprob, Tpar, Tlens, Tdot, Tu, Tp, Tθ, Trow}


PALCCache ?

These are already in ContState
z0u::Tu
z0p::Tp
τu::Tu
τp::Tp

rveltz · 2026-10-06T07:25:25Z

+    struct PALCContext{Tprob, Tpar, Tlens, Tdot, Tu, Tp, Tθ, Trow}
+        prob::Tprob
+        par::Tpar
+        lens::Tlens


this is in prob, isn't it?

rveltz · 2026-10-06T07:25:47Z

+        θ::Tθ
+        ds::Tp
+        row::Trow
+        c::Tp


what is this?

rveltz · 2026-10-06T07:27:22Z

Just out of curiosity, can you add tests for Fold/Hopf continuation please to see if the interface you propose is good enough?

ParadaCarleton · 2026-10-11T01:53:47Z

Thanks for the review. Plan, following your remarks:

State: u and p stay separate as ArrayPartition(x, [p]), and the solver accepts any AbstractArray x (not only vectors). To answer your question: NonlinearSolve returns an ArrayPartition when given one (checked with Float64 and Float32, dense Jacobian), so the scalar indexing and vcat go away.
Bordered solve: through BorderedLinearSolve.jl on top of LinearSolve.jl, replacing my hand-built bordered matrix and products. I checked it against a dense solve with a matrix J and a matrix-free J with GMRES, in Float64 and Float32 (errors at rounding level, element type kept).
Cleanup: Linearization -> BorderedJacobian, PALCContext -> PALCCache without the fields already in ContState/prob, a docstring for c, and Fold/Hopf continuation tests.

One catch: BorderedLinearSolve.jl is not registered yet (v0.1.0, requires LinearSolve 5 and SciMLBase 3), so BifurcationKit can only depend on it after a release. Should this PR wait for that, or should it land with BorderedLinearSolve as a weak dependency in the extension in the meantime?

…tems by BorderedLinearSolve.jl The unknown of PALC is ArrayPartition(x, [p]) for any AbstractArray state x, without vcat or scalar indexing. Its Jacobian is BorderedJacobian, a SciMLOperator given to NonlinearSolve as jac_prototype, whose linear systems are solved by BorderedLinearSolve.jl (keyword `bls`, default BorderingBLS with the `linsolve` of `alg` as the inner solver and one correction sweep) through a LinearSolveFunction in the Newton step of `alg`. This replaces the hand-built bordered matrix and products, and `corrector_jacobian`/`corrector_matrix`: a preconditioner now sees the operator of J alone, with the state in `A.u` and the parameters in `A.p`. PALCContext becomes PALCCache, which reads z0, τ0, ds from the continuation state and the lens from the problem, and documents `row` and `c`. BorderedLinearSolve.jl is a weak dependency (extension trigger) loaded from its git URL through [sources] until it is registered.

…e default corrector

…p from the problem's jet

…lem's jet (R01), not finite differences

…e as the bordered jvp With a matrix-free `J` (a function `dx -> J dx`), PALC's corrector gives NonlinearSolve the product with the whole bordered Jacobian `[J ∂ₚF; row' c]` (`palc_jvp`) and solves it with the iterative linear solver of `alg`, whose `precs` reads the iterate through `corrector_jacobian` and `corrector_matrix`. BorderedLinearSolve.jl is used only when `J` is a matrix. This narrows the use of BorderedLinearSolve (decision 2 of the redesign) to the matrix case. NonlinearSolve rebinds only its own JacobianOperator to each iterate, and BorderedLinearSolve has no API to swap the inner operator of its cache while keeping that cache, so a matrix-free `J` inside BLS would rebuild the preconditioner at every iterate and `jacobian_reuse` would have no effect (the precs test failed). The unknown stays `ArrayPartition(x, [p])`. `J` and `∂ₚF` are computed once per iterate (`iterate_jacobian!`), however many products the Krylov solver takes.

…p_min or p_max A corrector whose solution lies past the bound was rejected as not converged, so PALC halved the step toward the bound without ever reaching it. The corrector now reports convergence from its residual alone; a predictor past the bound is still corrected at the bound by Natural.

…es past p_min or p_max" This reverts commit 10f3a59.

…n or p_max is corrected at that bound from where it ended, as a predictor past the bound is

…r p_max is converged The corrector no longer rejects a solution past the bound; PALC corrects it again at the bound (previous commit). Rejecting it made PALC halve the step and reach the bound only through a later predictor past it. Tested on the circle x² + a² = 1 below its centre, where the corrector raises a past a predictor that is short of p_max: the branch ends on p_max exactly.

…ringBLS clashes with BifurcationKit's

rveltz · 2026-10-20T16:57:59Z

did you see: https://docs.sciml.ai/NonlinearSolve/stable/native/solvers/#NonlinearSolveBase.ArcLengthContinuation ?

ParadaCarleton added 5 commits October 1, 2026 17:55

NonlinearSolveCorrector: reuse_jacobian carries the Jacobian across c…

6026bef

…orrectors via reinit!(...; reuse_jacobian = true)

NonlinearSolveCorrector: matrix-free Jacobians, Jacobian reuse across…

be38741

… steps, bounds for algorithms that support them

Test NonlinearSolveCorrector with a projecting algorithm when Nonline…

756dbf0

…arSolve has project_bounds

Test the projecting corrector with bounds_handling = BoundsProjection…

e4910ce

…(), the keyword NonlinearSolve settled on

rveltz reviewed Oct 4, 2026

View reviewed changes

ParadaCarleton force-pushed the nonlinearsolve-corrector branch from aaec980 to e4910ce Compare October 5, 2026 00:07

ParadaCarleton added 2 commits October 5, 2026 13:01

Merge remote-tracking branch 'upstream/master' into nonlinearsolve-co…

8f136ae

…rrector

test: PALC's type now carries the corrector's type parameter

fa78753

rveltz reviewed Oct 6, 2026

View reviewed changes

Comment thread ext/NonlinearSolveExt/NonlinearSolveExt.jl Outdated

θ::Tθ

ds::Tp

row::Trow

c::Tp

rveltz Oct 6, 2026

Copy link
Copy Markdown

Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

what is this?

ParadaCarleton added 11 commits October 10, 2026 22:51

Test NonlinearSolveCorrector on PALC with a matrix state

0e71865

Test NonlinearSolveCorrector in Fold and Hopf continuation against th…

dc31044

…e default corrector

Test that the bordered Jacobian of NonlinearSolveCorrector takes ∂F/∂…

ff0bb44

…p from the problem's jet

NonlinearSolveCorrector: ∂F/∂p of the bordered Jacobian from the prob…

ce6893c

…lem's jet (R01), not finite differences

Revert "NonlinearSolveCorrector: accept a PALC corrector that converg…

50b6b76

…es past p_min or p_max" This reverts commit 10f3a59.

PALC: a corrector from an interior predictor that converges past p_mi…

c3063db

…n or p_max is corrected at that bound from where it ended, as a predictor past the bound is

NonlinearSolveCorrector docs: import BorderedLinearSolve, whose Borde…

2b135a2

…ringBLS clashes with BifurcationKit's

ParadaCarleton mentioned this pull request Oct 11, 2026

Register v0.1.0 in General SciML/BorderedLinearSolve.jl#5

Open

NonlinearSolveCorrector docs and comments: state what the corrector does

f26ae26

ParadaCarleton force-pushed the nonlinearsolve-corrector branch from e7a7cbd to f26ae26 Compare October 24, 2026 01:43

Conversation

ParadaCarleton commented Oct 2, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

rveltz commented Oct 4, 2026

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

rveltz commented Oct 6, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

ParadaCarleton commented Oct 11, 2026

Uh oh!

rveltz commented Oct 20, 2026

Uh oh!

This branch has not been deployed

Reviewers

Uh oh!

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

ParadaCarleton commented Oct 2, 2026 •

edited

Loading

rveltz commented Oct 6, 2026 •

edited

Loading