Generalized Minimal Residual Method

Generalized minimal residual method (GMRES) is an iterative Krylov subspace solver in numerical linear algebra and scientific computing, introduced by Yousef Saad and Martin H. Schultz in 1986, for approximating the solution of a large, sparse, general non-symmetric (or non-Hermitian) system of linear equations Ax = b where A ∈ ℝ^n×n is invertible. While the conjugate gradient method solves symmetric positive-definite systems and MINRES solves symmetric indefinite systems using short three-term recurrences from the Lanczos algorithm, the Faber–Manteuffel theorem (1984) proves that no universal short-recurrence optimal Krylov method exists for general non-symmetric matrices. GMRES overcomes this barrier by combining Arnoldi orthogonalization with progressive QR decomposition to find the iterate x_k in the affine Krylov subspace x₀ + K_k(A, r₀) that minimizes the Euclidean residual norm ‖b − Ax_k‖₂.

Because GMRES minimizes the true 2-norm of the residual over an expanding nested sequence of subspaces K₁ ⊂ K₂ ⊂ ⋯ ⊂ K_k, its residual norm ‖r_k‖₂ is guaranteed to decrease monotonically (‖r_k+1‖₂ ≤ ‖r_k‖₂) without the erratic residual spikes or division-by-zero breakdowns that can afflict biorthogonal methods such as BiCG and CGS. Consequently, preconditioned GMRES serves as the workhorse linear solver inside computational fluid dynamics, electromagnetic scattering, reactor neutronics, and Jacobian-free Newton–Krylov implementations of Newton's method.

Given an initial guess x₀ ∈ ℝⁿ with initial residual r₀ = b − Ax₀ and norm β = ‖r₀‖₂, the order-k Krylov subspace is defined as:

K_k(A, r₀) = span{r₀, Ar₀, A²r₀, …, A^k−1r₀}.

Because the monomial vectors A^jr₀ rapidly align with the dominant eigenvector of A (as in power iteration), forming a basis directly from {A^jr₀} would produce a nearly singular matrix. GMRES instead runs the Arnoldi iteration, which applies modified Gram–Schmidt orthogonalization (or Householder reflectors) to construct an orthonormal basis Q_k = [q₁, q₂, …, q_k] ∈ ℝ^n×k starting from q₁ = r₀ / β. At iteration j, the vector w = Aq_j is orthogonalized against all previous vectors q₁, …, q_j via h_ij = q_i^Tw and w ← w − ∑_i=1^j h_ijq_i, followed by setting h_j+1,j = ‖w‖₂ and q_j+1 = w / h_j+1,j.

After k steps, the Arnoldi vectors satisfy the Arnoldi relation:

AQ_k = Q_k+1 H̃_k,

where Q_k+1 ∈ ℝ^n×(k+1) has orthonormal columns and H̃_k ∈ ℝ^(k+1)×k is an upper Hessenberg matrix (having non-zero entries only on and above the main diagonal plus the first subdiagonal h_j+1,j).

Every candidate iterate in x₀ + K_k(A, r₀) can be written as x_k = x₀ + Q_ky for a low-dimensional coefficient vector y ∈ ℝ^k. Substituting the Arnoldi relation and noting that r₀ = βq₁ = Q_k+1(βe₁) yields:

‖b − Ax_k‖₂ = ‖r₀ − AQ_ky‖₂ = ‖Q_k+1(βe₁ − H̃_ky)‖₂ = ‖βe₁ − H̃_ky‖₂,

where the final equality holds because Q_k+1 has orthonormal columns and therefore preserves Euclidean norms. This reduces an n-dimensional minimization problem (where n may be 10⁶ to 10⁹) to a tiny (k + 1) × k upper Hessenberg least-squares problem (where k is typically 20 to 100).

Rather than refactoring H̃_k from scratch at every step, GMRES updates the QR decomposition of H̃_k incrementally: applying the previous k − 1 Givens plane rotations to the newly appended k-th column of H̃_k leaves only a single subdiagonal entry h_k+1,k, which is zeroed by one new Givens rotation G_k. As a remarkable byproduct, the magnitude |g_k+1| of the bottom component of the rotated right-hand side vector g = G_k⋯G₁(βe₁) equals the residual norm ‖b − Ax_k‖₂ at step k. Thus GMRES monitors convergence at every step without assembling the vector x_k = x₀ + Q_ky_k until the stopping tolerance is satisfied.

Memory Growth and Restarted GMRES(m)

Because Arnoldi orthogonalization must retain all k basis vectors q₁, …, q_k, the storage cost grows as O(nk) words and the orthogonalization work grows as O(nk²) flops after k steps. To cap memory and arithmetic overhead on large problems, Saad and Schultz introduced restarted GMRES, denoted GMRES(m): after a fixed cycle length m (commonly m = 30 to 100), if convergence has not yet occurred, the algorithm computes x_m, discards the basis Q_m, sets x₀ ← x_m, and begins a fresh m-step cycle. If A + A^T is positive-definite, GMRES(m) is guaranteed to converge for every restart parameter m ≥ 1; for strongly indefinite or non-normal matrices, small restart lengths can stall unless augmented by deflated invariant subspaces (GMRES-DR) or strong preconditioning.

Spectral vs. Non-Normal Convergence Behavior and Preconditioning

If A = V Λ V⁻¹ is diagonalizable with eigenvalues {λ₁, …, λ_n}, the residual at step k satisfies the polynomial bound:

‖r_k‖₂ / ‖r₀‖₂ ≤ κ₂(V) · min_{p∈Pk, p(0)=1} max_1≤i≤n |p(λ_i)|,

where P_k is the space of polynomials of degree at most k satisfying p(0) = 1. When A is close to normal (κ₂(V) ≈ 1) and its eigenvalues are tightly clustered away from the origin in the complex plane, a low-degree polynomial achieves tiny values across the entire spectrum, yielding rapid convergence. For highly non-normal matrices (such as convection-dominated flow operators where κ₂(V) ≫ 1), convergence is governed by the ε-pseudospectrum and the field of values of A rather than eigenvalues alone.

In practice, GMRES is almost always paired with a preconditioner M ≈ A—such as incomplete LU decomposition (ILU), algebraic multigrid (AMG), or domain decomposition—solving either the left-preconditioned system M⁻¹Ax = M⁻¹b or the right-preconditioned system AM⁻¹u = b (with x = M⁻¹u, which preserves the unpreconditioned residual norm ‖b − Ax_k‖₂). When the preconditioner itself varies from step to step (for example, an inner iterative solver), Saad's Flexible GMRES (FGMRES) stores the preconditioned vectors z_j = M_j⁻¹q_j separately to maintain convergence.

What is a "lucky breakdown" in GMRES?

If the subdiagonal entry h_k+1,k = 0 at step k < n, the Krylov subspace K_k(A, r₀) has become an invariant subspace of A (AK_k ⊆ K_k). Far from being a failure, this implies that the exact solution x^* = A⁻¹b lies inside x₀ + K_k(A, r₀), so GMRES terminates at step k with zero residual ‖r_k‖₂ = 0.

Why is right preconditioning often preferred over left preconditioning in GMRES?

Left preconditioning minimizes the preconditioned residual ‖M⁻¹(b − Ax_k)‖₂, which can distort stopping criteria if M is poorly conditioned. Right preconditioning solves AM⁻¹u = b and minimizes the true physical residual ‖b − Ax_k‖₂ directly.

1. Saad, Y., & Schultz, M. H. (1986). GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM Journal on Scientific and Statistical Computing, 7(3), 856–869. https://doi.org/10.1137/0907058
2. Faber, V., & Manteuffel, T. (1984). Necessary and sufficient conditions for the existence of a conjugate gradient method. SIAM Journal on Numerical Analysis, 21(2), 352–362. https://doi.org/10.1137/0721026
3. Saad, Y. (1993). A flexible inner-outer preconditioned GMRES algorithm. SIAM Journal on Scientific Computing, 14(2), 461–469. https://doi.org/10.1137/0914028
4. Greenbaum, A., Pták, V., & Strakoš, Z. (1996). Any nonincreasing convergence curve is possible for GMRES. SIAM Journal on Matrix Analysis and Applications, 17(3), 465–469. https://doi.org/10.1137/S0895479894275030