Conjugate Gradient Methods

Why This Matters

Many problems in ML reduce to solving large linear systems $Ax = b$ where $A$ is symmetric positive definite (SPD). Examples: solving the normal equations in least squares, computing natural gradient steps, solving preconditioned systems in second-order optimization, and inverting kernel matrices.

Direct methods (Cholesky factorization) cost $O(n^3)$ and require forming $A$ explicitly. Conjugate gradient solves $Ax = b$ using only matrix-vector products $Av$, never forming $A$. For sparse or structured matrices where $Av$ is cheap, CG is dramatically faster. It converges in at most $n$ iterations (exact arithmetic) but typically reaches a good approximate solution in far fewer.

Mental Model

Steepest descent moves in the direction of the negative gradient at each step. The problem: it revisits the same directions repeatedly, zigzagging toward the solution. CG eliminates this waste by choosing directions that are $A$-conjugate: each new direction is "orthogonal" with respect to the $A$-inner product to all previous directions. This means each direction extracts maximum progress without undoing work from previous steps.

Formal Setup

Definition

A-Conjugate Directions

Vectors $d_0, d_1, \ldots, d_{k}$ are A-conjugate (or conjugate with respect to $A$) if $d_i^\top A d_j = 0$ for all $i \neq j$. When $A = I$, this reduces to ordinary orthogonality.

Definition

Krylov Subspace

The Krylov subspace of order $k$ is:

$\mathcal{K}_k(A, r_0) = \text{span}\{r_0, Ar_0, A^2 r_0, \ldots, A^{k-1} r_0\}$

where $r_0 = b - Ax_0$ is the initial residual. CG finds the solution in $x_0 + \mathcal{K}_k(A, r_0)$ that minimizes the $A$-norm of the error.

The CG Algorithm

Given SPD matrix $A$, right-hand side $b$, initial guess $x_0$:

1. Set $r_0 = b - Ax_0$ (initial residual) and $d_0 = r_0$ (initial search direction).
2. For $k = 0, 1, 2, \ldots$:
   • Step length: $\alpha_k = r_k^\top r_k / (d_k^\top A d_k)$
   • Update solution: $x_{k+1} = x_k + \alpha_k d_k$
   • Update residual: $r_{k+1} = r_k - \alpha_k A d_k$
   • Conjugate direction coefficient: $\beta_{k+1} = r_{k+1}^\top r_{k+1} / (r_k^\top r_k)$
   • New search direction: $d_{k+1} = r_{k+1} + \beta_{k+1} d_k$

Each iteration requires one matrix-vector product ($Ad_k$), two inner products, and three vector updates. Storage: only $x_k$, $r_k$, $d_k$, and the product $Ad_k$.

Main Theorems

Theorem

CG Finite Termination

Statement

The conjugate gradient algorithm terminates in at most $n$ iterations with the exact solution $x^* = A^{-1}b$. Equivalently, $r_n = 0$.

Intuition

The directions $d_0, \ldots, d_{n-1}$ are A-conjugate and span $\mathbb{R}^n$. Each iteration solves the problem exactly in one new dimension. After $n$ dimensions, the entire space is covered.

Proof Sketch

The A-conjugacy of the search directions $d_0, \ldots, d_{k-1}$ can be proved by induction using the three-term recurrence. Since $A$ is $n \times n$, at most $n$ A-conjugate directions can exist. The minimization property ensures $x_k$ is the minimizer of $\frac{1}{2}x^\top A x - b^\top x$ over $x_0 + \mathcal{K}_k(A, r_0)$. At $k = n$, $\mathcal{K}_n(A, r_0) = \mathbb{R}^n$ (for generic $r_0$), so $x_n = x^*$.

Why It Matters

Finite termination distinguishes CG from generic iterative methods (which converge asymptotically but never reach the exact solution). In practice, finite termination is destroyed by floating-point errors, but it explains why CG converges so much faster than steepest descent.

Failure Mode

In floating-point arithmetic, rounding errors cause loss of A-conjugacy among the search directions. After about $n$ iterations, the directions are no longer conjugate and the residual may stagnate. Preconditioning and reorthogonalization are standard remedies.

Theorem

CG Convergence Rate

Statement

After $k$ iterations:

$\|x_k - x^*\|_A \leq 2 \left(\frac{\sqrt{\kappa} - 1}{\sqrt{\kappa} + 1}\right)^k \|x_0 - x^*\|_A$

where $\|v\|_A = \sqrt{v^\top A v}$ is the $A$-norm and $\kappa = \lambda_{\max}(A) / \lambda_{\min}(A)$.

Intuition

The convergence rate depends on $\sqrt{\kappa}$, not $\kappa$ itself. This is a quadratic improvement over steepest descent, which depends on $(\kappa - 1)/(\kappa + 1)$. For a matrix with condition number $\kappa = 10{,}000$, steepest descent reduces the error by a factor of $0.9999$ per step; CG reduces it by $0.98$ per step.

Proof Sketch

CG minimizes $\|x_k - x^*\|_A$ over all $x_0 + \mathcal{K}_k(A, r_0)$. This is equivalent to minimizing $\|p_k(A)(x_0 - x^*)\|_A$ over all polynomials $p_k$ of degree $k$ with $p_k(0) = 1$. The Chebyshev polynomials solve this min-max problem on $[\lambda_{\min}, \lambda_{\max}]$, giving the stated bound.

Why It Matters

This bound explains why preconditioning is so important: replacing $A$ by $M^{-1}A$ (where $M \approx A$ is cheap to invert) reduces $\kappa$ and accelerates convergence. It also explains why CG is fast for matrices with clustered eigenvalues, even if the condition number is large.

Failure Mode

The bound is pessimistic when eigenvalues are clustered. If $A$ has only $r$ distinct eigenvalues, CG converges in exactly $r$ steps (exact arithmetic). The Chebyshev bound captures worst-case behavior over $[\lambda_{\min}, \lambda_{\max}]$ and does not exploit clustering.

Nonlinear Conjugate Gradient

CG can be extended to minimize a general smooth function $f(x)$, not just the quadratic $\frac{1}{2}x^\top Ax - b^\top x$. Replace the residual $r_k$ with $-\nabla f(x_k)$ and compute $\alpha_k$ by line search.

The key difference is the choice of $\beta_k$:

• Fletcher-Reeves: $\beta_k = \|\nabla f_k\|^2 / \|\nabla f_{k-1}\|^2$
• Polak-Ribiere: $\beta_k = \nabla f_k^\top (\nabla f_k - \nabla f_{k-1}) / \|\nabla f_{k-1}\|^2$

Polak-Ribiere is generally preferred because it automatically resets to steepest descent when progress stalls ($\beta_k \approx 0$ when consecutive gradients are similar). Fletcher-Reeves can generate poor search directions after an imprecise line search.

Neither variant retains the $n$-step finite termination property for nonquadratic objectives. Periodic restarts (every $n$ or $\sqrt{n}$ steps) are common.

Common Confusions

Watch Out

CG is not gradient descent with momentum

Both use information from previous iterations, but the mechanisms differ. CG chooses A-conjugate directions (which are optimal for quadratics). Momentum methods use an exponential moving average of past gradients. For a quadratic objective, CG converges in at most $n$ steps; gradient descent with momentum does not have this property.

Watch Out

CG requires symmetry and positive definiteness

The standard CG algorithm solves $Ax = b$ only when $A$ is SPD. For non-symmetric or indefinite systems, use GMRES, BiCGSTAB, or MINRES. For symmetric indefinite systems, MINRES is the appropriate Krylov method.

Summary

• CG solves $Ax = b$ for SPD $A$ using only matrix-vector products
• Converges in at most $n$ iterations (exact arithmetic), often much fewer
• Error bound depends on $\sqrt{\kappa}$: quadratically better than steepest descent
• Preconditioning reduces $\kappa$ and accelerates convergence
• Nonlinear CG (Fletcher-Reeves, Polak-Ribiere) extends to general smooth optimization

Exercises

ExerciseCore

Problem

A $1000 \times 1000$ SPD matrix has condition number $\kappa = 100$. Using the CG convergence bound, how many iterations are needed to reduce the $A$-norm error by a factor of $10^{-6}$?

Hint

Reveal Solution

ExerciseAdvanced

Problem

Prove that if $A$ has exactly $r$ distinct eigenvalues, CG converges in at most $r$ iterations (exact arithmetic).

Hint

Reveal Solution

References

Canonical:

• Hestenes & Stiefel, Methods of Conjugate Gradients for Solving Linear Systems (1952)
• Trefethen & Bau, Numerical Linear Algebra, Lectures 38-39
• Nocedal & Wright, Numerical Optimization, Chapter 5

Current:

• Shewchuk, An Introduction to the Conjugate Gradient Method Without the Agonizing Pain (1994), a widely used tutorial

Next Topics

Understanding CG provides the foundation for preconditioned iterative methods and second-order optimization in large-scale ML.