NLA MT25, GMRES algorithm

Minimise the residual in the Krylov subspace:

What GMRES computes. GMRES (Generalised Minimal RESidual) finds at each iteration $k = 1, 2, \ldots$, the element of the Krylov subspace $\mathcal K _ k(A, b)$ that minimises the residual

The iterations are: extend the Krylov subspace by one dimension, solve the (slightly larger) projected least-squares problem, get a new iterate $x _ k$. The residual is non-increasing in $k$, and convergence is determined by ∆gmres-convergence.

The algorithm, at iteration $k$:

1, Extend Arnoldi by one step to obtain $Q _ k \in \mathbb R^{n \times k}$, $Q _ {k+1} \in \mathbb R^{n \times (k+1)}$ with orthonormal columns spanning $\mathcal K _ k(A, b)$ and $\mathcal K _ {k+1}(A, b)$ respectively, and $\tilde H _ k \in \mathbb R^{(k+1) \times k}$ upper hessenberg, satisfying the Arnoldi decomposition

Cost: one matrix-vector product with $A$, plus $O(nk)$ for orthogonalising against the existing basis. Cumulative cost over $k$ steps: $O(nk^2)$.

2, Reduce the residual minimisation to a small Hessenberg least-squares problem. Writing $x = Q _ k y$ for $y \in \mathbb R^k$, and using that $Q _ {k+1}$ has orthonormal columns:

Since the first column of $Q _ {k+1}$ is $q _ 1 = b/\|b\| _ 2$ (from Arnoldi), $Q^\top _ {k+1} b = \|b\| _ 2 e _ 1$ exactly, and so there is no need to ever compute $Q _ {k+1}^\top b$ as an inner product. The problem then becomes

3, Solve the Hessenberg least-squares problem via $k$ Givens rotations to get the QR factorisation of $\tilde H _ k$, then a triangular back-substitution. Because $\tilde H _ k$ is upper Hessenberg (only one subdiagonal), this is much cheaper than a general $(k+1) \times k$ least-squares problem. Moreover, as $k$ increases, $\tilde H _ {k+1}$ only adds one new column and one new row to $\tilde H _ k$, so the previous Givens factorisation can be updated incrementally.

Cost: $O(k)$ per iteration using incremental update. Cumulative cost over $k$ steps: $O(k^2)$.

4, Calculate residual. After step $k$, the GMRES residual norm equals

i.e. the absolute value of the last entry of the rotated RHS $g = G _ k \cdots G _ 1 (\|b\| _ 2 e _ 1) \in \mathbb R^{k+1}$. This is available at $O(1)$ cost per step, so there is no need to compute $y$ or $x _ k$ to check whether there is convergence.

5, Recover the iterate (only when stopping). Once you have decided to terminate at step $k$, back-solve $R _ k y = g _ {1 : k}$ (cost $O(k^2)$), and form $x _ k = Q _ k y$ (cost $O(nk)$).

Total cost over $k$ iterations:

• $O(nk^2)$ for Arnoldi
• $O(k^3)$ for cumulative Hessenberg least-squares solves, plus $k$ matrix-vector multiplications.

The memory cost is $O(nk)$ for storing $Q _ {k+1}$.

We aim to minimise the residual in the Krylov subspace:

Write $x = Q _ k y$. Then

This is minimised when $ \vert \vert \tilde H _ k y - Q^\top _ {k+1}b \vert \vert _ 2$ is minimised. This a Hessenberg least-squares problem and can be solved efficiently via QR.

Triangularise $\tilde H _ k$. Being upper Hessenberg, $\tilde H _ k$ has a single nonzero subdiagonal, so one Givens rotation per column clears it: $G _ j$ acts on rows $j, j+1$ to zero the subdiagonal entry in column $j$. After $k$ rotations,

Rotate the RHS. Apply the same rotations to the right-hand side, $g = G _ k \cdots G _ 1 (\|b\| _ 2 e _ 1) \in \mathbb R^{k+1}$, and split $g = \begin{bmatrix} g _ {1:k} \\ g _ {k+1} \end{bmatrix}$. The rotations are orthogonal, so they preserve the objective:

Solve and read off. The first term is killed by the triangular solve $R _ k y = g _ {1:k}$ (back-substitution); the second is irreducible. So the minimiser is $y = R _ k^{-1} g _ {1:k}$, and the minimal residual norm is $ \vert g _ {k+1} \vert $, the last entry of $g$, which needs no $y$ to evaluate.

Incremental update (the point). From step $k-1$ to $k$, $\tilde H _ k$ only appends one column and one subdiagonal entry to $\tilde H _ {k-1}$. So: apply the stored rotations $G _ 1, \ldots, G _ {k-1}$ to the new column, compute one new rotation $G _ k$ to zero $h _ {k+1,k}$, and apply $G _ k$ to $g$ (updating its last two entries and exposing the new $g _ {k+1}$). Nothing already computed is redone.

Cost. $O(k)$ per step (apply $\le k$ rotations to one length-$(k+1)$ column), hence $O(k^2)$ over $k$ steps, negligible beside Arnoldi’s $O(nk^2)$. The residual $ \vert g _ {k+1} \vert $ is available for free each step, so convergence is monitored without forming $y$ or $x _ k$; the back-solve $R _ k y = g _ {1:k}$ and lift $x _ k = Q _ k y$ are done once, only on stopping.

Why Givens, not Householder. The Hessenberg structure leaves a single entry per column to zero, so a local two-row Givens rotation is exactly the right tool, and it is what makes the append-one-column update cheap as $k$ grows.

Suppose:

• $A$ is a diagonalisable matrix, so that $A = X \Lambda X^{-1}$

Then the $k$-th GMRES iterate $x _ k$ satisfies

Recall that $x _ k \in \mathcal K _ k(A, b)$ implies $x _ k = p _ {k-1}(A)b$ where $p _ {k-1}$ is a polynomial of degree at most $k-1$. Hence the GMRES solution is

As $A$ is diagonalisable, we have

Hence overall

GMRES converges fast on matrices $A$ for which there exists a low-degree polynomial $p$ with $p(0) = 1$ and $ \vert p \vert $ small at every eigenvalue of $A$. Two common situations where this happens might be:

1. Eigenvalues clustered away from $0$. If $\lambda(A) \supset \overline{B(c, r)}$ for some centre $c$ with $ \vert c \vert > r$, then $p(z) = (1 - z/c)^k$ satisfies $p(0) = 1$ and $ \vert p(\lambda) \vert \le (r/ \vert c \vert )^k \to 0$ geometrically.
2. Few distinct eigenvalues. If $A$ has $k \ll n$ distinct eigenvalues $\lambda _ 1, \ldots, \lambda _ k$, set $p(z) = \prod^k _ {i = 1} \frac{(z - \lambda _ i)}{(0 - \lambda _ i)}$, which satisfies $p(0) = 1$ and $ \vert p(\lambda _ i) \vert = 0$ at every eigenvalue of $A$. So GMRES converges exactly in $\le k$ steps.

The condition number of the eigenvector matrix also matters. For normal $A$ it equals $1$; for non-normal $A$ it can amplify the bound substantially.

• In the first example, convergence is very fast because there are many eigenvalues clustered away from $0$.
• In the second example, convergence is very slow because there are many distinct eigenvalues all around $0$.

See also ∆gmres-fast-convergence-examples.

• Applying $M$ to a vector should be easy, since GMRES requires matrix-vector multiplications
• At least one of:
  1. $MA$ has clustered eigenvalues away from $0$
  2. $MA$ has a small number of distinct nonzero eigenvalues
  3. $MA$ is well-conditioned $\kappa _ 2(MA) = O(1)$

Setup: for a sparse $A$ (e.g. from a discretised PDE) the exact LU factorisation suffers fill-in — the elimination update $a _ {ij} \leftarrow a _ {ij} - \ell _ {ik} u _ {kj}$ creates nonzeros in positions where $A _ {ij} = 0$, so the true $L, U$ are far denser than $A$ (an $O(n)$-nonzero $A$ can give $O(n^2)$-nonzero factors): too expensive to form, store, and apply.

Incomplete: constrain $L, U$ to a prescribed sparsity pattern and discard any fill-in landing outside it. The simplest choice, ILU(0), keeps exactly the sparsity pattern of $A$ — run the elimination but zero out every update that would fall in a structurally-zero position. (Richer variants allow controlled extra fill: ILU($k$) by level-of-fill, ILUT by a drop tolerance — trading accuracy for density.)

The factorisation is therefore inexact: $\hat L \hat U = A + E$, where $E$ collects the dropped fill, so $\hat L \hat U \approx A$ but $\ne A$. Set

Then $MA = (A + E)^{-1} A \approx I$ (spectrum clustered near $1$ when $\|E\|$ is small) — the approximation comes entirely from the incomplete-factorisation step, not the inverse. Applying $M$ to a vector is two sparse triangular solves, $O(\mathrm{nnz})$ rather than $O(n^2)$, which is exactly the “easy to apply” property a preconditioner needs.

Then if $M$ is nonsingular, $M \mathcal A$ has eigenvalues $\{1, \frac 1 2 (1 \pm \sqrt 5)\}$ — only three distinct values, and so by ∆gmres-fast-convergence-examples (case 2) GMRES has three-step convergence on this preconditioned system.

The ideal preconditioner would be $M = A^{-1}$, since then

Algorithm:

1. Stop GMRES after $k _ \max$ iterations to get an approximate solution $\hat x _ 1$.
2. Solve $A\tilde x = b - A \tilde x _ 1$ via GMRES (i.e. fresh linear system with a new right-hand side, namely the current residual)
3. Return $\hat x _ 1 + \tilde x$ as the improved solution

If one restart isn’t enough, repeat.

Why useful: GMRES costs $O(nk^2)$ over $k$ iterations and stores $O(nk)$ for the Arnoldi basis $Q _ {k+1}$. Both costs grow in $k$, so for ill-conditioned or large problems the iteration becomes expensive and infeasible in memory before convergence. Restarting caps $k$ at a fixed budget $k _ \max$.

Trade-off: This means that GMRES at iteration $k$ no longer gives the exact solution to $x = \text{argmin} _ {x \in \mathcal K _ k}\|Ax - b\| _ 2$, which can cause issues with the convergence analysis.

Work in the affine space $x _ k = x _ 0 + \mathcal K _ k(A, r _ 0)$ where $r _ 0 = b - A x _ 0$ is the initial residual. Writing $x _ k = x _ 0 + z$, the minimisation becomes $\min _ {z \in \mathcal K _ k(A, r _ 0)} \|A(x _ 0 + z) - b\| _ 2 = \min _ z \|A z - r _ 0\| _ 2$ — i.e. plain GMRES applied to $A z = r _ 0$. The Arnoldi process is started from $r _ 0/\|r _ 0\| _ 2$ rather than $b/\|b\| _ 2$.

For a non-normal $A$, the eigenvector matrix $X$ in $A = X \Lambda X^{-1}$ can have arbitrarily large $\kappa _ 2(X)$ — Jordan-block-like matrices have ill-conditioned eigenvectors. Even when the polynomial factor $\max _ {z \in \lambda(A)} \vert p(z) \vert $ is small, the bound $\kappa _ 2(X) \cdot \max \vert p \vert $ can be vacuous. Refined “pseudospectral” GMRES bounds (non-examinable) avoid this dependence.

• Translate the residual-minimisation into a polynomial-approximation problem. Every $x _ k \in \mathcal K _ k(A, b)$ has $x _ k = p _ {k-1}(A) b$, so $Ax _ k - b = (A p _ {k-1}(A) - I) b = -p(A)b$ where $p \in \mathcal P _ k$ with $p(0) = 1$.
• Diagonalise $A = X \Lambda X^{-1}$, giving $p(A) = X p(\Lambda) X^{-1}$. Submultiplicativity then bounds $\|p(A)\| _ 2 \le \kappa _ 2(X) \max _ {z \in \lambda(A)} \vert p(z) \vert $.
• Take the min over $p \in \mathcal P _ k$ with $p(0) = 1$. Because GMRES chooses the best $x _ k$ in the Krylov subspace, it also achieves this $\min _ p$ bound, giving the result.