diff --git a/Project.toml b/Project.toml index f78592cb2..94a1ef203 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "NonlinearSolve" uuid = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" -version = "4.21.1" +version = "4.22.0" authors = ["SciML"] [deps] @@ -98,7 +98,7 @@ NLSolvers = "0.5" NLsolve = "4.5" NaNMath = "1" NonlinearProblemLibrary = "0.1.2" -NonlinearSolveBase = "2.35" +NonlinearSolveBase = "2.36" NonlinearSolveFirstOrder = "2" NonlinearSolveQuasiNewton = "1.12" NonlinearSolveSpectralMethods = "1.6" diff --git a/docs/src/native/solvers.md b/docs/src/native/solvers.md index dfb1f9b05..8d976794b 100644 --- a/docs/src/native/solvers.md +++ b/docs/src/native/solvers.md @@ -54,6 +54,7 @@ LimitedMemoryBroyden ```@docs HomotopySweep +KantorovichHomotopy ArcLengthContinuation HomotopyPolyAlgorithm ``` diff --git a/lib/NonlinearSolveBase/Project.toml b/lib/NonlinearSolveBase/Project.toml index ffd7b030a..3d0e42046 100644 --- a/lib/NonlinearSolveBase/Project.toml +++ b/lib/NonlinearSolveBase/Project.toml @@ -1,6 +1,6 @@ name = "NonlinearSolveBase" uuid = "be0214bd-f91f-a760-ac4e-3421ce2b2da0" -version = "2.35.0" +version = "2.36.0" authors = ["Avik Pal and contributors"] [deps] diff --git a/lib/NonlinearSolveBase/src/NonlinearSolveBase.jl b/lib/NonlinearSolveBase/src/NonlinearSolveBase.jl index 502559db6..30452c77e 100644 --- a/lib/NonlinearSolveBase/src/NonlinearSolveBase.jl +++ b/lib/NonlinearSolveBase/src/NonlinearSolveBase.jl @@ -65,6 +65,7 @@ include("timer_outputs.jl") include("tracing.jl") include("wrappers.jl") include("polyalg.jl") +include("kantorovich_homotopy.jl") include("homotopy_sweep.jl") include("arclength.jl") include("homotopy_polyalg.jl") @@ -118,7 +119,7 @@ export DescentResult, SteepestDescent, NewtonDescent, DampedNewtonDescent, Dogle export NonlinearSolvePolyAlgorithm -export HomotopySweep, ArcLengthContinuation, HomotopyPolyAlgorithm +export HomotopySweep, KantorovichHomotopy, ArcLengthContinuation, HomotopyPolyAlgorithm export NonlinearVerbosity diff --git a/lib/NonlinearSolveBase/src/homotopy_polyalg.jl b/lib/NonlinearSolveBase/src/homotopy_polyalg.jl index b71e16bad..d3047d6bb 100644 --- a/lib/NonlinearSolveBase/src/homotopy_polyalg.jl +++ b/lib/NonlinearSolveBase/src/homotopy_polyalg.jl @@ -28,25 +28,27 @@ cost of solving an ``(n+1)``-dimensional corrector system per step. ### Warm handoff -When a [`HomotopySweep`](@ref) stage fails *partway* along the span, everything it -accepted before the failure is genuine converged path: its last accepted iterate is a -solution of ``H(u, λ) = 0`` at some ``λ`` strictly between `λspan[1]` and the failure -point. With `warm_handoff = true` (the default), the next stage is first attempted on -the remaining stretch — the problem is rebuilt with `u0` set to that last accepted -iterate and `λspan` shrunk to `(λ_h, λspan[2])` — instead of redoing the -already-conquered prefix from a cold start at `λspan[1]`. +When a natural-parameter stage ([`HomotopySweep`](@ref) or +[`KantorovichHomotopy`](@ref)) fails *partway* along the span, everything it accepted +before the failure is genuine converged path: its last accepted iterate is a solution +of ``H(u, λ) = 0`` at some ``λ`` strictly between `λspan[1]` and the failure point. With +`warm_handoff = true` (the default), the next stage is first attempted on the remaining +stretch — the problem is rebuilt with `u0` set to that last accepted iterate and +`λspan` shrunk to `(λ_h, λspan[2])` — instead of redoing the already-conquered prefix +from a cold start at `λspan[1]`. The handoff λ is deliberately *backed off* from the failure: `λ_h` is placed 5% of the -span width behind the sweep's last accepted ``λ``. The sweep typically dies at a fold, -where the path turns vertical in ``λ``; a warm stage seeded right at the fold starts -with its initial pure-λ tangent nearly orthogonal to the true path direction and pays -for it in rejected steps (measured on the cubic S-curve: arclength warm-started at the -fold costs *more* residual calls than a full cold run, while backing off 5% costs -~15–25% less). The handed-over `u0` is the last accepted iterate — off-path at `λ_h` -by the backoff distance — and the stage's own λ-fixed anchor solve at `λ_h` pulls it -back onto the path for a few warm Newton iterations. Because the stages measure their -step-size *caps* (`max_step_factor`, and a sweep's fixed-size `nsteps`) as fractions -of the span width, the warm attempt rescales those caps by +span width behind the natural-parameter stage's last accepted ``λ``. Such a stage +typically dies at a fold, where the path turns vertical in ``λ``; a warm stage seeded +right at the fold starts with its initial pure-λ tangent nearly orthogonal to the true +path direction and pays for it in rejected steps (measured on the cubic S-curve: +arclength warm-started at the fold costs *more* residual calls than a full cold run, +while backing off 5% costs ~15–25% less). The handed-over `u0` is the last accepted +iterate — off-path at `λ_h` by the backoff distance — and the stage's own λ-fixed anchor +solve at `λ_h` pulls it back onto the path for a few warm Newton iterations. Because the +stages measure their step-size *caps* (`max_step_factor`, and a natural-parameter +stage's fixed-size `nsteps`) as fractions of the span width, the warm attempt rescales +those caps by `full_width / remaining_width` (capping the fraction at 1, i.e. at an absolute step of the remaining width) so a user-tightened absolute cap survives the span shrink — the initial step factor is left span-relative, since starting small right behind the @@ -56,13 +58,13 @@ polyalgorithm moves on, so enabling the handoff never costs robustness relative `warm_handoff = false` — only, in that rare double-failure case, the extra warm attempt. -The handoff only engages when the sweep made real progress (the backed-off `λ_h` lies -strictly past `λspan[1]`); a sweep that failed at the `λspan[1]` anchor itself, or -within the backoff width of it, leaves the fallback stages with the current cold -full-range behavior. A warm-handoff success is returned as a solution of the -*original* problem (same `prob`, same `u` type); the stage's solution of the shrunken -problem is attached as `original` only when `store_original = Val(true)` (see below) — -by default it is dropped so the returned solution stays concretely typed. +The handoff only engages when the natural-parameter stage made real progress (the +backed-off `λ_h` lies strictly past `λspan[1]`); a stage that failed at the `λspan[1]` +anchor itself, or within the backoff width of it, leaves the fallback stages with the +current cold full-range behavior. A warm-handoff success is returned as a solution of +the *original* problem (same `prob`, same `u` type); the stage's solution of the +shrunken problem is attached as `original` only when `store_original = Val(true)` (see +below) — by default it is dropped so the returned solution stays concretely typed. ### Arguments @@ -72,10 +74,11 @@ by default it is dropped so the returned solution stays concretely typed. ### Keyword Arguments - `warm_handoff`: when `true` (default), a stage following a partway-failed - [`HomotopySweep`](@ref) first attempts the remaining `(λ_h, λspan[2])` stretch - from the sweep's last accepted iterate (with `λ_h` backed off 5% of the span from - the failure), falling back to the cold full-range attempt only if that fails. - `false` recovers the plain try-each-stage-cold behavior. + [`HomotopySweep`](@ref) or [`KantorovichHomotopy`](@ref) first attempts the remaining + `(λ_h, λspan[2])` stretch from the natural-parameter stage's last accepted iterate + (with `λ_h` backed off 5% of the span from the failure), falling back to the cold + full-range attempt only if that fails. `false` recovers the plain + try-each-stage-cold behavior. - `store_original`: whether a warm-handoff success stores its shrunken-problem stage solution in the returned solution's `original` field. Default `Val(false)` keeps the returned solution concretely typed (the stage solution the handoff produces infers @@ -129,9 +132,9 @@ end # rescaled unboundedly vs 146 capped at the remaining width). The *initial* step # factor is NOT rescaled for the same reason — starting small right behind the fold # is measurably cheaper, and adaptive growth recovers the size within a few accepted -# steps. A sweep stage's fixed-size `nsteps` is rescaled like the cap: `nsteps` of -# the *remaining* span would shrink the absolute increment (and `adaptive = false` -# sweeps cannot recover from that). Stages of unknown type get the shrunken problem +# steps. A natural-parameter stage's fixed-size `nsteps` is rescaled like the cap: +# `nsteps` of the *remaining* span would shrink the absolute increment (and a fixed-step +# sweep cannot recover from that). Stages of unknown type get the shrunken problem # unchanged — no field contract to rescale against. _rescale_step_caps(stage, scale) = stage function _rescale_step_caps(stage::ArcLengthContinuation, scale) @@ -145,6 +148,13 @@ function _rescale_step_caps(stage::HomotopySweep, scale) msf = stage.max_step_factor return @set stage.max_step_factor = min(msf * scale, oneunit(msf)) end +function _rescale_step_caps(stage::KantorovichHomotopy, scale) + if stage.nsteps !== nothing + stage = @set stage.nsteps = max(1, ceil(Int, stage.nsteps / scale)) + end + msf = stage.max_step_factor + return @set stage.max_step_factor = min(msf * scale, oneunit(msf)) +end function CommonSolve.solve( prob::SciMLBase.HomotopyProblem{uType, iip}, alg::HomotopyPolyAlgorithm, @@ -154,8 +164,8 @@ function CommonSolve.solve( throw(ArgumentError("HomotopyPolyAlgorithm requires at least one algorithm")) nstages = length(alg.algs) λ0, λ1 = prob.λspan - # Handoff seed `(u_last, λ_h)` from the most recent partway-failed sweep stage, - # or `nothing` when none is available: `u_last` is the sweep's last accepted + # Handoff seed `(u_last, λ_h)` from the most recent partway-failed natural-parameter + # stage, or `nothing` when none is available: `u_last` is the stage's last accepted # iterate and `λ_h` its λ backed off by 5% of the span (see the loop below). handoff = nothing for (i, stage) in enumerate(alg.algs) @@ -186,21 +196,21 @@ function CommonSolve.solve( # through to the cold full-range attempt so the handoff never costs # robustness relative to the plain staged behavior. end - sol = if stage isa HomotopySweep + sol = if stage isa Union{HomotopySweep, KantorovichHomotopy} csol, λ_last = _homotopy_sweep_solve(prob, stage, args...; kwargs...) if λ_last !== nothing - # The sweep dies where the path turns hard (a fold), and a fallback - # seeded exactly there starts with its initial pure-λ tangent nearly - # orthogonal to the path — measured on the cubic S-curve, arclength - # warm-started AT the fold costs more residual calls than a full - # cold run, while 5% of the span behind it costs ~15–25% less. So - # the handoff λ is backed off by span/20 from the last accepted - # point; the stage's own λ-fixed anchor solve pulls `u_last` (5% of - # the span away in λ) back onto the path at `λ_h` for a few - # Newton iterations. The handoff only engages when the backed-off λ - # is still strictly past `λspan[1]` — a sweep that died within the - # backoff width of the anchor leaves the fallback no cheaper start - # than its own anchor solve. + # A natural-parameter stage dies where the path turns hard (a fold), + # and a fallback seeded exactly there starts with its initial pure-λ + # tangent nearly orthogonal to the path — measured on the cubic + # S-curve, arclength warm-started AT the fold costs more residual calls + # than a full cold run, while 5% of the span behind it costs ~15–25% + # less. So the handoff λ is backed off by span/20 from the last accepted + # point; the stage's own λ-fixed anchor solve pulls `u_last` (5% of the + # span away in λ) back onto the path at `λ_h` for a few Newton + # iterations. The handoff only engages when the backed-off λ is still + # strictly past `λspan[1]` — a stage that died within the backoff width + # of the anchor leaves the fallback no cheaper start than its own anchor + # solve. backoff = (oftype(λ_last, λ1) - oftype(λ_last, λ0)) / 20 if abs(λ_last - oftype(λ_last, λ0)) > abs(backoff) handoff = (csol.u, λ_last - backoff) diff --git a/lib/NonlinearSolveBase/src/homotopy_sweep.jl b/lib/NonlinearSolveBase/src/homotopy_sweep.jl index 5f95b8101..b0a1a02af 100644 --- a/lib/NonlinearSolveBase/src/homotopy_sweep.jl +++ b/lib/NonlinearSolveBase/src/homotopy_sweep.jl @@ -380,6 +380,137 @@ function _effort_wants_shrink(nit::Int, budget::Int) return nit >= 0 && nit <= budget && 4 * nit >= 3 * budget end +function _homotopy_corrector!(cache, ::HomotopySweep, ::Type{T}) where {T} + sol = CommonSolve.solve!(cache) + return sol, zero(T), zero(T), false, false +end + +function _homotopy_corrector!( + cache::NonlinearSolvePolyAlgorithmCache, ::KantorovichHomotopy, ::Type{T} + ) where {T} + sol = CommonSolve.solve!(cache) + return sol, zero(T), zero(T), false, false +end + +function _homotopy_corrector!( + cache, alg::KantorovichHomotopy, ::Type{T} + ) where {T} + # No-init algorithms do not expose intermediate residuals or a step interface. + if !hasfield(typeof(cache), :fu) || cache.retcode == ReturnCode.InitialFailure + sol = CommonSolve.solve!(cache) + return sol, zero(T), zero(T), false, false + end + + residualnormprev = L2_NORM(get_fu(cache)) + first_Θ = zero(T) + rejected_Θ = zero(T) + has_Θ = false + contraction_rejected = false + Θreject = T(alg.Θreject) + + while not_terminated(cache) + CommonSolve.step!(cache) + residualnorm = L2_NORM(get_fu(cache)) + Θ = iszero(residualnormprev) ? zero(T) : T(residualnorm / residualnormprev) + if !has_Θ + first_Θ = Θ + has_Θ = true + end + if !contraction_rejected && !(Θ <= Θreject) + rejected_Θ = Θ + contraction_rejected = true + end + residualnormprev = residualnorm + end + + sol = CommonSolve.solve!(cache) + return sol, first_Θ, rejected_Θ, has_Θ, contraction_rejected +end + + +function _kantorovich_step_factor( + alg::KantorovichHomotopy, Θ, has_Θ::Bool, ::Type{T} + ) where {T} + has_Θ && (!isfinite(Θ) || Θ < 0) && return T(alg.qmin) + Θeff = has_Θ ? max(T(Θ), T(alg.Θmin)) : T(alg.Θmin) + g(x) = sqrt(one(T) + 4 * x) - one(T) + q = T(alg.γ) * (g(T(alg.Θbar)) / g(Θeff))^(inv(T(alg.predictor_order))) + return clamp(q, T(alg.qmin), T(alg.qmax)) +end + +function _accepted_step_size( + alg::HomotopySweep, dλ, streak::Int, θ, nit::Int, budget::Int, + min_dλ, max_dλ, expand_factor, ::Type{T}, first_Θ, has_Θ + ) where {T} + alg.adaptive || return dλ, streak + if _effort_wants_shrink(nit, budget) + abs(dλ) / 2 >= min_dλ && (dλ = dλ / 2) + return dλ, 0 + end + + streak += 1 + corrector_cheap = nit >= 0 && nit <= 2 + if streak >= alg.expand_threshold && + (θ === nothing || θ <= T(alg.expand_quality) || corrector_cheap) + growth = _effort_growth_factor(nit, budget, expand_factor) + if growth > 1 + grown = growth * dλ + dλ = abs(grown) > abs(max_dλ) ? max_dλ : grown + streak = 0 + end + end + return dλ, streak +end + +function _accepted_step_size( + alg::KantorovichHomotopy, dλ, streak::Int, θ, nit::Int, budget::Int, + min_dλ, max_dλ, expand_factor, ::Type{T}, first_Θ, has_Θ + ) where {T} + q = _kantorovich_step_factor(alg, first_Θ, has_Θ, T) + corrector_cheap = nit >= 0 && nit <= 2 + if q > 1 && !(θ === nothing || θ <= T(alg.expand_quality) || corrector_cheap) + q = one(T) + end + proposed = q * dλ + if abs(proposed) > abs(max_dλ) + proposed = max_dλ + elseif abs(proposed) < min_dλ + proposed = sign(dλ) * min_dλ + end + return proposed, 0 +end + +function _rejected_step_size( + alg::HomotopySweep, dλ, attempted_dλ, min_dλ, + first_Θ, rejected_Θ, has_Θ, contraction_rejected, inner_success, + ::Type{T} + ) where {T} + alg.adaptive || return dλ, false + proposed = dλ / 2 + return proposed, abs(proposed) >= min_dλ +end + +function _rejected_step_size( + alg::KantorovichHomotopy, dλ, attempted_dλ, min_dλ, + first_Θ, rejected_Θ, has_Θ, contraction_rejected, inner_success, + ::Type{T} + ) where {T} + q = if inner_success && contraction_rejected && isfinite(rejected_Θ) && rejected_Θ >= 0 + _kantorovich_step_factor(alg, rejected_Θ, true, T) + else + T(alg.qmin) + end + proposed = q * attempted_dλ + return proposed, abs(proposed) >= min_dλ && !iszero(proposed) +end + +_strict_contraction_rejection(::HomotopySweep, contraction_rejected) = false +function _strict_contraction_rejection( + alg::KantorovichHomotopy, contraction_rejected + ) + return alg.strict && contraction_rejected +end + # Full-fidelity standalone solve at a fixed λ, with the user's original kwargs (full # iteration budget, full tolerances) outside the tracking cache. Used only where the # tracking cap / loose tracking tolerance are active but must not bind: the λspan[1] @@ -396,7 +527,8 @@ function _sweep_exempt_solve( end function CommonSolve.solve( - prob::SciMLBase.HomotopyProblem, alg::HomotopySweep, args...; kwargs... + prob::SciMLBase.HomotopyProblem, + alg::Union{HomotopySweep, KantorovichHomotopy}, args...; kwargs... ) sol, _ = _homotopy_sweep_solve(prob, alg, args...; kwargs...) return sol @@ -410,7 +542,7 @@ end # fallback cold at `λspan[1]`. function _homotopy_sweep_solve( prob::SciMLBase.HomotopyProblem{uType, iip}, - alg::HomotopySweep, args...; kwargs... + alg::Union{HomotopySweep, KantorovichHomotopy}, args...; kwargs... ) where {uType, iip} λ0, λ1 = prob.λspan λ = float(λ0) @@ -420,7 +552,7 @@ function _homotopy_sweep_solve( dλ = alg.nsteps === nothing ? λT(alg.initial_step_factor) * span : span / alg.nsteps min_dλ = alg.min_dλ === nothing ? sqrt(eps(λT)) : λT(alg.min_dλ) max_dλ = λT(alg.max_step_factor) * span # carries span's sign, like dλ - expand_factor = λT(alg.expand_factor) + expand_factor = alg isa HomotopySweep ? λT(alg.expand_factor) : one(λT) abs(dλ) > abs(max_dλ) && (dλ = max_dλ) u = copy(prob.u0) budget, cap_active = _tracking_budget(alg.tracking_maxiters, prob.kwargs, kwargs) @@ -553,7 +685,8 @@ function _homotopy_sweep_solve( # won the previous solve (the anchor's full-ladder run discovers the winner) # instead of re-failing the cheaper ladder members on every warm-started step. reinit_retaining!(cache, guess) - last_sol = CommonSolve.solve!(cache) + last_sol, first_Θ, rejected_Θ, has_Θ, contraction_rejected = + _homotopy_corrector!(cache, alg, λT) if next_λ == λend if cap_active && !SciMLBase.successful_retcode(last_sol) @@ -571,6 +704,10 @@ function _homotopy_sweep_solve( last_sol = _sweep_exempt_solve( prob, alg.inner, retry_guess, next_λ, args...; kwargs... ) + first_Θ = zero(λT) + rejected_Θ = zero(λT) + has_Θ = false + contraction_rejected = false elseif tol_active && SciMLBase.successful_retcode(last_sol) # The landing is exempt from the loose tracking tolerance: the loose # cache solve above did the bulk of the convergence, now re-polish at @@ -583,7 +720,9 @@ function _homotopy_sweep_solve( end end - if SciMLBase.successful_retcode(last_sol) + inner_success = SciMLBase.successful_retcode(last_sol) + strict_rejection = _strict_contraction_rejection(alg, contraction_rejected) + if inner_success && !strict_rejection # The secant prediction error θ (relative to the recent solution movement) # is a cheap local error estimate: it grows with the path's curvature times # dλ², so a large θ means the path is turning and a stale tangent would @@ -615,52 +754,32 @@ function _homotopy_sweep_solve( λ_prev = λ λ = next_λ λ == λend && break - if alg.adaptive - nit = last_sol.stats === nothing ? -1 : Int(last_sol.stats.nsteps) - if _effort_wants_shrink(nit, budget) - # Proactive shrink on a straining success (see - # `_effort_wants_shrink`); the floor guard keeps the increment - # from dropping below what bisection itself may reach. - abs(dλ) / 2 >= min_dλ && (dλ = dλ / 2) - streak = 0 - else - streak += 1 - # Growth requires a streak of successes (the classic heuristic) - # plus evidence the corrector has headroom: either a small relative - # prediction error (the quality gate) or a corrector that converged - # almost immediately. The iteration count covers paths that flatten - # exponentially, where θ stays at a constant mediocre value while - # the absolute corrections — and hence the corrector work — become - # negligible. The growth factor itself is scaled by the corrector's - # effort (see `_effort_growth_factor`), so the gates veto growth - # and the effort bands size it. The streak is NOT reset on a vetoed - # or held expansion, so growth resumes on the first step whose - # evidence recovers. - corrector_cheap = nit >= 0 && nit <= 2 - if streak >= alg.expand_threshold && - (θ === nothing || θ <= λT(alg.expand_quality) || corrector_cheap) - g = _effort_growth_factor(nit, budget, expand_factor) - if g > 1 - grown = g * dλ - dλ = abs(grown) > abs(max_dλ) ? max_dλ : grown - streak = 0 - end - end - end - end - elseif alg.adaptive && abs(dλ) / 2 >= min_dλ - dλ = dλ / 2 # bisect; retry from the same λ (do not advance) - streak = 0 - # a rejected step is evidence against the tangent: bisection retries (and - # the steps right after) warm-start constantly until quality re-accumulates - trust = 0 + nit = last_sol.stats === nothing ? -1 : Int(last_sol.stats.nsteps) + dλ, streak = _accepted_step_size( + alg, dλ, streak, θ, nit, budget, min_dλ, max_dλ, + expand_factor, λT, first_Θ, has_Θ + ) else - # on failure: u is the last converged iterate (λ<λ1); resid is from the failed step (advisory) - return build_solution_less_specialize( - prob, alg, u, last_sol.resid; - retcode = last_sol.retcode, original = last_sol, - store_original = alg.store_original - ), λ + attempted_dλ = next_λ - λ + proposed_dλ, can_retry = _rejected_step_size( + alg, dλ, attempted_dλ, min_dλ, first_Θ, rejected_Θ, + has_Θ, contraction_rejected, inner_success, λT + ) + if can_retry + dλ = proposed_dλ + streak = 0 + # a rejected step is evidence against the tangent: bisection retries (and + # the steps right after) warm-start constantly until quality re-accumulates + trust = 0 + else + # on failure: u is the last converged iterate (λ<λ1); resid is from the failed step (advisory) + retcode = strict_rejection ? ReturnCode.ConvergenceFailure : last_sol.retcode + return build_solution_less_specialize( + prob, alg, u, last_sol.resid; + retcode, original = last_sol, + store_original = alg.store_original + ), λ + end end end diff --git a/lib/NonlinearSolveBase/src/kantorovich_homotopy.jl b/lib/NonlinearSolveBase/src/kantorovich_homotopy.jl new file mode 100644 index 000000000..fcb4818a9 --- /dev/null +++ b/lib/NonlinearSolveBase/src/kantorovich_homotopy.jl @@ -0,0 +1,178 @@ +@doc raw""" + KantorovichHomotopy(; inner = nothing, nsteps = nothing, + initial_step_factor = 0.1, min_dλ = nothing, max_step_factor = 1.0, + qmin = 1 // 5, qmax = 5, Θmin = 1 // 8, Θreject = 0.95, + Θbar = 0.5, γ = 0.95, strict = true, predictor = :constant, + predictor_order = nothing, expand_quality = 0.25, + tracking_maxiters = 10, tracking_abstol = nothing, maxsteps = 10000, + store_original = Val(false)) + +Natural-parameter continuation for a [`SciMLBase.HomotopyProblem`](@ref), with step +sizes chosen from the observed contraction of the inner nonlinear corrector. This is +the Newton--Kantorovich path-following controller described in Section 5.1.3 of +Deuflhard, *Newton Methods for Nonlinear Problems*. + +For each accepted continuation point, the solver measures residual contraction ratios + +```math +\Theta_k = \frac{\lVert H(u^{k + 1}, \lambda)\rVert} + {\lVert H(u^k, \lambda)\rVert} +``` + +during the corrector. If ``\Theta_0`` is the first available ratio, the next parameter +increment is multiplied by + +```math +q = \operatorname{clamp}\left( + \gamma \left[\frac{g(\bar\Theta)}{g(\max(\Theta_0, \Theta_{min}))}\right]^{1/p}, + q_{min}, q_{max}\right), \qquad g(x) = \sqrt{1 + 4x} - 1. +``` + +Here ``p`` is the predictor order: by default `1` for `predictor = :constant` and `2` +for `predictor = :secant`. The constant predictor matches `ImplicitDiscreteSolve`'s +controller and is the default. A corrector whose contraction exceeds `Θreject` is rejected +when `strict = true`, even if the inner solver eventually converged, and is retried with +the smaller increment prescribed by the same formula. + +The driver shares the cache reuse, secant trust monitoring, prediction-quality growth +gate, interior iteration cap, optional loose tracking tolerance, analytic/sparse +Jacobian forwarding, and endpoint polishing of [`HomotopySweep`](@ref). Compared with +`HomotopySweep`, only the parameter-step controller differs: `HomotopySweep` uses +success streaks and coarse corrector-effort bands, while `KantorovichHomotopy` uses the +measured contraction rate after every corrector. + +Keyword arguments: + + - `inner`: the inner nonlinear algorithm. `nothing` selects NonlinearSolve's default + polyalgorithm. A polyalgorithm contains separate corrector sequences whose residual + contractions are not comparable across rungs, and algorithms without an iterative + cache do not expose intermediate residuals, so those solves use `Θmin`. Pass a single + cache-based iterative algorithm such as `NewtonRaphson()` to activate the + measured-contraction controller. + - `nsteps`: optional number of equal divisions used only to choose the initial + parameter increment. Subsequent increments remain adaptive. + - `initial_step_factor`: initial increment as a fraction of the `λspan` width when + `nsteps` is not supplied. + - `min_dλ`, `max_step_factor`: minimum absolute increment and maximum increment as a + fraction of the span width. `min_dλ = nothing` resolves to + `sqrt(eps(typeof(λ)))`. + - `qmin`, `qmax`: lower and upper bounds for the step-size multiplier. + - `Θmin`: contraction-rate floor used when a corrector converges before a ratio can be + measured or when the measured contraction is smaller than the floor. + - `Θbar`: target corrector contraction rate. + - `Θreject`: contraction rate above which a corrector is rejected in strict mode. + - `γ`: safety factor in the Kantorovich step-size formula. + - `strict`: reject converged correctors containing any contraction ratio greater than + `Θreject`. With `false`, every converged corrector is accepted and the ratio only + controls the following increment. + - `predictor`: `:secant` or `:constant`, with the same trust-monitored behavior as + [`HomotopySweep`](@ref). + - `predictor_order`: exponent denominator ``p`` in the controller formula. `nothing` + selects `2` for the secant predictor and `1` for the constant predictor. + - `expand_quality`: maximum relative secant-prediction error that permits `q > 1`. + A corrector taking at most two iterations also permits growth. `Inf` disables this + gate. + - `tracking_maxiters`, `tracking_abstol`, `maxsteps`, `store_original`: identical to + the corresponding [`HomotopySweep`](@ref) options. + +This algorithm follows `λ` monotonically and therefore cannot round a fold. Combine it +with [`ArcLengthContinuation`](@ref) in a [`HomotopyPolyAlgorithm`](@ref) when a fold +must be traversed. +""" +@concrete struct KantorovichHomotopy <: AbstractNonlinearSolveAlgorithm + inner + nsteps + initial_step_factor + min_dλ + max_step_factor + qmin + qmax + Θmin + Θreject + Θbar + γ + strict::Bool + predictor::Symbol + predictor_order::Int + expand_quality + tracking_maxiters + tracking_abstol + maxsteps::Int + store_original <: Val +end + +function KantorovichHomotopy(; + inner = nothing, nsteps = nothing, initial_step_factor = 0.1, + min_dλ = nothing, max_step_factor = 1.0, qmin = 1 // 5, qmax = 5, + Θmin = 1 // 8, Θreject = 0.95, Θbar = 0.5, γ = 0.95, + strict = true, predictor = :constant, predictor_order = nothing, + expand_quality = 0.25, tracking_maxiters = 10, tracking_abstol = nothing, + maxsteps = 10000, store_original::Val = Val(false) + ) + nsteps !== nothing && nsteps < 1 && + throw(ArgumentError("KantorovichHomotopy `nsteps` must be ≥ 1, got $nsteps")) + !(0 < initial_step_factor <= 1) && throw( + ArgumentError( + "KantorovichHomotopy `initial_step_factor` must be in (0, 1], " * + "got $initial_step_factor" + ) + ) + min_dλ !== nothing && min_dλ <= 0 && throw( + ArgumentError("KantorovichHomotopy `min_dλ` must be positive, got $min_dλ") + ) + !(0 < max_step_factor <= 1) && throw( + ArgumentError( + "KantorovichHomotopy `max_step_factor` must be in (0, 1], " * + "got $max_step_factor" + ) + ) + !(0 < qmin < 1) && + throw(ArgumentError("KantorovichHomotopy `qmin` must be in (0, 1), got $qmin")) + qmax < 1 && + throw(ArgumentError("KantorovichHomotopy `qmax` must be ≥ 1, got $qmax")) + !(0 < Θmin <= Θbar < Θreject < 1) && throw( + ArgumentError( + "KantorovichHomotopy requires 0 < Θmin ≤ Θbar < Θreject < 1; " * + "got Θmin = $Θmin, Θbar = $Θbar, Θreject = $Θreject" + ) + ) + !(0 < γ < 1) && + throw(ArgumentError("KantorovichHomotopy `γ` must be in (0, 1), got $γ")) + predictor in (:secant, :constant) || throw( + ArgumentError( + "KantorovichHomotopy `predictor` must be :secant or :constant, " * + "got :$predictor" + ) + ) + predictor_order = something(predictor_order, predictor === :secant ? 2 : 1) + predictor_order < 1 && throw( + ArgumentError( + "KantorovichHomotopy `predictor_order` must be ≥ 1, got $predictor_order" + ) + ) + !(expand_quality > 0) && throw( + ArgumentError( + "KantorovichHomotopy `expand_quality` must be positive, got $expand_quality" + ) + ) + tracking_maxiters !== nothing && tracking_maxiters < 1 && throw( + ArgumentError( + "KantorovichHomotopy `tracking_maxiters` must be ≥ 1, got " * + "$tracking_maxiters" + ) + ) + tracking_abstol !== nothing && !(tracking_abstol > 0) && throw( + ArgumentError( + "KantorovichHomotopy `tracking_abstol` must be positive, got " * + "$tracking_abstol" + ) + ) + maxsteps < 1 && + throw(ArgumentError("KantorovichHomotopy `maxsteps` must be ≥ 1, got $maxsteps")) + + return KantorovichHomotopy( + inner, nsteps, initial_step_factor, min_dλ, max_step_factor, qmin, qmax, + Θmin, Θreject, Θbar, γ, strict, predictor, predictor_order, + expand_quality, tracking_maxiters, tracking_abstol, maxsteps, store_original + ) +end diff --git a/test/Core/homotopy_sweep_tests__item22.jl b/test/Core/homotopy_sweep_tests__item22.jl index 4d5db546e..32a1f633a 100644 --- a/test/Core/homotopy_sweep_tests__item22.jl +++ b/test/Core/homotopy_sweep_tests__item22.jl @@ -38,7 +38,7 @@ prob = HomotopyProblem(H!, [4.0], [4.0]) # By default the drivers drop `original`, so their returned solution is concretely typed # (its `original` slot is `Nothing`, not `Any`), and the full `solve` no longer infers as # `Any`. -for alg in (HomotopySweep(), ArcLengthContinuation()) +for alg in (HomotopySweep(), KantorovichHomotopy(), ArcLengthContinuation()) sol = solve(prob, alg) @test SciMLBase.successful_retcode(sol) @test sol.u[1] ≈ 2.0 atol = 1.0e-6 @@ -52,9 +52,11 @@ end # The `store_original` field is carried and defaults to `Val(false)`. @test HomotopySweep().store_original === Val(false) +@test KantorovichHomotopy().store_original === Val(false) @test ArcLengthContinuation().store_original === Val(false) @test HomotopyPolyAlgorithm().store_original === Val(false) @test HomotopySweep(; store_original = Val(true)).store_original === Val(true) +@test KantorovichHomotopy(; store_original = Val(true)).store_original === Val(true) @test ArcLengthContinuation(; store_original = Val(true)).store_original === Val(true) @test HomotopyPolyAlgorithm(; store_original = Val(true)).store_original === Val(true) diff --git a/test/Core/kantorovich_homotopy_tests__item1.jl b/test/Core/kantorovich_homotopy_tests__item1.jl new file mode 100644 index 000000000..788a08a9a --- /dev/null +++ b/test/Core/kantorovich_homotopy_tests__item1.jl @@ -0,0 +1,193 @@ +using NonlinearSolve +using CommonSolve +using SciMLBase +using Test + +alg = KantorovichHomotopy() +@test alg.predictor === :constant +@test alg.predictor_order == 1 +@test alg.strict +@test alg.qmin == 1 // 5 +@test alg.qmax == 5 +@test KantorovichHomotopy(; predictor = :secant).predictor_order == 2 +@test KantorovichHomotopy(; predictor = :secant, predictor_order = 3).predictor_order == 3 +@test NonlinearSolveBase._kantorovich_step_factor(alg, NaN, true, Float64) == 0.2 + +@test_throws ArgumentError KantorovichHomotopy(; nsteps = 0) +@test_throws ArgumentError KantorovichHomotopy(; initial_step_factor = 0) +@test_throws ArgumentError KantorovichHomotopy(; min_dλ = 0) +@test_throws ArgumentError KantorovichHomotopy(; max_step_factor = 2) +@test_throws ArgumentError KantorovichHomotopy(; qmin = 1) +@test_throws ArgumentError KantorovichHomotopy(; qmax = 0.9) +@test_throws ArgumentError KantorovichHomotopy(; Θmin = 0.6, Θbar = 0.5) +@test_throws ArgumentError KantorovichHomotopy(; Θbar = 0.96) +@test_throws ArgumentError KantorovichHomotopy(; γ = 1) +@test_throws ArgumentError KantorovichHomotopy(; predictor = :quadratic) +@test_throws ArgumentError KantorovichHomotopy(; predictor_order = 0) +@test_throws ArgumentError KantorovichHomotopy(; expand_quality = 0) +@test_throws ArgumentError KantorovichHomotopy(; tracking_maxiters = 0) +@test_throws ArgumentError KantorovichHomotopy(; tracking_abstol = 0) +@test_throws ArgumentError KantorovichHomotopy(; maxsteps = 0) + +H(u, p, λ) = [u[1]^3 - (1 + λ)] +prob = HomotopyProblem(H, [1.0]; λspan = (0.0, 1.0)) +sol = solve(prob, KantorovichHomotopy(; inner = NewtonRaphson())) +@test SciMLBase.successful_retcode(sol) +@test sol.u[1] ≈ cbrt(2) atol = 1.0e-9 + +soldefault = solve(prob, KantorovichHomotopy()) +@test SciMLBase.successful_retcode(soldefault) +@test soldefault.u[1] ≈ cbrt(2) atol = 1.0e-9 + +solnocache = solve(prob, KantorovichHomotopy(; inner = SimpleNewtonRaphson())) +@test SciMLBase.successful_retcode(solnocache) +@test solnocache.u[1] ≈ cbrt(2) atol = 1.0e-9 + +function Hiip!(du, u, p, λ) + du[1] = u[1]^2 - (1 + λ) + return nothing +end +probiip = HomotopyProblem{true}(Hiip!, Float32[1]; λspan = (0.0f0, 1.0f0)) +soliip = solve(probiip, KantorovichHomotopy(; inner = NewtonRaphson())) +@test SciMLBase.successful_retcode(soliip) +@test eltype(soliip.u) == Float32 +@test soliip.u[1] ≈ sqrt(2.0f0) atol = 2.0f-5 + +Hdecreasing(u, p, λ) = [u[1] - (1 + λ)] +probdecreasing = HomotopyProblem(Hdecreasing, [2.0]; λspan = (1.0, 0.0)) +soldecreasing = solve( + probdecreasing, KantorovichHomotopy(; inner = NewtonRaphson()) +) +@test SciMLBase.successful_retcode(soldecreasing) +@test soldecreasing.u[1] ≈ 1.0 atol = 1.0e-10 + +λcalls = Float64[] +function Htracked(u, p, λ) + push!(λcalls, λ) + return [u[1] - λ] +end +probtracked = HomotopyProblem(Htracked, [0.0]; λspan = (0.0, 1.0)) +soltracked = solve( + probtracked, + KantorovichHomotopy(; + inner = NewtonRaphson(), initial_step_factor = 0.1, max_step_factor = 1.0 + ) +) +@test SciMLBase.successful_retcode(soltracked) +λattempts = [λcalls[i] for i in eachindex(λcalls) if i == 1 || λcalls[i] != λcalls[i - 1]] +q_expected = 0.95 * + ((sqrt(1 + 4 * 0.5) - 1) / (sqrt(1 + 4 * (1 / 8)) - 1)) +@test λattempts[1:3] ≈ [0.0, 0.1, 0.1 * (1 + q_expected)] atol = 1.0e-12 + +struct SlowContractionAlgorithm <: SciMLBase.AbstractNonlinearAlgorithm end + +mutable struct SlowContractionCache{U, P, A, S} <: + NonlinearSolveBase.AbstractNonlinearSolveCache + u::U + fu::U + prob::P + alg::A + stats::S + nsteps::Int + maxiters::Int + force_stop::Bool + retcode::SciMLBase.ReturnCode.T +end + +function SciMLBase.__init( + prob::SciMLBase.AbstractNonlinearProblem, alg::SlowContractionAlgorithm, args...; + maxiters = 3, kwargs... + ) + u = copy(prob.u0) + fu = prob.f(u, prob.p) + solved = iszero(NonlinearSolveBase.L2_NORM(fu)) + return SlowContractionCache( + u, fu, prob, alg, nothing, 0, maxiters, solved, + solved ? ReturnCode.Success : ReturnCode.Default + ) +end + +function SciMLBase.reinit!(cache::SlowContractionCache, u0; kwargs...) + copyto!(cache.u, u0) + cache.fu = cache.prob.f(cache.u, cache.prob.p) + cache.nsteps = 0 + cache.force_stop = iszero(NonlinearSolveBase.L2_NORM(cache.fu)) + cache.retcode = cache.force_stop ? ReturnCode.Success : ReturnCode.Default + return cache +end + +function CommonSolve.step!(cache::SlowContractionCache; kwargs...) + cache.nsteps += 1 + cache.u .-= 0.01 .* cache.fu + cache.fu = cache.prob.f(cache.u, cache.prob.p) + if cache.nsteps == cache.maxiters + cache.u .-= cache.fu + cache.fu = cache.prob.f(cache.u, cache.prob.p) + cache.retcode = ReturnCode.Success + cache.force_stop = true + end + return cache +end + +function CommonSolve.solve!(cache::SlowContractionCache) + while NonlinearSolveBase.not_terminated(cache) + CommonSolve.step!(cache) + end + cache.retcode == ReturnCode.Default && (cache.retcode = ReturnCode.MaxIters) + return SciMLBase.build_solution( + cache.prob, cache.alg, copy(cache.u), copy(cache.fu); retcode = cache.retcode + ) +end + +corrector_prob = NonlinearProblem((u, p) -> [u[1] - 1], [0.0]) +corrector_cache = init(corrector_prob, SlowContractionAlgorithm(); maxiters = 2) +corrector_sol, first_Θ, rejected_Θ, has_Θ, contraction_rejected = + NonlinearSolveBase._homotopy_corrector!(corrector_cache, alg, Float64) +@test SciMLBase.successful_retcode(corrector_sol) +@test has_Θ +@test first_Θ ≈ 0.99 +@test rejected_Θ ≈ 0.99 +@test contraction_rejected + +slowprob = HomotopyProblem((u, p, λ) -> [u[1] - λ], [0.0]; λspan = (0.0, 0.1)) +strictsol = solve( + slowprob, + KantorovichHomotopy(; + inner = SlowContractionAlgorithm(), initial_step_factor = 1.0, + min_dλ = 1.0e-3, strict = true + ) +) +@test strictsol.retcode == ReturnCode.ConvergenceFailure +@test strictsol.u == [0.0] + +nonstrictsol = solve( + slowprob, + KantorovichHomotopy(; + inner = SlowContractionAlgorithm(), initial_step_factor = 1.0, + min_dλ = 1.0e-3, strict = false + ) +) +@test SciMLBase.successful_retcode(nonstrictsol) +@test nonstrictsol.u ≈ [0.1] atol = 1.0e-12 + +foldtarget = 2.1038034 +Hfold(u, p, λ) = [u[1]^3 - 3 * u[1] - (-3 + 6λ)] +foldprob = HomotopyProblem(Hfold, [-foldtarget]; λspan = (0.0, 1.0)) +foldstage = KantorovichHomotopy(; + inner = NewtonRaphson(), min_dλ = 1.0e-2 +) +foldstage_sol = solve(foldprob, foldstage; maxiters = 10) +@test !SciMLBase.successful_retcode(foldstage_sol) + +foldpoly_sol = solve( + foldprob, + HomotopyPolyAlgorithm( + (foldstage, ArcLengthContinuation(; inner = SimpleNewtonRaphson())); + store_original = Val(true) + ); + maxiters = 10 +) +@test SciMLBase.successful_retcode(foldpoly_sol) +@test foldpoly_sol.u[1] ≈ foldtarget atol = 1.0e-4 +@test foldpoly_sol.original !== nothing +@test first(foldpoly_sol.original.prob.λspan) > first(foldprob.λspan) diff --git a/test/runtests.jl b/test/runtests.jl index e67c49a79..837ab996e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -87,6 +87,7 @@ else @time @safetestset "HomotopySweep secant predictor reduces corrector work" include("Core/homotopy_sweep_tests__item20.jl") @time @safetestset "HomotopySweep regrows the step after bisecting a hard region" include("Core/homotopy_sweep_tests__item21.jl") @time @safetestset "Continuation drivers return concretely-typed solutions (no Any-typed original)" include("Core/homotopy_sweep_tests__item22.jl") + @time @safetestset "Kantorovich homotopy controller" include("Core/kantorovich_homotopy_tests__item1.jl") @time @safetestset "ArcLengthContinuation construction + defaults" include("Core/arclength_tests__item1.jl") @time @safetestset "ArcLengthContinuation happy path (fold-free, matches sweep)" include("Core/arclength_tests__item2.jl") @time @safetestset "ArcLengthContinuation rounds a fold (non-monotone λ)" include("Core/arclength_tests__item3.jl")