2
$\begingroup$

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{equation} \begin{bmatrix} 1 & -0.25 & 0\\ 0 & 1 & -0.25\\ 0 & 0 & 1 \end{bmatrix} \vec x = \begin{bmatrix} 1\\ 1\\ 1 \end{bmatrix} \end{equation}

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

\begin{equation} \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} \end{equation}

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?

$\endgroup$
  • $\begingroup$ 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.) $\endgroup$ – Christian Clason May 4 '13 at 8:57
  • $\begingroup$ 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? $\endgroup$ – Wolfgang Bangerth May 4 '13 at 13:19
  • $\begingroup$ @WolfgangBangerth Oops, I meant after the first "complete" run of GMRES, so there are N=3 orthogonal vectors computed (updated question to clarify this). $\endgroup$ – helloworld922 May 4 '13 at 18:47
  • $\begingroup$ @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. $\endgroup$ – helloworld922 May 4 '13 at 18:49
  • 3
    $\begingroup$ 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. $\endgroup$ – Christian Clason May 4 '13 at 20:53
1
$\begingroup$

To summarize Christian's comments and some digging through my code, my understanding of GMRES was correct, I had a code issue involving just the QR decomposition so it was a rather odd situation that just happened to converge with more iterations.

$\endgroup$
  • 2
    $\begingroup$ Glad it's sorted out! Since this question is not likely to be of help for future readers, I'm proposing to close it to keep it from being bumped. $\endgroup$ – Christian Clason Jun 4 '13 at 7:37
  • $\begingroup$ ... and it has now been bumped by Community because there is no accepted answer. Added vote to close. $\endgroup$ – Jed Brown Aug 4 '13 at 15:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.