Newton root finding

Deep dive

How it works

In: residual F (and derivative or Jacobian J), start x₀ (plus a bracket in 1D), tolerances · Out: x with F(x) ≈ 0 · State: current x (and bracket); optionally a factored J reused across steps · Cost: one F and J evaluation plus one linear solve per iteration; about 5-6 iterations from a decent guess

1D safeguarded Newton (in the spirit of Numerical Recipes' rtsafe):

1. Start from a bracket [a, b] with f(a)·f(b) < 0 and x inside it.
2. Compute the Newton step x_new = x − f/f′.
3. If x_new leaves the bracket, or the step did not at least halve compared with two steps ago, take the bisection point (a + b)/2 instead.
4. Evaluate f(x_new) and shrink the bracket to the half that keeps the sign change.
5. Stop when |step| ≤ xtol·max(1, |x|) or f = 0.

Systems: compute F and J, solve J·Δ = −F by LU, then backtrack t from 1 until the sufficient-decrease (Armijo) condition holds on φ(x) = ½‖F(x)‖². Near the root t = 1 is accepted and quadratic convergence returns. If J is singular or the line search fails, the iterate is near a local minimum of ‖F‖ that is not a root; restart from elsewhere or switch to a trust-region (Levenberg-Marquardt) step.

Code

/// Safeguarded Newton on a bracket [a, b] with f(a) f(b) <= 0. `fd` returns (f(x), f'(x)).
pub fn newton_bracketed(fd: impl Fn(f64) -> (f64, f64), mut a: f64, mut b: f64, xtol: f64, max_iter: usize) -> Option<(f64, usize)> {
    let (fa, _) = fd(a);
    let (fb, _) = fd(b);
    if fa == 0.0 { return Some((a, 0)); }
    if fb == 0.0 { return Some((b, 0)); }
    if fa * fb > 0.0 { return None; }
    if fa > 0.0 { std::mem::swap(&mut a, &mut b); } // orient: f(a) < 0 < f(b)
    let mut x = 0.5 * (a + b);
    let (mut dx_old, mut dx) = ((b - a).abs(), (b - a).abs());
    let (mut f, mut df) = fd(x);
    for it in 1..=max_iter {
        let newton_out = ((x - b) * df - f) * ((x - a) * df - f) > 0.0; // step lands outside [a, b]
        let too_slow = (2.0 * f).abs() > (dx_old * df).abs();
        dx_old = dx;
        if newton_out || too_slow {
            dx = 0.5 * (b - a);
            x = a + dx;
        } else {
            dx = f / df;
            x -= dx;
        }
        if dx.abs() <= xtol * x.abs().max(1.0) { return Some((x, it)); }
        (f, df) = fd(x);
        if f == 0.0 { return Some((x, it)); }
        if f < 0.0 { a = x; } else { b = x; }
    }
    None
}

/// Solve the dense system J d = -F by Gaussian elimination with partial pivoting. None if singular.
fn solve_neg(mut j: Vec<Vec<f64>>, f: &[f64]) -> Option<Vec<f64>> {
    let n = f.len();
    let mut r: Vec<f64> = f.iter().map(|v| -v).collect();
    for k in 0..n {
        let p = (k..n).max_by(|&x, &y| j[x][k].abs().total_cmp(&j[y][k].abs()))?;
        if j[p][k].abs() < 1e-300 { return None; }
        j.swap(k, p); r.swap(k, p);
        for i in k + 1..n {
            let l = j[i][k] / j[k][k];
            for c in k..n { j[i][c] -= l * j[k][c]; }
            r[i] -= l * r[k];
        }
    }
    for i in (0..n).rev() {
        let s: f64 = (i + 1..n).map(|c| j[i][c] * r[c]).sum();
        r[i] = (r[i] - s) / j[i][i];
    }
    Some(r)
}

/// Forward-difference Jacobian, step sqrt(eps) * max(|x_j|, 1).
pub fn fd_jacobian(f: &impl Fn(&[f64]) -> Vec<f64>, x: &[f64], fx: &[f64]) -> Vec<Vec<f64>> {
    let n = x.len();
    let mut jac = vec![vec![0.0; n]; fx.len()];
    let mut xp = x.to_vec();
    for c in 0..n {
        let h = f64::EPSILON.sqrt() * x[c].abs().max(1.0);
        xp[c] = x[c] + h;
        let h = xp[c] - x[c]; // the step actually representable
        let fp = f(&xp);
        for r in 0..fx.len() { jac[r][c] = (fp[r] - fx[r]) / h; }
        xp[c] = x[c];
    }
    jac
}

pub struct NewtonStats { pub iterations: usize, pub f_evals: usize, pub residual: f64 }

/// Damped Newton with Armijo backtracking on 0.5 * ||F||^2.
pub fn newton_system(f: impl Fn(&[f64]) -> Vec<f64>, jac: Option<&dyn Fn(&[f64]) -> Vec<Vec<f64>>>,
                     x: &mut Vec<f64>, ftol: f64, max_iter: usize) -> Result<NewtonStats, NewtonStats> {
    let norm2 = |v: &[f64]| v.iter().map(|a| a * a).sum::<f64>();
    let mut fx = f(x);
    let mut evals = 1;
    for it in 0..max_iter {
        let phi = 0.5 * norm2(&fx);
        if phi.sqrt() <= ftol { return Ok(NewtonStats { iterations: it, f_evals: evals, residual: phi.sqrt() }); }
        let j = match jac { Some(jf) => jf(x), None => { evals += x.len(); fd_jacobian(&f, x, &fx) } };
        let Some(d) = solve_neg(j, &fx) else { break };
        let mut t = 1.0;
        loop {
            let xt: Vec<f64> = x.iter().zip(&d).map(|(a, b)| a + t * b).collect();
            let ft = f(&xt);
            evals += 1;
            // J d = -F, so the directional derivative of phi is -2 phi: Armijo with alpha = 1e-4
            if 0.5 * norm2(&ft) <= (1.0 - 2e-4 * t) * phi || t < 1e-10 {
                *x = xt; fx = ft; break;
            }
            t *= 0.5;
        }
    }
    let r = norm2(&fx).sqrt();
    if r <= ftol { Ok(NewtonStats { iterations: max_iter, f_evals: evals, residual: r }) }
    else { Err(NewtonStats { iterations: max_iter, f_evals: evals, residual: r }) }
}

Complexity and constants

• 1D: each step costs one f and f′; near the root digits double, so from 2 correct digits about 3 more steps reach 16.
• Measured with the code above: plain Newton for √2 from x₀ = 1 has errors 8.6·10⁻², 2.5·10⁻³, 2.1·10⁻⁶, 1.6·10⁻¹², 0; newton_bracketed finds √2 on [0, 2] in 6 steps, cos x = x on [0, 1] in 4, and the root −1.7693 of the cycling cubic on [−3, 1] in 7. newton_system on circle ∩ (y = eˣ − 1) from (10, 10) needs 17 iterations (the line search cuts early steps), 21 F evaluations with an analytic J and 55 with finite differences; the 50-variable Broyden tridiagonal system from x = −1 converges in 5 iterations to 2·10⁻¹⁵.
• Systems: a dense J costs n extra F evaluations by finite differences (or n forward-mode directional derivatives) plus an O(n³) solve; reusing J for several steps (chord / simplified Newton, as in stiff ODE solvers) trades quadratic for linear convergence at much lower cost.
• Finite-difference J with step √ε·|x|: derivatives accurate to about 8 digits, which still gives fast (superlinear, close to quadratic until 10⁻⁸ relative) convergence.

Efficiency

• Cost: one F and one derivative per 1D step, and for systems a dense J (n extra F evaluations by finite differences) plus an O(n³) solve. Measured with the code above: plain Newton for √2 reaches 1.6·10⁻¹² in 4 steps; newton_system on circle ∩ (y = eˣ − 1) from (10, 10) takes 17 iterations, 21 F evaluations with an analytic J and 55 with finite differences. Bound: cost of F and J, then the linear solve.
• Levers:
  • analytic or autodiff J instead of finite differences — 21 vs 55 F evaluations (measured above) — systems with a few variables
  • chord (simplified) Newton, reuse a factored J — much cheaper steps, linear instead of quadratic convergence (page above) — stiff ODE solvers
  • bracket plus bisection fallback — guaranteed convergence, 6 steps for √2 on [0, 2] (measured above) — unknown start
  • Newton-Krylov with matrix-free J·v — no dense J or O(n³) solve — large PDE systems
  • hardware estimate plus one Newton step — about 23 bits from a 12-bit rsqrtps estimate (page above) — reciprocals and square roots
• Naive → standard → state of the art:
  • bisection · one bit per step, about 52 steps to full f64 precision
  • safeguarded Newton or Newton with backtracking · 4-6 steps (measured above)
  • Newton-Krylov, or chord Newton with a reused factorization · cost per step far below a dense solve

Implementation notes

• Stopping: combine a step test |Δ| ≤ xtol·(|x| + xtol) with a residual test; f alone can be tiny at a far-from-root point when f is flat.
• Scaling: scale variables and equations to similar magnitudes; Newton is invariant to linear transformations in exact arithmetic but line search and stopping tests are not.
• Multiple roots: if convergence is linear with ratio ≈ (m−1)/m, use x ← x − m·f/f′ or apply Newton to f/f′.
• Halley / higher order: cubic convergence with f″; worth it only when f″ is nearly free (e.g. polynomials, special functions).
• Fast inverse sqrt / reciprocal: hardware estimate (rsqrtps, ~12 bits) plus one Newton step y ← y(1.5 − 0.5·x·y²) gives about 23 bits.
• Complex roots and polynomials: deflation is unstable; use companion-matrix eigenvalues or Aberth/Durand-Kerner for all roots.

Pitfalls

• Unsafeguarded Newton on x³ − 2x + 2 from x₀ = 0 cycles 0 → 1 → 0 forever.
• f′(x) = 0 at an iterate (local extremum): division by zero or a huge jump; the bracket check catches it.
• Finite-difference step too small (catastrophic cancellation) or too large (truncation): use √ε scaling with the actual representable step, as in fd_jacobian.
• Line search on ‖F‖² can stall in a local minimum that is not a root; detect and restart, or use homotopy/continuation from an easy problem.
• Inconsistent analytic Jacobian (a sign error in one entry): convergence degrades to linear or fails; test J against finite differences once.

Variants and alternatives

In production

Libraries: SciPy optimize.brentq, newton, root (hybr = MINPACK Powell hybrid, lm), Boost.Math newton_raphson_iterate / halley_iterate, GSL root solvers, PETSc SNES, SUNDIALS KINSOL, Ceres Solver (Levenberg-Marquardt). Electronics: implicit solvers in every circuit simulator (SPICE Newton iterations). Engineering: FEM codes. Robotics and animation: inverse-kinematics solvers (Jacobian-based Newton and damped least squares, Buss).

Implementations

Sources

• William H. Press, Saul A. Teukolsky, William T. Vetterling, Brian P. Flannery, Numerical Recipes: The Art of Scientific Computing, 3rd ed., Cambridge University Press, 2007, chapter 9. book
• J. E. Dennis Jr., Robert B. Schnabel, Numerical Methods for Unconstrained Optimization and Nonlinear Equations, SIAM Classics, 1996. book
• C. T. Kelley, Solving Nonlinear Equations with Newton's Method, SIAM, 2003. book
• Richard P. Brent, Algorithms for Minimization without Derivatives, Prentice-Hall, 1973.
• Jorge Nocedal, Stephen J. Wright, Numerical Optimization, 2nd ed., Springer, 2006, chapter 11. book
• Samuel R. Buss, Introduction to Inverse Kinematics with Jacobian Transpose, Pseudoinverse and Damped Least Squares Methods, 2004. pdf