Newton, quasi-Newton and second-order methods

Gradient descent uses only the slope of the objective, so on ill-conditioned problems it zigzags and needs thousands of steps (Lesson 1). Second-order methods also use curvature, the Hessian or an approximation to it, to take steps scaled to the shape of the objective. Newton's method converges in a handful of iterations near a solution; quasi-Newton methods such as BFGS and L-BFGS approach that speed using only gradients; Gauss-Newton and Levenberg-Marquardt exploit the structure of least-squares problems; and for large models, Hessian-vector products and structured approximations such as K-FAC and Shampoo bring curvature information within reach. This lesson explains these methods, why they are fast, what they cost, and where each is used in machine learning.

Newton's method

Newton's method minimizes the local quadratic model of \(f\) at each step:

On a quadratic it reaches the minimum in one step, whatever the conditioning, because the Hessian exactly undoes the different curvatures in different directions. On a general smooth, strongly convex function, it converges quadratically once close to the solution: the number of correct digits roughly doubles each iteration.

import numpy as np

rng = np.random.default_rng(2)
n, d = 1000, 20
X = rng.normal(size=(n, d)) * np.linspace(0.3, 3.0, d)       # badly scaled features
w_true = rng.normal(size=d)
y = (rng.random(n) < 1 / (1 + np.exp(-X @ w_true))).astype(float)
lam = 1e-3

def loss_grad_hess(w):
    z = X @ w
    p = 1 / (1 + np.exp(-z))
    loss = np.mean(np.logaddexp(0, z) - y * z) + lam / 2 * w @ w
    grad = X.T @ (p - y) / n + lam * w
    hess = (X.T * (p * (1 - p))) @ X / n + lam * np.eye(d)
    return loss, grad, hess

# Reference optimum from many Newton steps
w_star = np.zeros(d)
for _ in range(50):
    _, g, H = loss_grad_hess(w_star)
    w_star -= np.linalg.solve(H, g)
f_star = loss_grad_hess(w_star)[0]

# Newton method
w = np.zeros(d)
print("Newton")
for k in range(1, 8):
    _, g, H = loss_grad_hess(w)
    w = w - np.linalg.solve(H, g)
    print(f"  step {k}: f - f* = {loss_grad_hess(w)[0] - f_star:.2e}")

# Gradient descent with step 1/L
L = np.linalg.eigvalsh(X.T @ X / n).max() / 4 + lam
w = np.zeros(d)
for k in range(1, 5001):
    w -= loss_grad_hess(w)[1] / L
    if k in (10, 100, 1000, 5000):
        print(f"gradient descent step {k:5d}: f - f* = {loss_grad_hess(w)[0] - f_star:.2e}")

Newton
  step 1: f - f* = 1.69e-01
  step 2: f - f* = 6.94e-02
  step 3: f - f* = 2.28e-02
  step 4: f - f* = 4.36e-03
  step 5: f - f* = 2.58e-04
  step 6: f - f* = 1.23e-06
  step 7: f - f* = 3.15e-11
gradient descent step    10: f - f* = 8.86e-02
gradient descent step   100: f - f* = 8.81e-03
gradient descent step  1000: f - f* = 8.51e-06
gradient descent step  5000: f - f* = 3.09e-14

The output shows the two phases of Newton's method: a few damped-looking steps far from the solution, then quadratic convergence, with the error falling from \(10^{-4}\) to \(10^{-6}\) to \(10^{-11}\). Gradient descent makes steady linear progress, at a rate set by the condition number.

Making Newton robust

Pure Newton steps can fail far from a solution: if the Hessian is not positive definite (non-convex regions), the step may point uphill; if the quadratic model is poor, the step may overshoot. Practical implementations add:

• A line search on the Newton direction (damped Newton), accepting a step size \(t \le 1\) that decreases \(f\) sufficiently.
• Hessian modification: adding \(\lambda I\) to make it positive definite, or using a trust region that limits the step to where the model is trusted.

The cost

Each Newton step needs the Hessian (\(d^2\) entries) and a linear solve (about \(d^3/3\) operations by Cholesky factorization). For a logistic regression with a few thousand features this is cheap; for a neural network with \(10^9\) parameters, the Hessian alone would need \(10^{18}\) entries. Second-order methods for large problems must avoid forming the Hessian.

Quasi-Newton methods: BFGS and L-BFGS

Quasi-Newton methods build an approximation to the inverse Hessian from the gradients they have already computed. After each step \(s_k = x_{k+1} - x_k\) with gradient change \(y_k = \nabla f_{k+1} - \nabla f_k\), they update the approximation \(H_k\) so that it satisfies the secant condition \(H_{k+1}y_k = s_k\): it maps the observed change in gradient to the observed step, as the true inverse Hessian would for a quadratic. The BFGS update (Broyden, Fletcher, Goldfarb, Shanno) does this while keeping \(H\) symmetric and positive definite:

BFGS converges superlinearly near a solution, needing only gradients, but stores a dense \(d \times d\) matrix. Limited-memory BFGS (L-BFGS) stores only the last \(m\) pairs \((s_k, y_k)\), typically 5 to 20, and applies the approximate inverse Hessian with a two-loop recursion in \(O(md)\) operations. L-BFGS is the standard choice for smooth, deterministic problems of moderate to large size: logistic regression with many features, conditional random fields, physics-informed fitting, and the final refinement of some scientific models.

import numpy as np
from scipy.optimize import minimize, least_squares

# Quasi-Newton: L-BFGS needs only gradients, and few of them, on the Rosenbrock function
def rosen(x):
    return np.sum(100 * (x[1:] - x[:-1] ** 2) ** 2 + (1 - x[:-1]) ** 2)

def rosen_grad(x):
    g = np.zeros_like(x)
    g[:-1] = -400 * x[:-1] * (x[1:] - x[:-1] ** 2) - 2 * (1 - x[:-1])
    g[1:] += 200 * (x[1:] - x[:-1] ** 2)
    return g

x0 = np.full(50, -1.2)
for method in ("CG", "L-BFGS-B", "BFGS"):
    r = minimize(rosen, x0, jac=rosen_grad, method=method, options={"maxiter": 20000, "gtol": 1e-8})
    print(f"{method:9s}: iterations {r.nit:5d}, gradient evaluations {r.njev:5d}, f = {r.fun:.1e}")

# Gauss-Newton / Levenberg-Marquardt for nonlinear least squares: fit y = a exp(-b t) + c
rng = np.random.default_rng(4)
t = np.linspace(0, 10, 60)
y = 4.0 * np.exp(-0.7 * t) + 0.5 + rng.normal(0, 0.05, t.size)
def residuals(p):
    a, b, c = p
    return a * np.exp(-b * t) + c - y
fit = least_squares(residuals, x0=[1.0, 0.1, 0.0], method="lm")
print(f"Levenberg-Marquardt: a, b, c = {np.round(fit.x, 3)} after {fit.nfev} function evaluations")

CG       : iterations   374, gradient evaluations   625, f = 8.9e-17
L-BFGS-B : iterations   267, gradient evaluations   320, f = 9.6e-10
BFGS     : iterations   336, gradient evaluations   383, f = 1.3e-19
Levenberg-Marquardt: a, b, c = [4.004 0.7   0.495] after 8 function evaluations

Quasi-Newton methods need accurate gradients. With the noisy minibatch gradients of deep learning, the secant pairs are dominated by noise and the curvature estimates become unreliable, which is one reason deep networks are trained with first-order methods such as Adam (Lesson 5) rather than L-BFGS, although full-batch L-BFGS works well for small networks and for fine-tuning in some scientific applications.

Nonlinear least squares: Gauss-Newton and Levenberg-Marquardt

Many fitting problems minimize a sum of squared residuals \(f(x) = \frac{1}{2}\sum_i r_i(x)^2 = \frac{1}{2}\|r(x)\|^2\). The Hessian is \(J^TJ + \sum_i r_i\nabla^2 r_i\), where \(J\) is the Jacobian of the residuals. Near a good fit the residuals are small, so the second term can be dropped, giving the Gauss-Newton step:

which needs only first derivatives of the residuals. Levenberg-Marquardt adds a damping term, \((J^TJ + \lambda I)p = -J^Tr\), increasing \(\lambda\) when a step fails (moving toward gradient descent) and decreasing it when steps succeed (moving toward Gauss-Newton). It is the standard method for curve fitting, camera calibration, bundle adjustment in computer vision (AI-350), and parameter estimation in engineering models.

Curvature without the Hessian

For large models, Newton-type methods work with Hessian-vector products \(\nabla^2 f(x)\,v\), which cost about as much as one or two gradient evaluations: by automatic differentiation (Lesson 6, the "Pearlmutter trick"), or by finite differences of gradients. With them:

• Newton-CG (truncated Newton) solves the Newton system approximately with the conjugate gradient method, using only Hessian-vector products, and stops early.
• Power iteration or the Lanczos method estimates the largest Hessian eigenvalues, used to monitor sharpness during training (Lesson 9).

import numpy as np

# Newton-CG: solve H p = -g using only Hessian-vector products, never forming H
rng = np.random.default_rng(5)
n, d = 2000, 300
X = rng.normal(size=(n, d)) / np.sqrt(d)
y = (rng.random(n) < 0.5).astype(float)
lam = 1e-2
w = rng.normal(size=d) * 0.1

def grad(w):
    p = 1 / (1 + np.exp(-X @ w))
    return X.T @ (p - y) / n + lam * w

def hvp(w, v, eps=1e-5):
    # finite-difference Hessian-vector product: two gradient evaluations
    return (grad(w + eps * v) - grad(w - eps * v)) / (2 * eps)

def conjugate_gradient(Av, b, iters=50, tol=1e-10):
    x = np.zeros_like(b); r = b.copy(); p = r.copy(); rs = r @ r
    for k in range(iters):
        Ap = Av(p)
        alpha = rs / (p @ Ap)
        x += alpha * p; r -= alpha * Ap
        rs_new = r @ r
        if np.sqrt(rs_new) < tol:
            return x, k + 1
        p = r + rs_new / rs * p; rs = rs_new
    return x, iters

g = grad(w)
step, k = conjugate_gradient(lambda v: hvp(w, v), -g)
p = 1 / (1 + np.exp(-X @ w))
H = (X.T * (p * (1 - p))) @ X / n + lam * np.eye(d)
exact = np.linalg.solve(H, -g)
print(f"CG used {k} Hessian-vector products; relative difference from exact Newton step {np.linalg.norm(step - exact) / np.linalg.norm(exact):.1e}")
print(f"full Hessian for d = {d}: {d * d * 8 / 1e6:.2f} MB; for d = 1e9 parameters it would need {1e18 * 2 / 1e18:.0f} EB at 16-bit")

CG used 6 Hessian-vector products; relative difference from exact Newton step 7.0e-09
full Hessian for d = 300: 0.72 MB; for d = 1e9 parameters it would need 2 EB at 16-bit

Structured approximations for deep learning

Several methods approximate curvature with structures that are cheap to store and invert:

• Diagonal approximations: Adam and RMSProp (Lesson 5) can be viewed as using a diagonal preconditioner built from squared gradients.
• Natural gradient methods precondition with the Fisher information matrix, which for many losses equals the Gauss-Newton matrix; they make steps invariant to how the model is parameterized.
• K-FAC (Kronecker-factored approximate curvature) approximates each layer's Fisher block as a Kronecker product of two small matrices, one from the layer's inputs and one from its output gradients, which can be inverted cheaply.
• Shampoo and its successors keep a preconditioner for each dimension of each weight matrix (left and right factors) and apply their inverse roots. Variants have won optimizer benchmarks for training neural networks, with a modest extra cost per step that is repaid by fewer steps.

These methods occupy a middle ground: more expensive per step than Adam, far cheaper than Newton, and often faster in wall-clock time for large models when implemented efficiently and distributed across devices.

Choosing a method

• Small or medium smooth problems with exact gradients: Newton (with line search) or L-BFGS.
• Nonlinear least squares: Levenberg-Marquardt or trust-region Gauss-Newton.
• Large convex problems with structure (sparsity, constraints): proximal and first-order methods (Lesson 7), possibly with L-BFGS for smooth parts.
• Deep networks with minibatch gradients: Adam-family methods (Lesson 5) by default; structured second-order methods such as Shampoo when the infrastructure supports them.