SimpleSolvers

const ARMIJO_τ_DEMAND_FRACTION

The largest fraction of the demanded decrease that the round-off resolution $\tau$ is allowed to amount to, used by armijo_ulps to cap $\tau$ in low precision. Its value is 0.01, i.e. $\tau$ may distort the SufficientDecreaseCondition by at most one percent.

const BACKTRACKING_SHRINK_MIN

Lower bound on the factor by which a rejected step is shrunk by the interpolation in Backtracking: the new trial step is confined to $[\mathrm{BACKTRACKING\_SHRINK\_MIN}\cdot\alpha, p\alpha]$ (see [1, §3.5], [2, Alg. A6.3.1]). Its value is 0.1.

const DEFAULT_ARMIJO_p

Constant used in Backtracking. Its value is 0.5

This is the default for the constant $p$ by which α is decreased if the SufficientDecreaseCondition and the CurvatureCondition are not satisfied.

const DEFAULT_ARMIJO_α₀

The default starting value for $\alpha$ used in Backtracking. Its value is 1.0.

const DEFAULT_ARMIJO_τ_ULPS

The nominal number of units in the last place (ulps) of $\varphi(0)$ taken as the round-off resolution $\tau$ of a merit function, i.e. $\tau = \mathrm{DEFAULT\_ARMIJO\_τ\_ULPS}\cdot\mathrm{ulp}(\varphi(0))$ (see armijo_tolerance). Its value is 4.

Use armijo_ulps rather than this constant: it caps the nominal value at what the element type can actually support, which matters in Float16.

$\tau$ is used for three things, and it is worth keeping them apart:

1. it slackens the SufficientDecreaseCondition inside Backtracking to $\varphi(\alpha) \leq \min\{\varphi(0),\ \varphi(0) + c_1\alpha\varphi'(0) + \tau\}$. The $\min$ is what keeps the allowance honest: it may reduce the decrease demanded, but it can never license a step whose merit exceeds $\varphi(0)$;
2. it fixes the smallest informative trial step $\alpha_\mathrm{min}$ below which the condition could only be decided by rounding (see backtracking_αmin);
3. every LinesearchMethod uses it to decide whether an accepted step's decrease was genuine (LINESEARCH_DECREASED) or within the noise (LINESEARCH_FLOOR, see LinesearchOutcome).

Set τ_ulps = 0 in Backtracking to recover the exact condition for (1) and (2).

const DEFAULT_BRACKETING_k

Gives the default ratio by which the bracket is increased if bracketing was not successful. See bracket_minimum.

Default constant. Number of maximum iterations for bracket_minimum, bracket_minimum_with_fixed_point and bracket_root.

const DEFAULT_BRACKETING_s

Gives the default initial width of the interval (the bracket). Used for bracket_minimum, bracket_minimum_with_fixed_point and bracket_root.

The default number of iterations before the Jacobian is refactored when constructing a quasi-Newton method via QuasiNewton.

Backtracking shrink factor used by the PicardSolver residual safeguard.

const DEFAULT_WOLFE_c₁

A constant $c_1$ that is used in the SufficientDecreaseCondition (the Armijo condition):

\[f(\alpha) \leq f(\alpha_0) + c_1 \alpha f'(\alpha_0).\]

const DEFAULT_WOLFE_c₂

The constant used in the second Wolfe condition (the CurvatureCondition). According to [1, 3] we should have

\[c_2 \in (c_1, 1),\]

where $c_1$ is the constant specified by DEFAULT_WOLFE_c₁.

Furthermore [1] recommend $c_2 = 0.9$; in [3] the authors write: "it is common to set $c_2=0.1$ when approximate line search is used with the conjugate gradient method and to 0.9 when used with Newton's method." We use $c_2 = 0.9$ as default.

const DEFAULT_WOLFE_αmax

Default upper bound on the step length for the bracketing phase of StrongWolfe. Its value is 65536.0.

A factor by which s is reduced in each bracketing iteration (see bracket_minimum_with_fixed_point).

Factor by which the trust-region radius is expanded on a very good step ($\rho > 3/4$ at the boundary); the default of the dogleg_radius_expand field of Options, which the solver actually reads.

Default initial trust-region radius for the DogLegSolver; the default of the dogleg_radius_initial field of Options, which the solver actually reads.

Default maximum trust-region radius ($\hat\Delta$ in [1, Alg. 4.1]) for the DogLegSolver; the radius is never expanded beyond this. The default of the dogleg_radius_max field of Options, which the solver actually reads.

Factor by which the trust-region radius is shrunk on a poor step ($\rho < 1/4$); the default of the dogleg_radius_shrink field of Options, which the solver actually reads.

Minimum ρ (actual/predicted reduction) for a step to be accepted ($\eta$ in [1, Alg. 4.1]).

Upper ρ threshold above which the trust-region radius may be expanded.

Lower ρ threshold below which the trust-region radius is shrunk.

const MAX_STALLS

The default number of consecutive stalled steps after which a NonlinearSolver gives up; the default of the max_stalls field of Options. Its value is 2.

A step is stalled when it does not move the iterate while the residual is not small (see stalled_step), i.e. when the merit $\|F\|^2$ cannot be reduced along the current direction. One stalled step is not conclusive, because the next step is guaranteed to be attempted under better conditions: a stall forces a fresh Jacobian immediately (see maybe_refactorize! and needs_refresh) rather than waiting for the next refactorize multiple, and the DogLegSolver additionally resets a collapsed trust-region radius. A second consecutive stall is therefore one that a freshly evaluated Jacobian did not fix, which is conclusive — and it is conclusive for every refactorize, not just refactorize = 1.

Set max_stalls = typemax(Int) to restore the previous behaviour of running all the way to max_iterations.

Threshold for the maximum size a static matrix should have. See _static.

Encompasses the NoLinearProblem and the LinearProblem. Subtyped from AbstractProblem, coming from GeometricBase.

AbstractNonlinearSolverCache

An abstract type that comprises e.g. the NonlinearSolverCache and the DogLegCache.

Backtracking <: LinesearchMethod

Keys

The keys are:

• α₀=1.0: the initial step size $\alpha$. This is decreased iteratively by a factor $p$ until the SufficientDecreaseCondition is satisfied.
• c₁=0.0001: the constant $c_1$ in the SufficientDecreaseCondition (Armijo condition). Also see DEFAULT_WOLFE_c₁.
• c₂=0.9: the constant on whose basis the CurvatureCondition is tested. We should have $c_2\in(c_1, 1).$ The closer this constant is to 1, the easier it is to satisfy the CurvatureCondition.
• p=0.5: an upper bound on the factor by which $\alpha$ is decreased in every step until the stopping criterion is satisfied. The actual factor is chosen by interpolation and confined to $[$ BACKTRACKING_SHRINK_MIN $\cdot\alpha, p\alpha]$, so the trial sequence is never longer than the plain $\alpha \gets p\alpha$ ladder.
• τ_ulps=armijo_ulps(T, c₁) (4 in Float64 and Float32, less in Float16): the round-off resolution of the merit, in units in the last place of $\varphi(0)$. It slackens the SufficientDecreaseCondition (never past $\varphi(0)$), fixes $\alpha_\mathrm{min}$, and separates a genuine decrease from one within the noise. A value larger than armijo_ulps(T, c₁) is capped to it, since above that $\tau$ would swamp the decrease the condition demands. See DEFAULT_ARMIJO_τ_ULPS.

Implementation

The algorithm starts by setting

\[\begin{aligned} \varphi_0 &\gets \varphi(0),\\ d_0 &\gets \varphi'(0), \end{aligned}\]

where $\varphi$ is of type LinesearchProblem. Unless $\varphi_0$ and $d_0$ are finite with $d_0 < 0$ the search is abandoned at once — no $\alpha$ can satisfy the SufficientDecreaseCondition along a direction that is not decreasing, so shrinking $\alpha$ would only waste merit evaluations to find that out.

Otherwise it sets the round-off resolution $\tau$ (armijo_tolerance) and the smallest informative step $\alpha_\mathrm{min}$ (backtracking_αmin), and shrinks the trial step by backtracking_interpolation until one of the following happens:

1. the SufficientDecreaseCondition is satisfied — the step is accepted, and reported as a genuine decrease only if $\varphi(\alpha) \leq \varphi_0 - \tau$;
2. two consecutive trials return $\varphi(\alpha) = \varphi_0$ bit-exactly — the trial point no longer differs from the base point in floating point, so no smaller step can either;
3. $\alpha \leq \alpha_\mathrm{min}$ — a smaller step could only be judged by rounding;
4. the linesearch_max_iterations budget of Options is spent.

Cases 2–4 are distinguished in the returned LinesearchStatus: a merit that does not vary by more than $\tau$ has reached its round-off floor (LINESEARCH_FLOOR, benign and not improvable by any line search), whereas one that does vary contradicts $d_0 < 0$ (LINESEARCH_EXHAUSTED, a genuine inconsistency). See LinesearchOutcome.

The CurvatureCondition is not used to terminate the iteration — it cannot be honoured by shrinking alone — it is only checked afterwards to emit a warning (see curvature_diagnostic).

Extended help

Sometimes the parameters $p$ and $c_1$ have different names such as $\tau$ and $c$. Note that our $\tau$ is something else entirely (the round-off resolution above).

$\alpha_\mathrm{min}$ is a factor $2\,$ τ_ulps above the step $\alpha^*$ at which the condition degenerates into a test decided by rounding — provided backtracking_αmin's upper clamp at $\sqrt{\mathrm{eps}(T)}$ is inactive, which it is for a merit of ordinary steepness in double precision. Where the clamp binds (a very flat merit, or any merit in Float16) the search does trial steps below $\alpha^*$. That is deliberate and harmless: the $\min$ in the SufficientDecreaseCondition means the test there reduces to $\varphi(\alpha) \leq \varphi_0$, i.e. plain monotonicity, and such an accept is reported as LINESEARCH_FLOOR rather than as a decrease.

BacktrackingCondition

Abstract type comprising the conditions that are used for checking step sizes for the backtracking line search (see Backtracking). This comprises the SufficientDecreaseCondition and the CurvatureCondition.

BierlaireQuadratic <: LinesearchMethod

Algorithm taken from [4].

Bisection <: LinesearchMethod

See bisection for the implementation of the algorithm.

Extended help

When invoked with a single trial step α (i.e. solve(ls, α)), the bracket is always lower-anchored at $\alpha = 0$ — the only point where a genuine descent direction is guaranteed to have a decreasing merit ($\varphi'(0) < 0$), which one-sided rightward bracketing requires. The caller's α is then folded in via one extra derivative evaluation (see issue #164):

• if $\varphi'(\alpha) \geq 0$ then α overshot the minimum and $[0, \alpha]$ already brackets a stationary point, so it is handed straight to bisection with no bracketing loop;
• otherwise $\alpha$ still lies on the descent side, so the bracket is grown outward from $0$ with the initial step seeded from $|\alpha|$ — clamped between DEFAULT_BRACKETING_s and 1 so a large α does not over-coarsen the search and a tiny α does not crawl — rather than the fixed default step.

This keeps the safe $\alpha = 0$ anchor while letting the caller's α set the search scale and, when it overshoots, serve directly as the upper bracket bound.

BracketMinimumCriterion <: BracketingCriterion

The criterion used for bracket_minimum. It checks whether $y(c)$ is greater than or equal to $y(b)$ (i.e. checks whether we are passed the minimum). Compare this with BracketRootCriterion.

Functor

bc = BracketMinimumCriterion()

yc = .1
yb = .2

bc(yb, yc)

# output

false

BracketRootCriterion <: BracketingCriterion

The criterion used for bracket_root. It checks whether there is a sign change between $b$ and $c$ (i.e. checks whether there is a root between those two points). Compare this with BracketMinimumCriterion.

Functor

bc = BracketRootCriterion()

yc = .1
yb = -.2

bc(yb, yc)

# output

true

BracketingCriterion

Abstract type for the criteria used while bracketing. It determines when a bracket has been found. The two concrete subtypes are BracketMinimumCriterion (used by bracket_minimum) and BracketRootCriterion (used by bracket_root).

CurvatureCondition <: BacktrackingCondition

The second of the Wolfe conditions [1]. The first one is the SufficientDecreaseCondition.

This encompasses the standard curvature condition and the strong curvature condition. This can be specified via the mode keyword.

With the standard curvature condition we check:

\[f'(\alpha) ≥ c_2 d,\]

where $c_2$ is the associated hyperparameter and $d$ is the derivative at $\alpha_0$. Further note that $f'(\alpha_0)$ and $d$ should both be negative.

With the strong curvature condition we check:

\[|f'(\alpha)| ≤ c_2 |d|.\]

Constructor

CurvatureCondition(c, d₀, D, Val(:Standard))
CurvatureCondition(c, d₀, D, Val(:Strong))

Here D has to be a function computing the derivative of the objective. The mode is passed as a Val (defaulting to Val(:Standard)) so that it is encoded in the type and dispatch — and hence inference — is stable without relying on constant propagation of a Symbol keyword. The other inputs are numbers.

DogLeg(refactorize=1)

Powell's dogleg method [5].

Like Newton, the refactorize parameter determines after how many steps the Jacobian is re-evaluated and refactored (see factorize!). The default refactorize = 1 re-evaluates and refactorizes the Jacobian on every step; refactorize > 1 reuses the Jacobian (and its factorization) in between, giving a quasi-Newton-style dogleg method.

DogLegCache

Like NonlinearSolverCache but storing two directions (callable with direction₁ and direction₂).

DogLegSolver

The NonlinearSolver for the DogLeg method.

Note that the DogLeg solver_step! is a trust-region method: it chooses the step length via the trust-region radius, not a line search. Unlike the NewtonSolver — but like the PicardSolver — no linesearch keyword is accepted (passing one is an error rather than being silently ignored).

Gradient

Abstract type. structs that are derived from this need an associated functor that computes the gradient of a function (in-place).

Examples

Examples include:

• GradientFunction
• GradientAutodiff
• GradientFiniteDifferences

GradientAutodiff <: Gradient

A struct that realizes Gradient by using ForwardDiff.

Keys

The struct stores:

• F: a function that has to be differentiated.
• ∇config: result of applying ForwardDiff.GradientConfig.

Constructors

GradientAutodiff(F, x::AbstractVector)
GradientAutodiff{T}(F, nx::Integer)

Functor

The functor does:

grad(g, x) = ForwardDiff.gradient!(g, grad.F, x, grad.∇config)

GradientFiniteDifferences <: Gradient

A struct that realizes Gradient by using finite differences.

Keys

The struct stores:

• F: a function that has to be differentiated.
• ϵ: small constant on whose basis the finite differences are computed.
• e: auxiliary vector used for computing finite differences. It's of the form $e_1 = \begin{bmatrix} 1 & 0 & \cdots & 0 \end{bmatrix}^T$.
• tx: auxiliary vector used for computing finite differences. It stores the offset in the x vector.

Constructor(s)

GradientFiniteDifferences{T}(F, nx::Integer; ϵ)

By default for ϵ is default_ϵ(T).

Functor

The functor does (for grad(g, x)):

for j in eachindex(x,g)
    ϵⱼ = grad.ϵ * abs(x[j]) + grad.ϵ
    fill!(grad.e, 0)
    grad.e[j] = 1
    grad.tx .= x .- ϵⱼ .* grad.e
    f1 = grad.F(grad.tx)
    grad.tx .= x .+ ϵⱼ .* grad.e
    f2 = grad.F(grad.tx)
    g[j] = (f2 - f1) / (2ϵⱼ)
end

GradientFunction <: Gradient

A struct that realizes a Gradient by explicitly supplying a function.

Keys

The struct stores:

• F: a function that has to be differentiated.
• ∇F!: a function that can be applied in place.

Functor

The functor does:

grad(g, x) = grad.∇F!(g, x)

Hessian

Abstract type. structs derived from this need an associated functor that computes the Hessian of a function (in-place).

Also see Gradient.

Implementation

When a custom Hessian is implemented, a functor is needed:

function (hessian::Hessian)(h::AbstractMatrix, x::AbstractVector) end

Examples

Examples include:

• HessianFunction
• HessianAutodiff

HessianAutodiff <: Hessian

A struct that realizes Hessian by using ForwardDiff.

Keys

The struct stores:

• F: a function that has to be differentiated.
• Hconfig: result of applying ForwardDiff.HessianConfig.

Constructors

HessianAutodiff{T}(F, Hconfig)
HessianAutodiff(F, x::AbstractVector)
HessianAutodiff{T}(F, nx::Integer)

Functor

The functor does:

hes(H, x) = ForwardDiff.hessian!(H, hes.F, x, hes.Hconfig)

HessianFunction <: Hessian

A struct that realizes a Hessian by explicitly supplying a function.

Keys

The struct stores:

• H!: a function that can be applied in place.

Functor

The functor does:

hes(H, x) = hes.H!(H, x)

Jacobian

Abstract type. structs that are derived from this need an associated functor that computes the Jacobian of a function (in-place).

Implementation

When a custom Jacobian is implemented, a functor is needed:

function (j::Jacobian)(g::AbstractMatrix, x::AbstractVector) end

Examples

Examples include:

• JacobianFunction
• JacobianAutodiff
• JacobianFiniteDifferences

Jacobian{T}(F, nx, ny; mode = :autodiff, kwargs...)

Construct a Jacobian of element type T for a function F mapping nx inputs to ny outputs, selecting the backend via the mode keyword:

• mode = :autodiff (default) builds a JacobianAutodiff (ForwardDiff).
• mode = :finitedifferences builds a JacobianFiniteDifferences; any remaining keyword arguments (e.g. ϵ) are forwarded to it.

The convenience forms Jacobian{T}(F, n), Jacobian(F, x) and Jacobian(F, x, y) forward here.

JacobianAutodiff <: Jacobian

A struct that realizes Jacobian by using ForwardDiff.

Keys

The struct stores:

• F: a function that has to be differentiated.
• Jconfig: result of applying ForwardDiff.JacobianConfig.
• ty: vector that is used for evaluating ForwardDiff.jacobian!

Constructors

JacobianAutodiff(F, x::AbstractVector)
JacobianAutodiff{T}(F, nx::Integer)

Functor

The functor does:

jac(J, x, params) = ForwardDiff.jacobian!(J, (y, x) -> jac.F(y, x, params), jac.ty, x, jac.Jconfig)

JacobianFiniteDifferences <: Jacobian

A struct that realizes Jacobian by using finite differences.

Keys

The struct stores:

• F: a function that has to be differentiated.
• ϵ: small constant on whose basis the finite differences are computed.
• f1: $f$ evaluated at $x - \epsilon_j e_j$ with $\epsilon_j = \epsilon|x_j| + \epsilon$ for all $j$.
• f2: $f$ evaluated at $x + \epsilon_j e_j$ with $\epsilon_j = \epsilon|x_j| + \epsilon$ for all $j$.
• e: auxiliary vector used for computing finite differences. It's of the form $e_1 = \begin{bmatrix} 1 & 0 & \cdots & 0 \end{bmatrix}^T$.
• tx: auxiliary vector used for computing finite differences. It stores the offset in the x vector.

Constructor(s)

JacobianFiniteDifferences{T}(F, nx::Integer, ny::Integer; ϵ)

By default for ϵ is default_ϵ(T).

Functor

The functor does:

for j in eachindex(x)
    ϵⱼ = jac.ϵ * abs(x[j]) + jac.ϵ
    fill!(jac.e, 0)
    jac.e[j] = 1
    jac.tx .= x .- ϵⱼ .* jac.e
    jac.F(jac.f1, jac.tx, params)
    jac.tx .= x .+ ϵⱼ .* jac.e
    jac.F(jac.f2, jac.tx, params)
    for i in eachindex(jac.f1)
        J[i,j] = (jac.f2[i] - jac.f1[i]) / (2ϵⱼ)
    end
end

JacobianFunction <: Jacobian

A struct that realizes a Jacobian by explicitly supplying a function taken from the NonlinearProblem.

Functor

f(y, x, params) = y .= [1. √2.; √2. 3.] * x
∇f(j, x, params) = j .= [1. √2.; √2. 3.]

jac = JacobianFunction(f, ∇f, Float64)
j = zeros(Float64, 2, 2)
x = ones(Float64, 2)

jac(j, x, NullParameters())

# output

2×2 Matrix{Float64}:
 1.0      1.41421
 1.41421  3.0

struct LU <: DirectMethod

A custom implementation of an LU solver, meant to solve a LinearProblem.

Routines that use the LU solver include factorize!, ldiv! and solve!.

Constructor

The constructor is called with either no argument:

LU()

# output

LU{Missing}(missing, true)

or with pivot and static as optional booleans:

LU(; pivot=true, static=true)

# output

LU{Bool}(true, true)

Note that if we do not supply an explicit keyword static, the corresponding field is missing (as in the first case). In that default case the cache matrix type is chosen by size via _static: a matrix whose leading dimension does not exceed N_STATIC_THRESHOLD yields a mutable static (MMatrix) cache, a larger one yields a plain Matrix. An explicit static=true/false forces the choice regardless of the matrix size.

Example

We use the LU together with solve to solve a linear system:

A = [1. 2. 3.; 5. 7. 11.; 13. 17. 19.]
v = rand(3)
ls = LinearProblem(A, v)

lu = LU()

solve(lu, ls) ≈ inv(A) * v

# output

true

Note that role of LinearProblem here.

LUSolverCache <: LinearSolverCache

The cache for the LU solver.

Keys

• A: the factorized matrix A,
• pivots: a vector of pivots used during factorization,
• perms: a vector of permutations used during factorization,
• info: stores an index regarding pivoting.

LinearProblem

A LinearProblem describes $Ax = y$, where we want to solve for $x$.

Keys

• A
• y

Constructors

A LinearProblem can be allocated by calling:

LinearProblem(A, y)
LinearProblem(A)
LinearProblem(y)
LinearProblem{T}(n, m)
LinearProblem{T}(n)

LinearProblem(A, y) stores copies of A and y, so the problem is ready to solve right after construction (and later mutations of the caller's arrays do not affect the stored copies):

A = [1. 2. 3.; 4. 5. 6.; 7. 8. 9.]
y = [1., 2., 3.]
ls = LinearProblem(A, y)

# output

LinearProblem{Float64, Vector{Float64}, Matrix{Float64}}([1.0 2.0 3.0; 4.0 5.0 6.0; 7.0 8.0 9.0], [1.0, 2.0, 3.0])

The size-only constructors (LinearProblem(A), LinearProblem(y), LinearProblem{T}(n[, m])) allocate the unspecified parts as NaNs; use update! to fill the system with values:

ls = LinearProblem(y)
update!(ls, A, y)

# output

LinearProblem{Float64, Vector{Float64}, Matrix{Float64}}([1.0 2.0 3.0; 4.0 5.0 6.0; 7.0 8.0 9.0], [1.0, 2.0, 3.0])

LinearSolver <: AbstractSolver

A struct that stores LinearSolverMethods (for example LU) and LinearSolverCaches (for example LUSolverCache). LinearSolvers are used to solve LinearProblems.

Constructors

LinearSolver(method, cache)
LinearSolver(method, A)
LinearSolver(method, ls::LinearProblem)
LinearSolver(method, x)

You can manually factorize by either calling factorize! or solve!.

LinearSolverCache

An abstract type that summarizes all the caches used for LinearSolvers. See e.g. LUSolverCache.

LinearSolverMethod <: SolverMethod

Summarizes all the methods used for solving linear systems of equations such as the LU method.

Extended help

The abstract type SolverMethod was imported from GeometricBase.

Linesearch

A struct that stores a LinesearchProblem, LinesearchMethod and Options.

Keys

• problem::LinesearchProblem
• method::LinesearchMethod
• config::Options

Constructors

The following constructors can be used:

Linesearch{T}(problem, method, config)
Linesearch(problem, method=Static(); kwargs...)
Linesearch(problem, method, config::Options)

LinesearchMethod{T} <: SolverMethod

Examples include Static, Backtracking, Bisection , BierlaireQuadratic and Quadratic. See these examples for specific information on linesearch algorithms.

Extended help

A LinesearchMethod is usually used in Linesearch (or with solve).

It is a subtype of SolverMethod (imported from GeometricBase) — line searches are one-dimensional subproblems used inside nonlinear solvers and optimizers, so (unlike a NonlinearSolverMethod) a LinesearchMethod is not itself a nonlinear-solver method.

The line search contract

Every method reached through solve or solve_with_status guarantees:

1. It never throws. A situation it cannot handle is reported, never raised — a line search must not abort the enclosing solve. Bracketing helpers signal failure with nothing (see bracket_minimum, triple_point_finder) and the method maps that onto a LinesearchOutcome.
2. It returns $\alpha > 0$. Never the $\alpha = 0$ anchor, which would freeze the outer iterate (x .+= 0 .* d), and never a negative step: $\alpha$ scales a direction that has already been chosen, so its sign is not the line search's to decide.
3. It reports through linesearch_warnings only — one message site and one verbosity policy for all methods (genuine failure at verbosity ≥ 1, rate limited; the benign round-off-floor and stationary outcomes at ≥ 2).
4. A non-finite or ascending anchor is reported, not assumed away — see check_anchor.
5. It terminates in a bounded number of merit evaluations, independently of the merit's scale. Multiplying $\varphi$ by a constant must not change the cost.

The two families

What is not standardised is the meaning of the input $\alpha$ and what each method guarantees about the step, because there are two distinct kinds:

• Condition-satisfying, $\alpha$-relative — Backtracking, StrongWolfe and trivially Static: "given the trial step $\alpha$, return a step satisfying a decrease condition". The result depends on the input $\alpha$.
• Minimising, $\alpha$-independent — Bisection, Quadratic and BierlaireQuadratic: "approximate the minimiser of $\varphi$ along the direction". The input $\alpha$ only seeds the bracketing (see issue #164), and no Wolfe condition is checked.

LinesearchOutcome

Why a LinesearchMethod stopped. Stored in a LinesearchStatus, which is returned by solve_with_status.

• LINESEARCH_DECREASED: a step was found that decreased the merit by more than the round-off allowance $\tau$. This is the only outcome that reports progress. Note what it does not claim: Backtracking and StrongWolfe additionally verify their Wolfe condition before returning, whereas the minimising searches (Bisection, Quadratic, BierlaireQuadratic) approximate the line minimiser and never test one. The common guarantee across all of them is the $\tau$-exceeding decrease.
• LINESEARCH_FLOOR: the merit has reached its round-off floor — no trial step changes it by more than $\tau$, so no line search can make progress here. The returned step is the smallest informative one. This is not an error and is only reported at verbosity ≥ 2: it is the expected final state of a converged solve, and when it is not — when the residual is still large — the outer iteration reports it as stagnation instead (see stalled_step and Options).
• LINESEARCH_EXHAUSTED: no acceptable step although the merit does vary by more than $\tau$. Either $\varphi'(0)$ is inconsistent with $\varphi$ (a stale or regularized Jacobian, an inexact linear solve, a non-smooth merit), or the linesearch_max_iterations budget of Options was spent.
• LINESEARCH_NO_DESCENT: $\varphi'(0) > 0$, or $\varphi(0)$/$\varphi'(0)$ is not finite. No $\alpha$ can satisfy the sufficient decrease condition.
• LINESEARCH_STATIONARY: $\varphi'(0) = 0$, e.g. a vanishing direction at an exact root. Benign — there is nothing to search for.
• LINESEARCH_UNKNOWN: the method does not report an outcome (the generic fallback of solve_with_status).

LinesearchProblem <: AbstractProblem

In practice LinesearchProblems are allocated by calling linesearch_problem.

Constructors

f(x) = x^2 - 1
g(x) = 2x
δx(x) = - g(x) / 2
x₀ = 3.
_f(α,_) = f(compute_new_iterate(x₀, α, δx(x₀)))
_d(α,_) = g(compute_new_iterate(x₀, α, δx(x₀)))
ls_obj = LinesearchProblem{typeof(x₀)}(_f, _d)

# output

LinesearchProblem{Float64, typeof(_f), typeof(_d)}(_f, _d)

LinesearchStatus{T}

The step length returned by a line search together with the diagnostics needed to tell progress from stagnation. Obtained from solve_with_status; compare this to NonlinearSolverStatus, which plays the same role for the outer iteration.

The step length alone cannot express the difference: a tiny $\alpha$ may be the correct answer, or it may be all that is left after the merit turned out to be irreducible. See LinesearchOutcome.

Keys

• α: the returned step length (the value solve returns),
• outcome::LinesearchOutcome,
• trials: the number of trial steps $\alpha > 0$ at which the method actually evaluated the problem in its own iteration — not the linesearch_max_iterations budget. That is the merit for every method except Bisection, which drives on the derivative it bisects. Evaluations spent inside a bracketing helper (bracket_minimum, triple_point_finder) are not included, so for the bracketing searches this is a lower bound on the total cost; for Backtracking and StrongWolfe it is exact, and every merit evaluation is either the $\alpha = 0$ anchor or a counted trial,
• φ₀, d₀: the merit and its derivative at the anchor $\alpha = 0$,
• φ: the merit at the returned step,
• τ: the round-off resolution of the merit (see armijo_tolerance), against which every method decides whether the decrease it achieved was genuine,
• αmin: the smallest step length that could still be decided by the merit rather than by rounding (see backtracking_αmin). This is a shrinking-ladder quantity and is therefore zero — meaning "not applicable" — for the minimising searches (Bisection, Quadratic, BierlaireQuadratic) and for StrongWolfe, which bracket rather than shrink.

LinesearchStatus(α, outcome=LINESEARCH_UNKNOWN)

Construct a LinesearchStatus that carries only the step length and the LinesearchOutcome; the remaining diagnostics are filled with NaN/zero. Used by the generic fallback of solve_with_status for methods that do not report them.

Newton(refactorize=1)

The Newton (and quasi-Newton) nonlinear solver method.

Constructors

Newton()

# output

Newton(1)

QuasiNewton()

# output

Newton(5)

NewtonSolver

A const derived from NonlinearSolver as NewtonSolver{T} = NonlinearSolver{T,Newton}.

Constructors

The NewtonSolver can be called with a NonlinearProblem or with a Callable.

See NewtonSolver(::AbstractVector{T}, ::Callable, ::AbstractVector{T}) where {T}.

F(y, x, params) = y .= sin.(x) ^ 2
x = ones(5)
y = zeros(5)

ns = NewtonSolver(x, F, y)
typeof(ns) <: NewtonSolver

# output

true

Keywords

• linear_solver_method: the method used to build the linear solver (see LinearSolver) that computes the direction of the solver step (see solver_step!),
• DF!: an in-place function computing the Jacobian,
• linesearch::Linesearch
• jacobian::Jacobian
• refactorize::Int: determines after how many steps the Jacobian is re-evaluated and refactored (see factorize!). refactorize > 1 gives a quasi-Newton method (see QuasiNewton),
• options_kwargs: see Options

NewtonSolver(x, F, y)

Keywords

• linear_solver_method
• DF!
• linesearch
• jacobian
• refactorize
• options_kwargs: see Options

A dummy linear system used for the fixed point iterator (Picard).

NonlinearProblem

A NonlinearProblem describes $F(x) = y$, where we want to solve for $x$ and $F$ is in nonlinear in general (also compare this to LinearProblem).

Keys

• F
• J::Union{Callable, Missing}: accessed by calling jacobian.

Constructors

We show an example for one particular constructor:

F(y, x, params) = y .= sin.(x) .^ 2
NonlinearProblem(F, zeros(3))

# output

NonlinearProblem{typeof(F), Missing}(F, missing)

NonlinearSolver

A struct that comprises Newton solvers (see Newton), the Picard solver (also known as fixed-point iteration; see Picard) and the Dogleg solver (see DogLeg).

Keys

• nonlinearproblem::NonlinearProblem: the system that has to be solved. This can be accessed by calling nonlinearproblem,
• linearproblem::LinearProblem,
• jacobian::Jacobian: the Jacobian is used to compute the direction in the solver step (see solver_step!). This can be accessed by calling jacobian,
• linearsolver::LinearSolver: the linear solver is used to compute the direction of the solver step (see solver_step!). This can be accessed by calling linearsolver,
• linesearch::Linesearch
• method::NonlinearSolverMethod: the solver method (e.g. Newton),
• cache::NonlinearSolverCache
• config::Options

NonlinearSolverCache <: AbstractNonlinearSolverCache

Derived from AbstractNonlinearSolverCache. Used in NonlinearSolver.

Keys

• x: the next iterate (or guess thereof),
• Δx: search direction. This is updated when calling solver_step! via the LinearSolver stored in the NewtonSolver,
• rhs: the right-hand-side (this can be accessed by calling rhs),
• y: the problem evaluated at x,
• j::AbstractMatrix: the Jacobian evaluated at x. Note that this is not of type Jacobian!

NonlinearSolverMethod <: SolverMethod

A supertype collecting all nonlinear solver methods, i.e. Newton, Picard and DogLeg.

Compare this with LinesearchMethod: both are subtypes of SolverMethod, but a LinesearchMethod describes a one-dimensional line search (used inside a solver step) whereas a NonlinearSolverMethod describes the outer nonlinear iteration itself.

NonlinearSolverState <: AbstractSolverState

The NonlinearSolverState to be used together with a NonlinearSolver.

Examples

julia> state = NonlinearSolverState(zeros(3))
NonlinearSolverState{Float64, Vector{Float64}, Vector{Float64}}(0, [NaN, NaN, NaN], [NaN, NaN, NaN], [NaN, NaN, NaN], [NaN, NaN, NaN], NaN, 0, false)

NonlinearSolverStatus

Stores absolute and successive residuals for x and f. It is used as a diagnostic tool in NewtonSolver.

Keys

• iterations: number of iterations
• stalls: number of consecutive stalled steps, see stalled_step and isstalled,
• rxₛ: successive residual in x,
• rfₐ: absolute residual in f,
• rfₛ: successive residual in f,
• x_converged::Bool
• f_converged::Bool
• f_increased::Bool
• stalled::Bool: the last step stalled, see stalled_step

Examples

x = [1., 2., 3., 4.]
state = NonlinearSolverState(x)
cache = NonlinearSolverCache(x, x)
config = Options()
NonlinearSolverStatus(state, config)

# output

i=   0,
rxₛ= NaN,
rfₐ= NaN,
rfₛ= NaN

Options

Examples

Options()

# output

                x_abstol = 4.440892098500626e-16
                x_reltol = 4.440892098500626e-16
                x_suctol = 4.440892098500626e-16
                f_abstol = 0.0
                f_reltol = 1.4901161193847656e-8
                f_suctol = 4.440892098500626e-16
                f_mindec = 0.0001
          f_abstol_break = Inf
       allow_f_increases = true
          min_iterations = 0
          max_iterations = 1000
         warn_iterations = 1000
linesearch_max_iterations = 60
              max_stalls = 2
              show_trace = false
             store_trace = false
          extended_trace = false
              show_every = 1
               verbosity = 1
      nan_max_iterations = 10
              nan_factor = 0.5
   regularization_factor = 0.0
   dogleg_radius_initial = 1.0
    dogleg_radius_shrink = 0.25
    dogleg_radius_expand = 2.0
       dogleg_radius_max = 100.0

Picard <: NonlinearSolverMethod

See PicardSolver.

PicardSolver(x, F)

Arguments

• x: the initial guess for the solution.
• F: the nonlinear function to solve.
• y

Keywords

• DF!: the Jacobian of F,
• jacobian: the Jacobian of F, defaults to JacobianAutodiff,
• options_kwargs: see Options.

Note that the Picard solver_step! is a residual-safeguarded fixed-point iteration and uses no line search, so — unlike the other solvers — no linesearch keyword is accepted (passing one is an error rather than being silently ignored).

Examples

F(y, x, params) = y .= sin.(x) .^ 2
x = zeros(2)
y = similar(x)

s = PicardSolver(x, F, y)
state = SolverState(s)

solve!(x, s, state)

# output

2-element Vector{Float64}:
 0.0
 0.0

Quadratic <: LinesearchMethod

Quadratic Polynomial line search based on the polynomial

\[p(α) = p_0 + p_1(\alpha - \alpha_0) + p_2(\alpha - \alpha_0)^2.\]

Performs multiple iterations in which all parameters $p_0$, $p_1$ and $p_2$ are adapted. We do not check the SufficientDecreaseCondition here. We instead repeatedly build new quadratic polynomials until a minimum is found (to sufficient accuracy). The iteration may also stop after it reaches the maximum number of iterations, the linesearch_max_iterations field of Options (see linesearch_iterations).

Keywords

• ε: A constant that checks the precision/tolerance.
• s: A constant that determines the initial interval for bracketing. By default this is DEFAULT_BRACKETING_s.
• s_reduction: A constant that determines the factor by which s is decreased in each new bracketing iteration.

Extended help

The quadratic method. Compare this to BierlaireQuadratic. The algorithm is adjusted from [6].

Static <: LinesearchMethod

The static method.

Keys

Keys include:

• α: equivalent to a step size. The default is 1.

Examples

Static()

# output

Static with α = 1.0.

StrongWolfe{T} <: LinesearchMethod

A line search that finds a step $\alpha$ satisfying the strong Wolfe conditions

\[\begin{aligned} f(\alpha) &\leq f(0) + c_1\,\alpha\,f'(0), &\text{(sufficient decrease / Armijo)}\\ |f'(\alpha)| &\leq c_2\,|f'(0)|, &\text{(strong curvature)} \end{aligned}\]

with $0 < c_1 < c_2 < 1$. It implements the bracketing line search of [1, Alg. 3.5 and 3.6 (zoom)]: a bracketing phase grows the step until an interval containing an acceptable point is found, then a zoom phase shrinks that interval (by bisection) until the strong Wolfe conditions hold.

Unlike Backtracking — which enforces only sufficient decrease, since the curvature condition cannot be honoured by shrinking alone — StrongWolfe actually enforces the curvature condition, at the cost of evaluating the derivative at each trial step. Use it when curvature control is genuinely required; Backtracking is cheaper otherwise.

Keys

• c₁ (default DEFAULT_WOLFE_c₁): the Armijo constant $c_1$.
• c₂ (default DEFAULT_WOLFE_c₂): the curvature constant $c_2$. We require $c_1 < c_2 < 1$.
• αmax (default DEFAULT_WOLFE_αmax): the largest step the bracketing phase will try.

SufficientDecreaseCondition <: BacktrackingCondition

The condition that determines if the change induced by $\alpha_k$ is big enough. This is used in Backtracking.

Example

c = SimpleSolvers.DEFAULT_WOLFE_c₁
f(x) = (x - 1.) ^ 2
xₖ = 0.
fₖ = f(xₖ)
dₖ = 2xₖ - 2.

sdc = SufficientDecreaseCondition(c, fₖ, dₖ, f)
sdc(1.9), sdc(2.)

# output

(true, false)

Extended help

We call the constant that pertains to the sufficient decrease condition $c$. This is typically called $c_1$ in the literature [1]. See DEFAULT_WOLFE_c₁ for the relevant constant

The optional keyword τ slackens the condition by an absolute amount:

\[f(\alpha) \leq \min\{f_0,\ f_0 + c\alpha{}d_0 + \tau\}.\]

It defaults to zero, i.e. the exact condition. Without it the accept/reject decision is taken by rounding alone as soon as $c\alpha|d_0|$ drops below one unit in the last place of $f_0$: the right-hand side then rounds back up to $f_0$ and the test degenerates to $f(\alpha) \leq f_0$, which a merit that has reached its round-off floor passes or fails at random. See armijo_tolerance and backtracking_αmin for how Backtracking chooses $\tau$ and derives a meaningful smallest step from it.

The $\min$ bounds the slackening: $\tau$ may reduce the decrease that is demanded, but it never accepts a step whose merit exceeds $f_0$. For $d_0 < 0$ — which both callers guarantee via check_anchor — the $\min$ is inactive wherever $f_0 + c\alpha{}d_0$ is representably below $f_0$, so it changes nothing in double precision; it matters at low precision, where $\tau$ can exceed the demanded $c\alpha|d_0|$ outright.

update!(ls::LinearProblem, A, b)

Set the rhs vector to b and the matrix stored in the LinearProblem ls to A.

update!(state, x, y)

Update x̄, ȳ, x and y.

Examples

julia> f(y, x, params) = y .= sin.(x .- .5) .^ 2
f (generic function with 1 method)

julia> x = ones(1) / 4
1-element Vector{Float64}:
 0.25

julia> y = zero(x); f(y, x, NullParameters())
1-element Vector{Float64}:
 0.06120871905481365

julia> state = NonlinearSolverState(x)
NonlinearSolverState{Float64, Vector{Float64}, Vector{Float64}}(0, [NaN], [NaN], [NaN], [NaN], NaN, 0, false)

julia> update!(state, x, y)
NonlinearSolverState{Float64, Vector{Float64}, Vector{Float64}}(0, [0.25], [NaN], [0.06120871905481365], [NaN], NaN, 0, false)

julia> x = ones(1) / 2
1-element Vector{Float64}:
 0.5

julia> f(y, x, NullParameters())
1-element Vector{Float64}:
 0.0

julia> update!(state, x, y)
NonlinearSolverState{Float64, Vector{Float64}, Vector{Float64}}(0, [0.5], [0.25], [0.0], [0.06120871905481365], NaN, 0, false)

The NonlinearSolverState stores the previous solution, the previous value, the current solution and the current value.

All of these are updated during one update! step (and initialized with NaNs).

ldiv!(x, lsolver, b)

Compute inv(cache(lsolver).A) * b by utilizing the factorization of the lu solver (see LU and LinearSolver) and store the result in x.

Examples

julia> A = [1.; 0.; 0.;; 0.; 2.; 0.;; 0.; 0.; 4.]
3×3 Matrix{Float64}:
 1.0  0.0  0.0
 0.0  2.0  0.0
 0.0  0.0  4.0

julia> b = [1., 1., 1.]
3-element Vector{Float64}:
 1.0
 1.0
 1.0

julia> s = LinearSolver(LU(), A); factorize!(s); x = zeros(3)
3-element Vector{Float64}:
 0.0
 0.0
 0.0

julia> ldiv!(x, s, b)
3-element Vector{Float64}:
 1.0
 0.5
 0.25

QuasiNewton(refactorize=5)

Convenience constructor for a Newton method whose Jacobian is only re-evaluated and refactored every refactorize iterations. Equivalent to Newton(refactorize) but with a quasi-Newton default (see DEFAULT_ITERATIONS_QUASI_NEWTON_SOLVER).

_static(A)

Determine whether the LUSolverCache for a default LU should store A as a mutable static matrix (MMatrix) or as a plain Matrix. Every matrix whose leading dimension is smaller than or equal to N_STATIC_THRESHOLD is stored as an MMatrix.

This is only consulted for the default LU() (i.e. LU{Missing}); an explicit static=true/false keyword overrides it. See the examples in factorize!.

absolute_tolerance(T)

Determine the absolute tolerance for a specific data type. This is used in the constructor of Options.

In comparison to default_tolerance, this should return a very small number, close to zero (i.e. not just machine precision).

Examples

julia> absolute_tolerance(Float64)
0.0

julia> absolute_tolerance(Float32)
0.0f0

alloc_g(x)

Allocate NaNs of the size of the gradient of f (with respect to x).

alloc_h(x)

Allocate NaNs of the size of the Hessian of f (with respect to x).

alloc_x(x)

Allocate NaNs of the size of x.

armijo_tolerance(φ₀, n)

The absolute round-off resolution $\tau = n\cdot\mathrm{ulp}(\varphi_0)$ of a merit function, where n is a number of units in the last place — armijo_ulps for the default. See DEFAULT_ARMIJO_τ_ULPS for what it is used for and backtracking_αmin for the step length derived from it.

This lives here rather than next to Backtracking because triple_point_finder needs it too, and bracketing/ is included before linesearch/.

armijo_ulps(T, c₁)
armijo_ulps(T)

The number of ulps of $\varphi(0)$ to use as the round-off resolution $\tau$ for element type T: the nominal DEFAULT_ARMIJO_τ_ULPS, capped at what T can support. The one-argument form uses DEFAULT_WOLFE_c₁.

$\tau$ has to satisfy two requirements that pull in opposite directions. To recognise a merit sitting at its round-off floor it must be at least an ulp or so of $\varphi(0)$. To leave the SufficientDecreaseCondition meaningful it must be far below the decrease that condition demands, which for the canonical $\|F\|^2$ merit of a Newton step ($\varphi'(0) = -2\varphi(0)$) is $2c_1\varphi(0)$ at $\alpha = 1$. Hence the cap

\[n \leq \frac{\mathtt{ARMIJO\_τ\_DEMAND\_FRACTION}\cdot 2c_1}{\mathrm{eps}(T)} .\]

The two requirements are compatible only while $\mathrm{eps}(T) \ll 2c_1$. They are, by a wide margin, in double precision (the cap is $\sim10^{10}$) and comfortably in single ($\sim17$), so the nominal 4 stands in both. They are not in Float16, where $\mathrm{eps}(T) = 9.8 \cdot 10^{-4}$ already exceeds $2c_1 = 2\cdot10^{-4}$: no value of τ_ulps above zero can satisfy both, so the cap resolves the conflict in favour of a meaningful condition and drops $\tau$ to about $2\cdot10^{-3}$ ulps — in effect the exact condition.

Nothing is lost by that. The floor is still detected, by two mechanisms that do not depend on $\tau$: at a trial step small enough that $\mathrm{fl}(\varphi(0) + c_1\alpha\varphi'(0))$ rounds back to $\varphi(0)$ the condition is $\varphi(\alpha) \leq \varphi(0)$, and Backtracking additionally stops on two consecutive bit-identical merits. What is gained is that a genuine decrease of one or two ulps — the smallest a Float16 merit can express — is reported as LINESEARCH_DECREASED instead of as LINESEARCH_FLOOR, which the outer iteration would otherwise count towards max_stalls.

Examples

julia> armijo_ulps(Float64), armijo_ulps(Float32)
(4.0, 4.0f0)

julia> armijo_ulps(Float16)
Float16(0.002075)

assess_convergence(rxₛ, rfₐ, rfₛ, config, state)

Assess convergence for status::NonlinearSolverStatus and return the triple (x_converged, f_converged, f_increased).

The successive-change criteria (in x and f) alone are not sufficient to declare convergence: a stalled step (e.g. an artificially tiny line-search step) makes the successive residuals rxₛ and rfₛ vanish even when the absolute residual rfₐ is large. We therefore require the residual to be small — written residual_small below — in addition to the successive-change criterion before reporting convergence. The residual passes the standard atol + rtol·‖F₀‖ test:

residual_small ⟺ rfₐ ≤ config.f_abstol + config.f_reltol * initial_residual(state),

with the absolute tolerance atol = config.f_abstol (defaulting to 0) and the relative tolerance rtol = config.f_reltol (defaulting to √eps(T)) applied to the initial residual ‖F(x₀)‖. Concretely:

• x_converged: rxₛ ≤ norm(solution(state)) * config.x_suctol and residual_small,
• f_converged: (rfₛ ≤ norm(value(state)) * config.f_suctol and residual_small) or rfₐ ≤ config.f_abstol,
• f_increased: norm(value(state)) > norm(previousvalue(state)).

This guards the successive-change criteria against stagnation: it is loose enough that a genuinely converged iterate satisfies it (the successive-change criteria still supply the tight, machine-precision accuracy) yet tight enough to reject a step that stalls near its initial residual (rfₐ ≈ ‖F(x₀)‖ ≫ f_reltol·‖F(x₀)‖). The relative term is what lets a well-scaled solve whose residual floors at a large absolute value (e.g. a large-magnitude or ill-conditioned F) still converge; it drops to zero (leaving the pure absolute f_abstol test) until the state has been initialized (initial_residual is NaN).

Also see meets_stopping_criteria.

backtracking_interpolation(φ₀, d₀, α, φα, αp, φp, p)

The next trial step of the safeguarded polynomial backtracking used by Backtracking (see [1, §3.5], [2, Alg. A6.3.1]).

α/φα is the trial step that was just rejected and αp/φp the one rejected before it (αp is NaN on the first backtrack). The model interpolates $\varphi(0)$, $\varphi'(0)$ and the rejected value(s) — a quadratic on the first backtrack, a cubic afterwards — and its minimiser is clamped to $[$ BACKTRACKING_SHRINK_MIN $\cdot\alpha, p\alpha]$.

The clamp is what makes this safe: an unclamped interpolant can return $\alpha$ itself (no progress at all), collapse to numerically zero, or be meaningless because the merit values it is built from are rounding noise. Because the upper bound is $p$, the trial sequence is pointwise never longer than the plain $\alpha \gets p\alpha$ ladder.

backtracking_αmin(c₁, d₀, τ)

The smallest step length for which the SufficientDecreaseCondition can still be decided by the merit rather than by rounding:

\[\alpha_\mathrm{min} = \frac{\tau}{c_1|\varphi'(0)|} .\]

Below $\alpha_\mathrm{min}$ the demanded decrease $c_1\alpha|\varphi'(0)|$ is smaller than the round-off resolution $\tau$, so a trial step carries no information. Writing $\tau = n\cdot\mathrm{ulp}(\varphi(0))$ (see armijo_tolerance) gives

\[\alpha_\mathrm{min} = 2n\,\alpha^*, \qquad \alpha^* = \frac{\mathrm{ulp}(\varphi(0))}{2c_1|\varphi'(0)|},\]

where $\alpha^*$ is the step below which $\mathrm{fl}(\varphi(0) + c_1\alpha\varphi'(0))$ rounds back up to $\varphi(0)$ and the condition degenerates to $\varphi(\alpha) \leq \varphi(0)$.

The result is clamped to $[\mathrm{eps}(T), \sqrt{\mathrm{eps}(T)}]$: the lower bound is the historical negligible-step floor, and the upper bound makes sure that a nearly flat but genuine merit (very small $|\varphi'(0)|$) is still searched — an unclamped $\alpha_\mathrm{min}$ grows without bound as $|\varphi'(0)| \to 0$ and would stop the search before it began.

bisection(f, αmin, αmax, params, config)

Perform bisection of f in the interval [αmin, αmax] with Options config.

The algorithm is repeated until a root is found (up to tolerance config.f_abstol which is determined by default_tolerance by default).

Extended help

The bisection algorithm divides an interval into equal halves until a root is found (up to a desired accuracy).

We first initialize:

\[\begin{aligned} \alpha_0 \gets & \alpha_\mathrm{min}, \\ \alpha_1 \gets & \alpha_\mathrm{max}, \end{aligned}\]

and then repeat:

\[\begin{aligned} & \alpha \gets \frac{\alpha_0 + \alpha_1}{2}, \\ & \text{if $f(\alpha_0)f(\alpha) > 0$} \\ & \qquad \alpha_0 \gets \alpha, \\ & \text{else} \\ & \qquad \alpha_1 \gets \alpha, \\ & \text{end} \end{aligned}\]

So the algorithm checks in each step where the sign change occurred and moves the $\alpha_0$ or $\alpha_1$ accordingly. The loop is terminated if config.linesearch_max_iterations is reached (by default 60 for Float64 in the Options struct, see linesearch_iterations); in that case a warning is emitted (at verbosity ≥ 1) and the best estimate found so far is returned.

bracket(f, x, bc, s, k, nmax)

Grow a bracket outward from x (in steps scaled by k, starting from s) until the BracketingCriterion bc is satisfied. Used by bracket_minimum and bracket_root.

Extended help

Before entering the main loop we check whether the criterion is already satisfied just to the left of a (at a - s). This early exit is only valid for the BracketRootCriterion, where it corresponds to a sign change in (a - s, b). For the BracketMinimumCriterion it would instead bracket a maximum rather than a minimum, so it is deliberately skipped.

Returns nothing when no bracket is found within nmax steps. A line search must be able to report an unbracketable merit rather than abort the enclosing solve, so this is a nothing rather than an error; see bracket_minimum.

bracket_minimum(f, x)

Move a bracket successively in the search direction (starting at x) and increase its size until a local minimum of f is found.

This is used in bisections when only one x is given (and not an entire interval).

This bracketing algorithm is taken from [3]. Also compare it to bracket_minimum_with_fixed_point.

Arguments

• f: the function to be bracketed,
• x: the starting point,
• s: by default DEFAULT_BRACKETING_s,
• k: by default DEFAULT_BRACKETING_k,
• nmax: by default DEFAULT_BRACKETING_nmax.

Extended help

For bracketing we need two constants $s$ and $k$ (see DEFAULT_BRACKETING_s and DEFAULT_BRACKETING_k).

Before we start the algorithm we initialize it, i.e. we check that we indeed have a descent direction:

\[\begin{aligned} & a \gets x, \\ & b \gets a + s, \\ & \mathrm{if} \quad f(b) > f(a)\\ & \qquad\text{Flip $a$ and $b$ and set $s\gets-s$.}\\ & \mathrm{end} \end{aligned}\]

The algorithm then successively computes:

\[c \gets b + s,\]

and then checks whether $f(c) \geq f(b)$ (also see BracketMinimumCriterion). If this is true it returns $(a, c)$ or $(c, a)$, depending on whether $a<c$ or $c<a$ respectively. If this is not satisfied $a,$ $b$ and $s$ are updated:

\[\begin{aligned} a \gets & b, \\ b \gets & c, \\ s \gets & sk, \end{aligned}\]

and the algorithm is continued. If we have not found a bracket after $n_\mathrm{max}$ iterations (see DEFAULT_BRACKETING_nmax) the algorithm terminates and returns nothing. The interval that is returned by bracket_minimum is then typically used as a starting point for bisection.

bracket_minimum_with_fixed_point(f, x, s, k, nmax)

Find a bracket while keeping the left side (i.e. x) fixed.

The algorithm is similar to bracket_minimum (also based on DEFAULT_BRACKETING_s and DEFAULT_BRACKETING_k) with the difference that for the latter the left side is also moving.

The function bracket_minimum_with_fixed_point is used as a starting point for Quadratic (adapted from [6]), as the coefficient $p_2$ of the fitted polynomial is:

\[p_2 = \frac{f(b) - f(a) - f'(a)b}{b^2},\]

where $b = \mathtt{bracket\_minimum\_with\_fixed\_point}(a)$. The right end b is grown outward (with the left end a held fixed) until f stops decreasing, i.e. until the turning point f(b) ≥ f(b_\mathrm{prev}) is reached, so that a minimum is bracketed in (a, b). (The earlier variant compared against the fixed anchor f(a) instead, which failed to bracket a minimum whose right tail stays below f(a).) The Quadratic caller guards the fitted curvature (p_2 ≤ 0 falls back to a bisection step), so f(b) > f(a) is no longer required.

Returns the bracket together with the function values at its endpoints, (a, b, f(a), f(b)) with a < b. The values are already computed during bracketing, so the caller (the Quadratic line search) does not have to re-evaluate f at the endpoints.

Returns nothing if no bracket is found within nmax steps — a line search must be able to report an unbracketable merit rather than abort the enclosing solve.

bracket_root(f, x)

Make a bracket for the function based on x (for root finding).

This is largely equivalent to bracket_minimum. See the end of that docstring for more information.

cache(ls)

Return the cache of the LinearSolver.

Examples

For the default LU(), a small matrix (leading dimension ≤ N_STATIC_THRESHOLD) is stored as a mutable static matrix (MMatrix):

julia> ls = LinearSolver(LU(), [1.0 2.0; 3.0 4.0]);

julia> cache(ls)
SimpleSolvers.LUSolverCache{Float64, StaticArraysCore.MMatrix{2, 2, Float64, 4}}([1.0 2.0; 3.0 4.0], [0, 0], [0, 0], 0)

Passing static=false forces a plain Matrix cache regardless of size:

julia> ls = LinearSolver(LU(; static=false), [1.0 2.0; 3.0 4.0]);

julia> cache(ls)
SimpleSolvers.LUSolverCache{Float64, Matrix{Float64}}([1.0 2.0; 3.0 4.0], [0, 0], [0, 0], 0)

change_precision(T, method::LinesearchMethod)

Return a copy of the LinesearchMethod method with its numeric fields converted to the element type T.

This is an internal helper used when constructing a Linesearch: the method's precision is adapted to the working precision T. It replaces a former misuse of Base.convert (which was ambiguous with Base and violated the convert contract by returning a differently-typed object).

check_anchor(φ₀, d₀, α)

Validate the $\alpha = 0$ anchor of a line search problem. Return a LinesearchStatus that the caller should return immediately, or nothing if the anchor is usable and the search may proceed.

This is the one definition of the anchor policy shared by every LinesearchMethod:

• $\varphi(0)$ or $\varphi'(0)$ not finite, or $\varphi'(0) > 0$, gives LINESEARCH_NO_DESCENT: no $\alpha$ can decrease the merit along this direction, so shrinking or bracketing would only spend merit evaluations to discover that. The caller's trial step α is handed back — never the $\alpha = 0$ anchor, which would freeze the outer iterate (x .+= 0 .* d).
• $\varphi'(0) = 0$ gives LINESEARCH_STATIONARY. For the $\|F\|^2$ merit of a NonlinearSolver this is the exact root ($F = 0 \Rightarrow \varphi'(0) = 0$ and the direction vanishes), so it is benign and every $\alpha$ is equivalent.

An ascent anchor arises in practice when the direction did not come from an exact, freshly factorized Newton solve — a stale Jacobian under refactorize > 1, a nonzero regularization_factor, or an inexact linear solve. The correct response is to refresh the Jacobian, which is why the line search reports the situation instead of trying to salvage a step from it, and solver_step! acts on the report: on LINESEARCH_NO_DESCENT it leaves the iterate where it is (moving along a direction that cannot decrease the merit would only make the retry start from a worse point) and records a stall, which forces a fresh Jacobian on the next step (see needs_refresh and maybe_refactorize!) and gives up after max_stalls if that does not help.

The step handed back is therefore still positive, as the contract requires — whether to use it is the caller's decision, not the line search's.

check_gradient([io], g)

Check norm, maximum value and minimum value of a vector.

Output is written to io (defaulting to stdout).

Examples

julia> g = [1., 1., 1., 2., 0.9, 3.];

julia> SimpleSolvers.check_gradient(g; digits=3)
norm(Gradient):               4.1
minimum(|Gradient|):          0.9
maximum(|Gradient|):          3.0

check_hessian([io], H)

Check the condition number, determinant, max and min value of the Hessian H.

Output is written to io (defaulting to stdout).

julia> H = [1. √2.; √2. 3.];

julia> SimpleSolvers.check_hessian(H)
Condition Number of Hessian: 13.9282
Determinant of Hessian:      1.0
minimum(|Hessian|):          1.0
maximum(|Hessian|):          3.0

check_jacobian([io], J)

Check the condition number, determinant, max and min value of the Jacobian J.

Output is written to io (defaulting to stdout).

julia> J = [1. √2.; √2. 3.];

julia> SimpleSolvers.check_jacobian(J)
Condition Number of Jacobian: 13.9282
Determinant of Jacobian:      1.0
minimum(|Jacobian|):          1.0
maximum(|Jacobian|):          3.0

clear!(ls)

Write NaNs into matrix(ls) and rhs(ls).

Here ls is a LinearProblem.

compute_new_iterate!(xₖ₊₁, xₖ, αₖ, pₖ)

Compute xₖ₊₁ based on a direction pₖ and a step length αₖ.

Extended help

In the case of vector spaces this function simply does:

xₖ = xₖ + αₖ * pₖ

For manifolds we instead perform a retraction [7].

curvature_diagnostic(status, ls, params)

Method-specific extra diagnostic emitted by linesearch_warnings at verbosity ≥ 2. The fallback does nothing; Backtracking checks the CurvatureCondition, which costs a derivative evaluation — a full Jacobian for the line search problem of a NonlinearSolver, hence the verbosity gate.

default_precision(T)

Compute the default precision used for e.g. BierlaireQuadratic.

Compare this to the default_tolerance used in Options.

Examples

julia> default_precision(Float64)
1.7763568394002505e-15

julia> default_precision(Float32)
9.536743f-7

julia> default_precision(Float16)
Float16(0.007812)

default_tolerance(T)

Determine the default tolerance for a specific data type. This is used in the constructor of Options.

Compare this to default_precision.

Examples

julia> default_tolerance(Float64)
4.440892098500626e-16

julia> default_tolerance(Float32)
2.3841858f-7

julia> default_tolerance(Float16)
Float16(0.001953)

default_ϵ(::Type{T})

The default step size on whose basis finite differences are computed, for the working precision T. Used by GradientFiniteDifferences and JacobianFiniteDifferences.

Its value is $8\sqrt{\varepsilon_T}$, where $\varepsilon_T$ is the machine epsilon of T. Being precision-aware (eps(T), not a baked-in Float64 epsilon) is essential for Float32 finite differences to be accurate.

Examples

julia> default_ϵ(Float64)
1.1920928955078125e-7

julia> default_ϵ(Float32)
0.0027621358f0

direction!(d, x, s, params, iteration; stalled=false)

Compute the Newton direction (for the NewtonSolver). stalled is forwarded to maybe_refactorize!; see needs_refresh.

direction(cache)

Return the direction (i.e. the step vector $\Delta{}x$) stored in a solver cache such as NonlinearSolverCache or DogLegCache.

directions!(s, x, params, iteration=1; force_refactorize=false)

Compute direction₁ and direction₂ for the DogLegSolver.

This is equivalent to direction! for the NewtonSolver.

Examples

julia> J = [0 1; -1 0];

julia> f(y, x, params) = y .= cos.(J * x .- 2.) .^ 2 / l2norm(sin.(x) .- 1.);

julia> x = zeros(2); y = similar(x); s = DogLegSolver(x, y; F = f);

julia> directions!(s, x, NullParameters());

julia> direction₁(cache(s))
2-element Vector{Float64}:
 -0.25513686072399455
  0.1601152321012896

julia> direction₂(cache(s))
2-element Vector{Float64}:
 -0.22882877718014286
  0.22882877718014288

Extended help

The Gauss-Newton direction (i.e. direction₂) is computed the usual way:

\[\mathbf{d}_2 = -\mathbf{J}^{-1} \mathbf{r}\]

where $\mathbf{J}$ is the Jacobian matrix and $\mathbf{r}$ is the residual vector. The steepest descent direction (taken from [1, Equation (11.46)]) is different:

\[\mathbf{d}_1 = -\frac{||\mathbf{J}^T\mathbf{r}||^2}{\mathbf{r}^T(\mathbf{J}\mathbf{J}^T)(\mathbf{J}\mathbf{J}^T)\mathbf{r}}\mathbf{J}^T\mathbf{r}.\]

The DogLegSolver then interpolates between these two directions (this interpolation is piecewise linear).

As for the (quasi-)NewtonSolver, the Jacobian is only re-evaluated and refactored every refactorize iterations (see DogLeg), and always on a fresh state or the first step (iteration ≤ 1), or when force_refactorize = true (used by solver_step! to recover from a collapsed trust-region radius, and after any step that did not move the iterate — see needs_refresh). In between, the stale Jacobian and its factorization are reused for both directions. The default refactorize = 1 refactorizes on every step.

direction₁(cache::DogLegCache)

Return the steepest descent direction.

See directions!.

direction₂(cache::DogLegCache)

Return the Newton direction.

See directions!.

dogleg_direction!(cache, Δ)

Compute the (piecewise-linear) dogleg step for trust-region radius Δ from the steepest-descent direction direction₁ and the Newton direction direction₂ (both already stored in cache), writing the result into direction(cache).

direction₁ and direction₂ do not depend on Δ, so this may be called repeatedly while shrinking Δ without recomputing (and refactorizing) the Jacobian.

factorize!(lsolver::LinearSolver, A)

Factorize the matrix A and store the result in cache(lsolver).A.

Note that calling cache on lsolver returns the instance of LUSolverCache stored in lsolver.

Examples

julia> A = [1. 2. 3.; 5. 7. 11.; 13. 17. 19.]
3×3 Matrix{Float64}:
  1.0   2.0   3.0
  5.0   7.0  11.0
 13.0  17.0  19.0

julia> x = zeros(3);

julia> lsolver = LinearSolver(LU(; static=false), x);

julia> factorize!(lsolver, A).cache.A
3×3 Matrix{Float64}:
 13.0        17.0       19.0
  0.0769231   0.692308   1.53846
  0.384615    0.666667   2.66667

julia> y = A * ldiv!(x, lsolver, ones(3));

julia> round.(y; digits = 10)
3-element Vector{Float64}:
 1.0
 1.0
 1.0

Here cache(lsolver).A stores the factorized matrix. If we call factorize! with two input arguments as above, the method first copies the matrix A into the LUSolverCache. We can equivalently also do:

julia> lsolver = LinearSolver(LU(), A);

julia> factorize!(lsolver).cache.A
3×3 StaticArraysCore.MMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
 13.0        17.0       19.0
  0.0769231   0.692308   1.53846
  0.384615    0.666667   2.66667

Note the difference between the output types of the two refactorized matrices: the default LU() chose a mutable static (MMatrix) cache because the matrix is small (see _static and N_STATIC_THRESHOLD), whereas LU(; static=false) forced a plain Matrix.

Also see ldiv! for how the refactorized matrix is used.

flag_stall!(state)

Record that the line search of the current step reported that it cannot make progress along the current direction — either the merit is at its round-off floor (isfloor) or the anchor is not a descent direction at all (LINESEARCH_NO_DESCENT, see LinesearchOutcome). The flag is OR-ed into the verdict of the next record_stall!, which clears it again, and it makes needs_refresh true for the next step.

This is how a line search that knows it cannot help reports one iteration earlier than the step-based diagnosis of stalled_step — which remains the primary mechanism, since it is the only one that also covers a Static step along an underflowed direction, a collapsed DogLegSolver trust-region radius, and a locally expanding PicardSolver map.

increase_iteration_number!(state)

To be used together with NonlinearSolverState.

initial_residual(state)

Return the initial residual ‖F(x₀)‖ recorded by initialize! (NaN if the state has not been initialized). This is the reference scale for the relative-residual convergence test in assess_convergence.

initialize!(ls, x)

Initialize the LinearProblem ls.

This uses clear!(::LinearProblem).

initialize!(cache, x)

Initialize the NonlinearSolverCache with NaNs.

isconverged(status)

Check if either x or f has converged.

The status is a NonlinearSolverStatus.

isfloor(status)

true if the line search could not find any step that changes the merit by more than the round-off allowance τ, i.e. the merit has reached its round-off floor. The outer iteration cannot make progress in this state no matter how the step is chosen — which is why a NonlinearSolver counts it as a stalled step (see record_stall!).

isstalled(status, config)

Check whether the iteration has stagnated: config.max_stalls consecutive steps stalled_step.

Mutually exclusive with isconverged — a stalled step is by definition one whose residual is not small, whereas both convergence branches require that it is.

A stagnated solve has reached the numerical floor of its residual and cannot improve it. Whether that counts as success is the caller's decision, which is why the status is queryable (see status): if status.rfₐ is acceptable to you, treat isstalled as success — and consider raising f_abstol above it, since the tolerance you asked for is not attainable.

issufficient(status)

true if the line search found a step with a genuine sufficient decrease, i.e. one that decreased the merit by more than the round-off allowance τ of status. Compare isfloor.

iterate_settled(rxₛ, config, state)

Return true when the last step did not move the iterate, rxₛ ≤ ‖x‖·x_suctol. Used by assess_convergence and stalled_step.

iteration_number(state)

Return the number of iterations taken so far, as counted by increase_iteration_number! on state::NonlinearSolverState. Compare this to stall_number.

jacobian(nlp)

Return the Jacobian function stored in the NonlinearProblem nlp.

jacobianmatrix(solver::NonlinearSolver)

Return the evaluated Jacobian (a matrix) stored in the cache (see NonlinearSolverCache) of solver.

Also see jacobian(::NonlinearProblem).

linearsolver(solver)

Return the linear part (i.e. a LinearSolver) of an NewtonSolver.

Examples

x = rand(3)
y = rand(3)
F(x) = tanh.(x)
F!(y, x, params) = y .= F(x)
s = NewtonSolver(x, y; F = F!)
linearsolver(s)

# output

LinearSolver{Float64, LU{Missing}, SimpleSolvers.LUSolverCache{Float64, StaticArraysCore.MMatrix{3, 3, Float64, 9}}}(LU{Missing}(missing, true), SimpleSolvers.LUSolverCache{Float64, StaticArraysCore.MMatrix{3, 3, Float64, 9}}([0.0 0.0 0.0; 0.0 0.0 0.0; 0.0 0.0 0.0], [0, 0, 0], [0, 0, 0], 0))

linesearch_iterations(T)

Determine the default number of trial steps a line search may take, i.e. the default of the linesearch_max_iterations field of Options.

This is deliberately not the same quantity as max_iterations, which bounds the outer nonlinear iteration (see meets_stopping_criteria). A one-dimensional search inside a single solver step needs a budget on the order of the mantissa width, not thousands of trials: a Backtracking ladder $\alpha \gets p\alpha$ starting at $\alpha_0 = 1$ reaches the negligible-step floor after $\lceil-\log_2\varepsilon\rceil$ halvings (52 in double precision, 24 in single), and a bisection needs the same count to exhaust the mantissa. We take that count plus a small margin; everything beyond it can only produce denormals.

The count is derived for the default shrink factor $p = 0.5$. A Backtracking built with a p close to 1 needs more trials in principle, though in the case that matters — a merit frozen at its round-off floor — the safeguarded interpolation shrinks by about a half per trial regardless of p, because the quadratic model through a frozen value has its minimiser at $\alpha/2$.

Quadratic and BierlaireQuadratic are bounded by the same field even though they fit a quadratic rather than shrink a step: they converge on their own ε tolerance long before the budget, which serves only as a backstop, so there is no reason for them to carry a separate knob.

Compare this to default_tolerance and absolute_tolerance.

Examples

julia> linesearch_iterations(Float64)
60

julia> linesearch_iterations(Float32)
31

linesearch_max_iterations(o::Options)

The maximum number of trial steps a line search may take within a single solver step. See linesearch_iterations; not to be confused with max_iterations, which bounds the outer nonlinear iteration.

linesearch_problem(nl::NonlinearSolver)

Build a line search problem based on a NonlinearSolver.

Examples

We show how to set up the LinesearchProblem for a simple example and compute $f^\mathrm{ls}(\alpha_0)$ and $\partial{}f^\mathrm{ls}/\partial\alpha(\alpha_0)$.

julia> F(y, x, params) = y .= (x .- 1.).^2;

julia> x = ones(3)/2; y = similar(x); nl = NewtonSolver(x, y; F = F);

julia> _params = NullParameters();

julia> direction!(nl, x, _params, 1)
3-element Vector{Float64}:
 0.25
 0.25
 0.25

julia> ls_prob = linesearch_problem(nl);

julia> state = NonlinearSolverState(x); update!(state, x, F(y, x, _params));

julia> params = (parameters = _params, x = state.x)
(parameters = NullParameters(), x = [0.5, 0.5, 0.5])

julia> ls_prob.F(0., params)
0.1875

julia> ls_prob.D(0., params)
-0.375

linesearch_problem(nlp, jacobian, cache)

Make a line search problem for a Newton solver (the cache here is an instance of NonlinearSolverCache).

Extended help

The line search closures evaluate the merit at trial steps α using private scratch buffers rather than the solver's shared cache. The shared buffers (solution/value/jacobianmatrix) are read by the solver after the line search returns (e.g. the next direction! step and the convergence check), so writing trial iterates into them would be an aliasing hazard. The line search therefore only reads the current direction from the shared cache and the current iterate from params.x; every write goes to a closure-owned buffer.

params may carry an optional φ₀ field holding the merit at the $\alpha = 0$ anchor. Every line search evaluates that anchor first, and it is exactly the residual the solver has already computed at the current iterate, so solver_step! passes it along and saves one F evaluation per solver step — the most expensive single operation for a large residual. A caller who drives solver_step! by hand from a state whose value is stale must not supply it.

linesearch_warnings(status, ls, params=NullParameters())

Report a LinesearchStatus obtained from solve_with_status. Compare this to nonlinear_solver_warnings. This is the only place where a line search emits log messages, so solve and solver_step! report identically.

Two things keep this quiet in normal use. LINESEARCH_FLOOR and LINESEARCH_STATIONARY are reported only at verbosity ≥ 2, because both are the expected final state of a converged solve — a residual that cannot be improved because it is already as small as the arithmetic allows. And the remaining outcomes are rate limited with maxlog, because a solve that cannot make progress asks the line search for an impossible decrease at every one of its iterations, which an unconditional warning turns into thousands of identical messages.

Whether an irreducible merit actually matters is the outer iteration's call, and nonlinear_solver_warnings makes it: it reports stagnation once, naming the residual that was achieved and the tolerance that was requested.

The messages themselves live in report_linesearch_status rather than here, which is a compile-time rather than a stylistic decision — see its docstring before merging them back.

lucache_eltype(T)

The element type used by the LUSolverCache for an input matrix of element type T. Linear solves are only supported for floating-point problems — real (AbstractFloat) or complex (Complex{<:AbstractFloat}) — so any other element type (e.g. an integer or rational matrix) is rejected here with a clear error rather than silently promoted. For a supported type the cache uses T unchanged.

max_stalls(o::Options)

The number of consecutive stalled steps after which a NonlinearSolver gives up. See MAX_STALLS and stalled_step.

maybe_refactorize!(s, x, params, iteration; force=false, stalled=false)

Re-evaluate the Jacobian at x, copy it into the LinearProblem (adding the diagonal regularization_factor), and refactorize the LinearSolver — but only on a refactorization step: a fresh state or the first step (iteration ≤ 1), every refactorize iterations (see Newton), when the previous step made no progress (stalled, see needs_refresh), or when forced (used by the DogLegSolver to recover from a collapsed trust-region radius). Otherwise the stale Jacobian and its factorization are reused (quasi-Newton). Returns the solver s.

stalled is what makes the quasi-Newton mode safe to combine with max_stalls: a step that did not move the iterate would otherwise rebuild the same direction from the same stale Jacobian on the next refactorize - 1 iterations and reproduce the same negligible step, so the solve would be given up on (see stalled_step) for a reason a fresh Jacobian could have fixed. Refreshing immediately means the second consecutive stall is one that a fresh Jacobian did not fix, which is the conclusive evidence max_stalls = 2 assumes. It is also the response check_anchor prescribes for an ascent anchor, which is a stale-Jacobian symptom.

meets_stopping_criteria(state, config)

Determines whether the iteration stops based on the current NonlinearSolverState.

The function meets_stopping_criteria returns true if one of the following is satisfied:

• the status::NonlinearSolverStatus is converged (checked with isconverged) and state.iterations ≥ config.min_iterations,
• the status has stagnated (checked with isstalled, i.e. config.max_stalls consecutive steps that did not move the iterate while the residual is not small) and state.iterations ≥ config.min_iterations,
• status.f_increased and config.allow_f_increases = false (i.e. f increased even though we do not allow it),
• state.iterations ≥ config.max_iterations,
• status.rfₐ > config.f_abstol_break (by default Inf). In theory this returns true if the residual gets too big.
• one of the residuals (rxₛ, rfₐ, rfₛ) is NaN (checked with havenan) and state.iterations ≥ 1,

So convergence is only one possible criterion for which meets_stopping_criteria. We may also satisfy a stopping criterion without having convergence!

Examples

In the following example we show that meets_stopping_criteria evaluates to true when used on a freshly allocated NonlinearSolverStatus:

julia> config = Options(verbosity=0);

julia> x = [NaN, 2., 3.]
3-element Vector{Float64}:
 NaN
   2.0
   3.0

julia> f = [NaN, 10., 20.]
3-element Vector{Float64}:
 NaN
  10.0
  20.0

julia> cache = NonlinearSolverCache(x, copy(x));

julia> state = NonlinearSolverState(x);

julia> update!(state, x, f); state.iterations += 1
1

julia> status = NonlinearSolverStatus(state, config);

julia> meets_stopping_criteria(state, config)
true

This obviously has not converged. To check convergence we can use assess_convergence. ```

method(ls)

Return the method (of type LinearSolverMethod) of the LinearSolver.

minimum_decrease_threshold(T)

The minimum value by which a function $f$ should decrease during an iteration.

The default value of $10^{-4}$ is often used in the literature [4], [1].

Examples

julia> minimum_decrease_threshold(Float64)
0.0001

julia> minimum_decrease_threshold(Float32)
0.0001f0

nan_recovery!(s, x, params)

Damp direction(cache(s)) by nan_factor until the trial iterate x + d has a finite residual (or the nan_max_iterations budget is exhausted). On return solution(cache(s)) and value(cache(s)) hold the last trial iterate and its residual. Used by the generic and Picard solver_step!s. Returns the solver s.

needs_refresh(state)

true when the previous step made no progress, i.e. when a stall has been flagged for the current step (flag_stall!) or the consecutive-stall counter is nonzero (stall_number).

solver_step! passes this to maybe_refactorize! as its stalled keyword, so a quasi-Newton solver rebuilds its Jacobian immediately after a step that did not move the iterate instead of waiting for the next refactorize multiple. Both sources are consulted because record_stall! consumes the flag into the counter once per iteration, and a caller who drives solver_step! by hand may never call it.

nonlinear_solver_warnings(status, config)

Report a NonlinearSolverStatus at the end of a solve!: the iteration count if it reached warn_iterations, stagnation at the residual floor (see isstalled and stalled_step), a disallowed residual increase, a residual beyond f_abstol_break, and NaNs. Compare this to linesearch_warnings, which does the same for the inner line search, and to print_status.

All messages except the iteration count and the two hard-failure ones are gated on config.verbosity ≥ 1.

nonlinearproblem(solver)

Return the NonlinearProblem contained in the NonlinearSolver. Compare this to linearsolver.

outcome(status)

The LinesearchOutcome stored in status::LinesearchStatus.

pivot_index(v, k)

Return the index (starting from k) of the entry of v with the largest absolute value.

This is used for pivoting in factorize!.

print_jacobian([io], J)

Display the Jacobian J as an aligned text/plain table.

Output is written to io (defaulting to stdout).

print_status(status, config)

Print the solver status if:

• config.verbosity $\geq1$ and one of the following three

1. the solver is converged,
2. status.iterations ≥ config.max_iterations,
3. status.iterations ≥ config.warn_iterations

• config.verbosity $>1.$

record_stall!(state, config)

Update the consecutive-stall counter of state::NonlinearSolverState: increment it when the last step stalled_step or the line search flagged a stall (see flag_stall!), and reset it to zero otherwise. The flag is cleared either way. Returns the new count (see stall_number).

This must be called exactly once per iteration — solve! does so right after update!. That is why the counter is not maintained inside assess_convergence or NonlinearSolverStatus: those are pure and are evaluated more than once per iteration, so incrementing there would double-count. A hand-rolled iteration that drives solver_step! directly and never calls record_stall! simply keeps the count at zero and behaves exactly as before.

report_linesearch_status(status, name, config)

Emit the messages for a LinesearchStatus; the reporting half of linesearch_warnings, whose docstring documents the verbosity and maxlog policy.

Implementation

This is a function barrier, and its signature is what makes it one. linesearch_warnings is called from solver_step! on every iteration of every solve, and takes a Linesearch — which carries the closure types of its LinesearchProblem — and a NamedTuple of parameters, so it is specialized once per problem a solver is built for. A message in its body is specialized with it, and all of the Base.CoreLogging and string-interpolation code that @warn expands to is re-inferred and re-codegen'd for each one, which on a caller that builds one solver per tableau dominates the cost of the whole solve.

Taking name and config, and nothing whose type can vary per solver, bounds the specializations of this function to one per element-type combination for the whole session. nonlinear_solver_warnings and print_status have the same shape for the same reason.

So: do not give this function a parameter whose type varies per solver, and do not move the messages back into linesearch_warnings. test/linesearch_tests.jl asserts both — the first from the method signature, which bounds the specialization set rather than sampling it, and the second by scanning the lowered code of each function for Base.CoreLogging.

The @noinline is a guard rather than the mechanism: Julia's inliner refuses a body this size anyway, but a future one that is more willing would undo the barrier, and nothing in the caller wants this inlined.

The element types are deliberately not tied together as LinesearchStatus{T}/Options{T}: this is a reporting path, and a precision mismatch anywhere upstream should not turn a diagnostic into a MethodError that replaces the problem being diagnosed. nonlinear_solver_warnings is written the same way.

residual_small(rfₐ, config, state)

Return true when the absolute residual rfₐ passes the standard $\mathrm{atol} + \mathrm{rtol}\cdot\|F(x_0)\|$ residual test,

\[r^f_a \leq \texttt{f\_abstol} + \texttt{f\_reltol}\cdot\|F(x_0)\|,\]

with $\|F(x_0)\|$ the initial_residual of state. This lets a large-magnitude or ill-conditioned solve converge once its residual is reduced by f_reltol from $\|F(x_0)\|$, while a step that stalls near $\|F(x_0)\|$ still fails. The relative term drops to zero until the state has been initialized (initial_residual is NaN), leaving the pure absolute f_abstol test.

This gate is shared by assess_convergence, which requires it in addition to a successive-change criterion, and by stalled_step, which requires its negation: a frozen iterate is convergence when the residual is small and stagnation when it is not. The two are therefore mutually exclusive by construction.

residuals(state)

Compute the residuals for state::NonlinearSolverState. The computed residuals are the following:

• rxₛ : successive residual (the norm of $x - \bar{x}$),
• rfₐ: absolute residual in $f$,
• rfₛ : successive residual (the norm of $y - \bar{y}$).

resolve_jacobian(F, DF!, jacobian, x, y)

Resolve the Jacobian for a nonlinear-solver constructor: an explicit DF! wins (wrapped as a JacobianFunction), otherwise an explicit jacobian, otherwise a lazily-built JacobianAutodiff. Building the autodiff Jacobian lazily avoids allocating a ForwardDiff config when either DF! or a jacobian is supplied.

rhs(cache)

Return the right-hand side of the equation, stored in cache::NonlinearSolverCache.

shift_χ_to_avoid_stalling(χ, a, b, c, ε)

Check whether b is closer to a or c and shift χ accordingly.

This is taken from [4].

solve!(x, s, state)

Solve the NonlinearProblem contained in the NonlinearSolver with the initial condition x.

You also have to supply a NonlinearSolverState.

solve!(x, ls::LinearSolver, A, b)

Solve the linear system described by:

\[ Ax = b,\]

and store it in x. Here $A$ and $b$ are provided as an input arguments.

Compare this to solve(::LinearSolver, ::AbstractVector).

solve!(x, ls::LinearSolver, b)

Solve the linear system described by:

\[ Ax = b,\]

and store it in x. Here $b$ is provided as an input argument and the factorized $A$ is stored in the LinearSolver ls (respectively its LinearSolverCache).

solve!(x, ls::LinearSolver, lsys::LinearProblem)

Solve the LinearProblem lsys with the LinearSolver ls and store the result in x.

Also see solve(::LU, ::AbstractMatrix, ::AbstractVector).

Examples

julia> x = zeros(3)
3-element Vector{Float64}:
 0.0
 0.0
 0.0

julia> A = [1.; 0.; 0.;; 0.; 2.; 0.;; 0.; 0.; 4.]
3×3 Matrix{Float64}:
 1.0  0.0  0.0
 0.0  2.0  0.0
 0.0  0.0  4.0

julia> b = ones(3)
3-element Vector{Float64}:
 1.0
 1.0
 1.0

julia> ls = LinearSolver(LU(), x);

julia> problem = LinearProblem(x); update!(problem, A, b);

julia> solve!(x, ls, problem)
3-element Vector{Float64}:
 1.0
 0.5
 0.25

solve!(ls::LinearSolver, args...)

Solve the LinearProblem with the LinearSolver ls.

solve(lu, A, b)

Solve the linear problem determined by A and b.

This is the most straightforward way to solve this system.

Examples

julia> A = [1.; 0.; 0.;; 0.; 2.; 0.;; 0.; 0.; 4.]
3×3 Matrix{Float64}:
 1.0  0.0  0.0
 0.0  2.0  0.0
 0.0  0.0  4.0

julia> b = ones(3)
3-element Vector{Float64}:
 1.0
 1.0
 1.0

julia> solve(LU(), A, b)
3-element StaticArraysCore.SizedVector{3, Float64, Vector{Float64}} with indices SOneTo(3):
 1.0
 0.5
 0.25

Compare this to solve!(::AbstractVector, ::LinearSolver, ::LinearProblem).

solve(ls::LinearSolver, args...)

Counterpart of solve! for a prebuilt LinearSolver: allocates (and returns) a fresh solution vector instead of writing into a caller-supplied one. Note that the solver's cache is still updated in place (the factorization is computed there).

Accepts the same trailing arguments as solve!(ls, args...): a LinearProblem, a matrix-vector pair A, b, or a bare right-hand side b (the latter uses the factorization already stored in ls).

solve(linesearch, α, params=NullParameters())

Solve the LinesearchProblem (contained in Linesearch) starting at α.

The argument params needs to be of an appropriate form expected by the respective LinesearchProblem.

See linesearch_problem.

solve(ls::Linesearch{T,<:Backtracking}, α, params)

Run the backtracking line search from the trial step α, report the outcome through linesearch_warnings and return the accepted step length.

Use solve_with_status to obtain the LinesearchStatus instead: a caller that has to tell "I found a decreasing step" from "the merit is at its round-off floor and nothing can decrease it" cannot do so from the step length alone.

solve(ls::Linesearch{T,<:BierlaireQuadratic}, α, params)

Fit successive quadratics through three bracketing points to approximate the line minimiser, report the outcome through linesearch_warnings and return the step length. See BierlaireQuadratic and solve_with_status.

solve(ls::Linesearch{T,<:Bisection}, α, params)

Bisect the derivative of the merit to approximate the line minimiser, report the outcome through linesearch_warnings and return the step length. See Bisection and solve_with_status.

solve(ls::Linesearch{T,<:Quadratic}, α, params)

Fit successive quadratics to approximate the line minimiser, report the outcome through linesearch_warnings and return the step length. See Quadratic and solve_with_status.

solve(ls::Linesearch{T,<:StrongWolfe}, α, params)

Run the strong-Wolfe bracketing line search starting from the trial step α, report the outcome through linesearch_warnings and return the accepted step length. See StrongWolfe and solve_with_status.

solve_with_status(ls, α, params=NullParameters())

Like solve, but return a LinesearchStatus — the step length plus the reason the search stopped — and emit no log messages. Use linesearch_warnings to report the status.

Only Backtracking reports a genuine LinesearchOutcome; the generic fallback calls solve and reports LINESEARCH_UNKNOWN, so a caller may use solve_with_status uniformly for every LinesearchMethod.

solver_step!(x, s, state, params)

Compute one step for solving the problem stored in an instance s of NonlinearSolver.

Examples

julia> f(y, x, params) = y .= sin.(x .- .5) .^ 2
f (generic function with 1 method)

julia> x = ones(3) / 4
3-element Vector{Float64}:
 0.25
 0.25
 0.25

julia> y = zero(x)
3-element Vector{Float64}:
 0.0
 0.0
 0.0

julia> s = NewtonSolver(x, similar(x); F = f);

julia> state = NonlinearSolverState(x); update!(state, x, f(y, x, NullParameters()));

julia> solver_step!(x, s, state, NullParameters())
3-element Vector{Float64}:
 0.37767096061051814
 0.37767096061051814
 0.37767096061051814

solver_step!(x, s::PicardSolver, state, params)

Take one fixed-point (Picard) step $x \gets x + \alpha d$ with the residual direction $d = -F(x)$ (see direction!).

Unlike a Newton/Gauss-Newton step, the Picard direction $d = -F(x)$ is not in general a descent direction for the merit $\varphi = \|F\|^2$, so applying the derivative-based (Wolfe) line search used by the other NonlinearSolvers is inappropriate (a directional derivative that is not negative makes the sufficient- decrease/curvature tests meaningless).

Instead the step is damped by a residual-monotonicity backtracking: starting from the full fixed-point step $\alpha = 1$ the step is halved until the residual norm does not increase, $\|F(x + \alpha d)\| \le \|F(x)\|$. This safeguard uses only function values and makes no descent assumption. If no positive $\alpha$ reduces the residual (the fixed-point map is locally expanding), the smallest trial step is taken and the convergence test — which requires a small residual, not merely a small step — correctly reports non-convergence instead of a false positive.

stall_number(state)

Return the number of consecutive stalled steps recorded in state::NonlinearSolverState by record_stall!. Compare this to iteration_number.

stalled_step(rxₛ, rfₐ, config, state)

Return true when the last step stalled: it left the iterate unchanged (see iterate_settled) while the residual is not small (see residual_small).

A stalled step is the failure mode that the residual gate in assess_convergence correctly refuses to call convergence, and that used to be invisible to the solver. The step length $\alpha\|d\|$ has dropped below the round-off level of $x$, so the merit $\|F\|^2$ cannot be reduced along the current direction — typically because the requested f_abstol lies below the round-off floor of the residual itself. Taking another step recomputes the same direction and the same negligible $\alpha$, so the iteration would spin all the way to max_iterations, asking the line search on every one of those steps to improve a residual that is already pure round-off noise.

meets_stopping_criteria therefore stops after config.max_stalls consecutive stalled steps (counted by record_stall!), and nonlinear_solver_warnings reports the achieved residual against the requested tolerance instead of the misleading "Solver took 1000 iterations.".

status(solver, state)

Return the NonlinearSolverStatus for the NonlinearSolverState state as assessed with the Options of solver.

solve! returns the solution x (updated in place), not a status, so this is how a caller inspects the outcome of a solve — in particular whether it converged (isconverged) or merely stagnated at the residual floor (isstalled). The state is the caller's own object (it is passed to solve!), so nothing has to be threaded back out of the solve.

Examples

julia> F(y, x, params) = y .= x .^ 2 .- 2;

julia> x = [1.0]; s = NewtonSolver(x, similar(x); F = F, verbosity = 0);

julia> state = SolverState(s);

julia> solve!(x, s, state);

julia> isconverged(status(s, state))
true

steplength(status)

The step length stored in status::LinesearchStatus, i.e. what solve returns.

trials(status)

The number of trial steps at which the merit was evaluated, stored in status::LinesearchStatus.

triple_point_finder(f, x)

Find three points a < b < c (strictly ordered in position) with f(a) ≥ f(b) and f(c) > f(b), so that a minimum is bracketed in (a, c). This is used for performing a quadratic line search (see BierlaireQuadratic). Returns a Symbol instead of a triple when no such triple exists — see the warning below.

Implementation

For δ we take DEFAULT_BRACKETING_s as default. For nmax we take DEFAULT_BRACKETING_nmax as default.

Examples

julia> f(x) = x ^ 2
f (generic function with 1 method)

julia> x = -1.
-1.0

julia> a, b, c = round10.(triple_point_finder(f, x))
(-0.37, 0.27, 1.55)

julia> round10.((f(a), f(b), f(c)))
(0.1369, 0.0729, 2.4025)

Extended help

The algorithm is taken from [4, Chapter 11.2.1].

trust_radius!(cache::DogLegCache, Δ)

Store the trust-region radius $\Delta$ in the DogLegCache so it carries over to the next outer solver step.

trust_radius(cache::DogLegCache)

Return the current trust-region radius $\Delta$ carried by the DogLegCache.

value!(y, nlp, x, params)

Evaluate the NonlinearProblem at x.

with_config(ls, config)

Return a Linesearch with the LinesearchProblem and the LinesearchMethod of ls, but with the Options config.

This is how a NonlinearSolver makes its line search share its options. A Linesearch built by Linesearch(problem, method) carries an Options of its own, constructed from nothing but defaults — so verbosity and linesearch_max_iterations would be configured twice, and verbosity = 0 on the solver would not silence the line search.

Linesearch is an immutable three-field wrapper, so rebuilding it is cheap: the problem (and hence its closures and scratch buffers) and the method are shared, not copied.

The Options element type is pinned to the Linesearch element type, so a mismatched config raises a MethodError rather than silently producing a broken object — the same guarantee the three-argument Linesearch constructor gives.