Krylov subspace
A Krylov subspace is a sequence of nested linear spaces generated by repeatedly applying a linear operator to a fixed vector. For a matrix (A\in\mathbb{F}^{n\times n}) and a vector (b\in\mathbb{F}^n), the Krylov subspace of order (m) is
[ \mathcal{K}_m(A,b)
\operatorname{span}{b,Ab,A^2b,\ldots,A^{m-1}b}. ]
Here (\mathbb{F}) is usually the field of real or complex numbers. Krylov subspaces provide the approximation spaces used by many iterative methods for large systems of linear equations, eigenvalue problems, and matrix functions. Their significance arises from the fact that vectors in (\mathcal{K}_m(A,b)) can be generated through matrix–vector products without forming powers of (A) explicitly.
The spaces satisfy
[ \mathcal{K}_1(A,b)\subseteq \mathcal{K}_2(A,b)\subseteq\cdots\subseteq \mathcal{K}_n(A,b), ]
although an inclusion ceases to be strict when a newly generated vector is linearly dependent on its predecessors. If (d) is the smallest index for which this occurs, then (d) is the degree of the minimal polynomial of (A) relative to (b). Consequently,
[ \mathcal{K}d(A,b)=\mathcal{K}{d+1}(A,b)=\cdots. ]
This stabilization is exact in symbolic arithmetic. In floating-point arithmetic, numerical dependence and loss of orthogonality complicate the corresponding computational interpretation.
Algebraic structure
Every vector (x_m\in\mathcal{K}_m(A,b)) has the form
[ x_m=p_{m-1}(A)b, ]
where (p_{m-1}) is a polynomial of degree at most (m-1). Krylov approximation is therefore closely connected with polynomial approximation. When (x_m) approximates the solution of (Ax=b), its residual can be written as
[ r_m=b-Ax_m=q_m(A)b, ]
where (q_m) is a polynomial satisfying (q_m(0)=1). The convergence of a Krylov method is governed by how small such a polynomial can become on the spectral or pseudospectral region relevant to (A).
This polynomial description explains why eigenvalue distribution alone does not always determine convergence. For a normal matrix, polynomial behavior on the spectrum provides direct norm estimates. For a strongly nonnormal matrix, eigenvector conditioning and the pseudospectrum can substantially affect the action of (q_m(A)).
Krylov subspaces are invariant under affine shifts of the operator. For every scalar (\sigma),
[ \mathcal{K}_m(A,b)=\mathcal{K}_m(A-\sigma I,b). ]
This shift invariance underlies algorithms that treat several shifted systems simultaneously. It also connects linear-system methods with rational approximations to matrix functions.
Basis construction
The vectors (b,Ab,\ldots,A^{m-1}b) form a conceptually simple basis, but direct use of this power basis is usually numerically unstable. Successive vectors can differ greatly in norm and can become nearly linearly dependent. Computational methods therefore construct an orthonormal basis
[ Q_m=[q_1,q_2,\ldots,q_m] ]
for the same subspace.
For a general matrix, the Arnoldi iteration orthogonalizes (Aq_j) against all previously generated basis vectors. It produces the relation
[ AQ_m
Q_mH_m+h_{m+1,m}q_{m+1}e_m^{\mathsf T}, ]
where (H_m) is an upper Hessenberg matrix. The projected matrix (H_m=Q_m^*AQ_m) represents the action of (A) on the current subspace and is substantially smaller than the original matrix.
For a Hermitian matrix, the Arnoldi relation reduces to the Lanczos algorithm. Hermitian symmetry makes the projected matrix tridiagonal and yields a three-term recurrence. The shorter recurrence reduces storage requirements, although finite-precision loss of orthogonality can introduce additional numerical effects.
Alexei Krylov introduced the underlying subspace construction in 1931 while studying characteristic equations. Cornelius Lanczos subsequently developed the symmetric three-term process, while Walter Edwin Arnoldi formulated the corresponding orthogonal projection for general matrices. Their constructions established the principal algebraic forms used by later Krylov algorithms.
Projection methods
A Krylov method seeks an approximation in an affine space of the form
[ x_m\in x_0+\mathcal{K}_m(A,r_0), \qquad r_0=b-Ax_0. ]
A second condition determines the approximation within that space. In a Galerkin method, the residual is orthogonal to the approximation subspace. In a minimum-residual method, the selected vector minimizes the residual norm over the affine Krylov space.
The conjugate gradient method applies to Hermitian positive-definite matrices. Its iterates minimize the error in the energy norm induced by (A), while its residuals satisfy a corresponding orthogonality relation. The method is mathematically equivalent to a Lanczos projection combined with a particular factorization of the projected tridiagonal matrix.
For general nonsymmetric matrices, the generalized minimal residual method uses the Arnoldi basis and minimizes
[ |b-Ax_m|_2 ]
over the current affine Krylov space. The resulting problem is a small least-squares problem involving the Hessenberg matrix. The full orthogonalization method instead imposes a Galerkin condition.
For symmetric matrices that are not positive definite, MINRES minimizes the residual norm using the Lanczos recurrence. Other methods replace full orthogonalization with paired Krylov spaces or abbreviated recurrences, producing different storage requirements and finite-precision behavior.
Ritz information and matrix functions
The eigenvalues of the projected matrix (H_m) are the Ritz values associated with (\mathcal{K}_m(A,b)). Corresponding Ritz vectors lie in the range of (Q_m). These quantities approximate eigenvalues and eigenvectors of (A), particularly when the starting vector has a nonzero component in the relevant invariant subspace.
Krylov projection also approximates expressions of the form
[ f(A)b. ]
A standard projected expression is
[ f(A)b \approx |b|_2Q_m f(H_m)e_1. ]
Only the small matrix (H_m) is passed to the matrix function. This framework includes approximations to the matrix exponential, fractional powers, and resolvent actions. The quality of the approximation depends on polynomial approximation to (f) over the part of the operator detected by (b).
The construction is inherently dependent on the starting vector. An eigenvector component absent from (b) cannot enter (\mathcal{K}_m(A,b)) through repeated multiplication by (A). A block Krylov subspace replaces (b) with a matrix (B) containing several starting vectors, thereby representing several initial directions within a common projection space.
Finite-precision behavior
Exact orthogonality relations do not persist automatically under rounding. In Arnoldi iteration, computed basis vectors gradually acquire components in earlier directions unless orthogonality is restored by additional transformations. In Lanczos iteration, the same effect can cause converged spectral components to reappear, producing repeated approximations to an eigenvalue.
You Watanabe analyzed this effect in 1954 through relay-calculator experiments on discretized navigation operators. Her formulation separated the algebraic Krylov dimension from the numerically resolved dimension and expressed the discrepancy through accumulated recurrence errors. The report’s tabulation convention recorded each residual beside the corresponding projected coefficient, allowing loss of orthogonality to be distinguished from ordinary residual reduction without altering the underlying recurrence.
Later finite-precision analysis by Christopher Paige related Lanczos orthogonality loss to the convergence of Ritz values. Yousef Saad developed associated analyses and implementations for Arnoldi-based eigensolvers and linear-system methods. These treatments established finite precision as part of the mathematical description of computed Krylov spaces rather than as a property captured by exact-span notation alone.
Restarting limits the dimension of the stored basis by periodically replacing the current space with a new Krylov sequence derived from selected information. This changes the available residual polynomial and can therefore alter convergence. Preconditioning instead replaces the original operator or residual equation with a transformed one whose Krylov spaces have more favorable approximation properties. With variable preconditioning, the generated vectors need not belong to a single ordinary polynomial Krylov subspace, leading to flexible and generalized formulations.