Iterated Gauss-Seidel GMRES

The MGS-GMRES algorithm of Saad and Schultz (1986) is an iterative method for approximately solving linear systems $A{\bf{x}}={\bf{b}}$, with initial guess ${\bf{x}}_0$ and residual ${\bf{r}}_0 = {\bf{b}} - A{\bf{x}}_0$. The algorithm employs the Arnoldi expansion of the Krylov basis vectors (columns of $V_k$). Paige and Strako\'{s} (2002) observe that this is equivalent to the $QR$ factorization of the matrix $B = [\: {\bf{r}}_0, AV_k\:]$ at each iteration. Despite an ${\cal O}(\epsilon)\kappa(B)$ loss of orthogonality, the modified Gram-Schmidt (MGS) formulation was shown to be backward stable in the seminal paper by Paige, et al. (2006). We present an iterated Gauss-Seidel formulation of the GMRES algorithm based on Ruhe (1983) and \'{S}wirydowicz et al. (2020) that achieves ${\cal O}(\epsilon) \: \|A{\bf{{v}}}_k\|_2 /h_{k+1,k}$ loss of orthogonality. By projecting the vector $A{\bf{v}}_k$ onto the orthogonal complement of the space spanned by the properly normalized Krylov vectors $\tilde{V}_k$ where $\tilde{V}_k^T\tilde{V}_k = I + L_k + L_k^T$, the loss of orthogonality is at most ${\cal O}(\epsilon)\kappa(B)$. For a broad class of matrices, a significant loss of orthogonality does not occur and the Arnoldi relative residual $\|\rho\:{\bf{e}}_1 - H_{k+1,k}{\bf{y}}_k\|_2/ \rho$, $\rho = \|{\bf{r}}_0\|_2$, no longer stagnates above machine precision for highly non-normal systems, where ${\bf{x}} = {\bf{x}}_0 + V_k{\bf{y}}_k$ is the approximate solution. The Krylov vectors remain linearly independent and the smallest singular value of $V_k$ remains close to one. We also demonstrate that Henrici's departure from normality of the matrix $T_k \approx (\:V_k^TV_k\:)^{-1}$ in the approximate projector $P = I - V_kT_kV_k^T$ is appropriate for detecting the loss of orthogonality.