SimpleNonlinearSolve.jl

SimpleNonlinearSolve

Small, dependency-light nonlinear solvers.

This subpackage provides straightforward nonlinear solver algorithms such as SimpleNewtonRaphson, SimpleBroyden, SimpleTrustRegion, and SimpleDFSane. They are intended for direct use on small problems, scalar problems, and allocation sensitive contexts where the full NonlinearSolve.jl stack is unnecessary.

Example

using SimpleNonlinearSolve, SciMLBase

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

SimpleNewtonRaphson(autodiff)
SimpleNewtonRaphson(; autodiff = nothing)

A low-overhead implementation of Newton-Raphson. This method is non-allocating on scalar and static array problems.

Keyword Arguments

• autodiff: determines the backend used for the Jacobian. Defaults to nothing (i.e. automatic backend selection). Valid choices include jacobian backends from DifferentiationInterface.jl.

SimpleGaussNewton(autodiff)
SimpleGaussNewton(; autodiff = nothing)

Alias for SimpleNewtonRaphson used for nonlinear least-squares problems. It uses the same low-overhead Newton implementation and Jacobian backend selection.

SimpleBroyden(; linesearch = nothing, alpha = nothing)

A low-overhead implementation of Broyden. This method is non-allocating on scalar and static array problems.

Keyword Arguments

• linesearch: nothing for no line search, or any LineSearch.AbstractLineSearchAlgorithm. Extra keyword arguments to solve are forwarded to LineSearch.init. For more options, use Broyden from NonlinearSolve.jl.
• alpha: Scale the initial jacobian initialization with alpha. If it is nothing, we will compute the scaling using 2 * norm(fu) / max(norm(u), true).

SimpleHalley(autodiff)
SimpleHalley(; autodiff = nothing)

A low-overhead implementation of Halley's Method.

Keyword Arguments

• autodiff: determines the backend used for the Jacobian. Defaults to nothing (i.e. automatic backend selection). Valid choices include jacobian backends from DifferentiationInterface.jl.

SimpleKlement()

A low-overhead implementation of Klement [4]. This method is non-allocating on scalar and static array problems.

SimpleTrustRegion(;
    autodiff = AutoForwardDiff(), max_trust_radius = 0.0,
    initial_trust_radius = 0.0, step_threshold = nothing,
    shrink_threshold = nothing, expand_threshold = nothing,
    shrink_factor = 0.25, expand_factor = 2.0, max_shrink_times::Int = 32,
    nlsolve_update_rule = Val(false)
)

A low-overhead implementation of a trust-region solver. This method is non-allocating on scalar and static array problems.

Keyword Arguments

• autodiff: determines the backend used for the Jacobian. Defaults to nothing (i.e. automatic backend selection). Valid choices include jacobian backends from DifferentiationInterface.jl.
• max_trust_radius: the maximum radius of the trust region. Defaults to max(norm(f(u0)), maximum(u0) - minimum(u0)).
• initial_trust_radius: the initial trust region radius. Defaults to max_trust_radius / 11.
• step_threshold: the threshold for taking a step. In every iteration, the threshold is compared with a value r, which is the actual reduction in the objective function divided by the predicted reduction. If step_threshold > r the model is not a good approximation, and the step is rejected. Defaults to 0.0001. For more details, see Rahpeymaii, F.
• shrink_threshold: the threshold for shrinking the trust region radius. In every iteration, the threshold is compared with a value r which is the actual reduction in the objective function divided by the predicted reduction. If shrink_threshold > r the trust region radius is shrunk by shrink_factor. Defaults to 0.25. For more details, see Rahpeymaii, F.
• expand_threshold: the threshold for expanding the trust region radius. If a step is taken, i.e step_threshold < r (with r defined in shrink_threshold), a check is also made to see if expand_threshold < r. If that is true, the trust region radius is expanded by expand_factor. Defaults to 0.75.
• shrink_factor: the factor to shrink the trust region radius with if shrink_threshold > r (with r defined in shrink_threshold). Defaults to 0.25.
• expand_factor: the factor to expand the trust region radius with if expand_threshold < r (with r defined in shrink_threshold). Defaults to 2.0.
• max_shrink_times: the maximum number of times to shrink the trust region radius in a row, max_shrink_times is exceeded, the algorithm returns. Defaults to 32.
• nlsolve_update_rule: If set to Val(true), updates the trust region radius using the update rule from NLSolve.jl. Defaults to Val(false). If set to Val(true), few of the radius update parameters – step_threshold = 0.05, expand_threshold = 0.9, and shrink_factor = 0.5 – have different defaults.

SimpleDFSane(;
    σ_min::Real = 1.0e-10, σ_max::Real = 1.0e10, σ_1::Real = 1.0,
    M::Union{Int, Val} = Val(10), γ::Real = 1.0e-4, τ_min::Real = 0.1, τ_max::Real = 0.5,
    nexp::Int = 2, η_strategy::Function = (f_1, k, x, F) -> f_1 ./ k^2
)

A low-overhead 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

• σ_min: the minimum value of the spectral coefficient σ_k which is related to the step size in the algorithm. Defaults to 1e-10.
• σ_max: the maximum value of the spectral coefficient σ_k which is related to the step size in the algorithm. Defaults to 1e10.
• σ_1: the initial value of the spectral coefficient σ_k which is related to the step size in the algorithm.. Defaults to 1.0.
• M: The monotonicity of the algorithm is determined by a this positive integer. A value of 1 for M would result in strict monotonicity in the decrease of the L2-norm of the function f. However, higher values allow for more flexibility in this reduction. Despite this, the algorithm still ensures global convergence through the use of a non-monotone line-search algorithm that adheres to the Grippo-Lampariello-Lucidi condition. Values in the range of 5 to 20 are usually sufficient, but some cases may call for a higher value of M. The default setting is 10.
• γ: a parameter that influences if a proposed step will be accepted. Higher value of γ will make the algorithm more restrictive in accepting steps. Defaults to 1e-4.
• τ_min: if a step is rejected the new step size will get multiplied by factor, and this parameter is the minimum value of that factor. Defaults to 0.1.
• τ_max: if a step is rejected the new step size will get multiplied by factor, and this parameter is the maximum value of that factor. Defaults to 0.5.
• nexp: the exponent of the loss, i.e. $f_k=||F(x_k)||^{nexp}$. The paper uses nexp ∈ {1,2}. Defaults to 2.
• η_strategy: function to determine the parameter η_k, which enables growth of $||F||^2$. Called as η_k = η_strategy(f_1, k, x, F) with f_1 initialized as $f_1=||F(x_1)||^{nexp}$, k is the iteration number, x is the current x-value and F the current residual. Should satisfy $η_k > 0$ and $∑ₖ ηₖ < ∞$. Defaults to $||F||^2 / k^2$.

SimpleLimitedMemoryBroyden(;
    threshold::Union{Val, Int} = Val(27), linesearch = nothing, alpha = nothing
)

A limited memory implementation of Broyden. This method applies the L-BFGS scheme to Broyden's method.

If the threshold is larger than the problem size, then this method will use SimpleBroyden.

Keyword Arguments:

• linesearch: nothing for no line search, or any LineSearch.AbstractLineSearchAlgorithm. Extra keyword arguments to solve are forwarded to LineSearch.init. For more options, use Broyden from NonlinearSolve.jl.
• alpha: Scale the initial jacobian initialization with alpha. If it is nothing, we will compute the scaling using 2 * norm(fu) / max(norm(u), true).

SimpleHomotopySweep(;
    inner = SimpleNewtonRaphson(), 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
)

Natural-parameter continuation solver for a SciMLBase.HomotopyProblem, the SimpleNonlinearSolve counterpart of HomotopySweep. The algorithm is the same — anchor solve at λspan[1], then predictor-corrector λ-stepping with the classic success/failure step control (failure halves the increment; expand_threshold consecutive successes grow it by expand_factor up to max_step_factor of the span, gated on the expand_quality secant-prediction error estimate), with a trust-monitored :secant warm-start predictor (:constant disables extrapolation) — but the driver is written in the direct, value-oriented SimpleNonlinearSolve style: each step calls solve on a freshly constructed (stack-allocated) inner problem instead of maintaining an inner-solver cache, and all sweep state is plain values.

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, and jac_prototype, sparsity, and colorvec are forwarded unchanged. The prototype is not eltype-promoted with λ: supply a prototype whose eltype matches the promoted residual eltype if λ's precision differs from u0's.

With a StaticArray (or scalar) u0 and a SimpleNonlinearSolve inner solver, the entire sweep is non-allocating, which also makes it the variant of choice inside hot loops, on GPUs, and for compilation-sensitive targets. For large mutable systems, prefer HomotopySweep, whose cached driver reuses the inner solver's workspace across steps.

Keyword arguments are identical to HomotopySweep except that inner defaults to SimpleNewtonRaphson() (a SimpleNonlinearSolve corrector) rather than the NonlinearSolve polyalgorithm. In particular tracking_maxiters caps the inner corrector's iterations for interior tracking steps only (never the λspan[1] anchor, never the final landing on λspan[2], and never overriding an explicit user-passed maxiters), tracking_abstol (default nothing = disabled) loosens the inner corrector's absolute tolerance on those same interior steps only — the anchor and the final landing always run at the user's full tolerances (each step here is a fresh standalone solve, so the exemption is simply not splicing the loose tolerance into the exempt solves — no re-polish is needed), and an explicit user-passed abstol or reltol (solve kwarg or problem kwarg) disables the loosening entirely — the success-side step growth is scaled by the corrector's iteration count (the AUTO-07p ADPTDS effort bands, including a proactive halving when a success nearly exhausts the budget), and maxsteps caps the total number of predictor-corrector attempts, returning ReturnCode.MaxIters with the last converged iterate when exceeded.

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 and original come from the most recent inner solve (unlike HomotopySweep, the ReturnCode.Stalled path reports the residual of the last successful step rather than nothing: every return path builds the same concrete solution type, which is part of what keeps the sweep allocation-free).