Gradient descent

The least-squares note ended with two problems of the closed form. It needs the $d \times d$ matrix $X^\top X$ and a solve that costs $O(d^3)$, out of reach for a million features; and it exists only for squared error, while the Huber loss, the absolute loss and every classification loss in this course have no formula for their minimizer. Both problems have one answer. Instead of solving for the point where the gradient vanishes, walk towards it: start anywhere, compute the gradient, take a small step against it, repeat. This note analyses that walk on least squares, where everything can be computed exactly, and then states what survives for general losses.

Notation. $Q(w)$ is the loss to minimize, here the mean squared error $\frac1n\lVert Xw - y\rVert^2$ of the least-squares note, with design matrix $X \in \mathbb R^{n \times d}$, and $\hat w$ its minimizer. $\nabla Q(w)$ is the gradient, $w^{(k)}$ the iterate after $k$ steps and $\eta > 0$ the step size or learning rate. For least squares the Hessian is the constant matrix $A = \frac2n X^\top X$; its smallest and largest eigenvalues are $\mu$ and $L$, and $\kappa = L/\mu$ is its condition number.

Which way is downhill?#

Near a point $w$, a differentiable function is close to its linear approximation:

$$Q(w + h) = Q(w) + \langle \nabla Q(w), h\rangle + o(\lVert h\rVert) .$$

Among all steps $h$ of a fixed small length $\varepsilon$, the change is dominated by $\langle \nabla Q(w), h\rangle$, and by the Cauchy–Schwarz inequality this is at least $-\varepsilon\lVert\nabla Q(w)\rVert$, with equality only for $h = -\varepsilon\,\nabla Q(w)/\lVert\nabla Q(w)\rVert$. The negative gradient is the direction of steepest descent, and a small enough step along it lowers $Q$ whenever the gradient is not zero. Gradient descent repeats such steps:

$$w^{(k+1)} = w^{(k)} - \eta\, \nabla Q\big(w^{(k)}\big) .$$

For least squares, the gradient from the normal-equations derivation, divided by $n$, is

$$\nabla Q(w) = \frac2n X^\top (Xw - y) ,$$

two matrix–vector products, $O(nd)$ operations, and only $O(\text{number of nonzeros})$ when $X$ is sparse. Nothing of size $d \times d$ is ever formed.

Worked example. The four points of the least-squares note, $x = (0, 1, 2, 3)$, $y = (1, 3, 2, 5)$, with columns $1$ and $x$. From $w^{(0)} = 0$ the residual is $-y$, so $\nabla Q = -\frac12 X^\top y = -\frac12 (11, 22) = (-5.5, -11)$, and $Q = (1 + 9 + 4 + 25)/4 = 9.75$. With $\eta = 0.2$, $w^{(1)} = (1.1, 2.2)$ and $Q = 4.91$; the next gradient is $(3.3, 7.7)$, $w^{(2)} = (0.44, 0.66)$ and $Q = 2.659$. The loss falls towards its minimum $2.7/4 = 0.675$ at $\hat w = (1.1, 1.1)$, but the slope coefficient went from 0 past its target to $2.2$ and back to $0.66$. The overshoot is the subject of the next section.

How large a step? The quadratic case#

For least squares the loss is exactly quadratic around its minimum. Its gradient vanishes at $\hat w$ and its Hessian is $A$ everywhere, so

$$Q(w) = Q(\hat w) + \tfrac12 (w - \hat w)^\top A\, (w - \hat w), \qquad \nabla Q(w) = A\,(w - \hat w) .$$

Write the error as $e^{(k)} = w^{(k)} - \hat w$. One step gives

$$e^{(k+1)} = e^{(k)} - \eta A\, e^{(k)} = (I - \eta A)\, e^{(k)} .$$

$A$ is symmetric with eigenvalues $\mu = \lambda_1 \le \dots \le \lambda_d = L$ and orthonormal eigenvectors $u_1, \dots, u_d$. In that basis $I - \eta A$ is diagonal, and the component of the error along $u_i$ is multiplied by $1 - \eta\lambda_i$ at every step, independently of the others. Three consequences:

• The largest safe step. All components shrink, from every starting point, if and only if $\lvert 1 - \eta\lambda_i\rvert < 1$ for all $i$, that is $0 < \eta < 2/L$ (assuming $\mu > 0$). At $\eta > 2/L$ the steepest component flips sign and grows at every step: descent diverges from almost every start, however close to the minimum; only a start whose steepest component is exactly zero escapes.
• The best constant step. The speed is set by the worst factor $\rho(\eta) = \max\big(\lvert 1 - \eta\mu\rvert, \lvert 1 - \eta L\rvert\big)$. The first term decreases in $\eta$ and the second eventually increases, so the maximum is smallest where $1 - \eta\mu = \eta L - 1$: $\eta^* = 2/(L + \mu)$, with $\rho^* = \dfrac{L - \mu}{L + \mu} = \dfrac{\kappa - 1}{\kappa + 1}$.
• The number of steps. Shrinking the error by a factor $\varepsilon$ takes $\ln(1/\varepsilon)/\ln(1/\rho^*)$ steps, and $\ln\frac{\kappa + 1}{\kappa - 1} \approx 2/\kappa$ for large $\kappa$, so about $\tfrac{\kappa}{2}\ln(1/\varepsilon)$: linear in the condition number.

Figure

How much one step shrinks the error

Worked example. For the four points, $A = \frac24 X^\top X = \begin{pmatrix} 2 & 3 \\ 3 & 7\end{pmatrix}$. Its trace is 9 and its determinant 5, so the eigenvalues are $(9 \pm \sqrt{61})/2$: $\mu = 0.595$ and $L = 8.405$, $\kappa = 14.1$. Steps above $2/L = 0.238$ diverge; the best step is $2/(L + \mu) = 2/9 = 0.222$, with factor $\rho^* = 0.868$, and a thousandfold reduction of the error takes $\ln 1000/\ln(1/0.868) = 49$ steps. The step $\eta = 1/L = 0.119$, the usual safe default, has the factor $1 - 1/\kappa = 0.929$ as a bound and needs 80 steps from $w = 0$. The data-features note’s example of two features on scales 1 and 100 has $\kappa \approx 10^4$, and the same thousandfold reduction needs about $\tfrac{10^4}{2}\ln 1000 \approx 35{,}000$ steps.

Why is the error so slow along $u_1$? Because $\mu$ is the curvature in the flattest direction: the gradient there is small, so the steps are small, while the step size cannot grow because of the steepest direction. The iterate zigzags across a narrow valley and creeps along its floor.

Scaling the features#

The eigenvalues of $A = \frac2n X^\top X$ depend on the units of the features. In the four-point example the off-diagonal entry of $A$ is $\frac2n\sum_i x_i = 3$, the coupling between intercept and slope created by the feature’s nonzero mean, and the diagonal entries 2 and 7 differ because $x$ has a larger second moment than the constant column. Centring $x$ zeroes the coupling; scaling it to unit variance equalizes the diagonal. For the standardized feature $(x - 1.5)/1.118$, $A = 2I$, $\kappa = 1$, every direction has the same curvature, and the gradient points straight at the minimum: the step $\eta = 2/(L + \mu) = 0.5$ reaches $\hat w$ in a single step.

Figure

Gradient descent on a stretched and on a round bowl

The effect on real data is large. On the wine data bundled with scikit-learn, predicting alcohol content from the other 12 measurements (whose standard deviations range from about 0.1 to 300), the condition number of the least-squares Hessian is $1.6 \cdot 10^8$ with raw features and 45 after standardization. Gradient descent with $\eta = 1/L$ reduces the excess loss $Q(w) - Q(\hat w)$, relative to its initial value, to $10^{-1}$, $9 \cdot 10^{-2}$, $2 \cdot 10^{-2}$ and $10^{-2}$ after 10, 100, 1,000 and 10,000 steps on the raw features; on the standardized features to $8 \cdot 10^{-3}$ after 10 steps, $2 \cdot 10^{-8}$ after 100 and machine precision after 1,000 (the Code section). Scaling cannot remove the remaining $\kappa = 45$, which comes from correlations between the features. That would need a rotation as well, which is what second-order methods effectively apply.

Standardization changes nothing about the model: a linear model on standardized features is a linear model on the original ones with rescaled weights and a shifted intercept. It changes only the geometry the optimizer walks through. The scaling statistics are computed on the training data, as the features note requires.

Beyond quadratics: what is guaranteed#

Most losses are not quadratic, but the analysis carries over through two properties. A differentiable $Q$ is $L$-smooth if its gradient is $L$-Lipschitz, $\lVert\nabla Q(u) - \nabla Q(v)\rVert \le L\lVert u - v\rVert$; for twice differentiable $Q$ this means that no eigenvalue of the Hessian exceeds $L$ in absolute value. It is $\mu$-strongly convex if $Q(v) \ge Q(u) + \langle\nabla Q(u), v - u\rangle + \frac\mu2\lVert v - u\rVert^2$; for twice differentiable $Q$, every Hessian eigenvalue is at least $\mu$. Least squares with independent columns is both, with the $\mu$ and $L$ of $A$.

The descent lemma. Smoothness bounds the function above by a parabola: $Q(v) \le Q(u) + \langle\nabla Q(u), v - u\rangle + \frac L2\lVert v - u\rVert^2$. Put $v = u - \eta\nabla Q(u)$:

$$Q\big(u - \eta\nabla Q(u)\big) \le Q(u) - \eta\Big(1 - \frac{L\eta}{2}\Big)\lVert\nabla Q(u)\rVert^2 .$$

Every step with $0 < \eta < 2/L$ decreases the loss, the same threshold as in the quadratic case, and the guaranteed decrease is largest at $\eta = 1/L$, where it equals $\lVert\nabla Q(u)\rVert^2/(2L)$. Three guarantees follow. They are the theory behind the practical routine, which needs none of the rates: compute the gradient, check that the loss falls at every step, scale the inputs, and keep the largest step that does not diverge.

• Nonconvex, $L$-smooth, $\eta = 1/L$. Summing the decreases over $K$ steps, $\frac{1}{2L}\sum_{k < K}\lVert\nabla Q(w^{(k)})\rVert^2 \le Q(w^{(0)}) - Q(w^{(K)}) \le Q(w^{(0)}) - \inf Q$, so the smallest gradient norm seen satisfies $\min_{k < K}\lVert\nabla Q(w^{(k)})\rVert^2 \le 2L\big(Q(w^{(0)}) - \inf Q\big)/K$. The gradient becomes small, nothing more: the limit can be a local minimum or a saddle point. This is the only general guarantee for neural networks, and it assumes smoothness, which networks with ReLU activations do not have everywhere.
• Convex, $L$-smooth, $\eta = 1/L$. $Q(w^{(k)}) - Q(\hat w) \le 2L\lVert w^{(0)} - \hat w\rVert^2/k$ (Bubeck, 2015, Theorem 3.3): the excess loss falls like $1/k$.
• $\mu$-strongly convex and $L$-smooth, $\eta = 2/(L + \mu)$. $Q(w^{(k)}) - Q(\hat w) \le \frac L2 \exp\big(-\frac{4k}{\kappa + 1}\big)\lVert w^{(0)} - \hat w\rVert^2$ (Bubeck, 2015, Theorem 3.12): linear convergence, with the number of steps proportional to $\kappa$, as on the quadratic.

When $X$ has dependent columns, least squares is convex but not strongly convex ($\mu = 0$): the minimizer is not unique, and the general convex theorem above gives only the $1/k$ rate. Least squares is still a fixed quadratic, though, and the eigenvector analysis of the quadratic case still applies: with $0 < \eta < 2/L$, the error components along eigenvectors with positive eigenvalues shrink geometrically, at a speed set by the smallest positive eigenvalue, and the component in the null space of $A$ does not change at all. Gradient descent started at $w = 0$ then converges geometrically to the minimum-norm solution, because every step adds a vector of the form $X^\top v$ and stays in the row space of $X$.

Choosing the step in practice#

$L$ is rarely known. Three practical answers:

• A grid on a log scale. Try $\eta = 1, 0.1, 0.01, \dots$ for a few hundred steps and keep the largest that does not diverge. A loss that grows from the first steps means $\eta$ is above $2/L$.
• Backtracking line search. At each step, start from a large $\eta$ and halve it until $Q(w - \eta g) \le Q(w) - \frac\eta2\lVert g\rVert^2$ with $g = \nabla Q(w)$ (Boyd and Vandenberghe, 2004, section 9.2, here with the constant 1/2). By the descent lemma every $\eta \le 1/L$ passes the test, so the search ends with $\eta \ge 1/(2L)$ and needs no knowledge of $L$. It also accepts longer steps when the gradient points along a flat direction. On the standardized wine problem, where $1/L = 0.108$, 100 backtracking steps from $\eta = 1$ accepted $0.0625$ twelve times, $0.125$ 41 times, $0.25$ nine times, $0.5$ twelve times and $1$ 26 times, and ended with a gradient norm of $2 \cdot 10^{-6}$, against $10^{-3}$ for 100 steps of the fixed $\eta = 1/L$.
• Stopping. Stop when the gradient norm or the relative decrease of $Q$ falls below a tolerance. In machine learning there is a better criterion than the training loss: stop when the validation loss stops improving. The number of steps then acts as a regularizer, as the training note shows.

Newton’s method uses the curvature too: $w \leftarrow w - (\nabla^2 Q)^{-1}\nabla Q$. On a quadratic it reaches $\hat w$ in one step from anywhere, because it solves the normal equations; in general it converges in a handful of steps near the minimum, at the price of a $d \times d$ solve per step. The logistic-regression note uses it where $d$ is small. Gradient descent is the method that scales.

Losses without a closed form#

Gradient descent needs only the gradient, so it works unchanged for the losses that the closed form could not handle. For the Huber loss with threshold $\delta$, the derivative of $\rho_\delta(r)$ is the residual clipped to $[-\delta, \delta]$, so

$$\nabla Q(w) = -\frac1n X^\top \operatorname{clip}(y - Xw,\, -\delta,\, \delta) .$$

The second derivative of $\rho_\delta$ is at most 1, so $Q$ is $L$-smooth with $L$ the largest eigenvalue of $\frac1n X^\top X$, and $\eta = 1/L$ is safe. On the ten points with one outlier from the least-squares note, 5,000 steps reach $(0.986, 0.533)$, the Huber line of that note’s figure. The absolute loss is not differentiable at zero residual; its subgradient $\operatorname{sign}(r)$ can replace the gradient, but subgradient steps must shrink over time to converge, a theme of the next note. The logistic regression and linear classifiers notes apply the same update to classification losses, backpropagation computes the gradient for neural networks, and gradient boosting takes the steps in a space of functions.

What changed#

	Normal equations	Gradient descent
What it needs	$X^\top X$ and a $d \times d$ solve	the gradient, $O(nd)$ per step
Cost	$O(nd^2 + d^3)$ once	$O(nd)$ times about $\kappa\ln(1/\varepsilon)$ steps
Losses	squared error only	any differentiable loss
Accuracy	exact up to rounding	stops at a tolerance
Hyperparameters	none	step size $\eta$, number of steps
Sensitivity to scaling	numerical accuracy ($\kappa^2$ if $X^\top X$ is formed)	speed ($\kappa$ steps); diverges for $\eta > 2/L$

Limitations#

• Speed depends on the condition number. Badly scaled or correlated features make descent crawl; scaling fixes the first, not the second.
• The step size must be chosen. Too small wastes steps, above $2/L$ diverges; a line search costs extra loss evaluations.
• Every step is a full pass over the data. With $n$ in the millions one step is expensive, and most of the information in it is redundant.
• Nonconvex losses. Only a small gradient is guaranteed; the result depends on the starting point.
• Non-differentiable losses need subgradients or smoothing, and constant steps no longer converge.

What you should remember#

• The negative gradient is the direction of steepest descent; gradient descent iterates $w \leftarrow w - \eta\nabla Q(w)$. For least squares $\nabla Q = \frac2n X^\top(Xw - y)$.
• On a quadratic with Hessian eigenvalues in $[\mu, L]$: converges iff $\eta < 2/L$; the best constant step $2/(L + \mu)$ shrinks the error by $(\kappa - 1)/(\kappa + 1)$ per step; about $\frac\kappa2\ln(1/\varepsilon)$ steps.
• Standardize features before gradient descent: centring decouples the intercept, scaling equalizes the curvatures. On the wine data it lowers $\kappa$ from $1.6 \cdot 10^8$ to 45.
• For $L$-smooth losses with $\eta = 1/L$: each step decreases $Q$ by at least $\lVert\nabla Q\rVert^2/(2L)$; the excess loss is $O(1/k)$ for convex, $O(e^{-k/\kappa})$ for strongly convex losses; nonconvex losses only get a small gradient.
• Backtracking line search finds a good step without knowing $L$.
• Gradient descent handles losses with no closed form, such as Huber, whose gradient is the clipped residual.

Code#

The worked example, step sizes on the raw and standardized problem, the wine data, backtracking and the Huber fit:

import numpy as np
from sklearn.datasets import load_wine

def gd(grad, w0, eta, steps):
    w = np.array(w0, float)
    for _ in range(steps):
        w = w - eta * grad(w)
    return w

# the four-point example: Q(w) = mean squared error, gradient (2/n) X^T (Xw - y)
x = np.array([0., 1, 2, 3]); y = np.array([1., 3, 2, 5])
for name, col in (('raw', x), ('standardized', (x - x.mean()) / x.std())):
    X = np.column_stack([np.ones(4), col])
    A = 2 / 4 * X.T @ X                                  # the Hessian of Q
    grad = lambda w: A @ w - 2 / 4 * X.T @ y
    w_hat = np.linalg.solve(X.T @ X, X.T @ y)
    mu, L = np.linalg.eigvalsh(A)
    print(name, A.ravel(), round(mu, 3), round(L, 3), round(L / mu, 2), round(2 / L, 3), round(2 / (L + mu), 3), round((L - mu) / (L + mu), 3))
    for eta in (0.05, 1 / L, 2 / (L + mu), 0.24):
        w, k = np.zeros(2), 0                            # steps until the distance to w-hat falls below 1/1000 of the start
        while np.linalg.norm(w - w_hat) > 1e-3 * np.linalg.norm(w_hat) and k < 1000 and np.linalg.norm(w) < 1e3:
            w, k = w - eta * grad(w), k + 1
        print('  ', round(eta, 3), k if k < 1000 and np.linalg.norm(w) < 1e3 else 'diverges')
# raw [2. 3. 3. 7.] 0.595 8.405 14.13 0.238 0.222 0.868
#    0.05 194 | 0.119 80 | 0.222 49 | 0.24 diverges
# standardized [2. 0. 0. 2.] 2.0 2.0 1.0 1.0 0.5 0.0
#    0.05 66 | 0.5 1 | 0.5 1 | 0.24 11

X = np.column_stack([np.ones(4), x])
w = np.zeros(2)
for k in range(3):                                       # the worked example: eta = 0.2 from w = 0
    print(k, w.round(2), (2 / 4 * X.T @ (X @ w - y)).round(2), round(np.mean((X @ w - y) ** 2), 3))
    w = w - 0.2 * 2 / 4 * X.T @ (X @ w - y)
# 0 [0. 0.] [ -5.5 -11. ] 9.75 | 1 [1.1 2.2] [3.3 7.7] 4.91 | 2 [0.44 0.66] [-2.64 -5.06] 2.659

# wine: predict alcohol from the other 12 measurements; raw against standardized features, eta = 1/L
D = load_wine().data; t, F = D[:, 0], D[:, 1:]
for name, M in (('raw', F), ('standardized', (F - F.mean(0)) / F.std(0))):
    Xw = np.column_stack([np.ones(len(t)), M])
    A = 2 / len(t) * Xw.T @ Xw
    ev = np.linalg.eigvalsh(A)
    Q = lambda w: np.mean((Xw @ w - t) ** 2)
    Q_hat, Q_0 = Q(np.linalg.lstsq(Xw, t, rcond=None)[0]), Q(np.zeros(13))
    w, out = np.zeros(13), []
    for k in range(1, 10_001):
        w -= 1 / ev[-1] * 2 / len(t) * Xw.T @ (Xw @ w - t)
        if k in (10, 100, 1000, 10_000): out.append(f'{(Q(w) - Q_hat) / (Q_0 - Q_hat):.0e}')
    print(name, f'{ev[-1] / ev[0]:.2g}', out)
# raw 1.6e+08 ['1e-01', '9e-02', '2e-02', '1e-02']: relative excess loss after 10, 100, 1000, 10000 steps
# standardized 45 ['8e-03', '2e-08', '7e-19', '7e-19']

# backtracking line search (halve eta until Q drops by at least eta/2 * |grad|^2) finds a step near 1/L by itself
Xs = np.column_stack([np.ones(len(t)), (F - F.mean(0)) / F.std(0)])
Q = lambda w: np.mean((Xs @ w - t) ** 2); grad = lambda w: 2 / len(t) * Xs.T @ (Xs @ w - t)
w, etas = np.zeros(13), []
for k in range(100):
    g, eta = grad(w), 1.0
    while Q(w - eta * g) > Q(w) - eta / 2 * g @ g:
        eta /= 2
    w, etas = w - eta * g, etas + [eta]
L = np.linalg.eigvalsh(2 / len(t) * Xs.T @ Xs)[-1]
w_fixed = gd(grad, np.zeros(13), 1 / L, 100)
print(round(1 / L, 3), {e: etas.count(e) for e in sorted(set(etas))}, f'{np.linalg.norm(grad(w)):.0e} {np.linalg.norm(grad(w_fixed)):.0e}')
# 0.108 {0.0625: 12, 0.125: 41, 0.25: 9, 0.5: 12, 1.0: 26} 2e-06 1e-03

# the Huber loss of the least-squares note has no closed form; its gradient is the clipped residual
xo = np.arange(10.); yo = np.array([1.3, 1.1, 2.2, 3.0, 2.7, 3.6, 3.5, 4.9, 13.0, 5.3]); Xo = np.column_stack([np.ones(10), xo])
huber_grad = lambda w: -Xo.T @ np.clip(yo - Xo @ w, -1, 1) / 10
L = np.linalg.eigvalsh(Xo.T @ Xo / 10)[-1]
print(gd(huber_grad, [0, 0], 1 / L, 5000).round(3))   # [0.986 0.533]

Next#

Each step of gradient descent touches all $n$ examples, and the gradient is an average of $n$ per-example gradients. An average can be estimated from a small random sample, and a cheap noisy step is often worth more than an exact expensive one. The next note makes that trade: stochastic gradients and minibatches, why constant steps then stop short of the minimum, and the momentum and per-coordinate step sizes (AdaGrad, RMSProp, Adam) that attack the condition number from this note.

The same update trains models across the courses: the generative models course fits distributions by gradient steps on a KL estimate, and a GAN alternates descent steps for two players.

References#

• E. Sokolov et al. Машинное обучение (HSE), lecture notes, lecture 3: training linear regression, gradient descent and its convergence.
• S. Bubeck. Convex optimization: algorithms and complexity. Foundations and Trends in Machine Learning, 8(3–4), 2015. Theorem 3.3 (smooth convex functions), Theorem 3.12 (smooth and strongly convex functions).
• S. Boyd, L. Vandenberghe. Convex Optimization. Cambridge University Press, 2004. Section 9.2 (descent methods and backtracking line search), Section 9.3 (gradient descent and its convergence), Section 9.5 (Newton’s method).
• I. Goodfellow, Y. Bengio, A. Courville. Deep Learning. MIT Press, 2016. Section 4.2 (poor conditioning), Section 4.3 (gradient-based optimization and the Hessian).
• M. P. Deisenroth, A. A. Faisal, C. S. Ong. Mathematics for Machine Learning. Cambridge University Press, 2020. Section 7.1 (optimization using gradient descent).