NonlinearSolve.jl Solvers

NewtonRaphson(;
    concrete_jac = nothing, linsolve = nothing, linesearch = missing,
    autodiff = nothing, vjp_autodiff = nothing, jvp_autodiff = nothing,
    forcing = nothing, jacobian_reuse = nothing,
)

An advanced NewtonRaphson implementation with support for efficient handling of sparse matrices via colored automatic differentiation and preconditioned linear solvers. Designed for large-scale and numerically-difficult nonlinear systems.

Keyword Arguments

• autodiff: determines the backend used for the Jacobian. Defaults to nothing which means that a default is selected according to the problem specification.
• concrete_jac: whether to build a concrete Jacobian. If a Krylov-subspace method is used, then the Jacobian will not be constructed. Defaults to nothing.
• linsolve: the LinearSolve.jl solver used for the linear solves within the Newton method. Defaults to nothing, which means it uses the LinearSolve.jl default algorithm choice.
• linesearch: the line search algorithm to use. Defaults to missing (no line search).
• vjp_autodiff: backend for computing Vector-Jacobian products.
• jvp_autodiff: backend for computing Jacobian-Vector products.
• forcing: Adaptive forcing term strategy for Newton-Krylov methods. When using an iterative linear solver (e.g., KrylovJL_GMRES()), this controls how accurately the linear system is solved at each iteration. Use EisenstatWalkerForcing2() for the classical Eisenstat-Walker adaptive forcing strategy. Defaults to nothing (fixed tolerance from the termination condition).
• jacobian_reuse: a JacobianReuse policy, true to force the default policy on, or false to force it off. Defaults to nothing, which reuses the Jacobian when length(u0) ≥ 16.

DFSane(;
    sigma_min = 1 // 10^10, sigma_max = 1.0e10, sigma_1 = 1, M::Int = 10,
    gamma = 1 // 10^4, tau_min = 1 // 10, tau_max = 1 // 2, n_exp::Int = 2,
    max_inner_iterations::Int = 100, eta_strategy = (fn_1, n, x_n, f_n) -> fn_1 / n^2
)

A low-overhead and allocation-free implementation of the df-sane method for solving large-scale nonlinear systems of equations. For in depth information about all the parameters and the algorithm, see La Cruz et al. [2].

Keyword Arguments

• sigma_min: the minimum value of the spectral coefficient σ which is related to the step size in the algorithm. Defaults to 1e-10.
• sigma_max: the maximum value of the spectral coefficient σₙ which is related to the step size in the algorithm. Defaults to 1e10.

For other keyword arguments, see RobustNonMonotoneLineSearch in LineSearch.jl.

Broyden(;
    max_resets::Int = 100, linesearch = nothing, reset_tolerance = nothing,
    init_jacobian::Val = Val(:identity), autodiff = nothing, alpha = nothing,
    update_rule = Val(:good_broyden)
)

An implementation of Broyden's Method [3] with resetting and line search.

Keyword Arguments

• max_resets: the maximum number of resets to perform. Defaults to 100.

• reset_tolerance: the tolerance for the reset check. Defaults to sqrt(eps(real(eltype(u)))).

• alpha: If init_jacobian is set to Val(:identity), then the initial Jacobian inverse is set to be (αI)⁻¹. Defaults to nothing which implies α = max(norm(u), 1) / (2 * norm(fu)).

• init_jacobian: the method to use for initializing the jacobian. Defaults to Val(:identity). Choices include:

  • Val(:identity): Identity Matrix.
  • Val(:true_jacobian): True Jacobian. This is a good choice for differentiable problems.

• update_rule: Update Rule for the Jacobian. Choices are:

  • Val(:good_broyden): Good Broyden's Update Rule
  • Val(:bad_broyden): Bad Broyden's Update Rule
  • Val(:diagonal): Only update the diagonal of the Jacobian. This algorithm may be useful for specific problems, but whether it will work may depend strongly on the problem

Klement(;
    max_resets = 100, linsolve = nothing, linesearch = nothing,
    alpha = nothing, init_jacobian::Val = Val(:identity),
    autodiff = nothing
)

An implementation of Klement [4] with line search, preconditioning and customizable linear solves. It is recommended to use Broyden for most problems over this.

Keyword Arguments

• max_resets: the maximum number of resets to perform. Defaults to 100.

• alpha: If init_jacobian is set to Val(:identity), then the initial Jacobian inverse is set to be αI. Defaults to 1. Can be set to nothing which implies α = max(norm(u), 1) / (2 * norm(fu)).

• init_jacobian: the method to use for initializing the jacobian. Defaults to Val(:identity). Choices include:

  • Val(:identity): Identity Matrix.
  • Val(:true_jacobian): True Jacobian. Our tests suggest that this is not very stable. Instead using Broyden with Val(:true_jacobian) gives faster and more reliable convergence.
  • Val(:true_jacobian_diagonal): Diagonal of True Jacobian. This is a good choice for differentiable problems.

LimitedMemoryBroyden(;
    max_resets::Int = 3, linesearch = nothing, threshold::Val = Val(10),
    reset_tolerance = nothing, alpha = nothing
)

An implementation of LimitedMemoryBroyden [5] with resetting and line search.

Keyword Arguments

• max_resets: the maximum number of resets to perform. Defaults to 3.
• reset_tolerance: the tolerance for the reset check. Defaults to sqrt(eps(real(eltype(u)))).
• threshold: the number of vectors to store in the low rank approximation. Defaults to Val(10).
• alpha: The initial Jacobian inverse is set to be (αI)⁻¹. Defaults to nothing which implies α = max(norm(u), 1) / (2 * norm(fu)).

HomotopySweep(;
    inner = nothing, nsteps = nothing, adaptive = true,
    initial_step_factor = 0.1, min_dλ = nothing, max_step_factor = 1.0,
    expand_factor = 2.0, expand_threshold = 2, expand_quality = 0.25,
    predictor = :secant, tracking_maxiters = 10, tracking_abstol = nothing,
    maxsteps = 10000, store_original = Val(false)
)

Natural-parameter continuation solver for SciMLBase.HomotopyProblem. The scalar continuation parameter $λ$ is swept across the problem's λspan. The sweep first solves the system at λspan[1] (for the canonical (0, 1) span, the simplified system — the form the homotopy is designed to make solvable from a cold start) starting from u0; each subsequent step fixes $λ$, predicts a warm start by extrapolating along the solution path, and corrects it by solving the resulting standard nonlinear system with inner.

The inner solver is initialized once and re-driven each step through the init/reinit!/solve! cache interface, so the continuation loop reuses the inner solver's workspace (Jacobian buffers, linear-solver storage) instead of reconstructing the solver every step; the sweep's own per-step state lives in a fixed set of preallocated buffers. The complete sweep also supports the standard init/reinit!/solve! interface: initialize once, optionally update u0 and p with reinit!, and call solve! repeatedly to reuse both the continuation buffers and the inner nonlinear-solver cache across homotopy problems of the same type. Solver options such as abstol can also be updated when the option was supplied to init and its type is unchanged.

When the inner solver is a polyalgorithm (as the default is), the warm-started tracking steps additionally arm best-subalgorithm retention on that cache (see NonlinearSolvePolyAlgorithm): each step resumes from the subalgorithm that produced the previous step's success instead of re-running the polyalgorithm ladder from its start, escalating (and eventually wrapping around) only when the retained subalgorithm fails. The cold anchor solve at λspan[1] always runs the full ladder — that is where the winning subalgorithm is discovered.

Optional derivative fields of the problem's NonlinearFunction (which SciMLBase.HomotopyProblem requires to follow the same λ-extended argument convention as the residual) are consumed by this solver: an analytic jac(u, p, λ) / jac(J, u, p, λ) is λ-fixed exactly like the residual and handed to the inner solver as a standard 2/3-argument Jacobian (so e.g. the default polyalgorithm selects its Jacobian-based members), and jac_prototype, sparsity, and colorvec are forwarded unchanged, enabling sparse Jacobian handling in every inner solve. The prototype is not eltype-promoted with λ: supply a prototype whose eltype matches the promoted residual eltype if λ's precision differs from u0's.

The step size is governed by the classic success/failure heuristic of predictor-corrector path tracking (see e.g. Timme, Mixed precision path tracking for polynomial homotopy continuation, Advances in Computational Mathematics 47, 2021): a failed corrector halves the λ increment and retries from the last accepted point, while expand_threshold consecutive accepted steps grow the increment by expand_factor, capped at max_step_factor of the span width. Expansion is additionally gated on the quality of the secant prediction (a Deuflhard-style local error estimate): the step only grows when the corrector's correction was small relative to how far the solution moved, so the increment does not balloon right before a sharp turn in the path — where an oversized step would be rejected only after the inner solver exhausts its iterations. On success the growth multiplier is additionally scaled by the corrector's iteration count (the AUTO-07p ADPTDS bands): a near-free corrector earns the full expand_factor, moderate effort earns milder growth, effort past a quarter of the iteration budget holds the increment, and a success that nearly exhausted the budget proactively halves it — the corrector working that hard on an accepted step is the earliest warning that the next full-size step will be rejected. This lets the sweep crawl through ill-conditioned regions of the path and accelerate back out of them while keeping trial-and-error rejections cheap.

Keyword arguments:

• inner: the inner nonlinear algorithm; nothing selects NonlinearSolve's default polyalgorithm (NOT a hardcoded Newton).
• nsteps: when given, the initial λ increment is the span width divided by nsteps instead of initial_step_factor. Required when adaptive = false (the steps are then fixed-size).
• adaptive: when true (default), a step whose inner solve fails to converge halves the λ increment and retries, down to a floor of min_dλ, and consecutive successes expand the increment as described above.
• initial_step_factor: the initial λ increment as a fraction of the λspan width; used when nsteps is not given.
• min_dλ: the smallest λ increment bisection may reach; nothing (default) resolves to sqrt(eps(typeof(λ))) at solve time, so the floor scales with precision.
• max_step_factor: the largest λ increment, as a fraction of the λspan width, that success expansion may reach. Must be in (0, 1]. Smaller values bound how far any single step can move along the path, which reduces the risk of the corrector converging to a different solution branch ("path jumping") on multi-branch problems.
• expand_factor: the multiplier applied to the λ increment after expand_threshold consecutive successful steps. Must be ≥ 1; 1 disables expansion.
• expand_threshold: the number of consecutive successful steps required before the increment is expanded. Must be ≥ 1. Larger values make regrowth more cautious after a bisection, avoiding repeated fail-shrink-regrow churn inside a hard region.
• expand_quality: expansion additionally requires the secant prediction's error ‖u - u_predicted‖ to be at most expand_quality times the scale of the recent solution movement. The error is measured against the prediction the secant would have made regardless of the predictor setting, so the gate is active for both predictors. A step whose corrector reports convergence within 2 iterations passes the gate outright — the warm start was deep inside the convergence basin, which is the strongest evidence that a larger step is affordable (and the relative error measure is uninformative on stretches where the path barely moves). Must be positive; Inf disables the gate, leaving the unconditional success/failure heuristic.
• predictor: :secant (default) extrapolates the initial guess for the next step linearly through the last two accepted points, so the corrector starts on the path tangent rather than at the previous solution; :constant warm-starts from the previous solution unchanged. The secant is trust-monitored: whenever its measured prediction error is no better than half that of the trivial constant prediction (as happens right after a sharp turn in the path, where the stale tangent points away from the path), or a step is rejected outright, subsequent steps fall back to the constant warm start until two consecutive accepted steps measure good secant quality again.
• tracking_maxiters: iteration cap for the inner solver on interior tracking steps (default 10, in the range used by MatCont, HomotopyContinuation.jl, and OpenModelica; nothing disables). A rejected step retries at half the increment from a warm start, so failing fast is far cheaper than exhausting the inner solver's full budget. Never applied to the λspan[1] anchor solve or the final step landing on λspan[2]; an explicit user-passed maxiters always wins.
• tracking_abstol: loose absolute tolerance for the inner corrector on interior tracking steps (default nothing = disabled: every step solves to the full tolerance). Interior iterates only serve as warm starts and secant-predictor history for the next step, so they need just enough accuracy to stay inside the next corrector's convergence basin — the reference trackers exploit exactly this (Bertini tracks at 1e-5/1e-6 and only the endgame runs tight; HomotopyContinuation.jl accepts on Newton contraction and polishes the endpoint). Values around 1e-4–1e-6 are good starting points when opting in. Never applied to the λspan[1] anchor solve (the one cold start, and the returned solution for a zero-width span) or to the final step landing on λspan[2], so the returned solution always satisfies the full tolerances: the landing runs on the loose cache first and is then re-polished at the full tolerance from that warm start (~1–2 extra corrector iterations). An explicit user-passed abstol or reltol (solve kwarg or problem kwarg) disables the loosening entirely. The looser interior iterates slightly degrade the secant/quality signals the step controller reads, which is why the default stays tight (opt-in loose, per the discussion in SciML/NonlinearSolve.jl#1020).
• maxsteps: a hard cap on the total number of predictor-corrector attempts (accepted steps plus bisection retries). Exceeding it returns a ReturnCode.MaxIters failure carrying the last converged iterate. Must be ≥ 1.
• store_original: whether to store the failing inner solve in the original field of the returned solution. Default Val(false) to keep the returned solution type concrete: the inner solve is a fresh solve(inner_prob, inner) whose return type inference gives up (its residual is a FixLambda wrapper), so storing it would pin the driver's returned solution's original type-slot to Any, forcing dynamic dispatch on every downstream field read of that solution. Set to Val(true) to keep the payload for debugging (the returned type is then no longer concrete). Mirrors NonlinearSolvePolyAlgorithm's option of the same name.

When the sweep cannot reach the end of λspan, the returned solution carries a failure retcode: its u is the last converged iterate (at some $λ$ short of λspan[2], or u0 itself if the initial λspan[1] anchor solve failed), while resid comes from the most recent inner solve. On the ReturnCode.Stalled and ReturnCode.MaxIters paths no new failed corrector is available, so the residual is from the last accepted point.

This is the embedding-homotopy / continuation analogue used to robustly initialize systems whose target form is hard to solve cold; it is unrelated to the polynomial HomotopyContinuationJL.

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 SciMLBase.HomotopyProblem, 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

\[\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

\[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. 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.
• 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 options.

This algorithm follows λ monotonically and therefore cannot round a fold. Combine it with ArcLengthContinuation in a HomotopyPolyAlgorithm when a fold must be traversed.

ArcLengthContinuation(;
    inner = nothing, initial_step_factor = 0.1,
    adaptive = true, min_ds = nothing, max_step_factor = 1.0,
    expand_factor = 2.0, expand_threshold = 2, max_angle = π / 6,
    predictor = :secant, autodiff = nothing, linsolve = nothing,
    tracking_maxiters = 10, maxsteps = 10000, theta = 0.5,
    store_original = Val(false)
)

Pseudo-arclength continuation solver for a SciMLBase.HomotopyProblem. Unlike HomotopySweep, which marches the scalar parameter $λ$ monotonically, this solver tracks the solution curve $H(u, λ) = 0$ parameterized by arclength $s$ in the augmented $(u, λ)$ space. Each step takes a predictor step along the path and corrects it by solving the augmented$(n+1)$-dimensional system

H(u, λ)                 = 0          # n equations
τ ⋅ ([u; λ] - x₀) - Δs  = 0          # Keller pseudo-arclength constraint

with the inner solver. Because $λ$ is a free variable of the corrector (not held fixed), the augmented Jacobian stays nonsingular at turning points (folds) where $∂H/∂u$ is singular and $λ$ is non-monotone along the path. This lets the solver round folds that defeat natural-parameter continuation — the canonical reason a sweep "fails to reach λ = 1" when a real solution at λ = 1 does exist but only on a branch reachable by going around a fold.

The target is the point on the curve where $λ = λspan[2]$: the solver follows the path until a step brackets that $λ$, then performs one final $λ$-fixed correction to land on it exactly.

Optional derivative fields of the problem's NonlinearFunction (which follow the same λ-extended argument convention as the residual) are consumed. An analytic jac(u, p, λ) / jac(J, u, p, λ) supplies the $∂H/∂u$ block of the augmented path Jacobian $[∂H/∂u | ∂H/∂λ]$; the missing $∂H/∂λ$ column is a scalar-parameter derivative, obtained as one forward-mode derivative of the residual in $λ$ at fixed u (through autodiff) — the packed system is never differentiated wholesale. The assembled path Jacobian drives the :tangent predictor and, extended by the analytically known θ-weighted Keller constraint row, gives every augmented corrector solve a full analytic $(n+1)×(n+1)$ Jacobian; the λ-fixed anchor and landing solves consume the jac exactly as HomotopySweep does. A jac_prototype (or matrix sparsity) is extended to the augmented shapes: one structurally dense $∂H/∂λ$ column for the predictor's $n×(n+1)$ system, plus the structurally dense constraint row for the corrector's bordered $(n+1)×(n+1)$ system. Sparse and structured prototypes are promoted to SparseMatrixCSC — a bordered Tridiagonal is no longer tridiagonal, and CSC is the general container the coloring and sparse linear-solve machinery handle; this requires SparseArrays to be loaded (structured prototypes fall back to a dense bordered prototype otherwise), and the prototype is not eltype-promoted with λ. A sparsity detector is forwarded unchanged (it detects the augmented pattern from the augmented residual itself). A user colorvec is forwarded to the predictor's system extended by one fresh color for the dense λ column; it is not forwarded to the corrector, whose dense constraint row admits no nontrivial column coloring — there the bordered prototype's value is chiefly the sparse linear solve. When no derivative fields are present, construction is identical to before.

⟨(u₁, λ₁), (u₂, λ₂)⟩_θ = (θ/n)⋅⟨u₁, u₂⟩ + (1 - θ)⋅λ₁λ₂,       n = length(u)

Keyword arguments:

• inner: the inner nonlinear algorithm used for both the initial on-curve correction and the augmented corrector; nothing selects NonlinearSolve's default polyalgorithm. A polyalgorithm inner runs the augmented corrector with best-subalgorithm retention (see NonlinearSolvePolyAlgorithm): after the first corrector solve discovers the winning subalgorithm, each warm-started corrector resumes from it instead of re-running the ladder from its start, escalating only when it fails. The λ-fixed anchor and landing solves always run the full ladder.
• initial_step_factor: the initial arclength step Δs as a fraction of the λspan width. Must be in (0, 1].
• adaptive: when true (default), a corrector failure halves Δs and retries from the last accepted point (down to a floor of min_ds), and expand_threshold consecutive successes grow Δs by expand_factor up to max_step_factor of the span.
• min_ds: the smallest arclength step bisection may reach; nothing (default) resolves to sqrt(eps(typeof(λ))).
• max_step_factor: the largest arclength step, as a fraction of the λspan width. Must be in (0, 1].
• expand_factor: the Δs growth multiplier after expand_threshold consecutive successful steps. Must be ≥ 1; 1 disables expansion.
• expand_threshold: consecutive successful steps required before Δs is expanded. Must be ≥ 1.
• max_angle: the curvature control (radians, in (0, π]). A step is rejected and Δs halved when the path direction turns by more than max_angle between the previous and current accepted segments; Δs is only allowed to grow when the turn is below max_angle / 3. Because the solution curve is smooth in arclength even at a fold (the tangent rotates continuously), bounding the per-step turn forces small steps through a turning point while permitting large steps on straight stretches — and it is what prevents the secant predictor from overshooting onto a different branch ("path jumping"). This is the analogue of OpenModelica's homotopy bend parameter.
• predictor: how the initial guess for each corrector is extrapolated along the path. :secant (default) uses the secant through the last two accepted points, bootstrapped from a pure-$λ$ step — derivative-free, but the bootstrap step cannot round a fold located within the very first step, and it is not curvature-checked. :tangent instead computes the true path tangent at the current point as the (oriented) null vector of the augmented Jacobian $[∂H/∂u | ∂H/∂λ]$, which stays well-defined at a fold (where the tangent is vertical in $λ$). It is obtained from the bordered linear solve $[J; τ_prevᵀ] t = e_{n+1}$ (one LU factorization per step, sparse Jacobians supported), falling back to a dense SVD null-space computation only when the bordered matrix is (near-)singular — e.g. exactly at a branch point, or when the previous tangent is orthogonal to the path. The tangent is a higher-order predictor and is accurate from the first step, so it curvature-checks every step and can round a fold at the very start; the cost is one Jacobian factorization per step (see autodiff). This is the Euler tangent predictor of the classic path trackers and of OpenModelica's global homotopy.
• autodiff: the automatic-differentiation backend (an ADTypes.AbstractADType) used to form the augmented Jacobian for the :tangent predictor and, when the problem supplies an analytic jac, to take the single $∂H/∂λ$ scalar derivative that completes it; nothing (default) selects AutoForwardDiff(). Unused by the :secant predictor on problems without an analytic jac.
• linsolve: the LinearSolve.jl algorithm for the :tangent predictor's bordered solve (the same knob the Newton descent methods expose); nothing (default) selects LinearSolve's default. Unused by the :secant predictor.
• tracking_maxiters: iteration cap for the augmented corrector solves (default 10, in the range used by MatCont, HomotopyContinuation.jl, and OpenModelica; nothing disables). A rejected step retries at half the arclength increment from a warm start, so failing fast is far cheaper than exhausting the inner solver's full budget. Never applied to the anchor or final λ-fixed landing solves; an explicit user-passed maxiters always wins. On success, step growth is additionally scaled by the corrector's iteration count (AUTO-style bands) alongside the bend-angle gate.
• maxsteps: a hard cap on the total number of predictor-corrector attempts (including bisection retries). Required because the path is not monotone in λ, so a sweep that never reaches the target — a closed loop, or a branch escaping to infinity — would otherwise not terminate. Exceeding it returns a ReturnCode.MaxIters failure.
• theta: the weight $θ$ of the arclength metric (see the note above). Must be in (0, 1); the default 0.5 weighs the (size-normalized) state block and the parameter equally. Larger theta emphasizes the state components, smaller theta emphasizes $λ$.
• store_original: whether to store the failing inner solve in the original field of the returned solution. Default Val(false) to keep the returned solution type concrete (the anchor/corrector inner solves' return type inference gives up, so storing one would pin the returned solution's original slot to Any). Set to Val(true) to keep the payload for debugging. Mirrors NonlinearSolvePolyAlgorithm's option of the same name.

When the solver cannot reach λspan[2], the returned solution carries a failure retcode and its u is the last converged curve point.

With the default :secant predictor the continuation is derivative-free in the predictor; the augmented corrector obtains the derivatives it needs through the inner solver's own differentiation, exactly as a standard NonlinearProblem would. The :tangent predictor additionally differentiates the homotopy through autodiff to build the augmented Jacobian.

HomotopyPolyAlgorithm(algs::Tuple; warm_handoff = true, store_original = Val(false))
HomotopyPolyAlgorithm(; inner = nothing, warm_handoff = true, store_original = Val(false))

A polyalgorithm for SciMLBase.HomotopyProblem: a container for a tuple of continuation algorithms that are tried in order until one returns a solution with a successful retcode. The first success is returned immediately — later stages never run. If every stage fails, the last stage's failed solution is returned, so its retcode (and original, when the stage attaches one) describe the most robust attempt.

This is the default algorithm for SciMLBase.HomotopyProblem: solve(prob) and solve(prob, nothing) route here.

The zero-argument form defaults to

HomotopyPolyAlgorithm((HomotopySweep(), ArcLengthContinuation()))

which encodes the natural escalation for homotopy solves: HomotopySweep is the cheap first attempt — natural-parameter continuation marches the scalar $λ$ monotonically across λspan, reusing one inner-solver cache across all steps, but it can never reverse $λ$ and therefore cannot follow a solution branch around a fold (turning point). When the sweep fails, ArcLengthContinuation takes over: it tracks the curve by pseudo-arclength in the augmented $(u, λ)$ space, so $λ$ is free to decrease along the path and folds that defeat the sweep are rounded — at the higher cost of solving an $(n+1)$-dimensional corrector system per step.

Warm handoff

When a natural-parameter stage (HomotopySweep or KantorovichHomotopy) 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 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 fold is measurably cheaper. Should the warm-started attempt fail anyway, the stage is retried cold on the original full-range problem before the polyalgorithm moves on, so enabling the handoff never costs robustness relative to warm_handoff = false — only, in that rare double-failure case, the extra warm attempt.

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

• algs: a tuple of continuation algorithms to try in order. Each stage must support solve(prob::SciMLBase.HomotopyProblem, alg, args...; kwargs...).

Keyword Arguments

• inner: an inner nonlinear algorithm threaded into both default stages (HomotopySweep and ArcLengthContinuation) as their corrector, so a HomotopyProblem is continued with that specific algorithm — e.g. one carrying a chosen autodiff backend. nothing (default) leaves each stage its own default inner. Ignored when algs is passed explicitly.
• warm_handoff: when true (default), a stage following a partway-failed HomotopySweep or KantorovichHomotopy 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 as Any). Pass Val(true) to recover the stage solution through original for introspection, at the cost of the returned solution's type no longer being concrete.

Example

using NonlinearSolve

alg = HomotopyPolyAlgorithm() # HomotopySweep, then ArcLengthContinuation on failure
alg = HomotopyPolyAlgorithm(
    (
        HomotopySweep(; inner = NewtonRaphson()),
        ArcLengthContinuation(; predictor = :tangent),
    )
)
alg = HomotopyPolyAlgorithm(; warm_handoff = false) # always restart stages cold

FastShortcutHomotopyPolyalg(
    ::Type{T} = Float64;
    autodiff = nothing, concrete_jac = nothing, linsolve = nothing,
    vjp_autodiff = nothing, jvp_autodiff = nothing,
    warm_handoff::Bool = true, store_original::Val = Val(false)
) where {T}

The recommended default HomotopyPolyAlgorithm for solving a SciMLBase.HomotopyProblem — e.g. a Modelica homotopy(actual, simplified) initialization system — by continuation. It is the homotopy analogue of FastShortcutNonlinearPolyalg: a fast HomotopySweep (natural-parameter continuation) escalating to a robust ArcLengthContinuation (pseudo-arclength) on failure, with a FastShortcutNonlinearPolyalg built from the requested autodiff threaded in as the inner corrector of both stages.

Solving a HomotopyProblem with a plain nonlinear algorithm instead fixes $λ$ at the target and solves only the actual system, which can converge to the wrong branch. This sweeps $λ$ from the simplified anchor to the actual system, tracking the intended branch.

Arguments

• T: the eltype of the initial guess, forwarded to the inner FastShortcutNonlinearPolyalg. Defaults to Float64.

Keyword Arguments

• autodiff, concrete_jac, linsolve, vjp_autodiff, jvp_autodiff: forwarded to the inner FastShortcutNonlinearPolyalg that corrects each continuation step — this is where the differentiation backend is chosen. A HomotopyProblem whose residual is not ForwardDiff-safe is solved by passing autodiff = AutoFiniteDiff().
• warm_handoff, store_original: forwarded to HomotopyPolyAlgorithm.

GaussNewton(;
    concrete_jac = nothing, linsolve = nothing, linesearch = missing,
    autodiff = nothing, vjp_autodiff = nothing, jvp_autodiff = nothing,
    jacobian_reuse = nothing
)

An advanced GaussNewton implementation with support for efficient handling of sparse matrices via colored automatic differentiation and preconditioned linear solvers. Designed for large-scale and numerically-difficult nonlinear systems.

Keyword Arguments

• jacobian_reuse: a JacobianReuse policy, true to force the default policy on, or false to force it off. Defaults to nothing, which reuses the Jacobian when length(u0) ≥ 16.

TrustRegion(;
    concrete_jac = nothing, linsolve = nothing,
    radius_update_scheme = nothing, subproblem = TrustRegionSubproblem.More,
    max_trust_radius::Real = 0 // 1,
    initial_trust_radius::Real = 0 // 1, step_threshold::Real = 1 // 10000,
    shrink_threshold::Real = 1 // 4, expand_threshold::Real = 3 // 4,
    shrink_factor::Real = 1 // 4, expand_factor::Real = 2 // 1,
    max_shrink_times::Int = 32,
    vjp_autodiff = nothing, autodiff = nothing, jvp_autodiff = nothing,
    jacobian_reuse = nothing,
)

An advanced TrustRegion implementation with support for efficient handling of sparse matrices via colored automatic differentiation and preconditioned linear solvers. Designed for large-scale and numerically-difficult nonlinear systems.

Keyword Arguments

• radius_update_scheme: the scheme used to update the trust region radius. Defaults to RadiusUpdateSchemes.More, or RadiusUpdateSchemes.Simple when subproblem = TrustRegionSubproblem.Dogleg. See RadiusUpdateSchemes for more details. For a review on trust region radius update schemes, see Yuan [6].
• subproblem: how the trust-region subproblem is solved. TrustRegionSubproblem.More (default) uses MoreTrustRegionDescent, which solves the subproblem nearly exactly via Moré's safeguarded iteration on the damping parameter (MINPACK lmpar) and is substantially more robust on ill-conditioned least-squares problems; TrustRegionSubproblem.Dogleg uses Dogleg. A custom AbstractDescentDirection can also be passed directly — e.g. MoreTrustRegionDescent(; scaling = TrustRegionScaling.Auto) engages Moré's column-norm scaling only on badly scaled Jacobians.
• jacobian_reuse: a JacobianReuse policy, true to force the default policy on, or false to force it off. Defaults to nothing, which reuses the Jacobian when length(u0) ≥ 16. A rejected step computed from a fresh Jacobian reuses that Jacobian at the unchanged state.

For the remaining arguments, see NonlinearSolveFirstOrder.GenericTrustRegionScheme documentation.

TrustRegionDogleg(; subproblem = TrustRegionSubproblem.Dogleg, kwargs...)

TrustRegion with the classical two-piece polygonal dogleg subproblem — the cheaper per-step choice for well-conditioned problems where the nearly-exact Moré solve is not needed. TrustRegion() itself defaults to TrustRegionSubproblem.More, which is more robust on ill-conditioned and rank-deficient problems.

All TrustRegion keyword arguments are forwarded; subproblem is fixed to TrustRegionSubproblem.Dogleg (a Dogleg descent instance is also accepted) and radius_update_scheme defaults to RadiusUpdateSchemes.Simple as with TrustRegion.

TrustRegionRobust(;
    subproblem = MoreTrustRegionDescent(; linsolve = RobustTrustRegionLinsolve()),
    kwargs...
)

TrustRegion with the conditioning-preserving subproblem formulation: for sparse-structured Jacobians the damped solves run on the rectangular augmented system [J; √λD] p = [-fu; 0] via LinearSolve.SparseColumnPivotedQRFactorization — a rank-revealing column-pivoted sparse QR — rather than the normal equations JᵀJ + λD², which squares the condition number. Dense and matrix-free Jacobians keep their default paths, which already avoid the normal equations (MINPACK lmpar and the augmented Krylov operator respectively).

Use this when the Jacobian is sparse and badly conditioned enough that forming JᵀJ loses digits — at the price of a sparse QR refactorization of the (m + n) × n system per damping-parameter trial.

All TrustRegion keyword arguments are forwarded; subproblem is fixed to a MoreTrustRegionDescent carrying the robust linsolve selection.

LevenbergMarquardt(;
    linsolve = nothing,
    damping_initial::Real = 1.0, α_geodesic::Real = 0.75, disable_geodesic = Val(false),
    damping_increase_factor::Real = 2.0, damping_decrease_factor::Real = 3.0,
    finite_diff_step_geodesic = 0.1, b_uphill::Real = 1.0, min_damping_D::Real = 1e-8,
    autodiff = nothing, vjp_autodiff = nothing, jvp_autodiff = nothing,
    jacobian_reuse = nothing
)

An advanced Levenberg-Marquardt implementation with the improvements suggested in Transtrum and Sethna [1]. Designed for large-scale and numerically-difficult nonlinear systems.

Keyword Arguments

• damping_initial: the starting value for the damping factor. The damping factor is inversely proportional to the step size. The damping factor is adjusted during each iteration. Defaults to 1.0. See Section 2.1 of Transtrum and Sethna [1].
• damping_increase_factor: the factor by which the damping is increased if a step is rejected. Defaults to 2.0. See Section 2.1 of Transtrum and Sethna [1].
• damping_decrease_factor: the factor by which the damping is decreased if a step is accepted. Defaults to 3.0. See Section 2.1 of Transtrum and Sethna [1].
• min_damping_D: the minimum value of the damping terms in the diagonal damping matrix DᵀD, where DᵀD is given by the largest diagonal entries of JᵀJ yet encountered, where J is the Jacobian. It is suggested by Transtrum and Sethna [1] to use a minimum value of the elements in DᵀD to prevent the damping from being too small. Defaults to 1e-8.
• disable_geodesic: Disables Geodesic Acceleration if set to Val(true). It provides a way to trade-off robustness for speed, though in most situations Geodesic Acceleration should not be disabled.
• jacobian_reuse: a JacobianReuse policy, true to force the default policy on, or false to force it off. Defaults to nothing, which reuses the Jacobian when length(u0) ≥ 16.

For the remaining arguments, see GeodesicAcceleration and NonlinearSolveFirstOrder.LevenbergMarquardtTrustRegion documentations.

PseudoTransient(;
    concrete_jac = nothing, linesearch = missing, alpha_initial = 1e-3,
    linsolve = nothing, mass_matrix = nothing,
    autodiff = nothing, jvp_autodiff = nothing, vjp_autodiff = nothing,
    jacobian_reuse = nothing
)

An implementation of PseudoTransient Method [7] that is used to solve steady state problems in an accelerated manner. It uses an adaptive time-stepping to integrate an initial value of nonlinear problem until sufficient accuracy in the desired steady-state is achieved to switch over to Newton's method and gain a rapid convergence. This implementation specifically uses "switched evolution relaxation" [8] SER method.

The damped Newton step solves $(J(u) + (1/α) M) δu = -F(u)$, i.e. one implicit-Euler step of the fictitious dynamics $M u' = -F(u)$ with pseudo-timestep α. By default M = I, which recovers the classical pseudo-transient continuation. Supplying a mass matrix M generalizes this to differential-algebraic (DAE) steady-state problems [7], where M is the (possibly singular, structured) mass matrix of the underlying M u' = -F(u) system. This makes the continuation topology-aware: components with large M entries are damped less, algebraic components (zero rows of M) are treated consistently with the DAE structure.

Keyword Arguments

• alpha_initial : the initial pseudo time step. It defaults to 1e-3. If it is small, you are going to need more iterations to converge but it can be more stable.
• mass_matrix : the mass matrix M used for damping, i.e. the descent solves $(J + (1/α) M) δu = -F$. Defaults to nothing, which uses M = I (identity damping, bit-for-bit identical to the classical method). If nothing and the problem's NonlinearFunction carries a non-identity mass_matrix, that mass matrix is used automatically. A diagonal M (e.g. Diagonal(...)) uses an efficient diagonal update; a general sparse/dense M is supported as well. Intended for square DAE-derived systems.
• jacobian_reuse: a JacobianReuse policy, true to force the default policy on, or false to force it off. Defaults to nothing, which reuses the Jacobian when length(u0) ≥ 16. The damped system is still rebuilt when the pseudo-timestep changes.

NonlinearSolvePolyAlgorithm(algs; start_index::Int = 1, store_original = Val(false))

A general way to define PolyAlgorithms for NonlinearProblem and NonlinearLeastSquaresProblem. This is a container for a tuple of algorithms that will be tried in order until one succeeds. If none succeed, then the algorithm with the lowest residual is returned.

If every stage natively supports box bounds, the polyalgorithm preserves the original bounded coordinates and projects an infeasible initial guess onto the box without mutating the supplied state. Otherwise, bound handling uses the variable transformation shared by algorithms without native bound support.

Arguments

• algs: a tuple of algorithms to try in-order! (If this is not a Tuple, then the returned algorithm is not type-stable).

Keyword Arguments

• start_index: the index to start at. Defaults to 1.
• store_original: Whether to store the winning sub-algorithm's solution in the original field of the returned solution. Default Val(false) to keep the return type simple (required for Enzyme AD compatibility). Set to Val(true) for debugging to inspect the sub-algorithm's solution.

Example

using NonlinearSolve

alg = NonlinearSolvePolyAlgorithm((NewtonRaphson(), Broyden()))

Best-subalgorithm retention (reinit!(cache; retain_best = true))

When the polyalgorithm's cache is reused across a sequence of warm-started solves (init/reinit!/solve!), each reinit! by default restarts the ladder at start_index, so every solve re-fails the same cheap subalgorithms before reaching the one that works. Passing retain_best = true to reinit! makes the next solve! start from the subalgorithm that produced the most recent success instead. Escalation is preserved: if that subalgorithm fails, the ladder continues upward as usual, and once the last subalgorithm fails it wraps around to the subalgorithms that were skipped at the start (they occasionally succeed where the retained one stagnates), so retention never tries fewer subalgorithms than a full ladder run. When no success has been recorded yet, retain_best = true starts the ladder at start_index as usual. The default (retain_best = false) is the status-quo full restart.

The wrap-around never goes below the algorithm's start_index, so retention never attempts a subalgorithm the algorithm itself excludes.

Retention periodically re-probes the skipped cheaper subalgorithms: when the retained subalgorithm sits above start_index, every RETAIN_REPROBE_INTERVAL-th (8th) retained reinit! starts one solve from start_index again. A cheap subalgorithm that failed once transiently (escalating best to an expensive one that then always succeeds) is therefore rediscovered instead of being locked out for the rest of the warm-started sequence, at a bounded cost of at most one status-quo-style ladder step per interval.

A retain_best = truereinit! also reinitializes the subalgorithm caches lazily: only the starting subalgorithm's cache is reinitialized up front, and each further one is reinitialized at the moment escalation reaches it. Since every subcache reinit! evaluates the residual at the new u0, this reduces the per-reinit! residual cost from one evaluation per subalgorithm to one per subalgorithm actually attempted (usually just the retained one). A consequence is that the shared stats of a solution produced after mid-solve escalation only reflect the subalgorithms run since the last deferred reinitialization, i.e. effectively the winning subalgorithm's own effort. Updating solver options in the same reinit! eagerly updates every subcache so that later escalation observes the new values.

FastShortcutNonlinearPolyalg(
    ::Type{T} = Float64;
    concrete_jac = nothing,
    linsolve = nothing,
    must_use_jacobian::Val = Val(false),
    prefer_simplenonlinearsolve::Val = Val(false),
    autodiff = nothing, vjp_autodiff = nothing, jvp_autodiff = nothing,
    jacobian_reuse = nothing,
    u0_len::Union{Int, Nothing} = nothing
) where {T}

A polyalgorithm focused on balancing speed and robustness. It first tries less robust methods for more performance and then tries more robust techniques if the faster ones fail.

Arguments

• T: The eltype of the initial guess. It is only used to check if some of the algorithms are compatible with the problem type. Defaults to Float64.

Keyword Arguments

• u0_len: The length of the initial guess. If this is nothing, then the length of the initial guess is not checked. If this is an integer and it is less than 25, we use jacobian based methods.
• jacobian_reuse: forwarded to each first-order method in the polyalgorithm.

FastShortcutNLLSPolyalg(
    ::Type{T} = Float64;
    concrete_jac = nothing,
    linsolve = nothing,
    autodiff = nothing, vjp_autodiff = nothing, jvp_autodiff = nothing,
    jacobian_reuse = nothing
)

A polyalgorithm focused on balancing speed and robustness. It first tries less robust methods for more performance and then tries more robust techniques if the faster ones fail.

Arguments

• T: The eltype of the initial guess. It is only used to check if some of the algorithms are compatible with the problem type. Defaults to Float64.
• jacobian_reuse: forwarded to each first-order method in the polyalgorithm.

RobustMultiNewton(
    ::Type{T} = Float64;
    concrete_jac = nothing,
    linsolve = nothing,
    autodiff = nothing, vjp_autodiff = nothing, jvp_autodiff = nothing,
    jacobian_reuse = nothing
)

A polyalgorithm focused on robustness. It uses a mixture of Newton methods with different globalizing techniques (trust region updates, line searches, etc.) in order to find a method that is able to adequately solve the minimization problem.

Basically, if this algorithm fails, then "most" good ways of solving your problem fail and you may need to think about reformulating the model (either there is an issue with the model, or more precision / more stable linear solver choice is required).

Arguments

• T: The eltype of the initial guess. It is only used to check if some of the algorithms are compatible with the problem type. Defaults to Float64.
• jacobian_reuse: forwarded to each first-order method in the polyalgorithm.

NonlinearSolveFirstOrder

First-order nonlinear and nonlinear least-squares solver algorithms.

This subpackage implements Newton, Gauss-Newton, trust-region, pseudo-transient, and related algorithms that are re-exported by NonlinearSolve.jl. Users typically load NonlinearSolve and pass these algorithms to solve; solver-package authors may depend on this package directly when they need the first-order implementations.

Example

using NonlinearSolveFirstOrder, SciMLBase

prob = NonlinearProblem((u, p) -> u^2 - p, 1.0, 2.0)
sol = solve(prob, NewtonRaphson())

NonlinearSolveQuasiNewton

Quasi-Newton nonlinear solver algorithms.

This subpackage implements quasi-Newton methods such as Broyden, Klement, and LimitedMemoryBroyden. These algorithms are re-exported by NonlinearSolve.jl, but the subpackage can also be loaded directly by packages that only need quasi-Newton solver implementations.

Example

using NonlinearSolveQuasiNewton, SciMLBase

prob = NonlinearProblem((u, p) -> u^2 - p, 1.0, 2.0)
sol = solve(prob, Broyden())

NonlinearSolveSpectralMethods

Spectral residual methods for nonlinear systems.

This subpackage provides DFSane and the generalized spectral residual implementation used by NonlinearSolve.jl. Use it for derivative-free nonlinear solves where a spectral residual method is appropriate.

Example

using NonlinearSolveSpectralMethods, SciMLBase

prob = NonlinearProblem((u, p) -> u^2 - p, 1.0, 2.0)
sol = solve(prob, DFSane())

QuasiNewtonAlgorithm(;
    linesearch = missing, trustregion = missing, descent, update_rule, reinit_rule,
    initialization, max_resets::Int = typemax(Int), name::Symbol = :unknown,
    max_shrink_times::Int = typemax(Int), concrete_jac = Val(false)
)

Nonlinear Solve Algorithms using an Iterative Approximation of the Jacobian. Most common examples include Broyden's Method.

Keyword Arguments

• trustregion: Globalization using a Trust Region Method. This needs to follow the NonlinearSolveBase.AbstractTrustRegionMethod interface.
• descent: The descent method to use to compute the step. This needs to follow the NonlinearSolveBase.AbstractDescentDirection interface.
• max_shrink_times: The maximum number of times the trust region radius can be shrunk before the algorithm terminates.
• update_rule: The update rule to use to update the Jacobian. This needs to follow the NonlinearSolveBase.AbstractApproximateJacobianUpdateRule interface.
• reinit_rule: The reinitialization rule to use to reinitialize the Jacobian. This needs to follow the NonlinearSolveBase.AbstractResetCondition interface.
• initialization: The initialization method to use to initialize the Jacobian. This needs to follow the NonlinearSolveBase.AbstractJacobianInitialization interface.

GeneralizedFirstOrderAlgorithm(;
    descent, linesearch = missing,
    trustregion = missing, autodiff = nothing, vjp_autodiff = nothing,
    jvp_autodiff = nothing, max_shrink_times::Int = typemax(Int),
    concrete_jac = Val(false), jacobian_reuse = nothing, name::Symbol = :unknown
)

This is a Generalization of First-Order (uses Jacobian) Nonlinear Solve Algorithms. The most common example of this is Newton-Raphson Method.

First Order here refers to the order of differentiation, and should not be confused with the order of convergence.

Keyword Arguments

• trustregion: Globalization using a Trust Region Method. This needs to follow the NonlinearSolveBase.AbstractTrustRegionMethod interface.
• descent: The descent method to use to compute the step. This needs to follow the NonlinearSolveBase.AbstractDescentDirection interface.
• max_shrink_times: The maximum number of times the trust region radius can be shrunk before the algorithm terminates.
• jacobian_reuse: a JacobianReuse policy for reusing the Jacobian across accepted steps. true forces the default policy on and false forces it off. Defaults to nothing, which resolves against length(u0) when the cache is built.

GeneralizedDFSane(; linesearch, sigma_min, sigma_max, sigma_1, name::Symbol = :unknown)

A generalized version of the DF-SANE algorithm. This algorithm is a Jacobian-Free Spectral Method.

Arguments

• linesearch: Globalization using a Line Search Method. This is not optional currently, but that restriction might be lifted in the future.
• sigma_min: The minimum spectral parameter allowed. This is used to ensure that the spectral parameter is not too small.
• sigma_max: The maximum spectral parameter allowed. This is used to ensure that the spectral parameter is not too large.
• sigma_1: The initial spectral parameter. If this is not provided, then the algorithm initializes it as sigma_1 = <u, u> / <u, f(u)>.

JacobianReuse(; max_age::Int = 10,
    max_residual_ratio::Real = 0.1)

Reuse a Jacobian across accepted nonlinear iterations. This turns a first-order method into an adaptive modified-Newton method: the current Jacobian is reused while the residual norm keeps contracting fast enough, subject to a maximum Jacobian age. Solvers of an unchanged concrete linear system also reuse its factorization; damped and matrix-free systems retain their own linear-solver update behavior.

max_age is the number of accepted steps a single Jacobian may serve. The Jacobian is refreshed when any of these conditions holds:

• it has served max_age accepted steps;
• the new residual norm is not strictly less than max_residual_ratio times the previous residual norm;
• a linear solve or globalization step fails with stale Jacobian information.

max_age = 0 disables reuse and recovers exact Newton steps; it is how jacobian_reuse = false is spelled, and max_age = 1 means the same thing. Setting max_residual_ratio = Inf selects purely periodic refreshes. The reuse state is reset by reinit!; retaining a Jacobian across separate nonlinear solves requires the manual step!(cache; recompute_jacobian = false) interface.

Pass jacobian_reuse = JacobianReuse() to NewtonRaphson, TrustRegion, or another first-order solver to force the policy on, and jacobian_reuse = false to force it off. The default, jacobian_reuse = nothing, decides from the problem size: reuse is enabled when length(u0) ≥ 16 and disabled below it. A matrix-free Jacobian operator is bound to the current iterate on every step, so there is nothing to reuse and reuse is switched off for it.

EisenstatWalkerForcing2(; η₀ = 0.5, ηₘₐₓ = 0.9, γ = 0.9, α = 2, safeguard = true, safeguard_threshold = 0.1)

Algorithm 2 from the classical work by Eisenstat and Walker (1996) as described by formula (2.6): ηₖ = γ * (||rₖ|| / ||rₖ₋₁||)^α

Here the variables denote: rₖ residual at iteration k η₀ ∈ [0,1) initial value for η ηₘₐₓ ∈ [0,1) maximum value for η γ ∈ [0,1) correction factor α ∈ [1,2) correction exponent

Furthermore, the proposed safeguard is implemented: ηₖ = max(ηₖ, γηₖ₋₁^α) if γηₖ₋₁^α > safeguard_threshold to prevent ηₖ from shrinking too fast.