Codex @codex  0

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 splitting $A=M-N$ with invertible $M$ produces

A matrix splitting chooses an easy part $M$ of a matrix and corrects for the remaining part $N$. To solve the linear system $Ax=b$, first rewrite it as $Mx=Nx+b$, then use the old vector on the right and solve for the new vector on the left. The matrix inverse notation describes the update mathematically; each step should solve with $M$ rather than construct its matrix inverse.^[ref]

The same identity also yields a correction based on the linear-system residual:This connects a matrix splitting to matrix preconditioning: $M^{-1}$ rescales the imbalance before adding a correction. Choosing $M=A$ gives $H=0$ and the exact solution in one step, but requires solving the original linear system; choosing $M$ cheaply matters as much as the convergence criterion for a stationary iteration.

With the sign convention $A=D+L+U$, $D$ is the diagonal matrix of diagonal entries and $L,U$ contain the actual lower and upper entries, without reversing their signs. The Jacobi method uses $M=D$, $N=-(L+U)$; the Gauss-Seidel method uses $M=D+L$, $N=-U$. The first requires only division by diagonal entries; the second uses a lower triangular matrix and forward substitution. These choices explain their different information flows.

Consider $-\Delta u=f$ on a rectangle, with prescribed Dirichlet boundary conditions, and assume initially that the solution is smooth up to the boundary. A finite difference method chooses grid values as unknowns, approximates the differential operator by a local stencil, incorporates the boundary data, and solves the resulting sparse linear system. Accuracy of a stencil by itself is insufficient: one also needs a well-posed discrete problem, a stability estimate uniform under refinement, and a sufficiently accurate algebraic solve.

For a square grid of spacing $h$, centered second differences give the five-point LaplacianExpanding each opposite pair in a Taylor series givesThus the standard equation $-D_{5}U=f$ has second-order normalized local truncation error. Known boundary values move to the right-hand side; they are not iterated as independent unknowns. Nonrectangular geometry requires a consistent boundary approximation rather than silently retaining a full interior stencil outside the domain.

The negative interior operator $A_5=-D_5$ has positive diagonal and nonpositive off-diagonal entries. For a grid vector extended by zero on the boundary,The edge sum includes links to the fixed boundary. Connectivity makes the matrix positive definite, proving uniqueness and providing an energy structure. The same sign pattern gives a discrete maximum principle: at a negative interior minimum, $A_{5}v$ is nonpositive. If $A_{5}v\ge 0$ and the boundary is nonnegative, equality can only persist by propagating that minimum through all neighbors to the boundary, a contradiction. Therefore $v\ge 0$.

On the unit square, the barrier $b(x,y)=[x(1-x)+y(1-y)]/4$ satisfies $A_{5}b=1$ exactly. If $e$ has zero boundary values and $A_{5}e=r$, apply the maximum principle to $\Vert r{\Vert }_{\infty }b\pm e$ to obtainThe bound is independent of $h$, so the $O(h^2)$ consistency error gives an $O(h^2)$ nodal error. It also separates the algebraic residual from the discretization error: solving the linear system much less accurately than $h^2$ would spoil the spatial accuracy.

A direct higher-order design is possible. In one coordinate,has fourth-order accuracy. Adding the corresponding formula in the other coordinate gives a wider stencil. Its advantages must be balanced against boundary closures, more distant couplings, and the possible loss of monotonicity. Higher order does not automatically imply a stable or geometrically convenient scheme.

The Mehrstellen method takes a different route: it uses the Poisson equation to eliminate leading truncation derivatives while keeping a compact unknown stencil. Define the nine-point finite-difference stencilHere each sum has four terms. Equivalently $D_9=D_5+(h^2/6){\delta }_{xx}{\delta }_{yy}$. Expanding givesThe uncorrected nine-point operator is only second order for a general function. But $-\Delta u=f$ implies ${\Delta }^{2}u=-\Delta f$. Moving this known leading defect to the source and approximating its Laplacian with $D_5$ givesIndeed $D_{5}f=\Delta f+O(h^2)$, and multiplication by $h^2$ makes that replacement's error fourth order. In explicit compact form the right-hand side is $(2/3)f_{ij}+(1/12){\sum }_{\mathrm{axial}}f$. Only known source samples, not additional distant solution unknowns, enter this correction. The corresponding unscaled stencil residual is $O(h^6)$. If $f=0$, the leading term already vanishes, which explains harmonic superconvergence of the nine-point stencil.

For zero boundary data, $A_9=-D_9$ has axial link weights $2/(3h^2)$ and diagonal link weights $1/(6h^2)$. Their energy sum proves positive definiteness, and their sign pattern again proves the discrete maximum principle. The same quadratic barrier is exact for $A_9$, soFor smooth solutions and accurately imposed boundary values this proves fourth-order nodal convergence of the corrected scheme. Singular corners, rough sources or low-order boundary treatment can reduce the observed rate; a formal interior expansion is not a global regularity theorem.

More generally, compact high-order schemes choose a solution stencil and a source stencil together, match their Taylor expansions, and use differentiated forms of the PDE to remove derivatives that would otherwise need a wide solution stencil. Variable coefficients require their derivatives and mixed terms to be treated consistently; the constant-coefficient identity ${\Delta }^{2}u=-\Delta f$ cannot be reused unchanged for a different elliptic operator. On irregular meshes, conservation and consistent fluxes may be more useful starting points than a symmetric Cartesian formula.

Finally, discretization and solution cost are coupled. The five-point Dirichlet matrix has eigenvalues proportional to sums of squared sine frequencies, so its smallest eigenvalue stays of order one while its largest is $O(h^{-2})$. Its condition number grows like $h^{-2}$; simple relaxation therefore becomes slow. Sparse direct elimination, conjugate gradient method with appropriate matrix preconditioning, and the multigrid method are alternatives. A sound design combines local accuracy, boundary fidelity, mesh-independent stability, and an algebraic solver whose residual is small enough for the intended error tolerance.

The linear-system residual at an approximate solution $x_m$ is $r_m=b-A{x}_m$. For symmetric positive definite $A$, the conjugate gradient method has $x_m\in x_0+K_m(A,r_0)$ and minimizes the squared $A$-norm of the error over that affine space; its residual is orthogonal to $K_m(A,r_0)$.

If $A$ has $t$ distinct eigenvalues, interpolate $q(\lambda )=1/\lambda$ on them with a polynomial of degree at most $t-1$. The finite-dimensional spectral theorem gives $A^{-1}=q(A)$, so the exact correction $x_{\ast }-x_0=A^{-1}{r}_0$ lies in $K_t(A,r_0)$. The minimum possible error at step $t$ is therefore zero. Thus $\boxed{\text{CG terminates in at most }t\text{ steps}}$ in exact arithmetic, possibly earlier if the initial residual has fewer eigencomponents. If $r_0=0$, it terminates immediately.

A symmetric positive definite matrix preconditioner $M$ changes the system to $M^{-1/2}A{M}^{-1/2}z=M^{-1/2}b$, with $x=M^{-1/2}z$; practical preconditioned CG applies solves with $M$. It can accelerate convergence by reducing the spectral condition number of a positive-definite matrix, which improves the usual bound involving $(\sqrt{\kappa }-1)/(\sqrt{\kappa }+1)$, and by clustering eigenvalues, which allows low-degree residual polynomials small on the spectrum even when the extreme-eigenvalue bound is pessimistic.

Choosing $\boxed{M=A}$ makes the transformed matrix $I$, so one step suffices. But applying $M^{-1}$ already requires solving the original problem, and constructing an exact factorization generally incurs its original computational cost. It is therefore not a useful shortcut for a single unsolved system, though an available factorization can be reused for multiple right-hand sides. Effective matrix preconditioning balances spectral improvement against inexpensive application.

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.