PALC: optional step control by the corrector's quality and orientation check (off by default) by ParadaCarleton · Pull Request #347 · bifurcationkit/BifurcationKit.jl

ParadaCarleton · 2026-10-03T21:29:05Z

This adds two optional PALC settings for step control, following Allgower and Georg:

step_control = CorrectorQuality(; max_distance = nothing, max_contraction = 0.6, max_growth = 1.1, η = 0.1), which, following Allgower and Georg (Numerical Continuation Methods, 1990) and Gambit's logit tracer, rejects a step (retried at half ds) when the first Newton step is longer than max_distance or a later one contracts by more than max_contraction, and grows ds after an accepted step by a factor set by those same ratios instead of by the Newton iteration count, with nothing switching a test off.
orientation_check = true rejects a step whose new tangent has a negative dot product with the previous tangent (Allgower and Georg 6.1.10), testing that the continuation has not turned back on itself or taken a step so long that the branch bends past a right angle. The distance threshold is in the state's units, which is why that test is nothing unless the caller supplies one, while the contraction rate and the orientation have no units. test/continuation/palc_step_control.jl covers the rules, both settings through the folds of a cubic with the Secant and Bordered predictors in Float32 and Float64, and the orientation test, and News.md has an entry.

Cost of the orientation test (both settings are off by default):

A Secant tangent: the chord, two norms and one dot product per step, no factorization.
With orientation_check or max_turn set, a Bordered tangent is oriented by the sign of the determinant of the bordered matrix (Allgower and Georg, (2.1.7)). The sign comes from the LU factorization that solves for the tangent, at O(n) extra, so a step adds no factorization; the start computes one extra Bordered tangent.
That needs a MatrixBLS solver; with a matrix-free or iterative solver the continuation throws at initialization.

…rQuality) and orientation check, off by default

… makes an angle with the tangent of more than max_angle (off by default)

…Bordered tangent at the new point is computed for the test and reused by the next predictor); max_angle = :derived takes 2 max_distance / |ds| up to max_angle_cap

…ngent is oriented by the determinant of the bordered matrix and a step is rejected when the new tangent has a negative dot product with the previous one

…to the Newton iteration count; a debug record tells whether max_growth capped the growth of ds

… the bound by Newton's method with the corrector quality tests and the orientation test

…son, so a Float type without Int conversion of sign (not Float32) works

rveltz · 2026-10-06T06:48:25Z

 # we also extend this for ContIterable
 getdot(it::ContIterable) = getdot(it.alg)
-getθ(it::ContIterable) = getθ(it.alg)
+getθ(it::ContIterable) = convert(eltype(it), getθ(it.alg))


use oftype

Done in d290a4e: getθ(it::ContIterable) = oftype(it.contparams.ds, getθ(it.alg)) (Palc.jl:114).

rveltz · 2026-10-06T06:49:45Z

    # we want to start at (u0, p0), not at (u1, p1)
    _copyto!(state.z, state.z_old)
+    if alg.orientation_check
+        initial_orientation!(state, iter, alg)


name too generic, only called when orientation check

Renamed init_orientation_check! in d290a4e. It runs only when tracks_turn(alg) holds (orientation_check or a max_turn).

rveltz · 2026-10-06T06:51:54Z

+
+The corrector of a step whose predictor lies beyond `p_min` or `p_max`: the predictor is put on the bound and Newton's method solves `F(x, p) = 0` there (the `Natural` corrector), with the tests of `control` on the lengths of its steps, taken in the norm of the arc length constraint with no parameter component. Returns the solution and its quality like `_newton_palc`.
+"""
+function _newton_at_bound(iter::AbstractContinuationIterable,


I'd say this function is not need. A callback to newton would suffice

You are right. _newton_at_bound is deleted in d290a4e. The corrector quality tests and the orientation test run in a Newton callback that wraps callback_newton, so newton_palc is unchanged by them.

rveltz · 2026-10-06T07:01:56Z

    stepsizecontrol::Bool = true
+    "factor by which step size control multiplies `ds` after the current step. Set by the corrector when the algorithm decides it (for example `PALC` with `step_control`), `nothing` otherwise."
+    step_factor::Union{Nothing, T} = nothing
+    "tangent at the new solution, which the corrector computed to measure the turn of the tangent (`PALC` with `max_angle` and a `Bordered` predictor). The next predictor takes it instead of computing it again. `nothing` otherwise."


this seems too specific to PALC except for

next_tangent
Proposal:

add a field cache::Tcache defaulting to nothing to ContState

define a cache with step_factor, orientation

Done in d290a4e as you proposed. ContState has cache::Tcache = nothing; PALC fills it with a PALCCache (step_factor, orientation, and the data of the corrector quality tests). next_tangent stays a ContState field since it is not specific to PALC. Multiple and AutoSwitch carry the inner PALC cache (6690142).

rveltz · 2026-10-06T07:03:09Z

Hi,

Thank you for this production!

I think you need to take advantage of callbacks for newton and newton_palc instead of modifying the code of newton. If for every new feature, we need to change newton, this is not practical.
CorrectorQuality
- In the doc string, I read Gambit's logit tracer (PathTracer::TracePath, src/solvers/path/path.cc) What is this?
- It is written after Allgower and Georg, can it be more specific? Which algorithm / chapter ?
In the doc string of PALC, I read "Orientation test of Gambit's path tracer (PathTracer::TracePath, src/solvers/path/path.cc):. Same as above.

On a side note:

bordered_determinant_sign is way too expansive. Think about the case of a billion unknowns, or when the jacobian is matrix free. Is it really useful?

…outside stops the corrector, not converged, before F is evaluated there, so no clamped point is reported as converged

ParadaCarleton · 2026-10-11T14:21:35Z

Pushed 479fb7a. newton_palc clamped p into [p_min, p_max] at every iterate, so a corrector whose solution lies past a bound could not converge, the step halved toward the bound without reaching it, and a long clamped corrector could converge on another sheet of the curve and report fake folds. In the legacy run of test/continuation/palc_step_control.jl (below), step 6 runs 15 clamped iterations from u ≈ -0.92 to u ≈ 0.92, skipping the sheet between the folds, and both reported folds lie on the upper sheet, where the curve has none. Now an iterate that leaves [p_min, p_max] stops the corrector, not converged, before F is evaluated there, and the step is halved; the new tests check that the branch visits the middle sheet, that F is never called past p_max, and every other test directory gives the same results as before, with one known pre-existing edge (the finite-difference ∂F/∂p still evaluates F at p + ϵ, just past p_max for an iterate within ϵ of it).

F(u, p) = @. u^3 - u - (p.λ - 1)
prob = BifurcationProblem(F, [-1.3247179], (λ = 0.0,), (@optic _.λ))
opts = ContinuationPar(ds = 0.05, dsmin = 1e-6, dsmax = 0.5, a = 0.5, p_min = 0.0, p_max = 2.0, η = 150.0, tol_stability = 1e-10, dsmin_bisection = 1e-16, tol_bisection_eigenvalue = 1e-16, tol_param_bisection_event = 1e-16, max_steps = 200, detect_fold = true, detect_bifurcation = 0, newton_options = NewtonPar(tol = sqrt(eps())))
br = continuation(prob, PALC(), opts)
any(point -> abs(point.x[1]) < 1 / sqrt(3), br.sol)  # false before 479fb7ad: the branch never enters |u| < 1/√3

…tely typed and the corrector's keyword call dispatches statically

ParadaCarleton · 2026-10-20T15:46:38Z

PALC's corrector! rebound control to nothing during bisection, so it was typed Union{Nothing, CorrectorQuality} and the keyword call to the corrector dispatched at runtime, boxing the iterator, the corrector and the returned solution on every step. This commit splits corrector! into one call per step-control case, so control is concretely typed in _corrector! and the call dispatches statically. Behaviour is unchanged: the same corrector runs with the same arguments. In one downstream continuation workload this removed about 3.4 GB of 20.8 GB of allocations per solve.

rveltz · 2026-10-20T17:09:43Z

    state.itlinear  = sol.itlineartot
    # record previous solution
-    if converged(sol)
+    if converged(state)


Why did you change this?

That change is gone. It was there so that a step refused by the orientation test (converged corrector, tangent pointing back) was recorded as not converged and kept state.z_old. Since d290a4e the orientation test runs in the Newton callback: a refusal stops Newton with converged = false, so sol.converged already says it and _update_field_but_not_solution! is the same as on master. Contbase.jl now differs from master by _step_growth (and the step_size_control! signature), and _step_growth forms the Int ratio (Nmax - itnewton) / Nmax in the type of ds, because the Int quotient is a Float64 and promoted the step of a Float32 continuation.

rveltz · 2026-10-20T17:20:35Z

Thank you!

I thought PALC was already taking into account the clamp?
Did you see my comments above?

…ewton callback wrapping callback_newton _newton and newton_palc return a plain NonLinearSolution again. PALCStepCheck, built per corrector, rejects a step by returning false and leaves the corrector's quality in the PALCCache that ContState.cache carries. Docstrings state the tests and cite Allgower and Georg (6.1.10), Program 1 and (2.1.7); the determinant's O(n^3) cost is documented.

…e an iterate outside [p_min, p_max] refuses the long step that would jump to the upper sheet

…te after a failed line search and stops at an iterate outside [p_min, p_max]; the Bordered tangent takes dFdl from R01 at its point

…rg's step length adaptation (6.1.10): a step whose corrected tangent turns by max_turn is refused and a smaller turn slows the next step

ParadaCarleton · 2026-10-26T20:31:29Z

It does clamp, and the clamp is the problem: PALC clamps only p and leaves x where the unclamped step put it. The resulting point lies off the tangent line, so the corrector can converge onto another branch or sheet. This happens in two places, the predictor in corrector! and each iterate in newton_palc. I moved both fixes out of this PR into PALC: keep the state on the tangent line when p reaches p_min or p_max #363, whose description now explains each case.
Yes, sorry for the delay; I'll answer them inline.

…arch halves a trial past a bound; test wrappers and Float32 Bordered orientation Claude-Session: https://claude.ai/code/session_0195j23m7MRde8V6ApNmKjBR

With `by_determinant = true` and a `MatrixBLS`, `bordered_tangent` builds the bordered matrix once (`bordered_matrix`, shared with `MatrixBLS`), factorises it with `lu`, solves for the tangent and reads the determinant sign from the same factorisation (pivot signs times permutation parity). A sparse matrix takes the sign from its UMFPACK factorisation through `logabsdet`. The orientation test of a `Bordered` tangent then costs no factorisation beyond the tangent's own.

ParadaCarleton · 2026-10-31T08:06:15Z

You were right that bordered_determinant_sign was a second dense O(n^3) factorisation per step. With orientation_check = true or a max_turn (both off by default) and a Bordered tangent, bordered_tangent factorises the bordered matrix once with lu, solves for the tangent with it and reads the determinant sign from the same factors (pivot signs times the row-permutation parity), O(n) extra. It needs a MatrixBLS solver; any other solver throws at initialization. Without those options nothing changes. The Gambit references are gone; the docstrings cite Allgower and Georg (2.1.7, 6.1.10).

This is one commit (cd74ffa4 on my fork), not pushed here. If you agree with the approach, I will push it and update the PR description.

Why the determinant stays: orienting by the dot product with the previous tangent cannot see a single step across a turn of more than π/2 when the chord still looks forward; the determinant sign does. A chord-angle max_turn does not close it either (a true turn of 2 radians gives a chord angle of 1 radian, which max_turn = 1.2 accepts).

ParadaCarleton · 2026-10-31T08:07:46Z

Pushed the single-LU orientation as fab2994 and added its cost to the PR description.

…to tol The loops stop at res <= tol but the verdict required res < tol, so a solve that stopped exactly at tol (including tol = 0 at an exact root) was reported not converged.

… deceleration but rejects no step on it For a Newton callback that judges contraction itself, in its own norm. Claude-Session: https://claude.ai/code/session_014FMUtEXDn4dox6driwcGkm

…n, turn or tangent cosine refuses the step Each test compared with >=, > or <, which a NaN fails, so a NaN passed as an accepted step. Claude-Session: https://claude.ai/code/session_013WmUuZon4D1nexEEjGV1WM

…of reaching the orientation test NaN >= tol is false, so the closing check took a NaN corrector as converged and computed its bordered tangent; LAPACK then threw on the NaN Jacobian. Claude-Session: https://claude.ai/code/session_013WmUuZon4D1nexEEjGV1WM

… stops continuation With tol = 0 two Newton steps of length zero gave the contraction rate 0 / 0. With rejects_contraction = false nothing rejected it, so the deceleration and then ds went NaN. abs(NaN) <= dsmin is false, so every failed step halved ds forever. A step of length zero now has rate zero, and the dsmin test is written as a negation so that a NaN ds stops. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_0195j23kWb8rVLQf4ZxZ1o1R

… tol newton_palc accepts a point with `residual <= tol`, but PALCStepCheck ran the orientation test only when `residual < tol`. At tol = 0, a corrector that lands exactly on a root was accepted without the test, so a step onto a reversed tangent passed. The `<` dates from 5b0f67c (as `>= tol` returns early); af75aa3 kept it while refusing NaN. The guard now uses Newton's own predicate `_converged(residual <= tol, verdict)`, which still skips a NaN residual. Claude-Session: https://claude.ai/code/session_014FMUtEXDn4dox6driwcGkm

ParadaCarleton force-pushed the palc-acceptance branch 2 times, most recently from cfe0341 to 7f26c34 Compare October 4, 2026 05:20

PALC: optional step control by the quality of the corrector (Correcto…

3792c71

…rQuality) and orientation check, off by default

ParadaCarleton force-pushed the palc-acceptance branch from 7f26c34 to 3792c71 Compare October 5, 2026 17:57

ParadaCarleton changed the title ~~PALC: optional Allgower-Georg step acceptance tests and step-size rule (off by default)~~ PALC: optional step control by the corrector's quality and orientation check (off by default) Oct 5, 2026

ParadaCarleton force-pushed the palc-acceptance branch from cc200dd to 3792c71 Compare October 5, 2026 18:02

ParadaCarleton added 9 commits October 5, 2026 13:01

Merge remote-tracking branch 'upstream/master' into palc-acceptance

d1c4711

Merge branch 'float32-disruption' into palc-acceptance

8679685

PALC: max_angle replaces orientation_check: reject a step whose chord…

53aa5cf

… makes an angle with the tangent of more than max_angle (off by default)

PALC: max_angle measures the angle between consecutive tangents (the …

ef90ab2

…Bordered tangent at the new point is computed for the test and reused by the next predictor); max_angle = :derived takes 2 max_distance / |ds| up to max_angle_cap

PALC: orientation_check is Gambit's orientation test: the Bordered ta…

ac59d2c

…ngent is oriented by the determinant of the bordered matrix and a step is rejected when the new tangent has a negative dot product with the previous one

CorrectorQuality: max_growth = nothing keeps the tests and leaves ds …

d2b9bf7

…to the Newton iteration count; a debug record tells whether max_growth capped the growth of ds

PALC: the step whose predictor crosses p_min or p_max is corrected at…

5b0f67c

… the bound by Newton's method with the corrector quality tests and the orientation test

PALC: initial_orientation! takes the sign of the alignment by compari…

16bfd55

…son, so a Float type without Int conversion of sign (not Float32) works

Merge branch 'float32-disruption' into palc-acceptance

df5ba35

rveltz reviewed Oct 6, 2026

View reviewed changes

PALC: newton_palc no longer clamps p into [p_min, p_max]: an iterate …

479fb7a

…outside stops the corrector, not converged, before F is evaluated there, so no clamped point is reported as converged

PALC corrector!: one call per step-control case, so control is concre…

abb192b

…tely typed and the corrector's keyword call dispatches statically

rveltz reviewed Oct 20, 2026

View reviewed changes

ParadaCarleton added 4 commits October 21, 2026 10:56

palc_step_control test: the legacy run passes the middle sheet becaus…

9c44e73

…e an iterate outside [p_min, p_max] refuses the long step that would jump to the upper sheet

Merge master into palc-acceptance: newton_palc keeps master's recompu…

3b03f6a

…te after a failed line search and stops at an iterate outside [p_min, p_max]; the Bordered tangent takes dFdl from R01 at its point

CorrectorQuality: max_turn, the tangent-turn test of Allgower and Geo…

4233774

…rg's step length adaptation (6.1.10): a step whose corrected tangent turns by max_turn is refused and a smaller turn slows the next step

CorrectorQuality: keep the four-argument constructor, max_turn off

4d0d420

ParadaCarleton mentioned this pull request Oct 24, 2026

PALC: keep the state on the tangent line when p reaches p_min or p_max #363

Open

ParadaCarleton added 2 commits October 29, 2026 23:50

PALC: Multiple and AutoSwitch carry the inner PALC cache; the line se…

6690142

…arch halves a trial past a bound; test wrappers and Float32 Bordered orientation Claude-Session: https://claude.ai/code/session_0195j23m7MRde8V6ApNmKjBR

ParadaCarleton and others added 8 commits November 5, 2026 19:08

Merge upstream master

02d351e

Newton, PALC, MoorePenrose and Multiple converge on a residual equal …

9137ffc

…to tol The loops stop at res <= tol but the verdict required res < tol, so a solve that stopped exactly at tol (including tol = 0 at an exact root) was reported not converged.

CorrectorQuality: rejects_contraction = false keeps the contraction's…

24a8578

… deceleration but rejects no step on it For a Newton callback that judges contraction itself, in its own norm. Claude-Session: https://claude.ai/code/session_014FMUtEXDn4dox6driwcGkm

CorrectorQuality and orientation_check: a NaN Newton step, contractio…

8d0cec0

…n, turn or tangent cosine refuses the step Each test compared with >=, > or <, which a NaN fails, so a NaN passed as an accepted step. Claude-Session: https://claude.ai/code/session_013WmUuZon4D1nexEEjGV1WM

Merge tol-inclusive-pr (bifurcationkit#373) into palc-acceptance

19e7157

Conversation

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

Uh oh!

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

Uh oh!

ParadaCarleton commented Oct 11, 2026

Uh oh!

ParadaCarleton commented Oct 20, 2026

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

rveltz commented Oct 20, 2026

Uh oh!

ParadaCarleton commented Oct 26, 2026

Uh oh!

ParadaCarleton commented Oct 31, 2026

Uh oh!

ParadaCarleton commented Oct 31, 2026

Uh oh!

This branch has not been deployed

Reviewers

Uh oh!

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

ParadaCarleton commented Oct 3, 2026 •

edited

Loading