Gradient-Based Methods

Chapter 19 left us with a skeleton rather than a method. We know that from a current iterate \(\mathbf{x}_k\) we should move along a descent direction \(\mathbf{d}_k\) — one satisfying \(\nabla F(\mathbf{x}_k)^{\top}\mathbf{d}_k < 0\) — by a step length \(\alpha_k > 0\) found by a line search. What the framework never said is which descent direction to choose, and that single unanswered question is the whole content of this chapter.

Writing the direction this way is the organising idea of the chapter. Every method below is a choice of the matrix \(\mathbf{B}_k\) — a model of the curvature at \(\mathbf{x}_k\). Note that whenever \(\mathbf{B}_k\) is positive definite, the descent test is automatic, since \(\nabla F^{\top}\mathbf{d}_k = -\nabla F^{\top}\mathbf{B}_k^{-1}\nabla F < 0\) for \(\nabla F \ne \mathbf{0}\). Positive definiteness of the curvature model is therefore not a technicality; it is what guarantees the method keeps going downhill at all.

The simplest possible choice is to ignore curvature entirely: set \(\mathbf{B}_k = \mathbf{I}\), so that \(\mathbf{d}_k = -\nabla F(\mathbf{x}_k)\). This is the method of steepest descent (also called gradient descent or the Cauchy method, after its originator in 1847). It moves along the direction in which \(F\) decreases fastest locally — and that word locally is the source of all its trouble.

When \(F\) is quadratic, \(q(\mathbf{x}) = \tfrac12\mathbf{x}^{\top}\mathbf{A}\mathbf{x} - \mathbf{b}^{\top}\mathbf{x}\), the exact step length can be written down in closed form. Along the line \(\mathbf{x}_k - \alpha\mathbf{g}_k\) with \(\mathbf{g}_k = \nabla q(\mathbf{x}_k)\), the one-dimensional function \(\phi(\alpha)\) is a parabola in \(\alpha\), and setting \(\phi'(\alpha) = 0\) gives:

An exact line search has a curious and revealing consequence. Because \(\phi'(\alpha_k) = 0\) means \(\nabla F(\mathbf{x}_{k+1})^{\top}\mathbf{d}_k = 0\), consecutive steepest-descent directions are orthogonal. Each step therefore turns a perfect right angle from the last, and the trajectory zigzags. On a well-conditioned bowl (near-circular contours) the zigzag is harmless. On an ill-conditioned problem — long, thin, valley-like contours — it is catastrophic: the method bounces from wall to wall of the valley while creeping along its floor.

This behaviour can be quantified. For a quadratic with Hessian \(\mathbf{A}\), define the condition number \(\kappa = \lambda_{\max}/\lambda_{\min}\), the ratio of the largest to smallest eigenvalue — geometrically, how elongated the contour ellipses are. Steepest descent with exact line search obeys the classical error bound:

The rate is linear, and the factor tells the whole story. For \(\kappa = 1\) the factor is \(0\) — one step lands exactly on the minimum. For \(\kappa = 10\) the factor is \(9/11 \approx 0.82\), tolerable. For \(\kappa = 1000\) it is \(0.998\), and the method effectively stalls. Since real engineering problems routinely have \(\kappa\) in the thousands, plain steepest descent is rarely the final answer — but it is the honest baseline against which everything else is measured.

Before refining the direction, we must settle the step. Chapter 18 gave exact methods — golden section, Fibonacci, quadratic interpolation — for minimising \(\phi(\alpha) = F(\mathbf{x}_k + \alpha\mathbf{d}_k)\). But an exact line search inside every iteration of a multivariable method is wasteful: we are solving a subproblem to high precision only to throw the point away and rebuild a new direction. Modern practice uses an inexact line search — accept any \(\alpha\) that makes enough progress.

What is "enough"? Simply requiring \(F(\mathbf{x}_{k+1}) < F(\mathbf{x}_k)\) is not enough — a sequence of ever-tinier decreases can converge to a non-stationary point. The Armijo condition (sufficient decrease) fixes this by demanding a decrease proportional to what the slope promised:

Armijo alone still permits absurdly small steps, which also satisfy it. Two remedies exist. The first is backtracking: start optimistically at \(\alpha = 1\) and shrink by a factor \(\rho\) until Armijo passes. Because we test large steps first, the accepted step is never needlessly small, and no second condition is needed.

The second remedy adds an explicit curvature condition, forbidding steps that stop while the slope is still steeply negative. Together with Armijo it forms the Wolfe conditions, the standard requirement quoted in convergence theorems and — crucially for Section 20-5 — the condition that keeps quasi-Newton matrices positive definite.

Steepest descent uses only the linear term of the Taylor model. Newton's method uses the whole quadratic model and minimises it exactly. Recall the second-order expansion from Chapter 19: \(F(\mathbf{x}_k + \mathbf{p}) \approx F(\mathbf{x}_k) + \nabla F_k^{\top}\mathbf{p} + \tfrac12\mathbf{p}^{\top}\mathbf{H}_k\mathbf{p}\). Differentiate this model with respect to \(\mathbf{p}\), set the result to zero, and the Newton direction falls out immediately.

Two properties follow at once. First, on a quadratic function the model is not an approximation but the function itself, so Newton's method lands on the exact minimiser in a single step from any starting point. Second, on a general smooth \(F\), the model becomes ever more accurate as the iterates approach \(\mathbf{x}^\star\), and the method converges quadratically — the number of correct digits roughly doubles each iteration. Near the solution, three or four iterations is typical.

The price is threefold, and it is steep. Each iteration needs the \(n \times n\) matrix of second partial derivatives (\(O(n^2)\) entries, often unavailable analytically), then a linear solve (\(O(n^3)\) work). Worse, the method can fail. If \(\mathbf{H}_k\) is not positive definite — entirely possible far from the solution — the Newton direction need not even be a descent direction, and the method may march confidently toward a saddle point or a maximum. If \(\mathbf{H}_k\) is singular, the step is undefined.

Both defects are repaired by the same two ideas, giving the modified Newton method used in practice. First, damp the step: instead of accepting \(\alpha_k = 1\) blindly, run a line search on the Newton direction, so the method cannot overshoot far from the solution. (Near the solution the line search accepts \(\alpha_k = 1\) automatically, preserving quadratic convergence.) Second, modify the matrix: replace \(\mathbf{H}_k\) by \(\mathbf{H}_k + \mu\mathbf{I}\) with \(\mu > 0\) large enough to make it positive definite.

The shift \(\mu\mathbf{I}\) raises every eigenvalue by \(\mu\), so a large enough \(\mu\) guarantees positive definiteness and hence a descent direction. Notice what the two extremes of \(\mu\) give: as \(\mu \to 0\) we recover the pure Newton step, and as \(\mu \to \infty\) the direction tends toward \(-\nabla F_k/\mu\), a short steepest-descent step. The parameter therefore interpolates continuously between the two methods — cautious and gradient-like when far away, aggressive and Newton-like when close. This is exactly the idea behind the Levenberg–Marquardt method familiar from nonlinear least squares and curve fitting.

Newton's speed comes from curvature; its cost comes from computing curvature. Quasi-Newton methods ask the obvious question: can we get the curvature information for free, from the gradients we are already computing? The answer is yes, and it produces the most successful family of general-purpose optimizers ever devised.

The insight is that gradients differences already contain second-derivative information. Define the step and the gradient change over one iteration:

For a quadratic with Hessian \(\mathbf{A}\), these satisfy \(\mathbf{y}_k = \mathbf{A}\mathbf{s}_k\) exactly. So we demand that our approximation \(\mathbf{B}_{k+1}\) reproduce this same relationship — the secant condition, the multivariable generalisation of the secant method's slope estimate.

The secant condition is \(n\) equations for the \(n(n+1)/2\) unknowns of a symmetric matrix, so it does not determine \(\mathbf{V}_{k+1}\) — it leaves freedom, and different ways of using that freedom give different methods. The standard resolution is to take the previous matrix and add the smallest correction that satisfies the condition while preserving symmetry and positive definiteness. Two rank-two corrections have proved decisive.

Both updates preserve symmetry, both satisfy the secant condition, and both preserve positive definiteness — but only if the curvature condition \(\mathbf{s}_k^{\top}\mathbf{y}_k > 0\) holds. This is precisely what the second Wolfe condition of Section 20-3 guarantees, which is why the two topics belong in the same chapter: the line-search rule is not decoration, it is what keeps the update legal.

Started from \(\mathbf{V}_0 = \mathbf{I}\), a quasi-Newton method therefore takes a steepest-descent step first and grows into a Newton-like method as it accumulates curvature. On a quadratic in \(n\) variables the accumulated \(\mathbf{V}_n\) equals \(\mathbf{A}^{-1}\) exactly and the method terminates in \(n\) steps. On general functions it converges superlinearly — slower than Newton's quadratic rate, but with no Hessian and no \(O(n^3)\) solve. BFGS is the default choice in virtually all modern software: it is markedly more robust than DFP under inexact line searches, and its self-correcting behaviour recovers gracefully from a bad curvature estimate that leaves DFP stuck.

Quasi-Newton methods buy speed with \(O(n^2)\) storage — a full matrix. When \(n\) is very large (think \(10^6\) variables in a finite-element or machine-learning model), that matrix cannot be stored at all. Conjugate gradient (CG) methods achieve finite termination on quadratics while storing only a handful of vectors.

The idea attacks the zigzag directly. Steepest descent zigzags because its successive directions are merely orthogonal — each new step partially undoes the progress of the last. The cure is to make directions orthogonal in the geometry of the problem itself, measured through \(\mathbf{A}\) rather than through the identity. Directions \(\mathbf{d}_0, \mathbf{d}_1, \dots\) are \(\mathbf{A}\)-conjugate if

The payoff is the expanding-subspace property: minimising along a set of \(\mathbf{A}\)-conjugate directions in turn means that a minimisation along \(\mathbf{d}_k\) never spoils the minimisations already achieved along \(\mathbf{d}_0,\dots,\mathbf{d}_{k-1}\). Each direction settles one "coordinate" of the problem permanently. Since \(n\) conjugate directions span \(\mathbb{R}^n\), the method reaches the exact minimiser of a quadratic in at most \(n\) steps. Conjugacy is, in effect, the geometrically correct notion of a right angle for this problem — it is orthogonality after the elongated contours have been mentally stretched back into circles.

Remarkably, such directions can be generated on the fly, without ever knowing \(\mathbf{A}\), by adding a single correction term to the negative gradient:

Two classical formulas for \(\beta_k\) are in common use. They are identical on a quadratic with exact line search and differ on general functions:

On non-quadratic functions, Polak–Ribière is usually faster in practice, because when progress stalls the numerator \(\nabla F_k - \nabla F_{k-1}\) becomes small, \(\beta_k \to 0\), and the method automatically restarts itself with a fresh steepest-descent direction. (For theoretical guarantees one uses \(\beta_k = \max(\beta_k^{\mathrm{PR}}, 0)\).) In either case, deliberate restarting every \(n\) iterations is standard practice, since the conjugacy generated from a nonquadratic function degrades as the iterations accumulate.

All methods in this chapter share the same stopping logic. Terminate when the gradient is small, \(\lVert\nabla F(\mathbf{x}_k)\rVert \le \varepsilon\) — the direct test of the first-order condition — and in practice back it up with a check that the iterates and function values have stopped moving, \(\lVert\mathbf{x}_{k+1} - \mathbf{x}_k\rVert \le \varepsilon_x\) and \(|F(\mathbf{x}_{k+1}) - F(\mathbf{x}_k)| \le \varepsilon_F\), plus an iteration cap. Note carefully what none of them can promise: all converge to a stationary point, and only convexity (Chapter 19) makes that point a global minimum.

Two closing cautions. First, scaling matters more than the algorithm. A problem whose variables are measured in millimetres and megawatts has an enormous condition number for a purely clerical reason, and simply rescaling the variables to comparable magnitudes often does more good than switching methods. Second, all of these are local methods — they find the valley they were dropped into. For non-convex problems with many valleys, multiple starting points, or the global metaheuristics of Part 6, are the answer.

For each, compute gradients (and Hessians where needed), apply the relevant iteration, and state the result. Difficulty rises down the list.