GMRES

Every computer simulation — weather forecast, crash test, fluid dynamics — eventually reduces to solving a linear system $Ax = b$. When the matrix $A$ has a million rows, a straightforward Gaussian-elimination approach would need $10^{18}$ arithmetic operations. That is roughly one million years on today's fastest chip.

GMRES (Generalized Minimal RESidual), introduced by Yousef Saad and Martin Schultz in 1986, sidesteps this wall. Instead of factoring $A$ it builds, step by step, a small Krylov subspace

$$\mathcal{K}_k(A,r_0) = \mathrm{span}\{r_0,\, Ar_0,\, A^2 r_0,\, \dots,\, A^{k-1}r_0\}$$

where $r_0 = b - Ax_0$ is the initial residual. At each iteration $k$ it finds the vector $x_k$ inside that subspace that minimises $\|b - Ax_k\|_2$ — the residual norm. In many practical problems only a few hundred iterations are needed, even when $n$ is in the millions.

The key insight: you never need to store or factor $A$. You only need to multiply $A$ by a vector — a single sparse operation that costs $\mathcal{O}(n)$ when $A$ has $\mathcal{O}(n)$ non-zeros. That transforms an astronomically hard direct problem into something a laptop can track.

The demo below runs GMRES on a randomly generated $6 \times 6$ non-symmetric system $Ax = b$. Watch how the residual $\|r_k\| = \|b - Ax_k\|$ drops with each new iteration.

<!-- {{c_html_intro}} -->
<div class="controls">
  <button id="btn-step" type="button">{{btn_step}}</button>
  <button id="btn-run" type="button">{{btn_run}}</button>
  <button id="btn-new" type="button" class="ghost">{{btn_new}}</button>
</div>
<div id="info" class="info"></div>
<canvas id="chart" width="480" height="200" aria-label="{{chart_aria}}"></canvas>
<div id="table-wrap">
  <table id="iter-table">
    <thead><tr><th>{{th_iter}}</th><th>{{th_residual}}</th><th>{{th_bar}}</th></tr></thead>
    <tbody id="tbody"></tbody>
  </table>
</div>

/* {{c_css_intro}} */
* { box-sizing: border-box; }
body { font-family: system-ui, sans-serif; color: #222; margin: 0; padding: .5rem; }
.controls { display: flex; gap: .5rem; flex-wrap: wrap; margin-bottom: .6rem; }
button { font: 600 14px system-ui; padding: .4rem .85rem; border: 1px solid #1d3557;
         background: #1d3557; color: #fff; border-radius: 8px; cursor: pointer; }
button.ghost { background: #fff; color: #1d3557; }
button:disabled { opacity: .45; cursor: not-allowed; }
.info { font-size: .85rem; color: #555; min-height: 1.3em; margin-bottom: .4rem; }
canvas { display: block; width: 100%; max-width: 480px; border: 1px solid #d0d8e0;
         border-radius: 8px; background: #f8fafc; }
#table-wrap { margin-top: .5rem; max-height: 160px; overflow-y: auto; }
table { width: 100%; border-collapse: collapse; font-size: .82rem; }
th { background: #eef2f6; padding: .28rem .5rem; text-align: left; position: sticky; top: 0; }
td { padding: .22rem .5rem; border-bottom: 1px solid #eef2f6; }
td.bar-cell { width: 50%; }
.bar-inner { height: 10px; border-radius: 4px; background: #1d6fa5; transition: width .25s; }
tr.last-row { background: #f0f8f0; }

// Code not found

Notice how each iteration never increases the residual — that is the guarantee of the minimisation. A well-conditioned system collapses in a handful of steps; a poorly conditioned one converges more slowly. Press New system to generate a fresh random matrix and observe how the convergence curve changes.

GMRES belongs to the family of Krylov-subspace methods. Its theoretical backbone is the Arnoldi process, which builds an orthonormal basis $\{q_1, q_2, \dots, q_k\}$ for $\mathcal{K}_k(A, r_0)$ via a modified Gram-Schmidt procedure. At each step $k$ the algorithm solves a small least-squares problem of size $k$ to find the residual-minimising coefficient vector.

• Termination guarantee. In exact arithmetic GMRES terminates in at most $n$ steps, because $\mathcal{K}_n$ spans all of $\mathbb{R}^n$. In practice machines stop after dozens or hundreds of steps, not millions.
• Restarted GMRES. Storing the growing basis costs $\mathcal{O}(kn)$ memory. To cap this, practitioners restart every $m$ iterations (written GMRES($m$)), trading some convergence speed for bounded storage.
• Preconditioning is the real art. Multiplying by a matrix $M \approx A^{-1}$ that is cheap to apply transforms $Ax = b$ into $M^{-1}Ax = M^{-1}b$, clustering the eigenvalues and often cutting the iteration count from thousands to tens. The theory of why clustering helps is tied to the Chebyshev polynomial approximation of $A^{-1}$ on the eigenvalue spectrum.
• No convergence guarantee without structure. Unlike conjugate gradient, which is optimal for symmetric positive-definite systems, GMRES can stagnate on pathological non-symmetric matrices. Preconditioning is not optional for hostile problems.

The cost per iteration is $\mathcal{O}(kn)$ for the Gram-Schmidt orthogonalisation plus $\mathcal{O}(\mathrm{nnz}(A))$ for the matrix-vector product, where $\mathrm{nnz}$ is the number of non-zero entries.

Any time science or engineering must solve a large sparse linear system, GMRES (or a variant) is likely running under the hood:

• Computational fluid dynamics (CFD). Discretising the Navier-Stokes equations on a mesh produces a non-symmetric system; GMRES with an ILU preconditioner is a workhorse solver.
• Finite-element structural analysis. Car-crash simulations, bridge-load calculations and turbine-blade stress all reduce to $Ax = b$ with $n$ in the millions.
• Machine learning. Gaussian process regression, Bayesian neural networks and kernel methods require solving linear systems whose matrices are dense but structured — iterative solvers exploit that structure.
• Quantum chemistry. Computing molecular orbitals and reaction energies involves solving eigenvalue problems iteratively, often using GMRES-like Krylov steps inside an outer loop.
• Graphics and image processing. Solving the Poisson equation for lighting simulation or image in-painting amounts to a large sparse system, ideal for Krylov methods.

Whenever you see a technical report cite "PETSc", "Trilinos", "SciPy sparse.linalg" or "MATLAB backslash on a sparse matrix", there is a strong chance GMRES (or the closely related conjugate gradient) is doing the work.

GMRES is a quiet workhorse. You will not see it in a headline, but it is running every time a supercomputer models the atmosphere, every time a finite-element solver checks whether a bridge will stand, and every time a chip-design tool verifies that a circuit is correct.

Its elegance is that it never asks you to factor $A$ — it only asks you to multiply by $A$, one vector at a time. From those multiplications it builds the smallest subspace where the answer lives and finds the best approximation inside it. In exact arithmetic it always wins in $n$ steps; in practice, with a good preconditioner, it wins in far fewer.

Iterative solvers like GMRES are a reminder that in numerical computing, as in life, the direct route is often impossible — but a well-chosen sequence of small, improving steps can take you remarkably far.