newtonsolve(rf!, x0::AbstractVector, [cache::NewtonCache]; tol=1.e-6, maxiter=100, [logger])

Solve the nonlinear equation (system) r(x)=0 using the newton-raphson method by calling the mutating residual function rf!(r, x), with signature rf!(r::T, x::T)::T where T<:AbstractVector x0 is the initial guess and the optional cache can be preallocated by calling NewtonCache(x0). Note that x0 is not modified, unless aliased to getx(cache). tol is the tolerance for norm(r) and maxiter the maximum number of iterations. An AbstractLogger, see e.g. [StandardLogger](@ref), can be passed to log the iterations.

returns x, drdx, converged::Bool

drdx is the derivative of r wrt. x at the returned x.

newtonsolve(rf, x0::T; tol=1.e-6, maxiter=100, [linsolver], [logger]) where {T <: 
    Union{Number, StaticArrays.SVector, Tensors.Vec, Tensors.SecondOrderTensor}}

Solve the nonlinear equation (system) r(x)=0 using the newton-raphson method by calling the residual function r=rf(x), with signature rf(x::T)::T x0::T is the initial guess, tol the tolerance form norm(r), and maxiter the maximum number of iterations.

Optional keyword arguments:

• A non-standard linsolver can optionally be specified, please see Linear solvers for more information.
• An AbstractLogger, see e.g. [StandardLogger](@ref), can be passed to log the iterations

returns: x, drdx, converged::Bool

drdx is the derivative of r wrt. x at the returned x.

function NewtonCache(x::AbstractVector; [linsolver])

Create the cache used by the newtonsolve and linsolve!. Only a copy of x will be used.

A special linsolver can optionally be given, please see Linear solvers for more information.

getx(cache::NewtonCache)

Extract out the unknown values. This can be used to avoid allocations when solving defining the initial guess.

NoLogger()

No-op logger that doesn't do any logging. This is used by default if no logger is passed, and logging is not enabled via Preferences.jl (see Newton.logging_mode).

StandardLogger(x)

The standard logger to show the evolution of residual norms and changes, measured by norm, of the unknowns between iterations.

FullLogger(x)

A logger that logs everything known to the newton iterations,

• The full residual in each iteration
• The full unknowns in each iteration
• The full jacobian in each iteration

The report_logger will still just print the same info as the StandardLogger, but additional information can be obtained by directly accessing the field. This is intended for interactive debugging, and the fields may change.

Newton.logging_mode(; enable=true)

Helper to turn on (enable=true) or off (enable=false) logging of iterations in Newton.jl. This changes the default logger from NoLogger to StandardLogger.

Newton.jl comes with the following linear solvers

• StandardLinsolver (default)
• UnsafeFastLinsolver
• RecursiveFactorizationLinsolver

Newton.StandardLinsolver()

This is the default linear solver, which gives safe operations and don't require any special packages to be loaded.

UnsafeFastLinsolver()

This is a special linear solver, which calculates the inverse recursively by using the analytical inverses of 2x2 and 3x3 matrices. This gives exceptional performance for small matrices, but suffers from numerical errors and can be sensitive to badly conditioned matrices. When using this method, it may be advisable to (adaptively) try a slower method if the newton iterations fail to converge.

RecursiveFactorizationLinsolver()

This linear solver utilizes the LU decomposition in RecursiveFactorization.jl, which gives faster performance than the StandardLinsolver. While not as fast as UnsafeFastLinsolver, it is always accurate. Is available via an extension, requiring the user to load RecursiveFactorization.jl separately.

Newton.inv!(A::Matrix, cache::NewtonCache)

In-place inverse, which, depending on the linsolver in cache, can be much more efficient than inv(A). However, note that A will be used as workspace and values should not be used after calling Newton.inv!. In some cases, A will become its inverse, and the output aliased to A. This behavior is not true in general, and should not be relied upon.

sinv(a::SMatrix{d, d}) where {d}

Fast, but numerically unstable implementation of the inverse of a statically sized matrix. Beneficial up to d ≈ 50, but can give large floating point errors for badly conditioned and/or large matrices.

About 4 to 5 timers faster than StaticArrays for sizes ∈ [5, 20], and at least twice as fast up to 50 according to benchmarks on macbook with M3 processor.

sinv!(K::Matrix)

Invert K in-place using the unsafe static implementation up to a size of 20x20, and fall back to generic LinearAlgebra.inv

ad_newtonsolve(rf, x0, rf_dual_args::Tuple; kwargs...)

Solve rf(x, y₁(z), y₂(z), ...) = 0 to find x(y₁(z), y₂(z), ...), given the initial guess x0 as a non-dual number, and the dual input numbers rf_dual_args = (y₁(z), y₂(z), ...) where all are derivatives wrt. the same variable z. Return x of Dual type seeded such that it corresponds to the derivative dx/dz = ∂x/∂yᵢ ⋆ dyᵢ/dz where ⋆ is the appropriate contraction.

Implementation:

Uses the adjoint, i.e. dr/dyᵢ = 0 = ∂r/∂x ⋆ ∂x/∂yᵢ + ∂r/∂yᵢ ⇒ ∂x/∂yᵢ = -[∂r/∂x]⁻¹ ⋆ ∂r/∂yᵢ, such that we avoid doing newton iterations with dual numbers.

ad_newtonsolve(rf, x0, rf_args::Tuple; kwargs...)

If rf_args do not contain dual numbers, the standard newtonsolver is just called on f(x) = rf(x, y₁, y₂, ...), and the solution x is returned. This allows writing generic code where the gradient is sometimes required, but not always.

Example

using Newton, Tensors, ForwardDiff, BenchmarkTools
rf(x::Vec, a::Number) = a * x - (x ⋅ x) * x
function myfun!(outputs::Vector, inputs::Vector)
    x0 = ones(Vec{2}) # Initial guess
    a = inputs[1] + 2 * inputs[2]
    x, converged = ad_newtonsolve(rf, x0, (a,))
    outputs[1] = x ⋅ x
    outputs[2] = a * x[1]
    return outputs 
end
out = zeros(2); inp = [1.2, 0.5]
ForwardDiff.jacobian(myfun!, out, inp)

gives

2×2 Matrix{Float64}:
 1.0      2.0
 1.57321  3.14643

linsolve!(K::AbstractMatrix, b::AbstractVector, cache::NewtonCache)

Solves the linear equation system Kx=b, mutating both K and b. b is mutated to the solution x.

linsolve!(linsolver, K::AbstractMatrix, b::AbstractVector, cache::NewtonCache)

The default implementation will call this signature, which should be overloaded for a different linsolver passed to cache upon construction.

linsolve(linsolver, K, b)

Solve the equation (system) b = K ⋆ x where ⋆ is the appropriate contraction between K and x, and b and x have the same type and size. The following combinations of K and x are supported.

extract_submatrix(::Type{SMatrix{d1, d2}}, m::SMatrix, start_row, start_col)

Efficiently extract s::SMatrix{d1, d2} such that s == m[start_row:(start_row + d1 - 1), start_col:(start_col + d2 - 1)]

join_submatrices(a11, a12, a21, a22)

Efficiently join the submatrices to return a::SMatrix = [a11 a12; a21 a22].

reset_logger!(logger::AbstractLogger)

Reset the logger at the beginning of the newton iterations.

update_logger!(logger::AbstractLogger, r_norm, r, x, drdx)

Update the logger for the current iteration by passing the residual norm, r_norm, the full residual, r, the unknowns, x, and the jacobian, drdx.

report_logger(logger::AbstractLogger)

Report (print) the information from the logger.