GMRES: Solving the Unsolvable One Step at a Time

Introduction

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.

Try It

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.

<!-- GMRES demo — residual convergence visualisation -->
<div class="controls">
  <button id="btn-step" type="button">Step</button>
  <button id="btn-run" type="button">Run all</button>
  <button id="btn-new" type="button" class="ghost">New system</button>
</div>
<div id="info" class="info"></div>
<canvas id="chart" width="480" height="200" aria-label="Residual convergence chart for GMRES"></canvas>
<div id="table-wrap">
  <table id="iter-table">
    <thead><tr><th>k</th><th>‖r_k‖</th><th>Relative size</th></tr></thead>
    <tbody id="tbody"></tbody>
  </table>
</div>

/* Layout and chart styles */
* { 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; }

/* GMRES with Arnoldi + Givens rotations (exact least-squares each step) */
const i18nFmt = (s, v) => s.replace(/%(\w+)%/g, (_, k) => String(v[k] ?? ''));

/* Generate a random diagonally-dominant n×n matrix and RHS vector */
function randMatrix(n, seed) {
  let s = seed || 42;
  const rng = () => { s = (s * 1664525 + 1013904223) & 0xffffffff; return (s >>> 0) / 0xffffffff; };
  const A = [], b = [];
  for (let i = 0; i < n; i++) {
    A.push([]);
    for (let j = 0; j < n; j++) {
      /* Diagonal dominance ensures convergence for the demo */
      const v = rng() * 2 - 1;
      A[i].push(i === j ? n + rng() * 2 : v * 0.8);
    }
    b.push(rng() * 4 - 2);
  }
  return { A, b };
}

/* Basic dense vector/matrix utilities */
const dot = (a, b) => a.reduce((s, v, i) => s + v * b[i], 0);
const norm = v => Math.sqrt(dot(v, v));
const scale = (v, c) => v.map(x => x * c);
const sub = (a, b) => a.map((x, i) => x - b[i]);
const add = (a, b) => a.map((x, i) => x + b[i]);
const matvec = (A, v) => A.map(row => dot(row, v));
const zeros = n => Array(n).fill(0);

/* One step of the Arnoldi process — extends the Krylov basis by one vector */
function arnoldiStep(A, Q, k) {
  let q = matvec(A, Q[k]);
  const h = zeros(k + 2);
  for (let j = 0; j <= k; j++) {
    h[j] = dot(q, Q[j]);
    q = sub(q, scale(Q[j], h[j]));
  }
  h[k + 1] = norm(q);
  const qNext = h[k + 1] > 1e-14 ? scale(q, 1 / h[k + 1]) : zeros(A.length);
  return { h, qNext };
}

/* Apply accumulated Givens rotations to new Hessenberg column, update QR */
function applyGivens(H, cs, sn, k) {
  for (let i = 0; i < k; i++) {
    const tmp = cs[i] * H[i] + sn[i] * H[i + 1];
    H[i + 1] = -sn[i] * H[i] + cs[i] * H[i + 1];
    H[i] = tmp;
  }
  const r = Math.sqrt(H[k] * H[k] + H[k + 1] * H[k + 1]);
  cs[k] = H[k] / r; sn[k] = H[k + 1] / r;
  H[k] = cs[k] * H[k] + sn[k] * H[k + 1]; H[k + 1] = 0;
  return r;
}

const N = 6;
let state;

function initState(seed) {
  const { A, b } = randMatrix(N, seed);
  const r0 = b.slice(); /* Start from x0 = 0, so initial residual r0 = b */
  const beta = norm(r0);
  const Q = [scale(r0, 1 / beta)];
  state = {
    A, b, beta,
    Q, /* Orthonormal Krylov basis vectors Q[0..k] */
    Hfull: [], /* Upper Hessenberg matrix columns (not used for solution here, only convergence tracking) */
    cs: [], sn: [], /* Cosine/sine pairs for Givens QR of the Hessenberg matrix */
    g: [beta], /* Rotated RHS vector; |g[k+1]| = residual norm at step k */
    residuals: [beta],
    done: false,
    k: 0,
    seed
  };
}

function stepGMRES() {
  if (state.done) return;
  const { A, Q, cs, sn, g, k } = state;
  const { h, qNext } = arnoldiStep(A, Q, k);
  state.Hfull.push(h.slice());
  Q.push(qNext);
  const col = h.slice(0, k + 2);
  applyGivens(col, cs, sn, k);
  g.push(0);
  g[k + 1] = -sn[k] * g[k];
  g[k] = cs[k] * g[k];
  const res = Math.abs(g[k + 1]);
  state.residuals.push(res);
  state.k++;
  if (res < 1e-10 || state.k >= N) state.done = true;
}

/* Draw residual convergence curve on the canvas (log scale) */
function draw() {
  const cvs = document.getElementById('chart');
  const ctx = cvs.getContext('2d');
  const W = cvs.width, H = cvs.height;
  ctx.clearRect(0, 0, W, H);

  const resids = state.residuals;
  const maxR = Math.max(...resids);
  const logMax = Math.log10(maxR + 1e-15);
  const logMin = Math.min(Math.log10(1e-11), Math.log10(Math.min(...resids)) - 0.5);
  const pad = { l: 52, r: 14, t: 16, b: 36 };
  const cw = W - pad.l - pad.r, ch = H - pad.t - pad.b;

  /* Horizontal grid lines at integer powers of 10 */
  ctx.strokeStyle = '#d5dfe8'; ctx.lineWidth = 1;
  for (let exp = Math.ceil(logMin); exp <= Math.ceil(logMax) + 1; exp++) {
    const y = pad.t + ch * (1 - (exp - logMin) / (logMax - logMin));
    if (y < pad.t || y > pad.t + ch) continue;
    ctx.beginPath(); ctx.moveTo(pad.l, y); ctx.lineTo(pad.l + cw, y); ctx.stroke();
    ctx.fillStyle = '#888'; ctx.font = '10px system-ui'; ctx.textAlign = 'right';
    ctx.fillText('10' + exp, pad.l - 4, y + 3);
  }

  /* Blue polyline connecting (k, log‖r_k‖) points */
  ctx.strokeStyle = '#1d6fa5'; ctx.lineWidth = 2;
  ctx.beginPath();
  resids.forEach((r, i) => {
    const x = pad.l + (i / (N)) * cw;
    const logR = Math.log10(Math.max(r, 1e-15));
    const y = pad.t + ch * (1 - (logR - logMin) / (logMax - logMin));
    i === 0 ? ctx.moveTo(x, y) : ctx.lineTo(x, y);
  });
  ctx.stroke();

  /* Circle markers at each data point */
  ctx.fillStyle = '#1d6fa5';
  resids.forEach((r, i) => {
    const x = pad.l + (i / (N)) * cw;
    const logR = Math.log10(Math.max(r, 1e-15));
    const y = pad.t + ch * (1 - (logR - logMin) / (logMax - logMin));
    ctx.beginPath(); ctx.arc(x, y, 4, 0, 2 * Math.PI); ctx.fill();
  });

  /* Axis labels */
  ctx.fillStyle = '#555'; ctx.font = '11px system-ui'; ctx.textAlign = 'center';
  ctx.fillText('Iteration k', pad.l + cw / 2, H - 4);
  ctx.save(); ctx.translate(12, pad.t + ch / 2); ctx.rotate(-Math.PI / 2);
  ctx.fillText('‖r_k‖ (log scale)', 0, 0); ctx.restore();
}

function updateTable() {
  const tbody = document.getElementById('tbody');
  tbody.innerHTML = '';
  const maxR = state.residuals[0] || 1;
  state.residuals.forEach((r, i) => {
    const tr = document.createElement('tr');
    if (i === state.residuals.length - 1) tr.className = 'last-row';
    const pct = Math.max(2, (r / maxR) * 100).toFixed(1);
    tr.innerHTML = `<td>${i}</td><td>${r.toExponential(3)}</td><td class="bar-cell"><div class="bar-inner" style="width:${pct}%"></div></td>`;
    tbody.appendChild(tr);
  });
  tbody.lastElementChild && tbody.lastElementChild.scrollIntoView({ block: 'nearest' });
}

function updateInfo() {
  const el = document.getElementById('info');
  const k = state.k;
  const r = state.residuals[k];
  if (state.done && r < 1e-10) {
    el.textContent = i18nFmt('Converged in %k% iteration(s) — residual below 1e-10.', { k });
  } else if (state.done) {
    el.textContent = i18nFmt('Finished %k% iteration(s) — maximum Krylov dimension reached.', { k });
  } else {
    el.textContent = i18nFmt('After iteration %k%: residual = %r%', { k, r: r.toExponential(3) });
  }
  document.getElementById('btn-step').disabled = state.done;
  document.getElementById('btn-run').disabled = state.done;
}

function refresh() { draw(); updateTable(); updateInfo(); }

let runTimer = null;
document.getElementById('btn-step').addEventListener('click', () => {
  stepGMRES(); refresh();
});
document.getElementById('btn-run').addEventListener('click', () => {
  if (runTimer) return;
  /* Advance one iteration every 350 ms until convergence */
  runTimer = setInterval(() => {
    if (state.done) { clearInterval(runTimer); runTimer = null; return; }
    stepGMRES(); refresh();
  }, 350);
});
document.getElementById('btn-new').addEventListener('click', () => {
  clearInterval(runTimer); runTimer = null;
  initState(Math.floor(Math.random() * 9000) + 100);
  refresh();
});

/* Initialise with seed 42 and render */
initState(42);
refresh();

<!-- GMRES demo — residual convergence visualisation -->
<div class="controls">
  <button id="btn-step" type="button">Step</button>
  <button id="btn-run" type="button">Run all</button>
  <button id="btn-new" type="button" class="ghost">New system</button>
</div>
<div id="info" class="info"></div>
<canvas id="chart" width="480" height="200" aria-label="Residual convergence chart for GMRES"></canvas>
<div id="table-wrap">
  <table id="iter-table">
    <thead><tr><th>k</th><th>‖r_k‖</th><th>Relative size</th></tr></thead>
    <tbody id="tbody"></tbody>
  </table>
</div>

/* Layout and chart styles */
* { 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; }

/* GMRES with Arnoldi + Givens rotations (exact least-squares each step) */
const i18nFmt = (s, v) => s.replace(/%(\w+)%/g, (_, k) => String(v[k] ?? ''));

/* Generate a random diagonally-dominant n×n matrix and RHS vector */
function randMatrix(n, seed) {
  let s = seed || 42;
  const rng = () => { s = (s * 1664525 + 1013904223) & 0xffffffff; return (s >>> 0) / 0xffffffff; };
  const A = [], b = [];
  for (let i = 0; i < n; i++) {
    A.push([]);
    for (let j = 0; j < n; j++) {
      /* Diagonal dominance ensures convergence for the demo */
      const v = rng() * 2 - 1;
      A[i].push(i === j ? n + rng() * 2 : v * 0.8);
    }
    b.push(rng() * 4 - 2);
  }
  return { A, b };
}

/* Basic dense vector/matrix utilities */
const dot = (a, b) => a.reduce((s, v, i) => s + v * b[i], 0);
const norm = v => Math.sqrt(dot(v, v));
const scale = (v, c) => v.map(x => x * c);
const sub = (a, b) => a.map((x, i) => x - b[i]);
const add = (a, b) => a.map((x, i) => x + b[i]);
const matvec = (A, v) => A.map(row => dot(row, v));
const zeros = n => Array(n).fill(0);

/* One step of the Arnoldi process — extends the Krylov basis by one vector */
function arnoldiStep(A, Q, k) {
  let q = matvec(A, Q[k]);
  const h = zeros(k + 2);
  for (let j = 0; j <= k; j++) {
    h[j] = dot(q, Q[j]);
    q = sub(q, scale(Q[j], h[j]));
  }
  h[k + 1] = norm(q);
  const qNext = h[k + 1] > 1e-14 ? scale(q, 1 / h[k + 1]) : zeros(A.length);
  return { h, qNext };
}

/* Apply accumulated Givens rotations to new Hessenberg column, update QR */
function applyGivens(H, cs, sn, k) {
  for (let i = 0; i < k; i++) {
    const tmp = cs[i] * H[i] + sn[i] * H[i + 1];
    H[i + 1] = -sn[i] * H[i] + cs[i] * H[i + 1];
    H[i] = tmp;
  }
  const r = Math.sqrt(H[k] * H[k] + H[k + 1] * H[k + 1]);
  cs[k] = H[k] / r; sn[k] = H[k + 1] / r;
  H[k] = cs[k] * H[k] + sn[k] * H[k + 1]; H[k + 1] = 0;
  return r;
}

const N = 6;
let state;

function initState(seed) {
  const { A, b } = randMatrix(N, seed);
  const r0 = b.slice(); /* Start from x0 = 0, so initial residual r0 = b */
  const beta = norm(r0);
  const Q = [scale(r0, 1 / beta)];
  state = {
    A, b, beta,
    Q, /* Orthonormal Krylov basis vectors Q[0..k] */
    Hfull: [], /* Upper Hessenberg matrix columns (not used for solution here, only convergence tracking) */
    cs: [], sn: [], /* Cosine/sine pairs for Givens QR of the Hessenberg matrix */
    g: [beta], /* Rotated RHS vector; |g[k+1]| = residual norm at step k */
    residuals: [beta],
    done: false,
    k: 0,
    seed
  };
}

function stepGMRES() {
  if (state.done) return;
  const { A, Q, cs, sn, g, k } = state;
  const { h, qNext } = arnoldiStep(A, Q, k);
  state.Hfull.push(h.slice());
  Q.push(qNext);
  const col = h.slice(0, k + 2);
  applyGivens(col, cs, sn, k);
  g.push(0);
  g[k + 1] = -sn[k] * g[k];
  g[k] = cs[k] * g[k];
  const res = Math.abs(g[k + 1]);
  state.residuals.push(res);
  state.k++;
  if (res < 1e-10 || state.k >= N) state.done = true;
}

/* Draw residual convergence curve on the canvas (log scale) */
function draw() {
  const cvs = document.getElementById('chart');
  const ctx = cvs.getContext('2d');
  const W = cvs.width, H = cvs.height;
  ctx.clearRect(0, 0, W, H);

  const resids = state.residuals;
  const maxR = Math.max(...resids);
  const logMax = Math.log10(maxR + 1e-15);
  const logMin = Math.min(Math.log10(1e-11), Math.log10(Math.min(...resids)) - 0.5);
  const pad = { l: 52, r: 14, t: 16, b: 36 };
  const cw = W - pad.l - pad.r, ch = H - pad.t - pad.b;

  /* Horizontal grid lines at integer powers of 10 */
  ctx.strokeStyle = '#d5dfe8'; ctx.lineWidth = 1;
  for (let exp = Math.ceil(logMin); exp <= Math.ceil(logMax) + 1; exp++) {
    const y = pad.t + ch * (1 - (exp - logMin) / (logMax - logMin));
    if (y < pad.t || y > pad.t + ch) continue;
    ctx.beginPath(); ctx.moveTo(pad.l, y); ctx.lineTo(pad.l + cw, y); ctx.stroke();
    ctx.fillStyle = '#888'; ctx.font = '10px system-ui'; ctx.textAlign = 'right';
    ctx.fillText('10' + exp, pad.l - 4, y + 3);
  }

  /* Blue polyline connecting (k, log‖r_k‖) points */
  ctx.strokeStyle = '#1d6fa5'; ctx.lineWidth = 2;
  ctx.beginPath();
  resids.forEach((r, i) => {
    const x = pad.l + (i / (N)) * cw;
    const logR = Math.log10(Math.max(r, 1e-15));
    const y = pad.t + ch * (1 - (logR - logMin) / (logMax - logMin));
    i === 0 ? ctx.moveTo(x, y) : ctx.lineTo(x, y);
  });
  ctx.stroke();

  /* Circle markers at each data point */
  ctx.fillStyle = '#1d6fa5';
  resids.forEach((r, i) => {
    const x = pad.l + (i / (N)) * cw;
    const logR = Math.log10(Math.max(r, 1e-15));
    const y = pad.t + ch * (1 - (logR - logMin) / (logMax - logMin));
    ctx.beginPath(); ctx.arc(x, y, 4, 0, 2 * Math.PI); ctx.fill();
  });

  /* Axis labels */
  ctx.fillStyle = '#555'; ctx.font = '11px system-ui'; ctx.textAlign = 'center';
  ctx.fillText('Iteration k', pad.l + cw / 2, H - 4);
  ctx.save(); ctx.translate(12, pad.t + ch / 2); ctx.rotate(-Math.PI / 2);
  ctx.fillText('‖r_k‖ (log scale)', 0, 0); ctx.restore();
}

function updateTable() {
  const tbody = document.getElementById('tbody');
  tbody.innerHTML = '';
  const maxR = state.residuals[0] || 1;
  state.residuals.forEach((r, i) => {
    const tr = document.createElement('tr');
    if (i === state.residuals.length - 1) tr.className = 'last-row';
    const pct = Math.max(2, (r / maxR) * 100).toFixed(1);
    tr.innerHTML = `<td>${i}</td><td>${r.toExponential(3)}</td><td class="bar-cell"><div class="bar-inner" style="width:${pct}%"></div></td>`;
    tbody.appendChild(tr);
  });
  tbody.lastElementChild && tbody.lastElementChild.scrollIntoView({ block: 'nearest' });
}

function updateInfo() {
  const el = document.getElementById('info');
  const k = state.k;
  const r = state.residuals[k];
  if (state.done && r < 1e-10) {
    el.textContent = i18nFmt('Converged in %k% iteration(s) — residual below 1e-10.', { k });
  } else if (state.done) {
    el.textContent = i18nFmt('Finished %k% iteration(s) — maximum Krylov dimension reached.', { k });
  } else {
    el.textContent = i18nFmt('After iteration %k%: residual = %r%', { k, r: r.toExponential(3) });
  }
  document.getElementById('btn-step').disabled = state.done;
  document.getElementById('btn-run').disabled = state.done;
}

function refresh() { draw(); updateTable(); updateInfo(); }

let runTimer = null;
document.getElementById('btn-step').addEventListener('click', () => {
  stepGMRES(); refresh();
});
document.getElementById('btn-run').addEventListener('click', () => {
  if (runTimer) return;
  /* Advance one iteration every 350 ms until convergence */
  runTimer = setInterval(() => {
    if (state.done) { clearInterval(runTimer); runTimer = null; return; }
    stepGMRES(); refresh();
  }, 350);
});
document.getElementById('btn-new').addEventListener('click', () => {
  clearInterval(runTimer); runTimer = null;
  initState(Math.floor(Math.random() * 9000) + 100);
  refresh();
});

/* Initialise with seed 42 and render */
initState(42);
refresh();

Result

Console

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.

The Real Complexity

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.

Where It Matters

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.

Conclusion

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.

Comments

Loading comments...

https://www.kipuhub.com/en/article/gmres/Content licensed under CC BY-NC 4.0.

Cookies & privacy

GMRES

Finding the needle in a million-dimensional haystack