Conjugate Gradient Methods

The problem CG is built for

CG is not a general-purpose optimizer. It is specialized for a very important class of problems:

1) Solve a linear system

2) where A is symmetric positive-definite (SPD):

• •Symmetric: $A = A^\top$
• •Positive-definite: $\mathbf{z}^\top A\mathbf{z} > 0$ for all nonzero z

This single assumption (SPD) is what makes CG fast, stable, and geometric.

Quadratic optimization ↔ linear systems

The key equivalence is that solving $A\mathbf{x}=\mathbf{b}$ is the same as minimizing a strictly convex quadratic:

Compute the gradient (showing the work):

For symmetric $A$,

and

So

Setting $\nabla f(\mathbf{x}^*)=\mathbf{0}$ gives $A\mathbf{x}^* = \mathbf{b}$. So the minimizer and the linear-system solution coincide.

Why gradient descent isn’t enough

If you already know gradient descent, you might propose:

That works, but on ill-conditioned problems it “zig-zags” along narrow valleys (elliptical level sets). Convergence becomes painfully slow because the best step size is constrained by the largest eigenvalue while progress is governed by the smallest.

CG fixes this by changing the question from “which way is steepest right now?” to:

“Which direction lets me make progress without undoing what I already fixed?”

Two geometric lenses (keep both in mind)

CG is easiest to understand using two complementary pictures:

1) Residual view (linear algebra):

• •Define the residual rₖ as

• •Note that $\mathbf{r}_k = -\nabla f(\mathbf{x}_k)$.

2) Energy/inner-product view (geometry):

• •SPD $A$ defines an inner product

• •And an induced norm $\|\mathbf{u}\|_A = \sqrt{\mathbf{u}^\top A\mathbf{u}}$.

CG chooses directions that are orthogonal under this $A$-inner-product (called A-conjugate), which is the right notion of “independent directions” for this quadratic.

Anchor visualization (static SVG)

Suggested SVG anchor for this node: /images/cg-geometry.svg

What it should depict (single figure with labeled layers):

• •2D ellipsoidal level sets of $f(\mathbf{x})$
• •iterates x₀, x₁, x₂…
• •residual vectors r₀, r₁ at those points (steepest descent directions)
• •search directions p₀, p₁, p₂ (CG directions)
• •annotations showing:
• •rᵢ ⟂ rⱼ (Euclidean orthogonality of residuals)
• •pᵢ ⟂ₐ pⱼ (A-orthogonality / conjugacy)

Even if learners skip interactivity, this one image should connect “ellipses” ↔ “residuals” ↔ “conjugate directions.”

The key idea: don’t reuse the same curvature direction

For the quadratic $f(\mathbf{x})=\tfrac12\mathbf{x}^\top A\mathbf{x}-\mathbf{b}^\top\mathbf{x}$, curvature is encoded by $A$. If you move along direction p, the second derivative along that line is $\mathbf{p}^\top A\mathbf{p}$. So if two directions are “independent with respect to curvature,” they should satisfy:

This is A-conjugacy (also called A-orthogonality).

Why it matters: if you minimize the quadratic along p₀ and then later move along some other direction p₁ that is A-conjugate to p₀, you do not spoil the minimization you already achieved along p₀ in the quadratic’s geometry.

CG as “Gram–Schmidt, but with A in the middle”

In Euclidean space, we make directions orthogonal using dot products $\mathbf{u}^\top\mathbf{v}$. For CG, the orthogonality condition is instead $\mathbf{u}^\top A\mathbf{v}$.

Conceptually, CG is building an A-orthogonal basis of directions

and doing an exact line minimization along each.

The Krylov subspace view (the “growth” story)

A second major structural fact is that CG iterates live in expanding Krylov subspaces.

Let the initial residual be r₀ = b − Ax₀.

Define the Krylov subspace of order k:

CG chooses xₖ from the affine space

This matters because:

• •Each iteration adds one new dimension of freedom.
• •CG chooses the best point in that expanding space under the quadratic objective.

This is the conceptual reason for the celebrated finite-step property:

In exact arithmetic, CG finds the exact solution in at most n iterations (n = dimension), because after n steps the Krylov subspace can span the whole space.

In practice, roundoff and “near dependence” mean you often stop earlier based on tolerance.

Residual orthogonality (a different orthogonality)

A surprising and useful property (for exact arithmetic) is:

Residuals are mutually orthogonal in the standard dot product, even though search directions are A-orthogonal.

This is one of the best places to use visualization.

Interactive visualization plan (explicitly targets the weak spot)

Interactive canvas suggestion (two linked panels):

Panel A: Level sets + directions (geometry)

• •Show ellipsoidal contours of $f$ in 2D.
• •Display xₖ, rₖ, and pₖ.
• •Toggle overlays:

1) “Residual orthogonality”: draw r₀, r₁, r₂ from a common origin and show right-angle markers for rᵢ ⟂ rⱼ.

2) “A-conjugacy”: draw pᵢ and show right-angle markers in the A-metric (see below).

Panel B: Krylov subspace growth (algebra)

• •Show basis vectors being added:
• •start with span{r₀}
• •then span{r₀, Ar₀}
• •then add A²r₀, etc.
• •In 2D, show the subspace as a line that rotates/expands into the plane.
• •Show that xₖ is the minimizer restricted to that subspace (a projection onto the solution in the A-norm).

How to visualize A-orthogonality concretely:

A-orthogonality is Euclidean orthogonality after a change of variables. Since A is SPD, write $A = C^\top C$ (e.g., Cholesky). Then

So pᵢ ⟂ₐ pⱼ means transformed vectors Cpᵢ are orthogonal in Euclidean space. The canvas can show both p and Cp.

A small but important mental model

CG is not “steepest descent with a fancy step size.” It is closer to:

• •Build a low-dimensional model space where the solution lives (Krylov)
• •Find the best point in that space (exact minimization)
• •Expand the space by one dimension with a direction that is independent under A

That combination produces fast progress on SPD problems.

Start from the simplest objects

We track three sequences:

• •iterate xₖ (current solution estimate)
• •residual rₖ = b − Axₖ (negative gradient)
• •search direction pₖ (A-conjugate directions)

Initialize:

The first direction equals steepest descent. After that, directions become “corrected” to maintain conjugacy.

Step 1: exact line search to get αₖ

We move along pₖ:

Choose $\alpha_k$ to minimize $f(\mathbf{x}_k + \alpha\mathbf{p}_k)$.

Compute it explicitly. Let $\phi(\alpha)=f(\mathbf{x}_k+\alpha\mathbf{p}_k)$.

Start expanding:

Differentiate w.r.t. α:

Group terms:

Recall $A\mathbf{x}_k-\mathbf{b}=-\mathbf{r}_k$:

Set $\phi'(\alpha)=0$:

In standard CG (with pₖ built from rₖ), one can show $\mathbf{p}_k^\top\mathbf{r}_k = \mathbf{r}_k^\top\mathbf{r}_k$, giving the commonly used formula:

Step 2: update residual cheaply

Rather than recompute rₖ₊₁ = b − Axₖ₊₁ from scratch, use:

So

This highlights the computational core: each iteration needs one matrix-vector product $A\mathbf{p}_k$.

Step 3: build the next search direction via βₖ

We want pₖ₊₁ to be in the span of the new residual plus the previous direction, but also maintain A-conjugacy:

Choose $\beta_k$ so that

Compute:

So

With additional CG identities, this becomes the well-known Fletcher–Reeves-style expression:

(For SPD quadratics and exact arithmetic, these are consistent. In numerical practice, variants exist.)

Putting it together (vanilla CG algorithm)

Given SPD A, b, and x₀:

1) r₀ = b − Ax₀

2) p₀ = r₀

3) For k = 0, 1, 2, … until converged:

• •$\alpha_k = (\mathbf{r}_k^\top\mathbf{r}_k)/(\mathbf{p}_k^\top A\mathbf{p}_k)$
• •xₖ₊₁ = xₖ + αₖ pₖ
• •rₖ₊₁ = rₖ − αₖ Apₖ
• •$\beta_k = (\mathbf{r}_{k+1}^\top\mathbf{r}_{k+1})/(\mathbf{r}_k^\top\mathbf{r}_k)$
• •pₖ₊₁ = rₖ₊₁ + βₖ pₖ

What convergence means here

CG is typically stopped when the residual is small:

Because rₖ is the gradient (up to sign), small residual also means near-stationary for the quadratic.

A note on complexity and when CG shines

CG is attractive when:

• •A is large and sparse
• •you can compute $A\mathbf{v}$ fast without forming dense matrices
• •you cannot afford $O(n^3)$ factorization

Per iteration cost is dominated by one $A\mathbf{p}_k$ product plus a handful of dot products and axpy operations.

Memory is small: a few vectors (often 4–6) of length n.

This is why CG is the workhorse for PDE discretizations, large least-squares normal equations (carefully), and many kernel / Gaussian process subproblems.

Why preconditioning exists

In exact arithmetic, CG solves in ≤ n steps. So why do people obsess over convergence speed?

Because in realistic problems:

• •n may be millions
• •roundoff means you won’t see exact n-step termination
• •you stop when residual is “good enough,” long before n

So the practical question becomes:

How many iterations to reach tolerance?

This is governed largely by the condition number

When κ(A) is large, the ellipsoids are extremely stretched, and CG must work hard to correct error components along small-eigenvalue directions.

A classic bound (informal but insightful) is that the A-norm error decreases like:

You don’t need to memorize this; the message is:

• •convergence depends on √κ
• •reducing κ can massively reduce iterations

Preconditioning: change the problem without changing the solution

A preconditioner is an SPD matrix M that approximates A but is easier to invert (or solve with).

We want to solve

but instead run CG on a transformed system with better conditioning.

There are multiple equivalent views.

View 1: Left-preconditioning

Solve

This does not preserve symmetry in the usual Euclidean inner product, but it can be handled carefully (see View 3).

View 2: Symmetric preconditioning (the clean SPD story)

Use $M = C^\top C$ (Cholesky of M). Define $\mathbf{y} = C\mathbf{x}$.

Start with:

Substitute $\mathbf{x}=C^{-1}\mathbf{y}$:

Multiply by $C^{-\top}$:

Now the matrix

is SPD, and hopefully has a much smaller condition number than A.

Run standard CG on $\tilde{A}\mathbf{y}=\tilde{\mathbf{b}}$ and map back $\mathbf{x}=C^{-1}\mathbf{y}$.

View 3: “Preconditioned inner product” (most intuitive)

Preconditioned CG (PCG) can be seen as running CG where the notion of orthogonality is changed by M.

Instead of using raw residual rₖ, define a “preconditioned residual” zₖ:

Think: zₖ is what you get after cheaply solving $M\mathbf{z}_k=\mathbf{r}_k$.

Then PCG uses dot products like $\mathbf{r}_k^\top\mathbf{z}_k$ to measure progress.

PCG algorithm (practical form)

Given SPD A, SPD M, b, x₀:

• •r₀ = b − Ax₀
• •Solve $M\mathbf{z}_0 = \mathbf{r}_0$
• •p₀ = z₀

For k = 0, 1, 2, …:

• •$\alpha_k = \dfrac{\mathbf{r}_k^\top\mathbf{z}_k}{\mathbf{p}_k^\top A\mathbf{p}_k}$
• •xₖ₊₁ = xₖ + αₖ pₖ
• •rₖ₊₁ = rₖ − αₖ Apₖ
• •Solve $M\mathbf{z}_{k+1} = \mathbf{r}_{k+1}$
• •$\beta_k = \dfrac{\mathbf{r}_{k+1}^\top\mathbf{z}_{k+1}}{\mathbf{r}_k^\top\mathbf{z}_k}$
• •pₖ₊₁ = zₖ₊₁ + βₖ pₖ

Compared to CG, you add one extra operation: apply $M^{-1}$ (implemented as solving with M).

What makes a “good” preconditioner?

A good M satisfies a three-way tradeoff:

Common choices:

• •Jacobi / diagonal: $M = \mathrm{diag}(A)$ (cheap, modest gains)
• •Incomplete Cholesky (IC): sparse approximation of Cholesky (often excellent for PDEs)
• •SSOR / Gauss–Seidel variants: classical iterative preconditioners
• •Multigrid preconditioners: often near-optimal for elliptic PDEs

Visualization: preconditioning as “making ellipses rounder”

Preconditioning’s geometric effect is easier than its algebra:

• •Original problem has level sets shaped by A (possibly very elongated ellipses).
• •After the change of variables induced by M, those ellipses become closer to circles.

Interactive canvas suggestion (add-on mode):

• •Show original ellipsoids in x-space.
• •Let user choose M (diagonal, IC(0), or a toy exact M=A).
• •Show transformed coordinates y where $\tilde{A}=C^{-\top}AC^{-1}$.
• •Show how κ changes and how many CG steps are needed.

This makes the “conditioning → iteration count” link tangible.

Connection to least squares and normal equations (careful)

Many problems arrive as least squares:

Normal equations yield SPD system:

You can run CG on $A=B^\top B$, but note:

• •conditioning squares: $\kappa(B^\top B)=\kappa(B)^2$
• •can be significantly worse numerically

Often better: use methods like LSQR / LSMR (Krylov methods tailored to least squares), or precondition aggressively.

What you should leave with

CG is a bridge between optimization and linear algebra:

• •optimization perspective: exact line searches on a quadratic in cleverly chosen directions
• •linear algebra perspective: Krylov method minimizing error in A-norm
• •numerical perspective: preconditioning is everything