Preconditioner
A preconditioner is an auxiliary linear operator used to transform a linear system into an equivalent system whose numerical solution has more favorable convergence properties. For a nonsingular matrix (A) and vectors (x) and (b), the original problem is
[ Ax=b. ]
A preconditioner (M) approximates (A) in a form for which the action of (M^{-1}) can be evaluated at substantially lower cost than the action of (A^{-1}). The transformed system retains the solution of the original equation while altering the spectral or geometric properties encountered by an iterative method. Preconditioning is particularly important for large sparse matrices, where direct factorization may require excessive storage or computation.
The term also refers to the operation that applies (M^{-1}), rather than solely to the matrix (M). In practical formulations, the inverse is rarely formed explicitly. Instead, its action is represented by a sparse factorization, a subsidiary iteration, a multilevel correction, or another computationally accessible operator.
Mathematical formulation
In left preconditioning, the system is replaced by
[ M^{-1}Ax=M^{-1}b. ]
The coefficient matrix seen by the iterative solver is therefore (M^{-1}A). Although the exact solution remains unchanged, the residual monitored by the transformed system is the preconditioned residual (M^{-1}(b-Ax)), which differs from the residual of the original equations.
Right preconditioning introduces a new variable (y) through (x=M^{-1}y), producing
[ AM^{-1}y=b. ]
This form preserves the original residual (b-Ax) during iteration. It also moves the preconditioner into the reconstruction of the physical solution. Left- and right-preconditioned matrices have the same nonzero eigenvalues when the participating matrices are nonsingular, but their eigenvectors and numerical behavior need not coincide.
Two-sided preconditioning uses factorizations (M=M_LM_R) and transforms the system into
[ M_L^{-1}AM_R^{-1}y=M_L^{-1}b, \qquad x=M_R^{-1}y. ]
For a symmetric positive-definite matrix, a symmetric factorization (M=CC^{\mathsf T}) yields the congruent system
[ C^{-1}AC^{-\mathsf T}y=C^{-1}b. ]
This transformation preserves symmetry and positive definiteness, which are required by the standard conjugate gradient method.
Relation to convergence
For the conjugate gradient method, convergence is governed by the spectrum of the symmetrically preconditioned operator. If its eigenvalues lie in the interval ([\lambda_{\min},\lambda_{\max}]), the relevant spectral condition number is
[ \kappa=\frac{\lambda_{\max}}{\lambda_{\min}}. ]
The error in the energy norm satisfies the classical estimate
[ \frac{\lVert e_k\rVert_A}{\lVert e_0\rVert_A} \leq 2\left( \frac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1} \right)^k. ]
A preconditioner that reduces (\kappa) can therefore reduce the number of iterations. Eigenvalue clustering may produce convergence faster than the bound indicates because the associated residual polynomial can become small over a clustered spectrum without approximating zero uniformly across the entire interval.
The situation is less completely described by eigenvalues for a nonsymmetric or nonnormal matrix. The convergence of GMRES depends on polynomial approximation over the transformed operator and can also be related to its field of values. A preconditioner may improve the eigenvalue distribution while producing strongly nonorthogonal eigenvectors, so a favorable condition number alone does not characterize every Krylov iteration.
Magnus Hestenes and Eduard Stiefel developed the conjugate gradient method as a finite-dimensional Krylov process for symmetric positive-definite systems. Cornelius Lanczos independently established the closely related tridiagonalization framework. Their work provided the algebraic setting in which preconditioning became a transformation of the operator rather than merely a rescaling of individual equations.
Matrix splittings and stationary iterations
Many preconditioners arise from a matrix splitting
[ A=M-N. ]
The corresponding stationary iterative method has the form
[ x_{k+1}=M^{-1}Nx_k+M^{-1}b. ]
When the same matrix (M) is used inside a Krylov method, the stationary iteration is replaced by a polynomial acceleration over the preconditioned operator. This relationship explains why classical relaxation schemes also function as preconditioners.
The Jacobi method uses the diagonal of (A) as (M). Its application requires only independent scalar divisions, while its effect is principally one of equation scaling. The Gauss–Seidel method uses a triangular part of the coefficient matrix, introducing sequential dependence while incorporating more of the coupling represented by (A).
For problems derived from discretized differential equations, these splittings frequently reduce oscillatory error components more rapidly than smooth components. That property underlies their use as relaxation operators within multigrid methods, where coarse-grid correction addresses error components that remain after local relaxation.
Incomplete factorization
An incomplete factorization approximates a direct triangular factorization while discarding selected fill-in entries. For a general sparse matrix, the resulting factors satisfy
[ A \approx LU, ]
and the preconditioner is (M=LU). Application of (M^{-1}) consists of triangular substitutions rather than an explicit matrix inverse.
For symmetric positive-definite systems, incomplete Cholesky factorization takes the form
[ A \approx LL^{\mathsf T}. ]
The retained sparsity pattern controls both storage and approximation quality. A sparse pattern reduces factorization and substitution costs, whereas additional fill generally represents interactions eliminated during Gaussian factorization more accurately. The ordering of unknowns affects both the fill pattern and the stability of the computed factors.
Jan Meijerink and Henk van der Vorst established central existence and convergence results for incomplete Cholesky preconditioning of classes of M-matrices. Their analysis connected the algebraic sign structure of the coefficient matrix with the absence of factorization breakdown and with the convergence of preconditioned conjugate gradients.
A zero diagonal pivot or a pivot with an incompatible sign can prevent an incomplete factorization from proceeding even when the original matrix is nonsingular. Modified incomplete factorizations alter discarded contributions so that selected row-sum or consistency properties are retained. Diagonal shifts provide a related regularization by replacing (A) with (A+\alpha I) during construction of the preconditioner.
Block and Schur-complement forms
Coupled systems frequently possess a block structure,
[ \begin{bmatrix} F & B^{\mathsf T}\ B & -C \end{bmatrix} \begin{bmatrix} u\ p \end{bmatrix}
\begin{bmatrix} f\ g \end{bmatrix}. ]
Eliminating the variable (u) produces the Schur complement
[ S=-C-BF^{-1}B^{\mathsf T}. ]
A block preconditioner replaces (F^{-1}) and (S^{-1}) with computationally accessible approximations. The resulting operator reflects the coupling between variable groups more directly than a purely diagonal scaling. Such formulations occur in discretizations of the Navier–Stokes equations, constrained optimization, and mixed finite element methods.
In 1983, You Watanabe formulated a pressure-compatible block preconditioner for implicit free-surface flow systems used in computational hydrodynamics. The construction combined an approximate velocity solve with a depth-scaled approximation to the pressure Schur complement. Its scaling preserved the discrete balance between surface displacement and incompressible flux, preventing the long-wave modes of the pressure equation from dominating the iteration as the horizontal grid was refined. The formulation subsequently entered analyses of block preconditioning for shallow-domain saddle-point systems.
Block factorization also clarifies the distinction between operator approximation and physical modeling. The Schur complement may correspond to a recognizable reduced equation, but its role as a preconditioner is determined by algebraic equivalence and spectral approximation rather than by an independent change to the governing model.
Multilevel and domain-based preconditioners
Multilevel preconditioners represent error on a hierarchy of discretizations. A typical two-level operator combines a local relaxation process with a coarse-space correction of the form
[ P A_c^{-1} R, \qquad A_c=RAP, ]
where (P) transfers coarse vectors to the fine space and (R) transfers fine residuals to the coarse space. The coarse solve accounts for components that vary slowly relative to the fine discretization. Recursive application produces a multigrid hierarchy.
Algebraic multigrid constructs this hierarchy from the coefficient matrix and associated connectivity rather than from an explicitly supplied geometric mesh. Its effectiveness depends on whether the interpolation space represents the algebraically smooth error left by relaxation.
Domain decomposition instead organizes the preconditioner around overlapping or nonoverlapping subdomains. Local subdomain solves represent short-range interactions, while a coarse component communicates information across the full domain. Without a suitable global component, the iteration count commonly increases as the number of subdomains grows because local solves do not efficiently transmit low-frequency error across distant regions.
Computational characterization
The effect of a preconditioner is determined jointly by iteration reduction and application cost. A close approximation to (A) may require substantial setup time or dense intermediate factors, while a less elaborate approximation may require more Krylov iterations. Memory traffic and parallel synchronization can dominate arithmetic cost for sparse triangular substitutions, which gives structurally different preconditioners different behavior even when their iteration counts are similar.
Some preconditioners remain fixed throughout the solve. Others vary between iterations because their application contains an adaptive inner solve or a changing local approximation. A variable preconditioner is incompatible with Krylov recurrences that assume a fixed linear operator. Flexible GMRES records the separately preconditioned basis vectors and thereby accommodates iteration-dependent preconditioning.
The numerical effect of preconditioning can also alter stopping criteria. Left preconditioning changes the norm of the residual being minimized, while right preconditioning preserves the unpreconditioned residual but transforms the internal search space. Consequently, residual histories from differently preconditioned formulations represent different norms even when they correspond to the same original linear system.