1 Problem setting and Krylov subspaces

Krylov subspace methods solve linear algebra problems by searching for an approximate solution within a sequence of nested subspaces generated from a matrix acting on a starting vector. The key idea is that repeated matrix–vector products can reveal useful directions in which the error decreases, while the computation is confined to a smaller space where projections are inexpensive.

1.1 Linear systems and residuals

Consider a linear system \[ Ax=b, \] where \(A\) is an \(n\times n\) matrix (or a linear operator) and \(b\) is a given vector. An iterative method produces a sequence of approximations \(x_0, x_1, \dots\). The quality of an iterate is measured by its residual \[ r_k=b-Ax_k. \] Many Krylov methods build the next approximation so that the residual becomes small in a prescribed sense (e.g., minimized norm over a subspace).

1.2 Definition of the Krylov subspace \(K_m(A,v)\)

A Krylov subspace of dimension \(m\) associated with \(A\) and a vector \(v\) is defined as \[ K_m(A,v)=\text{span}\{v,Av,A^2v,\dots,A^{m-1}v\}. \] In practice, \(v\) is often chosen as the initial residual \(r_0=b-Ax_0\) or as a related vector after preconditioning.

1.2.1 Krylov basis vectors \(v, Av, \dots, A^{m-1}v\)

The subspace is generated by successive applications of the operator to the starting vector. Although these raw vectors are generally not orthogonal, algorithms typically replace them by an orthonormal (or bi-orthogonal) basis to obtain stable projections and reduced problems of manageable size.

1.2.2 Dimension, rank, and breakdown concepts

The subspace dimension can be smaller than \(m\) if the vectors become linearly dependent. This can occur because of exact arithmetic relationships or numerical effects. In algorithmic terms, this manifests as breakdown (a step where a needed vector direction vanishes or cannot be normalized). Near-breakdown can also arise when the vectors become almost dependent, degrading numerical behavior.

1.3 Polynomial viewpoint and moment interpretation

Krylov methods can be interpreted through polynomials in \(A\). If \(x\) is sought in the affine space \(x_0+K_m(A,r_0)\), then the error \(e_m=x-x_\ast\) (with \(x_\ast\) the exact solution) is often expressible via a polynomial filter: \[ e_m = p_m(A)e_0, \] for some polynomial \(p_m\) satisfying \(p_m(0)=1\). The residual similarly relates to such polynomials, which leads to approximation-theoretic interpretations of convergence.

1.4 Connections to projections and reduced models

Rather than solving the full problem in \(n\)-dimensional space, Krylov methods solve a smaller projected system or enforce conditions on the residual via projection. A typical pattern is:

  1. build a basis \(V_m\) for \(K_m(A,v)\),
  2. form a reduced operator \(A_m = V_m^\ast A V_m\) (or a related two-sided form),
  3. compute an approximate solution in the span of the basis by imposing a Galerkin or least-squares condition.

This yields reduced computational work per iteration, particularly when the reduced system is small relative to \(n\).

2 Building Krylov bases in practice

Constructing a Krylov basis is central to the effectiveness of the method. Since orthogonality properties strongly influence numerical stability and the structure of the reduced problem, the choice of basis-building process shapes both speed and accuracy.

2.1 Orthogonalization strategies

Most modern Krylov solvers rely on orthogonalizing the iteratively generated vectors to maintain a well-conditioned basis.

2.1.1 Arnoldi process for general matrices

For a general (possibly nonsymmetric) operator \(A\), the Arnoldi process builds an orthonormal basis \(V_m=[v_1,\dots,v_m]\) for \(K_m(A,v_1)\) and produces an upper Hessenberg matrix \(H_m\) such that \[ AV_m = V_{m+1}\,\bar{H}_m. \] Here, the relation captures how the operator maps the subspace into the next one, enabling projection-based residual computations.

2.1.2 Lanczos process for Hermitian/symmetric matrices

When \(A\) is Hermitian (or real symmetric), the Lanczos process produces a basis with far simpler structure: the reduced matrix becomes tridiagonal. This is advantageous for storage and speed, and it underpins methods like Conjugate Gradient and MINRES. As in Arnoldi, the next basis vector is created by applying \(A\) and orthogonalizing against prior vectors, but symmetry allows recurrence relations that avoid full orthogonalization.

2.2 Bi-orthogonal and two-sided variants

Some Krylov methods employ two distinct subspaces built from \(A\) and \(A^\ast\). This is particularly important for nonsymmetric problems where single-sided orthogonality may not produce desirable residual properties.

2.2.1 Bi-orthogonal Arnoldi basics

In bi-orthogonal Arnoldi, one constructs bases \(V_m\) and \(W_m\) such that \(W_m^\ast V_m = I\). The reduced operator then takes the form \(W_m^\ast A V_m\), which governs the projected problem. This two-sided construction supports algorithms like BiCG-type methods.

2.2.2 Stabilization considerations

Two-sided processes are more sensitive to numerical errors because orthogonality can be lost independently in both left and right bases. Stabilization may involve controlled reorthogonalization, careful normalization, or algorithmic variants designed to reduce the propagation of rounding errors.

2.3 Numerical stability and loss of orthogonality

In exact arithmetic, orthogonality properties hold by construction. In floating-point arithmetic, rounding can erode orthogonality, affecting convergence behavior and sometimes causing spurious stagnation.

2.3.1 Reorthogonalization heuristics

A common response is reorthogonalization: during basis construction, additional passes of orthogonalization are applied when orthogonality tests indicate degradation. Practical implementations often use modified Gram–Schmidt with selective or repeated orthogonalization to balance cost and robustness.

2.3.2 Impact on convergence and accuracy

When basis vectors lose orthogonality, the reduced problem no longer accurately represents the projected dynamics, and residual minimization may fail to behave as expected. Symptoms include irregular residual histories, unexpectedly slow progress, or convergence to an incorrect target in ill-conditioned cases.

3 Krylov subspace methods for linear systems

Many Krylov methods differ primarily in how they define projection conditions and how they generate short or long recurrences. The choice is linked to whether the operator is symmetric, nonsymmetric, or indefinite.

3.1 Minimal residual and Galerkin-type frameworks

Projection frameworks specify conditions under which an approximate iterate is chosen.

3.1.1 Projection conditions (Galerkin vs least-squares)

A Galerkin condition typically enforces orthogonality between the residual and the test space, such as \[ V_m^\ast r_m = 0. \] Least-squares formulations choose an iterate that minimizes the residual norm over the subspace. These two ideas lead to different algorithms and different stability profiles.

3.1.2 Residual minimization in subspaces

In minimal residual methods, the iterate is determined by solving a small optimization problem in the reduced space. This approach often yields strong convergence properties, though it may require more orthogonalization and storage than methods based on short recurrences.

Generalized Minimal Residual (GMRES) is a widely used method for nonsymmetric systems, designed around residual norm minimization.

3.2.1 Restarted GMRES (GMRES(m))

GMRES builds a Krylov basis and at each step computes the best residual within the expanding subspace. However, because memory and orthogonalization cost grow with subspace dimension, practical solvers often use restarts: after \(m\) inner iterations, the method discards the basis and begins again from the latest approximation. This is known as restarted GMRES.

3.2.2 Preconditioned GMRES

With preconditioning, the method applies to a transformed system that ideally clusters eigenvalues and improves conditioning. Preconditioning can be used in left, right, or split forms, each affecting how the residual is defined and how the reduced problems are assembled.

3.3 Conjugate Gradient (CG) and short-recurrence behavior

Conjugate Gradient targets symmetric positive definite (SPD) systems and is notable for short recurrences that reduce per-iteration cost.

3.3.1 CG in exact arithmetic

For SPD matrices, CG generates search directions that are mutually \(A\)-conjugate, enabling efficient updates without full basis storage. In exact arithmetic, the method converges in at most \(n\) steps, and the error after \(k\) iterations can be characterized by minimizing an \(A\)-norm error over a Krylov subspace.

3.3.2 CG for SPD and SPD-with-preconditioning

In practice, CG is run on preconditioned systems to improve convergence. When the preconditioner preserves symmetry and positive definiteness in the appropriate sense, CG-like guarantees extend, and the algorithm benefits from the same short-recurrence efficiencies.

3.4 BiCG, BiCGSTAB, and nonsymmetric Krylov solvers

Nonsymmetric linear systems often motivate methods that use two-sided information to maintain residual relations.

3.4.1 BiCG iterations and shadow residuals

Biconjugate Gradient (BiCG) uses a pair of Krylov subspaces generated by \(A\) and \(A^\ast\). It maintains bi-orthogonality between the residual and a “shadow” residual associated with the adjoint problem, producing recurrences that avoid full basis storage.

3.4.2 BiCGSTAB smoothing and variants

BiCGSTAB (Biconjugate Gradient Stabilized) introduces a stabilization step intended to reduce oscillations common in BiCG. By blending iterates from the BiCG process with additional polynomial smoothing, BiCGSTAB often yields more reliable residual behavior for challenging nonsymmetric problems.

3.5 MINRES and other symmetric-indefinite solvers

For symmetric matrices that are not positive definite, MINRES is a common choice.

3.5.1 MINRES projection rationale

MINRES seeks an approximation that minimizes the residual norm (or an equivalent quantity) over a Krylov subspace, while being compatible with symmetric indefinite operators. It relies on Lanczos-type tridiagonalization and solves reduced problems aligned with the minimization goal.

3.5.2 Handling indefinite operators

Indefiniteness can break assumptions required by CG, but MINRES remains applicable because its formulation does not require \(A\) to be positive definite. Convergence still depends on spectral features such as how eigenvalues are distributed and how well the Krylov subspace approximates the needed components.

4 Preconditioning and acceleration

Preconditioning changes the problem into a form where iterative methods converge more quickly, often by improving conditioning and altering spectral properties relevant to Krylov approximations.

4.1 The role of preconditioners

A preconditioner \(M\) approximates the action of \(A\) (or its inverse) and is used to solve a transformed system that is mathematically equivalent under appropriate forms.

4.1.1 Left vs right preconditioning

Left preconditioning typically transforms the system to \(M^{-1}Ax = M^{-1}b\), affecting how residuals relate to the original problem. Right preconditioning instead rewrites the unknown, such as \(A M^{-1} y = b\) with \(x=M^{-1}y\). Each choice influences monitoring and sometimes practical stability and stopping rules.

4.1.2 Effect on the spectrum and subspace geometry

Although preconditioning does not guarantee uniform improvement in every case, it often improves eigenvalue distribution in a way that makes polynomial approximations easier. Geometrically, it can reduce the angle between the current residual and directions where the Krylov subspace is effective.

4.2 Common preconditioner types

Preconditioners range from algebraic approximations to multilevel schemes.

4.2.1 Incomplete factorizations (IC/ILU) overview

Incomplete LU (ILU) and related incomplete factorizations attempt to factor \(A\) while dropping selected entries to limit fill-in. The quality of the approximation can strongly influence convergence, with parameter choices controlling the balance between accuracy and computational cost.

2.2.2 Multigrid and coarse-grid correction (high level)

Multigrid methods use a hierarchy of discretizations to eliminate errors at multiple scales. While classical multigrid can be used as a solver, its components also serve as preconditioners, accelerating Krylov methods by reducing low-frequency error quickly.

2.2.3 Block and domain decomposition approaches

For systems with natural block structure or spatial decomposition, block preconditioners or domain decomposition techniques apply approximate solves on subproblems. These can be especially useful in parallel computing contexts, where localized operations reduce communication overhead.

4.3 Practical guidelines for choosing preconditioners

Selection is driven by cost, robustness, and the ability to build and apply the preconditioner efficiently.

4.3.1 Trade-offs: setup cost vs iteration savings

Some preconditioners require substantial setup, such as building ILU factors or constructing multilevel operators. If setup is expensive, total runtime can be dominated by preprocessing; conversely, weak preconditioners may lead to many Krylov iterations. Effective choices minimize end-to-end cost rather than iteration count alone.

4.3.2 Tuning and monitoring performance

Monitoring typically includes tracking residual reduction per iteration, total work per solve, and sensitivity to tolerance and problem size. If performance stalls, one may adjust preconditioner parameters, change restart length (for restarted GMRES), or switch solver families.

5 Convergence theory and error bounds

Convergence analysis describes why Krylov methods tend to reduce error and how problem structure influences the rate.

5.1 Residual polynomials and best-approximation ideas

Because residuals in many Krylov methods correspond to polynomials applied to \(A\), convergence relates to how well those polynomials approximate zero on relevant spectral sets.

5.1.1 Chebyshev-type intuition (outline)

For certain matrix classes, the best polynomial approximation to zero on an interval can be linked to Chebyshev polynomials. This yields intuition that when the spectrum lies in a favorable region (e.g., clustered or contained in a bounded interval), residual norms can decrease rapidly.

5.1.2 Field of values considerations (for nonsymmetric cases)

For nonsymmetric problems, eigenvalues alone may not determine behavior. The field of values (numerical range) can provide bounds on how polynomials behave under \(A\), guiding expectations about convergence even when the operator is non-normal.

5.2 Spectral properties and cluster effects

Spectral distribution influences how quickly the Krylov subspace captures error components.

5.2.1 Eigenvalue clustering and faster convergence

When eigenvalues cluster, polynomial filters can annihilate clustered components effectively. This often results in faster convergence because fewer effective polynomial degrees are needed to reduce residual in the dominant eigendirections.

5.2.2 Jordan structure and non-normality

If \(A\) is defective or highly non-normal, eigenvectors may be poorly conditioned and transient growth can occur. Jordan structure affects the minimal polynomial and therefore the degree needed for exact annihilation in exact arithmetic, while non-normality can cause slower or more erratic convergence.

5.3 Stopping criteria

Because computations are approximate, practical stopping rules determine when to halt iterations.

5.3.1 Relative residual norms

A common approach uses a tolerance on the relative residual norm, such as \(\|r_k\|/\|b\|\) or \(\|r_k\|/\|r_0\|\). This ties termination to how close the current iterate is to satisfying the linear system.

5.3.2 Normwise vs componentwise stopping

Normwise stopping focuses on overall magnitude. Componentwise checks can be stricter in cases where some entries are important, but they are more sensitive to scaling and may require problem-specific thresholds.

5.4 Breakdown and near-breakdown

Some Krylov methods rely on divisions or normalizations that can fail if certain vectors become zero or nearly so.

5.4.1 Happy breakdown in certain methods

In exact arithmetic, a breakdown can sometimes coincide with having already reached an exact solution (a “happy breakdown”). This occurs when the constructed residual or correction direction vanishes for reasons aligned with success.

5.4.2 Remedies and safeguards

Safeguards include monitoring norms used in normalization, switching to an alternative method if instability is detected, or using regularization and reorthogonalization strategies. In production settings, robust fallbacks are important because exact arithmetic assumptions rarely hold.

Krylov methods extend naturally from linear systems to eigenvalue, singular value, and related approximation tasks by reusing projection ideas.

6.1 Krylov subspace methods for eigenvalue computation

Eigenvalue solvers often compute approximate eigenpairs from projected operators derived from Krylov bases.

6.1.1 Rayleigh–Ritz procedure

Given a basis \(V_m\), one forms a reduced eigenproblem \(V_m^\ast A V_m y = \theta y\). The approximate eigenvectors in the full space are then \(V_m y\), and the scalars \(\theta\) approximate true eigenvalues. Accuracy improves as \(m\) increases and as the basis captures invariant subspaces.

6.1.2 Implicitly restarted eigensolvers (overview)

To control cost, eigensolvers frequently use restart mechanisms that retain approximate eigencomponents while discarding less useful basis directions. Such approaches combine Arnoldi or Lanczos steps with periodic compression to keep storage and operations manageable.

6.2 Arnoldi-based eigenvalue extraction

For nonsymmetric operators, Arnoldi is often the backbone for computing eigenvalues.

6.2.1 Harmonic Ritz values (high-level)

Harmonic Ritz extraction targets eigenvalues near a specified shift by projecting a related problem using the subspace basis and its companion relations. This can be advantageous when one seeks interior eigenvalues rather than those at the spectrum’s extremes.

6.3 Singular value problems and subspace reuse

Singular values involve the eigenvalues of \(A^\ast A\) or \(A A^\ast\), but Krylov methods typically avoid forming these products explicitly.

6.3.1 Lanczos bidiagonalization overview

Lanczos bidiagonalization builds two sequences of vectors—one for \(A\) action and one for \(A^\ast\) action—producing a bidiagonal reduced matrix. The singular values can then be approximated from that reduced problem.

6.4 Deflation and subspace augmentation

Deflation removes components of the solution associated with already-identified eigenmodes or nearly invariant subspaces. In practice, subspace augmentation enlarges the Krylov space with previously computed vectors to accelerate subsequent solves or eigenvalue iterations.

7 Implementation considerations

Real performance depends on how algorithms are coded: how operators are applied, how vectors are stored, and how reduced problems are solved and updated.

7.1 Data structures and operator interfaces

Implementations typically treat \(A\) as a linear operator with a matrix-vector product routine rather than an explicit dense matrix. This enables Krylov methods to operate on sparse matrices and even on implicit operators arising from discretized physics or data models.

7.2 Complexity and memory usage

Cost includes both arithmetic and communication, especially in parallel environments.

7.2.1 Per-iteration costs by method family

Methods like GMRES require orthogonalization against a growing basis, making per-iteration cost depend on the subspace dimension. CG and Lanczos-based methods usually have cheaper per-step updates due to short recurrence relations. Two-sided methods such as BiCG maintain additional vectors and may increase overhead.

7.2.2 Storage growth with subspace size

Long-recurrence methods store basis vectors to reconstruct the projected problem and update residuals. Restarted variants trade additional iterations for reduced storage. Efficient implementations use compressed representations or careful reuse of memory buffers.

7.3 Practical aspects of matrix-vector products

In many applications, matrix-vector products dominate runtime. Efficient sparse storage formats, optimized kernels, and minimizing data movement strongly affect overall speed. Preconditioning also requires its own operator applications, so performance engineering often focuses on both \(A\) and \(M^{-1}\).

7.4 Handling stopping, restarts, and limits

Termination and restart policies influence both correctness and runtime.

7.4.1 Restart frequency selection

For restarted GMRES, the inner restart size \(m\) controls memory and effectiveness. Too small \(m\) can lead to repeated loss of progress, while too large \(m\) can make each cycle expensive. Selection often depends on problem size, desired tolerance, and hardware constraints.

7.4.2 Safeguards for numerical issues

Common safeguards include checking for breakdown indicators, monitoring orthogonality measures when orthogonalization is selective, and ensuring that reductions and updates do not produce NaNs or infinities. Robust error handling is important for large-scale runs.

8 Worked examples and comparison

Worked examples illustrate how different Krylov solvers behave under matrix type and how empirical evaluation guides method selection.

8.1 Small test problems (illustrative)

Small problems enable direct observation of convergence behavior, residual curves, and breakdown risks. Typical illustrative matrices vary in symmetry, definiteness, and non-normality to demonstrate how solver assumptions affect performance.

8.2 Comparing Krylov methods by matrix type

Solver choice often follows matrix structure.

8.2.1 SPD vs nonsymmetric vs symmetric-indefinite

For SPD systems, CG (with appropriate preconditioning) usually provides efficient short-recurrence convergence. Nonsymmetric problems are often addressed with GMRES or BiCGSTAB variants, while symmetric-indefinite systems are frequently handled by MINRES. In each case, the method’s projection framework aligns with the operator’s structure.

8.3 Interpreting iteration histories

Residual history plots help diagnose issues beyond final convergence.

8.3.1 Residual curves and stagnation patterns

Rapid initial decrease followed by plateau can indicate an ineffective preconditioner or loss of orthogonality. Oscillatory patterns may signal nonsymmetric dynamics or stabilization needs. Sudden drops can indicate successful polynomial annihilation of prominent error components.

8.4 Benchmarks and empirical evaluation (how-to)

Empirical evaluation typically compares total runtime, number of matrix-vector products, and achieved residual tolerances. For fairness, experiments account for preprocessing cost (e.g., ILU setup), communication overhead, and differences in restart sizes or stopping rules. Reproducible runs and standardized test cases help isolate algorithmic effects.

9 Summary and further reading

Krylov subspace methods form a versatile toolbox for large-scale linear algebra, with families tailored to symmetry, indefiniteness, and nonsymmetric behavior.

9.1 Key takeaways by method family

CG and Lanczos-based methods exploit symmetry and can provide efficient short recurrences for SPD or symmetric-indefinite problems. GMRES and minimal-residual approaches handle nonsymmetric systems by minimizing residual norms in projected spaces. BiCG and BiCGSTAB extend Krylov ideas using two-sided information and stabilization.

9.2 Suggested references and standard texts

Standard references include monographs and lecture notes on iterative methods, Krylov subspace approximations, and numerical linear algebra. These typically cover projection theory, stability considerations, preconditioning strategies, and practical solver selection.

9.3 Glossary of common terms

Key terms include Krylov subspace, residual, projection (Galerkin/least-squares), Arnoldi/Lanczos processes, preconditioner, restart, residual minimization, breakdown, and deflation.