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?
gmresfunction finds the exact solution for your example at iteration 3 with residual 0.) $\endgroup$