# A numerical GMRES example

I'm having trouble understanding how GMRES works. I've read the part in Saad's book and a few others but still I am confused. Can someone provide me a numerical example to understand it better? Or if that's not possible maybe a walk through the algorithm step by step?

• It would help if you first explain the parts of GMRES that you understand. Then, you can isolate the parts of the algorithm that you don't understand. – Paul May 20 '16 at 21:53
• In particular, do you understand what the algorithm does mathematically? (There's two levels to GMRES: A high-level "step x solves this subproblem" and a low level "to solve this subproblem efficiently, carry out these steps". I at least find it helpful to tackle one after the other, instead of jumping straight into the (pseudo-)code.) – Christian Clason May 20 '16 at 23:33
• +1 on Christian's comment. I think it's unfortunate that even many standard texts on the topic (GMRES, CG etc) go only as far as to state the "physical" algorithm, rather than describe the intuition behind "why" these steps are taken. – Richard Zhang May 20 '16 at 23:59
• Thank you for you answers. I understand that the method is trying to minimize the residual on a Krylov Subspace that is generated by the Arnoldi method, in a very theoretical way, I'm not sure if that's what you mean by mathematically @ChristianClason. Like Richard said I don't understand the intuition behind it. – System May 21 '16 at 8:07
• @System That is exactly what I meant. – Christian Clason May 21 '16 at 14:53

GMRES is indeed one of the hardest to understand Krylov methods. As you correctly state, the algorithm computes in each step $m$ a new approximation $x^m$ to the solution of $Ax=b$ as a minimizer of $\|b-Ax\|_2$ over the Krylov space $K_m(A,b) = \mathrm{span}\{b,Ab,A^2b,\dots,A^{m-1}b\}$. (I'm assuming for simplicity that $x^0 =0$). Clearly, if $A\in\mathbb{R}^{n\times n}$, then $x^n$ is the desired solution. The whole point of GMRES is to compute each minimizer $x^{m}$ via the Arnoldi process by reusing as much information from the previous steps as possible.

The derivation proceeds in two steps. First, recall that the Arnoldi process (which is an algorithm analogous to Gram-Schmidt to compute an orthonormal basis of $K_m(A,b)$) yields after $m$ steps the orthonormal vectors $q_1,\dots,q_{m+1}$ with $q_1 = b/\beta$, $\beta = \|{b}\|$, and scalars $h_{ij}$, $1\leq i\leq j+1\leq m+1$ such that $$AQ_m = Q_m H_m + h_{m+1,m}q_{m+1}e_m^T =: Q_{m+1}\bar H_m,$$ where $Q_{m+1}=[q_1,\dots,q_{m+1}] \in \mathbb{R}^{n\times (m+1)}$, $\bar H_m = (h_{ij})_{i,j} \in \mathbb{R}^{(m+1)\times m}$, and $e_m$ is the $m$th unit vector. Since every $x\in K_m(A,b)$ can be written as a linear combination of $q_1,\dots,q_m$, i.e., $x = Q_m\xi$ for some $\xi\in\mathbb{R}^m$, we have that \begin{aligned}[t] b-Ax &= b - AQ_m \xi\\ &= \beta q_1 - Q_{m+1}\bar H_m \xi \\ &= Q_{m+1}(\beta e_1 - \bar H_m\xi). \end{aligned} Since $Q_{m+1}$ is unitary, we hence have that $$\min_{x\in K_m(A,b)} \|{b-Ax}\|_2^2 = \min_{\xi\in\mathbb{R}^m} \|{\beta e_1 - \bar H_m\xi}\|_2^2.$$ The right-hand side is a (typically much) smaller least squares problem that can be solved by QR decomposition of $\bar H_m$: If $\bar H_m = GR$, where $G$ is unitary (I'm not using $Q$ here to avoid confusion) and $R$ is upper triangular, then $$\|\bar H_m\xi - \beta e_1\|^2_2 = \|\beta G^T e_1 -R\xi\|^2_2.$$ Since $\bar H_m\in\mathbb{R}^{(m+1)\times m}$, the matrix $R$ has the form $R=\begin{pmatrix}R_m\\0\end{pmatrix}$ with $R_m\in\mathbb{R}^{m\times m}$ invertible. Furthermore, we can partition $G^T e_1 =: (d_m,\rho_m)\in\mathbb{R}^{m+1}$ with $d_m\in\mathbb{R}^m$ and $\rho_m \in \mathbb{R}$. Thus, $$\|\beta G^T e_1 -R\xi\|^2_2 = \|R_m\xi - \beta d_m\|_2^2 + |\beta \rho_m|^2.$$ The minimum is therefore attained at $\xi^m := \beta R_m^{-1} d_m$, for which the residual is $\beta |\rho_m|$. This residuum -- which is used as a stopping criterion -- can therefore be evaluated without computing $x^m$. Furthermore, $\bar H_m$ already has upper Hessenberg form, hence this can be done quite fast by a single sweep of Givens rotations (whose product -- in matrix form -- is exactly $G^T$). The approximation $x^m$ is then simply obtained from $x^m =Q_m\xi^m$ (but this only needs to be done at the very end).

This is exactly what the GMRES algorithm does, but instead of starting from scratch each time, it alternates a step of the Arnoldi process and an update of the QR factorization to add a column: In step $m$,

1. A new set of $q_{m+1}$ and $h_{1,m},\dots,h_{m+1,m}$ are computed using the Arnoldi process. (This gives a new column of $Q$ as well as of $\bar H$.)

2. The extended matrix $\bar H_{m+1}$ is brought to upper triagonal form by applying $m+1$ Givens rotations to the new column -- of which the first $m$ have already been obtained during the previous steps, and only the last one has to be computed to zero out the bottom-most element. (This gives a new column of $R$ in $\bar H = GR$; note that the new rotation leaves the previous columns unchanged -- this is the key of GMRES.)

3. These Givens rotations are also applied to the right-hand side $e_1\in \mathbb{R}^{m+1}$ to obtain the vector $(d^m,\rho_m)$. Again, only the newest rotation needs actually be applied, because the older ones were already applied in the previous iterations.

4. If we are satisfied with the current relative residuum $|\rho_m|$ (or have reached the maximum number of iterations), the final iterate is computed via $x^m = \beta Q_m R_m^{-1}d^m$ (which is cheap since $R_m$ is small and upper triangular).

Since the steps become increasingly more expensive with $m$ (more data has to be stored, and the number of arithmetic operations grows quadratically), the algorithm is usually restarted after $M$ steps: $x^M$ is computed as above, and the procedure starts over with $m=0$ and a new initial point $x^0=x^M$.

Going into more detail here would get a bit long, but let me know and I'll try to fill in any gaps.

• Thank you for your answer. It's very informative and I think I get the gist of it now. In step 3 how is ρm obtained and why does it turn out to be the residual? – System May 25 '16 at 16:46
• I'm glad it's useful. I've added a few details about the solution of the small least squares problem using QR factorization; let me know if it's still less than clear! – Christian Clason May 25 '16 at 19:33

The intuition behind all Krylov's subspace methods is the following. Given a square matrix $A$ and a compatible vector $b$, there exists a unique monic polynomial $p$ such that $p(A)b = 0$. Mathematically, this is a consequence of the axiom of choice and Cayley's theorem which ensures that $q(A) = 0$ where $q$ is the characteristic polynomial of $A$. If $A$ is nonsingular, then $p(A)b=0$ can be transformed in $A w(A)b = b$ where $w$ is a polynomial closely related to $p$. In short, the solution $x = A^{-1}b$ can be written in the form $x = w(A)b$. Typically, $w$ has degree $n-1$, and so it is not suitable for practical computations. But this expression for $x$ is the reason why we seek approximations of similiar type, i.e. $x_k \in K_k(A,b)$.

• Why the downvote? I provided the information which is not in Saad's book which the OP has read. – Carl Christian May 27 '16 at 12:49
• Very useful comment, I think this is the most imptortant idea about Krylov subspace. By the way, do you have the matrices use in Saad's examples? – sunshine Dec 14 '19 at 1:15