Newton's Method

NewtonsMethod(;
    max_iters = 1,
    update_j = UpdateEvery(NewNewtonIteration),
    krylov_method = nothing,
    convergence_checker = nothing,
    verbose = Silent(),
    line_search = false,
)

Solve f(x) = 0 by iterating x[n+1] = x[n] - j(x[n]) \ f(x[n]), where j(x) = f'(x) is the Jacobian. Called via solve_newton!(method, cache, x, f!, j!, prepare_for_f!); x is modified in-place from its initial guess.

Keyword Arguments

• max_iters: maximum Newton iterations (default 1)
• update_j: UpdateSignalHandler controlling when the Jacobian is recomputed (see Jacobian update strategies below)
• krylov_method: a KrylovMethod to solve the linear system approximately. If nothing, uses direct ldiv! (see Krylov variant below)
• convergence_checker: a ConvergenceChecker that can terminate early based on x[n] and Δx[n]. Without one, always runs max_iters iterations; if verbose is set, a warning is printed when convergence has not been reached by max_iters.
• verbose: Verbose() to print ‖x‖ and ‖Δx‖ each iteration
• line_search: a LineSearch instance to apply backtracking (halving up to 5×) when the residual norm does not decrease or becomes NaN. Default nothing (no line search).

Jacobian update strategies

The update_j parameter accepts any UpdateSignalHandler:

• UpdateEvery(NewNewtonIteration) — fresh Jacobian every iteration (default)
• UpdateEvery(NewNewtonSolve) — reuse across iterations within one solve (the chord method: j(x[n]) ≈ j(x₀))
• UpdateEvery(NewTimeStep) — reuse across solves within a timestep

External signals can also update the Jacobian between solves via update!(method, cache, signal, j!).

Krylov variant

When krylov_method is set, Δx[n] is computed approximately — this is a Newton-Krylov method. If the Krylov method additionally uses a Jacobian-free JVP (see ForwardDiffJVP), neither j_prototype nor j! need to be specified (Jacobian-free Newton-Krylov). When both a JVP and j are provided, j serves as a left preconditioner.

Notes on j_prototype (in allocate_cache)

j_prototype should support ldiv! directly (e.g., a pre-factorized matrix or LinearOperator). Dense matrices are accepted for convenience but trigger an lu factorization on every solve — suitable only for testing. Note that Krylov.jl does not support dense-matrix preconditioners; when using a Jacobian-free JVP, j_prototype must be ldiv!-compatible.

LineSearch(; norm = norm)

Applies a simple backtracking line search to a Newton update.

The line search assumes that x has already been updated by the full Newton step x ← x - Δx, and that f initially contains f(x_old). The full Newton step is evaluated first. If the residual norm does not decrease or becomes non-finite, the step is repeatedly reduced by moving x back toward the previous iterate along the Newton direction. At most five halvings are performed.

The residual norm is measured using norm(f), which can be customized by providing an alternative norm function.

The update is performed in-place. The final (x, f) correspond to the lowest-residual iterate tried; in particular, a non-finite iterate is never returned when a finite one was found (falling back to the previous iterate if no step is finite).

KrylovMethod(;
    type = Val(GmresWorkspace),
    jacobian_free_jvp = nothing,
    forcing_term = ConstantForcing(0),
    args = (),
    kwargs = (; memory = 20),
    solve_kwargs = (;),
    disable_preconditioner = false,
    verbose = Silent(),
    debugger = nothing,
)

Iterative linear solver for Newton's method: finds Δx[n] such that ‖f(x[n]) - j(x[n]) * Δx[n]‖ ≤ rtol[n] * ‖f(x[n])‖, where rtol[n] is controlled by the ForcingTerm. Called via solve_krylov!(method, cache, Δx, x, f!, f, n, prepare_for_f!, j = nothing); Δx is modified in-place. Allocate cache with allocate_cache(method, x_prototype).

This is a wrapper around Krylov.jl solvers. By default, GMRES is used with a Krylov subspace of size 20.

Keyword Arguments

• type: Krylov solver type, wrapped in Val (default Val(GmresWorkspace)).
• jacobian_free_jvp: a JacobianFreeJVP for matrix-free operation (default nothing → uses j directly)
• forcing_term: a ForcingTerm setting rtol[n] (default ConstantForcing(0) → exact solve)
• args, kwargs: forwarded to the Krylov.KrylovSolver constructor (default args = (), kwargs = (; memory = 20) → GMRES subspace size 20)
• solve_kwargs: forwarded to Krylov.solve!
• disable_preconditioner: if true, skip preconditioning even when j is available (default false)
• verbose: Verbose() to print the Krylov residual each iteration
• debugger: a KrylovMethodDebugger run before each Krylov solve
• preconditioner: a custom left preconditioner M (e.g., a matrix-free preconditioner or block-diagonal operator) supporting ldiv! (default nothing)

Operator construction

The solver operates on a LinearOperator opj representing j(x[n]):

• With jacobian_free_jvp: opj evaluates mul!(jΔx, opj, Δx) via the JVP (e.g., finite-difference or AD), so no explicit Jacobian is needed.
• Without: opj wraps j directly, so mul! reduces to mul!(jΔx, j, Δx).

Preconditioning

When disable_preconditioner is false and a preconditioner is provided, it is used as the left preconditioner M. Otherwise, when both a jacobian_free_jvp and an explicit j are provided (and disable_preconditioner is false), j is used as a left preconditioner M. The solver calls ldiv!(Δx′, M, f′) (not mul!), so M is treated as an approximation of j rather than as its inverse. If no preconditioner or j is available, or preconditioning is disabled, M = I.

Tolerances

atol is fixed to 0 so the convergence criterion remains purely relative: ‖r‖ ≤ rtol * ‖f‖. A nonzero atol would add a constant floor that prevents the forcing term from driving the residual to zero, breaking the convergence guarantees of the Newton-Krylov method.

Convergence failures

A failed Krylov solve (most often a singular or inconsistent Jacobian) is not raised or warned about: warning from inside the hot, GPU-dispatched Newton loop would be both noisy and problematic on accelerators. Instead, Δx is set to the least-squares solution returned by Krylov.jl, and the failure is reported only at @debug level. To diagnose suspected non-convergence, set verbose = Verbose() (to print the Krylov residual each iteration) or attach a KrylovMethodDebugger such as PrintConditionNumber.

Extensibility

All constructor and solver arguments can be overridden via args, kwargs, and solve_kwargs, so any Krylov.jl feature not explicitly covered by this wrapper remains accessible.

ForcingTerm

Abstract type for the relative tolerance schedule rtol[n] used by KrylovMethod inside Newton-Krylov iterations. Called via get_rtol!(method, cache, f, n), which returns rtol[n]. Allocate cache with allocate_cache(method, x_prototype).

See ConstantForcing, EisenstatWalkerForcing, and Eisenstat & Walker (1996) for convergence guarantees.

ConstantForcing(rtol)

A ForcingTerm that returns the fixed value rtol ∈ [0, 1) on every Newton iteration.

Convergence properties

• rtol > 0: linear convergence with asymptotic rate ≤ rtol
• rtol = 0: exact Krylov solve → quadratic Newton convergence

Smaller rtol gives faster asymptotic convergence but increases the risk of oversolving (spending Krylov iterations on accuracy that Newton discards).

EisenstatWalkerForcing(;
    initial_rtol = 0.5,
    γ = 1,
    α = 2,
    min_rtol_threshold = 0.1,
    max_rtol = 0.9,
)

Adaptive ForcingTerm ("Choice 2" from Eisenstat & Walker, 1996) that automatically tightens rtol[n] as ‖f(x[n])‖ decreases, balancing convergence speed against oversolving risk.

Keyword Arguments

• initial_rtol ∈ [0, 1): tolerance for the first Newton iteration
• γ ∈ [0, 1]: scaling factor for the tolerance update
• α ∈ (1, 2]: convergence-order exponent — larger means faster convergence but higher oversolving risk
• min_rtol_threshold ∈ [0, 1): safeguard against tolerance decreasing too quickly
• max_rtol ∈ [0, 1): upper bound on rtol[n]

Notes

This is "Choice 2" (not "Choice 1") because it only requires ‖f(x[n])‖ to compute rtol[n], whereas "Choice 1" also needs the final Krylov residual.

JacobianFreeJVP

Abstract type for matrix-free Jacobian-vector product strategies. Subtypes compute j(x) * Δx without forming j explicitly, using only function evaluations of f. Called via jvp!(method, cache, jΔx, Δx, x, f!, f, prepare_for_f!); jΔx is modified in-place. Allocate cache with allocate_cache(method, x_prototype).

See ForwardDiffJVP.

ForwardDiffJVP(; default_step = ForwardDiffStepSize3(), step_adjustment = 1)

A JacobianFreeJVP that approximates the Jacobian-vector product via first-order forward differences:

`j(x) * Δx ≈ (f(x + ε Δx) - f(x)) / ε`,

where ε = step_adjustment * default_step(Δx, x).

Keyword Arguments

• default_step: a ForwardDiffStepSize (default ForwardDiffStepSize3)
• step_adjustment: multiplicative scaling factor for ε (default 1)

ForwardDiffStepSize

Abstract type for step-size strategies used by ForwardDiffJVP. Subtypes are callable: ε = step_size(Δx, x) returns the base step ε before any step_adjustment scaling.

See ForwardDiffStepSize1, ForwardDiffStepSize2, ForwardDiffStepSize3.

ForwardDiffStepSize1()

A ForwardDiffStepSize derived from a truncation-vs-roundoff error analysis of the forward difference approximation j(x) * Δx ≈ (f(x + ε Δx) - f(x)) / ε. Not commonly used with Newton-Krylov methods in practice, but provides intuition for setting step_adjustment in a ForwardDiffJVP.

Reference: Oregon State roundoff/truncation notes.

Returns

The optimal step size that minimizes the error upper bound is

`ε = step_adjustment * sqrt(eps(FT)) / ‖Δx‖`,

where step_adjustment = 2 * sqrt(S * R) (default 1). Increase step_adjustment when f is very smooth (S ≫ 1) or has large roundoff (R ≫ 1). For a central difference approximation, the sqrt becomes a cube root (generally, an n-th root for order n - 1).

Derivation

Forward difference error decomposition

The first-order Taylor expansion of f(x + ε Δx) around x is

`f(x + ε Δx) = f(x) + j(x)(ε Δx) + e_trunc(x, ε Δx)`,

where j(x) = f'(x). In floating point we can only evaluate f̂(x), with

`f(x) = f̂(x) + e_round(x)`.

Substituting and rearranging gives the approximation error

`‖error‖ = ‖e_trunc(x, ε Δx) - e_round(x + ε Δx) + e_round(x)‖ / ε`.

Applying the triangle inequality and approximating ‖e_round(x + ε Δx)‖ ≈ ‖e_round(x)‖ for small ε:

`‖error‖ ≤ (‖e_trunc(x, ε Δx)‖ + 2 ‖e_round(x)‖) / ε`.

Bounding truncation error

From Taylor's theorem for multivariate vector-valued functions (proof):

`‖e_trunc(x, ε Δx)‖ ≤ sup_{x̂ ∈ X} ‖f''(x̂)‖ / 2 · ‖ε Δx‖²`.

Defining the smoothness parameter S = ‖f(x)‖ / sup ‖f''(x̂)‖ (default S ≈ 1; larger values indicate a small Hessian relative to f):

`‖e_trunc(x, ε Δx)‖ ≤ ε² / (2S) · ‖Δx‖² · ‖f(x)‖`.

Bounding roundoff error

Assuming componentwise roundoff |e_round(x)[i]| ≤ R · eps(FT) · |f(x)[i]| (default R ≈ 1):

`‖e_round(x)‖ ≤ R · eps(FT) · ‖f(x)‖`.

Optimal step size

Substituting both bounds into the error bound:

`‖error‖ ≤ ε/(2S) · ‖Δx‖² · ‖f(x)‖ + 2R · eps(FT) · ‖f(x)‖ / ε`.

Minimizing over ε (set derivative to zero) gives ε = 2√(SR) · √eps(FT) / ‖Δx‖, i.e., step_adjustment = 2√(SR).

ForwardDiffStepSize2()

A ForwardDiffStepSize from Knoll & Keyes, "Jacobian-free Newton–Krylov methods: a survey of approaches and applications". This is the step size used by the Fortran package NITSOL:

`ε = √(eps(FT) * (1 + ‖x‖)) / ‖Δx‖`.

ForwardDiffStepSize3()

A ForwardDiffStepSize from Knoll & Keyes, "Jacobian-free Newton–Krylov methods: a survey of approaches and applications". This is the average step size obtained from element-wise forward differences:

`ε = √eps(FT) · Σᵢ(1 + |xᵢ|) / (length(x) · ‖Δx‖)`.

This is the default step size used by ForwardDiffJVP.

ConvergenceChecker(;
    norm_condition,
    component_condition,
    condition_combiner,
    norm = LinearAlgebra.norm,
)

Checks whether a sequence val[0], val[1], val[2], ... has converged to some limit L, given the errors err[iter] = val[iter] .- L. This is done by calling is_converged!(::ConvergenceChecker, cache, val, err, iter), where val = val[iter] and err = err[iter]. If the value of L is not known, err can be an approximation of err[iter]. The cache for a ConvergenceChecker can be obtained with allocate_cache(::ConvergenceChecker, val_prototype), where val_prototype is similar to val and err.

A ConvergenceChecker can perform two types of checks–-it can check whether norm(val) and norm(err) satisfy some ConvergenceCondition, and it can check whether all the components of abs.(val) and abs.(err) individually satisfy some ConvergenceCondition. These two checks can be combined with either & or |. If one of the checks is not needed, the corresponding ConvergenceCondition can be set to nothing.

Instead of LinearAlgebra.norm, norm can be set to anything that will convert val and err to non-negative scalar values.

ConvergenceCondition

An abstract type for objects that can check whether a sequence of non-negative scalar values val[0], val[1], val[2], ... has converged to some limit L, given the errors err[iter] = |val[iter] - L|.

Every subtype of ConvergenceCondition must define has_converged(::ConvergenceCondition, cache, val, err, iter). The cache, which is set to nothing by default, may be used to store information from previous iterations that is useful for determining convergence. In order to have access to a cache of some particular type, a subtype of ConvergenceCondition should define cache_type(::ConvergenceCondition, ::Type{FT}). To specify on which iterations this cache should be updated, it should define needs_cache_update(::ConvergenceCondition, iter). To specify how the cache should be update on those iterations, it should define updated_cache(::ConvergenceCondition, cache, val, err, iter).

Although cache_type can call promote_type to prevent potential type instability errors, this should be avoided in order to ensure that users write type-stable code.

MaximumError(max_err)

Checks whether err[iter] ≤ max_err. Since err[iter] ≥ 0, this can only be true if max_err ≥ 0.

MaximumRelativeError(max_rel_err)

Checks whether err[iter] ≤ max_rel_err * val[iter]. Since err[iter] ≥ 0 and val[iter] ≥ 0, this can only be true if max_rel_err ≥ 0.

MaximumErrorReduction(max_reduction)

Checks whether err[iter] ≤ max_reduction * err[0] for all iter ≥ 1. Since err[iter] ≥ 0, this can only be true if max_reduction ≥ 0. Also, it must be the case that max_reduction ≤ 1 in order for the sequence to not diverge (i.e., to avoid err[iter] > err[0]).

MinimumRateOfConvergence(rate, order = 1)

Stop iterating when convergence is too slow. Reports has_converged = true (i.e., "stop the iteration") when err[iter] ≥ rate * err[iter - 1]^order for any iter ≥ 1. In other words, if the error fails to decrease fast enough relative to the previous iteration, the solver gives up.

Typical values: rate ∈ (0, 1], order ≥ 1.

MultipleConditions(condition_combiner = all, conditions...)

Checks multiple ConvergenceConditions, combining their results with either all or any.

UpdateSignalHandler

Decide whether to update a value upon receiving an appropriate UpdateSignal. Query via needs_update!(::UpdateSignalHandler, ::UpdateSignal).

UpdateEvery(update_signal_type)

Perform the update whenever needs_update! receives an UpdateSignal whose type is update_signal_type or a subtype of it. The subtype matching means a broad handler like UpdateEvery(WithDSS) fires on any EndOfStage / EndOfStep event as well.

UpdateEveryN(n, update_signal_type, reset_signal_type = Nothing)

Perform the update every n-th time needs_update! receives an UpdateSignal of type update_signal_type. The internal counter cycles through 0, 1, …, n-1, and the update is performed when the counter is 0. If reset_signal_type is specified, the counter is also reset to 0 whenever needs_update! receives a signal of that type.

UpdateEveryDt(dt)

Perform the update whenever needs_update! receives a NewTimeStep signal and the elapsed time since the last update is at least dt. The first NewTimeStep always performs the update.

UpdateSignal

A signal that gets passed to an UpdateSignalHandler whenever a certain operation is performed.

NewTimeStep(t)

The signal for a new time step at time t.

NewNewtonSolve()

Signal emitted at the start of each Newton solve (once per implicit Runge-Kutta stage).

NewNewtonIteration()

The signal for a new iteration of Newton's method.

WithDSS

Broadest DSS-site signal type. Use UpdateEvery(WithDSS) to fire a hook at every dss! call (including pre-implicit and post-initialize_imp!).

EndOfStage

Signal type for state-ready-for-tendency-evaluation moments (post-Newton in implicit stages, right after DSS in explicit-only stages, and end of step). A subset of WithDSS.

EndOfStep

Signal type for the end of a full time step. A subset of EndOfStage.

KrylovMethodDebugger

Abstract type for diagnostic hooks run before each Krylov solve. Called via print_debug!(method, cache, j, M). Allocate cache with allocate_cache(method, x_prototype).

See PrintConditionNumber.

PrintConditionNumber()

A KrylovMethodDebugger that prints cond(j) and, when a preconditioner M is available, cond(M⁻¹ j) (the effective condition number seen by the Krylov solver).