Preconditioning

Intuition Beginner

The Krylov solvers of the last few units — conjugate gradient for symmetric problems, GMRES for general ones — are fast only when the matrix is well behaved. A well-behaved matrix is one whose stretching is roughly the same in every direction, so the bowl you are minimising is nearly round. When the matrix stretches space far more in one direction than another, the bowl is a long thin canyon, and even a clever solver needs many steps to crawl along it.

Preconditioning fixes the matrix before the solver ever runs. The idea is to find a second matrix, call it $M$, that is close to the original $A$ but much cheaper to undo. You then hand the solver the reweighted problem instead of the raw one. If $M$ really is close to $A$, the reweighted matrix is close to doing nothing at all, so the solver sees an almost-round bowl and finishes in a handful of steps.

There is a balance to strike. A preconditioner that exactly equals $A$ would make the reweighted problem perfectly round, but undoing it would be as hard as solving the original system, so nothing is gained. A preconditioner that is too crude is cheap to undo but barely rounds the bowl, so the step count stays high. The art is to pick an $M$ that is cheap to apply and still close enough to $A$ to round the problem out.

This single trick is what turns the Krylov methods from textbook curiosities into the workhorses behind weather models, structural engineering, and circuit simulation. The solver supplies the engine; the preconditioner supplies the smooth road.

Visual Beginner

The picture is the same long thin canyon seen twice: once as the raw problem the solver dreads, and once after preconditioning has reshaped it into a round bowl the solver finishes quickly.

Read the table top to bottom. The left column is the raw system, where the matrix stretches one direction far more than another and the solver takes many steps. The right column is the preconditioned system, where the reshaped matrix stretches every direction about equally and the same solver finishes fast.

feature	raw system $Ax=b$	preconditioned $M^{-1}Ax=M^{-1}b$
shape of the bowl	long thin canyon	nearly round
spread of stretching	very uneven	nearly equal
solver steps needed	many	few
cost per step	one multiply by $A$	one multiply by $A$ plus one undo of $M$

The takeaway: a Krylov solver is fast on round problems and slow on stretched ones. Preconditioning multiplies in a cheap approximate undo of the matrix, reshaping the stretched problem into a round one before the solver starts. The extra work is one undo of $M$ per step; the payoff is far fewer steps, and on hard problems that trade is hugely in your favour.

Worked example Beginner

Take a diagonal stand-in for a stretched matrix, $$ A = \begin{pmatrix} 100 & 0 \ 0 & 1 \end{pmatrix}, \qquad b = \begin{pmatrix} 100 \ 1 \end{pmatrix}, $$ whose exact answer is $x=(1,1)$. The matrix stretches the first direction $100$ times more than the second, so its condition number — the ratio of the largest stretch to the smallest — is $100/1=100$. A Krylov solver is slow when this ratio is large.

Choose the simplest preconditioner: the diagonal of $A$ itself, $M=\left(\begin{array}{cc}100& 0\\ 0& 1\end{array}\right)$, which here equals $A$. Undoing $M$ is one division per entry, very cheap. The reshaped matrix is $$ M^{-1}A = \begin{pmatrix} 1/100 & 0 \ 0 & 1 \end{pmatrix}\begin{pmatrix} 100 & 0 \ 0 & 1 \end{pmatrix} = \begin{pmatrix} 1 & 0 \ 0 & 1 \end{pmatrix}. $$ The reshaped matrix is the identity. Its condition number is $1/1=1$, the best possible, down from $100$. The reshaped right-hand side is $M^{-1}b=(1,1)$, and the reshaped problem $Ix=(1,1)$ is solved in one step: $x=(1,1)$, the exact answer.

What this tells us: dividing each row by its own diagonal entry — the Jacobi preconditioner — collapsed a condition number of $100$ down to $1$, so the solver finished immediately instead of crawling. The undo cost just two divisions. On a real matrix the diagonal is only an approximation, so the reshaped condition number will not fall all the way to $1$, but the same mechanism shrinks it and speeds the solver up.

Check your understanding Beginner

Formal definition Intermediate+

Let $A\in {\mathbb{F}}^{n\times n}$ ($\mathbb{F}\in \{\mathbb{R},\mathbb{C}\}$) be nonsingular and $b\in {\mathbb{F}}^n$. A preconditioner is a nonsingular matrix $M$, chosen so that (i) systems $Mz=r$ are cheap to solve relative to $Ax=b$, and (ii) the preconditioned operator is better conditioned or has a more clustered spectrum than $A$. One never forms $M^{-1}$ explicitly; "applying $M^{-1}$" means solving one system with $M$.

Left, right, and split preconditioning. There are three algebraically equivalent placements of $M$ ^{[Trefethen, L. N. & Bau, D. — Numerical Linear Algebra (SIAM, 1997)]}: $$ \underbrace{M^{-1}A,x = M^{-1}b}{\text{left}}, \qquad \underbrace{A M^{-1}u = b,\ \ x = M^{-1}u}{\text{right}}, \qquad \underbrace{L^{-1} A L^{-},w = L^{-1}b,\ \ x = L^{-}w}_{\text{split, } M = L L^}. $$ Left preconditioning changes the residual the Krylov method sees from $b-Ax$ to $M^{-1}(b-Ax)$; right preconditioning leaves the true residual unchanged and is preferred when a residual-based stopping test must track $\Vert b-Ax\Vert$. Split preconditioning, available when $M = LL^$isaHermitianpositive-definitefactorisation,applieshalfof$M^{-1}$oneachsidesothatthepreconditionedoperator$L^{-1}AL^{-*}$inheritsthesymmetryanddefinitenessof$A$ — the requirement that makes conjugate gradient applicable.

The spectrum-clustering goal. The convergence theory of 43.07.04 and 43.07.03 shows that the Krylov error after $k$ steps is governed by a min-max polynomial on the spectrum of the operator the method runs on: for CG the energy-norm error obeys $\Vert e_k{\Vert }_A/\Vert e_0{\Vert }_A\le {\min }_{p(0)=1,\mathrm{deg}p\le k}{\max }_{\lambda \in \sigma }|p(\lambda )|$. A degree-$k$ polynomial pinned to $p(0)=1$ is small on the spectrum exactly when the eigenvalues lie in a few tight clusters, or in a short interval bounded away from $0$. Preconditioning is the deliberate construction of an $M$ for which $\sigma (M^{-1}A)$ is so distributed — ideally clustered near $1$, since $M^{-1}A\approx I$ when $M\approx A$. The relevant figure of merit is the condition number $\kappa (M^{-1}A)$ for symmetric problems (replacing $\kappa (A)$ in the Chebyshev bound) and the spectral/field-of-values distribution of $M^{-1}A$ for nonsymmetric ones.

Preconditioned conjugate gradient (PCG). When $A\succ 0$ and $M\succ 0$, $M^{-1}A$ is generally not symmetric, yet CG can still be applied — by running it in the $M$-inner product $⟨u,v{⟩}_M=u^{\ast }Mv$, in which $M^{-1}A$ is self-adjoint. The resulting algorithm applies $M^{-1}$ exactly once per iteration, to the residual, producing the preconditioned residual $z_k=M^{-1}{r}_k$: $$ \alpha_k = \frac{r_k^* z_k}{p_k^* A p_k},\quad x_{k+1} = x_k + \alpha_k p_k,\quad r_{k+1} = r_k - \alpha_k A p_k,\quad z_{k+1} = M^{-1}r_{k+1},\quad \beta_k = \frac{r_{k+1}^* z_{k+1}}{r_k^* z_k},\quad p_{k+1} = z_{k+1} + \beta_k p_k. $$ This is the CG recurrence of 43.07.04 with every Euclidean residual inner product $r^{\ast }r$ replaced by $r^{\ast }z=r^{\ast }{M}^{-1}r$. Preconditioned GMRES is simpler: it runs the GMRES of 43.07.03 on $M^{-1}A$ (left) or $A{M}^{-1}$ (right) directly, since GMRES needs no symmetry.

Classical and modern preconditioners. The Jacobi preconditioner is $M=D=\mathrm{diag}(A)$, undone by one division per entry. The Gauss-Seidel and symmetric SOR (SSOR) preconditioners are built from the stationary splittings $A=D-L-U$ of 43.07.01: SSOR uses $M=(D-\omega L)D^{-1}(D-\omega U)/(2-\omega )$ up to scaling, a Hermitian positive-definite factorisation suitable for PCG. Incomplete LU (ILU) and incomplete Cholesky (IC) compute approximate factors $A\approx \bar{L}\bar{U}$ by performing Gaussian elimination but discarding fill-in outside a prescribed sparsity pattern, so $M=\bar{L}\bar{U}$ keeps $A$'s sparsity. Multigrid, domain decomposition, and sparse approximate inverse preconditioners are treated at high level in the Advanced results.

Counterexamples to common slips

• The preconditioned operator is not symmetric in the Euclidean sense. For $A\succ 0$ and $M\succ 0$ the matrix $M^{-1}A$ is self-adjoint in the $M$-inner product, not in the standard one. PCG works because it tacitly runs in $⟨\cdot ,\cdot {⟩}_M$; applying ordinary CG to the explicitly formed $M^{-1}A$ would be inconsistent.
• Left and right preconditioning are not interchangeable for stopping. They generate the same iterates in exact arithmetic up to the change of variable, but the residual a left-preconditioned method monitors is $M^{-1}(b-Ax)$, not $b-Ax$. A stopping test calibrated to the true residual must use right preconditioning or unwind the preconditioned residual.
• $M$ must be definite to precondition CG. A poor indefinite $M$ can destroy the positive-definiteness that CG requires, causing the denominator $p_k^{\ast }A{p}_k$ or the inner product $r_k^{\ast }{z}_k$ to lose its sign. SSOR and incomplete Cholesky are chosen precisely to keep $M\succ 0$.
• A clustered spectrum, not a small norm, is the goal. Making $\Vert M^{-1}A-I\Vert$ small is sufficient but not necessary; what the polynomial bound rewards is a few tight eigenvalue clusters. A preconditioner that leaves a handful of outliers but clusters the bulk can still give fast convergence, since the residual polynomial spends a few roots on the outliers.

Key theorem with proof Intermediate+

The signature result is that preconditioned conjugate gradient is exactly ordinary conjugate gradient applied to the symmetrically split operator $\hat{A}=L^{-1}A{L}^{-\ast }$, so the whole convergence theory of 43.07.04 transfers verbatim with $\kappa (A)$ replaced by $\kappa (M^{-1}A)$; and that the split, left, and $M$-inner-product forms produce the same iterates.

Theorem (PCG equivalence and the preconditioned Chebyshev bound). Let $A\succ 0$ and $M\succ 0$ with a Hermitian factorisation $M = LL^$.Set$\hat A = L^{-1}AL^{-}$,$\hat b = L^{-1}b$,and$\hat x = L^x$.Then$\hat A \succ 0$issimilarto$M^{-1}A$(so$\sigma(\hat A) = \sigma(M^{-1}A)$),andtheiterates$\hat x_k$producedbyordinaryCGon$\hat A\hat x = \hat b$correspond,under$x_k = L^{-}\hat x_k$, exactly to the PCG iterates of the Formal definition. Consequently $$ |x_k - x_\star|A \le 2\left(\frac{\sqrt{\hat\kappa} - 1}{\sqrt{\hat\kappa} + 1}\right)^{k}|x_0 - x\star|A, \qquad \hat\kappa = \kappa(M^{-1}A) = \frac{\lambda{\max}(M^{-1}A)}{\lambda_{\min}(M^{-1}A)}. $$ ^{[Saad, Y. — Iterative Methods for Sparse Linear Systems (2nd ed.)]}

Proof. First, $\hat{A}=L^{-1}A{L}^{-\ast }$ is Hermitian since $A=A^{\ast }$, and for $w\ne 0$, $w^{\ast }\hat{A}w=(L^{-\ast }w{)}^{\ast }A(L^{-\ast }w)>0$ because $A\succ 0$ and $L^{-\ast }$ is nonsingular, so $\hat{A}\succ 0$. The similarity $L^{\ast }(M^{-1}A)L^{-\ast }=L^{\ast }(L{L}^{\ast }{)}^{-1}A{L}^{-\ast }=L^{\ast }{L}^{-\ast }{L}^{-1}A{L}^{-\ast }=L^{-1}A{L}^{-\ast }=\hat{A}$ shows $\hat{A}$ and $M^{-1}A$ are similar, hence share a spectrum, which is real and positive because $\hat{A}\succ 0$.

Now run ordinary CG of 43.07.04 on $\hat{A}\hat{x}=\hat{b}$ from ${\hat{x}}_0=L^{\ast }{x}_0$. Its quantities are ${\hat{r}}_k=\hat{b}-\hat{A}{\hat{x}}_k$, ${\hat{p}}_k$, with the recurrences ${\hat{\alpha }}_k=({\hat{r}}_k^{\ast }{\hat{r}}_k)/({\hat{p}}_k^{\ast }\hat{A}{\hat{p}}_k)$, ${\hat{x}}_{k+1}={\hat{x}}_k+{\hat{\alpha }}_k{\hat{p}}_k$, ${\hat{r}}_{k+1}={\hat{r}}_k-{\hat{\alpha }}_k\hat{A}{\hat{p}}_k$, ${\hat{\beta }}_k=({\hat{r}}_{k+1}^{\ast }{\hat{r}}_{k+1})/({\hat{r}}_k^{\ast }{\hat{r}}_k)$, ${\hat{p}}_{k+1}={\hat{r}}_{k+1}+{\hat{\beta }}_k{\hat{p}}_k$. Introduce the untransformed variables $x_k=L^{-\ast }{\hat{x}}_k$, $r_k=L{\hat{r}}_k$, $p_k=L^{-\ast }{\hat{p}}_k$, and $z_k=M^{-1}{r}_k$. Then ${\hat{r}}_k=L^{-1}{r}_k$, so ${\hat{r}}_k^{\ast }{\hat{r}}_k=r_k^{\ast }{L}^{-\ast }{L}^{-1}{r}_k=r_k^{\ast }(L{L}^{\ast }{)}^{-1}{r}_k=r_k^{\ast }{M}^{-1}{r}_k=r_k^{\ast }{z}_k$. Likewise ${\hat{p}}_k^{\ast }\hat{A}{\hat{p}}_k={\hat{p}}_k^{\ast }{L}^{-1}A{L}^{-\ast }{\hat{p}}_k=p_k^{\ast }A{p}_k$, giving ${\hat{\alpha }}_k=(r_k^{\ast }{z}_k)/(p_k^{\ast }A{p}_k)={\alpha }_k$. The iterate update ${\hat{x}}_{k+1}={\hat{x}}_k+{\hat{\alpha }}_k{\hat{p}}_k$ becomes $x_{k+1}=x_k+{\alpha }_k{p}_k$ on multiplying by $L^{-\ast }$; the residual update ${\hat{r}}_{k+1}={\hat{r}}_k-{\hat{\alpha }}_k\hat{A}{\hat{p}}_k$ becomes, on multiplying by $L$, $r_{k+1}=r_k-{\alpha }_{k}A{L}^{-\ast }{\hat{p}}_k=r_k-{\alpha }_{k}A{p}_k$. The coefficient ${\hat{\beta }}_k=(r_{k+1}^{\ast }{z}_{k+1})/(r_k^{\ast }{z}_k)={\beta }_k$, and ${\hat{p}}_{k+1}={\hat{r}}_{k+1}+{\hat{\beta }}_k{\hat{p}}_k$ becomes $p_{k+1}=L{\hat{r}}_{k+1}\cdot$ — precisely $p_{k+1}=z_{k+1}+{\beta }_k{p}_k$ after multiplying by $L^{-\ast }$ and using $L^{-\ast }{\hat{r}}_{k+1}=L^{-\ast }{L}^{-1}{r}_{k+1}=M^{-1}{r}_{k+1}=z_{k+1}$. These are exactly the PCG recurrences, so the iterates coincide.

Finally, ordinary CG on $\hat{A}$ obeys the Chebyshev bound of 43.07.04 in the $\hat{A}$-norm: $\Vert {\hat{x}}_k-{\hat{x}}_{\star }{\Vert }_{\hat{A}}\le 2((\sqrt{\hat{\kappa }}-1)/(\sqrt{\hat{\kappa }}+1){)}^k\Vert {\hat{x}}_0-{\hat{x}}_{\star }{\Vert }_{\hat{A}}$ with $\hat{\kappa }=\kappa (\hat{A})=\kappa (M^{-1}A)$. The $\hat{A}$-norm pulls back to the $A$-norm: $\Vert \hat{x}-{\hat{x}}_{\star }{\Vert }_{\hat{A}}^2=(\hat{x}-{\hat{x}}_{\star }{)}^{\ast }\hat{A}(\hat{x}-{\hat{x}}_{\star })=(x-x_{\star }{)}^{\ast }L{L}^{-1}A{L}^{-\ast }{L}^{\ast }(x-x_{\star })=(x-x_{\star }{)}^{\ast }A(x-x_{\star })=\Vert x-x_{\star }{\Vert }_A^2$, using $\hat{x}-{\hat{x}}_{\star }=L^{\ast }(x-x_{\star })$. Substituting gives the stated bound. $\square$

Bridge. This equivalence is the foundational reason preconditioning costs nothing in theory and everything in practice: it builds toward the convergence and design theory of the Advanced results by showing that the entire CG apparatus of 43.07.04 — the energy-norm optimality, the conjugacy, the residual-polynomial reformulation — survives intact under preconditioning, with $\kappa (A)$ simply overwritten by $\kappa (M^{-1}A)$. This is exactly the spectrum-reshaping principle made quantitative: the Chebyshev bound that turned $\kappa$ into $\sqrt{\kappa }$ now runs on the preconditioned condition number, so the central insight is that a preconditioner's only job is to shrink or cluster $\sigma (M^{-1}A)$, and the $\sqrt{\hat{\kappa }}$ count generalises the $\sqrt{\kappa }$ count of unpreconditioned CG. The split form $\hat{A}=L^{-1}A{L}^{-\ast }$ is dual to the $M$-inner-product form: one symmetrises by changing the operator, the other by changing the geometry, and both yield identical iterates. The same construction appears again in 43.07.03, where preconditioned GMRES runs the Arnoldi-based minimisation on $M^{-1}A$ or $A{M}^{-1}$ and the field of values of the preconditioned operator replaces that of $A$. Putting these together, the bridge is that preconditioning is a similarity transform chosen to improve a spectrum, executed implicitly so that the iterate sequence is exactly that of the unpreconditioned method on the better-conditioned operator.

Exercises Intermediate+

Advanced results Master

Theorem 1 (spectral equivalence and mesh-independent convergence). Let $\{A_h\}$ be a family of SPD matrices arising from discretising an elliptic operator on a mesh of size $h$, with $\kappa (A_h)=O(h^{-2})\to \infty$ as $h\to 0$. A preconditioner family $\{M_h\}$ is spectrally equivalent *to $\{A_h\}$ if there are constants $0<c_1\le c_2$, independent of $h$, with $c_1{u}^{\ast }{M}_{h}u\le u^{\ast }{A}_{h}u\le c_2{u}^{\ast }{M}_{h}u$ for all $u$. Then $\kappa (M_h^{-1}{A}_h)\le c_2/c_1$ uniformly in $h$, and PCG converges to a fixed tolerance in a number of iterations bounded independently of $h$.* The bound is the condition-number estimate of Exercise 4 made uniform; a preconditioner achieving it is called optimal, and multigrid and certain domain-decomposition preconditioners achieve it for the model elliptic problems ^{[Greenbaum, A. — Iterative Methods for Solving Linear Systems (SIAM, 1997)]}.

Theorem 2 (existence of incomplete factorisations for M-matrices). Let $A$ be a symmetric M-matrix (symmetric positive-definite with nonpositive off-diagonal entries), and fix a sparsity pattern $P\supseteq \{(i,i)\}$. The incomplete Cholesky factorisation $A = \bar L\bar L^ - R$,computedbyCholeskyeliminationwithallfill-inoutside$P$settozero,existswith$\bar L$havingrealpositivediagonal,andthepreconditioner$M = \bar L\bar L^$ is symmetric positive-definite. The Meijerink-van der Vorst theorem guarantees every pivot stays positive during incomplete elimination because dropping a nonpositive fill entry only increases the remaining diagonal of an M-matrix; the same argument gives ILU($0$) existence for nonsymmetric M-matrices. This is the theoretical license for the ICCG method that made preconditioned CG the standard sparse SPD solver ^{[Meijerink, J. A. & van der Vorst, H. A. — An Iterative Solution Method for Linear Systems of Which the Coefficient Matrix is a Symmetric M-Matrix]}.

Theorem 3 (incomplete-factorisation hierarchy and fill control). ILU($0$) retains exactly the sparsity of $A$, discarding all fill; ILU($p$) retains fill of level at most $p$, where the level of an entry records how many elimination steps generated it; ILUT($\tau ,\ell$) drops entries below a relative threshold $\tau$ and keeps at most $\ell$ per row. Increasing $p$ or lowering $\tau$ makes $M=\bar{L}\bar{U}$ a more accurate approximation of $A$, clustering $\sigma (M^{-1}A)$ more tightly near $1$ and reducing the iteration count, at the cost of denser factors and more expensive applies. The hierarchy makes explicit the central tradeoff: the preconditioner's quality (how clustered it renders the spectrum) is bought with its construction and per-application cost, and the optimum minimises total time, not iteration count ^{[Saad, Y. — Iterative Methods for Sparse Linear Systems (2nd ed.)]}.

Theorem 4 (multigrid, domain decomposition, and sparse approximate inverse as preconditioners). Three modern families supply preconditioners whose applies are matrix-free or embarrassingly parallel. Multigrid uses a hierarchy of grids with smoothing (a stationary iteration of 43.07.01) on each level and coarse-grid correction; one V-cycle is a linear operator $M^{-1}$ that, for the model Poisson problem, is spectrally equivalent to $A$, giving an optimal $O(n)$ preconditioner. Domain decomposition (additive Schwarz) splits the domain into overlapping subdomains, solves local problems, and sums the corrections: $M^{-1}={\sum }_i{R}_i^{\ast }{A}_i^{-1}{R}_i$, with a coarse space added to bound the condition number independently of the number of subdomains. Sparse approximate inverse (SPAI) computes a sparse $M\approx A^{-1}$ directly by minimising $\Vert I-MA{\Vert }_F$ column by column, a set of independent small least-squares problems of 43.04.01, yielding a preconditioner whose apply is a single sparse matrix-vector product and so is fully parallel. Each trades a different resource — grid hierarchy, subdomain solves, or up-front least-squares work — for a clustered preconditioned spectrum ^{[Saad, Y. — Iterative Methods for Sparse Linear Systems (2nd ed.)]}.

Synthesis. Preconditioning is one operation — replace $A$ by a spectrally improved similar operator $M^{-1}A$ that is implicit, never formed — viewed under one invariant, the spectrum of the operator the Krylov method actually runs on, and the foundational reason it works is that the residual-polynomial convergence theory of 43.07.04 and 43.07.03 depends on $A$ only through that spectrum, so reshaping it is the entire lever. This is exactly the Chebyshev and min-max polynomial machinery of those units read backward: where CG turned $\kappa$ into $\sqrt{\kappa }$ and rewarded clustered spectra, preconditioning manufactures the small $\kappa (M^{-1}A)$ or the tight clusters that the polynomial bound rewards, and the central insight is that the $\sqrt{\kappa (M^{-1}A)}$ iteration count of preconditioned CG generalises the $\sqrt{\kappa (A)}$ count of plain CG with $A$ overwritten by the preconditioned operator. The split preconditioner $L^{-1}A{L}^{-\ast }$ is dual to the $M$-inner-product formulation — change the operator or change the geometry, identical iterates — and this duality is the bridge from the stationary splittings of 43.07.01, which were too weak to solve but are exactly strong enough to precondition: Jacobi, Gauss-Seidel, and SSOR all return here as $M$, their spectral radius creeping to one as solvers being precisely the defect a preconditioner need not cure, only approximate.

Putting these together, the chapter closes a circle. The stationary methods of 43.07.01 supplied the splittings; the Krylov methods of 43.07.02, 43.07.03, and 43.07.04 supplied the optimal polynomial acceleration whose rate is set by the spectrum; and preconditioning fuses the two, using a cheap splitting or incomplete factorisation or multilevel operator to reshape the spectrum so the optimal acceleration needs only a handful of steps. The condition number $\kappa$ of 43.01.02 is the quantity every member of the family manipulates, and the design space — Jacobi to incomplete Cholesky to multigrid — is a ladder of cost-versus-clustering tradeoffs, with the optimal preconditioner of Theorem 1 the one whose $\kappa (M^{-1}A)$ is bounded independently of the problem size while its apply stays cheap, converting the divergent $O(h^{-1})$ iteration count of the raw discretised operator into a constant.

Full proof set Master

Proposition 1 (split symmetrisation preserves definiteness and spectrum). For $A\succ 0$ and $M = LL^ \succ 0$,thesplitoperator$\hat A = L^{-1}AL^{-}$isHermitianpositive-definiteandsimilarto$M^{-1}A$,so$\sigma(M^{-1}A) = \sigma(\hat A) \subset (0,\infty)$.

Proof. ${\hat{A}}^{\ast }=(L^{-1}A{L}^{-\ast }{)}^{\ast }=L^{-1}{A}^{\ast }{L}^{-\ast }=L^{-1}A{L}^{-\ast }=\hat{A}$, Hermitian. For $w\ne 0$, $w^{\ast }\hat{A}w=(L^{-\ast }w{)}^{\ast }A(L^{-\ast }w)>0$ since $A\succ 0$ and $L^{-\ast }$ nonsingular, so $\hat{A}\succ 0$. The similarity $L^{\ast }(M^{-1}A)L^{-\ast }=L^{\ast }(L{L}^{\ast }{)}^{-1}A{L}^{-\ast }=L^{-1}A{L}^{-\ast }=\hat{A}$ gives $\sigma (M^{-1}A)=\sigma (\hat{A})$, real and positive by $\hat{A}\succ 0$. $\square$

Proposition 2 (PCG equals CG on the split operator). The PCG recurrences of the Formal definition produce iterates $x_k = L^{-}\hat x_k$,where$\hat x_k$aretheordinaryCGiterateson$\hat A\hat x = \hat b$with$\hat b = L^{-1}b$,$\hat x_0 = L^x_0$.

Proof. This is the Key theorem's computation: substituting ${\hat{r}}_k=L^{-1}{r}_k$, ${\hat{p}}_k=L^{\ast }{p}_k$, $z_k=M^{-1}{r}_k$ into the CG recurrences for $\hat{A}$ reproduces the PCG coefficients ${\alpha }_k=(r_k^{\ast }{z}_k)/(p_k^{\ast }A{p}_k)$, ${\beta }_k=(r_{k+1}^{\ast }{z}_{k+1})/(r_k^{\ast }{z}_k)$ and the updates $x_{k+1}=x_k+{\alpha }_k{p}_k$, $r_{k+1}=r_k-{\alpha }_{k}A{p}_k$, $p_{k+1}=z_{k+1}+{\beta }_k{p}_k$, using $r_k^{\ast }{z}_k={\hat{r}}_k^{\ast }{\hat{r}}_k$ and $p_k^{\ast }A{p}_k={\hat{p}}_k^{\ast }\hat{A}{\hat{p}}_k$. The map ${\hat{x}}_k\mapsto L^{-\ast }{\hat{x}}_k$ is the change of variable. $\square$

Proposition 3 (preconditioned Chebyshev bound). Under the hypotheses of Proposition 1, $\Vert x_k-x_{\star }{\Vert }_A\le 2((\sqrt{\hat{\kappa }}-1)/(\sqrt{\hat{\kappa }}+1){)}^k\Vert x_0-x_{\star }{\Vert }_A$ with $\hat{\kappa }=\kappa (M^{-1}A)$.

Proof. CG on $\hat{A}\succ 0$ obeys the Chebyshev bound of 43.07.04 in the $\hat{A}$-norm with $\kappa (\hat{A})=\kappa (M^{-1}A)=\hat{\kappa }$ (Proposition 1). The norm identity $\Vert \hat{x}-{\hat{x}}_{\star }{\Vert }_{\hat{A}}^2=(L^{\ast }(x-x_{\star }){)}^{\ast }\hat{A}(L^{\ast }(x-x_{\star }))=(x-x_{\star }{)}^{\ast }L{L}^{-1}A{L}^{-\ast }{L}^{\ast }(x-x_{\star })=\Vert x-x_{\star }{\Vert }_A^2$ transports the bound from the $\hat{A}$-norm to the $A$-norm via Proposition 2. $\square$

Proposition 4 (spectral-equivalence condition-number bound). *If $c_1{u}^{\ast }Mu\le u^{\ast }Au\le c_2{u}^{\ast }Mu$ for all $u$ with $M\succ 0$, $A\succ 0$, then $\kappa (M^{-1}A)\le c_2/c_1$.*

Proof. The eigenvalues of $M^{-1}A$ are the stationary values of the generalised Rayleigh quotient $R(u)=(u^{\ast }Au)/(u^{\ast }Mu)$ (the pencil $Au=\lambda Mu$). The hypothesis bounds $c_1\le R(u)\le c_2$, so by Courant-Fischer 01.01.14 for the pencil, ${\lambda }_{\min }(M^{-1}A)\ge c_1$ and ${\lambda }_{\max }(M^{-1}A)\le c_2$; with $\sigma \subset (0,\infty )$ (Proposition 1) this gives $\kappa ={\lambda }_{\max }/{\lambda }_{\min }\le c_2/c_1$. $\square$

Proposition 5 (low-rank preconditioner gives finite GMRES termination). If $M^{-1}A=I+R$ with $\mathrm{rank}(R)=r$ and $M^{-1}A$ nonsingular, preconditioned GMRES reaches the exact solution in at most $r+1$ steps in exact arithmetic.

Proof. $R$ has $0$ as an eigenvalue of multiplicity $\ge n-r$, so $I+R$ has at most $r+1$ distinct eigenvalues, all nonzero ($M^{-1}A$ nonsingular). Let ${\nu }_1,\dots ,{\nu }_t$ ($t\le r+1$) be the distinct eigenvalues and $q(z)={\prod }_{j=1}^t(z-{\nu }_j)$, $p(z)=q(z)/q(0)$ — well defined since $0\notin \{{\nu }_j\}$. Then $p(0)=1$, $\mathrm{deg}p=t\le r+1$, and $p$ vanishes on $\sigma (M^{-1}A)$, so $p(M^{-1}A)=0$ on a diagonalisable operator. The residual-polynomial bound of 43.07.03 gives $\Vert r_t\Vert \le \Vert p(M^{-1}A){\hat{r}}_0\Vert =0$, so GMRES terminates by step $t\le r+1$. $\square$

Proposition 6 (incomplete Cholesky exists for symmetric M-matrices). Let $A$ be a symmetric M-matrix and $P$ a sparsity pattern containing the diagonal. The incomplete Cholesky factorisation with fill restricted to $P$ produces $\bar{L}$ with positive diagonal, and $M = \bar L\bar L^ \succ 0$.*

Proof sketch. Incomplete Cholesky performs the recursion $a_{ij}\leftarrow a_{ij}-{\sum }_{k<i}{\bar{l}}_{ik}{\bar{l}}_{jk}$ but sets to zero any update at a position outside $P$. For an M-matrix the off-diagonal entries are $\le 0$ and the Schur complements of the exact factorisation stay M-matrices with positive diagonal (the exact pivots are positive by SPD-ness). Dropping a fill entry at position $(i,j)\notin P$ removes a nonpositive contribution $-{\bar{l}}_{ik}{\bar{l}}_{jk}$ from an off-diagonal and a corresponding nonnegative amount from the diagonal of the remaining Schur complement, so each incomplete pivot is at least as large as the corresponding exact pivot, hence remains positive. The induction (Meijerink-van der Vorst) shows every diagonal of $\bar{L}$ is real and positive, so $\bar{L}$ is nonsingular and $M=\bar{L}{\bar{L}}^{\ast }$ is Hermitian positive-definite, with $u^{\ast }Mu=\Vert {\bar{L}}^{\ast }u{\Vert }_2^2>0$ for $u\ne 0$. $\square$

Connections Master

• The conjugate gradient method 43.07.04 is the solver preconditioning most directly serves: preconditioned CG is, by the Key theorem, ordinary CG run on the symmetrically split operator $L^{-1}A{L}^{-\ast }$ or equivalently in the $M$-inner product, so the energy-norm optimality, the conjugacy of search directions, and the Chebyshev convergence bound all carry over with $\kappa (A)$ replaced by $\kappa (M^{-1}A)$. The clustered-spectrum acceleration that unit identifies — that the count of distinct or well-separated eigenvalues governs the true rate — is exactly what a preconditioner engineers, which is why incomplete-Cholesky-preconditioned CG (ICCG) is the standard solver for the sparse SPD systems of discretised elliptic PDEs.

• GMRES 43.07.03 is the nonsymmetric counterpart: preconditioned GMRES runs the Arnoldi-based minimal-residual minimisation on $M^{-1}A$ (left) or $A{M}^{-1}$ (right), and the diagonalisable convergence bound and field-of-values analysis of that unit apply to the preconditioned operator. The any-curve theorem there warns that for non-normal $M^{-1}A$ the spectrum alone need not predict convergence, so preconditioner design for GMRES targets the field of values or pseudospectrum, not merely the eigenvalues; the flexible variant FGMRES accommodates a preconditioner that changes from step to step, such as an inner Krylov solve.

• Stationary iterative methods: Jacobi, Gauss-Seidel, SOR 43.07.01 supply the classical preconditioners: the splitting $A=M-N$ that was a slow solver on its own — its iteration matrix $M^{-1}N$ having spectral radius creeping to one as the grid refines — becomes the preconditioner $M$ whose only job is to approximate $A$ cheaply, not to converge alone. Jacobi ($M=D$), Gauss-Seidel, and the symmetric SSOR factorisation built from $A=D-L-U$ all reappear here, and the SSOR construction yields the symmetric positive-definite $M$ that preconditioned CG requires, fusing the stationary and Krylov families that the chapter developed separately.

• Conditioning and condition numbers of problems 43.01.02 furnishes the figure of merit the whole subject optimises: preconditioning is the deliberate manipulation of the condition number $\kappa (M^{-1}A)$, and the spectral-equivalence theory measures a preconditioner by whether $\kappa (M^{-1}A)$ is bounded independently of the discretisation parameter. The least-squares machinery of 43.04.01 underlies the sparse-approximate-inverse preconditioner, whose columns solve independent small least-squares problems $\min \Vert I-MA{\Vert }_F$.

Historical & philosophical context Master

The idea of transforming a linear system to improve the convergence of an iteration predates the Krylov methods: David Young's 1950 thesis on successive over-relaxation, and the splitting framework of Richard Varga's 1962 Matrix Iterative Analysis, established that the spectral radius of $M^{-1}N$ controls a stationary iteration, which is the conceptual seed of using $M$ as a preconditioner rather than as a standalone solver. The word preconditioning in its modern sense, and the recognition that it should be paired with conjugate gradients rather than with a stationary iteration, crystallised in the 1970s.

The decisive step was the 1977 paper of Johannes Meijerink and Henk van der Vorst, which introduced the incomplete Cholesky factorisation as a preconditioner, proved its existence and stability for symmetric M-matrices, and demonstrated that incomplete-Cholesky-preconditioned conjugate gradient (ICCG) solved the sparse SPD systems of discretised PDEs far faster than any stationary method ^{[Meijerink, J. A. & van der Vorst, H. A. — An Iterative Solution Method for Linear Systems of Which the Coefficient Matrix is a Symmetric M-Matrix]}. ICCG converted conjugate gradients from a finite direct method that had disappointed in finite-precision arithmetic — the revival narrative of 43.07.04 — into the dominant iterative solver for elliptic problems. The multigrid method of Achi Brandt (1977), originally a standalone solver, was soon recognised as a spectrally optimal preconditioner; the additive-Schwarz domain-decomposition framework was placed on a rigorous condition-number footing through the 1980s and 1990s by Olof Widlund, Maksymilian Dryja, and collaborators. Yousef Saad's textbook consolidated the algebra of left, right, and split preconditioning, the ILU hierarchy, and approximate-inverse methods into the unified account standard today ^{[Saad, Y. — Iterative Methods for Sparse Linear Systems (2nd ed.)]}, and Anne Greenbaum supplied the spectral-equivalence convergence theory that distinguishes optimal from merely useful preconditioners ^{[Greenbaum, A. — Iterative Methods for Solving Linear Systems (SIAM, 1997)]}.

Bibliography Master

@article{MeijerinkVanDerVorst1977,
  author  = {Meijerink, J. A. and van der Vorst, H. A.},
  title   = {An Iterative Solution Method for Linear Systems of Which the Coefficient Matrix is a Symmetric {M}-Matrix},
  journal = {Mathematics of Computation},
  volume  = {31},
  number  = {137},
  year    = {1977},
  pages   = {148--162}
}

@book{Saad2003precond,
  author    = {Saad, Yousef},
  title     = {Iterative Methods for Sparse Linear Systems},
  edition   = {2},
  publisher = {Society for Industrial and Applied Mathematics},
  year      = {2003}
}

@book{Greenbaum1997precond,
  author    = {Greenbaum, Anne},
  title     = {Iterative Methods for Solving Linear Systems},
  series    = {Frontiers in Applied Mathematics},
  publisher = {SIAM},
  year      = {1997}
}

@book{Varga1962,
  author    = {Varga, Richard S.},
  title     = {Matrix Iterative Analysis},
  publisher = {Prentice-Hall},
  year      = {1962}
}

@article{Brandt1977,
  author  = {Brandt, Achi},
  title   = {Multi-Level Adaptive Solutions to Boundary-Value Problems},
  journal = {Mathematics of Computation},
  volume  = {31},
  number  = {138},
  year    = {1977},
  pages   = {333--390}
}

@book{ElmanSilvesterWathen2014,
  author    = {Elman, Howard C. and Silvester, David J. and Wathen, Andrew J.},
  title     = {Finite Elements and Fast Iterative Solvers},
  edition   = {2},
  publisher = {Oxford University Press},
  year      = {2014}
}

@article{Benzi2002,
  author  = {Benzi, Michele},
  title   = {Preconditioning Techniques for Large Linear Systems: A Survey},
  journal = {Journal of Computational Physics},
  volume  = {182},
  number  = {2},
  year    = {2002},
  pages   = {418--477}
}

@book{TrefethenBau1997precond,
  author    = {Trefethen, Lloyd N. and Bau, David},
  title     = {Numerical Linear Algebra},
  publisher = {SIAM},
  address   = {Philadelphia},
  year      = {1997}
}