Week 8: Newton and Quasi-Newton

Newton Step#

Gradient descent uses a first-order Taylor expansion: $f(x + \Delta x) \approx f(x) + \nabla f(x)^T \Delta x$. The model is linear in $\Delta x$; the step direction is $-\nabla f(x)$, independent of curvature.

Newton's method uses a second-order Taylor expansion:

$f(x + \Delta x) \approx f(x) + \nabla f(x)^T \Delta x + \tfrac{1}{2} \Delta x^T \nabla^2 f(x) \Delta x$

Minimising this quadratic model in $\Delta x$ (setting derivative to zero) gives the Newton step:

$\nabla^2 f(x) \, \Delta x_{\text{nt}} = -\nabla f(x)$

The Newton step $\Delta x_{\text{nt}}$ is the solution to an $n \times n$ linear system — it incorporates both the gradient direction and the local curvature (Hessian). When $\nabla^2 f(x) \succ 0$, the step is guaranteed to be a descent direction.

Worked 1D example. $f(x) = \tfrac{1}{2}a x^2$ with $a > 0$. The Newton step from $x_0$ solves $f''(x_0) \Delta x = -f'(x_0)$, i.e., $a \cdot \Delta x = -a x_0$, giving $\Delta x = -x_0$. Hence $x_1 = x_0 + \Delta x = 0$ — the exact minimiser in one step. For quadratic functions with $a > 0$, Newton converges in one iteration from any starting point.

Worked 2D example. $f(x, y) = x^2 + 10y^2$. The Hessian is $\nabla^2 f = \begin{pmatrix} 2 & 0 \\ 0 & 20 \end{pmatrix}$. The condition number is $20/2 = 10$.

Starting from $(10, 1)$, the gradient is $(20, 20)$ and GD with $\eta = 0.05$ takes many iterations. Newton: $\Delta x = -(2)^{-1} \cdot 20 = -10$, $\Delta y = -(20)^{-1} \cdot 20 = -1$, so $x_1 = (0, 0)$ — exact in one step. Newton is scale-invariant: the Hessian automatically corrects for differing curvatures along each coordinate.

Quadratic Convergence Intuition#

Newton exhibits quadratic convergence near a minimiser: there exists a constant $M$ such that

$\|x_{k+1} - x^*\| \leq M \|x_k - x^*\|^2$

The error squares at each iteration. If $\|x_k - x^*\| = 10^{-2}$, the next iterate has error $\approx 10^{-4}$, then $10^{-8}$, and so on — doubling the number of correct digits per step. Gradient descent, by contrast, is linear: $\|x_{k+1} - x^*\| \leq \rho \|x_k - x^*\|$ with $\rho < 1$, gaining a fixed number of digits per step.

Requirements for quadratic convergence:

• $\nabla^2 f$ is Lipschitz continuous near $x^*$
• $\nabla^2 f(x^*) \succ 0$ (the optimum is non-degenerate)
• The initial point is "close enough" to $x^*$

The last condition is the catch. Far from the optimum, the Hessian may be indefinite or poorly conditioned, and the pure Newton step can increase the objective. This is why damping and line search are essential in practice.

Damping and Line Search#

Damped Newton modifies the pure step with a line search:

$x_{k+1} = x_k + \alpha_k \Delta x_{\text{nt}}, \quad \alpha_k \in (0, 1]$

where $\alpha_k$ is chosen via backtracking to satisfy the Armijo condition:

$f(x_k + \alpha \Delta x_{\text{nt}}) \leq f(x_k) + c \alpha \nabla f(x_k)^T \Delta x_{\text{nt}}$

with $c \in (0, 0.5)$. When $\nabla^2 f(x_k) \succ 0$, $\nabla f(x_k)^T \Delta x_{\text{nt}} = -\nabla f(x_k)^T (\nabla^2 f(x_k))^{-1} \nabla f(x_k) < 0$, so the directional derivative is negative — a descent direction exists, and backtracking finds a stepsize. When close to $x^*$, $\alpha_k = 1$ is accepted (full Newton step), and the quadratic convergence kicks in.

When the Hessian is not PD. If $\nabla^2 f(x_k)$ has negative or zero eigenvalues, the Newton direction may point uphill. Practical Newton solvers modify the Hessian:

• Add a multiple of identity: $\nabla^2 f + \mu I$ (Levenberg–Marquardt).
• Use a modified Cholesky that perturbs negative eigenvalues to positive.
• Fall back to gradient descent when the Newton decrement is too large.

BFGS and L-BFGS#

The Hessian is expensive: storing $n^2$ entries and solving an $n \times n$ linear system costs $O(n^3)$ per iteration. Quasi-Newton methods avoid the Hessian entirely, building an approximation $B_k \approx \nabla^2 f(x_k)$ from gradient changes observed during optimisation.

BFGSBroyden–Fletcher–Goldfarb–Shanno (quasi-Newton) update. Given $s_k = x_{k+1} - x_k$ (step) and $y_k = \nabla f(x_{k+1}) - \nabla f(x_k)$ (gradient change), the BFGS formula updates the inverse Hessian approximation $H_k \approx (\nabla^2 f(x_k))^{-1}$:

$H_{k+1} = \left(I - \frac{s_k y_k^T}{y_k^T s_k}\right) H_k \left(I - \frac{y_k s_k^T}{y_k^T s_k}\right) + \frac{s_k s_k^T}{y_k^T s_k}$

This is a rank-2 update — computationally cheap ($O(n^2)$). The direction is $d_k = -H_k \nabla f(x_k)$. BFGS preserves positive definiteness of $H_k$ when $y_k^T s_k > 0$ (the curvature condition), which holds under convexity or with a suitable line search.

L-BFGS (Limited-memory BFGS). Storing the $n \times n$ matrix $H_k$ is still too expensive for large $n$. L-BFGS stores only the last $m$ pairs $(s_k, y_k)$ (typically $m = 10$–$20$), and computes $H_k \nabla f(x_k)$ implicitly via a two-loop recursion in $O(mn)$ time and $O(mn)$ memory. For $n = 10^7$ (a large deep net), $m = 10$ pairs cost $10^8$ doubles (~800 MB) — manageable; $n^2 = 10^{14}$ is not.

vs Adam#

| Property | Newton / L-BFGS | AdamAdaptive Moment Estimation (optimizer) | |---|---|---| | Uses curvature | Full or approximated Hessian | Diagonal (element-wise) second-moment estimate | | Per-iteration cost | $O(n^2)$ (BFGS) or $O(mn)$ (L-BFGS) | $O(n)$ | | Convergence near optimum | Quadratic (Newton) or superlinear (BFGS) | Linear (no acceleration without momentum tuning) | | Works well when | $n \lesssim 10^4$ (Newton) or $n \lesssim 10^7$ (L-BFGS), accurate gradients | $n$ is huge, gradients are stochastic, nonconvex | | Sensitive to | Ill-conditioning, nondifferentiability | Learning rate, $\beta_1, \beta_2$ |

In practice, L-BFGS is competitive with Adam for moderate-scale problems where gradients are computed accurately (full batch or large minibatch). Deep nets with stochastic gradients and $n > 10^7$ still favour AdamAdaptive Moment Estimation (optimizer) or SGDStochastic Gradient Descent with momentum because L-BFGS's curvature approximation degrades under gradient noise. However, for scientific computing, logistic regression, and moderate-scale ML, L-BFGS is often the faster and more reliable choice.

Browser lab#

Compare Newton's method vs gradient descent on a well-conditioned quadratic ($f_1$) and an ill-conditioned one ($f_2$). Count iterations to reach a tolerance.

import numpy as np
import matplotlib.pyplot as plt

def newton_quadratic(Q, x0, tol=1e-8, max_iter=50):
    """Newton on f(x) = 0.5 * x^T Q x."""
    x = x0.copy().astype(float)
    path = [x.copy()]
    for _ in range(max_iter):
        g = Q @ x
        if np.linalg.norm(g) < tol:
            break
        dx = np.linalg.solve(Q, -g)
        x = x + dx
        path.append(x.copy())
    return np.array(path)

def gd_quadratic(Q, x0, eta, tol=1e-8, max_iter=500):
    """Gradient descent on f(x) = 0.5 * x^T Q x."""
    x = x0.copy().astype(float)
    path = [x.copy()]
    for _ in range(max_iter):
        g = Q @ x
        if np.linalg.norm(g) < tol:
            break
        x = x - eta * g
        path.append(x.copy())
    return np.array(path)

# Well-conditioned: Q1 = diag(1, 1)
Q1 = np.diag([1.0, 1.0])
# Ill-conditioned: Q2 = diag(1, 100)
Q2 = np.diag([1.0, 100.0])

x0 = np.array([10.0, 1.0])

# Newton
path_n1 = newton_quadratic(Q1, x0)
path_n2 = newton_quadratic(Q2, x0)

# GD with optimal step size for quadratics
eta1 = 2.0 / (1.0 + 1.0)  # 2/(L+mu)
eta2 = 2.0 / (1.0 + 100.0)

path_g1 = gd_quadratic(Q1, x0, eta1)
path_g2 = gd_quadratic(Q2, x0, eta2)

print(f"Well-conditioned: Newton {len(path_n1)-1} iters, GD {len(path_g1)-1} iters")
print(f"Ill-conditioned: Newton {len(path_n2)-1} iters, GD {len(path_g2)-1} iters")

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))

def plot_path(ax, paths, labels, colors, title):
    for path, lbl, clr in zip(paths, labels, colors):
        ax.plot(path[:, 0], path[:, 1], f"{clr}o-", markersize=3, linewidth=1, label=lbl)
    ax.plot(0, 0, "kx", markersize=10)
    ax.set_xlabel("x1")
    ax.set_ylabel("x2")
    ax.set_title(title)
    ax.legend()
    ax.set_aspect("equal")

plot_path(ax1, [path_n1, path_g1], [f"Newton ({len(path_n1)-1})", f"GD ({len(path_g1)-1})"],
          ["b", "r"], "Well-conditioned: $Q = \mathrm{diag}(1,1)$")
plot_path(ax2, [path_n2, path_g2], [f"Newton ({len(path_n2)-1})", f"GD ({len(path_g2)-1})"],
          ["b", "r"], "Ill-conditioned: $Q = \mathrm{diag}(1,100)$")

plt.tight_layout()
plt.show()