---
address:
- ', , , '
- ', '
- ', , '
- ', , '
author:
- 'Zeyu LIAO\*'
- Ken HAYAMI
- Keiichi MORIKUNI
- 'Jun-Feng YIN'
bibliography:
- 'wileyNJD-VANCOUVER.bib'
nocite: '[@*]'
title: A Stabilized GMRES Method for Solving Underdetermined Least Squares Problems
---

Introduction {#sec1}
============

Consider solving the inconsistent underdetermined least squares problem $$\label{eq3}
\min_{x\in \mathbb{R}^n}\|b-Ax\|_2,\qquad A\in \mathbb{R}^{m\times n},\qquad b\in \mathbb{R}^{m},\qquad b\notin \rm{\mathcal{R}}(\it{A}),\qquad m<n,$$ where $A$ is ill-conditioned and may be rank-deficient. Here, $\rm{\mathcal{R}}$$(A)$ denotes the range space of $A$. Such problems may occur in ill-posed problems where $b$ is given by an observation which contains noise. The least squares problem $(\ref{eq3})$ is equivalent to the normal equations $$\label{eq4}
A^{\mathsf{T}}Ax=A^{\mathsf{T}}b.$$

The standard direct method for solving the least squares problem $(\ref{eq3})$ is to use the QR decomposition. However, when $A$ is large and sparse, iterative methods become necessary. The CGLS [@HS] and LSQR [@lsqr] are mathmetically equivalent to applying the conjugate gradient (CG) method to $(\ref{eq4})$. The convergence of these methods deteriorates for ill-conditioned problems and they require reorthogonalization [@Hayami10] to improve the convergence. Here, we say $(\ref{eq3})$ is ill-conditioned if the condition number $\kappa_2(A)=\|A\|_2\|A^\dag\|_2\gg 1$, where $A^\dag$ is the pseudoinverse of $A$. The LSMR [@lsmr] applies MINRES [@paige1975] to $(\ref{eq4})$.

Hayami et al. [@Hayami10] proposed preconditioning the $m\times n$ rectangular matrix $A$ of the least squares problem by an $n\times m$ rectangular matrix $B$ from the right and the left, and using the generalized minimal residual (GMRES) method[@saad1986] for solving the preconditioned least squares problems (AB-GMRES and BA-GMRES methods, respectively). For ill-conditioned problems, the AB-GMRES and BA-GMRES were shown to be more robust compared to the preconditioned CGNE and CGLS, respectively. Note here that the BA-GMRES works with Krylov subspaces in $n$-dimensional space, whereas the AB-GMRES works with Krylov subspaces in $m$-dimensional space. Since $m<n$ in the underdetermined case, the AB-GMRES works in a smaller dimensional space than the BA-GMRES and should be more computationally efficient compared to the BA-GMRES for each iteration. Moreover, the AB-GMRES has the advantage that the weight of the norm in $(\ref{eq3})$ does not change for arbitrary $B$. Thus, we mainly focus on using the AB-GMRES to solve the underdetermined least squares problem $(\ref{eq3})$. Morikuni [@Morikuni13] showed that the AB-GMRES may fail to converge to a least squares solution in finite-precision arithmetic for inconsistent problems. We will review this phenomenon. The GMRES applied to inconsistent problems was also studied in other papers[@bw; @cr; @REICHEL05; @Morikuni15; @MorikuniRozloznik2018SIMAX].

In this paper, we first analyze the deterioration of convergence of the AB-GMRES. To overcome the deterioration, we use the normal equations of the upper triangular matrix arising in the AB-GMRES to change the inconsistent subproblem to a consistent one. In finite precision arithmetic, forming the normal equations for the subproblem will not square its condition number as would be predicted by theory. In the ill-conditioned case, the tiny singular values are shifted upwards due to rounding errors. In finite precision arithmetic, applying the standard Cholesky decomposition to the normal equations will result in a well-conditioned lower triangular matrix, which will ensure that the forward and backward substitutions work stably, and overcome the problem. Numerical experiments on a series of ill-conditioned Maragal matrices[@florida] show that the proposed method converges to a more accurate approximation than the original AB-GMRES. The method can also be used to solve general inconsistent singular systems.

The rest of the paper is organized as follows. In Section 2, we briefly review the AB-GMRES and a related theorem. In Section 3, we demonstrate and analyze the deterioration of the convergence. In Section 4, we propose and present a stabilized GMRES method and explain a regularization effect of the method based on the normal equations for ill-conditioned problems. In Section 5, numerical results for the underdetermined case and the square case are presented. In Section 6, we conclude the paper.

All the experiments in this paper were done using MATLAB R2017b in double precision, unless specified otherwise (where we extended the arithmetic precision by using the Multiprecision Computing Toolbox for MATLAB [@mptfm]), and the computer uesd was Alienware 15 CAAAW15404JP with CPU Inter(R) Core(TM) i7-7820HK (2.90GHz).

Deterioration of convergence of AB-GMRES for inconsistent problems
==================================================================

In this section, we review previous results. First, we introduce the right-preconditioned GMRES (AB-GMRES), which is the basic algorithm in this paper. Then, we show the phenomenon that the convergence of the AB-GMRES deteriorates for inconsistent problems. Finally, we cite a related theorem to analyze the deterioration.

AB-GMRES method
---------------

AB-GMRES for least squares problems applies GMRES to $\min_{u\in \mathbb{R}^{m}}\|b-ABu\|_2$ with $x=Bu$, where $B\in \mathbb{R}^{n\times m}$. Let $x_0$ be the initial solution (in all our numerical experiments, we set $x_0=0$), and $r_0=b-Ax_0.$ Then, AB-GMRES searches for $u$ in the Krylov subspace $\mathcal{K}_i(AB,r_0)=\rm span\it\{r_0,ABr_0,\dots,(AB)^{i-1}r_0\}$. The algorithm is given in Algorithm \[AL1\][@Hayami10]. Here, $H_{i+1,i}=(h_{pq})\in \mathbb{R}^{(i+1)\times i}$ and $e_1=(1,0,\dots,0)^{\mathsf{T}}\in \mathbb{R}^{i+1}.$

Choose $x_0\in \mathbb{R}^{n}$,$r_0=b-Ax_0$,$v_1=r_0/\|r_0\|_2$ $w_i=ABv_i$ $h_{i,j}=w_i^{{\mathsf{T}}}v_j$, $w_i=w_i-h_{j,i}v_j$ $h_{i+1,i}=\|w_i\|_2$, $v_{i+1}=w_i/h_{i+1,i}$ Compute $y_i\in \mathbb{R}^i$ which minimizes $\|r_i\|_2=\|\|r_0\|_2e_1-H_{i+1,i}y_i\|_2$ $x_i=x_0+B[v_1, v_2, \dots, v_i]y_i$, $r_i=b-Ax_i$ stop

\[AL1\]

To find $y_i\in \mathbb{R}^i$ that minimizes $\|r_i\|_2=\|\|r_0\|_2e_1  -H_{i+1,i}y_i\|_2$ in Algorithm \[AB-GMRES method\], the standard approach computes the QR decomposition of $H_{i+1,i}$ $$\label{EQ3}
H_{i+1,i}=Q_{i+1}R_{i+1,i},\qquad Q_{i+1}\in \mathbb{R}^{{(i+1)}\times{(i+1)}},\qquad R_{i+1,i}=\left(                 %左括号
  \begin{array}{ccc}   %该矩阵一共3列，每一列都居中放置
    R_i\\  %第一行元素
    0^{\mathsf{T}}\\  %第二行元素
  \end{array}
\right)\in \mathbb{R}^{{(i+1)}\times{i}}, \qquad R_{i}\in \mathbb{R}^{{i}\times{i}},$$ where $Q_{i+1}$ is an orthogonal matrix and $R_i$ is an upper triangular matrix. Then, backward substitution is used to solve a system with the coefficient matrix $R_i$ as follows

$$\|r_i\|_2=\min_{y_i\in \mathbb{R}^i}\| Q_{i+1}^{\mathsf{T}}\beta e_1-R_{i+1,i}y_i\|_2,$$

where $$\beta=\|r_0\|_2,\qquad
Q_{i+1}^{\mathsf{T}}\beta e_1=\left(                 %左括号
  \begin{array}{ccc}   %该矩阵一共3列，每一列都居中放置
    t_i\\  %第一行元素
    \rho_{i+1}\\  %第二行元素
  \end{array}
\right),\qquad t_i\in \mathbb{R}^i,\qquad \rho_{i+1}\in \mathbb{R},\qquad y_i=R_i^{-1}t_i,$$ $$x_i=V_iy_i=V_i(R_i^{-1}t_i),\qquad V_i=[v_1, v_2, \dots, v_i]\in \mathbb{R}^{n\times i},\qquad V_i^{\mathsf{T}}\it{V_i}=I,$$ where $\rm{I}$ is the identity matrix.

Note the following theorem.

\[th1\] (Corollary 3.8 of Hayami et al.[@Hayami10]) If $\rm{\mathcal{R}}$$(A)=\rm{\mathcal{R}}$$(B^{\mathsf{T}})$ and $\rm{\mathcal{R}}$$(A^{\mathsf{T}})=\rm{\mathcal{R}}$$(B)$, then AB-GMRES determines a least squares solution of $\min_{x\in \mathbb{R}^n}\|b-Ax\|_2$ for all $b\in \mathbb{R}^m$ and for all $x_0\in \mathbb{R}^n$ without breakdown.

Here, breakdown means $h_{i+1, i}=0$ in Algorithm \[AB-GMRES method\]. See Appendix B of [@Morikuni15].

In fact, if $\rm{\mathcal{R}}$$(A^{\mathsf{T}})=\rm{\mathcal{R}}$$(B)$ and $x_0\in \rm{\mathcal{R}}$$(A^{\mathsf{T}})$, the solution is a minimum-norm solution since $x=Bu\in \rm{\mathcal{R}}$$(A^{\mathsf{T}})=\rm{\mathcal{N}}(A)^{\bot}$, where $\rm{\mathcal{N}}$$(A)$ is the null space of $A$.

From now on, we use AB-GMRES to solve $(\ref{eq3})$ with $B=A^{\mathsf{T}}$ and $x=Bu$, which means using the Krylov subspace $\mathcal{K}_i(AA^{\mathsf{T}},r_0)=\langle r_0,AA^{\mathsf{T}}r_0,\dots,(AA^{\mathsf{T}})^{i-1}r_0\rangle$ to approximate $u$. Hence, Theorem \[th1\] guarantees the convergence in exact arithmetic even in the inconsistent case. However, in finite precision arithmetic, AB-GMRES may fail to converge to a least squares solution for inconsistent problems, as shown later.

AB-GMRES for inconsistent problems
----------------------------------

In this section, we perform experiments to show that the convergence of AB-GMRES deteriorates for inconsistent problems. Experiments were done on the transpose of the matrix Maragal$\_3$ [@florida], denoted by Maragal$\_3$T etc. Table \[tb1\] gives the information on the Maragal matrices, including the density of nonzero entries, rank and condition number. Here, the rank and condition number were determined by using the MATLAB functions `spnrank` [@if] and `svd`, respectively.

[r|r|r|r|r|r]{}

& & & & &\

------------------------------------------------------------------------

Maragal$\_$3T & 858& 1682 & 1.27 & 613& 1.10$\times 10^{3}$\
Maragal$\_$4T & 1027& 1964 & 1.32 & 801&9.33$\times 10^{6}$\
Maragal$\_$5T & 3296& 4654 & 0.61 & 2147& 1.19$\times 10^{5}$\
Maragal$\_$6T & 10144& 21251 & 0.25 & 8331 &2.91$\times 10^{6}$\
Maragal$\_$7T & 26525& 46845 & 0.10 & 20843&8.91$\times 10^{6}$\

\[tb1\]

  -----------------------------------------------------------------------------------------------------------------------------------------------------------
   ![$\kappa_2$($R_i$) and relative residual norm versus the number of iterations for Maragal$\_$3T.[]{data-label="lllw1"}](lw1new.eps "fig:"){width="12cm"}
  -----------------------------------------------------------------------------------------------------------------------------------------------------------

Figure \[lllw1\] shows the relative residual norm $\|A^{\mathsf{T}}r_i\|_2/\|A^{\mathsf{T}}b\|_2$ and $\kappa_2$($R_i$) versus the number of iterations for AB-GMRES with $B=A^{\mathsf{T}}$ for Maragal$\_3$T, where $r_i=b-Ax_i$, and the vector $b$ was generated by the MATLAB function `rand` which returns a vector whose entries are uniformly distributed in the interval $(0,1)$. Here $\kappa_2$($R_i$)=$\kappa_2$($H_{i+1,i})$ holds from $(\ref{EQ3})$. The value of $\kappa_2(R_i$) was computed by the MATLAB function `cond`. The relative residual norm $\|A^{\mathsf{T}}r_i\|_2/\|A^{\mathsf{T}}b\|_2$ decreased to $10^{-8}$ until the 525th iteration, and then increased sharply. The value of `cond(R_i)` started to increase rapidly around iterations 450–550. This observation shows that $R_i$ becomes ill-conditioned before convergence. Thus, AB-GMRES failed to converge to a least squares solution. This phenomnenon was observed by Morikuni[@Morikuni13].

The reason why $R_i$ becomes ill-conditioned before convergence in the inconsistent case will be explained by a theorem in the next subsection.

GMRES for inconsistent problems {#sec2333}
-------------------------------

Brown and Walker [@bw] introduced an effective condition number to explain why GMRES fails to converge for inconsistent least squares problems $$\label{eq5}
\it{\min_{x\in \mathbb{R}^m}\|b-\widetilde{A}x\|_2},$$ where $\it{\widetilde{A}\in\mathbb{R}^{m\times m}}$ is singular, in the following Theorem \[th2\].

Let $b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}$ denote the orthogonal projection of $b$ onto $\rm{\mathcal{R}(\widetilde{\it{A}})}$. Assume $\rm{\mathcal{N}}$$(\widetilde{A})=$$\rm{\mathcal{N}}$$(\widetilde{\it{A}}^{\mathsf{T}})$ and grade$(\widetilde{A},b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}})=k$. Here, grade$(\widetilde{A},\widetilde{b})$ for $\widetilde{A}\in \mathbb{R}^{m\times m}$, $\widetilde{b}\in \mathbb{R}^m$ is defined as the minimum $k$ such that $\mathcal{K}_{k+1}(\widetilde{A}, \widetilde{b})=\mathcal{K}_{k}(\widetilde{A}, \widetilde{b})$. Then, dim($\mathcal{K}_k(\widetilde{\it{A}},b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}))=$ dim($\mathcal{K}_{k+1}(\widetilde{\it{A}},b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}))=$dim($\widetilde{\it{A}}\mathcal{K}_{k}(\widetilde{\it{A}},b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}))=$dim($\widetilde{\it{A}}\mathcal{K}_{k+1}(\widetilde{\it{A}},b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}))=k$ (See Appendix \[ap1\]). Since $\rm{\mathcal{N}}$$(\widetilde{\it{A}})=$$\rm{\mathcal{N}}$$(\widetilde{\it{A}}^{\mathsf{T}})$, we obtain $\widetilde{A}b|_{\mathcal{R}(\widetilde{A})}=\widetilde{A}b$ and dim($\widetilde{A}\mathcal{K}_{k+1}(\widetilde{A},b))=$dim$(\widetilde{A}\mathcal{K}_{k+1}(\widetilde{A},b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}))=k$. If $b\notin \mathcal{R}(\widetilde{A})$ and dim($\widetilde{A}\mathcal{K}_{k}(\widetilde{A},b))=k$, dim($\mathcal{K}_{k+1}(\widetilde{A},b))=k+1$ (See Appendix \[ap2\]).

Let $x_0$ be the initial solution and $r_0=b-\widetilde{A}x_0.$ In the inconsistent case, a least squares solution is obtained at iteration $k$, and at iteration $k+1$ breakdown occurs because of dim($\widetilde{A}\mathcal{K}_{k+1}(\widetilde{A},r_0))< $ dim($\mathcal{K}_{k+1}(\widetilde{A},r_0))$, i.e. rank deficiency of $\it{\min_{z\in \mathcal{K}_{k+1}(\widetilde{A},r_0)}\|b-\widetilde{A}(x_0+z)\|_2}=\it{\min_{z\in \mathcal{K}_{k+1}(\widetilde{A},r_0)}\|r_0-\widetilde{A}z\|_2}$[@bw]. This case is also called the hard breakdown[@REICHEL05].

However, even if $\rm{\mathcal{N}}$$(\widetilde{A})=$$\rm{\mathcal{N}}$$(\widetilde{A}^{\mathsf{T}})$, when (\[eq5\]) is inconsistent, the least squares problem $\it{\min_{z\in \mathcal{K}_{i}(\widetilde{A},r_0)}\|r_0-\widetilde{A}z\|_2}$ may become ill-conditioned as shown below.

\[th2\][@bw] Assume $\rm{\mathcal{N}}$$(\widetilde{A})=$$\rm{\mathcal{N}}$$(\widetilde{A}^{\mathsf{T}})$, and denote the least squares residual of *(\[eq5\])* by $r^*$, the residual at the $(i-1)$st iteration by $r_{i-1}$. If $r_{i-1}\neq r^*$, then $$\kappa_2(A_i)\geq \frac{\|A_i\|_2}{\|\bar{A_i}\|_2}\frac{\|r_{i-1}\|_2}{\sqrt{\|r_{i-1}\|_2^2-\|r^*\|_2^2}},$$ where $A_i\equiv \widetilde{A}|_{\mathcal{K}_i(A,r_0)}$and $\bar{A_i}\equiv \widetilde{A}|_{\mathcal{K}_i(A,r_0)+\rm{span}\{\it{r}^*\}}$. Here, $\widetilde{A}|_S$ is the restriction of $\widetilde{A}$ to a subspace $S\subseteq \mathbb{R}^{m}$.

Theorem \[th2\] implies that GMRES suffers ill-conditioning for $b\notin$ $\rm{\mathcal{R}}$$(\widetilde{A})$ as $\|r_i\|$ approaches $\|r^*\|$. We can apply Theorem \[th2\] to AB-GMRES for least-squares problems by setting $\widetilde{A}\equiv AA^{\mathsf{T}}$. Theorem \[th2\] also implies that even if we choose $B$ as $A^{\mathsf{T}}$, which satisfies the conditions in Theorem \[th1\], AB-GMRES still may not converge numerically because of the ill-conditioning of $R_i$, losing accuracy in the solution computed in finite-precision arithmetic when $r_{i-1}$ approaches $r^*$.

Analysis of the deterioration of convergence {#sec3}
============================================

In this section, we illustrate, the deterioration of convergence of GMRES through numerical experiments. There are two points to note in this section. The first point is that the condition number of $R_i$ tends to become very large as the iteration proceeds for inconsistent problems. Due to $H_{i+1,i}=Q_{i+1}R_{i+1,i}$, the condition number of $H_{i+1,i}$ is the same as that of $R_i$, and will also become very large. The second point is as follows. Since $y_i=R_i^{-1}t_i$, $y_i$ is obtained by applying backward substitution to the triangular system $R_iy_i=t_i$. When the triangular system becomes ill-conditioned, backward substitution becomes numerically unstable, and fails to give an accurate solution $y_i$.

Figure \[lllw1\] shows that at step 550 the relative residual norm suddenly increases. To understand this increase, observe the singular values of $R_{550}$.

  --------------------------------------------------------------------------------------------------------------------------------------------------------------------
   ![Singular value distribution of $R_{550}$ for Maragal$\_$3T in double and quadruple precision arithmetic.[]{data-label="nd"}](lw3bignew.eps "fig:"){width="12cm"}
  --------------------------------------------------------------------------------------------------------------------------------------------------------------------

The left of Figure \[nd\] shows the singular values of $R_{550}$ which were computed in double precision arithmetic. The smallest singular value of $R_{550}$ is $3.21\times 10^{-14}$, which means that the triangular matrix $R_{550}$ is very ill-conditioned and nearly singular in double precision arithmetic.

The right of Figure \[nd\] shows the singular values of $R_{550}$ which were computed in quadruple precision arithmetic using the Multiprecision Computing Toolbox for MATLAB [@mptfm]. The smallest singular value of $R_{550}$ is $5.39\times 10^{-15}$. Since quadruple precision is more accurate, from now on, we mainly show singular value distributions computed in quadruple precision.

  ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
   ![$\kappa_2(R_i)$, $\|y_i\|_2$, and $\|t_i-R_iy_i\|_2/\|t_i\|_2$ versus the number of iterations for Maragal$\_$3T.[]{data-label="lllw5"}](lrenew1.eps "fig:"){width="12cm"}
  ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

Figure \[lllw5\] shows $\kappa_2(R_i)$, $\|y_i\|_2$, and the relative residual norm $\|t_i-R_iy_i\|_2/\|t_i\|_2$ versus the number of iterations for AB-GMRES. The relative residual norm increases only gradually when the condition number of $R_i$ is less than $10^{8}$. When the condition number of $R_i$ becomes larger than $10^{10}$, the relative residual norm starts to increase sharply. This observation shows that when the condition number of $R_i$ becomes very large, the backward substitution will fail to give an accurate $y_i$. As a result, we would not get an accurate $x_i$, and the convergence of AB-GMRES would deteriorate.

Stabilized GMRES method
=======================

In this section, we first propose and present a stabilized GMRES method. Then, we explain its regularization effect comparing it with other regularization techniques.

The stabilized GMRES {#sec41}
--------------------

In order to overcome the deterioration of convergence of GMRES for inconsistent systems, we propose solving the normal equations $$\label{eq419}
R_i^{\mathsf{T}}R_iy_i=R_i^{\mathsf{T}}t_i$$instead of $
R_iy_i=t_i,
$ which we will call the stabilized GMRES. This makes the system consistent, and stabilizes the process, as will be shown in the following.

One may also consider using the normal equations of $H_{i+1,i}$. However, before breakdown, we use AB-GMRES, which means we do not have to store $H_{i+1,i}$. We only store $R_i$ and update it in each iteration, which is cheaper.

  -------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
   ![Comparison of the standard AB-GMRES with stabilized and TSVD stabilized AB-GMRES with $\mu=10^{-8}$ for Maragal$\_$3T.[]{data-label="lllw2"}](lw2bignew.eps "fig:"){width="12cm"}
  -------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

Figure \[lllw2\] shows the relative residual norm $\|A^{\mathsf{T}}r\|_2/\|A^{\mathsf{T}}r_0\|_2$ versus the number of iterations for the standard AB-GMRES and stabilized AB-GMRES with $B=A^{\mathsf{T}}$ for Maragal$\_3$T. The stabilized method reaches the relative residual norm level of $10^{-11}$ which improves a lot compared to the standard method. The method which we used for solving the normal equations (\[eq419\]) is the standard Cholesky decomposition. We replace line 8 of Algorithm \[AL1\] by Algorithm \[AL2\].

We first checked that the method works for the standard Cholesky decomposition coded by ourselves. Later we applied the backslash function of Matlab to (\[eq419\]) to speed up. We checked that in the backslash, the Cholesky decomposition method `chol` is used until the GMRES residual norm stagnates at a small level as seen in Figure \[lllw2\]. In order to continue with further GMRES iterations, the `chol` is automatically switched to the `ldl`, which works even for singular systems.

Compute the QR decomposition of $H_{i+1,i}=Q_{k+1}R_{i+1,i}.$ $R_{i+1,i}=\left(                 %左括号
  \begin{array}{ccc}   %该矩阵一共3列，每一列都居中放置
    R_i\\  %第一行元素
    0^{\mathsf{T}}\\  %第二行元素
  \end{array}
\right)$, $Q_{i+1}^{\mathsf{T}}\beta e_1=\left(                 %左括号
  \begin{array}{ccc}   %该矩阵一共3列，每一列都居中放置
    t_i\\  %第一行元素
    \rho_{i+1}\\  %第二行元素
  \end{array}
\right)$,$\widetilde{R_i}=R_i^{\mathsf{T}}R_i,\qquad \widetilde{t_i}=R_i^{\mathsf{T}}t_i$. Compute the Cholesky decomposition of $\widetilde{R_i}=LL^{\mathsf{T}}$. Solve $Lz_i=\widetilde{t_i}$ by forward substitution. Solve $L^{\mathsf{T}}y_i=z_i$ by backward substitution.

In spite of the above mentioned merits of stabilization, solving the normal equations in AB-GMRES is expensive. Actually, we only need the stabilized AB-GMRES when $R_i$ becomes ill-conditioned. Thus, we can speed up the process by switching AB-GMRES to stabilized AB-GMRES only when $R_i$ becomes ill-conditioned. The condition number of an incrementaly enlarging triangular matrix can be estimated by techniques in [@tebbens2014]. In this paper, we adopt the switching strategy by monitoring the relative residual norm $\|A^{\mathsf{T}}r_i\|_2/\|A^{\mathsf{T}}r_0\|_2$. Let ATR($i$)=$\|A^{\mathsf{T}}r_i\|_2/\|A^{\mathsf{T}}r_0\|_2$ for the $i$th iteration. When ATR($v$)/ $\min_{i=1, 2, \dots, v-1}$ATR($i$)&gt;10, we judge that a jump in relative residual norm has occured, and we switch AB-GMRES to stabilized AB-GMRES at the $v$th iteration.

Motivated by the stabilized AB-GMRES, we also applied the truncated singular value decomposition (TSVD) stabilization method and compared it with the stabilized AB-GMRES. The method modifies $R_i$ by truncating singular values smaller than $\mu$. More specifically, let $R_{i}=U\Sigma V^{\mathsf{T}}$ be the SVD of $R_i$, where the columns of $U=[u_1, u_2, \dots, u_i]$ and $V=[v_1, v_2, \dots, v_i]$ are the left and right singular vectors, respectively, and the diagonal entries of $\Sigma=$ diag$(\sigma_1, \sigma_2, \dots, \sigma_i)$ are the singular values of $R_i$ in discending order $\sigma_1\geq \sigma_2\geq \cdots\geq \sigma_i$. Then, the TSVD approximates $R_i\simeq \sum_{j=1}^k \sigma_j u_j v_j^{\mathsf{T}}$ with $k$ such that $\sigma_{k+1}\leq \mu\sigma_1\leq \sigma_k$ and $y_i= R_i^{-1}t_i\simeq \sum_{j=1}^k\frac{1}{\sigma_j}v_ju_j^{\mathsf{T}}t_i$.

When $\mu= 10^{-13}, 10^{-12}, \dots, 10^{-4}$, the method converges but when $\mu$ is smaller than $10^{-13}$ or larger than $10^{-4}$, it diverges and is similar to the original AB-GMRES. Numerical experiments showed that $\mu=\sqrt{\epsilon}\simeq 10^{-8}$, where $\epsilon$ is the machine epsilon (about $10^{-16}$ in double presion arithmetic), gave the best result among $\mu = 10^{-1}, 10^{-2},\dots, 10^{-16}$ in terms of the relavtive residual as shown in Figure \[lllw2\] for the problem Maragal$\_$3T. The convergence behaviour of the TSVD stabilization method is similar to the stabilized AB-GMRES method, which suggests that eliminating tiny singular values which are less than $10^{-8}$ is effctive for sovling problem (\[eq3\]). However, the TSVD method requires computing the truncated singular value decomposition of $R_i$, and requires choosing the value of the threshold parameter $\mu$, whereas the stabilized AB-GMRES does not require either of them.

[c|rrrrr]{}

matrix & Maragal$\_$3T & Maragal$\_$4T & Maragal$\_$5T & Maragal$\_$6T & Maragal$\_$7T\

------------------------------------------------------------------------

iter. &531 & 465 & 1110 & 2440 & 1864\
standard AB-GMRES &1.05$\times 10^{-8}$ & 2.09$\times 10^{-7}$ & 5.35$\times 10^{-6}$ & 8.26$\times 10^{-6}$ & 4.53$\times 10^{-6}$\

------------------------------------------------------------------------

iter. &552 & 598 & 1226 & 3002 & 2459\
stabilized AB-GMRES &5.99$\times 10^{-12}$ & 5.59$\times 10^{-8}$ & 4.22$\times 10^{-6}$ & 3.88$\times 10^{-6}$& 2.80$\times 10^{-7}$\

\[tb10\]

Table \[tb10\] gives more results for the Maragal matrices. The table shows that the stabilized AB-GMRES is more accurate than the standard AB-GMRES. This seems paradoxical, since forming the normal equations whose coefficient matrix $R_i^{\mathsf{T}}R_i$ would square the condition number compared to $R_i$, which would make the ill-conditioned problem even worse. Why can the stabilized AB-GMRES give a more accurate solution? We will explain why the stabilized AB-GMRES works in the next subsection.

Why the stabilized GMRES method works
-------------------------------------

Consider solving $R_iy_i=t_i, R_i\in \mathbb{R}^{i\times i}, t_i\in \mathbb{R}^{i}$ by solving the normal equations (\[eq419\]), which, in theory, squares the condition number and makes the problem become harder to solve numerically. However, in finite precision arithmetic, the condition number of the normal equations is not neccessarily squared. We will continue to illustrate the phenomenon by using the example in Section \[sec3\].

We used the MATLAB function `svd` in quadruple precision arithmetic [@mptfm] to calculate the singular values. The smallest singular value of $R_{550}$ is $ 5.39\times 10^{-15}$, so its square is $2.91\times 10^{-29}$.

  ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
   ![Singular values $\sigma_i($fl$_d(R_{550}^{\mathsf{T}}R_{550}))$, $i=1, 2, \dots, 550$ in quadruple precision arithmetic.[]{data-label="lll7"}](lll8big.eps "fig:"){width="12cm"}
  ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

Let fl$(\cdot)$ denote the evaluation of an expression in floating point arithmetic and fl$_d(\cdot)$ and fl$_q(\cdot)$ denote the result in double precision arithmetic and quadruple precision arithimetic, respectively. Figure \[lll7\] shows that, numerically, the smallest singular value of fl$_d(R_{550}^{\mathsf{T}}R_{550})$ is $ 7.21\times 10^{-14}$, which is much larger than $2.91\times 10^{-29}$. Further, the Cholesky factor $L$ of fl$_d(R_{550}^{\mathsf{T}}R_{550})~=~LL^{\mathsf{T}}$ computed in double precision precision arithmetic has the smallest singular value $3.50\times 10^{-7}$, which is also larger than $\sqrt{2.91\times 10^{-29}}=5.39\times 10^{-15}$. Thus, the triangular systems $Lz_i=\widetilde{t_i}$ and $L^{\mathsf{T}}y_i=z_i$ are better-conditioned than $R_iy_i=t_i$, which will ensure the stability of the forward and backward substitutions and succeeds in obtaining a much more accurate solution than the standard approach.

  ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
   ![Singular values $\sigma_i($fl$_d(R_{550}^{\mathsf{T}}R_{550}))$, $\sigma_i(R_{550})^2$, $\sigma_i$(fl$_d(R_{610}^{\mathsf{T}}R_{610}))$, and $\sigma_i(R_{610})^2$ in quadruple precision arithmetic.[]{data-label="lll8"}](lm2new.eps "fig:"){width="12cm"}
  ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

The left of Figure \[lll8\] compares the singular values $\sigma_i($fl$_d(R_{550}^{\mathsf{T}}R_{550}))$ and $\sigma_i(R_{550})^2, i=1, 2, \dots, 550$. The first to the 549th singular values of fl$_d(R_{550}^{\mathsf{T}}R_{550})$ and the corresponding $\sigma(R_{550})^2$ are almost the same, while the last one is different. What will happen when $R_i$ contains a cluster of small singular values?

The upper triangular matrix $R_{610}$ contains a cluster of small singular values. The right of Figure \[lll8\] compares the singular values $\sigma_i($fl $_q(R_{610}^{\mathsf{T}}R_{610}))$ and $\sigma_i(R_{610})^2$. The larger singular values are the same as the ‘exact’ values, while the smaller singular values become larger than the ‘exact’ ones.

Experiment results show that finite precision arithmetic has the effect of shifting the tiny singular value upwards. That is the reason why the normal equations (\[eq419\]) help to reduce the condition number and makes the problem become better-conditioned.

Next, we computed the multiplication $R_{550}^{\mathsf{T}}R_{550}$ in quadruple precision arithmetic and observed that the smallest singular values of $R_{550}^{\mathsf{T}}R_{550}$ conincided with the squared singular values $\sigma_i(R_{550})^2$ (blue circle symbol) in the left of Figure \[lll8\], unlike in double precision computation. Since the maximum of the elements of | fl$_q(R_{550}^{\mathsf{T}}R_{550})$ $-$ fl$_d(R_{550}^{\mathsf{T}}R_{550})$ | is approximately $8.16\times 10^{-12}$ , double precision arithmetic contains error of the order of $10^{-12}$. Thus, double precision arithmetic has an effect of regularizing the matrix $R_{550}^{\mathsf{T}}R_{550}$, since double precision matrix multiplication is not accurate enough to keep all the information.

Quadruple precision
-------------------

  -------------------------------------------------------------------------------------------------------------------------------------------
   ![Effect of the stabilized method in quadruple precision arithmetic for Maragal$\_$3T.[]{data-label="un1"}](un1.eps "fig:"){width="12cm"}
  -------------------------------------------------------------------------------------------------------------------------------------------

In order to see the effect of the machine precision $\epsilon$ on the convergence of the AB-GMRES, we compared the stabilized AB-GMRES with the AB-GMRES in quadruple precision arithmetic for the problem Maragal$\_$3T in Figure \[un1\]. For both methods, the relative residual norm reached a smaller level of $10^{-16}$ compared to $10^{-12}$ and $10^{-8}$, respectively, for double precision arithmetic in Figure \[lllw2\]. The curve of the relative residual norm became smoother compared to double precision. As seen in Figure \[un1\], the relative residual norm of the AB-GMRES method jumped to $10^{-1}$ after reaching $10^{-16}$, whereas the relative residual norm of the stabilized GMRES stayed around $10^{-16}$.

When the stabilized GMRES method works
--------------------------------------

[\[sec43\]]{} Motivated by the Läuchli matrix [@Higham], we consider solving the following EP (equal projection) problem $A_3x=(1,0,0)^\mathsf{T}$, where $A_3$ is null space symmetric, that is $\rm{\mathcal{N}}$$(A_3)=$$\rm{\mathcal{N}}$$(A_3^{\mathsf{T}})$ with null space $\rm{\mathcal{N}}$$(A_3)=\rm{span}\it\{(1,-1,1)^\mathsf{T}\}.$ $$\label{exep}
  A_3x=\left(
            \begin{array}{ccccc}
               \frac{\sqrt{2}}{2} &\qquad & \frac{\sqrt{2}}{2}- \frac{\sqrt{6\epsilon}}{6} &\qquad&- \frac{\sqrt{6\epsilon}}{6}\\
                \frac{\sqrt{2}}{2}&\qquad & \frac{\sqrt{2}}{2}+ \frac{\sqrt{6\epsilon}}{6}&\qquad& \frac{\sqrt{6\epsilon}}{6} \\
               0 &\qquad& \frac{\sqrt{6\epsilon}}{3} &\qquad&\frac{\sqrt{6\epsilon}}{3}\\
             \end{array}
           \right)x= \left(
             \begin{array}{c}
1\\
0\\
0\\
   \end{array}
           \right),$$ where $\epsilon$ is the machine epsilon.

Apply GMRES with $x_0=0$ to (\[exep\]). Let $R_s\in \mathbb{R}^{s\times s}$ be the upper triangular matrix obtained at the $s$th iteration of GMRES. In the second iteration, after applying the Givens rotation to $H_{3, 2}$, we obtain the following: $$\label{excep1}
R_2= \left(
             \begin{array}{cc}
               1 & 1 \\
               0 & \sqrt{\epsilon} \\
             \end{array}
           \right),\qquad R_2^{\mathsf{T}}R_2= \left(
             \begin{array}{cc}
               1 & 1 \\
               1 & 1+\epsilon \\
             \end{array}
           \right)\simeq \left(
             \begin{array}{cc}
               1 & 1 \\
               1 & 1\\ \end{array}
           \right).$$ Thus, there is a risk that the stabilized GMRES will give a numerically singular matrix $R_2^{\mathsf{T}}R_2$ in finite precision arithmetic for nonsingular $R_2$. We will analyze this phenominon.

We define the following.

$O(\epsilon)$ denotes that there exists a constant $c$ independent of $\epsilon$, such that $-c\epsilon<O(\epsilon)<c\epsilon$. Also, let $$\mathbb{O}(\epsilon)=\left(  \begin{array}{c}
                O(\epsilon)\\
                O(\epsilon)\\
                \vdots\\
              O(\epsilon) \\
             \end{array}\right)\in \mathbb{R}^n ,\quad \mathcal{O}(\epsilon)= [\mathbb{O}(\epsilon), \mathbb{O}(\epsilon), \cdots, \mathbb{O}(\epsilon)]\in \mathbb{R}^{n\times n}.$$ We assume that the basic arithmetic operations op $=$ $+, -, *, /$ satisfy fl$(x$ op $y)= (x $ op $y)(1+O(\epsilon))$ as in [@nh].

Note also that the following hold from [@nh]. Let $x, y\in \mathbb{R}^n, A,B\in \mathbb{R}^{n\times n},$ and $$|x|=\left(  \begin{array}{c}
                |x_1|\\
                |x_2|\\
                \vdots\\
              |x_n| \\
             \end{array}\right)\quad \rm for\quad \it x=\left(\begin{array}{c}
                x_1\\
                x_2\\
                \vdots\\
              x_n \\
             \end{array}\right),$$ $$|A|=\left(  \begin{array}{cccc}
                |a_{11}|&|a_{12}|&\cdots&|a_{1n}|\\
               |a_{21}|&    |a_{22}|        &   \cdots      &|a_{2n}|\\
                  \vdots        &    \vdots        & \ddots&     \vdots      \\
              |a_{n1}|& |a_{n2}| &  \cdots &|a_{nn}| \\
             \end{array}\right)$$ for $A=(a_{pq})$. Then

$$\begin{aligned}
 \rm f\!l\it(x^\mathsf{T}y)&=x^\mathsf{T}y+O(n\epsilon)|x|^{\mathsf{T}}|y|=x^{\mathsf{T}}y+O(n\epsilon),\\
\rm f\!l\it(Ax)&=Ax+\mathbb{O}(n\epsilon)|A||x|=Ax+\mathbb{O}(n\epsilon),\\
\rm f\!l\it(AB)&=AB+O(n\epsilon)|A||B|=AB+\mathcal{O}(n\epsilon).\end{aligned}$$

Note also that the following theorem holds from Theorem 8.10 of [@nh].

\[th333\] Let $T=(t_{pq})\in \mathbb{R}^{n\times n}$ be a triangular matrix and $b\in \mathbb{R}^{n}.$ Then, the computed solution $\hat{x}$ obtained from substitution applied to $Tx=b$ satisfies $$\hat{x} = x+O(n^2\epsilon)M(T)^{-1}|b|.$$ Here, $M(T)=(m_{ij})$ is the comparison matrix such that

$$m_{ij}=\left\{
             \begin{array}{cc}
               |t_{ij}|, & i=j, \\
               -|t_{ij}|, & i\neq j. \\
             \end{array}
           \right.$$

\
Further, we define the following.

Assume $\|A\|_2=O(1).$ We say $A\in \mathbb{R}^{n\times n}$ is numerically nonsingular if and only if $$\label{EQ14}
\rm f\!l\it (Ax)=\mathbb{O}(\epsilon)\quad\Rightarrow\quad x=\mathbb{O}(\epsilon).$$

Note that this definition of numerical nonsingularity agrees with that of numerical rank due to the following.

Let the SVD of $A=U\Sigma V^{\mathsf{T}}$ where $U, V$ are orthogonal matrices and $\Sigma=\rm diag\it (\sigma_1, \sigma_2, \dots, \sigma_n).$ Here, $\|A\|_2=\sigma_1=O(1).$ If the numerical rank of $A$ is $r<n$, there is a $\sigma_i=O(\epsilon),$ $r+1\leq i\leq n.$ Then, $Ax=U\Sigma V^{\mathsf{T}}x=O(\epsilon)$ admits $x'=V^{\mathsf{T}}x=(x_1', x_2', \dots, x_n')^{\mathsf{T}}$ such that $x_i'=O(1)$, and hence $x=\mathbb{O}(1).$ Thus, $A$ is numerical singular. Then, the following theorem holds.

\[th55\] Let $R_s=(r_{pq})\in \mathbb{R}^{s\times s}$ be an upper-triangular matrix and $$R_{s+1}=\left(
             \begin{array}{cc}
               R_s & d \\
              0^{\mathsf{T}} & r_{s+1,s+1} \\
             \end{array}
           \right)\in \mathbb{R}^{(s+1)\times(s+1)}.$$ Assume $R_s$ is numerically nonsingular, and $R_s=\mathcal{O}(1), R_s^{-1}=\mathcal{O}(1), M(R_s)^{-1}=\mathcal{O}(1), d=\mathbb{O}(1)$ and $O(s)=O(s^2)=O(1).$ Then, the following holds: $$\rm f\!l\it(R_{s+1}^{\mathsf{T}}R_{s+1})\rm~is~numerically~nonsingular \it\quad\Longleftrightarrow\quad \rm f\!l \it (r^2_{s+1,s+1})> \rm f\!l \it(d^{\mathsf{T}}d)O(\epsilon).$$

See Appendix \[ap3\].

Theorem \[th55\] gives the necessary and sufficient condition so that the stabilized GMRES works at the $(s+1)$st iteration, i.e. $ R_{s+1}^{\mathsf{T}}R_{s+1}$ is numerically nonsingular.

  -----------------------------------------------------------------------------------------------------------------------------------------------
   ![$r_{i,i}^2$ $(i=s+1)$ and $d^{\mathsf{T}}d$ in stabilized AB-GMRES for Maragal$\_$3T.[]{data-label="t4dr"}](lt4dr.eps "fig:"){width="12cm"}
  -----------------------------------------------------------------------------------------------------------------------------------------------

The difficulty in solving $R_iy_i=t_i$ by backward substitution is not because the diagonals of $R_i$ are tiny. The reason is that $R_i$ has tiny singular values. However, the exceptional example (\[excep1\]) exists where the stabilized AB-GMRES does not work. The condition fl$(r^2_{s+1, s+1})$ $>$ fl$(d^{\mathsf{T}}d)O(\epsilon)$ in Theorem \[th55\] excludes such exceptions.

Figure \[t4dr\] shows $r^2_{s+1, s+1}$ and $d^{\mathsf{T}}d$ together with the convergence of the AB-GMRES and that of the stabilized AB-GMRE for Maragal$\_$3T. The figure shows that upto 613 iterations, the conditions in Theorem \[th55\] are satisfied, and $ R_{s+1}^{\mathsf{T}}R_{s+1}$ is numerically nonsingular, so that the stabilized AB-GMRES works.

Comparison with Tikhonov regularization method
----------------------------------------------

Another approach to stabilize the AB-GMRES would be to apply Tikhonov regularization. There are two methods to implement it. The first method is to solve the following square system: $$\label{eq66}
(R_i^{\mathsf{T}}R_i+\lambda I)y_i=R_i^{\mathsf{T}}t_i,\qquad \lambda\geq 0$$ using the Cholesky decomposition.

The second method is to solve the regularized least suqares problem $$\label{eqz16}
  \min_{y_i\in\mathbb{R}^i}\left|\left|  \left( \begin{array}{c}
t_i\\
0\\
   \end{array}\right)-\left(\begin{array}{c}
R_i\\
\sqrt{\lambda} I\\
   \end{array} \right)y_i\right|\right|_2$$ using the QR decomposition.

These two methods are equivalent mathematically. However, they are not equivalent numerically. The behavior of the first method is similar to the stabilized AB-GMRES. Table \[tb33\] shows that AB-GMRES combined with the first method converges better when $\lambda=10^{-16}$ than when $\lambda=10^{-14}$. This method can be used to shift upwards the small singular values, but is less acurrate compared to the stabilized AB-GMRES (cf. Table \[tb10\]).\

  -------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
   ![Relative residual norm for the regularized AB-GMRES using (\[eqz16\]) versus number of iterations for different $\lambda$ for Maragal$\_$3T.[]{data-label="llwre999"}](llwre999.eps "fig:"){width="12cm"}
  -------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

Table \[tb33\] also shows that the second method is even more accurate compared with the stabilized AB-GMRES method. There is no need to form the normal equations, so that less information is lost due to rounding error. However, one needs to choose an appropriate value for the regularization parameter $\lambda$. Figure \[llwre999\] shows the relative residual norm $\| A^\mathsf{T} r_i \|_2 / \| A^\mathsf{T} r_0\|_2$ for the regularized AB-GMRES using (\[eqz16\]) versus the number of iterations for different values of $\lambda$ for Maragal$\_$3T. According to Figure \[llwre999\], $\lambda=10^{-16}$ was optimal among $10^{-12}, 10^{-14}, 10^{-16}, 10^{-18}$, so we recommend this value in practice.

We here note the following.

Let $\sigma_1\geq \sigma_2\geq\cdots\geq\sigma_i$ be the the singular values of $R_i$. Then, the singular values of $$R_i'=\left( \begin{array}{c}
R_i\\
\sqrt{\lambda}I\\
   \end{array}\right)$$ are given by $\sqrt{\sigma_1^2+\lambda}\geq\sqrt{\sigma_2^2+\lambda}\geq\cdots\geq\sqrt{\sigma_i^2+\lambda}.$

See Appendix \[ap4\].

Then, let $$\kappa\equiv\kappa_2(R_i)=\frac{\sigma_1}{\sigma_i},\quad \kappa'^2\equiv\kappa_2(R_i')^2=\frac{\sigma_1^2+\lambda}{\sigma_1^2/\kappa^2+\lambda}=1+\frac{\sigma_1^2(1-1/\kappa^2)}{\sigma_1^2/\kappa^2+\lambda}.$$ Since $\kappa\geq 1,  \rm d\it \kappa'/\rm d\it \lambda \leq 0$ for $\lambda\geq 0$ and $\kappa'(\lambda=0)=\kappa,  \kappa'(\lambda=+\infty)=1.$ Note also that $$\lambda=\frac{\sigma_1^2[1+(\kappa'/\kappa)^2]}{\kappa'^2-1}$$ Therefore, for instance, if $\kappa\gg 1$ and we want $\kappa'=\sqrt{\kappa}$, $$\lambda=\frac{\sigma_1^2(1+1/\kappa)}{\kappa-1}\simeq \frac{\sigma_1^2}{\kappa}.$$ For example, if $\kappa=10^{16}$ and we want $\kappa'=10^8$, we should choose $\lambda\simeq \sigma_1^2\times 10^{-16}.$ For Maragal$\_$3T, the largest singular value $\sigma_1$ is about 12.64, so that we can estimate a reasonable value of $\lambda \simeq 1.60 \times 10^{-14}$. However, this estimation assumes $\kappa'=\sqrt{\kappa}$, and needs an extra cost for computing $\sigma_1$. See [@brezinski2009error] for other estimation techniques for the regularization parameter.

[c|rrrrr]{}

matrix & Maragal$\_$3T & Maragal$\_$4T & Maragal$\_$5T & Maragal$\_$6T & Maragal$\_$7T\

------------------------------------------------------------------------

iter. &552 & 597 & 1304 & 2440 & 1864\
method (\[eq66\]) $\lambda=10^{-14}$ &5.08$\times 10^{-11}$ & 5.57$\times 10^{-8}$ & 1.05$\times 10^{-5}$ & 8.26$\times 10^{-6}$ & 4.53$\times 10^{-6}$\

------------------------------------------------------------------------

iter. &570 & 598 & 1226 & 2440 & 1864\
method (\[eq66\]) $\lambda=10^{-16}$ &5.80$\times 10^{-12}$ & 5.59$\times 10^{-8}$ & 4.22$\times 10^{-6}$ & 8.26$\times 10^{-6}$& 4.53$\times 10^{-6}$\

------------------------------------------------------------------------

iter. &553 & 547 & 1261 & 2937 & 2475\
method (\[eqz16\]) $\lambda=1.6\times 10^{-14}$ &7.54$\times 10^{-11}$ & 5.59$\times 10^{-8}$ & 1.15$\times 10^{-5}$ & 9.12$\times 10^{-6}$& 2.78$\times 10^{-7}$\

------------------------------------------------------------------------

iter. &551 & 547 & 1262 & 3037 & 2475\
method (\[eqz16\]) $\lambda=10^{-16}$ &3.37$\times 10^{-12}$ & 5.59$\times 10^{-8}$ & 5.64$\times 10^{-7}$ & 1.91$\times 10^{-6}$& 2.78$\times 10^{-7}$\

\[tb33\]

Comparisons with other methods
==============================

Underdetermined inconsistent least squares problems
---------------------------------------------------

First, we compared the stabilized AB-GMRES with the range restricted AB-GMRES (RR-AB-GMRES) [@NEUMAN2], where the Krylov subspace for the RR-AB-GMRES with $B=A^{\mathsf{T}}$ is $\textit{K}_i(AA^{\mathsf{T}},AA^{\mathsf{T}}r_0)$, AB-GMRES with $B=A^{\mathsf{T}}$, BA-GMRES with $B=A^{\mathsf{T}}$, LSQR[@lsqr] and LSMR[@lsmr]. All programs for iterative methods were coded according to the algorithms in [@NEUMAN2; @Hayami10; @lsqr; @lsmr]. Each method was terminated at the iteration step which gives the minimum relative residual norm within $m$ iterations, where $m$ is the number of the rows of the matrix. No restarts were used for GMRES. Experiments were done for rank-deficient matrices whose information is given in Table 1. Here, we have deleted the zero rows and columns of the test matrices beforehand. The elements of $b$ were randomly generated using the MATLAB function `rand`. Each experiment was done 10 times for the same right hand side $b$ and the average of the CPU times are shown. Symbol - denotes that $\|A^{\mathsf{T}}r_i\|_2/\|A^{\mathsf{T}}r_0\|_2$ did not reach $10^{-8}$ within $20n$ iterations.

[c|rrrrr]{}

matrix & Maragal$\_$3T & Maragal$\_$4T & Maragal$\_$5T & Maragal$\_$6T & Maragal$\_$7T\

------------------------------------------------------------------------

iter. &531 & 465 & 1110 & 2440 & 1864\
standard AB-GMRES &1.05$\times 10^{-8}$ & 2.09$\times 10^{-7}$ & 5.35$\times 10^{-6}$ & 8.26$\times 10^{-6}$ & 4.53$\times 10^{-6}$\

------------------------------------------------------------------------

iter. &552 & 598 & 1226 & 3002 & 2459\
stabilized AB-GMRES &5.99$\times 10^{-12}$ & 5.59$\times 10^{-8}$ & 4.22$\times 10^{-6}$ & 3.88$\times 10^{-6}$& 2.80$\times 10^{-7}$\

------------------------------------------------------------------------

iter. &553 & 565 & 1223 & 2374 & 2474\
RR-AB-GMRES &2.57$\times 10^{-11}$ & 5.59$\times 10^{-8}$ & 3.62$\times 10^{-6}$ & 1.63$\times 10^{-5}$& 2.78$\times 10^{-7}$\

------------------------------------------------------------------------

iter. &562 & 626 & 1263 & 4373 & 5658\
BA-GMRES &2.88$\times 10^{-14}$ & 7.92$\times 10^{-11}$ & 2.29$\times 10^{-12}$ & 5.12$\times 10^{-11}$& 2.03$\times 10^{-10}$\

------------------------------------------------------------------------

iter. &1682 & 2375 & 4576 & 151013 & 97348\
LSQR &5.64$\times 10^{-14}$ & 2.77$\times 10^{-10}$ & 1.11$\times 10^{-11}$ & 5.87$\times 10^{-10}$& 1.33$\times 10^{-9}$\

------------------------------------------------------------------------

iter. &1654 &2308 & 4273 & 127450 & 70242\
LSMR &5.51$\times 10^{-14}$ & 3.00$\times 10^{-10}$ & 3.25$\times 10^{-11}$ & 4.16$\times 10^{-10}$& 9.95$\times 10^{-10}$\

\[tb22\]

[c|rrrrr]{}

matrix & Maragal$\_$3T & Maragal$\_$4T & Maragal$\_$5T & Maragal$\_$6T & Maragal$\_$7T\

------------------------------------------------------------------------

iter. &- & - & - & - & -\
standard AB-GMRES &- & - & - & - & -\

------------------------------------------------------------------------

iter. &546 (526) & - & - & - & -\
stabilized AB-GMRES &2.01 & - & - & -& -\

------------------------------------------------------------------------

iter. &545 & - & - & - & -\
RR-AB-GMRES &1.84 & - & - & -& -\

------------------------------------------------------------------------

iter. &530 & 608 & 1232 & 3623 & 5001\
BA-GMRES &2.10 & 3.19 & 4.25$\times 10^{1}$ &1.81$\times 10^{3}$& 9.20$\times 10^{3}$\

------------------------------------------------------------------------

iter. &1465 & 2120 & 4032 & 101893 & 54444\
LSQR &1.27$\times 10^{-1}$ & 2.56$\times 10^{-1}$ & 1.49 & 2.93$\times 10^{2}$& 4.33$\times 10^{2}$\

------------------------------------------------------------------------

iter. &1456 &1989 & 4013 & 54017 & 31206\
LSMR &1.25$\times 10^{-1}$ & 2.37$\times 10^{-1}$ & 1.49 & 1.50$\times 10^{2}$& 2.23$\times 10^{2}$\

\[ntb1196\]

Table \[tb22\] shows that the stabilized AB-GMRES is generally more accurate than the RR-AB-GMRES. The stabilized AB-GMRES took more iterations to attain the same order of the smallest residual norm than the RR-AB-GMRES. Table \[tb22\] also shows that for the same underdetermined least squares problems, the BA-GMRES was the best in terms of the attainable smallest relative residual norm and that the LSQR and LSMR are comparable to the BA-GMRES, but require less CPU time according to Tabel \[ntb1196\].

Inconsistent systems with highly ill-conditioned square coefficient matrices
----------------------------------------------------------------------------

The stabilized AB-GMRES is not restricted to solving underdetermined problems but can also be applied to solving the least squares problem $\min_{x\in \mathbb{R}^n} \|b-Ax\|_2$, where $A\in \mathbb{R}^{n\times n}$ is a highly ill-conditioned square matrix. Thus, we also test on square matrices of different kinds. Table \[tb94\] gives the information of the matrices.

These matrices are all numerically singular. We generated the right-hand side $b$ by the MATLAB function `rand`, so that the systems are generically inconsistent. We compared the stabilized AB-GMRES with the standard AB-GMRES, RR-AB-GMRES, BA-GMRES with $B=A^{\mathsf{T}}$, LSMR [@lsmr], and LSQR [@lsqr]. Table \[tb95\] gives the smallest relative residual norm and the number of iterations. Table \[tb695\] gives the CPU times in seconds required to obtain relative residual norm $\|A^{\mathsf{T}}r_i\|_2/\|A^{\mathsf{T}}r_0\|_2<10^{-8}$. The switching strategy which was introduced in Section \[sec41\] was used for the stabilized AB-GMRES when measuing CPU times. The number of iterations when switching occurred is in brackets.

[r|r|r|r|r|r]{}

& & & & &\

------------------------------------------------------------------------

Harvard500 & 500 & 1.05 &170 &1.30$\times 10^2$ & web connectivity\
netz4504 & 1961 & 0.13 & 1342&3.41$\times 10^1$ & 2D/3D finite element problem\
TS & 2142 & 0.99 & 2140 &3.52$\times 10^3$ &counter example problem\
grid2$\_$dual & 3136 & 0.12 & 3134 &8.58$\times 10^3$ &2D/3D finite element problem\
uk & 4828 & 0.06 & 4814 & 6.62$\times 10^3$&undirected graph\
bw42 & 10000 & 0.05 & 9999 &2.03$\times 10^{3}$ & partial differential equation[@bw]\

\[tb94\]

[c|rrrrrr]{}

matrix & Harvard500 & netz4504 &TS & grid2$\_$dual & uk & bw42\

------------------------------------------------------------------------

iter. &104 & 144 & 1487 & 3134 & 4620 &715\
standard AB-GMRES &9.38$\times 10^{-9}$ & 4.51$\times 10^{-10}$ &1.56$\times 10^{-9}$ &5.98$\times 10^{-10}$& 1.35$\times 10^{-9}$ & 8.06$\times 10^{-8}$\

------------------------------------------------------------------------

iter. &175 & 201 & 1617 & 3135 & 4779& 788\
stabilized AB-GMRES &4.53$\times 10^{-14}$ & 1.51$\times 10^{-14}$ & 1.54$\times 10^{-9}$ & 1.14$\times 10^{-9}$&6.81$\times 10^{-10}$& 1.66$\times 10^{-7}$\

------------------------------------------------------------------------

iter. &135 & 200 & 1652 & 3134 & 4706&1163\
RR-AB-GMRES &7.78$\times 10^{-14}$ & 3.36$\times 10^{-14}$ & 4.56$\times 10^{-9}$ & 6.52$\times 10^{-8}$& 8.33$\times 10^{-8}$&1.56$\times 10^{-5}$\

------------------------------------------------------------------------

iter. &139 & 194 & 1628 & 3134 & 4724&1520\
BA-GMRES &1.91$\times 10^{-15}$ & 7.27$\times 10^{-16}$ & 8.43$\times 10^{-13}$ & 1.23$\times 10^{-13}$& 6.94$\times 10^{-14}$& 1.97$\times 10^{-11}$\

------------------------------------------------------------------------

iter. &391 & 198 & 6047 & 12549 & 6249&1256\
LSQR &3.59$\times 10^{-15}$ & 5.86$\times 10^{-16}$ & 1.96$\times 10^{-12}$ & 2.51$\times 10^{-13}$& 6.56$\times 10^{-14}$&1.59$\times 10^{-11}$\

------------------------------------------------------------------------

iter. &338 &195 & 6219 & 12497 & 6199& 1212\
LSMR &2.01$\times 10^{-15}$ & 5.97$\times 10^{-16}$ & 1.25$\times 10^{-12}$ & 2.34$\times 10^{-13}$& 7.35$\times 10^{-14}$&1.60$\times 10^{-11}$\

\[tb95\]

[r|rr]{}

& &\

------------------------------------------------------------------------

standard GMRES & 147& 8.08$\times 10^{-9}$\
stabilized GMRES & 219& 1.94$\times 10^{-11}$\
RR-GMRES & 220& 3.13$\times 10^{-11}$\

\[tbbw42\]

[c|rrrrrr]{}

matrix & Harvard500 & netz4504 &TS & grid2$\_$dual & uk & bw42\

------------------------------------------------------------------------

iter. &104 & 134 & 1411 & 3134 & 4583 &-\
standard AB-GMRES &4.72$\times 10^{-2}$ & 1.87$\times 10^{-1}$ &2.14$\times 10$ &2.16$\times 10^{2}$& 6.93$\times 10^{2}$ & -\

------------------------------------------------------------------------

iter. &104 & 134 & 1531 (182) & 3134 & 4679 (4199)& -\
stabilized AB-GMRES &4.78$\times 10^{-2}$ & 1.89$\times 10^{-1}$ & 8.19$\times 10$ & 2.21$\times 10^{2}$&1.93$\times 10^{3}$& -\

------------------------------------------------------------------------

iter. &114 & 153 & 1530 & - & -&-\
RR-AB-GMRES &6.42$\times 10^{-2}$ & 2.62$\times 10^{-1}$ & 2.68$\times 10$ & -& -&-\

------------------------------------------------------------------------

iter. &103 & 131 & 1379 & 3134 & 4562&738\
BA-GMRES &5.48$\times 10^{-2}$ & 1.72$\times 10^{-1}$ & 2.06$\times 10$ & 2.44$\times 10^{2}$& 7.55$\times 10^{2}$& 2.33$\times 10$\

------------------------------------------------------------------------

iter. &222 & 134 & 4239 & 11802 & 5948&913\
LSQR &5.63$\times 10^{-3}$ & 6.61$\times 10^{-3}$ & 7.86$\times 10^{-1}$ & 1.15& 8.65$\times 10^{-1}$&3.12$\times 10^{-1}$\

------------------------------------------------------------------------

iter. &215 &132 & 3913 & 11746 & 5898& 655\
LSMR &5.34$\times 10^{-3}$ & 6.42$\times 10^{-3}$ & 7.04$\times 10^{-1}$ & 1.15& 8.42$\times 10^{-1}$& 2.32$\times 10^{-1}$\

\[tb695\]

Table \[tb95\] shows that for most problems the BA-GMRES was the best in terms of accuracy of relative residual norm. The LSQR and LSMR are similar and are comparable to the BA-GMRES, beacuse they all change the inconsistent problem into a consistent problem. The LSQR and LSMR are more suitable for large and sparse problems compared to the BA-GMRES because they require less CPU time and memory.

For Harvard500 and bw42, the AB-GMRES could only converge to the level of $10^{-9}$ regarding the relative residual norm, while the stabilized AB-GMRES converged to the level of $10^{-14}$. The stabilized AB-GMRES was robust in the sense that it could continue to compute even when the upper triangular matrix $R_i$ became seriously ill-conditioned, and the relative residual norm did not increase sharply towards the end, but just stagnated at a low level, just like for consistent problems. Comparing the CPU time in Tabel \[tb695\], LSMR was the fastest. The stabilized AB-GMRES was usually faster than BA-GMRES. Thus, our stabilization method also makes AB-GMRES stable for highly ill-conditioned inconsistent systems with square coefficient matrices.

The coefficient matrix $A$ of bw42 is singular and satisfies $\rm{\mathcal{N}}$$(A)=$$\rm{\mathcal{N}}$$(A^{\mathsf{T}})$. The problem comes from a finite-difference discritization of a PDE with periodic boundary condition (Experiment 4.2 in Brown and Walker[@bw] with the original $b$). Since the matrix is range symmetric, the GMRES, RR-GMRES, and stabilized GMRES can be directly applied to $Ax=b$ (See [@bw] Theorem 2.4, [@hayami2011g] Theorem 2.7, and [@calvetti2000g] Theorem 3.2.) as shown in Table \[tbbw42\]. The stabilized GMRES gave the relative residual norm 1.94$\times 10^{-11}$ for bw42 at the 219th iteration, similar to the BA-GMRES.

Concluding Remarks
==================

We proposed a stabilized AB-GMRES method for ill-conditioned underdetermined and inconsistent least squares problems. It shifts upwards the tiny singular values of the upper triangular matrix appearing in AB-GMRES, making the process more stable, giving better convergence, and more accurate solutions compared to AB-GMRES. The method is also effective for making AB-GMRES stable for inconsistent least squares problems with highly ill-conditioned square coefficient matrices.

Acknowledgments {#acknowledgments .unnumbered}
===============

Ken Hayami was supported by JSPS KAKENHI Grant Number 15K04768.\
Keiichi Morikuni was supported by JSPS KAKENHI Grant Number 16K17639 and Hattori Hokokai Foundation.\
Jun-Feng Yin was supported by the National Natural Science Foundation of China (No. $11971354$).

Proof of statement in section 2.3 {#ap1}
=================================

\[lemma1\] Assume $\rm{\mathcal{N}}$$(\widetilde{A})$ $\cap$ $\rm{\mathcal{R}}$$(\widetilde{\it{A}})=\{0\}$, and grade$(\widetilde{A},b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}})=k$. Then, $\mathcal{K}_{k+1}(\widetilde{\it{A}},b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}})$ = $\widetilde{A}\mathcal{K}_{k}(\widetilde{\it{A}},b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}})$ holds.

Note that $$\widetilde{\it{A}}\mathcal{K}_{k}(\widetilde{\it{A}},b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}})=\rm{span}\it\{\widetilde{\it{A}}b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}, \widetilde{\it{A}}^2b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}, \cdots,\widetilde{\it{A}}^k b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}\}\subseteqq \rm{span}\it\{b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}},\widetilde{\it{A}}b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}},\cdots,\widetilde{A}^k b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}\}=\mathcal{K}_{k+1}(\widetilde{\it{A}},b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}).$$

grade$(\widetilde{\it{A}},b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}})=k$ implies that

$$\mathcal{K}_{k+1}(\widetilde{\it{A}},b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}})=\mathcal{K}_{k}(\widetilde{\it{A}},b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}})= \rm{span}\it\{b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}},\widetilde{\it{A}}b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}},\cdots,\widetilde{\it{A}}^{k-1} b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}\}.$$

Hence, $$\widetilde{\it{A}}^k b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}=c_0b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}+c_1\widetilde{\it{A}}b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}+\cdots+c_{k-1}\widetilde{\it{A}}^{k-1} b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}},\quad c_i\in \mathbb{R}, i=0, 1, 2, \cdots, k-1.$$

If $c_0=0,$ $$\widetilde{\it{A}}^k b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}=c_1\widetilde{\it{A}}b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}+c_2\widetilde{\it{A}}^2b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}+\cdots+c_{k-1}\widetilde{\it{A}}^{k-1} b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}.$$

Hence, $$c_1\widetilde{\it{A}}b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}+c_2\widetilde{\it{A}}^2b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}+\cdots+c_{k-1}\widetilde{\it{A}}^{k-1} b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}-\widetilde{\it{A}}^k b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}=\widetilde{\it{A}}(c_1b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}+\cdots+c_{k-1}\widetilde{\it{A}}^{k-2} b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}-\widetilde{\it{A}}^{k-1} b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}})=0.$$

Hence, $$c_1b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}+c_2\widetilde{\it{A}}^2b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}+\cdots+c_{k-1}\widetilde{\it{A}}^{k-2} b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}-\widetilde{\it{A}}^{k-1} b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}\in \rm{\mathcal{N}}(\widetilde{\it{A}})\cap \rm{\mathcal{R}}(\widetilde{\it{A}}) = \{0\}.$$

which implies

$$\widetilde{\it{A}}^{k-1} b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}=c_1b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}+c_2\widetilde{\it{A}}b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}+\cdots+c_{k-1}\widetilde{\it{A}}^{k-2} b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}.$$

Thus, $$\mathcal{K}_{k}(\widetilde{A},b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}})=\mathcal{K}_{k-1}(\widetilde{\it{A}},b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}),$$ which contradicts with grade$(\widetilde{A},b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}})=k$. Hence, $c_0\neq 0$, and

$$b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}=d_1\widetilde{\it{A}}b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}+d_2\widetilde{\it{A}}^2b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}+\cdots+d_{k-1}\widetilde{\it{A}}^{k-1} b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}+d_k\widetilde{\it{A}}^k b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}.$$

Hence, $$\mathcal{K}_{k+1}(\widetilde{A},b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}})=\rm{span}\it\{b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}},\widetilde{A}b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}},\cdots,\widetilde{A}^k b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}\}\subseteqq\rm{span}\it\{\widetilde{\it{A}}b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}, \widetilde{\it{A}}^2b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}, \cdots,\widetilde{A}^k b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}\}=\widetilde{A}\mathcal{K}_{k}(\widetilde{A},b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}).$$ Thus, $$\mathcal{K}_{k+1}(\widetilde{A},b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}})=\widetilde{A}\mathcal{K}_{k}(\widetilde{A},b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}).$$

\[corollary1\] Assume $\rm{\mathcal{N}}$$(\widetilde{A})$ $=$ $\rm{\mathcal{N}}$$(\widetilde{\it{A}}^{\mathsf{T}})$, and grade$(\widetilde{A},b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}})=k.$ Then, $\mathcal{K}_{k+1}(\widetilde{A},b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}})=\widetilde{A}\mathcal{K}_{k}(\widetilde{A},b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}})$ holds.

$\rm{\mathcal{N}}$$(\widetilde{A})$ $=$ $\rm{\mathcal{N}}$$(\widetilde{\it{A}}^{\mathsf{T}})$ implies that $$\rm{\mathcal{N}}(\widetilde{\it{A}})\cap \rm{\mathcal{R}}(\widetilde{\it{A}})=\rm{\mathcal{N}}(\widetilde{\it{A}}^{\mathsf{T}})\cap \rm{\mathcal{R}}(\widetilde{\it{A}})=\rm{\mathcal{R}}(\widetilde{\it{A}})^{\bot}\cap \rm{\mathcal{R}}(\widetilde{\it{A}})=\{0\}.$$ Hence, from Lemma \[lemma1\], Corollary \[corollary1\] holds.

Proof of statement in section 2.3 {#ap2}
=================================

\[lemma2\] Assume $\rm{\mathcal{N}}$$(\widetilde{A})$ $\cap$ $\rm{\mathcal{R}}$$(\widetilde{\it{A}})=\{0\}$, grade$(\widetilde{A},b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}})=k$, and $b\notin \mathcal{R}(\widetilde{A})$. Then, dim$(\mathcal{K}_{k+1}(\widetilde{A},b))=k+1$ holds.\

Let $c_0, c_1, \dots, c_k\in \mathbb{R}$ satisfy $$c_0b+c_1\widetilde{\it{A}}b+\cdots+c_{k}\widetilde{\it{A}}^{k} b=0.$$ Since $\rm{\mathcal{N}}$$(\widetilde{A})$ $\cap$ $\rm{\mathcal{R}}$$(\widetilde{\it{A}})=\{0\}$, $$b=b|_{\rm{\mathcal{R}}(\widetilde{\it{A}})}\oplus b|_{\rm{\mathcal{N}}(\widetilde{\it{A}})},$$ where $b|_{\rm{\mathcal{N}}(\widetilde{\it{A}})}$ denotes the orthogonal projection of $b$ onto $\mathcal{N}(\widetilde{\it{A}}).$ Hence, $$c_0b|_{\rm{\mathcal{N}(\widetilde{\it{A}})}}+c_0b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}+c_1\widetilde{\it{A}}b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}+\cdots+c_{k}\widetilde{\it{A}}^{k} b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}=0.$$

If $c_0\neq 0$ $$b|_{\rm{\mathcal{N}(\widetilde{\it{A}})}}=-b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}-\frac{c_1}{c_0}\widetilde{A}b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}-\cdots-\frac{c_{k}}{c_0}\widetilde{A}^{k} b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}\in \rm{\mathcal{R}(\widetilde{\it{A}})}.$$

Hence, $$b|_{\rm{\mathcal{N}(\widetilde{\it{A}})}}\in \rm{\mathcal{N}}(\widetilde{\it{A}})\cap \rm{\mathcal{R}}(\widetilde{\it{A}})=\{0\}.$$ Thus, $b|_{\rm{\mathcal{N}(\widetilde{\it{A}})}}=0,$ which contradicts $b\notin \mathcal{R}(\widetilde{\it{A}})$. Hence, we have $c_0= 0$, and $$c_1\widetilde{A}b+c_2\widetilde{A}^2b+\cdots+c_{k}\widetilde{A}^{k} b=c_1\widetilde{A}b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}+c_2\widetilde{A}^2b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}+\cdots+c_{k}\widetilde{A}^{k} b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}=0.$$

But, since $$\rm dim( \rm{span}\it\{\widetilde{\it{A}}b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}, \widetilde{A}^2b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}, \cdots,\widetilde{A}^k b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}\})=\rm dim( \rm\widetilde{\it{A}} {span}\it\{b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}},\widetilde{A}b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}\cdots,\widetilde{A}^{k-1} b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}\})=\rm dim(\it \widetilde{\it{A}}\mathcal{K}_{k}(\widetilde{A},b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}}))=k$$

holds from Lemma \[lemma1\], we have $c_1=c_2=\cdots=c_k=0,$ which implies $ \rm dim(\it \mathcal{K}_{k+1}(\widetilde{\it{A}},b))=k+1$.

Assume $\rm{\mathcal{N}}$$(\widetilde{\it{A}})=$ $\rm{\mathcal{N}}$$(\widetilde{A}^{\mathsf{T}})$, grade$(\widetilde{A},b|_{\rm{\mathcal{R}(\widetilde{\it{A}})}})=k$, and $b\notin \mathcal{R}(\widetilde{A})$. Then, dim$(\mathcal{K}_{k+1}(\widetilde{A},b))=k+1$ holds.\

$\rm{\mathcal{N}}$$(\widetilde{\it{A}})=$ $\rm{\mathcal{N}}$$(\widetilde{A}^{\mathsf{T}})$ implies $\rm{\mathcal{N}}$$(\widetilde{A})$ $\cap$ $\rm{\mathcal{R}}$$(\widetilde{\it{A}})=\{0\}$. Hence, the corollary follows from Lemma \[lemma2\].

Proof of Theorem 4 {#ap3}
==================

Note that $$R_{s+1}^{\mathsf{T}}R_{s+1}=\left(
             \begin{array}{cc}
               R_s & 0 \\
               d^{\mathsf{T}} & r_{s+1,s+1} \\
             \end{array}
           \right)\left(
             \begin{array}{cc}
               R_s &   d \\
            0^{\mathsf{T}}  & r_{s+1,s+1} \\
             \end{array}
           \right)=\left(
             \begin{array}{cc}
               R_s^{\mathsf{T}}R_s & R_s^{\mathsf{T}}d \\
               d^{\mathsf{T}}R_s & d^{\mathsf{T}}d+r_{s+1,s+1}^2 \\
             \end{array}
           \right).$$

Assume fl$(r^2_{s+1,s+1})$ $\leq$ fl$(d^{\mathsf{T}}d)O(\epsilon).$ Then, since\
$$\begin{aligned}
\rm f\!l\it (d^{\mathsf{T}}d)&=d^{\mathsf{T}}d+O(s\epsilon)d^{\mathsf{T}}d=(1+O(s\epsilon))d^{\mathsf{T}}d,\\
\rm f\!l\it(d^{\mathsf{T}}d+r^2_{s+1,s+1})&=(d^{\mathsf{T}}d+r^2_{s+1,s+1})(1+O(s\epsilon))=d^{\mathsf{T}}d(1+O(s\epsilon)),\end{aligned}$$ we have\
$$\label{eqapc1}
\rm f\!l\it(R_{s+1}^{\mathsf{T}}R_{s+1})=\left(
             \begin{array}{cc}
               R_s^{\mathsf{T}}R_s +O(s\epsilon)|R_s|^{\mathsf{T}}|R_s|& R_s^{\mathsf{T}}d+O(s\epsilon)|R_s|^{\mathsf{T}}|d| \\
               d^{\mathsf{T}}R_s+O(s\epsilon)|d|^{\mathsf{T}}|R_s| & d^{\mathsf{T}}d+O(s\epsilon)d^{\mathsf{T}}d\\
             \end{array}
           \right)=\left(
             \begin{array}{c}R_s^{\mathsf{T}}\\ d^{\mathsf{T}}\end{array}\right)\left(
             \begin{array}{cc}R_s & d\end{array}\right)+\mathcal{O}(s\epsilon),$$ since $R_s=\mathcal{O}(1)$ and $d=\mathbb{O}(1).$ Note $$\left(
             \begin{array}{cc}R_s & d\end{array}\right)\left(
             \begin{array}{c}-R_s^{-1}d\\ 1\end{array}\right)=-R_sR_s^{-1}d+d=0,$$ since $R_s$ is nonsingular.

Hence,\
$$\rm f\!l\it(\left(
             \begin{array}{cc}R_s & d\end{array}\right)\left(
             \begin{array}{c}-R_s^{-1}d\\ 1\end{array}\right))=\rm f\!l\it \{R_s  \rm f\!l\it(-R_s^{-1}d)+d\} = [\rm f\!l \it \{R_s\rm f\!l\it(-R_s^{-1}d)\}+d]\{1+O(\epsilon)\}.$$ Note here that $$\rm f\!l\it \{R_s\rm f\!l\it (-R_s^{-1}d)\}=R_s\rm f\!l\it(-R_s^{-1}d)+O(s\epsilon)|R_s||R_s^{-1}d|,$$ and $$\label{eqapc2}
\rm f\!l\it (-R_s^{-1}d)= -R_s^{-1}d+O(s^2\epsilon)M(R_s)^{-1}|d|$$ from Theorem \[th333\]. Hence, $$\rm f\!l\it (\left(
             \begin{array}{cc}R_s & d\end{array}\right)\left(
             \begin{array}{c}-R_s^{-1}d\\ 1\end{array}\right)) =
 O(s^2\epsilon)R_sM(R_s)^{-1}|d|+O(s\epsilon)|R_s||R_s^{-1}d| =\mathbb{O}(s^2\epsilon),$$ since $R_s^{-1}=\mathcal{O}(1)$ and $M(R_s)^{-1}=\mathcal{O}(1).$

Then, $$\rm f\!l\it(R_{s+1}^{\mathsf{T}}R_{s+1}\left(
             \begin{array}{c}-R_s^{-1}d\\ 1\end{array}\right)) = \rm f\!l\it(\{\left(
             \begin{array}{c}R_s^{\mathsf{T}}\\ d^{\mathsf{T}}\end{array}\right)\left(
             \begin{array}{cc}R_s & d\end{array}\right)+\mathcal{O}(s\epsilon)\}\left(
             \begin{array}{c}-R_s^{-1}d+\mathcal{O}(s^2\epsilon)M(R_s)^{-1}|d|\\ 1\end{array}\right))=\mathbb{O}(s^2\epsilon)=\mathbb{O}(\epsilon),$$ since (\[eqapc1\]), (\[eqapc2\]), and $O(s^2)=O(1).$ Since $\left(\begin{array}{c}-R_s^{-1}d\\ 1\end{array}\right)=\mathbb{O}(1),$ $R_{s+1}^{\mathsf{T}}R_{s+1}$ is numerically singular. By contraposition, ($\Leftarrow$) holds.

Assume $R_{s+1}^{\mathsf{T}}R_{s+1}$ is not numerically singular. Then, there exists a vector $\left(
                                                                                                        \begin{array}{c}
                                                                                                          z \\
                                                                                                          w \\
                                                                                                        \end{array}
                                                                                                      \right)\in \mathbb{R}^{s+1}$ such that $\left|\left(
                                                                                                        \begin{array}{c}
                                                                                                          z \\
                                                                                                          w \\
                                                                                                        \end{array}
                                                                                                      \right)\right|>\mathbb{O}(\epsilon),$ and $$\begin{aligned}
\rm f\!l\it \{R_{s+1}^{\mathsf{T}}R_{s+1}\left(
                                                                                                        \begin{array}{c}
                                                                                                          z \\
                                                                                                          w \\
                                                                                                        \end{array}
                                                                                                      \right)\}= R_{s+1}^{\mathsf{T}}\left(R_{s+1}\left(
                                                                                                        \begin{array}{c}
                                                                                                          z \\
                                                                                                          w \\
                                                                                                        \end{array}
                                                                                                      \right)+|R_{s+1}|\left|\left(
                                                                                                        \begin{array}{c}
                                                                                                          z \\
                                                                                                          w \\
                                                                                                        \end{array}
                                                                                                      \right)\right|O((s+1)\epsilon)\right)+\\ \left|R_{s+1}^{\mathsf{T}}\right|\left|R_{s+1}\left(
                                                                                                        \begin{array}{c}
                                                                                                          z \\
                                                                                                          w \\
                                                                                                        \end{array}
                                                                                                      \right)+|R_{s+1}|\left|\left(
                                                                                                        \begin{array}{c}
                                                                                                          z \\
                                                                                                          w \\
                                                                                                        \end{array}
                                                                                                      \right)\right|O((s+1)\epsilon)\right|O((s+1)\epsilon)=\mathbb{O}(\epsilon)\end{aligned}$$ assuming $O(s+1)=O(1).$\
Hence, $$\rm f\!l\it\{R_{s+1}^{\mathsf{T}}R_{s+1}\left(
                                                                                                        \begin{array}{c}
                                                                                                          z \\
                                                                                                          w \\
                                                                                                        \end{array}
                                                                                                      \right)\}=\left(
             \begin{array}{cc}
               R_s^{\mathsf{T}}R_s & R_s^{\mathsf{T}}d \\
               d^{\mathsf{T}}R_s & d^{\mathsf{T}}d+r_{s+1,s+1}^2 \\
             \end{array}
           \right)\left(
                                                                                                        \begin{array}{c}
                                                                                                          z \\
                                                                                                          w \\
                                                                                                        \end{array}
                                                                                                      \right)+\mathbb{O}(\epsilon)=\mathbb{O}(\epsilon).$$ Thus, $$\label{c1}
 R_s^{\mathsf{T}}R_sz+wR_s^{\mathsf{T}}d=\mathbb{O}(\epsilon),$$ $$\label{c2}
d^{\mathsf{T}}R_sz+(d^{\mathsf{T}}d+r_{s+1,s+1}^2)w=\mathbb{O}(\epsilon).$$ (\[c1\]) can be expressed as $R_s^{\mathsf{T}}(R_sz+wd)=\mathbb{O}(\epsilon).$ From Lemma \[lemma3\], $R_s^\mathsf{T}$ is numerically nonsingular, so that $$\label{c3}
R_sz+wd=\mathbb{O}(\epsilon).$$ Hence, from (\[c2\]), $d^{\mathsf{T}}R_sz+w(d^{\mathsf{T}}d+r_{s+1,s+1}^2)=d^{\mathsf{T}}(R_sz+wd)+wr_{s+1,s+1}^2=O(\epsilon).$ Thus, $wr_{s+1,s+1}^2=O(\epsilon)$. If $w=O(\epsilon),$ $R_sz=\mathbb{O}(\epsilon)$ from (\[c3\]). Since $R_s$ is numerically nonsingular, $z=\mathbb{O}(\epsilon),$ which contradicts with the assumption.

Hence, $|w|>O(\epsilon),$ so that $r_{s+1,s+1}^2=O(\epsilon),$ which gives $$\rm  f\!l\it(r_{s+1,s+1}^2) = O(\epsilon)\leq \rm f\!l\it( d^{\mathsf{T}}d)O(\epsilon).$$

\[lemma3\]

Let $n=O(1)$. If $A\in \mathbb{R}^{n\times n}$ is numerically nonsingular, and $A^{-1}=\mathcal{O}(1)$, then $A^{\mathsf{T}}$ is numerically nonsingular.

If $$\rm f\!l\it (A^{\mathsf{T}}x) = A^{\mathsf{T}}x+\mathbb{O}(n\epsilon)|A^{\mathsf{T}}||x|=\mathbb{O}(n\epsilon),$$ then $$\rm f\!l\it(x^{\mathsf{T}}A) = x^{\mathsf{T}}A+\mathbb{O}^{\mathsf{T}}(n\epsilon)=\mathbb{O}^{\mathsf{T}}(n\epsilon).$$ Thus, $$\rm f\!l(\it x^{\mathsf{T}}Ay) = \rm f\!l\it (x^{\mathsf{T}}A)y+O(n\epsilon)|\rm f\!l\it(x^{\mathsf{T}}A)||y|=O(n\epsilon)$$ holds for all $y=\mathbb{O}(1).$

For arbitrary $z=\mathbb{O}(1)\in \mathbb{R}^{n},$ let $$y=A^{-1}z=\mathbb{O}(1).$$ Then, $$\rm f\!l(\it Ay) = Ay+O(n\epsilon)|A||y|=z+O(n\epsilon)|A||y|.$$ Hence, $$z=\rm f\!l\it(Ay) + O(n\epsilon)|A||y| = \rm f\!l\it(Ay) + \mathbb{O}(n\epsilon).$$ Thus, we have $$\rm f\!l(\it x^{\mathsf{T}}z) = x^{\mathsf{T}}z + O(n\epsilon)|x|^{\mathsf{T}}|z| =\rm f\!l\it(x^{\mathsf{T}}Ay) + O(n\epsilon)= O(n\epsilon)$$ for arbitrary $z=\mathbb{O}(1)\in \mathbb{R}^n.$ Hence, $x=\mathbb{O}(\epsilon),$ so that $A^{\mathsf{T}}$ is numerically nonsingular.

Proof of Theorem 5 in section 4.5 {#ap4}
=================================

Let the singular value decomposition of $R_i$ be given by $R_i=U\Sigma V^{\mathsf{T}}\in \mathbb{R}^{i\times i},$ where $U, V$ are orthogonal matrices and $\Sigma=\rm diag\it(\sigma_1, \sigma_2, \dots, \sigma_i).$ Let $I_i\in \mathbb{R}^{i\times i}$ be the identity matrix. Then, we have $R_i'=\left( \begin{array}{c}
R_i\\
\sqrt{\lambda}I_i\\
   \end{array}\right)=U'\Sigma'V^{\mathsf{T}}$, where $U'=\left( \begin{array}{cc}
U & 0\\
0 & V
   \end{array}\right)$ and $\Sigma'=\left( \begin{array}{c}
\Sigma\\
\sqrt{\lambda}I_i\\
   \end{array}\right).$ Since $\Sigma'^{\mathsf{T}}\Sigma'=\Sigma^2+\lambda I_i=\rm diag\it (\sigma_1^2+\lambda, \sigma_2^2+\lambda, \dots, \sigma_i^2+\lambda),$ the singular values of $\left( \begin{array}{c}
R_i\\
\sqrt{\lambda}I_i\\
   \end{array}\right)$ are $\sqrt{\sigma_1^2+\lambda}\geq\sqrt{\sigma_2^2+\lambda}\geq\cdots\geq\sqrt{\sigma_i^2+\lambda}.$
