Classic Linear Solvers - A Review

Summary

This note casts the abstract iterative skeleton LinSolve_Project (Algorithm 17.2) in the language of classical numerical linear algebra: projection methods, orthogonal / Galerkin methods, conjugate-direction methods, and Krylov subspace methods. The main structural results (Theorem 18.2, Lemma 18.3, Theorem 18.4, Corollary 18.5) pin down exactly which properties of the abstract inference rule make the skeleton equal to Conjugate Gradients. It closes with preconditioning and the observation that, unlike quadrature, no single Gaussian prior is forced by these classical constraints.

Overview

The point of Ch. 18 is to connect the generic solver of The Linear Algebra Problem and Evaluation Strategies (Algorithm 17.2), whose inference rule in line 12 is left abstract, to established families of classical solvers. Doing so lets us later (Corollary 18.5) identify a requirement that any probabilistic solver must satisfy to be equivalent to CG. The taxonomy is nested:

$$\text{projection}\supset \text{orthogonal (Galerkin)}\supset \text{conjugate-direction}\supset \text{Krylov}\ni \text{CG}.$$

Main Content

Projection methods

Independent of the inference rule, Algorithm 17.2 selects its estimate as $x_i\in x_0+{\mathbb{K}}_i$ with the expanding subspace ${\mathbb{K}}_i=\mathrm{span}\{d_1,\dots ,d_i\}$ spanned by the search directions.

Projection method ^def-projection

An algorithm choosing estimates $x_i\in x_0+{\mathbb{K}}_i$ for some sequence of subspaces ${\mathbb{K}}_i$ is a projection method. Because the step size ${\alpha }_i$ (Eq. 17.3) minimises $f$ along $d_i$, the resulting gradient satisfies a constraint: there is a space ${\mathbb{L}}_i$ orthogonal to the gradient,

$$\nabla f(x_i)=r_i=A{x}_i-b\perp {\mathbb{L}}_i.\text{\hspace{1em}(18.1)}$$

If $\mathbb{K}=\mathbb{L}$ the method is an orthogonal projection method and (18.1) is the Galerkin condition. If $\mathbb{K}\ne \mathbb{L}$ the method is oblique (Petrov–Galerkin).

Symbols: ${\mathbb{K}}_i$ = space in which the estimate lives; ${\mathbb{L}}_i$ = space orthogonal to the residual; $r_i$ = residual/gradient at step $i$.

Conjugate directions

Conjugate directions (A-conjugacy) ^def-conjugate

Given symmetric $A\in {\mathbb{R}}^{N\times N}$, two vectors $v,w\in {\mathbb{R}}^N$ are $A$-conjugate if $v^{\top }Aw=0$. A solver for $Ax=b$ with symmetric $A$ is a conjugate directions method if it chooses projections $d_i$ that are pairwise $A$-conjugate.

Conjugate-direction methods enjoy:

• Linear consistency: they converge to the exact $x$ in at most $N$ steps (Nocedal & Wright, Thm. 5.1).
• The gradient $r_k$ after $k$ steps is orthogonal to all preceding directions, $r_k^{\top }{d}_i=0$ for $i<k$ (Nocedal & Wright, Thm. 5.2) — so they are orthogonal projection methods.

Theorem 18.2 (symmetric estimator ⇒ conjugate directions) ^thm-182

If $A$ is symmetric, and the inference rule in line 12 of Algorithm 17.2 produces a symmetric estimator $H_i$, then Algorithm 17.2 is a conjugate directions method.

(Proof by induction on the directions, using symmetry of $H_i$ in the key algebraic step; full proof in Ch. 22, reproduced in ^proof-182.)

This is the load-bearing bridge to the probabilistic story: symmetry of the inverse estimate is exactly what buys $A$-conjugacy of the search directions.

Krylov subspace methods

Krylov sequence and Krylov methods ^def-krylov

The Krylov sequence of order $i$ generated by $r_0$ and $A$ is

$${\mathbb{K}}_i={\mathcal{K}}^i(r_0,A)=\mathrm{span}\{r_0,A{r}_0,A^2{r}_0,\dots ,A^{i-1}{r}_0\}.\text{\hspace{1em}(18.2)}$$

A projection method with ${\mathbb{K}}_i$ equal to the Krylov sequence is a Krylov subspace method. Examples: CG, GMRES, BiCG, CGS, QMR. CG is closely tied to the Lanczos process for eigenvalues of symmetric matrices; GMRES to the more general Arnoldi process.

Because line 12 of Algorithm 17.2 is generic, the skeleton is not a Krylov method in general — §18.4 states what is needed to make it one.

Lemma 18.3 (Krylov characterisation) ^thm-183

Assumption (18.3) — namely that the search directions satisfy

$$d_i=-H_{i-1}{r}_{i-1}\in \mathrm{span}\{s_1,\dots ,s_{i-1},y_1,\dots ,y_{i-1},r_{i-1}\}\text{\hspace{1em}(18.3)}$$

— is equivalent to Algorithm 17.2 being a Krylov subspace method, i.e. to

$$s_i\in \mathrm{span}\{r_0,A{r}_0,A^2{r}_0,\dots ,A^{i-1}{r}_0\}.\text{\hspace{1em}(18.4)}$$

(Proof by a nearly trivial induction; Ch. 22.2. Here $s_i={\alpha }_i{d}_i$, $y_i={\alpha }_i{z}_i=A{s}_i$, and $r_{i-1}=r_0+{\sum }_{j<i}{y}_j$.)

Conjugate Gradients as the intersection

Theorem 18.4 (CG conditions) ^thm-184

If $A$ is spd, $H_i$ is symmetric for all $i\ge 0$, Assumption (18.3) holds, and Algorithm 17.2 does not terminate before step $k<N$, then the residuals are mutually orthogonal,

$$r_i\perp r_j\forall 0\le i\ne j\le k,$$

and there exist scalars ${\gamma }_i\in \mathbb{R}\setminus 0$ for $i<k$ such that line 12 can be written

$$d_i=-H_{i-1}{r}_{i-1}={\gamma }_i\left(-r_{i-1}+\frac{{\beta }_i}{{\gamma }_{i-1}}{d}_{i-1}\right),{\beta }_i:=\frac{r_{i-1}^{\top }{r}_{i-1}}{r_{i-2}^{\top }{r}_{i-2}}.\text{\hspace{1em}(18.5)}$$

Comparing with CG (Algorithm 16.1), the directions are identical up to rescaling by ${\gamma }_i$: $d_i^{\mathrm{CG}}={\gamma }_i{d}_i^{\text{Probabilistic}}$.

(Proof in Ch. 22.3; the full argument is developed in Conjugate Gradients as Probabilistic Inference.)

Corollary 18.5 (equivalence to CG) ^thm-185

Under the assumptions of Theorem 18.4 on the estimators $H_i$, Algorithm 17.2 is equivalent to the method of conjugate gradients: it produces the exact same sequence $\{x_0,x_1,\dots ,x_N\}$ of estimates as CG initialised at $x_0=H_{0}b$.

(The scaling ${\gamma }_i$ is cancelled by the step size ${\alpha }_i$ in line 7, so both algorithms produce identical estimates.)

CG is thus characterisable as the conjugate-direction method that is also a Krylov subspace method. Fig. 18.1 gives the geometric intuition: gradients sampled by CG are drawn (in effect) from a spherical Gaussian around the extremum and are dominated by the eigenvectors of the largest eigenvalues — the Lanczos process’s low-rank capture of the dominant structure that underlies conjugate gradients.

Preconditioning

Krylov subspaces are invariant under scaling (${\mathcal{K}}^m(\sigma r,\tau A)={\mathcal{K}}^m(r,A)$), translation (${\mathcal{K}}^m(r,A+\sigma I)={\mathcal{K}}^m(r,A)$), and orthonormal change of basis (${\mathcal{K}}^m(Ur,UA{U}^{\top })=U{\mathcal{K}}^m(r,A)$) (Eq. 18.6).

Preconditioned system ^def-precond

Given a non-singular preconditioner $C$, solve the transformed problem

$$\tilde{A}\tilde{x}=\tilde{b},\tilde{A}:=C^{-\top }A{C}^{-1},\tilde{x}:=Cx,\tilde{b}:=C^{-\top }b.\text{\hspace{1em}(18.7)}$$

A good $C$ makes CG converge much faster. Preconditioning is externalised from the solver by adding lines that solve $Ky=r_{k+1}$ with $K=C{C}^{\top }$ (Algorithm 18.1, pCG). The ideal but impractical choice is $K=A$ (one-step solution). Practical choices are incomplete/approximate decompositions of $A$.

In the probabilistic language: a measure $p_{\mathrm{CG}}(A)$ consistent with CG, altered so directions map to the transformed space,

$$d_i\in \mathrm{span}\{s_1,\dots ,s_{i-1},y_1,\dots ,y_{i-1},K^{-1}{r}_{i-1}\},\text{\hspace{1em}(18.8)}$$

becomes consistent with preconditioned CG. This modification will later (Corollary 19.17) correspond to a change of the prior of the probabilistic solver, with the ideal $K=A$ associated with a particularly “natural” prior.

Examples

Why linear algebra has no forced Gaussian prior (contrast with integration)

None of the restrictions above (projection, conjugacy, Krylov) implies a concrete — let alone Gaussian — form for $p(A)$. This differs sharply from integration, where each quadrature rule maps (up to one scalar) onto a specific Gaussian process prior. The reason: an $N\times N$ invertible matrix has a finite $N^2$ degrees of freedom; at each step the solver identifies $N$ of them perfectly while learning nothing about the other $N(N-1)$. Linear algebra is largely sophisticated book-keeping rather than inference of an intractable object — so uncertainty is a “corner case” that classical solvers legitimately ignore.

When uncertainty does matter

Uncertainty in linear solvers becomes prominent when: (i) the object is an infinite-dimensional operator (PDEs); (ii) one needs the matrix inverse itself more than the solution (e.g. Laplace approximations / inverse Hessians of deep nets, which have few prominent eigenvalues and are stopped very early); or (iii) projection observations are corrupted by noise (data subsampling — sometimes re-phrasable as a Schur-complement projection).

Connections

• Theorem 18.2 (symmetry ⇒ conjugacy) is the pivot reused in Probabilistic Linear Solvers - Algorithmic Scaffold (Theorem 19.15) once the estimator is Gaussian.
• The Krylov/Lanczos picture reappears in Uncertainty Calibration for Linear Solvers via Rayleigh regression.
• Preconditioning-as-prior (Corollary 19.17) is developed in Probabilistic Linear Solvers - Algorithmic Scaffold.

See Also

• The Linear Algebra Problem and Evaluation Strategies — defines Algorithm 17.2 that these results characterise.
• Conjugate Gradients as Probabilistic Inference — states Theorem 18.4 / Corollary 18.5 with full proof sketch.
• Probabilistic Linear Solvers - Algorithmic Scaffold — the Gaussian models satisfying these classical requirements.
• Gauss-Markov Processes and SDEs — Krylov/Lanczos parallels the Markovian state-space compression seen elsewhere.