When a linear system \(\bA\bx = \bb\) is nonsymmetric or indefinite, the quadratic energy minimization framework of the Conjugate Gradient method ceases to apply. By the celebrated Faber-Manteuffel theorem, it is mathematically impossible for any Krylov subspace method on general matrices to simultaneously maintain a short three-term recurrence and satisfy an optimal error minimization property.
The Generalized Minimal Residual (GMRES) method of Saad and Schultz resolves this dilemma by retaining strict mathematical optimality: at every step \(k\), GMRES finds the unique vector \(\bx^{(k)}\) within the affine Krylov subspace \(\bx^{(0)} + \mathcal{K}_k(\bA, \br^{(0)})\) that minimizes the Euclidean norm of the residual \(\|\bb - \bA\bx\|_2\). To achieve this, GMRES constructs and stores an explicit orthonormal basis for the expanding Krylov subspace via the Arnoldi process.
The Arnoldi Process
The Arnoldi iteration is a Modified Gram-Schmidt orthogonalization procedure that computes an orthonormal basis for the Krylov subspace \(\mathcal{K}_k(\bA, \br^{(0)})\).
Let \(\bA \in \fR^{n \times n}\) and let \(\br^{(0)} = \bb - \bA\bx^{(0)}\) be a nonzero initial residual with norm \(\beta = \|\br^{(0)}\|_2\). Set the first basis vector to \(\bv_1 = \br^{(0)} / \beta\). For \(j = 1, 2, ..., k\):
Compute the new Krylov vector: \(\bw = \bA \bv_j\).
For \(i = 1, 2, ..., j\), compute the projection coefficients and orthogonalize: \[
\begin{align}
h_{ij} = \bv_i^T \bw, \qquad \bw \leftarrow \bw - h_{ij} \bv_i.
\end{align}
\]
Compute the subdiagonal entry: \(h_{j+1,j} = \|\bw\|_2\).
If \(h_{j+1,j} = 0\), the iteration terminates (happy breakdown). Otherwise, set \(\bv_{j+1} = \bw / h_{j+1,j}\).
After \(k\) steps of the Arnoldi iteration, the matrix \(\bV_k = [\bv_1, ..., \bv_k] \in \fR^{n \times k}\) has orthonormal columns spanning \(\mathcal{K}_k(\bA, \br^{(0)})\), and satisfies the algebraic relation \[
\begin{align}
\bA \bV_k = \bV_{k+1} \bar{\bH}_k,
\end{align}
\] where \(\bar{\bH}_k \in \fR^{(k+1) \times k}\) is an upper Hessenberg matrix whose entries are the orthogonalization coefficients \(h_{ij}\). Furthermore, \(\bH_k = \bV_k^T \bA \bV_k \in \fR^{k \times k}\) is the square upper Hessenberg matrix formed by omitting the final row of \(\bar{\bH}_k\).
Reduction to Hessenberg Least Squares
The critical algorithmic insight of GMRES is that the \(n\)-dimensional residual minimization problem over the Krylov subspace can be mapped isometrically onto a small \((k+1) \times k\) linear least-squares problem.
Let any vector in the affine space \(\bx^{(0)} + \mathcal{K}_k(\bA, \br^{(0)})\) be parameterized as \(\bx = \bx^{(0)} + \bV_k \by\) for some vector \(\by \in \fR^k\). The residual vector satisfies \[
\begin{align}
\br(\bx) = \bb - \bA(\bx^{(0)} + \bV_k \by) = \br^{(0)} - \bA \bV_k \by = \beta \bv_1 - \bV_{k+1} \bar{\bH}_k \by = \bV_{k+1} (\beta \be_1 - \bar{\bH}_k \by).
\end{align}
\] Because the columns of \(\bV_{k+1}\) are orthonormal, \(\bV_{k+1}\) preserves the Euclidean norm. Therefore, \[
\begin{align}
\|\bb - \bA\bx\|_2 = \|\beta \be_1 - \bar{\bH}_k \by\|_2,
\end{align}
\] where \(\be_1 = (1, 0, ..., 0)^T \in \fR^{k+1}\). The GMRES iterate is \(\bx^{(k)} = \bx^{(0)} + \bV_k \by_k\), where \(\by_k\) uniquely solves the \((k+1) \times k\) upper Hessenberg least-squares problem \[
\begin{align}
\by_k = \arg\min_{\by \in \fR^k} \|\beta \be_1 - \bar{\bH}_k \by\|_2.
\end{align}
\]
Substituting the Arnoldi identity \(\bA \bV_k = \bV_{k+1} \bar{\bH}_k\) and expressing \(\br^{(0)} = \beta \bv_1 = \bV_{k+1} (\beta \be_1)\) gives \(\br(\bx) = \bV_{k+1} (\beta \be_1 - \bar{\bH}_k \by)\). Taking the Euclidean norm and applying the isometry property \(\|\bV_{k+1}\bz\|_2 = \|\bz\|_2\) (which holds because \(\bV_{k+1}^T \bV_{k+1} = \bI_{k+1}\)) establishes the result.
Solution via Givens Rotations
Because \(\bar{\bH}_k\) is upper Hessenberg, it has nonzero entries only on the main diagonal, upper triangular entries, and a single subdiagonal. We solve the least-squares problem incrementally using a sequence of \(2 \times 2\) Givens rotations \(\bG_1, \bG_2, ..., \bG_k\).
At step \(k\), applying the previous rotations \(\bG_1, ..., \bG_{k-1}\) transforms the first \(k-1\) entries of the new column of \(\bar{\bH}_k\). The subdiagonal entry at \((k+1, k)\) is then eliminated using a new Givens rotation \(\bG_k\) defined by \[
\begin{align}
\bG_k = \begin{pmatrix} c_k & s_k \\ -s_k & c_k \end{pmatrix}, \qquad c_k = \frac{h_{kk}}{\sqrt{h_{kk}^2 + h_{k+1,k}^2}}, \quad s_k = \frac{h_{k+1,k}}{\sqrt{h_{kk}^2 + h_{k+1,k}^2}}.
\end{align}
\] Applying \(\bG_k ... \bG_1\) transforms \(\bar{\bH}_k\) into an upper triangular matrix \(\bR_k\) with a zero final row, and transforms \(\beta \be_1\) into a vector \(\bg_{k+1} = (\gamma_1, ..., \gamma_k, \gamma_{k+1})^T\). The residual norm at iteration \(k\) is obtained immediately without forming \(\bx^{(k)}\): \[
\begin{align}
\|\br^{(k)}\|_2 = |\gamma_{k+1}|.
\end{align}
\]
(Happy breakdown vs. stagnation) If the subdiagonal entry \(h_{k+1,k}\) becomes zero during the Arnoldi process, the Krylov subspace has become an invariant subspace of \(\bA\). In this case, the least-squares problem has residual zero, and the GMRES iterate \(\bx^{(k)}\) is the exact solution to the linear system. This event is termed happy breakdown. Conversely, if the residual norm ceases to decrease across iterations, the method stagnates, which typically indicates poor conditioning or an inadequate preconditioner.
Restarted GMRES and Practical Considerations
Full GMRES requires storing \(k\) basis vectors of length \(n\), which incurs \(O(kn)\) memory storage and \(O(k^2 n)\) orthogonalization flops. To avoid memory exhaustion on large problems, the algorithm is restarted every \(m\) steps:
Run \(m\) steps of GMRES starting from \(\bx^{(0)}\) to produce \(\bx^{(m)}\).
Compute the residual \(\br^{(m)} = \bb - \bA\bx^{(m)}\).
Clear the Krylov basis vectors, set \(\bx^{(0)} \leftarrow \bx^{(m)}\), and repeat.
(The trade-off of restarting) Restarted GMRES(\(m\)) bounds memory usage to \(O(mn)\) and bounds orthogonalization cost per cycle. However, restarting discards the accumulated Krylov subspace and destroys the global optimality property. If \(m\) is chosen too small, the algorithm may stagnate completely.
Perform two steps of the Arnoldi iteration by hand for \[
\begin{align}
\bA = \begin{pmatrix} 1 & 1 & 0 \\ 0 & 2 & 1 \\ 0 & 0 & 3 \end{pmatrix}, \qquad \br^{(0)} = \begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix}.
\end{align}
\] Explicitly write down \(\bV_2\) and the \(3 \times 2\) upper Hessenberg matrix \(\bar{\bH}_2\).
Using the Hessenberg matrix from the previous exercise, construct the Givens rotations that triangularize \(\bar{\bH}_2\), and compute the second GMRES residual norm.
Explain why GMRES applied to a symmetric positive definite matrix produces mathematically identical iterates to the Conjugate Gradient method, but requires significantly more memory.
Compare Left Preconditioning (\(\bM^{-1}\bA\bx = \bM^{-1}\bb\)) and Right Preconditioning (\(\bA\bM^{-1}\by = \bb, \bx = \bM^{-1}\by\)). Why is Right Preconditioning preferred when monitoring true residual norms?