16 Gradient Descent
Almost every model you will ever fit is trained by one idea: compute the gradient, take a small step against it, repeat.
It is barely more than Section 12.4 plus a loop. What fills this chapter is everything that goes wrong — steps that overshoot, valleys that cause zig-zagging, gradients too expensive to compute exactly — and the fixes that turned a simple idea into the workhorse of machine learning.
16.1 Descent directions
A direction \(\mathbf{d}\) is a descent direction at \(\mathbf{x}\) if moving a little way along it decreases \(f\). By Equation 12.3 that means
\[ \nabla f(\mathbf{x})^\top\mathbf{d} < 0 \]
any direction making an angle greater than \(90°\) with the gradient. There are many. The obvious choice is the steepest one, \(\mathbf{d} = -\nabla f\), which is where the method gets its other name, steepest descent.
It is worth noticing early that “steepest” is a local claim, and locally steepest is not the same as fastest overall. That gap is the subject of half this chapter.
16.2 The gradient descent algorithm
\[ \mathbf{x}_{k+1} = \mathbf{x}_k - \eta\,\nabla f(\mathbf{x}_k) \tag{16.1}\]
The scalar \(\eta > 0\) is the step size, or in machine learning the learning rate. That is the whole algorithm.
Start with the easy case: a circular bowl, where the gradient always points straight at the minimum.
Now the ill-conditioned valley from Figure 12.2, where \(\kappa = 10\).
[,1] [,2]
[1,] 1.0000 3.00000
[2,] -0.8000 2.46000
[3,] 0.6400 2.01720
[4,] -0.5120 1.65410
[5,] 0.4096 1.35637The \(x\) coordinate flips sign every step while shrinking; the \(y\) coordinate crawls down. Steepest descent is wasting most of its effort going sideways. This is the picture to keep — nearly every refinement below exists to fix it.
16.3 Choosing a step size
The single most consequential hyperparameter, and the failure modes at each end are completely different.
For a quadratic with curvature \(L\) in the steepest direction, the update multiplies each component by \((1 - \eta\lambda_i)\). Stability requires \(|1-\eta\lambda_i| < 1\) for every eigenvalue, which gives
\[ \eta < \frac{2}{L} \tag{16.2}\]
Above that, the largest-curvature direction is amplified rather than damped and the whole thing diverges. Here \(L = 20\), so the bound is \(\eta < 0.1\).
[,1] [,2] [,3]
[1,] 19.000000 1.900000e+01 1.900000e+01
[2,] 8.427397 2.554182e+00 6.266758e+01
[3,] 4.012633 6.205068e-03 1.469772e+04
[4,] 1.787840 2.026467e-06 2.160228e+07Three regimes, in three columns:
| \(\eta\) | Behavior |
|---|---|
| \(0.01\) — too small | stable, but the flat direction barely moves |
| \(0.0909\) — optimal | fastest possible for this problem |
| \(0.11\) — too large | past \(2/L\); the loss grows without bound |
The optimal step size for a quadratic is
\[ \eta^* = \frac{2}{L + m} \]
which here is \(2/22 \approx 0.0909\) — just under the stability bound, and not a coincidence: the best you can do is damp the stiffest direction as hard as stability allows.
Equation 16.2 depends on \(L\), the largest Hessian eigenvalue, which you almost never know for a real model. In practice \(\eta\) is tuned, and the symptom of exceeding the bound is unmistakable: the loss increases, often to NaN within a few steps.
If your loss explodes, reduce the learning rate first. It is by far the most common cause.
16.4 Convergence
For a strongly convex, smooth function (Section 15.7) with optimal step size, the distance to the minimum shrinks by a constant factor each step:
\[ \|\mathbf{x}_k - \mathbf{x}^*\| \leq \left(\frac{\kappa - 1}{\kappa + 1}\right)^{k}\|\mathbf{x}_0 - \mathbf{x}^*\| \tag{16.3}\]
This is linear convergence: error falls geometrically, and the ratio is set entirely by the condition number.
[,1] [,2] [,3] [,4] [,5]
kappa 1 2.0000 10.0000 100.0000 1000.000
rate 0 0.3333 0.8182 0.9802 0.998
steps_for_1e_3 1 7.0000 35.0000 346.0000 3454.000The last row is the headline. At \(\kappa = 1\) convergence is immediate; at \(\kappa = 10\) it takes about 35 steps to gain three digits; at \(\kappa = 1000\), roughly 3,500. Conditioning, not dimension, is what makes a problem slow.
This is why feature standardization helps so reliably. Rescaling variables to comparable ranges reduces \(\kappa\), and Equation 16.3 turns that directly into fewer iterations.
16.5 Line search
Rather than fixing \(\eta\), choose it at each step by asking how far to go along the current direction:
\[ \eta_k = \arg\min_{\eta > 0} f(\mathbf{x}_k - \eta\nabla f(\mathbf{x}_k)) \]
Exact line search solves that one-dimensional problem. It sounds ideal and is usually not worth it — each step costs many function evaluations, and on the valley it still zig-zags, because the direction was the problem, not the distance.
In practice a backtracking line search is used: start with a generous \(\eta\), halve it until the step actually decreases \(f\) by a reasonable amount.
backtrack <- function(f, grad, x, eta0 = 1, c1 = 1e-4) {
g <- grad(x)
eta <- eta0
while (f(x - eta * g) > f(x) - c1 * eta * sum(g^2)) {
eta <- eta / 2
}
eta
}
fv <- function(p) f_valley(p[1], p[2])
c(
at_start = backtrack(fv, grad_valley, c(1, 3)),
near_min = backtrack(fv, grad_valley, c(0.05, 0.2))
)at_start near_min
0.0625 0.0625 The condition being tested is the Armijo rule: accept a step only if it delivers at least a small fraction \(c_1\) of the decrease the gradient promised. It is cheap, needs no knowledge of \(L\), and is what most serious optimizers use.
Both calls return \(0.0625\), comfortably under the stability bound \(2/L = 0.1\) — found by halving from an optimistic guess, without ever computing an eigenvalue. That is the appeal: the rule discovers a safe step from function values alone.
16.6 Momentum
Here is the fix for the zig-zag. Instead of stepping along the current gradient, accumulate a running average of past gradients and step along that:
\[ \begin{aligned} \mathbf{v}_{k+1} &= \beta\mathbf{v}_k + \nabla f(\mathbf{x}_k) \\ \mathbf{x}_{k+1} &= \mathbf{x}_k - \eta\,\mathbf{v}_{k+1} \end{aligned} \tag{16.4}\]
with \(\beta \in [0,1)\), typically \(0.9\). Setting \(\beta = 0\) recovers plain gradient descent.
The intuition is physical, and it is the reason for the name. A ball rolling down the valley keeps its momentum along the floor while the side-to-side pushes, alternating in sign, cancel out. Oscillating components average away; consistent components accumulate.
momentum_descent <- function(grad, x0, eta, beta, n) {
path <- matrix(0, nrow = n + 1, ncol = length(x0))
path[1, ] <- x0
v <- rep(0, length(x0))
for (k in seq_len(n)) {
v <- beta * v + grad(path[k, ])
path[k + 1, ] <- path[k, ] - eta * v
}
path
}
p_plain <- momentum_descent(grad_valley, c(1, 3), 0.02, 0, 40)
p_mom <- momentum_descent(grad_valley, c(1, 3), 0.02, 0.9, 40)
rbind(
plain = f_valley(p_plain[, 1], p_plain[, 2])[c(11, 21, 41)],
momentum = f_valley(p_mom[, 1], p_mom[, 2])[c(11, 21, 41)]
) [,1] [,2] [,3]
plain 3.978388 1.758295 0.343511399
momentum 2.879627 1.225979 0.001741788After 40 steps momentum has reduced the loss roughly two hundredfold further, at identical cost per step.
16.7 Stochastic gradient descent
For a loss summed over \(n\) observations,
\[ L(\boldsymbol{\theta}) = \frac{1}{n}\sum_{i=1}^{n}\ell_i(\boldsymbol{\theta}) \]
the exact gradient costs a pass over the entire dataset. With millions of observations that is one step per pass — unaffordable.
Stochastic gradient descent uses one randomly chosen observation instead:
\[ \boldsymbol{\theta}_{k+1} = \boldsymbol{\theta}_k - \eta\,\nabla\ell_{i_k}(\boldsymbol{\theta}_k) \]
The gradient is now wrong, but right on average — its expectation is the true gradient — so the steps are noisy yet drift the right way. In exchange each step costs \(1/n\) as much.
set.seed(11)
n <- 500
xd <- rnorm(n)
yd <- 2 + 3 * xd + rnorm(n, sd = 0.5)
Xd <- cbind(1, xd)
grad_i <- function(theta, i) {
r <- drop(Xd[i, ] %*% theta) - yd[i]
2 * r * Xd[i, ]
}
theta <- c(0, 0)
for (k in 1:2000) {
i <- sample(n, 1)
theta <- theta - 0.01 * grad_i(theta, i)
}
rbind(
sgd = round(theta, 4),
exact = round(coef(lm(yd ~ xd)), 4)
) xd
sgd 2.0227 2.9980
exact 2.0102 3.0154Close to the least squares solution, having touched 2,000 single observations rather than making 2,000 full passes.
SGD does not converge to the minimum and stay there — it rattles around in a neighborhood whose size is set by \(\eta\) and the gradient noise. Getting a precise answer requires decaying the learning rate, and the classical conditions are \(\sum_k \eta_k = \infty\) with \(\sum_k \eta_k^2 < \infty\): large enough to travel any distance, small enough that the noise dies out. A schedule like \(\eta_k \propto 1/k\) satisfies both.
16.8 Mini-batches
Between one observation and all of them sits the practical choice: average the gradient over a mini-batch of \(B\) observations.
| Batch size | Gradient noise | Cost per step | Steps per epoch |
|---|---|---|---|
| \(B = 1\) | highest | lowest | \(n\) |
| \(1 < B < n\) | moderate | moderate | \(n/B\) |
| \(B = n\) (full batch) | none | highest | 1 |
Averaging \(B\) independent gradients reduces the noise standard deviation by \(\sqrt{B}\) — the same \(1/\sqrt{N}\) from Section 22.4. Note what that means: going from \(B=1\) to \(B=100\) costs a hundred times more per step and buys only a tenfold noise reduction. Diminishing returns are built into the arithmetic, which is why mini-batches are typically 32 to 512 rather than as large as memory allows.
Some gradient noise also appears to help, by discouraging the optimizer from settling into sharp narrow minima.
16.9 Adaptive methods
Plain gradient descent uses one \(\eta\) for every coordinate, which is a poor fit when they have wildly different curvatures — precisely the ill-conditioned case. Adaptive methods give each coordinate its own effective step size, scaled by the size of the gradients it has seen.
| Method | Idea |
|---|---|
| AdaGrad | divide by the square root of the accumulated squared gradient |
| RMSProp | same, but on an exponentially decaying average, so it does not stall |
| Adam | RMSProp plus momentum, with bias correction |
Adam is the default in deep learning. It is momentum and per-coordinate scaling combined — an approximation to using curvature information without ever forming a Hessian.
Adaptive methods are effectively estimating a diagonal approximation to the Hessian from the gradients they have already computed. That places them between plain gradient descent, which ignores curvature entirely, and Newton’s method, which uses all of it and cannot be afforded.
16.10 Newton’s method
Gradient descent uses a linear model of \(f\). Newton’s method uses a quadratic one — the second-order Taylor expansion from Equation 11.5 — and jumps straight to its minimum:
\[ \mathbf{x}_{k+1} = \mathbf{x}_k - \bigl[\nabla^2 f(\mathbf{x}_k)\bigr]^{-1}\nabla f(\mathbf{x}_k) \tag{16.5}\]
Multiplying by \(\mathbf{H}^{-1}\) undoes the stretch of the level sets: it turns an elongated valley back into a circular bowl, and then walks straight to the bottom.
On a quadratic that means it finishes in one step, whatever the condition number.
[,1] [,2]
start 1 3
after_one_newton_step 0 0Exactly the minimum, from a starting point that took gradient descent dozens of steps to approach.
The catch is cost. The Hessian has \(p^2\) entries and solving with it costs \(O(p^3)\) (Section 5.14), so for \(p\) in the millions Newton’s method is unthinkable. It also requires \(\mathbf{H}\) to be positive definite; near a saddle it can step uphill.
16.11 Quasi-Newton methods
The compromise: build an approximation to \(\mathbf{H}^{-1}\) from the gradients you are already computing, without ever forming the Hessian.
BFGS is the standard, updating its estimate after each step from the change in gradient. L-BFGS (“limited memory”) stores only the last few updates instead of a full \(p \times p\) matrix, which makes it usable at large scale. R’s optim() provides both.
x y value gradient_evals
0 0 0 9 Essentially exact in a handful of gradient evaluations.
| Method | Uses | Cost per step | Steps needed |
|---|---|---|---|
| Gradient descent | \(\nabla f\) | \(O(p)\) | many, set by \(\kappa\) |
| Momentum / Adam | \(\nabla f\) + history | \(O(p)\) | fewer |
| L-BFGS | \(\nabla f\) + recent history | \(O(mp)\) | far fewer |
| Newton | \(\nabla f\), \(\nabla^2 f\) | \(O(p^3)\) | very few |
The choice tracks problem size. For a few thousand parameters — a GLM, a small likelihood — L-BFGS is excellent and is what optim() and glm() reach for. For a neural network with millions of parameters, only the \(O(p)\) methods are viable, which is why Adam and SGD dominate despite needing far more steps. It is a straight trade of steps against cost per step.
16.12 Summary
| Idea | Statement |
|---|---|
| Update | \(\mathbf{x} \leftarrow \mathbf{x} - \eta\nabla f\) |
| Stability | \(\eta < 2/L\) |
| Optimal step (quadratic) | \(\eta^* = 2/(L+m)\) |
| Convergence rate | \((\kappa-1)/(\kappa+1)\) per step |
| Momentum | average past gradients; cancels oscillation |
| SGD | one observation per step; right on average |
| Mini-batch | noise falls as \(1/\sqrt{B}\) |
| Newton | \(\mathbf{x} \leftarrow \mathbf{x} - \mathbf{H}^{-1}\nabla f\); one step on a quadratic |
| In R | optim(..., method = "BFGS") |
16.13 Exercises
1. For \(f(x) = x^2\), write the gradient descent update and find the \(\eta\) that converges fastest.
\(f' = 2x\), so \(x_{k+1} = x_k - 2\eta x_k = (1-2\eta)x_k\). Convergence needs \(|1-2\eta| < 1\), i.e. \(0 < \eta < 1\), and the factor is smallest — zero — at \(\eta = 0.5\).
[1] 0.32768 0.00000 -0.32768 -2.48832At \(\eta = 0.5\) it lands exactly on the minimum in one step. That is Newton’s method in disguise: \(H^{-1} = 1/2\), so the optimal step for a one-dimensional quadratic is the Newton step.
2. Run gradient descent on the valley with \(\eta = 0.05\) and count steps to reach \(f < 10^{-3}\).
[1] 44Slower than the optimal \(\eta = 0.0909\) would be. With \(\eta = 0.05\) the \(x\) factor is \(|1-1| = 0\) — that direction is solved immediately — but the \(y\) factor is \(1 - 0.1 = 0.9\), and the flat direction now sets the pace.
A step size can be excellent for one direction and poor for another. That is precisely what adaptive methods try to fix.
3. Show that momentum with \(\beta = 0\) is plain gradient descent.
With \(\beta = 0\), Equation 16.4 gives \(\mathbf{v}_{k+1} = \nabla f(\mathbf{x}_k)\), so the update becomes \(\mathbf{x}_{k+1} = \mathbf{x}_k - \eta\nabla f(\mathbf{x}_k)\) — exactly Equation 16.1.
4. Verify that Newton’s method solves any quadratic in one step, using \(f(x,y) = 3x^2 + 5y^2 - 2x + 4y\).
\(\nabla f = (6x - 2,\; 10y + 4)\) and \(\mathbf{H} = \operatorname{diag}(6, 10)\), constant.
One step from an absurd starting point lands exactly on \((1/3, -2/5)\), and the gradient there is zero. Newton’s method is exact on a quadratic because the quadratic model is the function — no approximation error to iterate away.
5. Compare SGD against full-batch gradient descent on the regression data above, at equal numbers of observation touches. Run SGD from several seeds. What do you notice?
grad_full <- function(theta) {
r <- drop(Xd %*% theta) - yd
2 * drop(crossprod(Xd, r)) / n
}
th_full <- c(0, 0)
for (k in 1:20) th_full <- th_full - 0.1 * grad_full(th_full)
sgd_run <- function(seed, eta_fn) {
set.seed(seed)
th <- c(0, 0)
for (k in 1:10000) {
th <- th - eta_fn(k) * grad_i(th, sample(n, 1))
}
th
}
exact <- unname(coef(lm(yd ~ xd)))
err <- function(th) sqrt(sum((th - exact)^2))
fixed <- sapply(1:3, sgd_run, eta_fn = function(k) 0.01)
decay_eta <- function(k) 0.05 / (1 + k / 500)
decay <- sapply(1:3, sgd_run, eta_fn = decay_eta)
round(c(
full_batch = err(th_full),
sgd_fixed = apply(fixed, 2, err),
sgd_decay = apply(decay, 2, err)
), 4)full_batch sgd_fixed1 sgd_fixed2 sgd_fixed3 sgd_decay1 sgd_decay2 sgd_decay3
0.0460 0.0780 0.0698 0.0955 0.0371 0.0086 0.0388 All four use the same budget: 10,000 observation touches, which is 20 full passes.
The instructive part is the spread. Full batch is deterministic — run it again and you get the identical answer. Fixed-step SGD gives a different answer every seed, scattered around the solution rather than converging to it, because a constant \(\eta\) leaves a noise floor it cannot get below. That is exactly the warning in Section 16.7, made concrete.
Decaying the step size shrinks the scatter, as the classical conditions promise.
Note what this exercise does not show: SGD beating full batch. With \(n = 500\) and a well-conditioned problem, 20 full passes is plenty, and SGD’s noise is a pure cost. Its advantage appears when \(n\) is large enough that even a handful of full passes is unaffordable — then the comparison is 10,000 noisy steps against two or three exact ones, and the noisy steps win easily. The method is not better; it is better per unit of compute, and only once that unit is scarce.
6. Increase the condition number to 100 by using \(f = 100x^2 + y^2\) and observe the effect on convergence.
predicted_rate observed_ratio
0.9801980 0.9802041 The per-step contraction matches Equation 16.3’s \((\kappa-1)/(\kappa+1) = 0.980\) — the loss is the squared distance, hence the square root. Ten times worse conditioning means roughly ten times as many iterations, which is the practical content of Equation 16.3.