Matrix preconditioning ($M^{-1/2}A{M}^{-1/2}$)

Matrix preconditioning transforms a linear system using an inexpensive auxiliary solve to improve its eigenvalues. For symmetric positive-definite matrices, a symmetric positive-definite matrix $M$ produces $M^{-1/2}A{M}^{-1/2}$, preserving the symmetric structure needed by the conjugate gradient method.

A useful matrix preconditioner combines a cheap solve with favorable transformed eigenvalues; neither property alone is enough.^[ref]

The aim is to make each correction more effective without making it expensive. The example under diagonal scaling for a two-dimensional linear system turns thin ellipses into nearly circular ones and makes preconditioned steepest descent effective.

A matrix preconditioner is an auxiliary invertible matrix whose matrix inverse is cheap to apply and whose use improves the eigenvalues of a linear system. In symmetric positive-definite matrix problems, a symmetric positive-definite matrix $M$ preserves that structure under the symmetric transformation $M^{-1/2}A{M}^{-1/2}$.

Applying $M^{-1}$ means solving $Mz=r$ for a supplied vector $r$. This notation does not require forming a general matrix inverse. A diagonal matrix gives componentwise divisions; a factored matrix gives auxiliary solves such as forward substitution and back substitution.^[ref]

Choosing $M=I$ makes the solve trivial but changes nothing; choosing $M=A$ gives the exact correction but asks us to solve the original linear system.^[ref]

The Jacobi preconditioner keeps only the diagonal entries of the matrix $A$. For a symmetric positive-definite matrix, these entries are positive, so its auxiliary solve is $z_i=r_i/a_{ii}$. It rescales coordinates with very different diagonal sizes, but cannot remove arbitrary off-diagonal coupling.^[ref]

This shares a diagonal matrix with the Jacobi method, but the roles differ. The Jacobi method specifies a whole stationary iterative method for a linear system; a Jacobi preconditioner supplies an auxiliary solve to another iterative method, including the Preconditioned conjugate gradient method.

Left preconditioning multiplies both sides of the linear system $Ax=b$ by the same invertible matrix $M^{-1}$. It leaves the unknown vector $x$ unchanged. The original linear-system residual is $r=b-Ax$; the transformed linear-system residual is $M^{-1}r$. Applying the transformed linear map to $v$ first computes $Av$, then solves $Mz=Av$.

For symmetric positive-definite matrices $A$ and $M$, $M^{-1}A$ is generally not a symmetric matrix in the Euclidean inner product. Its positive eigenvalues alone do not justify ordinary Euclidean conjugate gradient method recurrences. The remedy is symmetric positive-definite preconditioning, implemented through the Preconditioned conjugate gradient method.

Right preconditioning substitutes $x=M^{-1}u$ into the linear system $Ax=b$. Applying its transformed linear map first solves $Mz=v$, then computes $Az$. Recover the solution by solving $Mx=u$. Its linear-system residual $b-A{M}^{-1}u$ is exactly the original linear-system residual $b-Ax$, which is useful when monitoring the imbalance in the original equations.

The left preconditioning and right preconditioning transformed matrices are similar matrices:They therefore have the same eigenvalues, although their Euclidean norm condition numbers and residual histories need not agree. A similarity transformation preserves eigenvalues; a nonorthogonal change of coordinates need not preserve lengths.

Let $A$ and $M$ be real symmetric positive-definite matrices. The principal square root of a positive semidefinite matrix construction gives an invertible matrix for each. SetThen $Ax=b$ is equivalent to $By=c$. This is a matrix congruence, and $B$ is a symmetric positive-definite matrix: $B^T=B$, and for every nonzero vector $v$,since $S^{-1}v\ne 0$. This short proof establishes both requirements for applying the conjugate gradient method to the transformed linear system.

The left transformed matrix has the same eigenvalues becauseIts correct geometry uses the weighted inner product $⟨u,v{⟩}_M=u^{T}Mv$: $⟨u,M^{-1}Av{⟩}_M=u^{T}Av=⟨M^{-1}Au,v{⟩}_M$. Thus it is symmetric in this inner product, even when it is not a symmetric matrix in ordinary coordinates. The spectral condition number governing the transformed symmetric problem is ${\kappa }_2(B)$, not generally ${\kappa }_2(M^{-1}A)$.

The Preconditioned conjugate gradient method realizes this transformation using auxiliary solves $M{z}_k=r_k$, without constructing $S$ or its matrix inverse. With fixed $M$, initialize $r_0=b-A{x}_0$ and $p_0=z_0$; while $r_k\ne 0$, usethen solve $M{z}_{k+1}=r_{k+1}$ and setStop at a zero linear-system residual before forming another ratio. Positive definiteness makes the denominators positive before termination.^[ref]

An invertible matrix alone is insufficient: with $A=I$ and $M=\mathrm{diag}(-1,1)$, the left transformed matrix has eigenvalues $-1,1$, and $M$ has no real symmetric positive-definite matrix square root. Even a valid symmetric positive-definite matrix can worsen conditioning if its scaling is unsuitable, as the next worked example shows.

For a symmetric positive-definite matrix $A$, solving $A{x}_{\ast }=b$ is equivalent to minimizing the quadratic form with linear termWriting the error vector $e=x-x_{\ast }$ and expanding gives the exact energy norm identityThe cross term vanishes because $A{x}_{\ast }=b$. With $\eta =Se$ and $B=S^{-1}A{S}^{-1}$, the same quantity is ${\eta }^{T}B\eta /2$. Symmetric positive-definite preconditioning changes coordinates without changing the energy assigned to a corresponding error vector.

In an orthonormal eigenbasis of $A$, the level set $e^{T}Ae=t>0$ has semiaxis lengths $\sqrt{t/{\lambda }_i}$. Its longest-to-shortest axis ratio is $\sqrt{{\kappa }_2(A)}$. A large spectral condition number of a positive-definite matrix therefore appears geometrically as a thin ellipsoid. Good scaling makes the transformed ellipsoid rounder, allowing a simple search direction to improve all coordinates together.

The negative gradient of $F$ is the linear-system residual $r$. Relative to the norm $\Vert v{\Vert }_M=\sqrt{v^{T}Mv}$, a steepest direction is $z=M^{-1}r$: the Cauchy-Schwarz inequality applied to $r^{T}v=(S^{-1}r{)}^T(Sv)$ shows that the greatest decrease per unit $M$-length occurs along this direction. Equivalently, ordinary gradient descent in $y=Sx$ becomes preconditioned steepest descent in $x$.^[ref]

Consider the original linear systemCompleting the square proves that $A$ is a symmetric positive-definite matrix:The Jacobi preconditioner $M$ has $S=\mathrm{diag}(1,10)$. Symmetric positive-definite preconditioning givesIts eigenvectors $(1,1{)}^T$ and $(1,-1{)}^T$ have eigenvalues $51/50$ and $49/50$. HenceThe characteristic polynomial of $A$ is ${\lambda }^2-101\lambda +2499/25$. Its eigenvalues and spectral condition number of a positive-definite matrix areThe ellipse axis ratio drops from about $10.0020$ to $1.0202$. The transformation multiplies the second coordinate by ten, bringing the two energy scales together.

Both panels show the same energy level sets, expressed in error coordinates $e=x-x_{\ast }$ and $\eta =M^{1/2}e$. Their coordinate axes have equal scales. The marked displacement $(1,0.1{)}^T$ becomes $(1,1{)}^T$; its energy remains $e^{T}Ae/2=1.02$. These are exact contours of the specified quadratic form.

Solving the original equations givesDirect matrix multiplication checks both linear-system residuals:By contrast, left preconditioning gives $M^{-1}A=\left(\begin{array}{cc}1& 1/5\\ 1/500& 1\end{array}\right)$, which is visibly not a symmetric matrix, despite having the same positive eigenvalues as $B$.

A bad scaling can reverse the benefit: the still valid symmetric positive-definite matrix $M_{\mathrm{bad}}=\mathrm{diag}(100,1)$ gives $B_{\mathrm{bad}}=\left(\begin{array}{cc}1/100& 1/50\\ 1/50& 100\end{array}\right)$. Its Rayleigh quotients on the coordinate unit vectors are $0.01$ and $100$, so ${\kappa }_2(B_{\mathrm{bad}})\ge 10000$. Cheapness and positive definiteness do not guarantee improved conditioning.

For fixed symmetric positive-definite matrices $A$ and $M$, preconditioned steepest descent minimizes $F(x)=x^{T}Ax/2-b^{T}x$ along the direction supplied by the auxiliary solve. Given $x_k$, form the linear-system residual $r_k=b-A{x}_k$. If $r_k=0$, stop: $x_k=x_{\ast }$. Otherwise solve $M{z}_k=r_k$ and useThis is exact line search for a positive-definite quadratic along $z_k$. Indeed,whose derivative vanishes at the displayed positive step size. Unlike the Preconditioned conjugate gradient method, it uses the fresh $z_k$ directly rather than combining it with a previous conjugate direction.

For diagonal scaling for a two-dimensional linear system, start at $x_0=0$. Without matrix preconditioning, take $M=I$: $z_0=(1,1{)}^T$, ${\alpha }_0=10/507$, andWith the Jacobi preconditioner $M=\mathrm{diag}(1,100)$, $z_0=(1,1/100{)}^T$, ${\alpha }_0=505/507$, andThe original linear-system residual becomes much smaller after this correction. Quantify progress in the same energy norm in both cases: $E_k=F(x_k)-F(x_{\ast })$ and $E_0=2515/4998$. Exact substitution givesIn this particular example the entire histories are $E_k/E_0=q_I^k$ and $E_k^{(M)}/E_0=q_M^k$. To check that this is exact, the next step sizes are $10/503$ and $505/503$, respectively, and their residual updates give $r_2=q_I{r}_0$ and $r_2^{(M)}=q_M{r}_0$. Multiplying a linear-system residual by a nonzero scalar leaves the next step size unchanged; every pair of steps repeats the same directions at a smaller scale. The two first-step energy ratios then establish these formulas for all $k$.

The curves come from exact rational updates of the specified linear system, converted to decimal values for plotting. The left panel counts iterations. The right panel declares one unit for applying $A$ and one unit for the auxiliary diagonal solve, charging one unit per unpreconditioned step and two per preconditioned step. It excludes setup and vector operations, so it compares a declared cost model rather than measured runtime. Scaling is advantageous even with that extra charge.

These curves compare steepest descent method iterations. On this two-dimensional linear system, both the conjugate gradient method and the Preconditioned conjugate gradient method terminate within two steps in exact arithmetic. The plot therefore establishes no CG iteration-count advantage for this example; the benefit of scaling for larger problems depends on their eigenvalues and auxiliary-solve costs.

Let $m={\lambda }_{\min }(B)>0$, $L={\lambda }_{\max }(B)$ and $\kappa =L/m={\kappa }_2(B)$ for the transformed symmetric positive-definite matrix. Then preconditioned steepest descent satisfiesThis is a worst-case bound, not necessarily the actual rate. For our scaled example the bound is $(1/50{)}^2=0.0004$, slightly above its computed factor $q_M\approx 0.000384321291$. If the initial transformed error is an eigenvector, exact line search for a positive-definite quadratic removes it in one step, regardless of the full condition number.

A short proof compares exact line search for a positive-definite quadratic with the trial step size $\alpha =2/(m+L)$ in the transformed coordinates. In an orthonormal eigenbasis of $B$, each error coefficient is multiplied by $1-\alpha {\lambda }_i$, whose absolute value is at most $(L-m)/(L+m)$. Summing ${\lambda }_i$ times the squared coefficients bounds the trial energy by the square of this factor. The exact minimizing step size does at least as well. This also explains the terminating case $m=L$: the trial step removes the entire error in one iteration.