Gradient Descent Theory, in depth: steepest descent under a norm, the eigenvalue view of step size, preconditioning, momentum, PL and implicit bias

What a gradient step minimises

Fix a point x with gradient g. A gradient step is the exact minimiser of a simple model of f around x: the linear Taylor term plus a penalty on how far you move.

x_next = argmin_d  f(x) + <g, d> + (1 / (2 * eta)) * ||d||^2      =>   d = -eta * g

Two things are hidden in that line. First, eta is the inverse of the curvature you assume: the model says f bends upward like a quadratic of curvature 1/eta in every direction. If the true curvature in some direction exceeds 2/eta, the model is badly wrong there and the step overshoots. Second, the gradient is the steepest direction only for the Euclidean norm. Change the norm and the minimiser changes:

• With ||d||_P^2 = d^T P d for a positive definite P, the step is d = -eta * P^-1 g. That is preconditioned gradient descent.
• With P equal to the Hessian, the step is Newton's method: the model matches the true curvature everywhere, and on a quadratic one step with eta = 1 lands on the minimum.
• With the l-infinity norm (every coordinate may move by the same amount), the steepest direction is -sign(g), the update behind sign-based optimisers.
• With the l1 norm, the steepest direction moves only the coordinate with the largest gradient magnitude, which is greedy coordinate descent.

So 'which optimiser' is largely the question 'which norm do I measure step length in', and a good norm is one under which the curvature looks uniform. Equivalently, gradient descent is the forward Euler discretisation of gradient flow dx/dt = -grad f(x), and every step-size limit below is the stability limit of an explicit integrator on a stiff system.

Quadratics: gradient descent as a linear system

Take the quadratic f(x) = 1/2 x^T A x - b^T x with A symmetric positive definite, minimiser x* = A^-1 b, smallest eigenvalue mu and largest L. The error e = x - x* obeys an exact linear recurrence:

e_{k+1} = (I - eta * A) e_k
In the eigenbasis of A, component i is multiplied by (1 - eta * lambda_i) every step.
rho(eta) = max_i |1 - eta * lambda_i| = max(|1 - eta * mu|, |1 - eta * L|)

The iteration converges if and only if rho is below 1, which requires 0 < eta < 2 / L. The factor is V-shaped in lambda, so only the two ends of the spectrum matter: small steps make the mu component crawl, large ones make the L component oscillate and grow. Equalising the ends gives the best fixed step:

eta* = 2 / (mu + L)          rho* = (L - mu) / (L + mu) = (kappa - 1) / (kappa + 1)
iterations to shrink the error by eps  ~  ln(1/eps) / -ln(rho)  ~  (kappa / 2) * ln(1/eps)

The condition number kappa = L / mu is the only property of A that matters. Near a non-degenerate minimum every smooth function is a quadratic with A equal to its Hessian, so this recurrence describes gradient descent locally on any model.

Measured: step size against the spectrum

The script below runs plain gradient descent on A = diag(1, 10, 100) (so kappa = 100) from the all-ones vector and counts iterations until the error norm falls by a factor of one million. It also predicts the count from rho.

import numpy as np, math

def gd(grad, x0, lr, tol=1e-6, max_iter=100000):
    x = np.array(x0, dtype=float)
    r0 = np.linalg.norm(x)
    for k in range(1, max_iter + 1):
        x = x - lr * grad(x)
        n = np.linalg.norm(x)
        if not np.isfinite(n) or n > 1e12 * r0:
            return x, None                      # diverged
        if n <= tol * r0:
            return x, k
    return x, max_iter

lam = np.array([1.0, 10.0, 100.0]); A = np.diag(lam)
for lr in [1/100, 2/101, 1.9/100, 1.99/100, 2.01/100]:
    _, k = gd(lambda x: A @ x, np.ones(3), lr)
    rho = max(abs(1 - lr * l) for l in lam)
    print(lr, round(rho, 4), k)

Measured counts sit just below the bound because only part of the starting error lies in the slowest direction. The steps 1/L and 1.99/L give identical counts: one is limited by the mu end, the other by the L end oscillating at 0.99. And the 1 percent change from 1.99/L to 2.01/L separates converging from overflowing, which is why a learning-rate sweep ends in a cliff, not a slope.

Rotation, preconditioning and the condition number

Plain gradient descent is rotation invariant. Rotating the problem with a random orthogonal Q, so that B = Q A Q^T has the same eigenvalues but no axis-aligned structure, took 676 iterations at the optimal step against 681 for the diagonal version; the small difference comes from how the fixed starting vector projects onto the eigenvectors. Only the spectrum matters.

Preconditioning changes the spectrum. Running x -= eta * P^-1 grad f(x) is plain gradient descent on the substituted variable y = P^(1/2) x, whose Hessian is P^(-1/2) A P^(-1/2). With A diagonal and P equal to its diagonal, that Hessian is the identity and the measured run converged in 1 iteration. On the rotated B, the same recipe (P equal to the diagonal of B) only reduced kappa from 100 to 8.39, and the optimal step then needed 60 iterations. Per-coordinate scaling fixes curvature that is aligned with the coordinates and does nothing for curvature that is spread across them.

Standardising input features is a diagonal preconditioner that works because unscaled features create axis-aligned curvature. Adam and RMSProp estimate a diagonal preconditioner from squared gradients, as derived in the Adam optimizer article, and share its blind spot for rotated curvature. Newton methods attack the full matrix, at quadratic storage cost.

The lower bound and momentum

Is kappa iterations unavoidable? No. For methods whose iterates stay in the span of the gradients seen so far, Nesterov proved a lower bound of order sqrt(kappa) * ln(1/eps) iterations on quadratics, and momentum attains it. Polyak's heavy ball adds a fraction of the previous step:

x_{k+1} = x_k - eta * grad f(x_k) + beta * (x_k - x_{k-1})
eta = 4 / (sqrt(L) + sqrt(mu))^2     beta = ((sqrt(kappa) - 1) / (sqrt(kappa) + 1))^2
asymptotic rate = sqrt(beta) = (sqrt(kappa) - 1) / (sqrt(kappa) + 1)   -> 9/11 = 0.818 for kappa = 100

On the same kappa = 100 problem heavy ball needed 92 iterations against 681 for the best plain step. The pure rate predicts about 69; the gap is real and has a cause. At this tuning the extremal eigen-directions of the two-step iteration matrix have repeated eigenvalues, so their error decays like k times rho to the k rather than rho to the k, and the extra factor of k costs a couple of dozen iterations. This guarantee is for quadratics. On general strongly convex functions these parameters can fail to converge (Lessard, Recht and Packard showed a counterexample), which is why Nesterov's accelerated gradient, with its look-ahead gradient, is the method with the general guarantee.

Beyond strong convexity: PL and implicit bias

Strong convexity (mu above 0) is a strong assumption. Many objectives that converge linearly in practice violate it: an over-parameterised least squares problem has a whole affine space of minimisers, so its Hessian has zero eigenvalues. The condition that actually delivers a linear rate is the Polyak-Lojasiewicz (PL) inequality:

1/2 * ||grad f(x)||^2  >=  mu * (f(x) - f*)       for all x
with eta = 1/L:   f(x_k) - f*  <=  (1 - mu / L)^k * (f(x_0) - f*)

PL needs neither a unique minimiser nor convexity. For least squares f(w) = 1/2 ||X w - y||^2 with more columns than rows, it holds with mu equal to the smallest non-zero squared singular value of X.

The measured case: X is 20 by 50 with Gaussian entries, so infinitely many w fit y exactly. The squared singular values run from 5.54 to 133.5, an effective kappa of 24.1. Gradient descent from w = 0 with eta = 1/L drove the residual below 1e-8 in 434 iterations, and the solution it found was within 4.2e-9 of the minimum-norm solution pinv(X) @ y.

That second fact is implicit bias. Every gradient X^T (X w - y) lies in the row space of X, so iterates started at zero never gain a null-space component, and the only exact fit inside the row space is the minimum-norm one. Start from w0 and you get the fit closest to w0: in under-determined problems, initialisation is part of the model.

Nonconvex objectives: stationarity and saddles

Without convexity, the descent lemma still gives a guarantee with eta = 1/L: the smallest gradient norm seen in k steps satisfies min ||grad f||^2 <= 2 L (f(x_0) - f_low) / k. That is a guarantee of stationarity, not of a minimum, and stationary points include saddles.

The linear-system view says exactly how saddles behave. Near a saddle the Hessian has a negative eigenvalue lambda, so the component along it is multiplied by 1 + eta * |lambda| per step: it grows, but only from wherever it started. On f(x, y) = x^2/2 - y^2/2 + y^4/4 (saddle at the origin, minima at y = plus or minus 1) with eta = 0.1, starting at (1, y0):

The difference of 97 iterations matches ln(1e4) / ln(1.1) = 96.6: escape time is logarithmic in the initial offset. Lee, Simchowitz, Jordan and Recht proved that from a random initialisation gradient descent avoids strict saddles with probability one, but this experiment shows the catch: 'avoids' can still mean a long plateau. The geometric picture of saddles, plateaus and sharpness in deep networks is covered in the optimization landscape article.

Operational guidance

• Estimate L before you tune. Run a few iterations of power iteration on Hessian-vector products (or on X^T X for linear models) to estimate the top eigenvalue, then start the step-size sweep near 1/L and expect a cliff near 2/L.
• Read divergence by its shape. A loss that oscillates with growing amplitude means eta times the top curvature is above 2. A loss that falls steadily but slowly means the small end of the spectrum is the bottleneck, and a larger step will not help much; preconditioning or momentum will.
• Precondition by the structure you have. Standardise features, normalise activations and use per-parameter scaling when curvature is roughly axis-aligned; reach for full-matrix or block methods only when it is not.
• Fix initialisation when runs must agree; implicit bias makes the answer depend on it.

Failure modes

• Step size copied from a different problem. L changes with data scaling and architecture, so a borrowed learning rate can sit beyond 2/L.
• Assuming a bigger step is always faster. Past 2/(mu+L) the L end dominates and the rate gets worse, then collapses at 2/L.
• Trusting diagonal scaling on rotated curvature. Correlated features leave most of kappa in place, as the 100 to 8.39 measurement shows.
• Reading a small gradient as a minimum. Saddles and plateaus also have small gradients; check the loss and, where cheap, the sign of a curvature estimate.

Trade-offs

For small convex problems, prefer the second-order solvers in the convex optimization article; for large-scale training, momentum plus a diagonal preconditioner is the usual compromise.

What to do next

1. Run the script above, then change the middle eigenvalue and confirm the iteration counts do not move: only mu and L matter.
2. Rotate the problem with a random orthogonal matrix and compare diagonal scaling against the axis-aligned case on your own spectrum.
3. For one model you train, estimate the top Hessian eigenvalue with power iteration and compare 2/L with the learning rate where your sweep diverges.
4. Fit an under-determined least squares problem from two different starts and measure the distance between the two solutions.
5. Implement Nesterov's method next to heavy ball and compare them on a logistic regression, where the quadratic tuning formula no longer applies exactly.