On GMRES for singular EP and GP systems

TL;DR

This paper investigates the numerical behavior of GMRES for singular systems, identifying factors that cause accuracy deterioration and comparing it with RR-GMRES, supported by numerical experiments.

Contribution

It reveals how inconsistency, initial residuals, and principal angles impact GMRES accuracy for singular systems and compares GMRES with RR-GMRES.

Findings

GMRES accuracy deteriorates due to system inconsistency

Initial residuals' proximity to nullspace affects solution quality

Principal angles between ranges influence conditioning and accuracy

Abstract

In this contribution, we study the numerical behavior of the Generalized Minimal Residual (GMRES) method for solving singular linear systems. It is known that GMRES determines a least squares solution without breakdown if the coefficient matrix is range-symmetric (EP), or if its range and nullspace are disjoint (GP) and the system is consistent. We show that the accuracy of GMRES iterates may deteriorate in practice due to three distinct factors: (i) the inconsistency of the linear system; (ii) the distance of the initial residual to the nullspace of the coefficient matrix; (iii) the extremal principal angles between the ranges of the coefficient matrix and its transpose. These factors lead to poor conditioning of the extended Hessenberg matrix in the Arnoldi decomposition and affect the accuracy of the computed least squares solution. We also compare GMRES with the range restricted GMRES (RR-GMRES) method. Numerical experiments show typical behaviors of GMRES for small problems with EP and GP matrices.

Figures

Equations

$$A\mathbf{x}=\mathbf{b},$$

$$\Vert \mathbf{b}-A\mathbf{x}\Vert =\underset{\mathbf{u}\in {\mathbb{R}}^n}{\min }\Vert \mathbf{b}-A\mathbf{u}\Vert ,$$

$$\Vert \mathbf{b}-A{\mathbf{x}}_k\Vert =\underset{\mathbf{z}\in {\mathcal{K}}_k}{\min }\Vert \mathbf{b}-A({\mathbf{x}}_0+\mathbf{z})\Vert =\underset{\mathbf{z}\in {\mathcal{K}}_k}{\min }\Vert {\mathbf{r}}_0-A\mathbf{z}\Vert =\Vert {\mathbf{r}}_0-A{\mathbf{z}}_k\Vert$$

$$\Vert \mathbf{b}-A{\mathbf{x}}_k^{\mathrm{R}}\Vert$$

$$\Vert \mathbf{b}-A{\mathbf{x}}_k^{\mathrm{R}}{\Vert }^2$$

$$=\Vert (Q_{k+1}^{\mathrm{R}}{)}^{\mathsf{T}}{\mathbf{r}}_0-H_{k+1,k}^{\mathrm{R}}{\mathbf{y}}_k^{\mathrm{R}}{\Vert }^2+\Vert [\mathrm{I}-Q_{k+1}^{\mathrm{R}}(Q_{k+1}^{\mathrm{R}}{)}^{\mathsf{T}}]{\mathbf{r}}_0{\Vert }^2$$

$$=\underset{\mathbf{y}\in {\mathbb{R}}^k}{\min }\Vert (Q_{k+1}^{\mathrm{R}}{)}^{\mathsf{T}}{\mathbf{r}}_0-H_{k+1,k}^{\mathrm{R}}\mathbf{y}{\Vert }^2+\Vert [\mathrm{I}-Q_{k+1}^{\mathrm{R}}(Q_{k+1}^{\mathrm{R}}{)}^{\mathsf{T}}]{\mathbf{r}}_0{\Vert }^2.$$

$${\sigma }_k(H_{k+1,k})\le \underset{\mathbf{y}\in {\mathbb{R}}^{k-1}\backslash \{0\}}{\min }\frac{\Vert A{Q}_{k-1}^{\mathrm{R}}\mathbf{y}\Vert }{\Vert Q_{k-1}^{\mathrm{R}}\mathbf{y}\Vert }={\sigma }_{k-1}(H_{k,k-1}^{\mathrm{R}})$$

$$[U_1,V_2{]}^{-1}=\left[\begin{array}{cc}(V_1^{\mathsf{T}}{U}_1{)}^{-1}& \mathrm{O}\\ -V_2{U}_1(V_1^{\mathsf{T}}{U}_1{)}^{-1}& \mathrm{I}\end{array}\right][V_1,V_2{]}^{\mathsf{T}}.$$

$$A=U\left[\begin{array}{cc}{\Sigma }_{r}K& {\Sigma }_{r}L\\ \mathrm{O}& \mathrm{O}\end{array}\right]U^{\mathsf{T}}=V\left[\begin{array}{cc}K{\Sigma }_r& \mathrm{O}\\ M{\Sigma }_r& \mathrm{O}\end{array}\right]V^{\mathsf{T}}$$

$$A^{†}=U\left[\begin{array}{cc}{K}^{\mathsf{T}}{\Sigma }_r^{-1}& \mathrm{O}\\ L^{\mathsf{T}}{\Sigma }_r^{-1}& \mathrm{O}\end{array}\right]U^{\mathsf{T}}=V\left[\begin{array}{cc}{\Sigma }_r^{-1}{K}^{\mathsf{T}}& {\Sigma }_r^{-1}{M}^{\mathsf{T}}\\ \mathrm{O}& \mathrm{O}\end{array}\right]V^{\mathsf{T}}.$$

$$\displaystyle\mathcal{R}(U_{1})\cap\mathcal{R}(V_{2})=\{\boldsymbol{0}\}\Longleftrightarrow\operatorname{rank}([U_{1},V_{2}])=n\Longleftrightarrow V_{1}^{\mathsf{T}}\!U_{1}\mbox{($=K$) is nonsingular}$$

$$[V_1,V_2{]}^{\mathsf{T}}[U_1,V_2]=\left[\begin{array}{cc}{V}_1^{\mathsf{T}}{U}_1& \mathrm{O}\\ V_2^{\mathsf{T}}{U}_1& \mathrm{I}\end{array}\right].$$

$$A^{\#}=U\left[\begin{array}{cc}{K}^{-1}{\Sigma }_r^{-1}& K^{-1}{\Sigma }_r^{-1}{K}^{-1}L\\ \mathrm{O}& \mathrm{O}\end{array}\right]U^{\mathsf{T}}=V\left[\begin{array}{cc}{\Sigma }_r^{-1}{K}^{-1}& \mathrm{O}\\ M{K}^{-1}{\Sigma }_r^{-1}{K}^{-1}& \mathrm{O}\end{array}\right]V^{\mathsf{T}}.$$

$$\displaystyle\kappa(A\rvert_{\mathcal{K}_{k}})=\frac{\max_{\boldsymbol{z}\in\mathcal{K}_{k}\backslash\{\boldsymbol{0}\}}\|A\boldsymbol{z}\|\mbox{\Big{/}}\|\boldsymbol{z}\|}{\min_{\boldsymbol{z}\in\mathcal{K}_{k}\backslash\{\boldsymbol{0}\}}\|A\boldsymbol{z}\|\mbox{\Big{/}}\|\boldsymbol{z}\|}$$

$$\{\underset{\mathbf{z}\in {\mathcal{K}}_k\backslash \{0\}}{\max ,\min }\}\frac{\Vert A\mathbf{z}\Vert }{\Vert \mathbf{z}\Vert }=\{\underset{\mathbf{w}\in {\mathbb{R}}^k\backslash \{0\}}{\max ,\min }\}\frac{\Vert A{Q}_k\mathbf{w}\Vert }{\Vert Q_k\mathbf{w}\Vert }=\{\underset{\mathbf{w}\in {\mathbb{R}}^k\backslash \{0\}}{\max ,\min }\}\frac{\Vert H_{k+1,k}\mathbf{w}\Vert }{\Vert \mathbf{w}\Vert }.$$

$$\Vert H_{k+1,k}\Vert =\underset{\mathbf{z}\in {\mathcal{K}}_k\backslash \{0\}}{\max }\frac{\Vert A\mathbf{z}\Vert }{\Vert \mathbf{z}\Vert }\le \underset{\mathbf{z}\in \mathrm{span}\{{\mathbf{r}}_{\ast }\}\oplus \mathcal{R}(A)\backslash \{0\}}{\max }\frac{\Vert A\mathbf{z}\Vert }{\Vert \mathbf{z}\Vert }\le \underset{\mathbf{z}\in {\mathbb{R}}^n\backslash \{0\}}{\max }\frac{\Vert A\mathbf{z}\Vert }{\Vert \mathbf{z}\Vert }=\Vert A\Vert .$$

$${\sigma }_k(H_{k+1,k})=\underset{\mathbf{z}\in {\mathcal{K}}_k\backslash \{0\}}{\min }\frac{\Vert A\mathbf{z}\Vert }{\Vert \mathbf{z}\Vert }\ge \underset{\mathbf{z}\in \mathrm{span}\{{\mathbf{r}}_{\ast }\}\oplus \mathcal{R}(A)\backslash \{0\}}{\min }\frac{\Vert A\mathbf{z}\Vert }{\Vert \mathbf{z}\Vert }\ge \underset{z\in {\mathbb{R}}^n\backslash \{0\}}{\min }\frac{\Vert A\mathbf{z}\Vert }{\Vert \mathbf{z}\Vert }=0$$

$${\sigma }_k(H_{k+1,k})\ge \underset{\mathbf{z}\in \mathcal{R}(A)\backslash \{0\}}{\min }\frac{\Vert A\mathbf{z}\Vert }{\Vert \mathbf{z}\Vert },$$

$$\underset{\mathbf{z}\in \mathcal{R}(A)\backslash \{0\}}{\min }\frac{\Vert A\mathbf{z}\Vert }{\Vert \mathbf{z}\Vert }=\underset{\mathbf{z}\in \mathcal{R}(A^{\mathsf{T}})\backslash \{0\}}{\min }\frac{\Vert A\mathbf{z}\Vert }{\Vert \mathbf{z}\Vert }={\sigma }_r(A)>0$$

$$\kappa (A{|}_{\mathcal{R}(A)})=\frac{\Vert A\Vert }{{\min }_{\mathbf{z}\in \mathcal{R}(A)\backslash \{0\}}\Vert A\mathbf{z}\Vert /\Vert \mathbf{z}\Vert }=\kappa (A).$$

$${\sigma }_k(H_{k+1,k})\ge \underset{\mathbf{z}\in \mathrm{span}\{{\mathbf{r}}_{\ast }\}\oplus \mathcal{R}(A)\backslash \{0\}}{\min }\frac{\Vert A\mathbf{z}\Vert }{\Vert \mathbf{z}\Vert }=0.$$

$${\sigma }_k(H_{k+1,k})$$

$$=\frac{\Vert A({\mathbf{r}}_{k-1}-{\mathbf{r}}_{\ast })\Vert }{\Vert {\mathbf{r}}_{k-1}\Vert }\le \Vert A\Vert \frac{\Vert {\mathbf{r}}_{k-1}-{\mathbf{r}}_{\ast }\Vert }{\Vert {\mathbf{r}}_{k-1}\Vert }.$$

$${\sigma }_k(H_{k+1,k})=\underset{\mathbf{z}\in {\mathcal{K}}_k\backslash \{0\}}{\min }\frac{\Vert A\mathbf{z}\Vert }{\Vert \mathbf{z}\Vert }\le \frac{\Vert A{\mathbf{r}}_0\Vert }{\Vert {\mathbf{r}}_0\Vert }\le \frac{\Vert A\Vert \Vert {\mathbf{r}}_0{|}_{\mathcal{R}(A^{\mathsf{T}})}\Vert }{\Vert {\mathbf{r}}_0\Vert }=\frac{\Vert A\Vert \Vert {\mathbf{r}}_0{|}_{\mathcal{R}(A)}\Vert }{\Vert {\mathbf{r}}_0\Vert }.$$

$$A=\left[\begin{array}{cc}1& 0\\ 0& 0\end{array}\right],\mathbf{b}=\left[\begin{array}{c}1\\ \epsilon \end{array}\right],\epsilon >0.$$

$$Q_2=[{\mathbf{q}}_1,{\mathbf{q}}_2]=\frac{1}{\sqrt{1+{\epsilon }^2}}\left[\begin{array}{cc}1& \epsilon \\ \epsilon & -1\end{array}\right],H_{2,2}=\frac{1}{1+{\epsilon }^2}\left[\begin{array}{cc}1& \epsilon \\ \epsilon & {\epsilon }^2\end{array}\right].$$

$$A=\left[\begin{array}{cc}1& 0\\ 0& \delta \end{array}\right],\mathbf{b}=\left[\begin{array}{c}1\\ \epsilon \end{array}\right],\delta ,\epsilon >0$$

$$H_{2,2}=\frac{1}{1+{\epsilon }^2}\left[\begin{array}{cc}1+\delta {\epsilon }^2& (1-\delta )\epsilon \\ (1-\delta )\epsilon & \delta +{\epsilon }^2\end{array}\right],$$

$$A=\left[\begin{array}{cc}D& 0\\ 0& 0\end{array}\right]\in {\mathbb{R}}^{128\times 128},\mathbf{b}=\left[\begin{array}{c}\mathbf{\gamma }\\ \mathbf{\delta }\end{array}\right],$$

$${\sigma }_k(H_{k+1,k}^{\mathrm{R}})$$

$$\ge \underset{\mathbf{z}\in \mathcal{R}(A)\backslash \{0\}}{\min }\frac{\Vert A\mathbf{z}\Vert }{\Vert \mathbf{z}\Vert }=\underset{\mathbf{z}\in \mathcal{R}(A^{\mathsf{T}})\backslash \{0\}}{\min }\frac{\Vert A\mathbf{z}\Vert }{\Vert \mathbf{z}\Vert }={\sigma }_r(A)$$

Full text

1 Introduction

Consider solving linear systems of equations

[TABLE]

where $A\in\mathbb{R}^{n\times n}$ may be singular and $\boldsymbol{b}\in\mathbb{R}^{n}$ is not necessarily in $\mathcal{R}(A)=\{\boldsymbol{y}\in\mathbb{R}^{n}\mid\boldsymbol{y}=A\boldsymbol{x},\ \boldsymbol{x}\in\mathbb{R}^{n}\}$, the range of $A$. We say that $A\boldsymbol{x}=\boldsymbol{b}$ is consistent if $\boldsymbol{b}\in\mathcal{R}(A)$, and otherwise it is inconsistent. If eq. 1.1 is inconsistent, instead of eq. 1.1, it is natural to consider solving the least squares problem

[TABLE]

where $\|\cdot\|$ denotes the Euclidean norm. We call a minimizer $\boldsymbol{x}\in\operatorname*{argmin}_{\boldsymbol{u}\in\mathbb{R}^{n}}\|\boldsymbol{b}-A\boldsymbol{u}\|$ a least squares solution, which is not necessarily unique.

In order to analyze iterative methods for solving eq. 1.1 in terms of the spaces associated with $A$, we give some required definitions and notations. Let $\mathcal{N}(A)=\{\boldsymbol{x}\in\mathbb{R}^{n}\mid A\boldsymbol{x}=\boldsymbol{0}\}$ be the nullspace of $A$. Then, we have $\mathcal{N}(A^{\mathsf{T}})\oplus\mathcal{R}(A)=\mathcal{N}(A)\oplus\mathcal{R}(A^{\mathsf{T}})=\mathbb{R}^{n}$, $\dim\mathcal{N}(A^{\mathsf{T}})=\dim\mathcal{N}(A)$, and $\dim\mathcal{R}(A^{\mathsf{T}})=\dim\mathcal{R}(A)=\operatorname{rank}(A)$, where $\oplus$ denotes the direct sum of subspaces. Let $r=\operatorname{rank}(A)$ and denote the singular value decomposition (SVD) of $A$ by $U\Sigma V^{\mathsf{T}}$, where $U\in\mathbb{R}^{n\times n}$ and $V\in\mathbb{R}^{n\times n}$ are orthogonal matrices $U^{\mathsf{T}}\!U=UU^{\mathsf{T}}=V^{\mathsf{T}}V=VV^{\mathsf{T}}=\mathrm{I}$, $\mathrm{I}$ is the identity matrix, $\Sigma=\operatorname{diag}(\sigma_{1},\sigma_{2},\dots,\sigma_{r},0,0,\dots,0)\in\mathbb{R}^{n\times n}$, and $\sigma_{i}$ is the $i$th largest nonzero singular value of $A$. Let $U=[U_{1},U_{2}]$ and $V=[V_{1},V_{2}]$, where the columns of $U_{1}\in\mathbb{R}^{n\times r}$ and $U_{2}\in\mathbb{R}^{n\times(n-r)}$ form orthonormal bases of $\mathcal{R}(A)=\mathcal{R}(U_{1})$ and $\mathcal{R}(A)^{\perp}=\mathcal{N}(A^{\mathsf{T}})=\mathcal{R}(U_{2})$, respectively, and the columns of $V_{1}\in\mathbb{R}^{n\times r}$ and $V_{2}\in\mathbb{R}^{n\times(n-r)}$ form orthonormal bases of $\mathcal{N}(A)^{\perp}=\mathcal{R}(A^{\mathsf{T}})=\mathcal{R}(V_{1})$ and $\mathcal{N}(A)=\mathcal{R}(V_{2})$, respectively.

We recall the definitions of generalized inverses. We call a matrix $X\in\mathbb{R}^{n\times n}$ the Moore-Penrose generalized inverse (pseudoinverse) of $A\in\mathbb{R}^{n\times n}$ if $X$ satisfies the Penrose equations $AXA=A$, $XAX=X$, $(AX)^{\mathsf{T}}=AX$, and $(XA)^{\mathsf{T}}=XA$, denote it by $A^{\dagger}$, and have the identity $A^{\dagger}=V_{1}\Sigma_{r}^{-1}{U_{1}}^{\mathsf{T}}$, where $\Sigma_{r}=\operatorname{diag}(\sigma_{1},\sigma_{2},\dots,\sigma_{r})$. The condition number of $A$ is denoted by $\kappa(A)=\|A\|\|A^{\dagger}\|$ [5, Definition 1.4.2]. The smallest nonnegative integer $k$ such that $\operatorname{rank}(A^{k})=\operatorname{rank}(A^{k+1})$ is called the index of $A$ [10, Definition 7.2.1], and is denoted by $\operatorname{index}(A)$. In addition, $k\geq\mathrm{index}(A)\iff\mathcal{N}(A^{k})\oplus\mathcal{R}(A^{k})=\mathbb{R}^{n}$ [10, p. 121]. Let $\operatorname{index}(A)=1$ and $X\in\mathbb{R}^{n\times n}$ be such that $AXA=A$, $XAX=X$, and $AX=XA$. Then, $X$ is unique, and we call $X$ the group inverse of $A$ and denote it by $A^{\#}$. The group inverse of the matrix $A$ is on $\mathcal{R}(A)$ equal to the inverse of the restriction of $A$ to its range ${\cal R}(A)$ and admits $\mathcal{N}(A)$ as its nullspace [21, Theorem 2]. The group inverse can be characterized by the Jordan canonical form. Let $S$ be a nonsingular matrix such that $J=S^{-1}\!AS$ is the Jordan canonical form of $A$. Then, $A^{\#}=SJ^{\dagger}S^{-1}$ holds. In particular, we have $\mathcal{R}(A^{\#})=\mathcal{R}(A)$, $\mathcal{N}(A^{\#})=\mathcal{N}(A)$, and $A^{\#}\!A=AA^{\#}=P_{\mathcal{R}(A),\mathcal{N}(A)}$ [4]. Here, $P_{\mathcal{R}(A),\mathcal{N}(A)}$ denotes the projection onto $\mathcal{R}(A)$ along $\mathcal{N}(A)$.

Now, we express solutions of eqs. 1.1 and 1.2. The vector $\boldsymbol{x}_{*}=A^{\dagger}\boldsymbol{b}=V_{1}\Sigma_{r}^{-1}U_{1}^{\mathsf{T}}\boldsymbol{b}$ is called the minimum-norm least squares or pseudoinverse solution of eq. 1.1 or eq. 1.2, and it belongs to $\mathcal{R}(A^{\mathsf{T}})$. We next give the expressions of the residual of eq. 1.2. Denote the orthogonal projector onto $\mathcal{R}(A)$ by $P_{\mathcal{R}(A)}=U_{1}U_{1}^{\mathsf{T}}$ and that onto $\mathcal{N}(A^{\mathsf{T}})$ by $P_{\mathcal{N}(A^{\mathsf{T}})}=U_{2}U_{2}^{\mathsf{T}}$. If $\boldsymbol{v}\rvert_{\mathcal{S}}\in\mathbb{R}^{n}$ is the orthogonal projection of a vector $\boldsymbol{v}\in\mathbb{R}^{n}$ onto the subspace $\mathcal{S}\subseteq\mathbb{R}^{n}$, then for any $\boldsymbol{x}_{0}\in\mathbb{R}^{n}$ and any $\boldsymbol{b}\in\mathbb{R}^{n}$, the corresponding residual of eq. 1.2 is $\boldsymbol{r}_{0}=\boldsymbol{b}-A\boldsymbol{x}_{0}=\boldsymbol{b}\rvert_{\mathcal{N}(A^{\mathsf{T}})}+\boldsymbol{b}\rvert_{\mathcal{R}(A)}-A\boldsymbol{x}_{0}=\boldsymbol{r}_{*}+\boldsymbol{r}_{0}\rvert_{\mathcal{R}(A)}$, where $\boldsymbol{r}_{*}=\boldsymbol{b}\rvert_{\mathcal{N}(A^{\mathsf{T}})}=P_{\mathcal{N}(A^{\mathsf{T}})}\boldsymbol{b}\in\mathcal{N}(A^{\mathsf{T}})$ is the least squares residual and $\boldsymbol{b}\rvert_{\mathcal{R}(A)}=P_{\mathcal{R}(A)}\boldsymbol{b}$. In particular, we have $\boldsymbol{r}_{*}=\boldsymbol{0}$ for $\boldsymbol{b}\in\mathcal{R}(A)$. Assume $\operatorname{index}(A)=1$. It holds that $\boldsymbol{x}_{\#}=A^{\#}\boldsymbol{b}$ is a solution of $A\boldsymbol{x}=\boldsymbol{b}$ $\Longleftrightarrow$ $\boldsymbol{b}\in\mathcal{R}(A)$ [10], and it is the unique solution of $A\boldsymbol{x}=\boldsymbol{b}$ in $\mathcal{R}(A)$. Furthermore, the minimum Euclidean norm solution of eq. 1.2 satisfies $\boldsymbol{x}_{*}=P_{\mathcal{R(A^{\mathsf{T}})}}\boldsymbol{x}_{\#}$.

In this paper, we are interested in the numerical behavior of the Generalized Minimal Residual (GMRES) method [22] applied in particular to singular systems eq. 1.1. In section 2, we give some well-known conditions under which GMRES determines a solution without breakdown when applied to certain classes of singular matrices. We discuss also a relation to the range-restricted GMRES (RR-GMRES) method proposed in [7]. In section 3, we examine the conditioning of the coefficient matrix $A$ restricted to the Krylov subspaces that significantly influences the numerical behavior of GMRES. We consider first the case of EP (equal projection) or range-symmetric matrices and distinguish between the consistent and inconsistent cases showing that the consistent case is similar to the nonsingular case. Then we discuss the inconsistent EP case where GMRES suffers from an instability, since the convergence means ill-conditioned restriction of $A$ to the Krylov subspaces. In section 4 we study the case of group (GP) matrices. We show that the numerical behavior of GMRES applied to such problems depends substantially on the principal angles between the subspaces $\mathcal{R}(A)$ and $\mathcal{R}(A^{\mathsf{T}})$. Surprisingly, difficulties can be expected for non-range-symmetric problems even for consistent systems. In section 5, we conclude the paper.

2 GMRES methods and its convergence for singular systems

GMRES for the linear system eq. 1.1 with initial iterate $\boldsymbol{x}_{0}\in\mathbb{R}^{n}$, independent of any particular implementation of the algorithm, determines the $k$th iterate $\boldsymbol{x}_{k}$ over $\boldsymbol{x}_{0}+\mathcal{K}_{k}(A,\boldsymbol{r}_{0})$ that minimizes $\|\boldsymbol{b}-A\boldsymbol{x}_{k}\|$, where $\boldsymbol{r}_{0}=\boldsymbol{b}-A\boldsymbol{x}_{0}$ is the initial residual and $\mathcal{K}_{k}(A,\boldsymbol{r}_{0})=\mathrm{span}\{\boldsymbol{r}_{0},A\boldsymbol{r}_{0},\dots,A^{k-1}\boldsymbol{r}_{0}\}$ is the Krylov subspace of order $k$. Note that there exist $\boldsymbol{x}_{k}$ and hence $\boldsymbol{r}_{k}=\boldsymbol{b}-A\boldsymbol{x}_{k}$ for all $k\geq 0$ but $\boldsymbol{x}_{k}$ may not be unique in the singular case. Denote $\mathcal{K}_{k}(A,\boldsymbol{r}_{0})$ by $\mathcal{K}_{k}$ for simplicity. It is clear that $\mathcal{K}_{k}=\mathrm{span}\{\boldsymbol{r}_{0}\}+A\mathcal{K}_{k-1}\subset\mathrm{span}\{\boldsymbol{r}_{*}\}+\mathcal{R}(A)$ holds. If $\dim A\mathcal{K}_{k}=\dim\mathcal{K}_{k}$, then the problem

[TABLE]

has a unique solution $\boldsymbol{x}_{k}=\boldsymbol{x}_{0}+\boldsymbol{z}_{k}$ and hence $\boldsymbol{r}_{k}=\boldsymbol{b}-A\boldsymbol{x}_{k}\in\boldsymbol{r}_{0}+A\mathcal{K}_{k}$ is uniquely determined.

General studies on Krylov subspace methods in the singular case were done in [16], [24], [31], [23]. Particular studies on GMRES-type methods in the singular case were done in [7], [25], [20], [26], [13]. See [8], [9], [14] for GMRES on ill-posed linear systems, and [32] for GMRES with preconditioning. See also [15] for GMRES and [1] for GMRES with preconditioning in Hilbert spaces.

The GMRES method for solving ill-posed problems was studied in [8], [9], [14]. The preconditioned GMRES method applied to singular systems was considered in [32] . See also [15] for GMRES and [1] for the preconditioned GMRES method applied in Hilbert spaces.

In the nonsingular case, GMRES determines the solution of $A\boldsymbol{x}=\boldsymbol{b}$ for all $\boldsymbol{b}\in\mathbb{R}^{n}$ and for all $\boldsymbol{x}_{0}\in\mathbb{R}^{n}$ within $n$ iterations. In the singular case, GMRES may fail to determine a solution of eq. 1.1, and is said to break down at some step $k$ if $\dim A\mathcal{K}_{k}<\dim\mathcal{K}_{k}~{}\mbox{or}~{}\dim\mathcal{K}_{k}<k$ [6, p. 38]. Note that, in general, $\dim A\mathcal{K}_{k}\leq\dim\mathcal{K}_{k}\leq k$ holds for each $k$.

We give an explicit expression of the iterate $\boldsymbol{x}_{k}$ for GMRES using the Arnoldi decomposition $AQ_{k}=Q_{k+1}H_{k+1,k}$, $k=1,2,\dots$, where the columns of $Q_{k}=\left[\boldsymbol{q}_{1},\boldsymbol{q}_{2},\dots,\boldsymbol{q}_{k}\right]$ form an orthonormal basis of the Krylov subspace $\mathcal{K}_{k}$, and $H_{k+1,k}=(h_{i,j})\in\mathbb{R}^{(k+1)\times k}$ is an extended Hessenberg matrix. Then the iterate is given by $\boldsymbol{x}_{k}=\boldsymbol{x}_{0}+Q_{k}\boldsymbol{y}_{k}$ with $\boldsymbol{y}_{k}=\mathrm{arg\,min}_{\boldsymbol{y}\in\mathbb{R}^{k}}\|\beta\boldsymbol{e}_{1}-H_{k+1,k}\boldsymbol{y}\|$, where $\beta=\|\boldsymbol{r}_{0}\|$, $\boldsymbol{e}_{1}$ is the first column of the identity matrix and $\|\boldsymbol{b}-A\boldsymbol{x}_{k}\|=\|\boldsymbol{r}_{0}-AQ_{k}\boldsymbol{y}_{k}\|=\|\beta\boldsymbol{e}_{1}-H_{k+1,k}\boldsymbol{y}_{k}\|$.

It is clear that if $h_{i+1,i}\neq 0$ for $i=1,2,\dots,k-1$, the breakdown does not occur until step $k-1$ of GMRES with $\dim A\mathcal{K}_{i}=i$, or $\operatorname{rank}(H_{i+1,i})=i$, $i=1,2,\dots,k-1$. At breakdown of GMRES at step $k$ with $h_{k+1,k}=0$, one of the following cases holds [18, Appendix B] (cf. [6, Theorem 2.2]):

Case I.

$\dim A\mathcal{K}_{k+1}=k<\dim\mathcal{K}_{k+1}=k+1$, whereas $\operatorname{rank}(H_{k,k})=k-1$.

Case II.

$\dim A\mathcal{K}_{k}=k=\dim\mathcal{K}_{k+1}<k+1$, whereas $\operatorname{rank}(H_{k,k})=k$ (GMRES determines a solution of $A\boldsymbol{x}=\boldsymbol{b}$ at step $k$).

Here, $H_{k.k}=(h_{i,j})\in\mathbf{R}^{k\times k}$.

A variant of GMRES called the range restricted GMRES (RR-GMRES) method was proposed in [7]. RR-GMRES determines the $k$th iterate by minimizing the same objective function as GMRES over a different Krylov subspace

[TABLE]

It was shown in [11, Theorem A2] that if RR-GMRES applied to eq. 1.2 breaks down at step $m$ with $\operatorname{rank}(A)=m-1$ and $\dim A\mathcal{K}_{m}(A,A\boldsymbol{r}_{0})=m-1$, then it determines a solution of eq. 1.2. Here, RR-GMRES is said to break down if $\dim A\mathcal{K}_{k}(A,A\boldsymbol{r}_{0})<\dim\mathcal{K}_{k}(A,A\boldsymbol{r}_{0})$ or $\dim\mathcal{K}_{k}(A,A\boldsymbol{r}_{0})<k$.

We give an explicit expression of the RR-GMRES iterate $\boldsymbol{x}_{k}$ using the Arnoldi decomposition $AQ_{k}^{\mathrm{R}}=Q_{k+1}^{\mathrm{R}}H_{k+1,k}^{\mathrm{R}}$, $k=1,2,\dots$, where the columns of $Q_{k}^{\mathrm{R}}=[\boldsymbol{q}_{1}^{\mathrm{R}},\boldsymbol{q}_{2}^{\mathrm{R}},\dots,\boldsymbol{q}_{k}^{\mathrm{R}}]$ form an orthonormal basis of the Krylov subspace $\mathcal{K}_{k}(A,A\boldsymbol{r}_{0})$ with the initial vector $\boldsymbol{q}_{1}^{\mathrm{R}}=A\boldsymbol{r}_{0}/\|A\boldsymbol{r}_{0}\|$, and $H_{k+1,k}^{\mathrm{R}}=(h_{i,j}^{\mathrm{R}})\in\mathbb{R}^{(k+1)\times k}$ is an extended Hessenberg matrix. Then, the iterate is given by $\boldsymbol{x}_{k}^{\mathrm{R}}=\boldsymbol{x}_{0}+Q_{k}^{\mathrm{R}}\boldsymbol{y}_{k}^{\mathrm{R}}$ with $\boldsymbol{y}_{k}^{\mathrm{R}}=\mathrm{arg\,min}_{\boldsymbol{y}\in\mathbb{R}^{k}}\|(Q_{k+1}^{\mathrm{R}})^{\mathsf{T}}\boldsymbol{r}_{0}-H_{k+1,k}^{\mathrm{R}}\boldsymbol{y}\|$, where

[TABLE]

The last term is equal to the $k$th residual norm for the simpler GMRES method [29], which is not larger than the $k$th residual norm for RR-GMRES, i.e., $\|\boldsymbol{b}-A\boldsymbol{x}_{k}\|\leq\|\boldsymbol{b}-A\boldsymbol{x}_{k}^{\mathrm{R}}\|$. Note also that $\|H_{k,k-1}^{\mathrm{R}}\|=\|AQ_{k-1}^{\mathrm{R}}\|\leq\|AQ_{k}\|=\|H_{k+1,k}\|$ and

[TABLE]

leading to an interesting bound $\kappa(H_{k,k-1}^{\mathrm{R}})\leq\kappa(H_{k+1,k})$ for $k=2,\dots,n-1$, where $\sigma_{k}(\cdot)$ denotes the $k$th largest singular value of a matrix.

In the following, we present conditions under which GMRES determines a solution of eq. 1.1. We start with the observation that in the case of $\mathcal{N}(A)\cap\mathcal{R}(A)\neq\{\boldsymbol{0}\}$, GMRES breaks down and fails to determine a solution.

Proposition 2.1**.**

If $\boldsymbol{b}\in\mathcal{R}(A)$ and $\boldsymbol{0}\neq\boldsymbol{r}_{0}\in\mathcal{N}(A)\cap\mathcal{R}(A)$, then GMRES breaks down at step 1 without determining a solution of $A\boldsymbol{x}=\boldsymbol{b}$.

Proof.

Since $\boldsymbol{r}_{0}\neq\boldsymbol{0}$, we have $\dim\mathcal{K}_{1}=\dim\mathrm{span}\{\boldsymbol{r}_{0}\}=1$. Since $\boldsymbol{r}_{0}\in\mathcal{N}(A)$ gives $A\boldsymbol{r}_{0}=\boldsymbol{0}$, we have $\dim A\mathcal{K}_{1}=\dim\mathrm{span}\{A\boldsymbol{r}_{0}\}=\dim\mathrm{span}\{\boldsymbol{0}\}=0$. Hence, $\dim A\mathcal{K}_{1}<\dim\mathcal{K}_{1}$ holds. Therefore, GMRES breaks down at step 1 without determining a solution of $A\boldsymbol{x}=\boldsymbol{b}$. ∎

Similarly to GMRES, RR-GMRES also breaks down at step 1 without determining a solution of $A\boldsymbol{x}=\boldsymbol{b}$ if $\boldsymbol{b}\in\mathcal{R}(A)$ and $\boldsymbol{0}\neq\boldsymbol{r}_{0}\in\mathcal{N}(A)\cap\mathcal{R}(A)$. Therefore, we will restrict our attention to the cases of $\mathcal{N}(A)\cap\mathcal{R}(A)=\{\boldsymbol{0}\}$. The following statement holds.

Theorem 2.2** ([6, Theorem 2.6], [31, Theorem 3.2]).**

If $\mathcal{N}(A)\cap\mathcal{R}(A)=\{\boldsymbol{0}\}$, then GMRES determines a solution of $A\boldsymbol{x}=\boldsymbol{b}$ without breakdown for all $\boldsymbol{b}\in\mathcal{R}(A)$ and for all $\boldsymbol{x}_{0}\in\mathbb{R}^{n}$. The solution is $\boldsymbol{x}_{\#}+(\mathrm{I}-A^{\#}\!A)\boldsymbol{x}_{0}$.

The condition $\mathcal{N}(A)\cap\mathcal{R}(A)=\{\boldsymbol{0}\}$ is equivalent to $\operatorname{rank}([U_{1},V_{2}])=n$, or $\mathcal{R}(U_{1})\cap\mathcal{R}(V_{2})=\{\boldsymbol{0}\}$ [10, Lemma 7.2.1]. Then, we have $P_{\mathcal{R}(A),\mathcal{N}(A)}=A^{\#}\!A=ff[U_{1},\mathrm{O}][U_{1},V_{2}]^{-1}=\break U_{1}(V_{1}^{\mathsf{T}}U_{1})^{-1}V_{1}^{\mathsf{T}}$, because of [4, Exercise 30, p. 167] and

[TABLE]

Thus, $\kappa(A^{\#}\!A)=\kappa(V_{1}^{\mathsf{T}}U_{1})$ holds.

For a special class of singular matrices, GMRES determines a least squares solution.

Theorem 2.3** ([6, Theorem 2.4]).**

If $\mathcal{R}(A)=\mathcal{R}(A^{\mathsf{T}})$, then GMRES determines a solution of $\min_{\boldsymbol{x}\in\mathbb{R}^{n}}\|\boldsymbol{b}-A\boldsymbol{x}\|$ without breakdown for all $\boldsymbol{b}\in\mathbb{R}^{n}$ and for all $\boldsymbol{x}_{0}\in\mathbb{R}^{n}$.

A matrix $A\in\mathbb{R}^{n\times n}$ satisfying $\mathcal{N}(A)\cap\mathcal{R}(A)=\{\boldsymbol{0}\}$, is called a GP (group) matrix. A GP matrix satisfying in addition $\mathcal{R}(A^{\mathsf{T}})=\mathcal{R}(A)$, or equivalently $\mathcal{R}(U_{1})=\mathcal{R}(V_{1})$, is called an EP (equal projection) or range-symmetric matrix. Now, we characterize the GP and EP matrices in terms of their singular value decompositions. The matrix $A$ can be decomposed into

[TABLE]

with the identity $KK^{\mathsf{T}}+LL^{\mathsf{T}}=\mathrm{I}$, where $K=V_{1}^{\mathsf{T}}\!U_{2}$, $L=V_{1}^{\mathsf{T}}U_{1}$, and $M=V_{2}^{\mathsf{T}}U_{1}$, and the pseudoinverse matrix of $A$ can be decomposed into

[TABLE]

The equivalences for GP matrices

[TABLE]

follow from the equation

[TABLE]

The group inverse of a GP matrix $A$ can be decomposed into

[TABLE]

See [2, Section 1]. The conditioning of $K=V_{1}^{\mathsf{T}}\!U_{1}$ is independent of the conditioning of $A$ but it gives a difficulty in solving singular linear systems with GMRES. The EP case $\mathcal{R}(A^{\mathsf{T}})=\mathcal{R}(A)$ is equivalent to that $K=V_{1}^{\mathsf{T}}U_{1}$ is orthogonal, since $L=V_{1}^{\mathsf{T}}U_{2}=\mathrm{O}$.

Next, we characterize GP and EP matrices in terms of the principal angles. In the EP case $\mathcal{R}(A^{\mathsf{T}})=\mathcal{R}(A)$, the matrix $V_{1}^{\mathsf{T}}U_{1}$ is orthogonal and the cosines of the principal angles between $\mathcal{R}(A)$ and $\mathcal{R}(A^{\mathsf{T}})$ are all zero. In the GP case, since the columns of $U_{1}$ and $V_{1}$ form bases of $\mathcal{R}(A)$ and $\mathcal{R}(A^{\mathsf{T}})$, respectively, the cosines of the canonical angles between $\mathcal{R}(A)$ and $\mathcal{R}(A^{\mathsf{T}})$ are the singular values of $V_{1}^{\mathsf{T}}U_{1}$ [12, section 1.2]. Hence, the condition number of $V_{1}^{\mathsf{T}}\!U_{1}$ is related to the extremal principal angles.

Note that due to $\|V_{1}^{\mathsf{T}}U_{1}\|\leq 1$ all singular values of $V_{1}^{\mathsf{T}}U_{1}$ are less than or equal to 1 and the number of those equal exactly to 1 gives the dimension of $\mathcal{R}(U_{1})\cap\mathcal{R}(V_{1})$. So, if $\|V_{1}^{\mathsf{T}}U_{1}\|<1$, then $\mathcal{R}(U_{1})\cap\mathcal{R}(V_{1})=\mathcal{R}(A)\cap\mathcal{R}(A^{\mathsf{T}})=\{\boldsymbol{0}\}$. If a matrix $A\in\mathbb{R}^{n\times n}$ satisfies $\mathcal{R}(A)\cap\mathcal{R}(A^{\mathsf{T}})=\{\boldsymbol{0}\}$, it is called a disjoint range (DR) matrix [3].

3 GMRES and EP matrices

As was already noted, the GMRES iterate $\boldsymbol{x}_{k}=\boldsymbol{x}_{0}+\boldsymbol{z}_{k}$ solves the least squares problem eq. 2.1. Therefore, the restriction of $A$ to the Krylov subspace $\mathcal{K}_{k}\subseteq\mathbb{R}^{n}$ denoted by $A\rvert_{\mathcal{K}_{k}}$ plays an important role in the numerical behavior of GMRES. Indeed, the ill-conditioning of $A\rvert_{\mathcal{K}_{k}}$ was studied and its condition number

[TABLE]

was introduced by Brown and Walker in [6]. In practical computations, the iterate $\boldsymbol{x}_{k}$ is computed as $\boldsymbol{x}_{k}=\boldsymbol{x}_{0}+Q_{k}\boldsymbol{y}_{k}$, where the columns of $Q_{k}$ form an orthonormal basis of the Krylov subspace $\mathcal{K}_{k}$ and the vector $\boldsymbol{y}_{k}$ is a solution of the extended Hessenberg least squares problem $\min_{\boldsymbol{y}\in\mathbb{R}^{k}}\|\beta\boldsymbol{e}_{1}-H_{k+1,k}\boldsymbol{y}\|$ (see section 2). The accuracy of $\boldsymbol{x}_{k}$ is thus affected directly by the conditioning of the matrix $H_{k+1,k}$, whereas the identity $\kappa(H_{k+1,k})=\kappa(A\rvert_{\mathcal{K}_{k}})$ follows from the identities

[TABLE]

Next, we give bounds on the extremal singular values of $H_{k+1,k}$. The norm of the matrix $H_{k+1,k}$ can be always bounded above by that of $A$

[TABLE]

This approach cannot be used to bound the $k$th (or smallest) singular value of $H_{k+1,k}$ due to

[TABLE]

as the last equality holds for $A$ singular.

In the consistent case, the condition number is bounded by $\kappa(H_{k+1,k})\leq\kappa(A\rvert_{\mathcal{R}(A)})$ from $\mathcal{K}_{k}\subseteq\mathcal{R}(A)$ and

[TABLE]

where $A\rvert_{\mathcal{R}(A)}$ denotes the restriction of $A$ to the range $\mathcal{R}(A)$. If $A$ is an EP matrix $\mathcal{R}(A)=\mathcal{R}(A^{\mathsf{T}})=\mathcal{N}(A)^{\perp}$, then

[TABLE]

and

[TABLE]

Indeed, the consistent EP case is similar to the nonsingular case, and the condition number of the extended Hessenberg matrix $H_{k+1,k}$ is bounded by $\kappa(H_{k+1,k})\leq\kappa(A)$ (cf. [33, Remark 3.2, Theorem 3.6]). Consequently, the rank deficiency of the least squares problem eq. 2.1 cannot occur and GMRES will terminate if a solution is reached at some step with a degeneracy of the Krylov subspace at the next step.

In the inconsistent EP case, the equivalence $\mathcal{R}(A^{\mathsf{T}})=\mathcal{R}(A)$ $\Longleftrightarrow$ $\mathcal{N}(A^{\mathsf{T}})=\mathcal{N}(A)$ shows that the nonzero least squares residual $\boldsymbol{r}_{*}\in\mathcal{N}(A^{\mathsf{T}})$ belongs also to $\mathcal{N}(A)$ and

[TABLE]

It follows from eq. 2.1 that the residual $\boldsymbol{r}_{k-1}$ at step $k-1$ belongs to the Krylov subspace $\mathcal{K}_{k}$ and satisfies $\boldsymbol{r}_{k-1}-\boldsymbol{r}_{*}\in\mathcal{R}(A)$. In addition, due to $A\boldsymbol{r}_{*}=\boldsymbol{0}$ we have

[TABLE]

This result was derived in a somewhat different form in [6, Theorem 2.5]. It is clear that in the inconsistent case, the least squares problem eq. 2.1 becomes ill-conditioned as the GMRES iterate converges to a least squares solution. This situation is illustrated in Figure 3.1. Note that we also have

[TABLE]

This bound indicates that if the norm of $A\boldsymbol{r}_{0}$ is too small, then $H_{k+1,k}$ becomes ill-conditioned and the inaccuracy can be expected at all subsequent steps of GMRES. Finally, since a symmetric matrix is an EP matrix, the above discussion also covers the MINRES method [19] applied to symmetric singular systems.

The conditioning of the extended Hessenberg matrix $H_{k+1,k}$ for GMRES and its relation to the conditioning of $A$ are illustrated on small examples. First for simplicity, consider applying GMRES with $\boldsymbol{x}_{0}=\boldsymbol{0}$ to $A\boldsymbol{x}=\boldsymbol{b}$, where

[TABLE]

The matrix $A$ is EP, it has the range $\mathcal{R}(A)=\mathcal{R}(A^{\mathsf{T}})=\mathrm{span}\{[1,0]^{\mathsf{T}}\}$ and the nullspace $\mathcal{N}(A)=\mathcal{N}(A^{\mathsf{T}})=\mathrm{span}\{[0,1]^{\mathsf{T}}\}$, and its minimum nonzero singular value is $\sigma_{1}(A)=\break\min_{\boldsymbol{z}\in\mathcal{R}(A)\backslash\{\boldsymbol{0}\}}\|A\boldsymbol{z}\|/\|\boldsymbol{z}\|=1$. The first two steps of the Arnoldi process for $A$ and the initial vector $\boldsymbol{q}_{1}=\boldsymbol{b}/\|\boldsymbol{b}\|$ give the decomposition $AQ_{2}=Q_{2}H_{2,2}$, where

[TABLE]

Hence, $H_{2,2}$ is singular and we have $\sigma_{1}(A)=\sigma_{1}(H_{2,2})=1$. Solving $\min_{\boldsymbol{y}\in\mathbb{R}^{2}}\|\beta\boldsymbol{e}_{1}-H_{2,2}\boldsymbol{y}\|$, where $\beta=\sqrt{1+\varepsilon^{2}}$, we have the minimum-norm solution $\boldsymbol{y}_{2}=\break 1/\sqrt{1+\varepsilon^{2}}[1,\varepsilon]^{\mathsf{T}}$ and $\|\boldsymbol{y}_{2}\|=1$. Therefore, $\|\boldsymbol{x}_{2}\|=\|\boldsymbol{y}_{2}\|$ and thus the norm of the iterate is not large, even if $\varepsilon$ is very small. It is also clear that for $\varepsilon=0$ the system eq. 3.3 becomes consistent and then GMRES will deliver the minimum norm solution $\boldsymbol{x}_{*}=[1,0]^{\mathsf{T}}$ in one iteration.

For comparison, we also consider the $2\times 2$ nonsingular ill-posed linear system $A\boldsymbol{x}=\boldsymbol{b}$, where

[TABLE]

with the condition number $\kappa(A)=1/\delta$, where $\delta$ is a small scalar and the right-hand side $\boldsymbol{b}$ is contaminated by the error $[0,\varepsilon]^{\mathsf{T}}$. The exact solution is $[1,0]^{\mathsf{T}}$ and the error-contaminated solution is $[1,\varepsilon/\delta]^{\mathsf{T}}$. GMRES applied to $A\boldsymbol{x}=\boldsymbol{b}$ with $\boldsymbol{x}_{0}=\boldsymbol{0}$ gives the same $Q_{2}$ as the one for eq. 3.3 but the different Hessenberg matrix

[TABLE]

whose condition number is the same as that of $A$.

Thus, we have the upper bound $\|A^{-1}\boldsymbol{b}\|=\|\boldsymbol{x}_{2}\|=\|\boldsymbol{y}_{2}\|\leq\|A^{-1}\|\|\boldsymbol{b}\|=\|H_{2,2}^{-1}\|\|\boldsymbol{r}_{0}\|=\sqrt{1+\varepsilon^{2}}/\delta$. This means that if $\varepsilon\neq 0$ and $\delta$ is small, then the iterates computed by GMRES can be inaccurate. In the exactly singular case with $\delta=0$, this is not the case.

In the following numerical examples, we examine the accuracy of the GMRES iterate with respect to the degree of consistency of linear systems by using the test matrix and right hand side vectors

[TABLE]

where $D=\operatorname{diag}(10^{\frac{0}{63}},10^{\frac{-4}{63}},10^{\frac{-8}{63}},\dots,10^{-4})\in\mathbb{R}^{64\times 64}$, $\boldsymbol{\gamma}=[\gamma,\gamma,\dots,\gamma]^{\mathsf{T}}\in\mathbb{R}^{64}$ and $\boldsymbol{\delta}=\left[\delta,\delta,\dots,\delta\right]^{\mathsf{T}}\in\mathbb{R}^{64}$. Hence, $A$ has the condition number $10^{4}$, and $\boldsymbol{b}\not\in\mathcal{R}(A)$ $\Longleftrightarrow$ $\delta\neq 0$. The degree of inconsistency of the linear system $A\boldsymbol{x}=\boldsymbol{b}$ can be controlled by the ratio between $\gamma$ and $\delta$, as $\|\boldsymbol{b}\rvert_{\mathcal{R}(A)}\|=8\gamma$ and $\|\boldsymbol{b}\rvert_{\mathcal{N}(A^{\mathsf{T}})}\|=8\delta$ hold. Since $\mathcal{R}(A)=\mathcal{R}(A^{\mathsf{T}})$, GMRES should determine the least squares solution of $\min_{\boldsymbol{x}\in\mathbb{R}^{n}}\|\boldsymbol{b}-A\boldsymbol{x}\|$ for all $\boldsymbol{b}\in\mathbb{R}^{128}$ (Theorem 2.3). Throughout all our numerical experiments, we use GMRES and RR-GMRES with the Householder orthogonalization process [28] to ensure the best possible orthogonality among the Arnoldi basis vectors $\boldsymbol{q}_{1}$, $\boldsymbol{q}_{2}$, …, $\boldsymbol{q}_{k}$ and we compute the $k$th residual $\boldsymbol{r}_{k}=\boldsymbol{b}-A\boldsymbol{x}_{k}$ explicitly from $\boldsymbol{x}_{k}$ by solving the extended Hessenberg least squares problem $\min_{\boldsymbol{y}\in\mathbb{R}^{k}}\|\beta\boldsymbol{e}_{1}-H_{k+1,k}\boldsymbol{y}\|$ with the Matlab backslash solver, which utilizes column pivoting. A mathematically equivalent solution of $\min_{\boldsymbol{y}\in\mathbb{R}^{k}}\|\beta\boldsymbol{e}_{1}-H_{k+1,k}\boldsymbol{y}\|$ was presented in [17].

Figure 3.4 shows the relative residual norm $\|A^{\mathsf{T}}\boldsymbol{r}_{k}\|/\|A^{\mathsf{T}}\boldsymbol{b}\|$ versus the number of iterations of GMRES in the weakly inconsistent cases $(\gamma,\delta)=(1,0)$, $(1,10^{-12})$, $(1,10^{-8})$, and $(1,10^{-4})$ on the left, and in the strongly inconsistent cases $(\gamma,\delta)=(1,1)$, $(10^{-4},1)$, $(10^{-8},1)$, and $(10^{-12},1)$ on the right. Note that the relative residual norm $\|A^{\mathsf{T}}\boldsymbol{r}_{k}\|/\|A^{\mathsf{T}}\boldsymbol{b}\|$ is associated with the normal equations $A^{\mathsf{T}}\!A\boldsymbol{x}=A^{\mathsf{T}}\boldsymbol{b}$, which are mathematically equivalent to the linear least squares problem $\min_{\boldsymbol{x}\in\mathbb{R}^{n}}\|\boldsymbol{b}-A\boldsymbol{x}\|$. Similarly, Figures 3.4 and 3.4 show the relative residual error norm $\|\boldsymbol{r}_{k}-\boldsymbol{r}_{*}\|/\|\boldsymbol{r}_{k}\|$ (cf. equation eq. 3.2) and the extremal singular values of $A$ and $H_{k+1,k}$, respectively. If the inconsistency is small ($\delta\ll\gamma$), then GMRES is sufficiently accurate (Figure 3.2a); otherwise the relative residual norm $\|A^{\mathsf{T}}\boldsymbol{r}_{k}\|/\|A^{\mathsf{T}}\boldsymbol{b}\|$ stagnates before attaining the accuracy on the level of $u\kappa(A)$ (Figure 3.2b), where $u\simeq 1.1\cdot 10^{-16}$ is the unit roundoff. In contrast to the nonsingular case, GMRES deteriorates not only due to the condition number of $A$ but also due to the inconsistency measured here by $\delta>0$. For strongly inconsistent systems with $\delta\gg\gamma$, i.e., for $\boldsymbol{r}_{0}$ close to $\mathcal{N}(A)$, even though $\|A^{\mathsf{T}}\boldsymbol{r}_{k}\|/\|A^{\mathsf{T}}\boldsymbol{b}\|$ is large and stagnates, and $\|\boldsymbol{r}_{*}\|$ and hence $\|\boldsymbol{r}_{k}\|$ are large, the residual $\boldsymbol{r}_{k}$ approaches $\boldsymbol{r}_{*}$. Figure 3.3b shows that for strongly inconsistent problems, $H_{k+1,k}$ has a condition number significantly larger than $A$, tends to become more ill-conditioned in the subsequent steps, and becomes numerically rank-deficient with $u\|H_{k+1,k}\|\|H_{k+1,k}^{\dagger}\|\geq 1$ as the iteration proceeds. In particular, for $\gamma=0$ and $\delta=1$, GMRES breaks down at step $1$ but gives a least squares solution. Comparing Figures 3.2b and 3.3b, we see that in these cases the bound eq. 3.2 gives a reasonably good upper estimate for the smallest singular value of $H_{k+1,k}$. The behavior described by this example is illustrated on a practical problem from a discretization of a partial differential equation with periodic boundary conditions in [6, Experiment 4.2].

A remedy for the ill-conditioning occurring in GMRES due to inconsistency is to form the Krylov subspace $\mathcal{K}_{k}(A,A\boldsymbol{r}_{0})$ by starting with the initial vector $A\boldsymbol{r}_{0}$ in $\mathcal{R}(A)$ instead of $\boldsymbol{r}_{0}$ as is done in RR-GMRES. Note that, on the other hand, the RR-GMRES residual norm is always larger than or equal to the GMRES residual norm (see section 2). Similarly to the above, we show numerical results for RR-GMRES on the same inconsistent linear systems eq. 3.5. Figures 3.7, 3.7 and 3.7 show the same quantities as Figures 3.4, 3.4 and 3.4 for RR-GMRES. For any inconsistency parameter $\delta>0$, the condition number of $H_{k+1,k}^{R}$ is bounded above by the condition number of $A$ and RR-GMRES is sufficiently accurate, as

[TABLE]

for $\mathcal{R}(A)=\mathcal{R}(A^{\mathsf{T}})$. Thus, the accuracy of the RR-GMRES iterate is affected

only by the condition number of $A$, even though the inconsistency increases or $\boldsymbol{r}_{0}$ approaches $\mathcal{N}(A)$. Hence, for inconsistent problems with EP matrices RR-GMRES is a successful alternative to GMRES.

4 GMRES and GP matrices

We have shown in section 3 that the condition number $\kappa(A\rvert_{\mathcal{R}(A)})$ plays an important role in the behavior of GMRES and in the EP case we have $\kappa(A\rvert_{\mathcal{R}(A)})=\kappa(A)$. Thus for consistent problems with EP matrices, the condition number of $A$ represents an upper bound for the condition number of $H_{k+1,k}$ due to $\kappa(H_{k+1,k})\leq\kappa(A\rvert_{\mathcal{R}(A)})=\kappa(A)$ and the accuracy of the GMRES iterates is actually determined by the singular values of $A$. Consider now applying GMRES to $A\boldsymbol{x}=\boldsymbol{b}$, where $A$ is a GP matrix (Theorem 2.2). We will show that in the GP case $\kappa(A\rvert_{\mathcal{R}(A)})$ can be significantly larger than $\kappa(A)$ and thus the condition number $\kappa(H_{k+1,k})$ can become larger than $\kappa(A)$ even in the consistent case. The accuracy of the GMRES iterates can be then affected by the inaccurate solution of the extended Hessenberg least squares problem that can be ill-conditioned even if $A$ is well-conditioned.

According to Theorem 2.2, GMRES in the consistent GP case determines $\boldsymbol{x}_{\#}+(\mathrm{I}-A^{\#}\!A)\boldsymbol{x}_{0}$. The vector $\boldsymbol{x}_{\#}$ that belongs to $\mathcal{R}(A)$ can be related to the vector $\boldsymbol{x}_{*}$ that belongs to $\mathcal{R}(A^{\mathsf{T}})$ as follows:

[TABLE]

which follows from the identity $\boldsymbol{x}_{*}=P_{\mathcal{R}(A^{\mathsf{T}})}\boldsymbol{x}_{\#}=V_{1}V_{1}^{\mathsf{T}}U_{1}U_{1}^{\mathsf{T}}\boldsymbol{x}_{\#}$. Note that $\boldsymbol{x}_{\#}$ may have a large component in $\mathcal{N}(A)$, if the angle between $\mathcal{N}(A)$ and $\mathcal{R}(A)$ is small, which may affect the accuracy of GMRES iterates (see Figure 4.1), due to the ill-conditioning of the extended Hessenberg matrix $H_{k+1,k}$.

In the consistent case, the extremal singular values of $H_{k+1,k}$ can be bounded as

[TABLE]

[TABLE]

Consequently, $\kappa(H_{k+1,k})\leq\kappa(A)\kappa(V_{1}^{\mathsf{T}}U_{1})$ is related to the extremal principal angles between $\mathcal{R}(A)$ and $\mathcal{R}(A^{\mathsf{T}})$ (cf. [30, Theorem 2.1]). The lower bound (4.2) shows that, in the consistent case, the smallest singular value of $H_{k+1,k}$ can be smaller than the smallest nonzero singular value of $A$, depending on the smallest nonzero singular value of $V_{1}^{\mathsf{T}}U_{1}$. In addition, it is easy to see that $\sigma_{k}(H_{k+1,k})$ can be bounded by

[TABLE]

Here, the last inequality is implied by the splitting $A\boldsymbol{r}_{0}=A(\boldsymbol{r}_{0}\rvert_{\mathcal{N}(A)}+\boldsymbol{r}_{0}\rvert_{\mathcal{R}(A^{\mathsf{T}})})=A\boldsymbol{r}_{0}\rvert_{\mathcal{R}(A^{\mathsf{T}})}$. Although Brown and Walker mention in [6, p. 50] that the condition number of $A\rvert_{\mathcal{K}_{k}}$ cannot become arbitrarily large through an unfortunate

choice of $\boldsymbol{b}$ and $\boldsymbol{x}_{0}$, it is clear from eq. 4.3 that if the residual $\boldsymbol{r}_{0}$ has a very small component in $\mathcal{R}(A^{\mathsf{T}})$ then the condition number of $H_{k+1,k}$ can be large for any singular matrix $A$.

In the following, we give illustrative examples that lead to ill-conditioned extended Hessenberg matrix $H_{k+1,k}$ in GMRES. First, we consider GMRES with $\boldsymbol{x}_{0}=\boldsymbol{0}$ applied to $A\boldsymbol{x}=\boldsymbol{b}$, where

[TABLE]

The matrix $A$ has the following ranges and nullspaces

[TABLE]

It is clear that for $\varepsilon=0$ the matrix $A$ is DR (see section 2). In addition, the initial vector $\boldsymbol{r}_{0}$ satisfies $\boldsymbol{r}_{0}\in\mathcal{N}(A)\cap\mathcal{R}(A)$ and thus the GMRES method breaks down at step 1. Now suppose that $0<\varepsilon\ll 1$. Then, the matrix $A$ is GP but not EP, i.e., $\mathcal{R}(A^{\mathsf{T}})\neq\mathcal{R}(A)$. Since $\sigma_{1}(V_{1}^{\mathsf{T}}U_{1})=\varepsilon/\sqrt{1+\varepsilon^{2}}$, we have

[TABLE]

The smallest singular value of $H_{2,1}$ is significantly smaller than the smallest nonzero singular value of $A$, $\sigma_{1}(H_{2,1})=\varepsilon\ll\sqrt{1+\varepsilon^{2}}=\sigma_{1}(A)$. Indeed, the components of $H_{2,1}$ are $H_{2,1}=[\varepsilon,0]^{\mathsf{T}}$. Furthermore, by solving $\min_{\boldsymbol{y}\in\mathbb{R}^{1}}\|\beta\boldsymbol{e}_{1}-H_{2,1}\boldsymbol{y}\|$ with $\beta=1$, we have $\boldsymbol{y}_{1}=1/\varepsilon$, i.e., $\boldsymbol{y}_{1}$ has a large component. We see that $\boldsymbol{x}_{1}=Q_{1}\boldsymbol{y}_{1}$ and $\boldsymbol{x}_{1}=\boldsymbol{x}_{\#}=A^{\#}\boldsymbol{b}=[1/\varepsilon,0]^{\mathsf{T}}$ for $\boldsymbol{b}=[1,0]^{\mathsf{T}}$. Therefore, $\|\boldsymbol{x}_{1}\|$ becomes very large even if the condition number of $A$ and the norm of the right-hand side are small. Thus, the vector $\boldsymbol{x}_{\#}$ contains a large component in $\mathcal{N}(A)$, whereas $\boldsymbol{x}_{*}=1/(1+\varepsilon^{2})[\varepsilon,1]^{\mathsf{T}}$, and the inequalities (4.1) are satisfied.

In the following, we also consider the nonsingular ill-posed linear system $A\boldsymbol{x}=\boldsymbol{b}$, where

[TABLE]

where $\delta$ is a small parameter and the right-hand side $\boldsymbol{b}$ is contaminated by the error $[0,-\varepsilon]^{\mathsf{T}}$. The extremal singular values of $A$ are $\sigma_{1}(A)\simeq 1$ and $\sigma_{2}(A)\simeq\delta\varepsilon$. Indeed, the condition number of $A$ is bounded as

[TABLE]

Then, GMRES applied to $A\boldsymbol{x}=\boldsymbol{b}$ with the error-free right-hand side $\boldsymbol{b}=[1,0]^{\mathsf{T}}$ and $\boldsymbol{x}_{0}=\boldsymbol{0}$ gives the same Arnoldi vector $\boldsymbol{q}_{1}$ and the same extended Hessenberg matrix $H_{2,1}$ as the ones for eq. 4.4. GMRES applied to $A\boldsymbol{x}=\boldsymbol{b}$ with (4.5) give

[TABLE]

Hence, we have the upper bound $\|A^{-1}\boldsymbol{b}\|=\|\boldsymbol{x}_{2}\|=\|\boldsymbol{y}_{2}\|\leq\|A^{-1}\|\|\boldsymbol{b}\|=\|{H_{2,2}}^{-1}\|\|\boldsymbol{r}_{0}\|=\break\sqrt{1+\varepsilon^{2}}/(\delta\epsilon)$. This means that GMRES breaks down in the exact singular case, whereas, in the ill-posed case, GMRES does not break down but the iterates computed by GMRES can be inaccurate as $\delta$, $\epsilon$ go to zero. In both exactly singular case and nonsingular ill-posed case, the iterates computed by GMRES will be inaccurate due to the ill-conditioning of $H_{2,2}$.

The above behavior of GMRES for singular linear systems with respect to the condition number of $V_{1}^{\mathsf{T}}U_{1}$ is illustrated on numerical examples $A\boldsymbol{x}=\boldsymbol{b}$, where

[TABLE]

$D=\mathrm{diag}(d_{1},d_{2},\dots,d_{64})\in\mathbb{R}^{64\times 64}$ is a diagonal matrix whose diagonal entries have the so-called Strakoš distribution [27]

[TABLE]

and $\boldsymbol{f}=(f_{i})\in\mathbb{R}^{64}$ has the entries $f_{i}=10^{-(64-i)\rho/63}$, $j=1,2,\dots,64$. Note that $A$ is GP but not EP. This setting gives well-conditioned $A$ with $\kappa(A)=\sqrt{2/(10^{-2\rho}+1)}\simeq\sqrt{2}$ and ill-conditioned $V_{1}^{\mathsf{T}}U_{1}$ for $\kappa(V_{1}^{\mathsf{T}}U_{1})=10^{\rho}\sqrt{(10^{-2\rho}+1)/2}\simeq 10^{\rho}/\sqrt{2}$ for $\rho\gg 1$. Furthermore, the norms of vectors $\boldsymbol{f}$, $D\boldsymbol{f}$, $D^{2}\boldsymbol{f}$ decrease and this reduction is pronounced as the value of $\rho$ increases.

Figures 4.3 and 4.3 show the relative residual norm and the smallest singular value of $H_{k+1,k}$ and $H_{k+1,k}^{\mathrm{R}}$ versus the number of iterations for GMRES and RR-GMRES, respectively, applied to the above linear systems with $\rho=1$, $4$, $8$, and $12$. As the value of $\rho$ increases, the condition number of the extended Hessenberg matrix increases, and the accuracy of the relative residual for both GMRES and RR-GMRES is significantly lost. It is clear from our experiments that while RR-GMRES does help in the inconsistent EP case by starting with a vector in $\mathcal{R}(A)$ to construct the Krylov subspace, in the GP case both GMRES and RR-GMRES may not give accurate solutions, even in the consistent case, when the condition number of $V_{1}^{\mathsf{T}}U_{1}$ is large.

5 Conclusions

In this paper we have considered the behavior of the GMRES method for solving a linear system $A\boldsymbol{x}=\boldsymbol{b}$, where $A$ is singular. We have discussed two classes of singular matrices (EP and GP) satisfying the conditions under which GMRES converges to a least squares solution and to the group inverse solution, respectively. We have distinguished between the consistent and inconsistent cases and showed that the conditioning of the extended Hessenberg least squares problem can significantly affect the accuracy of approximate solutions computed by GMRES in finite precision arithmetic.

It appears that the consistent EP case is similar to the nonsingular case. The rank deficiency of the extended Hessenberg least squares problem does not occur and GMRES converges to the accurate approximate solution and terminates with a degeneracy of the Krylov space in the next step. If the coefficient matrix is EP, but system $A\boldsymbol{x}=\boldsymbol{b}$ is inconsistent, then despite the theoretical guarantee for convergence to the least squares solution, the extended Hessenberg least squares problem becomes seriously ill-conditioned and this may lead to very inaccurate approximate solutions computed by GMRES. This happens when the distance of the initial residual to the nullspace is too small or when the residual vector converges gradually to the least squares residual. For such cases, a remedy is to use RR-GMRES.

It is known that if the coefficient matrix is GP and the system is consistent, then theoretically GMRES converges to the solution given by the group inverse of $A$. We have shown, however, that the extended Hessenberg least squares problem can be ill-conditioned even in the consistent case. Indeed, the conditioning of the extended Hessenberg matrix $H_{k+1,k}$ in GMRES depends not only on the conditioning of the coefficient matrix $A$ (as it is in the case of consistent EP problems) but also on the smallest principal angle between the spaces $\mathcal{R}(A^{\mathsf{T}})$ and $\mathcal{R}(A)$ that can be quite large. In such cases, both GMRES and RR-GMRES may compute inaccurate approximate solutions.

We believe that under conditions guaranteeing the convergence of GMRES to the generalized least squares solution considered in [31], our results can be extended to singular systems with a general $\mathrm{index}(A)$. Note also that in this paper we assume only exact arithmetic and our considerations form a groundwork for future work on rounding error analysis.

We would like to point out that here we have considered the behavior of GMRES applied to exactly singular problems. The extension of our results to GMRES applied to almost singular (or numerically singular) linear systems is far from straightforward as it was also illustrated by our small examples.

Acknowledgments

We would like to thank the referees for their valuable comments.