Krylov Subspace: Definition, Applications, and Methods

• Krylov subspaces are subspaces generated by successive applications of a matrix (or linear operator) on a starting vector, useful for approximating solutions in various computational contexts.
• Primary methods such as Arnoldi and Lanczos are used to construct Krylov spaces, which help approximate solutions of large-scale linear systems, matrix functions, and optimization problems without explicitly forming high-dimensional matrices.
• Key techniques like preconditioning, rational extensions, and recycling enhance efficiency by reducing computational complexity, increasing accuracy, and enabling the reuse of prior information.

A Krylov subspace is the polynomially generated subspace associated with a matrix or linear operator and a starting vector. For $A\in\mathbb C^{n\times n}$ and $v\neq 0$, the order-$m$ Krylov subspace is

$$\mathcal K_m(A,v)=\operatorname{span}\{v,Av,A^2v,\ldots,A^{m-1}v\} =\{p(A)v:\deg p<m\}.$$

Its dimension is at most $\min(m,n)$ and may be smaller when the generated vectors become linearly dependent. Krylov methods exploit this low-dimensional structure to approximate solutions of linear systems, matrix functions, eigenpairs, nonlinear optimization steps, transfer operators, graph kernels, and quantum dynamics without explicitly forming matrix inverses, dense matrix functions, or full Hessians. The central computational pattern is the reduction of a large problem to a projected problem involving a small Hessenberg, tridiagonal, or structured matrix (1811.09025).

1. Definition, polynomial structure, and invariant subspaces

Given a linear system

$Ax=b$

and an initial approximation $x_0$, the initial residual is

$r_0=b-Ax_0.$

Krylov iterative methods search for approximations in the affine space

$x_m\in x_0+\mathcal K_m(A,r_0).$

Every correction in this space has the form

$z_m=p_{m-1}(A)r_0,$

where $v\neq 0$0 is a polynomial of degree at most $v\neq 0$1. The associated residual is

$v\neq 0$2

The subspaces are nested:

$v\neq 0$3

Although $v\neq 0$4 is generated by $v\neq 0$5 vectors, its actual dimension may be smaller. If

$v\neq 0$6

the sequence terminates. In exact arithmetic, this breakdown can indicate that an invariant subspace has been reached and, for some problems, that an exact solution is contained in the generated space (1811.09025).

The stabilized dimension is the degree of the minimal polynomial of $v\neq 0$7 relative to $v\neq 0$8. Equivalently, it is the smallest dimension of a subspace that contains $v\neq 0$9 and is invariant under $m$0. For a symmetric matrix and vector $m$1, this maximal space is

$m$2

where $m$3 is the smallest integer satisfying

$m$4

A truncated Krylov space contains the starting vector exactly but is generally only approximately invariant under the matrix (Jiang et al., 2024).

The polynomial interpretation is fundamental. Krylov methods do not approximate arbitrary vectors independently; they select a polynomial filter applied to the starting vector. This explains their sensitivity to eigenvalue distributions, nonnormality, conditioning, spectral clustering, and the starting vector. For matrix functions, the same principle gives

$m$5

where $m$6 is a Krylov basis and $m$7 is a small projected matrix (Frommer et al., 2022).

2. Arnoldi, Lanczos, and projected representations

For a general square matrix, Arnoldi orthogonalization constructs an orthonormal basis

$m$8

for $m$9. Starting with

$$\mathcal K_m(A,v)=\operatorname{span}\{v,Av,A^2v,\ldots,A^{m-1}v\} =\{p(A)v:\deg p<m\}.$$0

each iteration computes $$\mathcal K_m(A,v)=\operatorname{span}\{v,Av,A^2v,\ldots,A^{m-1}v\} =\{p(A)v:\deg p<m\}.$$1, orthogonalizes it against the previous basis vectors, and normalizes the result. The resulting Arnoldi relation is

$$\mathcal K_m(A,v)=\operatorname{span}\{v,Av,A^2v,\ldots,A^{m-1}v\} =\{p(A)v:\deg p<m\}.$$2

where $$\mathcal K_m(A,v)=\operatorname{span}\{v,Av,A^2v,\ldots,A^{m-1}v\} =\{p(A)v:\deg p<m\}.$$3 is upper Hessenberg. The Hessenberg structure expresses the inclusion

$$\mathcal K_m(A,v)=\operatorname{span}\{v,Av,A^2v,\ldots,A^{m-1}v\} =\{p(A)v:\deg p<m\}.$$4

For GMRES, an approximation is written as

$$\mathcal K_m(A,v)=\operatorname{span}\{v,Av,A^2v,\ldots,A^{m-1}v\} =\{p(A)v:\deg p<m\}.$$5

Since $$\mathcal K_m(A,v)=\operatorname{span}\{v,Av,A^2v,\ldots,A^{m-1}v\} =\{p(A)v:\deg p<m\}.$$6 with $$\mathcal K_m(A,v)=\operatorname{span}\{v,Av,A^2v,\ldots,A^{m-1}v\} =\{p(A)v:\deg p<m\}.$$7, residual minimization reduces to

$$\mathcal K_m(A,v)=\operatorname{span}\{v,Av,A^2v,\ldots,A^{m-1}v\} =\{p(A)v:\deg p<m\}.$$8

Thus a large least-squares problem in the original dimension becomes a small least-squares problem of dimension $$\mathcal K_m(A,v)=\operatorname{span}\{v,Av,A^2v,\ldots,A^{m-1}v\} =\{p(A)v:\deg p<m\}.$$9. GMRES minimizes the Euclidean residual norm over the affine Krylov space and satisfies

$\min(m,n)$0

For symmetric matrices, Arnoldi reduces to Lanczos. The projected matrix becomes tridiagonal, and the method admits a short three-term recurrence. Conjugate gradients applies this construction to symmetric positive-definite matrices and computes

$\min(m,n)$1

where

$\min(m,n)$2

The residuals are mutually orthogonal in exact arithmetic, while the search directions are $\min(m,n)$3-conjugate. Under the standard SPD assumptions, the classical estimate is

$\min(m,n)$4

Arnoldi requires increasing orthogonalization and storage. After $\min(m,n)$5 steps, the orthogonalization cost is approximately $\min(m,n)$6 and the basis requires $\min(m,n)$7 storage. Lanczos reduces these costs through short recurrences but requires symmetry or a corresponding structured inner product. Loss of orthogonality can occur in finite precision; modified Gram–Schmidt, reorthogonalization, Householder transformations, or structure-preserving reorthogonalization may be required (1811.09025).

3. Extensions: rational, extended, and structured Krylov spaces

Polynomial Krylov spaces use nonnegative powers of $\min(m,n)$8. Rational Krylov spaces replace monomials by rational functions with prescribed poles:

$\min(m,n)$9

A finite pole $Ax=b$0 corresponds to a shifted inverse

$Ax=b$1

whereas an infinite pole corresponds to multiplication by $Ax=b$2. Every finite pole therefore requires solving a shifted linear system (Buggenhout et al., 2018).

An extended Krylov space includes both positive and negative powers:

$Ax=b$3

It is a rational Krylov space with poles at $Ax=b$4 and $Ax=b$5. The negative-power component can improve approximation for functions with important behavior near the origin or across widely separated spectral scales, but requires linear solves with $Ax=b$6 or shifted variants; the inverse need not be formed explicitly.

For non-Hermitian matrices, rational Krylov methods can use distinct right and left spaces,

$Ax=b$7

with biorthogonal bases satisfying

$Ax=b$8

The oblique projection is

$Ax=b$9

Rational Arnoldi naturally produces a matrix pencil

$x_0$0

and the reduced eigenvalue problem is

$x_0$1

Biorthogonal extended Krylov methods can choose complementary Hessenberg and inverse-Hessenberg representations so that the projected operator is represented by a tridiagonal pencil:

$x_0$2

The pole sequence appears in ratios of subdiagonal and superdiagonal entries. This yields compact storage and short recurrences analogous to Lanczos, although oblique projection is more vulnerable to loss of biorthogonality (Buggenhout et al., 2018).

For Hamiltonian matrices, ordinary Euclidean projection generally fails to preserve Hamiltonian structure. The Hamiltonian Extended Krylov Subspace method instead constructs a $x_0$3-orthogonal basis $x_0$4 satisfying

$x_0$5

The projected matrix

$x_0$6

is Hamiltonian. HEKS combines positive and negative powers with structure-preserving projection and uses short recurrences involving at most five previously generated basis vectors. It is intended for approximating $x_0$7, particularly when inverse powers are beneficial, as in matrix-sign computations (Benner et al., 2022).

4. Preconditioning, augmentation, and recycling

Preconditioning changes the operator whose Krylov space is generated. Left preconditioning uses

$x_0$8

whereas right preconditioning uses

$x_0$9

Right-preconditioned GMRES minimizes the residual in the original equation, while the Krylov space is generated by $r_0=b-Ax_0.$0. Incomplete LU and incomplete Cholesky are examples of sparse preconditioners.

Augmented Krylov methods enlarge the search space with a fixed subspace $r_0=b-Ax_0.$1:

$r_0=b-Ax_0.$2

The augmentation space may contain approximate eigenvectors, previous solution components, or directions generated by earlier solves. In projected recycling methods, the Krylov process is generated using an operator such as

$r_0=b-Ax_0.$3

where $r_0=b-Ax_0.$4 projects onto the image of the augmentation space. In unprojected methods, the Krylov space remains

$r_0=b-Ax_0.$5

and the augmentation is incorporated through the residual-minimization equations. This separation can be useful when the augmentation space is inaccurate, because it prevents erroneous recycled information from determining the entire Krylov recurrence (Soodhalter, 2021).

Subspace recycling is particularly effective for sequences of related linear systems. For shifted systems

$r_0=b-Ax_0.$6

the polynomial Krylov spaces satisfy

$r_0=b-Ax_0.$7

when the right-hand side is fixed or suitably compatible. However, this shift invariance does not imply that the complete recycled or augmented solve is shift-independent. Distinct right-hand sides produce distinct residual directions, and

$r_0=b-Ax_0.$8

makes the effect of the recycled space shift-dependent. One scheme therefore uses a separate deflation space for each shift; another retains one recycled space, computes accurate corrections for a base system, and constructs approximate corrections for the remaining shifts recursively (Soodhalter et al., 2013).

For evolving finite-element meshes, algebraic vectors from one system may have different dimensions, node meanings, and sparsity structures in the next system. A recycled vector can instead be interpreted as the coefficient vector of a finite-element function, transferred to the new mesh, and re-expanded in the new basis. Approximate invariant subspaces associated with small eigenvalues are especially useful because they are often smooth, low-frequency modes and therefore more transferable. A warm-start Krylov–Schur process can improve a transferred subspace when geometric mapping is inaccurate (Bolten et al., 2020).

Proper truncation is essential. Retaining all previous Krylov vectors increases storage, orthogonalization, and reduced-system costs. Goal-oriented proper orthogonal decomposition can compress the accumulated vectors using a metric $r_0=b-Ax_0.$9, snapshot weights, and an energy criterion. With $x_m\in x_0+\mathcal K_m(A,r_0).$0, the resulting basis is adapted to the energy norm; with $x_m\in x_0+\mathcal K_m(A,r_0).$1, it is adapted to a quantity of interest. This approach combines POD compression, augmented CG, and a direct/iterative reduced-space strategy (Carlberg et al., 2015).

5. Matrix functions, numerical ranges, and applications

Krylov methods approximate matrix-function actions without explicitly forming $x_m\in x_0+\mathcal K_m(A,r_0).$2. The standard Arnoldi approximation is

$x_m\in x_0+\mathcal K_m(A,r_0).$3

This applies to functions represented through polynomial or contour approximations, including exponentials, fractional powers, inverse functions, logarithms, square roots, and sign functions.

For Laplace-transform functions,

$x_m\in x_0+\mathcal K_m(A,r_0).$4

the Arnoldi error can itself be represented as a Laplace-transform action on the next Arnoldi vector:

$x_m\in x_0+\mathcal K_m(A,r_0).$5

where

$x_m\in x_0+\mathcal K_m(A,r_0).$6

This identity supports a fixed-size restart method using quadrature, spline interpolation, and matrix-exponential actions on small Hessenberg matrices (Frommer et al., 2022).

For sequences of matrix-function applications, recycled FOM augments a newly generated Krylov space with a retained subspace:

$x_m\in x_0+\mathcal K_m(A,r_0).$7

The retained vectors commonly approximate eigenvectors associated with eigenvalues near singularities or branch points of $x_m\in x_0+\mathcal K_m(A,r_0).$8. A closed-form augmented FOM approximation evaluates $x_m\in x_0+\mathcal K_m(A,r_0).$9 on the projected augmented matrix and avoids quadrature. Randomized sketching can reduce the cost of full orthogonalization by replacing large inner products with approximate inner products in a lower-dimensional sketch. The resulting sketched-and-recycled FOM method is useful when orthogonalization, rather than matrix-vector multiplication, is the computational bottleneck (Burke et al., 2022, Burke et al., 2023).

The numerical range of a matrix is

$z_m=p_{m-1}(A)r_0,$0

If $z_m=p_{m-1}(A)r_0,$1 is an orthonormal Krylov basis, then

$z_m=p_{m-1}(A)r_0,$2

For normal matrices, $z_m=p_{m-1}(A)r_0,$3. With a random starting vector, the numerical range of an Arnoldi compression approximates the full numerical range with gap-independent estimates of order

$z_m=p_{m-1}(A)r_0,$4

for general normal matrices and

$z_m=p_{m-1}(A)r_0,$5

when the spectrum lies on the unit circle. Corresponding lower bounds of order $z_m=p_{m-1}(A)r_0,$6 and $z_m=p_{m-1}(A)r_0,$7 show that these dependencies are essentially optimal in general. For nonnormal matrices, the quality additionally depends on eigenbasis conditioning and on the distinction between the spectral convex hull and the numerical range (Chen et al., 2024).

Krylov spaces also support graph kernel computation. If $z_m=p_{m-1}(A)r_0,$8 is a sparse graph Laplacian and

$z_m=p_{m-1}(A)r_0,$9

block Krylov methods approximate $v\neq 0$00 without forming the full kernel matrix. Classical block Lanczos, global block Lanczos, sequential Lanczos, Chebyshev, and squared Chebyshev methods offer different tradeoffs among iteration count, memory, symmetry, and positive definiteness. Classical block Lanczos preserves positive definiteness of the collocation matrix under the stated assumptions; squared Chebyshev enforces positive semidefiniteness by squaring a polynomial approximation (Erb, 2023).

6. Optimization, dynamical systems, and quantum Krylov spaces

In second-order optimization, a Hessian-gradient Krylov space is

$v\neq 0$01

where $v\neq 0$02 and $v\neq 0$03. Krylov Cubic Regularized Newton minimizes the cubic Newton model within this space. Lanczos produces a tridiagonal projected Hessian using one gradient evaluation and Hessian-vector products, avoiding explicit Hessian formation. For convex objectives with Lipschitz Hessian, the method has convergence rate

$v\neq 0$04

independent of the ambient dimension in its iteration bound. If the Hessian-gradient pair generates an invariant Krylov space of dimension at most $v\neq 0$05, the truncation term vanishes and the full cubic-regularized Newton rate $v\neq 0$06 is recovered (Jiang et al., 2024).

For nonlinear stochastic dynamical systems, Krylov methods can be applied to a Perron–Frobenius operator acting on kernel mean embeddings. The operator propagates probability measures in an RKHS, while finite-data inner products are computed from kernel Gram matrices. Ordinary Arnoldi is appropriate when the transfer operator is bounded. When it may be unbounded, shift-invert Arnoldi constructs a Krylov space for

$v\neq 0$07

a bounded resolvent, and reconstructs

$v\neq 0$08

Maximum mean discrepancy provides a prediction and anomaly-detection metric in the embedded distribution space (Hashimoto et al., 2019).

In network analysis, the columns

$v\neq 0$09

are retained as normalized power-iteration vectors, and their rows are interpreted as node-specific Krylov subspace trajectories. These trajectories encode how node weights propagate across successive path lengths. Unlike eigenvector centrality, which uses primarily the final converged vector, trajectory-based methods retain transient information for node similarity, clustering, perturbation analysis, and exploratory node-importance measures (Frost, 2024).

Quantum Krylov methods construct

$v\neq 0$10

or generalized spaces generated by functions $v\neq 0$11. The projected Hamiltonian and overlap matrices satisfy the generalized eigenvalue problem

$v\neq 0$12

Finite-shot measurement makes $v\neq 0$13 nearly singular as the basis dimension grows. Gaussian-power bases,

$v\neq 0$14

suppress distant spectral components and reduce measurement instability. Quantum Krylov methods are used for ground-state and excited-state energies, molecular properties, reduced density matrices, and quantum simulation (Zhang et al., 2023, Oumarou et al., 9 Jan 2025).

In quantum dynamics, the Krylov space may be generated by a Hamiltonian acting on a state, by the Liouvillian

$v\neq 0$15

acting on an operator, or by a Lindbladian acting on a density matrix. Lanczos orthogonalization yields a tridiagonal representation. For an operator Krylov basis $v\neq 0$16,

$v\neq 0$17

The evolved operator is represented by amplitudes on a one-dimensional Krylov chain. The Krylov complexity is

$v\neq 0$18

while the variance and entropy characterize the spreading of the Krylov wavepacket. Linear Lanczos-coefficient growth can imply exponential Krylov complexity, but it is more directly associated with rapid operator growth or scrambling than with quantum chaos itself. Initial-operator dependence, finite-size saturation, inner-product choice, non-Hermiticity, and ultraviolet or infrared regularization must be considered when interpreting these quantities (Nandy et al., 2024).

Across these applications, the defining tradeoff is consistent: increasing the Krylov dimension improves representational capacity and spectral resolution but increases storage, orthogonalization, matrix-vector products, shifted solves, or measurement cost. Recycling, rational enrichment, structured projections, restarting, truncation, POD compression, and randomized sketching modify this tradeoff without changing the central principle: approximate a high-dimensional operation through information generated by repeated action of an operator on a carefully chosen starting object.