Tags:Arnoldi expansion, backward stability, Krylov subspace methods, loss of orthogonality and low-synchronization algorithms
Abstract:
The GMRES algorithm of Saad and Schultz (1986) for nonsymmetric linear systems relies on the Arnoldi expansion for the Krylov basis. The algorithm computes the $QR$ factorization of the matrix $B = [\: \vect{r}_0, AV_m\:]$. Despite an ${\cal O}(\eps)\kappa(B)$ loss of orthogonality, the modified Gram-Schmidt (MGS) formulation was shown to be backward stable in the seminal papers by Paige, et al. (2006) and Paige and Strako\v{s} (2002). Classical Gram-Schmidt (CGS) exhibits an ${\cal O}(\eps)\kappa^2(B)$ loss of orthogonality, whereas DCGS-2 (CGS with delayed reorthogonalization) reduces this to ${\cal O}(\eps)$ in practice (without a formal proof). We present a post-modern (viz not classical) GMRES algorithm based on Ruhe (1983) and the low-synch algorithms of Swirydowicz et al (2020) that achieves ${\cal O}(\eps) \|Av_k \|_2 / h_{k+1,k}$ loss of orthogonality. By projecting the vector $Av_m$ with Gauss-Seidel onto the orthogonal complement of the space spanned by the computed Krylov vectors $\bv_m$ where $\bv_m^T\bv_m = I + L_m + L_m^T$, we can further demonstrate that the loss of orthogonality closely follows ${\cal O}(\eps)$. For a broad class of matrices, unlike MGS-GMRES, significant loss of orthogonality does not occur and the relative residual no longer stagnates for highly non-normal systems. The Krylov vectors remain linearly independent and the smallest singular value of $\tv_m$ is close to one. We also demonstrate that Henrici's departure from normality of the lower triangular matrix $T_m \approx (\:V_m^TV_m\:)^{-1}$ in the Gram-Schmidt projector $P = I - V_mT_mV_m^T$ is an appropriate quantity for detecting the loss of orthogonality.