Gradient Descent

Gradient descent minimizes a differentiable function by starting from a guess and repeatedly stepping in the direction in which the function decreases fastest, the negative gradient. It needs only first derivatives and scales to millions of unknowns. In computer vision it minimizes variational energies for denoising and optical flow, fits models to pixel data, and, in stochastic form, trains neural networks.

Augustin-Louis Cauchy described the method in an 1847 note, motivated by solving for the six elements of an orbit without eliminating variables. As Lemaréchal (2012) recounts, Cauchy chose the step either small enough to decrease the function or by minimizing along the descent line, reduced a system of equations to minimizing their sum of squares, and argued convergence only informally.

“Gradient” here means the derivative of an objective with respect to the parameters, not the image gradient, the spatial derivative of intensity.

Problem

Find a minimizer $\boldsymbol{\theta}^\star = \arg\min_{\boldsymbol{\theta}} f(\boldsymbol{\theta})$ of a differentiable $f : \mathbb{R}^n \to \mathbb{R}$ when no closed form exists or it is too expensive. Linear least squares has a closed form, but factoring its normal equations for an image with a million unknowns is expensive; a gradient costs about as much as evaluating $f$.

Inputs and Outputs

• Input: $f$ and its gradient $\nabla f$, a starting point $\boldsymbol{\theta}_0$, a step-size rule, and a stopping criterion, usually a small gradient norm or an iteration budget.
• Output: a point whose gradient is approximately zero: an approximate global minimizer if $f$ is convex, an approximate stationary point otherwise.

Intuition

On a foggy hillside, the only information is the slope underfoot. Walking a short way downhill and re-reading the slope eventually reaches a valley floor. Too long a stride overshoots; in a long, narrow valley the steepest direction points mostly across it, so the walker zig-zags between the walls.

Algorithm

1. Choose a starting point $\boldsymbol{\theta}_0$.
2. Compute $\mathbf{g}_k = \nabla f(\boldsymbol{\theta}_k)$. Stop if $\lVert \mathbf{g}_k \rVert$ is below a tolerance.
3. Choose a step size $\eta_k > 0$, fixed or by a line search. In machine learning it is called the learning rate.
4. Set $\boldsymbol{\theta}_{k+1} = \boldsymbol{\theta}_k - \eta_k \mathbf{g}_k$ and return to step 2.

Mathematical Formulation

Why the negative gradient

To first order, $f(\boldsymbol{\theta} + \mathbf{v}) \approx f(\boldsymbol{\theta}) + \nabla f(\boldsymbol{\theta})^\top \mathbf{v}$. Among steps of fixed Euclidean length, the Cauchy–Schwarz inequality makes $\nabla f^\top \mathbf{v}$ most negative when $\mathbf{v}$ points along $-\nabla f$: the negative gradient is the steepest-descent direction for the Euclidean norm (Boyd and Vandenberghe, §9.4.1).

Step size on a quadratic

Linear least squares and the quadratic energies of variational vision have the form $f(\boldsymbol{\theta}) = \tfrac{1}{2}\boldsymbol{\theta}^\top H \boldsymbol{\theta} - \mathbf{b}^\top \boldsymbol{\theta}$ with a symmetric positive definite Hessian $H$ ($H = A^\top A$ for least squares). Since $\nabla f = H\boldsymbol{\theta} - \mathbf{b}$, the error $\mathbf{e}_k = \boldsymbol{\theta}_k - \boldsymbol{\theta}^\star$ evolves as

$$\mathbf{e}_{k+1} = (I - \eta H)\,\mathbf{e}_k ,$$

so along an eigenvector of $H$ with eigenvalue $\mu_i$ the error is multiplied by $1 - \eta\mu_i$ per step. With $m$ and $L$ the smallest and largest eigenvalues:

• If $0 < \eta < 2/L$, every factor has magnitude below one and the iteration converges.
• If $\eta > 2/L$, the component along the top eigenvector grows by $|1 - \eta L| > 1$ per step, flipping sign each time: the iteration diverges.
• The slowest factor is $\max(|1 - \eta m|, |1 - \eta L|)$. With $\eta = 1/L$ it is $1 - 1/\kappa$, where $\kappa = L/m$ is the condition number of $H$; the best constant step, $2/(m + L)$, makes the two terms equal and gives $(\kappa - 1)/(\kappa + 1)$. Either way, the number of steps needed grows in proportion to $\kappa$.

Convergence rates

For general $f$, $L$ is a Lipschitz constant of the gradient, $\lVert \nabla f(\mathbf{x}) - \nabla f(\mathbf{y}) \rVert \le L \lVert \mathbf{x} - \mathbf{y} \rVert$, and $f$ is called $L$-smooth.

• Convex and $L$-smooth. With $\eta = 1/L$ and iterations numbered from 1, $f(\boldsymbol{\theta}_t) - f^\star \le 2L\lVert \boldsymbol{\theta}_1 - \boldsymbol{\theta}^\star \rVert^2/(t - 1)$ (Bubeck, Theorem 3.3): the error decays as $O(1/t)$.
• Strongly convex and $L$-smooth. If the Hessian’s eigenvalues are also bounded below by $m > 0$, convergence is linear, a constant factor per step. With exact line search, $f(\boldsymbol{\theta}_k) - f^\star \le (1 - m/L)^k\,(f(\boldsymbol{\theta}_0) - f^\star)$, with a similar bound for backtracking (Boyd and Vandenberghe, §9.3.1); for $\eta = 1/L$, $\lVert \boldsymbol{\theta}_{t+1} - \boldsymbol{\theta}^\star \rVert^2 \le e^{-t/\kappa} \lVert \boldsymbol{\theta}_1 - \boldsymbol{\theta}^\star \rVert^2$ (Bubeck, Theorem 3.10).

The constant matters: from their examples, Boyd and Vandenberghe conclude that convergence can be very slow with condition numbers in the hundreds, and that at 1000 or more gradient descent is useless in practice (§9.3.2).

Line search

Since $L$ and $m$ are rarely known, practical codes choose $\eta_k$ at each step. Exact line search minimizes $f$ along the ray $\boldsymbol{\theta}_k - \eta \mathbf{g}_k$, which pays off only when that one-dimensional problem is cheap. Backtracking line search starts from a trial step and multiplies it by $\beta \in (0, 1)$ until the sufficient decrease (Armijo) condition holds,

$$f(\boldsymbol{\theta}_k - \eta\,\mathbf{g}_k) \le f(\boldsymbol{\theta}_k) - \alpha\,\eta\,\lVert \mathbf{g}_k \rVert^2, \qquad 0 < \alpha < \tfrac{1}{2},$$

that is, until the step achieves a fraction $\alpha$ of the decrease the linear model predicts (Boyd and Vandenberghe, §9.2). Larry Armijo (1966) proved convergence for a rule of this kind, halving the step until the decrease is at least half the linear prediction, assuming among other conditions a Lipschitz gradient and a unique minimizer that is the only stationary point. Since every accepted step decreases $f$, backtracking cannot diverge, and it needs no knowledge of $L$.

Parameters

• Step size $\eta$. For a fixed step, $1/L$ is the safe default; on a quadratic, a step above $2/L$ diverges from almost every start, and a tiny step is proportionally slow.
• Backtracking constants. Boyd and Vandenberghe give typical values of $\alpha$ from 0.01 to 0.3 and of $\beta$ from 0.1 to 0.8 (§9.2); in their examples the choice had a noticeable but not dramatic effect (§9.3.2).
• Starting point. Irrelevant to the limit for strictly convex $f$, decisive for nonconvex $f$; in image restoration the observed image is the natural start.

Complexity

Each iteration costs one gradient, one function evaluation per line-search trial, and $O(n)$ memory. For least squares with an $m \times n$ matrix, a gradient costs $O(mn)$, against $O(mn^2)$ for a direct solve; for image energies built from local differences it is proportional to the number of pixels.

Implementation

The first example denoises a $32 \times 32$ synthetic image $\mathbf{y}$ by minimizing a quadratic energy, the discrete analogue of the regularized energies of variational image denoising and Horn–Schunck flow:

$$f(\mathbf{u}) = \tfrac{1}{2}\lVert \mathbf{u} - \mathbf{y} \rVert^2 + \tfrac{\lambda}{2}\lVert D\mathbf{u} \rVert^2, \qquad \nabla f(\mathbf{u}) = (\mathbf{u} - \mathbf{y}) + \lambda D^\top D\,\mathbf{u},$$

where $D$ stacks horizontal and vertical forward differences. By Gershgorin’s theorem the eigenvalues of $D^\top D$ are at most 8, and a constant image has $D\mathbf{u} = \mathbf{0}$, so the Hessian’s eigenvalues lie in $[1, 1 + 8\lambda]$; we take $L = 1 + 8\lambda$, and $\kappa$ is just below it. The gradient was checked against finite differences.

import numpy as np

rng = np.random.default_rng(0)

# 32x32 synthetic image: a bright square and a darker disk, plus Gaussian noise.
n = 32
yy, xx = np.mgrid[0:n, 0:n]
clean = np.zeros((n, n))
clean[6:18, 6:18] = 1.0
clean[(yy - 22) ** 2 + (xx - 21) ** 2 < 49] = 0.6
noisy = clean + rng.normal(0, 0.2, size=clean.shape)

def dx(u): return np.diff(u, axis=1, append=u[:, -1:])   # forward differences,
def dy(u): return np.diff(u, axis=0, append=u[-1:, :])   # zero at the far border
def dxT(p): q = -p.copy(); q[:, 1:] += p[:, :-1]; return q   # adjoints
def dyT(p): q = -p.copy(); q[1:, :] += p[:-1, :]; return q

def f(u, lam):      # 1/2 ||u - noisy||^2 + lam/2 ||D u||^2
    return 0.5 * np.sum((u - noisy) ** 2) + 0.5 * lam * np.sum(dx(u) ** 2 + dy(u) ** 2)

def grad_f(u, lam):
    return (u - noisy) + lam * (dxT(dx(u)) + dyT(dy(u)))

def gd(lam, eta, iters):
    u = noisy.copy()
    for _ in range(iters):
        u = u - eta * grad_f(u, lam)
    return u

def gd_armijo(lam, iters, alpha=0.3, beta=0.5):
    u, halvings = noisy.copy(), 0
    for _ in range(iters):
        g, t = grad_f(u, lam), 1.0
        while f(u - t * g, lam) > f(u, lam) - alpha * t * np.sum(g * g):
            t *= beta
            halvings += 1
        u = u - t * g
    return u, halvings

lam = 1.0
L = 1 + 8 * lam                        # bound on the largest Hessian eigenvalue
f_star = f(gd(lam, 1 / L, 5000), lam)  # reference minimum
print(f"lam = {lam}, L = {L}: f(noisy) - f* = {f(noisy, lam) - f_star:.2e}")
for c in [0.1, 1.0, 1.9, 2.1]:
    u = gd(lam, c / L, 30)
    print(f"  eta = {c:3.1f}/L, 30 steps: f - f* = {f(u, lam) - f_star:.2e}")
u, h = gd_armijo(lam, 30)
print(f"  Armijo backtracking, 30 steps: f - f* = {f(u, lam) - f_star:.2e} ({h} halvings)")
rms = lambda u: np.sqrt(np.mean((u - clean) ** 2))
print(f"  RMS error: noisy {rms(noisy):.3f}, denoised {rms(gd(lam, 1 / L, 5000)):.3f}")

# Stronger smoothing makes the problem worse conditioned (kappa = 1 + 8 lam).
for lam in [1.0, 10.0, 100.0]:
    L, ref = 1 + 8 * lam, gd(lam, 1 / (1 + 8 * lam), 50000)
    u, k = noisy.copy(), 0
    while f(u, lam) - f(ref, lam) > 1e-6 * (f(noisy, lam) - f(ref, lam)):
        u, k = u - grad_f(u, lam) / L, k + 1
    print(f"lam = {lam:5.1f}: steps of size 1/L to reduce f - f* by 1e6: {k}")

Output:

lam = 1.0, L = 9.0: f(noisy) - f* = 7.77e+01
  eta = 0.1/L, 30 steps: f - f* = 3.83e+00
  eta = 1.0/L, 30 steps: f - f* = 1.65e-04
  eta = 1.9/L, 30 steps: f - f* = 1.15e-03
  eta = 2.1/L, 30 steps: f - f* = 2.14e+02
  Armijo backtracking, 30 steps: f - f* = 1.40e-11 (63 halvings)
  RMS error: noisy 0.195, denoised 0.122
lam =   1.0: steps of size 1/L to reduce f - f* by 1e6: 33
lam =  10.0: steps of size 1/L to reduce f - f* by 1e6: 236
lam = 100.0: steps of size 1/L to reduce f - f* by 1e6: 886

A step of $0.1/L$ is still far from the minimum, $1/L$ and $1.9/L$ converge, and $2.1/L$ diverges, ending above the noisy image’s energy. In this run backtracking accepted steps mostly alternating between $4.5/L$ and $1.125/L$; steps longer than $2/L$ are allowed because the Armijo rule only requires each step to decrease $f$ enough, and here the alternation converged much faster than a constant $1/L$.

Raising $\lambda$ from 1 to 100 raises $\kappa$ from about 9 to about 800, and the step count from 33 to 886. Along the slowest eigenvector the error in $f$ shrinks by $(1 - 1/\kappa)^2$ per step, so the worst case is about $(\kappa/2)\ln 10^6 \approx 5500$ steps. Our reading of the error recursion explains the gap: the initial error $\mathbf{y} - \mathbf{u}^\star$ is mostly high-frequency noise, along large-eigenvalue directions that are damped quickly, while the smooth, slowly converging components start close to the minimizer.

The second example fits a linear least-squares classifier, $\min_{\mathbf{w}} \tfrac{1}{2N}\lVert A\mathbf{w} - \mathbf{y} \rVert^2$, to the 357 images of 3s and 8s in scikit-learn’s 8×8 digit set, with labels $\pm 1$ and one column of $A$ per pixel plus a bias column. The four pixel scalings describe the same set of classifiers, so they share $f^\star$; only the optimization differs.

import numpy as np
from sklearn.datasets import load_digits

digits = load_digits()                   # 1797 images of 8x8 pixels, values 0..16
keep = np.isin(digits.target, [3, 8])
X = digits.data[keep]                    # one row of 64 pixel values per image
X = X[:, X.std(axis=0) > 0]              # drop pixels that are 0 in every image
y = np.where(digits.target[keep] == 3, 1.0, -1.0)
N = len(y)

def design(Xs):                          # pixels plus a constant column for the bias
    return np.column_stack([Xs, np.ones(N)])

def relative_gap(Xs, eta, iters=1000):
    """Run gradient descent on f(w) = 1/(2N) ||A w - y||^2 and return
    (f(w_k) - f*) / (f(0) - f*) after `iters` steps."""
    A = design(Xs)
    f = lambda w: 0.5 * np.mean((A @ w - y) ** 2)
    f_star = f(np.linalg.lstsq(A, y, rcond=None)[0])
    w = np.zeros(A.shape[1])
    f0 = f(w)
    for _ in range(iters):
        w = w - eta * A.T @ (A @ w - y) / N
    with np.errstate(over="ignore", invalid="ignore"):
        return (f(w) - f_star) / (f0 - f_star)

def L_and_kappa(Xs):
    e = np.linalg.eigvalsh(design(Xs).T @ design(Xs) / N)   # Hessian eigenvalues
    return e[-1], e[-1] / e[0]

scalings = {
    "raw, 0..16": X,
    "divided by 16": X / 16,
    "centered, / 16": (X - X.mean(axis=0)) / 16,
    "standardized": (X - X.mean(axis=0)) / X.std(axis=0),
}
for name, Xs in scalings.items():
    L, kappa = L_and_kappa(Xs)
    print(f"{name:15s} L = {L:7.1f}  kappa = {kappa:9.0f}  "
          f"relative gap after 1000 steps of 1/L: {relative_gap(Xs, 1 / L):.1e}")

# A step size that is safe for one scaling diverges for another.
eta = 1 / L_and_kappa(X / 16)[0]
print(f"1/L of 'divided by 16' used on raw pixels: {relative_gap(X, eta, iters=20):.1e} after 20 steps")

Output:

raw, 0..16      L =  3059.0  kappa =   1775492  relative gap after 1000 steps of 1/L: 9.9e-03
divided by 16   L =    12.9  kappa =   1920151  relative gap after 1000 steps of 1/L: 1.0e-02
centered, / 16  L =     1.0  kappa =    148445  relative gap after 1000 steps of 1/L: 2.6e-03
standardized    L =     8.6  kappa =       229  relative gap after 1000 steps of 1/L: 1.7e-07
1/L of 'divided by 16' used on raw pixels: 1.0e+91 after 20 steps

Dividing the pixels by 16 shrinks $L$ about 240-fold but leaves $\kappa$, and so the progress per step at $\eta = 1/L$, almost unchanged; its practical effect is that a learning rate tuned for one scale diverges on the other. Centering each pixel reduces $\kappa$ about thirteenfold, and standardizing brings it to 229, with a gap of about $10^{-7}$ instead of $10^{-2}$. Border pixels that are almost always dark have tiny variance and stretch the cost surface along their weights; standardizing removes the stretch.

Properties and Behavior

• Zig-zagging. With exact line search each step is orthogonal to the previous one, since the new gradient has no component along the line just minimized. On $\tfrac{1}{2}(x_1^2 + \gamma x_2^2)$ started from $(\gamma, 1)$, $f$ shrinks by exactly $\big((\gamma - 1)/(\gamma + 1)\big)^2$ per step (Boyd and Vandenberghe, §9.3.2).
• Scaling as preconditioning. Steepest descent in the norm $\lVert \mathbf{v} \rVert_P = \sqrt{\mathbf{v}^\top P \mathbf{v}}$ steps along $-P^{-1}\nabla f$, which is plain gradient descent after the change of variables $\bar{\boldsymbol{\theta}} = P^{1/2}\boldsymbol{\theta}$; a $P$ close to the Hessian makes the transformed problem well conditioned (Boyd and Vandenberghe, §9.4.1, §9.4.4). Standardizing pixels is a cheap diagonal version.
• Nonconvex objectives. For $L$-smooth $f$, a step of $1/L$ decreases $f$ by at least $\lVert \nabla f \rVert^2/(2L)$ (Bubeck, inequality 3.5; its derivation uses only smoothness). If $f$ is bounded below, the gradient therefore becomes small, but only near a stationary point: a local minimum, a saddle point, or even a maximum. Lee and colleagues (2016) proved that if every critical point is a local minimum or a strict saddle (a Hessian with a negative eigenvalue), gradient descent with a constant step below $1/L$ from a random start almost surely does not converge to a saddle point. Du and colleagues (2017) constructed smooth functions on which, from a natural random start, it needs exponentially many steps in the dimension to escape a chain of saddle points, while a randomly perturbed version needs polynomially many.

Limitations

• Slow on ill-conditioned problems. Poorly scaled pixels, strong regularization, or nearly collinear features make $\kappa$, and the step count, large.
• First-order only. For sums of squared residuals, the Gauss–Newton and Levenberg–Marquardt methods of nonlinear least squares use curvature and usually need far fewer steps.
• Full gradients over millions of training images are expensive.

Variants

• Preconditioned gradient descent uses the scaled step above, often with a diagonal $P$.
• Momentum methods add a fraction of the previous step to each update. Nesterov’s accelerated gradient is one: for smooth convex $f$ it improves the $O(1/t)$ rate to $O(1/t^2)$, and in the strongly convex case its step count grows with $\sqrt{\kappa}$ instead of $\kappa$ (Bubeck, §3.7).
• Projected gradient descent handles a convex constraint set, such as pixel values in $[0, 1]$, by projecting each iterate onto it (Bubeck, §3.1).
• Stochastic gradient descent estimates the gradient from a random mini-batch of images; with variants such as Adam, it trains neural networks, whose gradients come from backpropagation.
• Coarse-to-fine optimization minimizes an image energy on an image pyramid, initializing each level from the coarser one.