# Convergence of GMRES

From what I understand the GMRES method is (using Arnoldi Iterations/Modified Gram-Schmidt):

The first vector of the Krylov subspace span of A is the normalized vector $\frac{\vec b - A\vec x_0} {|| \vec b - A\vec x_0 ||}$

At each iteration i, calculate a single new orthonormal vector of the existing Krylov subspace. If the norm of that vector is 0 (or close to 0 in floating-precision calculations), then the subspace is "complete", and we should be able to calculate the converged solution. In theory, this occurs when we have N orthonormal vectors (where N is the size of A).

The solution can then be computed using $x = x_0 + V y$, where V is the orthonormal vectors computed above, and y minimizes $||H y - \beta e_0||$.

$\beta = ||\vec b - A\vec x_0||$

$e_0 = (1, 0, ..., 0)$

This can be solved in a few different ways, (say, the way Saad and Schultz proposed in their original paper using Givens rotations/QR decomposition interleaved with the Arnoldi iteration), and in theory x is converged to the solution.

Now I tried applying this to a sample test dataset:

$$\begin{bmatrix} 1 & -0.25 & 0\\ 0 & 1 & -0.25\\ 0 & 0 & 1 \end{bmatrix} \vec x = \begin{bmatrix} 1\\ 1\\ 1 \end{bmatrix}$$

Starting with $\vec x_0$ as a vector of 0's, I calculated 3 valid orthonormal vectors (calculated with double precision floats):

$$\begin{bmatrix} 0.57735026918962584\\ 0.57735026918962584\\ 0.57735026918962584 \end{bmatrix} \begin{bmatrix} -0.40824829046386368\\ -0.40824829046386368\\ 0.81649658092772548 \end{bmatrix} \begin{bmatrix} 0.70710678118655124\\ -0.7071067811865438\\ 2.5639502485114192e-15 \end{bmatrix}$$

As expected, the 4th vector calculated has a norm close to 0, and is not orthogonal to the other vectors in the Krylov subspace.

However, after the first "complete" run of the GMRES algorithm (so I have N=3 orthonormal vectors to work with), I calculated x as:

\begin{bmatrix} 1.3075447066949195\\ 1.2455012338896092\\ 1.0068140690151826 \end{bmatrix}

For reference, the "true" solution is:

\begin{bmatrix} 1.3125\\ 1.25\\ 1 \end{bmatrix}

This seems too large of a difference for me to be able to say that the difference is a result of floating point arithmetic.

Assuming the problem has converged (or stagnated), I would imagine re-running the GMRES would not significantly change the solution. However, running GMRES again would improve the solution significantly up to 3 or 4 times. Repeating the test for similar problems with different dimensions different dimensions and again about re-runs of GMRES improved the solution till it was closer to a small multiple of machine epsilon of the true solution (actual number of times varies slightly, no where near the order of the problem size).

Am I misunderstanding what convergence means in terms of the GMRES algorithm, or is there something else I'm missing in my implementation? Or is this really a result of imprecise arithmetic and extra iterations of the GMRES algorithm just happen to improve to solution?

• Without seeing your implementation, it's hard to guess what is going on. (I suspect inexact solution of the projected least-squares system.) Since the Arnoldi process in every iteration orthogonalizes against all previous vectors, it doesn't suffer from numerical loss of orthogonality like the CG method does. (And in fact, Matlab's gmres function finds the exact solution for your example at iteration 3 with residual 0.) – Christian Clason May 4 '13 at 8:57
• When you say "However, after the first iteration..." do you mean (i) after the first iteration of the GMRES method, or (ii) after running the entire GMRES method once for all $N=3$ steps? – Wolfgang Bangerth May 4 '13 at 13:19
• @WolfgangBangerth Oops, I meant after the first "complete" run of GMRES, so there are N=3 orthogonal vectors computed (updated question to clarify this). – helloworld922 May 4 '13 at 18:47
• @ChristianClason I'm not expecting anyone to debug the code, I was just curious as to what is typically considered an iteration in GMRES, and what convergence means in terms of GMRES. I can search through my code on my own time, and ask any additional questions I find in the future. – helloworld922 May 4 '13 at 18:49
• I wasn't inferring that, just pointing out that you are not misunderstanding what convergence means (GMRES should have found the exact solution, even - up to a very small error - in floating point arithmetic), and that it is most likely that there is indeed something missing (or incorrect) in your implementation. – Christian Clason May 4 '13 at 20:53