NLA MT25, Conjugate gradient method

Set $r _ 0 := b - A x _ 0$ (initial residual) and look for an approximate solution in the affine Krylov space $x _ 0 + \mathcal K _ k(A, r _ 0)$, i.e. $x = x _ 0 + Q _ k y$ where $\mathrm{range}(Q _ k) = \mathcal K _ k(A, r _ 0)$, and fix it by Galerkin orthogonality: force the residual $Ax - b$ to be orthogonal to $\mathrm{range}(Q _ k)$,

which using $A x _ 0 - b = -r _ 0$ rearranges to

With the ∆lanczos-decomposition $AQ _ k = Q _ k T _ k + q _ {k+1} [0, \ldots, 0, t _ {k+1, k}]$ giving $Q _ k^\top A Q _ k = T _ k$ (symmetric tridiagonal), this becomes

Why this is a good idea:

• For $A \succ 0$, Galerkin orthogonality is equivalent to minimising the $A$-norm error $\|x - x _ \ast\| _ A$ over $x _ 0 + \mathcal K _ k(A, r _ 0)$ (∆conjugate-gradient-motivation): every step returns the best approximation in a subspace that grows by one dimension, so the error never increases and the exact solution is reached in $\le n$ steps.
• Building $\mathcal K _ k(A, r _ 0)$ needs only matrix-vector products with $A$ — no $O(n^3)$ factorisation — so it is feasible for the large sparse $A$ on which direct methods are not.
• The projected system $T _ k y = Q _ k^\top r _ 0$ is small and tridiagonal, and the Lanczos three-term recurrence keeps memory at $O(n)$; convergence is a polynomial-approximation problem on $\lambda(A)$, hence fast (often $O(1)$ iterations) when the spectrum is clustered or well-conditioned (∆conjugate-gradient-convergence).

The lecture’s $x _ 0 = 0$ specialisation drops the affine offset: then $r _ 0 = b$, $\mathcal K _ k(A, r _ 0) = \mathcal K _ k(A, b)$, and $T _ k y = Q _ k^\top b$.

Given an initial guess $x _ 0$ and the corresponding initial residual $r _ 0 := b - A x _ 0$, the iterate $x = x _ 0 + Q _ k y$ (with $\mathrm{range}(Q _ k) = \mathcal K _ k(A, r _ 0)$) satisfies the Galerkin condition

equivalently

i.e. the residual $b - Ax$ is orthogonal to $\mathrm{range}(Q _ k) = \mathcal K _ k(A, r _ 0)$.

Why it matters:

• Selects a unique iterate: “$x \in x _ 0 + \mathcal K _ k(A, r _ 0)$” is $k$-dimensional and underdetermined; Galerkin supplies exactly $k$ scalar conditions for the $k$ unknowns $y$, collapsing the $n \times n$ problem to the square system $T _ k y = Q _ k^\top r _ 0$ — tridiagonal via the ∆lanczos-decomposition, so $O(k)$ to solve.
• Optimality in disguise: for $A \succ 0$, residual orthogonality to $\mathcal K _ k(A, r _ 0)$ is precisely the inner-product best-approximation condition, equivalent to $x = \mathrm{argmin} _ {x \in x _ 0 + \mathcal K _ k(A, r _ 0)} \|x - x _ \ast\| _ A$ (∆conjugate-gradient-motivation). So CG’s iterate is the best point in the affine subspace, not merely a feasible one.
• The projection-method template: the same “residual $\perp$ trial subspace” idea underlies GMRES/MINRES (Petrov-Galerkin: orthogonal to a different test space) and Rayleigh-Ritz for eigenproblems (∆rayleigh-ritz-algorithm); it originates in the finite element method.

The $x _ 0 = 0$ specialisation replaces $r _ 0$ with $b$ and the affine space with $\mathcal K _ k(A, b)$ throughout.

The CG algorithm minimises the $A$-norm of the error over the affine Krylov space,

Writing $x = x _ 0 + Q _ k y$, and using $A(x _ \ast - x _ 0) = b - A x _ 0 = r _ 0$, we have $x - x _ \ast = -(x _ \ast - x _ 0) + Q _ k y$ and so

using $(x _ \ast - x _ 0)^\top A = r _ 0^\top$ and $A$ symmetric. Differentiating in $y$, the minimum is at

which is exactly $Q _ k^\top (A Q _ k y - r _ 0) = 0$.

Another way to see this is via the results in Notes - Numerical Analysis HT24, Best approximation in inner product spaces^U: the minimum of $\|x - x _ \ast\| _ A$ over the affine space $x _ 0 + \mathrm{range}(Q _ k)$ is attained when $\langle q _ i, x _ \ast - x \rangle _ A = 0$ for every column $q _ i$ of $Q _ k$.

The $x _ 0 = 0$ specialisation replaces $r _ 0$ with $b$ and $x _ \ast - x _ 0$ with $x _ \ast$ throughout.

For a symmetric positive definite matrix $A \in \mathbb R^{n \times n}$, $A$ induces an inner product $\langle u, v \rangle _ A = u^\top A v$. Two vectors $u, v$ are $A$-conjugate when they’re orthogonal in this inner product, so that

Equivalently, $u$ and $v$ are $A$-conjugate if $A^{1/2} u$ and $A^{1/2} v$ are orthogonal in the ordinary sense.

Given $x _ 0 \in \mathbb R^n$ (take $x _ 0 = 0$ if no better guess is available), define $r _ 0 := b - A x _ 0$. The conjugate gradient method computes at each step $k$ an iterate $x _ k$ that is the unique element of the affine Krylov space $x _ 0 + \mathcal K _ k(A, r _ 0) = x _ 0 + \text{span}(r _ 0, A r _ 0, \ldots, A^{k-1} r _ 0)$ satisfying any of the following equivalent characterisations:

• ∆galerkin-condition: $Q _ k^\top (Ax _ k - b) = 0$, where $\mathrm{range}(Q _ k) = \mathcal K _ k(A, r _ 0)$ is the Lanczos basis. This says that the residual is orthogonal to the Krylov subspace.
• $A$-norm minimisation over the affine Krylov space: $x _ k = \text{argmin} _ {x \in x _ 0 + \mathcal K _ k(A, r _ 0)}\|x - x _ \ast\| _ A$ where $x _ \ast = A^{-1} b$.
• Tridiagonal projection: $x _ k = x _ 0 + Q _ k y$ with $T _ k y = Q _ k^\top r _ 0$, where $T _ k = Q _ k^\top A Q _ k$ is the symmetric tridiagonal Lanczos projection.

Algorithm:

• Initialisation:
  • $x _ 0 \in \mathbb R^n$ (arbitrary initial guess; take $x _ 0 = 0$ if none is available).
  • $r _ 0 = b - Ax _ 0$ (initial residual).
  • $p _ 0 = r _ 0$ (initial search direction).
• For $k = 0, 1, 2, \ldots$:
  • $\alpha _ k = \frac{\langle r _ k, r _ k \rangle}{\langle p _ k, Ap _ k \rangle}$ (step size, minimises $\|x - x _ \ast\| _ A$ along $p _ k$ from $x _ k$).
  • $x _ {k+1} = x _ k + \alpha _ k p _ k$ (update iterate).
  • $r _ {k+1} = r _ k - \alpha _ k A p _ k$ (update residual, coming directly from $r _ {k+1} = b - Ax _ {k+1}$).
  • $\beta _ k = \frac{\langle r _ {k+1}, r _ {k+1} \rangle}{\langle r _ k, r _ k \rangle}$ (Gram-Schmidt-like coefficient making the next $p _ {k+1}$ $A$-conjugate to $p _ k$, and in fact enforcing $A$-conjugacy with all previous $p _ j$, see ∆conjugate-gradient-maintains-conjugacy).
  • $p _ {k+1} = r _ {k+1} + \beta _ k p _ k$ (new search direction).
• …stop when $\|r _ k\|$ is small enough.

The lecture’s $x _ 0 = 0$ specialisation gives $r _ 0 = b$, so $\mathcal K _ k(A, r _ 0) = \mathcal K _ k(A, b)$ and the affine offset disappears.

Time:

• 1 matrix-vector product $Ap _ k$, $O(n^2)$
• 2 inner products $\langle r _ k, r _ k\rangle$, and $\langle p _ k, Ap _ k \rangle$, $O(n)$
• 3 axpys ($Ax + y$), $O(n)$
• Total per step: $O(n^2)$.

Memory: $O(n)$

CG minimises the $A$-norm of the error $\|x-x _ \ast\| _ A$ over the affine Krylov space $x _ 0 + \mathcal K _ k(A, r _ 0)$, where $x _ \ast = A^{-1} b$ and $r _ 0 = b - A x _ 0$. The consequence of $A$-conjugacy means that if the search directions $p _ 0, p _ 1, \ldots, p _ {k-1}$ are mutually $A$-conjugate, then minimising $\|x - x _ \ast\| _ A$ along each direction is independent. Hence the iterate

with each $\alpha _ i$ the 1D $A$-norm-minimum along $p _ i$, is the exact minimiser of $\|x - x _ \ast\| _ A$ over $x _ 0 + \text{span}(p _ 0, \ldots, p _ {k-1})$, and hence after $k = n$ steps this is the exact solution to $Ax = b$. So there is no need to revisit the previous directions.

Setup: Let $A$ be symmetric positive definite, $x _ 0$ an arbitrary initial guess, $r _ 0 = b - A x _ 0$, and $r _ k = b - Ax _ k$ the residual at step $k$. We have the following CG identities:

1. $r _ k = (-1)^k \|r _ k\| _ 2 q _ {k+1}$ for the Lanczos basis $\{q _ j\}$ of $\mathcal K _ k(A, r _ 0)$ (i.e. the CG residuals are scaled Lanczos vectors).
2. $\langle r _ i, r _ j \rangle = 0$ for $i \ne j$ (the residuals are mutually orthogonal)
3. $T _ k = Q _ k^\top A Q _ k$ is tridiagonal (∆arnoldi-tridiagonal-when-a-symmetric), so $\langle q _ i, q _ j\rangle _ A = 0$ whenever $ \vert i - j \vert > 1$.
4. The search directions satisfy $\text{span}(p _ 0, \ldots, p _ k) = \text{span}(r _ 0, \ldots, r _ k) = \mathcal K _ {k+1}(A, r _ 0)$ (induction on $p _ {k+1} = r _ {k+1} + \beta _ k p _ k$).

Proof: (by induction on $k$).

Base case $k = 0$: this is vacuous.

Inductive step: Assume $\langle p _ i, p _ j \rangle _ A = 0$ for all $i \ne j$ with $i, j \le k$. We show $\langle p _ {k+1}, p _ j \rangle _ A = 0$ for all $j \le k$.

• For $j = k$, this is enforced by the construction. The choice $\beta _ k = \langle r _ {k+1}, r _ {k+1} \rangle / \langle r _ k, r _ k \rangle$ is the Gram-Schmidt coefficient that makes $p _ {k+1} = r _ {k+1} + \beta _ k p _ k$ satisfy

(using the identity $\langle r _ {k+1}, p _ k \rangle _ A = -\langle r _ {k+1}, r _ {k+1} \rangle / \alpha _ k$ and $\langle p _ k, p _ k \rangle _ A = \langle r _ k, r _ k \rangle / \alpha _ k$ from the definitions of $\alpha _ k$ and the residual recurrence).

• For $j < k$, we expand

The second term vanishes by the inductive hypothesis, so that $\langle p _ k, p _ j \rangle _ A = 0$ since $j < k$. For the first term, $Ap _ {j} \in A \cdot \mathcal K _ {j+1}(A, r _ 0) \subseteq \mathcal K _ {j+2}(A, r _ 0) = \text{span}(r _ 0, \ldots, r _ {j+1})$ by (4). Since $j+1 \le k$, every $r _ i$ in this span has $i \le k$, and $r _ {k+1}$ is orthogonal to all of them by (2):

So $\langle r _ {k+1}, p _ j\rangle _ A = r^\top _ {k+1} A p _ j = 0$, and the first term vanishes too.

Hence $\langle p _ {k+1}, p _ j \rangle _ A = 0$ for all $j \le k$.

We prove three separate subclaims:

• (i) Residual orthogonality: $\langle r _ i, r _ j \rangle = 0$ for $i \ne j$
• (ii) $A$-conjugacy of search directions: $\langle p _ i, Ap _ j \rangle = 0$ for $i \ne j$
• (iii) Subspace identity: $\text{span}(p _ 0, \ldots, p _ {k-1}) = \text{span}(r _ 0, \ldots, r _ {k-1}) = \mathcal K _ k(A, r _ 0)$

Once these are proved, the main result follows since:

• By (ii) and (iii), $\{p _ 0, \ldots, p _ {k-1}\}$ is an $A$-conjugate basis of $\mathcal K _ k(A, r _ 0)$,
• For any $A$-conjugate basis, the unique $A$-norm minimiser over the affine space decomposes into 1D minima along each direction, and
• The algorithm’s $\alpha _ j$ is the 1D $A$-norm-error coefficient.

The fact error decreases monotonically and that the algorithm converges to $x _ \ast$ in at most $n$ steps is immediate from this main claim.

Step 1: Simultaneous induction proving (i) and (ii).

We prove by induction on $k$ the inductive hypothesis “$\mathcal H _ k$” that for all $0 \le j < k$,

• $\langle r _ k, r _ j \rangle = 0$,
• $\langle r _ k, p _ j \rangle = 0$,
• $\langle p _ k, Ap _ j \rangle = 0$.

Base case $\mathcal H _ 1$: prove $\langle r _ 1, r _ 0 \rangle = \langle r _ 1, p _ 0 \rangle = \langle p _ 1, Ap _ 0 \rangle = 0$.

Inductive step: Assume $\mathcal H _ k$ holds, now we aim to prove

• $\langle r _ {k+1}, r _ j \rangle = 0$,
• $\langle r _ {k+1}, p _ j \rangle = 0$,
• $\langle p _ {k+1}, Ap _ j \rangle = 0$.

for all $j = 0, \ldots, k$.

To see $\langle r _ {k+1}, p _ j \rangle = 0$, we split into two cases.

• Case $j = k$: Then $\langle r _ {k+1}, p _ k \rangle = \langle r _ k, p _ k \rangle - \alpha _ k \langle p _ k, Ap _ k \rangle$. Compute

from applying $\mathcal H _ k$’s identity $\langle r _ k, p _ {k-1}\rangle = 0$. With $\alpha _ k \langle p _ k, Ap _ k\rangle = \langle r _ k, r _ k \rangle$, this evaluates to $0$.

• Case $j < k$: Then $\langle r _ {k+1}, p _ j\rangle = \langle r _ k, p _ j \rangle$, and both terms are zero by $\mathcal H _ k$.

To see $\langle r _ {k+1}, r _ j \rangle = 0$, note that for any $j \le k$, $r _ j = p _ j - \beta _ {j-1} p _ {j-1}$ (with $\beta _ {-1} = 0$). Then

by the previous part.

To see $\langle p _ {k+1}, Ap _ j \rangle = 0$, we split into two cases.

• Case $j = k$: This is the defining property of $\beta _ k$. Substitute $p _ {k+1} = r _ {k+1} + \beta _ k p _ k$ and use $Ap _ k = (r _ k - r _ {k+1})/\alpha _ k$, together with the residual orthogonalities from the previous $\langle r _ {k+1}, r _ j \rangle = 0$ calculation:

With $\beta _ k = \langle r _ {k+1}, r _ {k+1} \rangle / \langle r _ k, r _ k \rangle$ and $\alpha _ k \langle p _ k, Ap _ k \rangle = \langle r _ k, r _ k \rangle$:

• Case $j < k$: Then $\langle p _ {k+1}, Ap _ j \rangle = \langle r _ {k+1}, Ap _ j \rangle + \beta _ k \langle p _ k, Ap _ j \rangle$. The second term vanishes by $\mathcal H _ k$. For the first, substitute $Ap _ j = (r _ j - r _ {j+1})/\alpha _ j$:

by the residual orthogonalities from the previous $\langle r _ {k+1}, r _ j \rangle = 0$ calculation, since both $j, j+1 \le k$.

This completes the induction, so that (i) and (ii) hold for all $k$.

Step 2: subspace identity (iii)

By induction on $k$, we aim to prove $\text{span}(p _ 0, \ldots, p _ {k-1}) = \text{span}(r _ 0, \ldots, r _ {k-1}) = \mathcal K _ k(A, r _ 0)$.

Base $k = 1$: $p _ 0 = r _ 0$, so all three are $\text{span}(r _ 0) = \mathcal K _ 1(A, r _ 0)$.

Inductive step: Assume the identity holds for $k$.

To see $\text{span}(r _ 0, \ldots, r _ k) = \mathcal K _ {k+1} (A, r _ 0)$, note that:

• $r _ k = r _ {k-1} - \alpha _ {k-1} A p _ {k-1}$ with $r _ {k-1} \in \mathcal K _ k(A, r _ 0)$ and $Ap _ {k-1} \in A \cdot \mathcal K _ k(A, r _ 0) \subseteq \mathcal K _ {k+1} (A, r _ 0)$, so $r _ k \in \mathcal K _ {k+1}(A, r _ 0)$.
• The $r _ 0, \ldots, r _ k$ are mutually orthogonal by (i). Assuming the algorithm has not terminated (else we are done), they are nonzero, and hence linearly independent. So $\dim \text{span}(r _ 0, \ldots, r _ k) = k+1$.
• $\dim \mathcal K _ {k+1}(A, r _ 0) \le k+1$. We have both an inclusion and a dimension match, so there is equality.

To see $\text{span}(p _ 0, \ldots, p _ k) = \text{span}(r _ 0, \ldots, r _ k)$, note that $p _ k = r _ k + \beta _ {k-1} p _ {k-1}$, so each side is contained in the other.

Step 3: the variational characterisation

The minimisation $x _ k = \text{argmin} _ {x \in x _ 0 + \mathcal K _ k(A, r _ 0)} \|x - x _ \ast\| _ A$ is, after the shift $\tilde x := x - x _ 0$, equivalent to

where $\tilde x _ \ast := x _ \ast - x _ 0$, so $A \tilde x _ \ast = b - A x _ 0 = r _ 0$.

We have established that $\{p _ 0, \ldots, p _ {k-1}\}$ is an $A$-conjugate basis of $\mathcal K _ k(A, r _ 0)$. By ∆a-conjugacy, the unique minimiser decomposes uniquely as

using $A \tilde x _ \ast = r _ 0$.

It remains to verify that the algorithm’s iterate $x _ k - x _ 0 = \sum^{k-1} _ {j = 0} \alpha _ j p _ j$ has $\alpha _ j = \gamma _ j$, equivalently that $\langle p _ j, r _ 0 \rangle = \langle r _ j, r _ j \rangle$.

Claim: $\langle p _ j, r _ 0 \rangle = \langle r _ j, r _ j \rangle$.

Proof: Iterating the residual recurrence $r _ {i+1} = r _ i - \alpha _ i A p _ i$ from $i = 0$ up to $i = j - 1$,

Taking the inner product with $p _ j$,

by the $A$-conjugacy claim in (ii). For the remaining term, $p _ j = r _ j + \beta _ {j-1} p _ {j-1}$ gives

by the claim in (i) that $\langle r _ j, p _ {j-1} \rangle = 0$.

Hence

so $x _ k - x _ 0 = \sum _ j \alpha _ j p _ j = \tilde x _ k$, i.e. $x _ k$ is the variational minimiser.

Step 4: monotonicity and termination.

The error decreases monotonically, $\|x _ {k+1} - x _ \ast\| _ A \le \|x _ k - x _ \ast\| _ A$, since $x _ 0 + \mathcal K _ k(A, r _ 0) \subseteq x _ 0 + \mathcal K _ {k+1}(A, r _ 0)$, so $x _ {k+1}$ is the minimiser over a larger affine subspace containing $x _ k$.

The algorithm also converges to $x _ \ast$ in at most $n$ steps in exact arithmetic, so that $x _ n = x _ \ast$, by the following: $\dim \mathcal K _ n(A, r _ 0) \le n$. If $\mathcal K _ n(A, r _ 0) = \mathbb R^n$, then $x _ 0 + \mathcal K _ n(A, r _ 0) = \mathbb R^n$, so $x _ n = \text{argmin} _ {x \in \mathbb R^n} \| x - x _ \ast \| _ A = x _ \ast$. If $\mathcal K _ n(A, r _ 0) \subsetneq \mathbb R^n$, then by (i), (ii) and (iii) the algorithm must have already terminated at some step $k < n$ with $r _ k = 0$, hence $A x _ k = b$, hence $x _ k = x _ \ast$.

Let $e _ k := x _ \ast - x _ k$ be the error after the $k$th CG iteration where $x _ \ast$ is the exact solution; in particular $e _ 0 = x _ \ast - x _ 0$. Then

The lecture’s $x _ 0 = 0$ specialisation gives $e _ 0 = x _ \ast$.

Setup: Because $A$ is symmetric, the spectral theorem (∆symmetric-eigenvalue-decomposition) gives $A = V \Lambda V^\top$ with $V$ orthogonal ($V^\top V = V V^\top = I$) and $\Lambda = \text{diag}(\lambda _ 1, \ldots, \lambda _ n)$ — symmetry alone makes $A$ orthogonally diagonalisable with real eigenvalues. Positive-definiteness then forces every $\lambda _ i > 0$, which is also what makes $\|z\| _ A := \sqrt{z^\top A z}$ a norm. Write $w := V^\top e _ 0$. Then:

so $\|e _ k\| _ A \le 2\, \|e _ 0\| _ A \big(\tfrac{\sqrt{\kappa _ 2(A)} - 1}{\sqrt{\kappa _ 2(A)} + 1}\big)^k$. $\square$

The steps:

• $(\star 1)$ CG optimality. By the variational characterisation of CG (∆conjugate-gradient-correctness-proof), $x _ k$ is the minimiser of $\|x - x _ \ast\| _ A$ over $x _ 0 + \mathcal K _ k(A, r _ 0)$, so $\|e _ k\| _ A = \min _ {x \in x _ 0 + \mathcal K _ k(A, r _ 0)} \|x - x _ \ast\| _ A$; divide by the fixed constant $\|e _ 0\| _ A$ ($x$ is now the minimisation variable, not the iterate $x _ k$).
• $(\star 2)$ Krylov $\to$ polynomial. Every $x \in x _ 0 + \mathcal K _ k(A, r _ 0)$ has the form $x = x _ 0 + p _ {k-1}(A) r _ 0$ for some $p _ {k-1}$ of degree $\le k-1$ (∆bite-krylov-polynomial-form). Then $x - x _ \ast = p _ {k-1}(A) r _ 0 - (x _ \ast - x _ 0) = p _ {k-1}(A) r _ 0 - e _ 0$, and the minimisation over $x$ becomes one over $p _ {k-1} \in \mathcal P _ {k-1}$.
• $(\star 3)$ Reindex the polynomial. Substitute $r _ 0 = A e _ 0$: $p _ {k-1}(A) r _ 0 - e _ 0 = \big(p _ {k-1}(A) A - I\big) e _ 0 = -p(A) e _ 0$ where $p(z) := 1 - z\, p _ {k-1}(z)$ (polynomials in $A$ commute). The sign does not affect the norm, and as $p _ {k-1}$ ranges over $\mathcal P _ {k-1}$, $p$ ranges over exactly the degree-$\le k$ polynomials with $p(0) = 1$.
• $(\star 4)$ Diagonalise (uses symmetry). $p(A) = V\, \text{diag}\big(p(\lambda _ 1), \ldots, p(\lambda _ n)\big)\, V^\top$. For any $z$, $\|z\| _ A^2 = z^\top A z = \sum _ i \lambda _ i (V^\top z) _ i^2$; with $z = p(A)\, e _ 0$ this gives $(V^\top z) _ i = p(\lambda _ i)\, w _ i$, and with $z = e _ 0$ it gives $(V^\top z) _ i = w _ i$ — hence the two sums.
• $(\star 5)$ Bound by the worst eigenvalue. Each weight $\lambda _ i w _ i^2 \ge 0$ (positive-definiteness, $\lambda _ i > 0$), so $\sum _ i \lambda _ i\, p(\lambda _ i)^2\, w _ i^2 \le \big(\max _ j p(\lambda _ j)^2\big) \sum _ i \lambda _ i\, w _ i^2$; take square roots — the $\big(\sum _ i \lambda _ i w _ i^2\big)^{1/2}$ cancels the denominator.
• $(\star 6)$ Relax to the interval. $\{\lambda _ j\} \subseteq [\lambda _ {\min}(A), \lambda _ {\max}(A)]$, and maximising over a larger set can only increase the value.
• $(\star 7)$ The assumed Chebyshev estimate (∆cg-chebyshev-bound-on-eigenvalue-interval).

For ease of notation, denote $a = \lambda _ {\min}$, $b = \lambda _ {\max}$, and $c := \tfrac{b + a}{b - a}$. Let

where $c _ k = 1 / T _ k(-c)$. Then $p(0) = 1$, and by ∆chebyshev-polynomials-growth ($ \vert T _ k \vert \le 1$ on $[-1, 1]$) we have

on $x \in [a, b]$, using the symmetry $ \vert T _ k(-y) \vert = \vert T _ k(y) \vert $ and the fact that $T _ k(c) > 0$ for $c > 1$ to drop the absolute value.

Switch to the $z$-parametrisation. From the Chebyshev definition,

We need to find the $z$ corresponding to $x = c$.

Solve for $z$. Multiplying $\tfrac 1 2 (z + z^{-1}) = c$ through by $2z$ gives the quadratic

Since $b > a > 0$ implies $c > 1$, both roots are real and positive. Their product is $z _ + z _ - = 1$ (the constant term of the monic quadratic), so they are reciprocals: $z _ - = 1/z _ +$, with $z _ + > 1 > z _ - > 0$.

Either root gives the same $T _ k(c)$. The map $z \mapsto \tfrac 1 2(z^k + z^{-k})$ is symmetric in $z \leftrightarrow z^{-1}$, so plugging in $z _ +$ or $z _ -$ yields the same value. We pick $z := z _ +$ (the root with $z > 1$) purely for convenience: that lets us write the clean lower bound $T _ k(c) \ge \tfrac 1 2 z^k$ below.

Simplify $z _ +$. With $c = \tfrac{b + a}{b - a}$,

so

Dividing numerator and denominator by $\sqrt a$ and using $\kappa _ 2(A) = b/a$ (SPD $A$),

$z > 1$ provided $\kappa _ 2(A) > 1$, i.e. $A$ is not a scalar multiple of $I$ (in the trivial case $\kappa _ 2(A) = 1$, both $z _ \pm = 1$ and the bound holds with equality, $ \vert p(x) \vert \le 1$).

Conclude the bound. Since $z > 1$,

so $1/T _ k(c) \le 2 z^{-k}$. Hence

When $\kappa _ 2(A)$ is large, so that there is a wide relative spread where $\lambda _ \min$ is much smaller than $\lambda _ \max$ in relative terms.

Setup: choose a matrix $M$ such that

and instead of solving $Ax = b$ directly, solve the preconditioned system

Since $M^\top A M$ is still symmetric positive definite (as long as $M$ is nonsingular and $A \succ 0$), CG can be applied to this modified system. For example, we could take $M = A^{-1/2}$.

We desire:

• $M$ approximates a square root of $A^{-1}$, e.g. perhaps we have an approximate Cholesky factorisation $A \approx R^\top R$, we may take $M = R^{-\top}$
• $M^\top A M$ is easy to multiply to a vector (e.g. perhaps sparse)
• $\kappa _ 2(M^\top A M)$ is small (e.g. $M^\top A M$ has clustered eigenvalues away from $0$)

All three pick out the same vector in $x _ 0 + \mathcal K _ k(A, r _ 0)$ — the first two say what the CG iterate is, the third says how it is actually computed:

• $A$-norm error minimiser in affine Krylov space: $x _ k = \arg\min _ {x \in x _ 0 + \mathcal K _ k(A, r _ 0)} \|x - x _ \ast\| _ A$.
• Galerkin orthogonality: $Q _ k^\top (A x _ k - b) = 0$ with $\mathrm{range}(Q _ k) = \mathcal K _ k(A, r _ 0)$ — the residual is orthogonal to the Krylov subspace.
• Tridiagonal projection: $x _ k = x _ 0 + Q _ k y$ with $T _ k y = Q _ k^\top r _ 0$, where $T _ k = Q _ k^\top A Q _ k$.

The lecture’s $x _ 0 = 0$ specialisation replaces $r _ 0$ with $b$ and the affine offset disappears.

Why the third matters (it is the computational one): substituting $x _ k = x _ 0 + Q _ k y$ into the Galerkin condition gives exactly $T _ k y = Q _ k^\top r _ 0$ (and minimising the $A$-norm error yields the same normal equations, ∆conjugate-gradient-motivation), so it is just Galerkin orthogonality rewritten in the $k$ unknowns $y$. Its point is efficiency: because $Q _ k$ is the Lanczos basis, $T _ k = Q _ k^\top A Q _ k$ is symmetric tridiagonal (∆lanczos-decomposition), so the projected $k \times k$ system costs only $O(k)$ to solve, and that tridiagonal structure is exactly what lets CG run as a short three-term recurrence in $O(n)$ memory instead of storing all of $Q _ k$. It is the Rayleigh-Ritz projection idea (∆rayleigh-ritz-algorithm) applied to $Ax = b$: project onto $\mathcal K _ k(A, r _ 0)$, solve the small compressed system, lift back via $x = x _ 0 + Q _ k y$.

Because CG minimises $\|x - x _ \ast\| _ A$ over the Krylov subspace, and $\|x\| _ A = \sqrt{x^\top A x}$ is only a norm when $A \succ 0$. For symmetric indefinite $A$, MINRES is used instead — it minimises the 2-norm of the residual rather than the $A$-norm of the error.

From the Chebyshev polynomial estimate. Setting $a = \lambda _ \min$, $b = \lambda _ \max$ and solving $\tfrac12(z + z^{-1}) = \tfrac{b+a}{b-a}$ for $z > 1$ gives $z = \tfrac{\sqrt{b/a} + 1}{\sqrt{b/a} - 1} = \tfrac{\sqrt\kappa + 1}{\sqrt\kappa - 1}$, and $T _ k(z) \ge \tfrac12 z^k$, so $1/T _ k(z) \le 2 z^{-k} = 2\left(\tfrac{\sqrt\kappa - 1}{\sqrt\kappa + 1}\right)^k$.

• Induct on $k$. Goal: $\langle p _ {k+1}, p _ j \rangle _ A = 0$ for all $j \le k$.
• Case $j = k$: enforced by the definition of $\beta _ k$. The choice $\beta _ k = \langle r _ {k+1}, r _ {k+1} \rangle / \langle r _ k, r _ k \rangle$ is exactly the Gram-Schmidt coefficient.
• Case $j < k$: expand $\langle p _ {k+1}, p _ j \rangle _ A = \langle r _ {k+1}, Ap _ j \rangle + \beta _ k \langle p _ k, Ap _ j \rangle$. Second term is zero by induction. For the first, observe $A p _ j \in \mathcal K _ {j+2}(A, r _ 0) = \mathrm{span}(r _ 0, \ldots, r _ {j+1})$, and $j + 1 \le k$, so $r _ {k+1}$ is orthogonal to every $r _ i$ in that span by residual orthogonality.

The whole argument runs on three structural facts: residual orthogonality, the four-term recurrence, and the Krylov-subspace identity $\mathrm{span}(p _ 0, \ldots, p _ j) = \mathcal K _ {j+1}(A, r _ 0)$.

Establish three facts via simultaneous / chained induction:

• (i) Residual orthogonality: $\langle r _ i, r _ j \rangle = 0$ for $i \ne j$.
• (ii) $A$-conjugacy of search directions: $\langle p _ i, A p _ j \rangle = 0$ for $i \ne j$ (compare ∆conjugate-gradient-maintains-conjugacy).
• (iii) Subspace identity: $\mathrm{span}(p _ 0, \ldots, p _ {k-1}) = \mathrm{span}(r _ 0, \ldots, r _ {k-1}) = \mathcal K _ k(A, r _ 0)$.

The main claim then follows from one-dimensional optimisation along each $A$-conjugate direction: an $A$-conjugate basis lets you decompose the $A$-norm minimisation over the affine space into independent 1D minima, and the algorithm’s $\alpha _ j$ is precisely that 1D coefficient. The third step verifies $\alpha _ j = \gamma _ j$ via the side-claim $\langle p _ j, r _ 0 \rangle = \langle r _ j, r _ j \rangle$ (derived by iterating the residual recurrence to write $r _ 0 = r _ j + A \sum _ {i < j} \alpha _ i p _ i$, then using $A$-conjugacy).

$\mathcal H _ k$: for all $0 \le j < k$,

• $\langle r _ k, r _ j \rangle = 0$ (residuals mutually orthogonal),
• $\langle r _ k, p _ j \rangle = 0$ (residual orthogonal to earlier search directions),
• $\langle p _ k, A p _ j \rangle = 0$ (search directions mutually $A$-conjugate).

All three are needed simultaneously: the inductive step for any one of the three identities at level $k+1$ uses the other two at level $k$.

• The whole optimisation framing assumes $\|\cdot\| _ A$ is a norm, which requires $A \succ 0$ (positive definite — not just symmetric).
• The decomposition of $A$-norm minimisation along an $A$-conjugate basis into independent 1D minima needs $\langle \cdot, \cdot \rangle _ A$ to be an inner product (i.e. $A \succ 0$, not merely $A \succeq 0$); otherwise an $A$-conjugate “basis” need not be linearly independent.
• Many CG identities expand $A p _ k$ via $r _ k - r _ {k+1} = \alpha _ k A p _ k$. Dividing by $\alpha _ k = \langle r _ k, r _ k \rangle / \langle p _ k, A p _ k \rangle$ implicitly assumes $\langle p _ k, A p _ k \rangle > 0$, again requiring $A \succ 0$.

Symmetry alone is not enough — for symmetric indefinite $A$, MINRES is used instead (it minimises the 2-norm of the residual, not the $A$-norm of the error).

• Rewrite the $A$-norm error minimisation as a polynomial-approximation problem: every $x \in x _ 0 + \mathcal K _ k(A, r _ 0)$ has the form $x _ 0 + p _ {k-1}(A) r _ 0$, so $x - x _ \ast = p _ {k-1}(A) r _ 0 - e _ 0 = -q(A) e _ 0$ with $q(0) = 1$ (using $r _ 0 = A e _ 0$). Hence $\|e _ k\| _ A / \|e _ 0\| _ A = \min _ {q \in \mathcal P _ k, q(0) = 1} \|q(A) e _ 0\| _ A / \|e _ 0\| _ A$.
• Diagonalise $A = V \Lambda V^\top$, then $\|q(A) e _ 0\| _ A^2 = \sum _ i \lambda _ i q(\lambda _ i)^2 (V^\top e _ 0) _ i^2 \le (\max _ j q(\lambda _ j))^2 \|e _ 0\| _ A^2$.
• Bound by relaxing $\max _ j q(\lambda _ j)$ to $\max _ {x \in [\lambda _ {\min}, \lambda _ {\max}]} q(x)$.
• Plug in the shifted-and-scaled Chebyshev polynomial (∆cg-chebyshev-bound-on-eigenvalue-interval), giving the $\sqrt\kappa$ rate.