RasmahRasmah

Nonlinear solvers

Overview

Root finding is the search for the inputs that make a function exactly zero: solve $f(x) = 0$. It shows up everywhere a model has an equilibrium — the displacement where internal forces balance external ones, the temperature where heat in equals heat out, the steady state of a reacting flow. Optimization finds a minimum; root finding finds a zero, and the two are cousins: the minimum of a smooth function is where its gradient is zero.

The classic idea is Newton's method. Start at a guess $x_0$, draw the straight line tangent to $f$ at that point, and jump to where that line crosses zero. Each step solves a linear problem $f'(x_k)\,\Delta = -f(x_k)$ and updates $x_{k+1} = x_k + \Delta$. When the guess is close, the error squares every step — roughly doubling the correct digits. The cost is that you need the derivative $f'$ (or, in many dimensions, the Jacobian) and a linear solve per step.

The unified solver

Every method in this chapter hangs off one entry point, solve_nonlinear, which dispatches on a method symbol:

Rasmah.NonlinearSolve.solve_nonlinear — Function
solve_nonlinear(f, x0; method=:newton, kwargs...) -> NonlinearResult

Solve the nonlinear root problem $f(x) = 0$ starting from x0, dispatching by method across Rasmah's nonlinear-solver layer.

Theory

A root of $f$ is a point $x^*$ with $f(x^*) = 0$. Newton's method (the default) linearizes the residual, $f(x) + J(x)\, \Delta x = 0$, and steps $x_{k+1} = x_k + \Delta x$; it converges quadratically near a root when the Jacobian $J$ is nonsingular. The other methods trade that local speed for robustness: quasi-Newton (Broyden, Anderson) approximates $J$ instead of reforming it, and continuation (pseudo-transient, homotopy) follows a smooth path from a known state to the solution.

Arguments

  • f: the residual f(x) -> vector, or a LinearResidual/InplaceResidual.
  • x0: the starting point.

Keyword options

  • method = :newton: the solver — :newton, :broyden, :anderson, :picard, :newton_krylov, :pseudotransient, :homotopy, :lm, :nlcg, :halley, :chebyshev, :householder, :mpria, … (:anderson/:picard use the fixed-point map G(x) = x - f(x)).
  • jac / jvp: analytic Jacobian / Jacobian-vector product (Newton; otherwise automatic differentiation is used).
  • linesearch, rtol, atol, maxiter: step-size rule and tolerances.

Returns

A NonlinearResult with solution, residual, residual_norm, converged, iterations, and reason.

Example

julia> using Rasmah

julia> r = solve_nonlinear(x -> [x[1]^2 - 4, x[2]^3 - 8], [1.0, 1.0]);

julia> round.(r.solution; digits = 6) == [2.0, 2.0]
true

See also: minimize, newton_solve.

source
solve_nonlinear(x -> [x[1]^2 - 4.0], [1.0]).solution
1-element Vector{Float64}:
 2.000000000000002

The square $x^2 = 4$ has two roots; started near $1$ the solver settles on $x = 2$. Every solver returns the same report card:

Rasmah.NonlinearSolve.NonlinearResult — Type
NonlinearResult{T}

Returned by every *_solve entry point. Fields:

  • solution: converged (or final) iterate.
  • residual: f(solution).
  • converged: whether a convergence criterion was met.
  • iterations: outer iterations taken.
  • function_evals: residual/function evaluations.
  • jacobian_evals: Jacobian assemblies (0 for matrix-free/JFNK).
  • linear_iterations: total inner Krylov/direct iterations.
  • residual_norm: norm(residual).
  • reason: :converged, :max_iter, :stagnation, or :singular.
  • trace: per-iterate solution snapshots when store_trace=true, else nothing.
source
r = solve_nonlinear(x -> [x[1]^2 - 4.0], [1.0])
r.converged, r.iterations, r.residual_norm
(true, 5, 8.881784197001252e-15)

Choosing a method

Newton's method is the default and the workhorse. The method symbol selects a specialist for when Newton's assumptions break:

  • :newton — Newton with a line search or dogleg trust region (default).
  • :broyden — quasi-Newton: approximate the Jacobian instead of assembling it, so expensive derivatives are avoided.
  • :anderson / :picard — fixed-point acceleration: treat $x = G(x)$ and mix the last few iterates to converge a contracting map fast.
  • :newton_krylov (JFNK) — matrix-free: never build the Jacobian, just multiply it against vectors, ideal for huge sparse systems.
  • :pseudotransient — march in a fake time step toward steady state, robust when the initial guess is far away.
  • :homotopy / :pseudo_arclength — continuation: track a family of solutions as a parameter changes, following the solution around turning points where plain Newton would jump or fail.
  • :deflation — find several roots by excluding the ones already found.

These are advanced levers; the reference lists the full menu. For most problems solve_nonlinear(f, x0) is all you need.

Root finding for whole solution families

When the root depends on a parameter $\theta$ — solve $f(\theta, x) = 0$ for $x$ at many values of $\theta$ — Rasmah can build a callable map $\theta \mapsto x^*(\theta)$ and differentiate through it using the implicit function theorem, without re-running the solver for each derivative:

Rasmah.NonlinearSolve.implicit_function — Function
implicit_function(f; solver=nothing, x0=nothing, backend=ad_backend(), kwargs...) -> ImplicitFunction

Build an ImplicitFunction for f(θ, x) = 0. Either pass a solver callable θ -> x*, or an initial guess x0 (then a Newton solver with an AD Jacobian is used; remaining kwargs forward to newton_solve).

source

This is what lets a downstream optimizer treat "the solution of my simulation" as a differentiable function of its inputs.

Next steps

Root finding and optimization share the same linear-algebra backbone. When the bottleneck is the linear solve inside each step, see Sparse linear algebra — and when the goal is to minimize something rather than hit zero, return to Optimization.