1 Introduction to Krylov Subspace Methods
1.1 The residual and minimal residual idea
For a linear system \(Ax=b\), the residual is \(r_k=b-Ax_k\), measuring how far an approximate solution \(x_k\) is from satisfying the equation. Krylov subspace methods build approximations so that the residual becomes small. The defining idea behind GMRES is to choose, at each iteration, the approximation from a growing subspace that minimizes a norm of the residual. This “minimal residual” viewpoint gives GMRES its robustness for nonsymmetric problems and its close connection to least-squares projection methods.
1.2 Krylov subspaces and why they matter
Krylov subspaces are generated from the initial residual \(r_0=b-Ax_0\). The \(k\)-th Krylov subspace is \[ \mathcal{K}_k(A,r_0)=\text{span}\{r_0, Ar_0, A^2r_0,\dots,A^{k-1}r_0\}. \] Any method that searches within \(\mathcal{K}_k\) produces iterates expressible via low-degree polynomials in \(A\) applied to \(r_0\). This structure lets the algorithm exploit matrix-vector products without requiring matrix factorization, and it provides a path for deriving practical projection and least-squares formulations.
1.3 Comparison with related iterative methods
GMRES belongs to a family of Krylov methods that includes variants tailored to symmetric or positive definite systems. Conjugate Gradient (CG) uses a short-recurrence framework and assumes symmetry and positive definiteness, whereas GMRES does not require these properties. Methods such as BiCG-type algorithms can handle nonsymmetric matrices but may produce less directly interpretable residual minimization. In practice, GMRES is often favored for general square nonsymmetric systems when reliable residual reduction is important, though it may cost more per iteration than short-recurrence alternatives.
2 Mathematical Formulation of GMRES
2.1 Problem setting: solving Ax = b
Consider a square linear system \(Ax=b\) with \(A\in\mathbb{C}^{n\times n}\) or \(\mathbb{R}^{n\times n}\). Given an initial guess \(x_0\), define \(r_0=b-Ax_0\). GMRES seeks iterates \(x_k\) of the form \[ x_k = x_0 + y_k, \] where \(y_k\in \mathcal{K}_k(A,r_0)\). The method focuses on how to choose \(y_k\) to make the residual small.
2.2 Residual norm minimization over Krylov subspaces
| GMRES chooses \(x_k\) to minimize \(\|r_k\|_2\) over all candidates in the Krylov subspace: |
|---|
\[
| x_k = \arg\min_{x_0 + y \,:\, y\in \mathcal{K}_k(A,r_0)} \|b-A(x_0+y)\|_2. |
|---|
\] Because the minimization occurs in a space of growing dimension, the residual norm is nonincreasing (in exact arithmetic) up to the dimension limit, making the “minimal residual” objective central to the method’s behavior.
2.3 Orthogonal projection viewpoint
The residual minimization can be recast as an orthogonal projection process. In the Arnoldi-based implementation, GMRES builds an orthonormal basis \(V_{k+1}\) and a related basis \(V_k\) spanning \(\mathcal{K}_k\). The computed iterate corresponds to choosing coefficients \(y_k\) so that the new residual is orthogonal to the expanded subspace generated by the basis. This projection interpretation explains why GMRES reduces the original problem to a smaller least-squares problem: the residual norm in the full space becomes the norm of a structured vector in a low-dimensional setting.
3 Arnoldi Process Underlying GMRES
3.1 Building an orthonormal Krylov basis
| GMRES typically relies on the Arnoldi process to produce an orthonormal basis for the Krylov subspace. Starting with \(v_1=r_0/\|r_0\|\), each new vector \(v_{j+1}\) is generated from \(Av_j\) and then orthogonalized against the existing basis vectors. In exact arithmetic, this yields |
|---|
\[ V_k=[v_1,\dots,v_k], \quad V_k^\ast V_k = I, \] with the columns spanning \(\mathcal{K}_k(A,r_0)\).
3.2 The upper Hessenberg relation
A key structural byproduct of Arnoldi is the relation \[ AV_k = V_{k+1}\,\bar{H}_k, \] where \(\bar{H}_k\) is an \((k+1)\times k\) upper Hessenberg matrix. This matrix encodes the projection coefficients from the orthogonalization step. The Hessenberg form is crucial: it makes the reduced problem inexpensive to solve compared with operating directly in \(\mathbb{C}^n\).
3.3 Least-squares subproblem structure
Using the Arnoldi relation, the GMRES residual minimization becomes a least-squares problem in the small space of dimension \(k\). The coefficients \(y_k\) are chosen so that the projected residual has minimal 2-norm, often expressed as \[
| \min_{z\in\mathbb{C}^k}\left\|\beta e_1 - \bar{H}_k z\right\|_2, |
|---|
\]
| where \(\beta=\|r_0\|\) and \(e_1\) is the first coordinate basis vector. Once \(z\) is found, the GMRES iterate is recovered by \(x_k=x_0+V_k z\). |
|---|
3.4 Numerical aspects of orthogonalization
In floating-point arithmetic, orthogonality can degrade as \(k\) grows. The choice of orthogonalization method (e.g., classical vs. modified Gram–Schmidt, and whether reorthogonalization is used) affects stability and attainable residual reduction. Since GMRES depends on accurate basis construction for its minimization property, numerical orthogonality directly influences whether the residual norm behaves as predicted and whether stagnation occurs earlier than expected.
4 Algorithmic Description
4.1 Step-by-step GMRES iteration
A typical GMRES iteration proceeds as follows:
| 1. Initialize: choose \(x_0\), compute \(r_0=b-Ax_0\), set \(\beta=\|r_0\|\). If \(\beta=0\), stop. |
|---|
- Normalize: set \(v_1=r_0/\beta\).
- For \(j=1,\dots,k\):
a. Compute \(w=A v_j\). b. Orthogonalize \(w\) against \(v_1,\dots,v_j\) to obtain the next basis vector component and the Hessenberg coefficients. c. Normalize to get \(v_{j+1}\) (if possible).
- Solve the reduced least-squares problem to obtain coefficients \(z\).
- Form \(x_k=x_0+V_k z\).
- Update the residual norm from the reduced system solution.
In efficient implementations, the least-squares solve is performed incrementally using QR-based updates rather than recomputing from scratch each time.
4.2 Implementation details: stopping criteria
Stopping criteria usually combine a relative residual test and absolute tolerance, such as \[
| \frac{\|b-Ax_k\|_2}{\|b\|_2} \le \text{tol} \quad \text{or} \quad \|b-Ax_k\|_2 \le \text{abs\_tol}. |
|---|
\]
| Because the residual computed via projections can differ slightly from the true residual in finite precision, robust codes may periodically recompute \(\|b-Ax_k\|_2\) explicitly or use safeguards when stagnation is suspected. |
|---|
4.3 Handling breakdown and stagnation
Two practical failure modes are:
- Breakdown: a new Arnoldi vector cannot be formed due to near-zero norm, which may occur if the Krylov subspace becomes invariant or if round-off causes loss of progress. In exact arithmetic, an early breakdown can indicate exact solvability within the current subspace.
- Stagnation: the residual norm stops decreasing meaningfully even though the method continues iterating. This can result from loss of orthogonality, poor conditioning, or matrix structures that limit the effectiveness of polynomial residual reduction.
Common responses include restarting, switching orthogonalization strategies, or applying preconditioning.
5 Computational Complexity and Storage
5.1 Cost per iteration
Each iteration requires one matrix-vector product \(A v_j\), which usually dominates when \(A\) is large and sparse. Additional cost comes from orthogonalization against the existing basis vectors, leading to a growing number of inner products and vector updates. The total arithmetic cost up to iteration \(k\) reflects both the per-iteration matvec cost and the cumulative orthogonalization work.
5.2 Memory requirements
GMRES stores the growing basis vectors \(V_k\) and related Hessenberg/QR data structures. This can be memory-intensive for large \(k\) because storing \(k\) vectors of length \(n\) scales like \(O(nk)\). The storage burden is one reason restarted GMRES is widely used in practice.
5.3 Effects of orthogonalization strategy
Better orthogonalization improves stability but increases overhead due to more inner products or reorthogonalization passes. Conversely, cheaper orthogonalization may reduce runtime but risk loss of the minimization property and earlier stagnation. The “best” strategy depends on precision, matrix size, and whether the environment favors memory bandwidth, latency, or floating-point throughput.
6 Convergence Behavior
6.1 Factors influencing convergence
Convergence rates depend on spectral properties of \(A\), how eigenvalues cluster, and how the initial residual aligns with invariant subspaces. The condition number of eigenvector matrices (for nonnormal matrices) can also matter: even when eigenvalues look favorable, nonnormality can make GMRES behave unpredictably. Preconditioning can reshape the effective spectrum seen by GMRES and thus improve residual reduction.
6.2 Residual polynomial interpretation
GMRES residuals can be expressed in terms of polynomials. With \(x_k=x_0+V_k z\), the residual norm satisfies \[ r_k = p_k(A)\,r_0 \]
| for some polynomial \(p_k\) of degree \(k\) with \(p_k(0)=1\), and GMRES selects the polynomial minimizing \(\|p_k(A)r_0\|_2\) over the Krylov subspace. This framework links convergence to how well low-degree polynomials can damp components associated with the spectrum (and, for nonnormal systems, other pseudospectral features). |
|---|
6.3 Practical observations vs. theory
Theoretical convergence bounds often rely on idealized assumptions or may be conservative. In real computations, observed behavior frequently reflects:
- finite precision effects that limit basis quality,
- imperfect preconditioners,
- model mismatch in matrix formation,
- and the presence of nearly invariant subspaces.
As a result, practitioners often tune restarting, tolerances, and preconditioner choices based on empirical performance while using theory to guide expectations.
7 Restarted GMRES (GMRES(m))
7.1 Motivation for restarting
Because memory and cost grow with the number of Krylov vectors, restarted GMRES limits the subspace dimension. After \(m\) inner iterations, the method discards the basis and begins again using the current approximate solution as the new initial guess. This caps storage at \(m\) vectors and keeps per-cycle cost controlled.
7.2 Trade-offs between iteration count and subspace size
Using a smaller restart parameter \(m\) reduces memory and orthogonalization overhead but can slow convergence because the method loses information accumulated in earlier cycles. Larger \(m\) increases work per cycle but can yield stronger residual reduction per cycle. The optimal \(m\) depends on problem size, matrix nonnormality, available memory, and preconditioning quality.
7.3 Choosing the restart parameter m
Common practice selects \(m\) from a moderate range (often problem-dependent heuristics such as 20–200 in many legacy setups). Criteria include:
- acceptable memory usage for storing \(m\) basis vectors,
- runtime balance between matvecs and orthogonalization,
- and observed convergence trends per cycle.
Monitoring residual decrease per restart cycle helps decide whether to increase or decrease \(m\).
8 Preconditioning for GMRES
8.1 Left, right, and split preconditioning
Preconditioning transforms the system into one with more favorable properties for iterative solution.
- Left preconditioning: solve \(M^{-1}Ax = M^{-1}b\). The GMRES residual norm is influenced by the left transformation.
- Right preconditioning: solve \(AM^{-1}y = b\), with \(x=M^{-1}y\). The residual is computed with respect to the original equation, which can be advantageous for stopping criteria.
- Split preconditioning: use different transformations on the left and right when appropriate, balancing spectrum and norm effects.
8.2 Designing effective preconditioners
An effective preconditioner approximates \(A\) in a way that clusters eigenvalues (or reduces pseudospectral spread) while being cheap to apply each iteration. For nonsymmetric problems, quality often depends on how the preconditioner handles nonnormality. Typical strategies include incomplete factorization approaches, algebraic multigrid adapted for nonsymmetric operators, and specialized approximate solvers.
8.3 Interaction of preconditioning and Krylov subspaces
Preconditioning changes the Krylov subspaces because the basis is built using matrix-vector products with the preconditioned operator. Consequently, the residual minimization polynomials act on the transformed system rather than the original one. This can significantly improve or occasionally worsen convergence depending on how accurately the preconditioner captures the operator’s action.
9 Variants and Enhancements
9.1 Flexible GMRES (FGMRES)
FGMRES allows the preconditioner to vary between iterations. This is useful when the application of the preconditioner is itself iterative, stochastic, nonlinear, or otherwise not fixed. Standard GMRES assumes a consistent linear operator in the Krylov construction; FGMRES relaxes that assumption by building the approximation space using the preconditioned directions while maintaining a proper residual minimization framework.
9.2 Block GMRES concepts
Block methods handle multiple right-hand sides simultaneously, constructing Krylov subspaces with matrix-valued bases. Block GMRES can be advantageous when several linear systems \(Ax=b^{(i)}\) share the same coefficient matrix, as it can exploit shared structure and improve robustness against rank deficiencies in the residuals. The computational pattern changes from vector operations to small-matrix operations.
9.3 Deflation and recycling subspaces
Deflation methods aim to remove the influence of troublesome modes that slow convergence. Recycling strategies preserve parts of the Krylov basis from previous solves or cycles to accelerate later iterations, particularly when solving a sequence of related systems, such as time-stepping problems. Deflation can be incorporated into restarted GMRES by augmenting the search space with approximate eigenvector information or selected directions extracted from prior runs.
9.4 Variant stopping and accuracy strategies
Enhancements may use inexact computations—for example, approximate preconditioner solves or relaxed orthogonalization thresholds—and adjust stopping rules accordingly. Some implementations use residual estimates from the reduced least-squares system, while others periodically recompute the true residual for reliability. Accuracy strategies balance performance against the risk of false convergence or excessive iterations.
10 Practical Guidance and Best Practices
10.1 Selecting initial guesses
A good initial guess can reduce the number of Krylov iterations needed. In problems where solutions are correlated across parameters or time steps, a previous solution can serve as \(x_0\). For multiple right-hand sides, scaling and approximate solves can also yield better starting points. Since GMRES builds the subspace from the initial residual, alignment between \(r_0\) and difficult modes strongly affects performance.
10.2 Monitoring residuals and detecting issues
Residual monitoring should include:
- checking whether the residual decreases steadily over iterations or cycles,
- observing if residual reduction stalls after orthogonality loss or restart,
- and verifying that the reported residual matches the true residual when troubleshooting.
When progress halts, interventions may include increasing \(m\), applying a stronger preconditioner, changing orthogonalization, or switching to a flexible or recycled strategy.
10.3 Numerical stability considerations
Stability hinges on maintaining an accurate orthonormal basis. Practical steps include using modified Gram–Schmidt with reorthogonalization, employing QR updates for the least-squares solve, and avoiding unnecessary loss of information in reductions. For very ill-conditioned problems, scaling of the system, careful choice of norms, and periodic recomputation of residuals may improve reliability.
11 Applications and Use Cases in Applied Mathematics
11.1 Nonsymmetric linear systems in modeling
Many models lead to nonsymmetric operators, for example when advection, damping, or non-self-adjoint effects are present. GMRES is suited to such systems because it does not rely on symmetry-based recurrences. Its projection-based residual minimization provides predictable behavior for a wide class of operators encountered in modeling and simulation.
11.2 Discretized PDEs and linear solvers
Discretization of partial differential equations often produces large sparse linear systems. For convection-dominated or coupled PDEs, nonsymmetry is common. GMRES is a standard choice when direct solvers are too expensive and when iterative refinement via preconditioning can make the Krylov process efficient.
11.3 Systems arising in scientific computing
GMRES also appears in areas such as computational fluid dynamics, electromagnetic simulations, and inverse problems, where repeated solution of linear systems is routine. In many workflows, combining GMRES with tailored preconditioners, restart parameters, and recycling strategies can turn a difficult sequence of solves into a practical computational pipeline.
12 Summary and Further Reading
12.1 Key takeaways
GMRES is an iterative method that searches for an approximate solution within a Krylov subspace and chooses the iterate that minimizes the residual norm. It relies on the Arnoldi process and reduces each iteration’s work to a small least-squares problem. While costs and storage grow with subspace dimension, restarted GMRES and preconditioning make the method practical for large nonsymmetric systems.
12.2 Suggested resources for deeper study
For deeper study, readers often consult texts covering Krylov subspace methods, numerical linear algebra, and iterative solver design. Key topics to pursue include Arnoldi decomposition, least-squares formulations, nonnormal convergence phenomena, and the theory and practice of preconditioning and restarting. Additional value comes from sources that emphasize implementation details, stability, and performance tuning on sparse systems.