Newton's Method

Newton's method, also called the Newton–Raphson method, is a foundational iterative algorithm in numerical analysis, scientific computing, and mathematical optimization that successively approximates the roots of a differentiable nonlinear function F(x) = 0 or the stationary points ∇f(x) = 0 of a twice-differentiable objective function f(x). Originating in algebraic root approximations described by Isaac Newton in 1669 and cast into modern derivative notation by Joseph Raphson in 1690 and Thomas Simpson in 1740, the method replaces a nonlinear problem at each iterate x_k with its local first-order Taylor linearization, solving for the step Δx_k that zeroes the tangent model.

When initialized within the basin of attraction of a simple root or a strongly convex local minimum, Newton's method exhibits local quadratic convergence: the number of accurate significant digits roughly doubles with every single iteration. This rapid asymptotic convergence distinguishes Newton's method from first-order methods such as gradient descent, which converge only linearly, and motivates both quasi-Newton approximations such as the Broyden–Fletcher–Goldfarb–Shanno algorithm and derivative-free alternatives such as the Nelder–Mead method.

Given a continuously differentiable scalar function g : ℝ → ℝ and an initial estimate x₀ near a root x satisfying g(x) = 0, expanding g in a first-order Taylor series around the current iterate x_k gives the tangent-line approximation:

g(x_k + Δx) ≈ g(x_k) + g′(x_k) Δx = 0.

Assuming the derivative g′(x_k) is non-zero, solving for the x-intercept of this tangent line yields the scalar Newton–Raphson recurrence:

x_k+1 = x_k − g(x_k) / g′(x_k).

Error Analysis and Quadratic Convergence Rate

Let e_k = x_k − x denote the error at step k, and suppose g is twice continuously differentiable near a simple root x where g(x) = 0 and g′(x) ≠ 0. By Taylor's theorem with Lagrange remainder around x_k, there exists a point ξ_k between x_k and x such that 0 = g(x) = g(x_k) − g′(x_k)e_k + (1/2)g″(ξ_k)e_k². Dividing by g′(x_k) and rearranging produces:

e_k+1 = e_k − g(x_k) / g′(x_k) = [g″(ξ_k) / (2g′(x_k))] e_k².

Thus, for x_k close to x, |e_k+1| ≈ C |e_k|² with asymptotic constant C = |g″(x) / (2g′(x^*))|. Once |e_k| ≤ 10⁻² (assuming C ≈ 1), three additional iterations reduce the error sequentially to 10⁻⁴, 10⁻⁸, and 10⁻¹⁶, reaching IEEE 754 double-precision machine accuracy in 4 to 6 steps. In floating-point hardware and math libraries, division-free Newton iterations compute reciprocals 1/a via x_k+1 = x_k(2 − ax_k) and square roots √a via Heron's formula x_k+1 = (1/2)(x_k + a/x_k).

If x is a multiple root of multiplicity m > 1 (so g(x) = g′(x) = ⋯ = g^(m−1)(x) = 0), standard Newton iteration degrades to linear convergence with rate 1 − 1/m unless modified to x_k+1 = x_k − m · g(x_k)/g′(x_k).

Systems of Nonlinear Equations and Kantorovich's Theorem

For a vector-valued mapping F : ℝⁿ → ℝⁿ, the first-order Taylor model around x_k ∈ ℝⁿ is F(x_k) + J_F(x_k) Δx_k = 0, where J_F(x_k) ∈ ℝ^n×n is the Jacobian matrix with entries [J_F]_ij = ∂F_i/∂x_j. Rather than inverting the Jacobian explicitly, each step solves the n × n linear system:

J_F(x_k) Δx_k = −F(x_k), x_k+1 = x_k + Δx_k,

using direct LU decomposition (for dense or moderate-size sparse problems) or an iterative Krylov solver. Leonid Kantorovich (1948) proved a landmark sufficient condition—the Newton–Kantorovich theorem—which guarantees both the existence of a root x^* and quadratic convergence solely from bounds on ‖J_F(x₀)⁻¹F(x₀)‖ and the Lipschitz constant of J_F at the initial point x₀, without assuming in advance that a root exists.

Newton's Method in Unconstrained and Constrained Optimization

To minimize a twice continuously differentiable objective function f : ℝⁿ → ℝ, one seeks a stationary point where the gradient vanishes: F(x) = ∇f(x) = 0. The Jacobian of the gradient is the symmetric Hessian matrix H(x) = ∇²f(x) ∈ ℝ^n×n with entries H_ij = ∂²f/(∂x_i∂x_j). At each iteration, Newton's method minimizes the local quadratic Taylor model m_k(p) = f(x_k) + ∇f(x_k)^Tp + (1/2)p^T∇²f(x_k)p, yielding the Newton equations:

∇²f(x_k) p_k = −∇f(x_k), x_k+1 = x_k + α_k p_k.

When ∇²f(x_k) is symmetric positive-definite, p_k is solved via Cholesky decomposition and is a descent direction (∇f(x_k)^Tp_k < 0). Furthermore, the pure Newton step (α_k = 1) is affine invariant: linearly rescaling the coordinate system x = Ay leaves the trajectory of iterates unchanged, making Newton's method immune to the coordinate ill-conditioning that slows gradient descent.

Far from a solution, the pure Newton step (α_k = 1) can overshoot, oscillate in a 2-cycle, diverge to infinity near inflection points where g′(x_k) ≈ 0, or step toward a saddle point or local maximum if the Hessian ∇²f(x_k) has negative eigenvalues. Modern production solvers globalize Newton's method via two complementary mechanisms:

1. Line-search damping and Levenberg–Marquardt regularization: If ∇²f(x_k) is not positive-definite, adding a positive shift ∇²f(x_k) + λ_kI restores positive-definiteness, and a backtracking Armijo–Wolfe line search selects step size α_k ∈ (0, 1] until f(x_k + α_kp_k) decreases sufficiently, switching to α_k = 1 near the minimum to preserve quadratic convergence.
2. Trust-region Newton–CG (Steihaug–Toint) and Jacobian-Free Newton–Krylov (JFNK): For large-scale PDE systems where forming or factoring J_F(x_k) is prohibitive, each outer Newton step solves J_F(x_k) Δx_k ≈ −F(x_k) inexactly using the generalized minimal residual method or the conjugate gradient method, evaluating matrix-vector products via directional finite differences J_F(x_k)v ≈ [F(x_k + εv) − F(x_k)]/ε or automatic differentiation.

What is the Gauss–Newton method for nonlinear least squares?

When minimizing a sum of squared residuals f(x) = (1/2)‖r(x)‖₂², the full Hessian is ∇²f(x) = J_r(x)^TJ_r(x) + ∑_i r_i(x)∇²r_i(x). The Gauss–Newton method drops the second-derivative term and solves the linear least-squares subproblem J_r(x_k) Δx_k ≈ −r(x_k) via QR decomposition at each step, achieving near-quadratic convergence when the optimal residual ‖r(x^*)‖₂ is small.

Why do Newton fractals appear in the complex plane?

Applying Newton's method to a complex polynomial such as z³ − 1 = 0 partitions the complex plane ℂ into basins of attraction for the three cube roots of unity. By Arthur Cayley's and Gaston Julia's analysis, the boundary between any two basins is simultaneously the boundary of all three basins, forming a self-similar Julia fractal where arbitrarily small perturbations in initial guess z₀ lead to different roots.

1. Kantorovich, L. V. (1948). Functional analysis and applied mathematics. Uspekhi Matematicheskikh Nauk, 3(6), 89–185. https://doi.org/10.6028/jres.049.014
2. Dennis, J. E., & Moré, J. J. (1977). Quasi-Newton methods, motivation and theory. SIAM Review, 19(1), 46–89. https://doi.org/10.1137/1019005
3. Knoll, D. A., & Keyes, D. E. (2004). Jacobian-free Newton–Krylov methods: A survey of approaches and applications. Journal of Computational Physics, 193(2), 357–397. https://doi.org/10.1016/j.jcp.2003.08.010
4. Nocedal, J., & Wright, S. J. (2006). Numerical Optimization (2nd ed.). Springer. https://doi.org/10.1007/978-0-387-40065-5