Right preconditioned GMRES for arbitrary singular systems
Kota Sugihara, Ken Hayami
TL;DR
This work addresses solving linear systems $A x=b$ where $A$ may be singular and non-range-symmetric, conditions under which GMRES can break down. By applying GMRES to the right-preconditioned system $A C A^{T} z=b$ with SPD $C$ (i.e., $B=C A^{T}$), the authors obtain a least-squares solution $x=C A^{T} z$ for arbitrary $A$ and $b$, while using a pseudoinverse with a threshold and reorthogonalization to stabilize severely ill-conditioned Hessenberg systems. Numerical experiments on GP and index-2 matrices show that the proposed AB-GMRES approach with pseudoinverse outperforms standard GMRES and left-preconditioned variants in both inconsistent and consistent scenarios, highlighting improved robustness and convergence. The study also discusses practical considerations, including the cost of pseudoinverse computations and the potential for more efficient preconditioners to enhance applicability in larger problems.
Abstract
Brown and Walker (1997) showed that GMRES determines a least squares solution of $ A x = b $ where $ A \in {\bf R}^{n \times n} $ without breakdown for arbitrary $ b, x_0 \in {\bf R}^n $ if and only if $A$ is range-symmetric, i.e. $ {\cal R} (A^{\rm T}) = {\cal R} (A) $, where $ A $ may be singular and $ b $ may not be in the range space ${\cal R} A)$ of $A$. In this paper, we propose applying GMRES to $ A C A^{\rm T} z = b $, where $ C \in {\bf R}^{n \times n} $ is symmetric positive definite. This determines a least squares solution $ x = CA^{\rm T} z $ of $ A x = b $ without breakdown for arbitrary (singular) matrix $A \in {\bf R}^{n \times n}$ and $ b \in {\bf R}^n $. To make the method numerically stable, we propose using the pseudoinverse with an appropriate threshold parameter to suppress the influence of tiny singular values when solving the severely ill-conditioned Hessenberg systems which arise in the Arnoldi process of GMRES when solving inconsistent range-symmetric systems. Numerical experiments show that the method taking $C$ to be the identity matrix and the inverse matrix of the diagonal matrix whose diagonal elements are the diagonal of $A A^{\rm T}$ gives a least squares solution even when $A$ is not range-symmetric, including the case when $ {\rm index}(A) >1$.
