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_GROW_MIN

Lower bound on the factor by which the expansion phase of Backtracking lengthens an accepted step: unless the model minimiser lies at least $\mathrm{BACKTRACKING\_GROW\_MIN}\cdot\alpha$, the trial step is kept and no further merit evaluation is spent (see backtracking_extrapolation). Its value is 2.0.

This is the counterpart of BACKTRACKING_SHRINK_MIN on the growing side, and it is what makes the expansion phase free for a well-scaled direction: a Newton or BFGS step is already at its model minimum, so the test fails and the search returns at once.

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_τ_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_BACKTRACKING_NEXPAND

The default cap on the number of expansion trials of Backtracking; each one costs a merit evaluation. Its value is 3.

const DEFAULT_BACKTRACKING_q

The default upper bound on the factor by which the expansion phase of Backtracking lengthens the step in one round, i.e. the counterpart of $p$ on the growing side. Its value is 10.0.

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.

const DEFAULT_LINESEARCH_αmax

The largest step length a LinesearchMethod will try unless a caller asks for less; its value is 65536.0 ($2^{16}$). See linesearch_αmax for the two ways the ceiling is set and method_αmax for the per-method field it defaults.

An absolute number is the right shape here, and this is the one place in the package where that is worth arguing. $\alpha$ scales a direction that has already been chosen, so it is a step-length fraction and of order one — the same reason the bracket-width tolerance of BierlaireQuadratic is absolute while its merit comparisons are not. A ceiling of $2^{16}$ is therefore four to five orders of magnitude of headroom above any step a well-scaled direction asks for, and still rules out the $\alpha \approx 4\cdot10^7$ that an unbounded bracketing search can reach on a merit whose minimiser is far away or whose fit is nearly flat.

The value is StrongWolfe's, which has carried exactly this field since before the other methods had one; DEFAULT_WOLFE_αmax is now defined as this constant so the two cannot drift apart.

Use default_linesearch_αmax rather than T(DEFAULT_LINESEARCH_αmax) to obtain it in a given precision: $2^{16}$ is above floatmax(Float16), so the plain conversion overflows to Inf and silently removes the ceiling in the precision the package otherwise takes most care over.

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. This is DEFAULT_LINESEARCH_αmax, which every method's ceiling now defaults to: the field was StrongWolfe's alone until the bracketing searches were found to extrapolate without one, and defining this as that constant is what keeps the two from drifting apart. The constructor obtains it through default_linesearch_αmax, which saturates it at floatmax(T) rather than overflowing to Inf in Float16.

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 F_STALL_FACTOR

The factor by which the residual has to drop for an iteration to count as progress; the default of the f_stall_factor field of Options. Its value is 0.5, i.e. progress means "the residual halved".

record_progress! keeps the residual $r^f_a$ of the last iteration that counted as progress, and counts the iterations since, so that iterations_since_progress measures how long the residual has been going nowhere. That one number is used twice: unconditionally by nonlinear_solver_warnings to explain a solve that spent its whole budget, and — only when f_stall_window > 0 — as the stopping criterion no_progress.

A factor closer to one (say 0.999, i.e. "any improvement of a tenth of a percent counts") makes the measurement more permissive: a residual creeping down that slowly keeps resetting the clock, so only one that is essentially flat is ever reported. A smaller factor demands a steeper descent of an iteration before it counts as progress, and so reports sooner.

Only $0 < f_{\mathrm{stall\,factor}} \le 1$ is meaningful. Nothing asserts it — Options validates none of its fields — but a factor above one would make a residual that grew count as progress and so let the reference climb, which is exactly the monotonicity record_progress! relies on to be immune to a residual that jumps around.

See F_STALL_WINDOW.

const F_STALL_REPORT_MINIMUM

The fewest iterations without progress that spent_without_progress will report on. Its value is 10.

Unlike f_stall_window this is not configurable, because it is not a policy: it is the point below which the proportion that predicate measures is not evidence of anything. At two iterations the proportion is satisfied by a single one, and a solve that gets nowhere for four iterations has not shown you anything a solve that got somewhere on the fifth would not.

It is the second of the two guards spent_without_progress carries. The first is isconverged — a residual that stopped improving because it was already small enough is success, not stagnation, and that is what excuses the short solve which converged on its successive-change criterion without halving its residual on the last step. This minimum then covers the solve that has not converged and is merely too short for the proportion to mean anything.

This bounds only what is said about a solve, never what is done with it — the stopping criterion is f_stall_window alone (see no_progress).

const F_STALL_WINDOW

The default number of iterations without progress after which a NonlinearSolver gives up; the default of the f_stall_window field of Options. Its value is 0, which disables the criterion (see no_progress).

This is the counterpart of max_stalls for a solve whose iterate keeps moving: the residual sits on a floor far above the requested tolerance while the step is nowhere near the round-off level of $x$, so neither stalled_step nor either convergence branch can fire and the iteration spends max_iterations in full. That floor is typically set by the problem itself — a model, discretisation or ansatz error that $F$ cannot resolve — rather than by round-off, so no $\mathrm{eps}(T)$-scaled tolerance can bound it.

It is off by default because the threshold is a policy, not a test, and a wrong one gives up on a healthy solve. An iteration converging linearly with rate $\rho$ improves by $\rho^W$ over a window $W$: at f_stall_factor = 0.5 and f_stall_window = 50 every $\rho > 2^{-1/50} \approx 0.986$ is abandoned, and a PicardSolver on a stiff problem is slower than that. There is no value that is right for every problem, which is why the diagnosis is unconditional (a solve that spent its budget without progressing says so, see nonlinear_solver_warnings) and only the stopping is opt-in: the report costs a solve that has already failed nothing, whereas the criterion can cost a solve that would have succeeded everything.

So set it once you have seen the report and know the floor is real — a window of a few tens of iterations then turns a full budget into a prompt answer.

HAS_PREALLOCATED_GETRF

Whether this Julia's LinearAlgebra.LAPACK.getrf! accepts a pre-allocated pivot vector.

getrf!(A, ipiv) arrived after the 1.10 LTS, and it is the whole reason factorize! can be allocation-free for LapackLU. Where it is missing, the one-argument getrf!(A) is used and its pivot vector copied into the cache, which costs one O(n) allocation per factorization — 3.3 kB at n = 384. The O(n^2) working matrix is reused on every version, and ldiv! is allocation-free on every version, because getrs! has always taken the pivot vector as an argument.

Since LapackLU is the default (see default_linear_solver_method), on the LTS that allocation lands once per factorization in a nonlinear solve. Closing it would mean calling getrf through a hand-written ccall to pass a pre-allocated pivot vector — the only place in the package reaching past the stdlib — for a benefit confined to a Julia version with a finite lifetime. Not worth the maintenance; when the compat floor reaches 1.11 this constant and its else branch delete cleanly.

This is a feature check rather than a VERSION comparison, so the exact release it arrived in does not have to be tracked here.

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.

const NLINESEARCH_OUTCOMES

The number of LinesearchOutcomes, i.e. the length of the tally a NonlinearSolverState keeps of the outcomes its line search reported (see record_linesearch!). A const rather than length(instances(LinesearchOutcome)) at each use, because it is a type parameter of that tally's MVector and so has to be known to the compiler.

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

const WARNING_COUNTS

How often each diagnosis has been reported so far, keyed by the Symbol passed to should_report!. Process-global, and guarded by WARNING_LOCK.

The number of keys is bounded by the number of message sites, so this does not grow without limit. It is reset by reset_warning_counts!.

Guards WARNING_COUNTS; see should_report!.

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 trial step $\alpha$ is not a key: it is the argument of solve, which is its only source. (A Backtracking used to carry an α₀ field for it, which the algorithm never read — see issue #174.) Neither is αmax: this search shrinks, so the trial step is already its ceiling and the expansion phase carries its own in $q^{\mathrm{nexpand}}$. A params.αmax still binds — on the trial step itself as well as on the expansion — since a caller that says no step above a given length is admissible means the first trial too. See SimpleSolvers.linesearch_αmax.

The keys are:

• 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.
• expand=false: whether the search may lengthen the trial step (the expansion phase described below). Off by default, so a Backtracking is the classical one-sided algorithm unless it is asked for. Setting it requires the merit to be evaluable — finite or not, but not throwing — out to $q^{\mathrm{nexpand}}\alpha$, since that is the largest step the phase can try.
• q=10.0: an upper bound on the factor by which $\alpha$ is increased in one expansion round — the counterpart of $p$. Only used when expand is set. See DEFAULT_BACKTRACKING_q.
• nexpand=3: the cap on the number of expansion trials, each of which costs one merit evaluation. It bounds the phase from within the linesearch_max_iterations of Options rather than beside it: whichever of the two is smaller applies, so the whole search still spends at most linesearch_max_iterations merit evaluations. Only used when expand is set. See DEFAULT_BACKTRACKING_NEXPAND.
• τ_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). With expand set it is satisfied more often, because the expansion phase moves the accepted step toward the line minimiser, so the diagnostic fires less.

The expansion phase

With expand = true and the first trial step accepted, the search may also lengthen it, by backtracking_extrapolation, while each longer trial still satisfies the SufficientDecreaseCondition and strictly improves the merit; at most nexpand such trials are made and the best step seen is returned. A shrunken step is never expanded again: once the ladder has backtracked, the longer steps are already known to fail.

This is the one place where the search leaves the interval $[0, \alpha]$ the caller offered. The largest step it can try is $q^{\mathrm{nexpand}}\alpha$ — a thousand times the trial step on the defaults — and a trial whose merit is not finite is simply rejected, at the cost of the one evaluation. A merit that throws outside its domain is the caller's to guard, and it is one reason the phase is opt-in.

The trials it spends come out of the linesearch_max_iterations budget of Options, not out of a second budget beside it, so termination case 4 above still bounds the whole search: the phase makes at most nexpand trials and at most as many as that budget has left.

This is what makes the search two-sided. A shrink-only search returns the trial step it was given whenever that step is acceptable, so on a direction whose natural scale is larger than the trial step it pins $\alpha$ at that ceiling on every iteration and the outer solve crawls — by two orders of magnitude in the DFP case of issue #174, where the direction wanted $\alpha \approx 11$ throughout.

The phase costs nothing where it can gain nothing. The model it extrapolates from is the same quadratic that backtracking_interpolation uses on the way down, built from $\varphi(0)$, $\varphi'(0)$ and the merit at the trial step — all three already known — so the decision whether to expand at all is free, and a direction that is already scaled like a Newton step (which is at its model minimum at $\alpha = 1$) fails the test and returns without a further merit evaluation. That matters because for the merit of a NonlinearSolver an evaluation is a full residual evaluation, the most expensive single operation of a solver step, which is also why the phase does not test the CurvatureCondition to decide when to stop growing: that would cost a full Jacobian per trial. Use StrongWolfe, which brackets on the derivative, where curvature control is genuinely required.

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].

Keywords

• ε: the bracket-width tolerance of the fit.
• ξ: the threshold below which $|\varphi'(\alpha_0)|$ counts as stationary.
• αmax: the largest step the triple-point bracketing will try, by default DEFAULT_LINESEARCH_αmax. Without it the bracketing doubles its increment until the merit stops falling, which for a nearly flat or distantly-minimised $\varphi$ is arbitrarily far; see linesearch_αmax, which is also how a caller imposes a smaller ceiling of its own.

Bisection <: LinesearchMethod

See bisection for the implementation of the algorithm.

Keywords

• αmax: the largest step the bracketing will try, by default DEFAULT_LINESEARCH_αmax. See linesearch_αmax, which is also how a caller imposes a smaller ceiling of its own.

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.

What it converges to is a minimum, not a maximum

Bisection drives on the sign of $\varphi'$, so it converges to whichever crossing the sign at its left endpoint selects — and that sign is invariant under the halving. From $\varphi'(\mathrm{lo}) < 0$ the interval shrinks onto a $-\to+$ crossing, a minimum; from $\varphi'(\mathrm{lo}) > 0$ onto a $+\to-$ crossing, a maximum.

Nothing in bracket_minimum rules the second out. It brackets a minimum in value, sampling $\varphi$ and never $\varphi'$, so on a non-convex ray its interval can enclose several stationary points and its left endpoint can sit past one of them. Landing on a maximum used not to be caught: the step was classified by the merit like any other step, so it was LINESEARCH_DECREASED when it happened to improve $\varphi$ and LINESEARCH_FLOOR when it did not — the same overclaim the previous section describes, reached by a different route, and GeometricOptimizers observed exactly it (LINESEARCH_FLOOR reported with $\varphi(1) = \varphi(0)$, and the step taken regardless).

The orientation is now checked and repaired: when $\varphi'(\mathrm{lo}) > 0$ the search bisects $[0, \mathrm{lo}]$ instead, which brackets a minimum by construction — the anchor check has established $\varphi'(0) < 0$ — and finds the earlier one, nearer the anchor. The check costs nothing on a ray that does not need it, because it reads a value the bisection computes anyway. See SimpleSolvers._bisect_for_minimum.

A bracket that fails is never a floor

Bisection drives on the sign of $\varphi'$, so it can only work on an interval whose endpoints straddle a sign change. When bracket_minimum hands it an interval that brackets a minimum in value but over which $\varphi'$ keeps its sign — a non-smooth or noisy merit, or a derivative inconsistent with it — there is nothing to bisect and _bisection_core reports BISECTION_NOBRACKET.

That case used to be folded into "converged", so the endpoint with the smallest $|\varphi'|$ was claimed as the line minimiser and, when it did not improve the merit, classified as LINESEARCH_FLOOR — which asserts that no line search can make progress along this direction and makes the outer iteration count the step towards max_stalls. A failed bracket establishes nothing of the kind. So the outcome is now classified by the merit alone: LINESEARCH_DECREASED when the returned step still beats $\varphi(0)$ by more than $\tau$, and LINESEARCH_EXHAUSTED when it does not. LINESEARCH_FLOOR is reachable only from a bisection that actually converged.

BisectionOutcome

Why _bisection_core stopped. This is what lets a caller tell "found the root" from "gave up" — a distinction a Bool cannot carry, and whose absence made an unbracketable derivative look like a located line minimiser (see Bisection).

• BISECTION_CONVERGED: a root was located, either because $|f(\alpha)| \leq$ f_abstol or because the bracket collapsed to x_suctol.
• BISECTION_NOBRACKET: the endpoint values share a sign, so the interval contains no root of odd multiplicity and bisection cannot start. The endpoint with the smallest $|f|$ is returned, but it is not a root — this is a failure to report, not a result.
• BISECTION_EXHAUSTED: the linesearch_max_iterations budget of Options was spent with the interval still straddling a sign change. Unlike BISECTION_NOBRACKET there is a root in the interval, so a larger budget would find it; the best estimate so far is returned.

This is internal: it is the return type of a private function, like the Symbol that _triple_point_core reports. Callers of a line search see a LinesearchOutcome instead.

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).

DogLegSolver(x, nlp::NonlinearProblem, y = zero(x))

Build a DogLegSolver for the NonlinearProblem nlp with the initial guess x. See NewtonSolver(::AbstractVector{T}, ::NonlinearProblem, ::AbstractVector{T}) where {T} for the rôle of y; as above, no linesearch keyword is accepted — the step length comes from the trust-region radius.

Keywords

• linear_solver_method
• jacobian: see resolve_jacobian,
• jacobian_prototype: see the NewtonSolver counterpart. Note that SparspakLU reports singularity too late for the DogLeg fallback; see its docstring,
• refactorize
• options_kwargs: see Options.

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.

struct LapackLU <: PivotedLUMethod

A LAPACK-backed LU solver, meant to solve a LinearProblem.

Where LU is a self-contained implementation that works for any number type and for static matrices, this method delegates the factorization to LAPACK's getrf. It is restricted to the element types LAPACK handles — Float32, Float64, ComplexF32 and ComplexF64 — and throws an ArgumentError naming the type when handed anything else. Any AbstractMatrix storage is accepted, but the cache always holds a plain Matrix, because that is what LAPACK can be pointed at.

Routines that use it are the same as for LU: factorize!, ldiv! and solve!, all shared with the other PivotedLUMethods through PivotedLUCache.

Constructor

LapackLU()

# output

LapackLU()

When to use which

LU is the better choice for very small systems, where its static-matrix cache avoids allocation altogether, and it is the only choice for element types LAPACK does not know about. LapackLU is the better choice everywhere else, and is the default for a dense matrix of a LAPACK element type — see default_linear_solver_method. Measured on an Apple M4 Max against OpenBLAS, in microseconds for factorize! including the copy-in:

LapackLU's triangular solve is the faster one too, by 3.5–4.5× across that whole range: 21.9 µs against LU's 77.3 µs at n = 384. See RecursiveLU for where the crossover in the factorization sits and why it depends on the BLAS in use.

The effect on a real problem is large. Measured from PoissonBrackets.jl, where a Newton step factorizes a dense $384 \times 384$ Jacobian, LU accounted for 74 % of the cost of one implicit time step — about 17 ms against 0.6 ms for the same factorization through LAPACK.

Allocation is not part of the trade-off. Like LU, and unlike a bare LinearAlgebra.lu!, both factorize! and ldiv! are allocation-free after the LinearSolver is built: the working matrix and the pivot vector are allocated once and reused. See PivotedLUCache.

The one exception is factorize! on a Julia too old for LAPACK.getrf!(A, ipiv) — the 1.10 LTS — where the pivot vector costs one O(n) allocation per call. The O(n^2) working matrix is reused on every version, and ldiv! is allocation-free on every version. See HAS_PREALLOCATED_GETRF.

Example

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

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

# output

true

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. The one deliberate exception is a caller error rather than a situation arising from the merit: a params.αmax that is not a usable ceiling raises an ArgumentError before any evaluation — see linesearch_αmax.
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; the benign round-off-floor and stationary outcomes at ≥ 2). And it reports there only when the user called it: a program calls solve_with_status, acts on the LinesearchStatus and sees no messages at all. A NonlinearSolver is such a program — see record_linesearch!. This is structural rather than a convention a method has to keep: a method implements solve_with_status, and solve is derived from it once, for all methods, as that call plus the report.
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.
6. It returns $\alpha \leq \alpha_\mathrm{max}$ — see linesearch_αmax for the two ways that ceiling is set and why one of them has to be per call.

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$. Backtracking shrinks it only, and therefore returns it unchanged whenever it is acceptable, unless its expand key is set; StrongWolfe brackets in both directions.
• 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 — Static, which evaluates no merit and so has established nothing, and any third-party method that chooses not to report one.

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: what the search cost, as the number of times it evaluated the problem — not the linesearch_max_iterations budget. That is the merit for every method except Bisection, which drives on the derivative it bisects and brackets on the merit, so its count is of both. For Backtracking and StrongWolfe it is exactly the number of trial steps $\alpha > 0$: every merit evaluation is either the $\alpha = 0$ anchor or a counted trial. For the searches that bracket it includes what the bracketing spent (bracket_minimum, triple_point_finder) — that is where those searches do their work, and on the path where a ceiling binds it is the whole of it, so a count omitting it reported one evaluation, or none at all, for a search of any size. One of the evaluations it then counts is the bracketing's own re-evaluation of the anchor, so for those methods the number is the cost rather than exactly the number of distinct positive steps,
• φ₀, 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. For a method that does not measure them — Static, or a third-party method whose solve_with_status has nothing more to say.

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}, ::NonlinearProblem, ::AbstractVector{T}) where {T} and 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
• jacobian_prototype: the matrix whose storage and, if sparse, whose sparsity pattern are adopted by the Jacobian, the LinearProblem and the LinearSolver cache. It is copied, so the caller's matrix is left alone. A SparseMatrixCSC here is what runs a sparse Jacobian through the solver, and it requires DF! (see checkjacobianprototype) and a pattern that includes the diagonal if regularization_factor is non-zero. It also selects the default linear_solver_method; see default_linear_solver_method,
• 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, nlp::NonlinearProblem, y = zero(x))

Build a NewtonSolver for the NonlinearProblem nlp with the initial guess x, assembling the Jacobian, the LinearProblem, the LinearSolver, the Linesearch and the NonlinearSolverCache.

y is a prototype for the residual $F(x)$: it supplies a size and an element type, and nothing that is computed from it survives (alloc_j turns it into a NaN matrix and the cache stores zero(y)). It is not, however, left alone — JacobianAutodiff keeps it as the buffer ForwardDiff writes the residual into, so a caller-supplied y is overwritten on every Jacobian evaluation. It defaults to zero(x), which is what a square system needs — and every system here is square, since the LinearSolver factorizes the Jacobian.

The Jacobian stored in nlp (if any) takes precedence over autodiff, exactly as the DF! keyword of NewtonSolver(::AbstractVector{T}, ::Callable, ::AbstractVector{T}) where {T} does — see resolve_jacobian.

Keywords

• linear_solver_method
• linesearch
• jacobian
• jacobian_prototype
• refactorize
• options_kwargs: see Options

Examples

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

NewtonSolver(x, nlp) isa NewtonSolver

# output

true

NewtonSolver(x, F, y)

Keywords

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

The Callable F (and the optional DF!) are wrapped in a NonlinearProblem, so this is NewtonSolver(::AbstractVector{T}, ::NonlinearProblem, ::AbstractVector{T}) where {T} with the problem built for the caller.

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, NaN, 0, [0, 0, 0, 0, 0, 0])

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,
• iterations_since_progress: iterations since the residual last dropped by f_stall_factor, see iterations_since_progress,
• 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
• not_progressing::Bool: the iteration is not getting anywhere, see no_progress
• ls_outcomes: how often the line search reported each LinesearchOutcome during the solve, indexed by linesearch_index; see linesearch_outcomes and linesearch_failures

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_stall_factor = 0.5
          f_abstol_break = Inf
       allow_f_increases = true
          min_iterations = 0
          max_iterations = 1000
         warn_iterations = 1000
linesearch_max_iterations = 60
              max_stalls = 2
          f_stall_window = 0
              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

PicardSolver(x, nlp::NonlinearProblem, y = zero(x))

Build a PicardSolver for the NonlinearProblem nlp with the initial guess x. See NewtonSolver(::AbstractVector{T}, ::NonlinearProblem, ::AbstractVector{T}) where {T} for the rôle of y; as above, no linesearch keyword is accepted.

Keywords

• jacobian: see resolve_jacobian,
• options_kwargs: see Options.

PivotedLUCache <: LinearSolverCache

The cache shared by every PivotedLUMethod — LapackLU and RecursiveLU.

Keys

• A: the working copy of the matrix, which the factorization overwrites in place,
• ipiv: the pivot vector the factorization fills,
• info: the index of the first zero pivot, or 0 on success,
• factorized: whether factorize! has run at all, which info == 0 cannot express.

The cache holds the pieces of the factorization rather than a LinearAlgebra.LU object, so that factorize! can be called any number of times without allocating: LU is immutable, so a fresh one would have to be built and boxed into a field on every refactorization, and its ipiv freshly allocated. Both are O(n)-to-O(1) costs against an O(n^3) factorization, but a nonlinear solve in a time-stepping loop refactorizes on every step of every step, and LU — the method these sit beside — allocates nothing at all. (On a Julia without LAPACK.getrf!(A, ipiv) — the 1.10 LTS — LapackLU fills ipiv from a per-call temporary instead, which costs one O(n) allocation per factorization; see HAS_PREALLOCATED_GETRF.)

Use factorization to get a LinearAlgebra.LU view of these pieces when one is actually wanted (for det, say).

The pivot vector is typed Vector{LinearAlgebra.BlasInt} rather than Vector{Int}: that is what getrf fills, and the two are not the same type under a 32-bit-integer BLAS. It is also what RecursiveLU writes into, whose ipiv is any AbstractVector{<:Integer} — so one cache type serves both, and the LAPACK triangular solve can be shared with it.

PivotedLUMethod <: DirectMethod

The methods that compute a partially-pivoted LU factorization in LAPACK's layout, and can therefore share a cache, a triangular solve and everything built on them: LapackLU and RecursiveLU.

They differ only in which kernel computes the factors — see _getrf! — and in which element types they accept. See PivotedLUCache.

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.
• αmax: the largest step the bracketing will try, by default DEFAULT_LINESEARCH_αmax. Without it the bracket grows outward until the merit stops falling, which for a nearly flat or distantly-minimised $\varphi$ is arbitrarily far; see linesearch_αmax, which is also how a caller imposes a smaller ceiling of its own.

Extended help

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

struct RecursiveLU <: PivotedLUMethod

A RecursiveFactorization.jl-backed LU solver, meant to solve a LinearProblem.

Only available once RecursiveFactorization is loaded; the constructor exists either way, and building a LinearSolver with it says what to load if it is not. It is a package extension because that dependency is heavy — LoopVectorization, Polyester, StrideArraysCore, TriangularSolve and VectorizedRNG — and the win is confined to a middle range of sizes.

Restricted to Float32 and Float64, and to Matrix storage. RecursiveFactorization has no complex support, so unlike LapackLU this covers two of the four BLAS element types; the cache constructor says so by name.

Everything except the factorization itself is shared with LapackLU — the cache, the triangular solve, singular_index and every solve! form — because RecursiveFactorization writes LAPACK-layout factors with LAPACK's pivot convention, so getrs applies unchanged. See PivotedLUMethod and PivotedLUCache.

Like LapackLU, factorize! and ldiv! are allocation-free once the LinearSolver exists, on every Julia version — the pivot vector is pre-allocated in the cache and RecursiveFactorization takes it as an argument.

Constructor

RecursiveLU()

When to use it

For a middle range of sizes, where a blocked pure-Julia kernel beats the BLAS's own but the O(n^3) term has not yet taken over. Measured on an Apple M4 Max, factorize! in microseconds including the copy-in:

Below n ≈ 10 the default LU() uses a static-matrix cache and allocates nothing, which is the better trade there. Above the crossover use LapackLU. For element types neither LAPACK nor RecursiveFactorization handles, LU remains the only option.

SparseDirectMethod <: DirectMethod

The direct methods that factorize a sparse matrix, keeping the ordering and symbolic factorization across refactorizations: UmfpackLU and SparspakLU.

Unlike a PivotedLUMethod these need the sparsity pattern up front, at LinearSolver construction, and they refuse a dense matrix. See SparseFactorizationCache.

SparseFactorizationCache <: LinearSolverCache

The cache shared by the SparseDirectMethods, UmfpackLU and SparspakLU.

Keys

• F: the backend's factorization object, which owns the ordering, the symbolic factorization and the numeric factors,
• n: the leading dimension, kept here because the two backends store it differently,
• info: 0 if the last factorization or solve found nothing wrong, non-zero otherwise; see singular_index,
• factorized: whether factorize! has run at all.

Unlike PivotedLUCache there is no working copy of the matrix: both backends take the SparseMatrixCSC as an argument to their refactorize call and read its nzval directly, so a copy here would be dead weight. The consequence is that the pattern is fixed at construction — which is exactly the contract a Newton loop wants, since reusing the ordering and symbolic factorization is where the saving is.

Neither backend is allocation-free, and that is inherent to them rather than to this wrapper. Measured on a periodic banded matrix at n = 384: UmfpackLU allocates ~374 kB per refactorization but 0 B per ldiv!; SparspakLU allocates ~11 kB per refactorization and ~10 kB per ldiv!.

struct SparspakLU <: SparseDirectMethod

A sparse LU solver backed by Sparspak.jl, meant to solve a sparse LinearProblem.

Only available once Sparspak is loaded; the constructor exists either way, and building a LinearSolver with it says what to load if it is not.

Constructor

SparspakLU()

Why this exists alongside UmfpackLU

Element types. Sparspak is generic in the element type where UMFPACK is not. Probed on a periodic banded matrix:

That last row is the point: a sparse solve over ℚ with no rounding at all, which nothing else here can do. It is also a pure-Julia stack, with no SuiteSparse binary. The Float32 row is not a typo either: UMFPACK converts a 32-bit matrix in lu/lu! but has no 32-bit solve, so UmfpackLU refuses those element types at construction rather than failing later inside ldiv!, and this is one of the two methods that cover them.

Not for speed on Float64. Its factorization is in fact the faster of the two, by 1.3–1.5×, but its triangular solve is about 9× slower, and a Newton loop does at least one solve per factorization. Measured at n = 384 on a periodic banded matrix: 59.5 µs factorize + 32.7 µs solve against UMFPACK's 76.2 + 3.5. Prefer UmfpackLU for the element types it handles.

Allocation

Neither factorize! nor ldiv! is allocation-free — about 11 kB and 10 kB respectively at n = 384 — and that is inside Sparspak rather than in this wrapper.

Singularity: reported late

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.

struct UmfpackLU <: SparseDirectMethod

A sparse LU solver backed by SuiteSparse's UMFPACK, meant to solve a sparse LinearProblem.

No extension and no new dependency: UMFPACK ships inside the SparseArrays standard library, which this package depends on for the sparse-Jacobian plumbing anyway. Restricted to the element types UMFPACK provides — Float64 and ComplexF64 — and it names the alternatives for anything else. That includes Float32 and ComplexF32: SuiteSparse converts them in lu and lu! but has no 32-bit solve, so a 32-bit factorization would be a ldiv! MethodError waiting to happen rather than a working narrow path. Use SparspakLU to keep a 32-bit matrix sparse, or LapackLU to densify it.

This is the method to reach for when the Jacobian is genuinely sparse and its element type is a standard one. It is also what default_linear_solver_method selects in that case.

Constructor

UmfpackLU()

Why sparse is worth it, and when it is not

Measured on an Apple M4 Max, periodic banded matrices of bandwidth 2, factorize! and ldiv! in microseconds against a dense LapackLU on the same matrix:

So the two are a wash around n = 64 and sparse wins by roughly 7× at n = 384 and 12× at n = 1024. A dense matrix handed to this method is an error rather than a conversion — see checksparse — because a SparseMatrixCSC with no structural zeros factorizes slower than LapackLU does.

Against SparspakLU

Sparspak's factorization is the faster of the two, by about 1.3–1.5×. Its triangular solve is about 9× slower, which more than reverses that in a Newton loop, where one factorization is followed by one or more solves. UmfpackLU is the better default; SparspakLU is for the element types UMFPACK cannot do at all.

Allocation

ldiv! is allocation-free. factorize! is not — lu! allocates about 374 kB at n = 384 inside SuiteSparse — and unlike the dense methods that cannot be fixed from here. See SparseFactorizationCache.

Singularity

UMFPACK reports singularity as a status, not as a pivot index, so singular_index is a flag: 0 on success, non-zero otherwise. factorize! records it and ldiv! raises SingularException, matching LapackLU's "report it when it is used" contract so that a quasi-Newton method that factorizes speculatively is not interrupted.

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, NaN, 0, [0, 0, 0, 0, 0, 0])

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

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, NaN, 0, [0, 0, 0, 0, 0, 0])

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).

_bisect_for_minimum(ls, lo, hi, params)

Bisect $\varphi'$ on [lo, hi] for a minimum, returning (α, outcome, n) as _bisection_core does without its ylo.

A bisection converges to whichever crossing the sign of $\varphi'$ at its left endpoint selects, and that sign is invariant under the halving (see _bisection_core): from $\varphi'(lo) < 0$ the interval shrinks onto a $- \to +$ crossing, which is a minimum, and from $\varphi'(lo) > 0$ onto a $+ \to -$ crossing, which is a maximum. Nothing in bracket_minimum rules the second out — it brackets a minimum in value, sampling $\varphi$ and never $\varphi'$, so on a non-convex ray its interval can enclose several stationary points and its left endpoint can sit past one of them.

So when $\varphi'(lo) > 0$ this bisects [0, lo] instead, which brackets a minimum by construction: check_anchor has established $\varphi'(0) < 0$, and $\varphi'(lo) > 0$ gives the opposite sign at the other end, so the crossing between them is a $- \to +$ one. It is also the earlier minimum, the one nearer the anchor, which is the one a line search wants.

The repair costs nothing on the path that does not need it: ylo is a value _bisection_core computes anyway, so a correctly oriented bracket is recognised without a single extra evaluation of $\varphi'$ — which for the $\|F\|^2$ merit of a NonlinearSolver is a full Jacobian. Only the pathological ray pays for a second bisection.

The same condition covers the case where [lo, hi] has no crossing at all and ascends at lo (BISECTION_NOBRACKET with ylo > 0): there too a minimum lies in $(0, lo)$ and there too the bisection of [0, lo] finds it. The retry's verdict then replaces the first one rather than being merged with it, unlike the negative-step retry in solve_with_status — the two disagree here because the first bracket was in the wrong place, which is precisely what ylo detected, and not because $\varphi'$ is inconsistent with $\varphi$.

Private; solve_with_status is the public entry point.

_bracket_core(f, x, bc, s, k, nmax, αmax)

The loop of bracket, returning (lo, hi, n, status) with n the number of evaluations of f it spent and status one of :ok, :capped or :unbracketable. Splitting the loop from the reporting is what bisection/_bisection_core and triple_point_finder/_triple_point_core already do, and here it is what lets a caller tell a bracket that ends at the ceiling αmax from one that satisfied the criterion: the first says the turning point lies beyond the largest step the caller allows, so the answer is αmax itself, and the interval is not worth fitting anything to.

n is reported for the same reason the status is: it is the cost of the bracketing, and for the searches that bracket it is most of the cost of the whole line search. Without it the trials of their LinesearchStatus could only ever be a lower bound — vacuous on the path where the bracketing is the search, which is exactly the path a ceiling produces.

Private: bracket and bracket_minimum are the public entry points.

_bracket_minimum_core(f, x, s, k, nmax, αmax)

bracket_minimum with the (lo, hi, n, status) return of _bracket_core, so that a caller can tell a bracket truncated at the ceiling αmax from one that found a turning point and can charge its cost n to the trials it reports. Private; see _bracket_core.

_bracket_minimum_with_fixed_point_core(f, x, s, k, nmax, αmax)

bracket_minimum_with_fixed_point with the number of evaluations of f it spent and the bracket's status appended — :ok, :capped or :unbracketable, as _bracket_core reports them. Private; the split exists so that Quadratic can tell a bracket that ends at the ceiling from one that found a turning point, and hand back the ceiling instead of fitting a polynomial to an interval over which the merit only falls.

_getrf!(method, cache)

Compute the LU factors of cache.A in place and record the zero-pivot index in cache.info.

This is the only thing that differs between the PivotedLUMethods; everything else — the cache, ldiv!, singular_index and every solve! form — is shared. See PivotedLUCache.

_pivoted_lu_cache(A)

Build a PivotedLUCache from A.

Any AbstractMatrix storage is accepted, but the cache always holds a plain Matrix, because that is what both backends can be pointed at. It is a copy rather than undef so that the single-argument factorize! has something to factorize — as it does for LU, whose cache is likewise seeded from A.

_solve_core!(x, s, state, params)

Run the nonlinear iteration and return the NonlinearSolverStatus it ends in, having already reported it through nonlinear_solver_warnings and print_status.

This is the shared body of solve!, which returns x, and of solve_with_status!, which returns the status. It exists so that the status is built once: the loop has to build one anyway to report on itself, and having solve_with_status! ask status for a second one afterwards spent five l2norm passes per solve reproducing a value it had just discarded — on exactly the path whose reason for existing is to keep per-solve work out of a caller's loop.

The two readings still cannot disagree, and for a stronger reason than before: nothing touches the state between the end of the loop and the caller, so a later status(s, state) rebuilds the same value from the same fields. That is asserted in the test suite.

Private: solve! and solve_with_status! are the public entry points.

_sparse_ldiv!(method, cache, x, b)

Solve with the factorization in cache, writing the result to x, and report a singular factorization as a SingularException.

The per-method half of ldiv! for a SparseDirectMethod; the shared half does the guards. The two backends differ in when they can tell that the matrix was singular, which is why this is not shared. See singular_index.

_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!.

Size is not the only condition: an MArray cannot setindex! a non-isbitstype element, so a BigFloat matrix gets a plain Matrix cache at every size. Without that, the default LU() — which default_linear_solver_method selects for exactly those element types — built an MMatrix cache for anything up to N_STATIC_THRESHOLD and then failed inside factorize! with StaticArrays' "setindex!() with non-isbitstype eltype is not supported", a long way from the choice that caused it.

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

add_to_diagonal!(A, α)

Add α to every diagonal entry of A, in place.

Returns immediately for α == 0, which is the default (SimpleSolvers.REGULARIZATION_FACTOR), so the regularization step of SimpleSolvers.maybe_refactorize! costs nothing unless it was asked for.

For a sparse matrix the diagonal has to be structurally present; a linear-index view into a SparseMatrixCSC, as the dense path uses, would both be wrong for a structural zero and cost O(log nnz) per entry. A missing diagonal entry is an error, because a regularized Jacobian needs somewhere to put the shift.

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_rhs(A)

Allocate a dense right-hand-side/solution vector of length size(A, 1) for the matrix A.

Deliberately dense even when A is not. The obvious spelling, alloc_x(A[:, 1]), gives a sparse vector for a SparseMatrixCSC, and NaN * 0 != 0, so the sparse broadcast has to store every entry anyway — a sparse vector with no structural zeros, which is strictly worse than the Vector it should have been. The right-hand side of a linear system is dense in every caller here.

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_extrapolation(φ₀, d₀, α, φα, q)

The next trial step of the expansion phase of Backtracking, or α itself to say that the step should not be lengthened.

α/φα is the step that was just accepted. The model is the same quadratic through $\varphi(0)$, $\varphi'(0)$ and $\varphi(\alpha)$ that backtracking_interpolation uses on the first backtrack, and its minimiser is

\[\alpha^\star = \frac{-\varphi'(0)\,\alpha^2}{2\big(\varphi(\alpha) - \varphi(0) - \varphi'(0)\alpha\big)} ,\]

clamped from above to $q\alpha$. A denominator that is negative or zero means the model is not convex — the merit fell at least as fast as its tangent, so it is still dropping steeply — and the step grows by the full factor $q$. A denominator that is not finite is a different thing entirely, namely no model at all, and returns α.

Everything the model needs has already been evaluated, so the decision costs no merit evaluation. That is what the lower bound BACKTRACKING_GROW_MIN is for: unless the step the search would actually try is at least that multiple of $\alpha$, α is returned unchanged and the search stops without spending one. The test is on the clamped step rather than on $\alpha^\star$ itself, which matters only for $q <$ BACKTRACKING_GROW_MIN: there the clamp, not the model, is what decides, and a convex model must not be allowed to buy a growth by $q$ that a non-convex one is refused. A direction scaled like a Newton step has $\alpha^\star \approx \alpha$ at $\alpha = 1$ and therefore pays nothing at all, and a merit sitting at its round-off floor ($\varphi(\alpha) \approx \varphi(0)$) gives $\alpha^\star \approx \alpha/2$, so the model declines to expand into rounding noise without needing a special case for it.

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, αmax)

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.

αmax bounds how far to the right the search may probe (see linesearch_αmax); a bracket truncated by it is returned like any other, and only _bracket_core distinguishes the two.

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,
• αmax: how far to the right the search may probe, by default Inf (see linesearch_αmax). A bracket truncated by it is returned like any other, so a caller that has to tell the two apart uses _bracket_minimum_core.

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, αmax)

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]): that line search fits the polynomial centred at the bracket's left endpoint $a$ (which may differ from the input x after the bracketer’s initial direction flip),

\[p(\alpha) = f(a) + f'(a)(\alpha - a) + p_2(\alpha - a)^2,\]

so interpolating $f$ at the right endpoint $b$ fixes the coefficient $p_2$ as

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

where $(a, b)$ is the bracket returned by $\mathtt{bracket\_minimum\_with\_fixed\_point}$. 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.

αmax bounds how far to the right the search may probe; see linesearch_αmax and _bracket_minimum_with_fixed_point_core, which is what tells a truncated bracket from a genuine one.

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)

capped_status(prob, params, αmax, φ₀, d₀, τ, n=0)

The LinesearchStatus of a search whose bracketing reached the ceiling αmax (see linesearch_αmax) with the merit still falling. The turning point then lies beyond the largest step the caller allows, so αmax is the best admissible step, and this is the shared definition of what to report about it.

It is deliberately not a failure and not a distinct LinesearchOutcome. The merit is evaluated at αmax and classified by exactly the rule every other returned step is classified by: LINESEARCH_DECREASED when it beats $\varphi(0)$ by more than the round-off allowance $\tau$, LINESEARCH_FLOOR when it does not. That is honest in both directions — on a compact merit, where this case arises, $\varphi(\alpha_\mathrm{max})$ is genuinely lower and the step genuinely decreases the merit — and a caller that wants to know whether its ceiling bound the search can compare the ceiling it supplied against steplength.

n is what the bracketing that reached the ceiling spent, which the caller gets from the bracketing core (see _bracket_core). It has to be passed in rather than assumed, because on this path the bracketing is the search: reporting only the single evaluation made here would make trials say a capped search cost one step whatever it actually cost.

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₀, α, αmax=Inf)

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. It also respects αmax (see linesearch_αmax), so that the ceiling holds on every return of a line search and not only on the ones that searched.

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

checkfactorized(lsolver::LinearSolver{T,<:PivotedLUMethod})

Throw an ArgumentError if factorize! has not been called on lsolver yet.

The PivotedLUMethod counterpart of the perms[1] == 0 guard in LU's ldiv!: without it, getrs would be handed an all-zero pivot vector and return garbage rather than complain.

checkfactorized(lsolver::LinearSolver{T,<:SparseDirectMethod})

Throw an ArgumentError if factorize! has not been called on lsolver yet.

The SparseDirectMethod counterpart of the PivotedLUMethod guard.

checkjacobianprototype(jacobian, jacobian_prototype)

Throw an ArgumentError if a sparse jacobian_prototype is paired with a Jacobian that cannot write into one.

A sparse prototype is a promise that the Jacobian will be assembled with a fixed sparsity pattern, and the two Jacobians the package builds for itself cannot keep it: JacobianAutodiff hands ForwardDiff a dense matrix and JacobianFiniteDifferences fills one column at a time, so both would write to structurally-zero positions. Caught here, at construction, rather than inside jacobian! on the first iteration.

Those two are named rather than JacobianFunction being allow-listed, so that a caller's own Jacobian subtype — which knows perfectly well how to assemble into a fixed pattern — is not refused for not being one of ours.

checkpattern(lsolver, A)

Throw a DimensionMismatch unless A matches the size the cache was built for.

The pattern itself is checked by the backend, which is where the useful error lives: both are told the pattern may not change, because the ordering and symbolic factorization the cache holds were computed for one.

checksparse(method, A)

Throw an ArgumentError unless A is a square SparseMatrixCSC.

The SparseDirectMethods need the sparsity pattern, and handing one a dense matrix is a mistake worth naming rather than silently converting: converting would build a SparseMatrixCSC with no structural zeros, whose factorization is slower than LapackLU's in every measurement.

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].

convert_αmax(T, αmax)

Convert the step-length ceiling αmax to the precision T, saturating a finite value at floatmax(T) instead of letting it overflow to Inf.

Inf is passed through unchanged: it is not an overflow but a statement — "no ceiling of my own" (see linesearch_αmax) — and turning it into floatmax(T) would impose one the caller declined. Used by change_precision for the αmax field of every method that has one.

copy_matrix!(dest, src)

Copy src into dest, preserving dest's storage.

The sparse-aware replacement for dest .= src. Sparse-to-sparse copies the stored values only, and requires the two patterns to be identical — a differing pattern is an error rather than a silent reallocation, because the solver cache built from dest holds an ordering and a symbolic factorization computed for dest's pattern, and quietly changing it underneath would either be wrong or turn every step into a fresh O(nnz log nnz) ordering.

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.

Reached only from a direct solve call, since that is the only caller of linesearch_warnings. A NonlinearSolver at verbosity = 2 therefore no longer pays for this once per iteration; to see it for a step of a solve, call the line search on that step's problem directly.

default_linear_solver_method(A)

The LinearSolverMethod to use for a system whose matrix looks like A.

This is what NewtonSolver and DogLegSolver fall back on when no linear_solver_method is given, and it dispatches on the Jacobian prototype rather than on the element type alone, because the right answer depends on the storage as much as on the number type:

A sparse matrix is never densified for you. UmfpackLU, the sparse default, solves Float64 and ComplexF64 systems only — UMFPACK converts a 32-bit matrix in lu/lu! but has no 32-bit solve, and it does not handle BigFloat or Rational at all. For every other element type this raises, naming the two things a caller might have meant: SparspakLU, which is generic in the element type and keeps the matrix sparse, or a dense method that discards the sparsity.

Densifying is a real answer — it is often the right one for a small matrix — but it throws away the structure the caller went to the trouble of building, and it is a decision that belongs to them rather than to a fallback. Choosing SparspakLU instead is not available to a default either: it lives in a package extension, so a default that reached for it would work or fail depending on what the caller had imported. So this asks. Once asked, both answers work: pass either as linear_solver_method.

LU used to be the default everywhere. It is still the one for a dense matrix of an element type LAPACK does not know, and the right choice for a very small system, where its static-matrix cache allocates nothing. (A Rational or Integer matrix reaches it and is then refused by lucache_eltype, which names the conversion to make — the package has no dense method for those at all, and saying so is better than a default that pretends otherwise.)

But its allocation-free MMatrix path stops at N_STATIC_THRESHOLD = 10 — and at any size for a non-isbitstype element type, see _static — and above that it is a scalar triple loop with no blocking. Measured on an Apple M4 Max it is 2× slower than LapackLU at n = 64 and 32× slower at n = 768, with its triangular solve a further 3.5–4.5× behind getrs throughout. A caller who did not know to pass linear_solver_method paid all of that; downstream, that was 74 % of an implicit time step.

RecursiveLU is never selected automatically: it lives in a package extension, its useful range depends on which BLAS is loaded, and it is not always installed. Choose it explicitly.

default_linesearch_αmax(T)

DEFAULT_LINESEARCH_αmax in the precision T, saturated at floatmax(T).

The saturation is the whole point of the function. 65536 exceeds floatmax(Float16) = 65504, so Float16(DEFAULT_LINESEARCH_αmax) is Inf — a Float16 line search built from the plain conversion would carry no ceiling at all, which is precisely the defect the field exists to fix, absent in the one precision where every other tolerance in this file is special-cased.

default_linesearch_αmax(Float16)

# output

Float16(6.55e4)

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> round.(direction₂(cache(s)); digits = 11)
2-element Vector{Float64}:
 -0.22882877718
  0.22882877718

The rounding is not decoration. direction₂ comes out of the LU solve, so its last bit carries the accumulation order of whichever BLAS kernel is underneath — and that varies by machine, not by Julia version. Mathematically the two components of this one are exact negatives; printed to sixteen significant digits they are -0.22882877718014286 and 0.22882877718014286 on aarch64 (identically on Julia 1.10, 1.12 and 1.13) and -0.22882877718014286 and 0.22882877718014288 on x86_64 OpenBLAS, which is what CI runs on. Asserting that last digit is asserting the noise rather than the result. direction₁ above needs no such treatment: it is mul! and norms, and never touches the factorization.

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.

dominant_linesearch_outcome(status, count_floor=true)

The non-isbenign LinesearchOutcome the line search reported most often during the solve, or nothing if it reported none. This is what nonlinear_solver_warnings names when it explains a solve that did not converge: a solve whose line search failed does so for one reason nearly every time, and naming that one reason is more use than a histogram.

Ties go to the outcome declared first in LinesearchOutcome, which orders them from the benign end (LINESEARCH_FLOOR, the merit is simply irreducible) towards the actionable one (LINESEARCH_NO_DESCENT, the direction is wrong and the Jacobian is the suspect) — so a tie is broken away from the more alarming diagnosis rather than towards it.

count_floor = false skips LINESEARCH_FLOOR, which is how a converged solve is asked whether anything went wrong: reaching the merit's round-off floor on the last step is how a solve converges, so for that question the floor is not a failure at all. It is the tie rule above that makes the distinction matter — a converged solve that floored once and was exhausted once would otherwise be explained by the floor, which is the half of it that is expected.

f_stall_factor(o::Options)

The factor by which the residual has to drop for an iteration to count as progress. See F_STALL_FACTOR and record_progress!.

f_stall_window(o::Options)

The number of iterations without progress after which a NonlinearSolver gives up (0 disables the criterion). See F_STALL_WINDOW and no_progress.

factorization(lsolver::LinearSolver{T,<:PivotedLUMethod})

A LinearAlgebra.LU view of the factorization held in the cache.

This wraps the cache's arrays rather than copying them, so it is only valid until the next factorize!. Neither ldiv! nor singular_index goes through it — building it allocates, and they are on the hot path — but it is what to reach for when a LinearAlgebra.Factorization is what you want, e.g. det(factorization(lsolver)).

factorize!(lsolver::LinearSolver{T,<:PivotedLUMethod}[, A])

Factorize in place in cache(lsolver).A, with whichever kernel the method selects — see _getrf!. With two arguments A is first copied into the cache; with one, whatever the cache already holds is factorized.

The factorization is not checked here. A singular matrix is reported when the factorization is used, by ldiv!, so that a caller that factorizes speculatively — as a quasi-Newton method does — is not interrupted by a matrix it may never solve with.

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.

factorize!(lsolver::LinearSolver{T,UmfpackLU}, A)

Refactorize A, reusing the ordering and symbolic factorization already in the cache.

SparseArrays' lu! reuses the ordering and symbolic factorization the cache holds, which is the SparseFactorizationCache contract, and raises ArgumentError: pattern of the matrix changed if the pattern does not match the one they were computed for — which is why the error worth reading comes from the backend and checkpattern only checks the size. Keeping the pattern fixed is still the caller's job, and the reason LinearSolver construction takes the prototype.

fill_nan!(A)

Fill A with NaNs.

For a sparse matrix only the stored entries are filled. That is not a shortcut: for a sparse Jacobian the sparsity pattern is structural information the linear solver depends on — the ordering and symbolic factorization in a SparseFactorizationCache were computed for one pattern — so clearing has to preserve it. fill!(A, NaN) would throw for a SparseMatrixCSC anyway, since NaN is not the structural zero.

Used wherever a cache or problem is initialized or cleared: clear!(::LinearProblem) and the initialize! methods of NonlinearSolverCache and DogLegCache.

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.

havenonfinite(status)

Check whether any of the three residuals of a NonlinearSolverStatus — rxₛ, rfₐ, rfₛ — is not finite, i.e. whether the iteration has left the region where the problem is representable. Used by meets_stopping_criteria to give up and by nonlinear_solver_warnings to say so.

The test is isfinite, not !isnan: a residual that has overflowed is as unusable as an undefined one, and the pure-NaN test used to miss it entirely — neither this predicate nor the rfₐ > f_abstol_break gate fired for an infinite residual, since f_abstol_break defaults to Inf and Inf > Inf is false, so such a solve ran its whole max_iterations budget with no diagnosis at all. The same widening was made to the solver-side guards; see nan_recovery!.

Note that a status is never converged by accident here — every comparison with NaN is false and no infinite residual passes residual_small — so this predicate decides when to stop and what to report, not whether the answer is good.

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.

isbenign(oc)

true for the LinesearchOutcomes that report no failure: LINESEARCH_DECREASED (a genuine decrease), LINESEARCH_STATIONARY (nothing to search for) and LINESEARCH_UNKNOWN (the method does not report one). The remaining three — LINESEARCH_FLOOR, LINESEARCH_EXHAUSTED and LINESEARCH_NO_DESCENT — are what linesearch_failures counts and what linesearch_warnings reports on.

LINESEARCH_FLOOR counts as a failure here even though it is the expected final state of a converged solve, because whether it matters is the outer iteration's call and this is how the outer iteration is told: the tally is read by nonlinear_solver_warnings, which names it only for a solve that did not converge.

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!).

isnotprogressing(status)

Check whether the iteration has been given up on for lack of progress: config.f_stall_window iterations without the residual dropping by config.f_stall_factor, see no_progress. Always false at the default f_stall_window = 0.

Mutually exclusive with isconverged, for the same reason isstalled is: the criterion requires the residual not to be small, whereas both convergence branches require that it is.

As with isstalled, whether this counts as failure is the caller's decision: the solve reached status.rfₐ and could not do better within the window it was given.

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.

An infinite step never counts as settled, even though Inf ≤ ‖x‖·x_suctol holds once ‖x‖ has overflowed too: an iterate that jumped to infinity has neither converged nor frozen, it has broken down, and that is what havenonfinite is for.

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.

iterations_since_progress(state)

Return the number of iterations since the residual last dropped by config.f_stall_factor, as recorded in state::NonlinearSolverState by record_progress!. Zero on a freshly initialized state.

This measures the failure mode stall_number cannot see. A stalled step is one that did not move the iterate; here the iterate moves perfectly normally — by far more than the round-off level of $x$ — while the residual descends towards a floor above the requested tolerance, so neither stalled_step nor either branch of assess_convergence can fire and the solve spends max_iterations in full. nonlinear_solver_warnings reports it, and no_progress stops it when the caller has opted in with f_stall_window.

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, LapackLU, SimpleSolvers.PivotedLUCache{Float64, Matrix{Float64}}}(LapackLU(), SimpleSolvers.PivotedLUCache{Float64, Matrix{Float64}}([NaN NaN NaN; NaN NaN NaN; NaN NaN NaN], [0, 0, 0], 0, false))

linesearch_failures(state)

The number of steps in this solve whose line search reported a failure, i.e. the tally of linesearch_outcomes summed over the outcomes that are not isbenign.

linesearch_failures(status)

The number of steps whose line search reported a failure, i.e. linesearch_outcomes summed over the outcomes that are not isbenign.

linesearch_index(oc)

The index of the LinesearchOutcome oc in a tally of length NLINESEARCH_OUTCOMES — the enum's value plus one, since the enum starts at zero and Julia indexes from one. See record_linesearch! and linesearch_outcomes.

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_outcomes(state)

Return the tally of LinesearchOutcomes reported during this solve, indexed by linesearch_index. Maintained by record_linesearch! and reset per solve by initialize!; NonlinearSolverStatus carries a copy of it out to the caller.

linesearch_outcomes(status)

Return the tally of LinesearchOutcomes the line search reported during the solve, indexed by linesearch_index. See record_linesearch! for why this, and not the log, is how a caller learns that the line search was rejected.

Examples

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

julia> x = [1.0];

julia> st = solve_with_status!(x, NonlinearProblem(F, zero(x)), Newton(); verbosity = 0);

julia> linesearch_outcomes(st)[linesearch_index(LINESEARCH_NO_DESCENT)]
0

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.

It may carry an optional αmax field too, the caller's ceiling on the step length; see linesearch_αmax, which documents why that one has to be per call. That one is read by the method rather than by these closures, but through the same hasproperty guard, resolved from the parameter type at compile time, so supplying neither costs nothing. params therefore has to answer property access throughout — a NamedTuple or any struct, which is what the required x and parameters fields already demanded.

linesearch_reason(status, config, oc=dominant_linesearch_outcome(status))

The clause nonlinear_solver_warnings appends to explain a failed solve in terms of what its line search reported, or "" when it reported nothing but success. oc is the outcome to explain — by default the dominant one, and for a converged solve the dominant one other than LINESEARCH_FLOOR, which is what made that message fire (see dominant_linesearch_outcome); the clause has to name the outcome the caller acted on, not a different one that happens to be as frequent.

This is the only thing said about the line search during a solve. A line search does not log from inside the iteration — it reports to the solver, which tallies the outcomes (see record_linesearch!) — so without this clause a solve that stagnated because every one of its steps was rejected would name the symptom and not the cause. A count is what makes it evidence: "the line search rejected 194 of 200 steps" is a diagnosis, "the line search failed once" is noise from the last step of an otherwise healthy solve.

Like no_progress_reason and linesearch_exhausted_reason, this is called from inside the @warn message so the string is built only for a message that is actually shown.

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 every LinesearchMethod reports identically.

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.

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, the tolerance that was requested, and what the line search reported along the way.

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.

linesearch_αmax(method, params)

The ceiling on the step length a LinesearchMethod may return: the smaller of the method's own method_αmax and the caller's params.αmax, if the caller supplied one. This is the single definition of that policy, as check_anchor is of the anchor policy, and every built-in method calls it once before it searches.

Why the caller's half has to be per call

A bracketing search grows its bracket outward until the merit stops falling, so where $\varphi$ is nearly flat — or its minimiser genuinely far away — the step it returns is bounded only by the bracketing budget: $2^{100}$ times the initial step, in principle. In a Euclidean problem that is self-correcting, since $\varphi$ grows with $\alpha$ and the search's own decrease test throws such a step out. On a compact manifold it is not: $\varphi$ is bounded there and can be genuinely lower at $\alpha = 10^9$ than at $\alpha = 0$, so nothing the search can measure calls the step too large. The merit is not a bound on the step, and the bound that does exist belongs to the caller — for a manifold solver, the $2\pi$ of a rotation, divided by the norm of the direction, which changes at every solver step.

Hence the two halves. method_αmax is a backstop that costs nothing and needs no knowledge of the problem; params.αmax is how a caller with a scale of its own imposes it, and it is read from params rather than taken as a keyword so that the extension point solve_with_status(ls, α, params) keeps its signature — a third-party method that does not know about the ceiling still compiles, and one that does needs only this call. It is the same channel params.φ₀ uses (see linesearch_problem), and the hasproperty guard is resolved at compile time, so a caller that supplies nothing pays nothing.

Example

solve_with_status(ls, one(T), (x = x, parameters = params, αmax = 2π / norm(direction)))

The other fields are the ones the merit closures of linesearch_problem read; αmax is additional and independent of them.

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.

lucache_matrix(static, A, Tf)

The working matrix for a LUSolverCache: an MMatrix if static, else a plain Matrix.

An explicit static = true for a non-isbitstype element type is an ArgumentError here rather than StaticArrays' setindex! failure later; the default never asks for it, see _static.

Note the Matrix in the non-static branch, rather than a broadcast that would preserve the input's storage. A sparse A would otherwise give a sparse cache, and factorize!'s scalar loops write to positions that are structurally zero — so it would fail deep inside the factorization, a long way from the cause. LapackLU densifies for the same reason. Use UmfpackLU or SparspakLU to actually exploit sparsity.

Densifying only ever happens because the caller asked for LU by name: default_linear_solver_method never selects a dense method for a sparse matrix, it raises instead.

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,
• the status is making no progress (checked with isnotprogressing, i.e. config.f_stall_window iterations without the residual dropping by config.f_stall_factor while it is not small) and state.iterations ≥ config.min_iterations; this is opt-in and never fires at the default f_stall_window = 0,
• 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 not finite (checked with havenonfinite) 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.