An incomplete generalized minimum backward perturbation algorithm for large nonsymmetric linear systems

Abstract

This paper gives the truncated version of the generalized minimum backward error algorithm (GMBACK)—the incomplete generalized minimum backward perturbation algorithm (IGMBACK) for large nonsymmetric linear systems. It is based on an incomplete orthogonalization of the Krylov vectors in question, and gives an approximate or quasi-minimum backward perturbation solution over the Krylov subspace. Theoretical properties of IGMBACK including finite termination, existence and uniqueness are discussed in details, and practical implementation issues associated with the IGMBACK algorithm are considered. Numerical experiments show that, the IGMBACK method is usually more efficient than GMBACK and GMRES, and IMBACK, GMBACK often have better convergence performance than GMRES. Specially, for sensitive matrices and right-hand sides being parallel to the left singular vectors corresponding to the smallest singular values of the coefficient matrices, GMRES does not necessarily converge, and IGMBACK, GMBACK usually converge and outperform GMRES.

Graphical abstract

Keywords

Nonsymmetric linear systems / Krylov subspace methods / minimum backward perturbation / incomplete orthogonalization process / GMBACK / GMRES

1 Introduction

In many applied sciences and engineering calculations, it is often necessary to solve large nonsymmetric sparse linear systems:

(1)\begin{array}{l}Ax=b,\end{array}

where A\in {\mathbb{R}}^{n\times n} is nonsingular, x,b\in {\mathbb{R}}^n. The Krylov subspace method is a very effective method for solving (1). Krylov subspace methods [20] such as the generalized minimum residual method (CMRES algorithm [4, 22]) often use the residual norm as a condition for determining the termination of the algorithm. If the approximate solution is accurate, then the residual norm is small. Conversely, a small residual norm does not mean that the approximate solution is accurate, especially when A is a ill-posed matrix. In order to overcome the shortcomings of residual norm as a termination condition, the idea of using backward perturbation norm [10] as the termination condition is proposed in reference [1]. Kasenally minimized the backward perturbation norm {\Vert {\mathrm{\Delta }}_A\Vert }_F on A by finding an approximate solution satisfying the perturbation equation when solving nonsymmetric linear systems (A-{\mathrm{\Delta }}_A)x_m=b, that is, \underset{x_m\in x_0+{\mathcal{K}}_m(A,r_0)}{min}{\Vert {\mathrm{\Delta }}_A\Vert }_F makes

(2)(A-{\mathrm{\Delta }}_A)x_m=b.

The generalized minimum backward perturbation method (CMBACK algorithm [13]) for solving large nonsymmetric linear systems is presented. In this notation, x_m represents the approximate solution of (1) with the form: x_m=x_0+t_m, where t_m\in {\mathcal{K}}_m(A,r_0), x_0 represents the original estimate of the approximate solution, denoted r_0:=b-A{x}_0. {\mathrm{\Delta }}_A and {\mathrm{\Delta }}_b represent perturbations to matrix A and vector b, respectively, forming a joint backward perturbation matrix [{\mathrm{\Delta }}_A,{\mathrm{\Delta }}_b]. In the late 1990s, Kasenally, Simoncini and Zhihao Cao popularized the GMBACK algorithm by perturbing 6 on top of 4, that is, \underset{x_m\in x_0+{\mathcal{K}}_m(A,r_0)}{min}{\Vert {\mathrm{\Delta }}_A,{\mathrm{\Delta }}_b\Vert }_F makes

\begin{array}{l}(A-{\mathrm{\Delta }}_A)x_m=(b+{\mathrm{\Delta }}_b).\end{array}

The minimum joint backward perturbation method for solving large nonsymmetric linear systems is obtained (Minpert algorithm [8, 14]). The GMRES algorithm can be regarded as minimizing the perturbation norm {\Vert {\mathrm{\Delta }}_b\Vert }_2 on b by finding an approximate solution satisfying the equation A{x}_m=b-{\mathrm{\Delta }}_b, that is, \underset{x_m\in x_0+{\mathcal{K}}_m(A,r_0)}{min}{\Vert {\mathrm{\Delta }}_b\Vert }_2 makes

\begin{array}{l}A{x}_m=b-{\mathrm{\Delta }}_b.\end{array}

Therefore, GMBACK and GMRES algorithms can be regarded as minimizing the perturbation norm on the coefficient matrix A and vector b, respectively. While the Minpert algorithm is minimizing the norm of joint backward perturbation matrix [{\mathrm{\Delta }}_A,{\mathrm{\Delta }}_b] on matrix A and vector b.

The GMBACK, GMRES, and Minpert algorithms are all a set of bases v_1,v_2,\dots ,v_m that utilize the Arnoldi process to generate {\mathcal{K}}_m(A,r_0). This means that these algorithms use long recursive formulas, which leads to a dramatic increase in computation and storage capacity with the number of steps, and the algorithms become unusable when the number of steps increases to a certain point. Therefore, these algorithms usually need to be restarted, but for difficult problems, even when combined with preprocessing techniques, the number of steps taken in a restart is still quite large to ensure convergence of the algorithm. In order to overcome the shortcomings of long recursive formulas, a popular technique is to resort to truncation strategies. The idea is to use only a few rather than all the previously calculated vectors to construct a new short recursive formula to calculate the following vectors, which greatly reduces the computation and storage capacity. The truncation form of CMRES algorithm is given in references [5, 6]: quasi-incomplete orthogonalization method (QGMRES) or incomplete generalized minimum Residue Method (ICMRES). Sun Lei provided the truncated form of the Minpert algorithm in reference [25]: incomplete minimum joint backward perturbation method (IMinpert algorithm). However, the truncated form of the CMIBACK algorithm has not yet been given. In this paper, by means of incomplete orthogonalization process [12, 20, 21], a group of bases v_1,v_2,\dots ,v_m of {\mathcal{K}}_m(A,r_0) will be generated. The truncation form of CMBACK algorithm is given: incomplete generalized minimum backward perturbation method (IGMBACK). As with GMBACK, the approximate solution of IGMBACK is generated by solving the minimum problem (2).

The structure of this paper is as follows: Section 2 gives the theoretical derivation process of the ICMBACK algorithm, followed by the IGMBACK algorithm, and some theoretical studies are made on the algorithm, including the finite termination of the algorithm, existence and uniqueness of the solution. Section 3 gives the execution of IGMBACK algorithm. In Section 4, we will show with numerical examples that IGMBACK is generally more efficient than GMBACK and GMRES, and IGMBACK and GMBACK generally converge better than GMRES. In particular, the restarted GMRES algorithm does not necessarily converge if the coefficient matrix is a sensitive matrix and the vector b on the right-hand side of equations is parallel to the left singular vector corresponding to the minimum singular value of A. However, IGMBACK and GMBACK algorithms generally converge and converge better than GMRES. Section 5 is the summary of the full text and the prospect of the work.

2 IGMBACK algorithm

2.1 Analysis of the backward perturbation matrix

The following lemma gives the general formula for the backward perturbed matrix {\mathrm{\Delta }}_A.

Lemma 2.1　 Suppose that an incomplete orthogonalization process of m steps has been performed and a set of basis vectors v_1,v_2,\dots ,v_m and v_{m+1} of {\mathcal{K}}_m(A,r_0) is obtained, denoted V_m=[v_1,v_2,\dots ,v_m], V_{m+1}=[v_1,v_2,\dots ,v_m,v_{m+1}]. At the same time, an upper Hessenberg matrix H_m\in {\mathbb{R}}^{(m+1)\times m} is obtained satisfying:

\begin{array}{l}A{V}_m=V_{m+1}{H}_m=V_m{\overline{H}}_m+h_{m+1,m}{v}_{m+1}{e}_m^T,\end{array}

where the last line of H_m is removed to get {\overline{H}}_m. The approximate solution of the equation (1) can be written as follows:

(3)\begin{array}{l}{x}_m=x_0+V_m{y}_m,y_m\in {\mathbb{R}}^m.\end{array}

Denoted \beta ={\Vert r_0\Vert }_2. Thus the set S=\left\{{\mathrm{\Delta }}_A\right\} of backward perturbation matrices satisfying \left(A-{\mathrm{\Delta }}_A\right)x_m=b can be expressed as:

(4)\begin{array}{l}S=\left\{V_{m+1}\left(H_m{y}_m-\beta e_1\right){\Vert x_m\Vert }_2^{-2}{x}_m^{\mathrm{T}}+R\left(I-w_m{w}_m^{\mathrm{T}}\right):\mathrm{\forall }R\in {\mathbb{R}}^{n\times n}\right\},\end{array}

where w_m=x_m{\Vert x_m\Vert }_2^{-1}.

Proof　Substituting A{V}_m=V_{m+1}{H}_m and (3) into \left(A-{\mathrm{\Delta }}_A\right)x_m=b, we get

\begin{array}{l}{\mathrm{\Delta }}_A{w}_m=V_{m+1}\left(H_m{y}_m-\beta e_1\right){\Vert x_m\Vert }_2^{-1}.\end{array}

Thus for any R\in {\mathbb{R}}^{n\times n}, we immediately obtain (4) (see [3]).

By Lemma 2.1, we can further obtain F- in S the backward perturbation matrix {\mathrm{\Delta }}_{min} with the minimum norm and its norm {\Vert {\mathrm{\Delta }}_{min}\Vert }_F.

Theorem 2.1　 The backward perturbed matrix {\mathrm{\Delta }}_{min} with the minimum norm of F in S can be expressed as

{\mathrm{\Delta }}_{min}:=V_{m+1}\left(H_m{y}_m-\beta e_1\right){\Vert x_m\Vert }_2^{-2}{x}_m^{\mathrm{T}}=-\frac{r_m{w}_m^{\mathrm{T}}}{{\Vert x_m\Vert }_2},

where r_m=b-A{x}_m. Thus the backward perturbation norm {\Vert {\mathrm{\Delta }}_{\text{min}}\Vert }_F=\frac{{\Vert r_m\Vert }_2}{{\Vert x_m\Vert }_2}.

Proof　Require F- in S the backward perturbation matrix with the minimum norm. According to Lemma 2.1, let

g\left(y_m,R\right)=\Vert V_{m+1}\left(H_m{y}_m-\beta e_1\right)\Vert x_m{\Vert {}_2^{-2}{x}_m^{\mathrm{T}}+R\left(I-w_m{w}_m^{\mathrm{T}}\right)\Vert }_F^2.

We can write g\left(y_m,R\right) as

\begin{array}{rl}g\left(y_m,R\right)=& \Vert \mathrm{v}\mathrm{e}\mathrm{c}\left\{V_{m+1}\left(H_m{y}_m-\beta e_1\right){\Vert x_m\Vert }_2^{-2}{x}_m^{\mathrm{T}}\right\}\\ & +\left\{\left(I-w_m{w}_m^{\mathrm{T}}\right)\otimes I\right\}\mathrm{v}\mathrm{e}\mathrm{c}(R){\Vert }_2^2,\end{array}

where \mathrm{v}\mathrm{e}\mathrm{c}(\cdot ) indicates that the columns of the matrix are straightened, \otimes denotes tensor product. From the above equation, the minimum value of g\left(y_m,R\right) is required. By the knowledge of optimization theory, we have {\mathrm{\nabla }}_{\mathrm{v}\mathrm{e}\mathrm{c}(\mathrm{R})}g\left(y_m,R\right)=0. Thus we have R\left(I-w_m{w}_m^{\mathrm{T}}\right)=0, and the conclusion of the theorem holds. □

2.2 Generation of IGMBACK algorithm

This section will derive the IGMBACK algorithm. Lemma 2.1 gives the general formula (4) for the backward perturbation matrix {\mathrm{\Delta }}_A. In the proof of Theorem 2.1, we establish the function g\left(y_m,R\right) according to the general formula (4). It can be seen that the only free variables in g\left(y_m,R\right) are R and y_m, so that the minimum value problem (2) can be transformed into \underset{y_m\in {\mathbb{R}}^m,R\in {\mathbb{R}}^{n\times n}}{min}g\left(y_m,R\right), i.e.,

(5)\begin{array}{l}\underset{y_m\in {\mathbb{R}}^m,R\in {\mathbb{R}}^{n\times n}}{min}\Vert V_{m+1}\left(H_m{y}_m-\beta e_1\right)\Vert x_m{\Vert {}_2^{-2}{x}_m^{\mathrm{T}}+R\left(I-w_m{w}_m^{\mathrm{T}}\right)\Vert }_F^2.\end{array}

According to Theorem 2.1, problem (5) can be further transformed into

(6)\begin{array}{l}\underset{y_m\in {\mathbb{R}}^m}{min}\Vert V_{m+1}\left(H_m{y}_m-\beta e_1\right)\Vert x_m{\Vert {}_2^{-2}{x}_m^{\mathrm{T}}\Vert }_F^2.\end{array}

Since V_{m+1}=\left[v_1,v_2,\dots ,v_m,v_{m+1}\right] is not a column orthogonal norm matrix, this can be computationally intensive, and since

(7)\begin{array}{rl}\Vert V_{m+1}\left(H_m{y}_m-\beta e_1\right)\Vert x_m{\Vert {}_2^{-2}{x}_m^{\mathrm{T}}\Vert }_F& ⩽{\Vert V_{m+1}\Vert }_F\frac{{\Vert H_m{y}_m-\beta e_1\Vert }_2}{{\Vert x_m\Vert }_2}\\ & =\sqrt{m+1}\frac{{\Vert H_m{y}_m-\beta e_1\Vert }_2}{{\Vert x_m\Vert }_2},\end{array}

we can choose y_m such that \frac{{\Vert H_m{y}_m-\beta e_1\Vert }_2}{{\Vert x_m\Vert }_2} is minimal. That is, problem (6) is approximated as

(8)\underset{y_m\in {\mathbb{R}}^m}{min}\frac{{\Vert H_m{y}_m-\beta e_1\Vert }_2^2}{{\Vert x_m\Vert }_2^2}.

For convenience, we define two matrices L_m\in {\mathbb{R}}^{(m+1)\times (m+1)} and G_m\in {\mathbb{R}}^{n\times (m+1)}, where

(9)\begin{array}{l}{L}_m=\left[\begin{array}{l}{H}_m-\beta e_1\end{array}\right],G_m=\left[\begin{array}{l}{V}_m{x}_0\end{array}\right].\end{array}

The following theorem will give the solution to the minimum value problem (8).

Theorem 2.2　 Assuming that the m-step incomplete orthogonalization process has been performed, the form of x_m in problem (8) is as in (3). Let {\left\{{\lambda }_i,{\mu }_i\right\}}_{i=1,2,\dots ,m+1} be the combination of all generalized eigenvalues and corresponding eigenvectors of \left(L_m^{\mathrm{T}}{L}_m,G_m^{\mathrm{T}}{G}_m\right), where {\lambda }_i⩾{\lambda }_{i+1},i=1,2,\dots ,m. If the last element of the vector {\mu }_{m+1} is not 0, i.e., {\mu }_{m+1,m+1}\ne 0, then the solution y_m of the minimum value problem (8) satisfies

(10)\left[\begin{array}{c}{y}_m\\ 1\end{array}\right]=\frac{1}{{\mu }_{m+1,m+1}}{\mu }_{m+1}

and

(11){\Vert {\mathrm{\Delta }}_{min}\Vert }_F⩽\sqrt{m+1}{\lambda }_{m+1}^{\frac{1}{2}}.

Proof　The preceding reasoning process eventually transforms problem (2) into the minimum value problem (8). By

\underset{y_m\in {\mathbb{R}}^m}{min}\frac{{\Vert H_m{y}_m-\beta e_1\Vert }_2^2}{{\Vert x_m\Vert }_2^2}=\underset{y_m\in {\mathbb{R}}^m}{min}\frac{{\Vert \left[\begin{array}{l}{H}_m-\beta e_1\end{array}\right]{\left[\begin{array}{l}{y}_m^{\mathrm{T}}1\end{array}\right]}^{\mathrm{T}}\Vert }_2^2}{{\Vert \left[\begin{array}{l}{V}_m{x}_0\end{array}\right]{\left[\begin{array}{l}{y}_m^{\mathrm{T}}1\end{array}\right]}^{\mathrm{T}}\Vert }_2^2}=\underset{\mu \in {\mathbb{R}}^{m+1}}{min}\frac{{\Vert L_m\mu \Vert }_2^2}{{\Vert G_m\mu \Vert }_2^2},

and the Courant−Fischer theorem [24], we have,

(12)\underset{\mu \in {\mathbb{R}}^{m+1}}{min}\frac{{\Vert L_m\mu \Vert }_2^2}{{\Vert G_m\mu \Vert }_2^2}={\lambda }_{m+1},

where {\lambda }_{m+1} is the minimum generalized eigenvalue of \left(L_m^{\mathrm{T}}{L}_m,G_m^{\mathrm{T}}{G}_m\right). Since {\mu }_{m+1,m+1}\ne 0, the solution y_m of (8) can be calculated by using equation (10). From (7), we soon obtain (11).

By Theorem 2.2, we derive the IGMBACK algorithm. In order to reduce the storage and computation capacity of the algorithm, the algorithm in this paper adopts the m step loop format. The main steps for \mathrm{I}\mathrm{G}\mathrm{M}\mathrm{B}\mathrm{A}\mathrm{C}\mathrm{K}(m) are given as follows.

Algorithm 1　Restart \mathrm{I}\mathrm{G}\mathrm{M}\mathrm{B}\mathrm{A}\mathrm{C}\mathrm{K}(m):

① Choose the initial estimate x_0 of (1) and the parameter q, where q satisfies 2⩽q⩽m, compute r_0:=b-A{x}_0 and v_1:=\frac{r_0}{\beta }, where \beta :={\Vert r_0\Vert }_2.

② Incomplete orthogonalization process (q) : Perform m steps of the incomplete orthogonalization process to compute V_{m+1} and H_m, for j= 1,2,\dots ,m.

i. Calculate {\hat{v}}_{j+1}:=A{v}_j.

ii. For i=max\{1,j-q+1\},\dots ,j-1,j, compute h_{ij}=v_i^{\mathrm{T}}{\hat{v}}_{j+1},{\hat{v}}_{j+1}:={\hat{v}}_{j+1}-h_{ij}{v}_i.

iii. Calculate h_{j+1,j}={\Vert {\hat{v}}_{j+1}\Vert }_2. If h_{j+1,j}\ne 0, compute v_{j+1}:=\frac{{\hat{v}}_{j+1}}{h_{j+1,j}}, otherwise stop.

③ Solving generalized eigenvalue problems.

(13)L_m^{\mathrm{T}}{L}_{m}u=\lambda G_m^{\mathrm{T}}{G}_{m}u,

where L_m and G_m are defined in (9).

④ The approximate solution x_m^{\mathrm{I}\mathrm{G}\mathrm{B}}:=x_0+V_m{y}_m is calculated from Equation (10) in Theorem 2.2.

⑤ Start over: Define \Vert \left(H_m{y}_m-\beta e_1\right)\Vert x_m{\Vert {}_2^{-2}{x}_m^{\mathrm{T}}\Vert }_F={\lambda }_{m+1}^{\frac{1}{2}}. If the termination condition is met, stop, otherwise redefine x_0:=x_m^{\mathrm{I}\mathrm{G}\mathrm{B}}, compute r_0:=b-A{x}_0,v_1:=\frac{r_0}{\beta }, where \beta :={\Vert r_0\Vert }_2, and return to step (2).

Denote

(14)P=L_m^{\mathrm{T}}{L}_m=\left[\begin{array}{cc}{H}_m^{\mathrm{T}}{H}_m& -\beta H_m^{\mathrm{T}}{e}_1\\ -\beta e_1^{\mathrm{T}}{H}_m& {\beta }^2\end{array}\right],Q=G_m^{\mathrm{T}}{G}_m=\left[\begin{array}{cc}{V}_m^{\mathrm{T}}{V}_m& V_m^{\mathrm{T}}{x}_0\\ x_0^{\mathrm{T}}{V}_m& x_0^{\mathrm{T}}{x}_0\end{array}\right].

Therefore, the generalized eigenvalue problem (13) can then be written as Pu=\lambda Qu.

The minimum value problem (12) can be rewritten as

(15){\lambda }_{m+1}^{\frac{1}{2}}=\underset{\mu \in {\mathbb{R}}^{m+1}}{min}{\left(\frac{{\Vert G_m\mu \Vert }_2}{{\Vert \mu \Vert }_2}\right)}^{-1}\frac{{\Vert L_m\mu \Vert }_2}{{\Vert \mu \Vert }_2}.

Using equation (15), we can give an estimate of the range of \Vert \left(H_m{y}_m-\beta e_1\right)\Vert x_m{\Vert {}_2^{-2}{x}_m^{\mathrm{T}}\Vert }_F={\lambda }_{m+1}^{\frac{1}{2}} in the termination condition of the algorithm. This is shown by the following theorem. Here x_m=x_0+V_m{y}_m is the approximate solution of the IGMBACK algorithm.

Theorem 2.3　 Assuming that m steps of incomplete orthogonalization have been performed, {\sigma }_1\left(G_m\right) denotes the maximum singular value of G_m. {\sigma }_{m+1}\left(L_m\right),{\sigma }_{m+1}\left(G_m\right) are the minimum singular values of L_m,G_m respectively. Let {\sigma }_{m+1}\left(G_m\right)>0, we have

{\sigma }_{m+1}\left(L_m\right){\sigma }_1^{-1}\left(G_m\right)⩽{\lambda }_{m+1}^{\frac{1}{2}}⩽{\sigma }_{m+1}\left(L_m\right){\sigma }_{m+1}^{-1}\left(G_m\right).

Proof　By employing equation (15), we have

{\lambda }_{m+1}^{\frac{1}{2}}⩾\underset{z_1\in {\mathbb{R}}^{m+1}}{min}{\left(\frac{{\Vert G_m{z}_1\Vert }_2}{{\Vert z_1\Vert }_2}\right)}^{-1}\underset{z_2\in {\mathbb{R}}^{m+1}}{min}\frac{{\Vert L_m{z}_2\Vert }_2}{{\Vert z_2\Vert }_2}={\sigma }_1^{-1}\left(G_m\right){\sigma }_{m+1}\left(L_m\right)

and

{\lambda }_{m+1}^{\frac{1}{2}}⩽\underset{z_1\in {\mathbb{R}}^{m+1}}{max}{\left(\frac{{\Vert G_m{z}_1\Vert }_2}{{\Vert z_1\Vert }_2}\right)}^{-1}\underset{z_2\in {\mathbb{R}}^{m+1}}{min}\frac{{\Vert L_m{z}_2\Vert }_2}{{\Vert z_2\Vert }_2}={\sigma }_{m+1}^{-1}\left(G_m\right){\sigma }_{m+1}\left(L_m\right).

Thus, the conclusion of the theorem holds.□

2.3 The theoretical research of IGMBACK algorithm

This section presents a theoretical study of the IGMBACK algorithm in terms of finite termination, existence and uniqueness of the solution.

2.3.1 Finite termination of the algorithm

Suppose that after m steps of incomplete orthogonalization, we get h_{m+1,m}=0, i.e., P is singular, then the vector v_{m+1} cannot be formed. When this happens, we call it a lucky break [22].

Let the column vectors of H_{m-1} be linearly independent and the F- backward perturbation matrix with minimum norm {\mathrm{\Delta }}_i generated in the IGMBACK algorithm satisfies {\mathrm{\Delta }}_i\ne 0 when . According to Theorem 2.1, a sufficient condition for {\mathrm{\Delta }}_{min}={\mathrm{\Delta }}_m=0 is H_m{y}_m-\beta e_1=0, then we have

{\mathrm{\Delta }}_{min}=0\iff h_{m+1,m}=0.

In fact, assume h_{m+1,m}=0, y_m={\overline{H}}_m^{-1}\beta e_1, then {\mathrm{\Delta }}_{min}=0, thus x_m=x_0+V_m{\overline{H}}_m^{-1}\beta e_1 is an exact solution to the system of equations (1). Conversely, if {\mathrm{\Delta }}_{\text{min}}=0, then H_m{y}_m-\beta e_1=0. Let y_m={\left[{\left(y_m^{(1)}\right)}^{\mathrm{T}},y_m^{(2)}\right]}^{\mathrm{T}}, where y_m^{(2)}\in \mathbb{R}. From the last line of the equation H_m{y}_m-\beta e_1=0, we get h_{m+1,m}{y}_m^{(2)}=0. If y_m^{(2)}=0, then H_{m-1}{y}_m^{(1)}-\beta e_1=0, and thus {\mathrm{\Delta }}_{m-1}=0, which contradicts the previous assumption and gives h_{m+1,m}=0.

We know that the algorithm can execute at most m=n steps before a lucky break occurs, so the following corollary holds.

Corollary 2.1　 For any A\in {\mathbb{R}}^{n\times n} and b\in {\mathbb{R}}^n, the IGMBACK algorithm can converge by up to n steps.

According to the above discussion, h_{i+1,i}\ne 0,i⩽m is assumed in the rest of this paper. That is, P is assumed to be non-singular.

2.3.2 Existence and uniqueness of solutions

There are two cases where the solution of the algorithm does not exist: one is when the generalized eigenvalue problem Pu=\lambda Qu has degenerate generalized eigenvalues, the other is when {\mu }_{m+1,m+1}=0 in (10).

Lemma 2.2　 Let the incomplete orthogonalization process in m steps have been performed, then, the necessary and sufficient condition for Q=G_m^{\mathrm{T}}{G}_m to be a singular matrix is x_0=0 or x_0\in {\mathcal{K}}_m\left(A,r_0\right).

Proof　According to the basic knowledge of algebra, Q=G_m^{\mathrm{T}}{G}_m is a singular matrix if and only if the rank m+1 of matrix G_m, where G_m is shown in (9). On the one hand, if Q=G_m^{\mathrm{T}}{G}_m is a singular matrix, then , i.e., the set of column vectors of G_m=\left[\begin{array}{l}{V}_m{x}_0\end{array}\right] is linearly related. The set of column vectors of V_m is linearly independent, so x_0 can be uniquely linearly tabulated by the set of column vectors of V_m, i.e., x_0\in {\mathcal{K}}_m\left(A,r_0\right). On the other hand, if x_0=0 or x_0\in {\mathcal{K}}_m\left(A,r_0\right), i.e., x_0 can be linearly tabulated by the set of column vectors of V_m, then the set of column vectors of G_m is linearly related. Thus , i.e., Q is a singular matrix.

By Lemma 2.2, Q is potentially singular. We choose x_0 in the actual implementation of the IGMBACK algorithm such that x_0\notin {\mathcal{K}}_m\left(A,r_0\right), which ensures that Q is a non-singular matrix. Thus, the generalized eigenvalue problem Pu=\lambda Qu belongs to the non-degenerate generalized eigenvalue problem [7], without degenerate generalized eigenvalues and eigenvectors. In general, the dimension m of {\mathcal{K}}_m\left(A,r_0\right) chosen in the IGMBACK algorithm is much smaller than the order n of A. Therefore, there is plenty of scope for choosing x_0 in the n-dimensional vector space {\mathbb{R}}^n such that x_0\notin {\mathcal{K}}_m\left(A,r_0\right). Throughout the rest of the paper, we assume that x_0\notin {\mathcal{K}}_m\left(A,r_0\right), i.e., that Q is a non-singular matrix.

The following theorem gives a sufficient condition for the IGMBACK algorithm to fail to form an approximate solution when the multiplicity of {\lambda }_{m+1} is 1.

Theorem 2.4　 Suppose \left\{{\lambda }_{m+1},{\mu }_{m+1}\right\} is the combination of the minimum generalized eigenvalue and the corresponding eigenvector of (P,Q), Q is a non-singular matrix, and the multiplicity of {\lambda }_{m+1} is 1. Let {\mu }_{m+1}={\left[ll{\hat{\mu }}^{\mathrm{T}}{\mu }_{m+1,m+1}\right]}^{\mathrm{T}}. Then

{\mu }_{m+1,m+1}=0\iff H_m^{\mathrm{T}}{H}_m\hat{\mu }={\lambda }_{m+1}{V}_m^{\mathrm{T}}{V}_m\hat{\mu },\hat{\mu }\perp \left({\lambda }_{m+1}{V}_m^{\mathrm{T}}{x}_0+H_m^{\mathrm{T}}\beta e_1\right).

Proof　 \Rightarrow Let {\mu }_{m+1,m+1}=0. Substituting the expressions for P,Q in (14) into Pu=\lambda Qu, we have

\left[\begin{array}{cc}{H}_m^{\mathrm{T}}{H}_m& -\beta H_m^{\mathrm{T}}{e}_1\\ -\beta e_1^{\mathrm{T}}{H}_m& {\beta }^2\end{array}\right]\left[\begin{array}{c}\hat{\mu }\\ 0\end{array}\right]={\lambda }_{m+1}\left[\begin{array}{cc}{V}_m^{\mathrm{T}}{V}_m& V_m^{\mathrm{T}}{x}_0\\ x_0^{\mathrm{T}}{V}_m& x_0^{\mathrm{T}}{x}_0\end{array}\right]\left[\begin{array}{c}\hat{\mu }\\ 0\end{array}\right].

From the first line of the matrix equation above, we get H_m^{\mathrm{T}}{H}_m\hat{\mu }={\lambda }_{m+1}{V}_m^{\mathrm{T}}{V}_m\hat{\mu }. From the second line we get -\beta e_1^{\mathrm{T}}{H}_m\hat{\mu }= {\lambda }_{m+1}{x}_0^{\mathrm{T}}{V}_m\hat{\mu }, i.e., ({\lambda }_{m+1}{V}_m^{\mathrm{T}}{x}_0+H_m^{\mathrm{T}}\beta e_1{)}^{\mathrm{T}}\hat{\mu }=0, that is,

\hat{\mu }\perp \left({\lambda }_{m+1}{V}_m^{\mathrm{T}}{x}_0+H_m^{\mathrm{T}}\beta e_1\right).

\Leftarrow Assume H_m^{\mathrm{T}}{H}_m\hat{\mu }={\lambda }_{m+1}{V}_m^{\mathrm{T}}{V}_m\hat{\mu }, and {\left({\lambda }_{m+1}{V}_m^{\mathrm{T}}{x}_0+H_m^{\mathrm{T}}\beta e_1\right)}^{\mathrm{T}}\hat{\mu }=0. Then {\left[{\hat{\mu }}^{\mathrm{T}},0\right]}^{\mathrm{T}} is the eigenvector of Pu=\lambda Qu, and the corresponding generalized eigenvalue is {\lambda }_{m+1}. Since Q is a non-singular matrix, Pu=\lambda Qu does not have degenerate generalized eigenvalues and eigenvectors, and since the multiplicity of {\lambda }_{m+1} is 1, k{\left[{\hat{\mu }}^{\mathrm{T}},0\right]}^{\mathrm{T}}(k\in \mathbb{R},k\ne 0) is the eigenvector of all Pu=\lambda Qu corresponding to {\lambda }_{m+1}. Therefore, there must be {\mu }_{m+1,m+1}=0.

The existence and uniqueness of the solution are discussed further below. For convenience, according to (14) and (8), denote

\begin{array}{rl}& P_m=P=L_m^{\mathrm{T}}{L}_m,Q_m=Q=G_m^{\mathrm{T}}{G}_m=\left[\begin{array}{l}{V}_m^{\mathrm{T}}{V}_m{V}_m^{\mathrm{T}}{x}_0\\ x_0^{\mathrm{T}}{V}_m{x}_0^{\mathrm{T}}{x}_0\end{array}\right],\\ & {\tilde{\mathrm{\Delta }}}_m=\underset{y_m\in {\mathbb{R}}^m}{min}\frac{\left[\begin{array}{l}{y}_m^{\mathrm{T}}1\end{array}\right]P_m{\left[\begin{array}{l}{y}_m^{\mathrm{T}}1\end{array}\right]}^{\mathrm{T}}}{\left[\begin{array}{l}{y}_m^{\mathrm{T}}1\end{array}\right]Q_m{\left[\begin{array}{l}{y}_m^{\mathrm{T}}1\end{array}\right]}^{\mathrm{T}}}.\end{array}

Let’s start with a lemma.

Lemma 2.3　 Suppose that an incomplete orthogonalization process of m+1 steps has been performed. {\left\{{\lambda }_i\right\}}_{i=1,2,\dots ,m+1} and {\left\{{\hat{\lambda }}_i\right\}}_{i=1,2,\dots ,m+2} represent all generalized eigenvalues of \left(P_m,Q_m\right) and \left(P_{m+1},Q_{m+1}\right) in non-increasing order, respectively, then

{\hat{\lambda }}_i⩾{\lambda }_i⩾{\hat{\lambda }}_{i+1},i=1,2,\dots ,m+1.

Proof　Let {\left\{{\lambda }_i,{\mu }_i\right\}}_{i=1,2,\dots ,m+1} be all generalized characteristic pairs of \left(P_m,Q_m\right). According to Subsection 2.3.1, if h_{m+2,m+1}\ne 0, then L_{m+1},L_m are non-singular matrices. Therefore, we can write Pu=\lambda Qu as {\left(L_m^{-1}\right)}^{\mathrm{T}}{Q}_m{L}_m^{-1} z_i={\xi }_i{z}_i, where z_i=L_m{\mu }_i,{\xi }_i is the reciprocal of {\lambda }_i. The size order of {\left\{{\lambda }_i\right\}}_{i=1,2,\dots ,m+1} determines \cdots ⩽{\xi }_{m+1}. By

K=\left[\begin{array}{c}{H}_{m+1}-\beta e_1\end{array}\right]=\left[\begin{array}{cc}{L}_m& {\hat{h}}_{m+1}\\ 0& h_{m+2,m+1}\end{array}\right]P_{m+2}(m+1,m+2),

where {\hat{h}}_{m+1}=\left[\begin{array}{l}{I}_{m+1}0\end{array}\right]H_{m+1}{e}_{m+1},P_{m+2}(m+1,m+2) represents the primary matrix after interchanging the m+2 and m+1 columns of the unit matrix I_{m+2}, one gets that

L_{m+1}^{-1}=P_{m+2}(m+1,m+2)\left[\begin{array}{cc}{L}_m^{-1}& -L_m^{-1}{\hat{h}}_{m+1}{h}_{m+2,m+1}^{-1}\\ 0& h_{m+2,m+1}^{-1}\end{array}\right].

Substituting the above equation into {\left(L_{m+1}^{-1}\right)}^{\mathrm{T}}{Q}_{m+1}{L}_{m+1}^{-1}, we have

{\left(L_{m+1}^{-1}\right)}^{\mathrm{T}}{Q}_{m+1}{L}_{m+1}^{-1}=\left[\begin{array}{cc}{\left(L_m^{-1}\right)}^{\mathrm{T}}{Q}_m{L}_m^{-1}& f\\ f^{\mathrm{T}}& g\end{array}\right],

where

\begin{array}{rl}& g=h_{m+2,m+1}^{-2}+2h_{m+2,m+1}^{-1}{d}_{m+1}{h}_{m+1}+d_{m+1}{Q}_m{d}_{m+1}^{\mathrm{T}},\\ & f=h_{m+2,m+1}^{-1}{\left(L_m^{-1}\right)}^{\mathrm{T}}\left(h_{m+1}-Q_m{L}_m^{-1}{\hat{h}}_{m+1}\right),\\ & h_{m+1}={\left[\begin{array}{l}{V}_m^{\mathrm{T}}{v}_{m+1}\\ x_0^{\mathrm{T}}{v}_{m+1}\end{array}\right]}_{(m+1)\times 1},d_{m+1}=-h_{m+2,m+1}^{-1}{\hat{h}}_{m+1}^{\mathrm{T}}{\left(L_m^{-1}\right)}^{\mathrm{T}}\end{array}.

Thus, according to Theorem IV.4.2 in reference [24], one obtains that

{\hat{\xi }}_{m+2}⩾{\xi }_{m+1}⩾{\hat{\xi }}_{m+1}⩾{\xi }_m⩾\cdots ⩾{\hat{\xi }}_2⩾{\xi }_1⩾{\hat{\xi }}_1,

where {\left\{{\hat{\xi }}_i\right\}}_{i=1,2,\dots ,m+2} is all the eigenvalues of {\left(L_{m+1}^{-1}\right)}^{\mathrm{T}}{Q}_{m+1}{L}_{m+1}^{-1}. Since {\xi }_i,{\hat{\xi }}_i are the reciprocals of {\lambda }_i and {\hat{\lambda }}_i, respectively, the conclusion of the lemma holds.

According to Lemma 2.3, we obtain the following conclusions about the convergence, existence and uniqueness of the solution.

Corollary 2.2　 Assuming that the incomplete orthogonalization process has been performed in m+1 steps, the following conclusion holds:

(a) {\tilde{\mathrm{\Delta }}}_{m+1}⩽{\tilde{\mathrm{\Delta }}}_m.

(b) If and {\mu }_{m+2,m+2}\ne 0, then there is a unique approximate solution to the IGMBACK algorithm.

(c) If {\hat{\lambda }}_{m+2}={\hat{\lambda }}_{m+1}, and the approximate solution of the IGMBACK algorithm exists, then the approximate solution of the IGMBACK algorithm is not unique, and {\tilde{\mathrm{\Delta }}}_{m+1}={\tilde{\mathrm{\Delta }}}_m, i.e., the IGMBACK algorithm stalls from step m to step m+1.

Let the minimum generalized eigenvalue {\lambda }_{m+1} of \left(P_m,Q_m\right) be of multiplicity not 1. Assume {\lambda }_{m-k+1}>{\lambda }_{m-k+2}=\cdots ={\lambda }_{m+1}, i.e., let {\lambda }_{m+1} be of k multiplicity, and the corresponding standard orthogonal eigenvectors are {\mu }_{m-k+2},{\mu }_{m-k+3},\dots ,{\mu }_{m+1}, respectively. Let

U_k=\left[{\mu }_{m-k+2}{\mu }_{m-k+3}\dots {\mu }_{m+1}\right]\in {\mathbb{R}}^{(m+1)\times k}.

Then for any z={\left[\begin{array}{l}{z}_1^{\mathrm{T}}{z}_2\end{array}\right]}^{\mathrm{T}}\in \mathrm{s}\mathrm{p}\mathrm{a}\mathrm{n}\left\{U_k\right\},z_2\in \mathbb{R}, if z_2\ne 0, the approximate solution of the IGMBACK algorithm takes the form

x_m=x_0+V_m\frac{z_1}{z_2}.

In the practical implementation of the algorithm, we often choose \frac{z_1}{z_2} such that {\Vert z_1\Vert }_2{\Vert }_2 is minimal to make the resulting approximate solution unique. Maximize \mathrm{s}\mathrm{p}\mathrm{a}\mathrm{n}\left\{U_k\right\} in all unit vectors containing span \{z:z={\left[\begin{array}{ll}{z}_1^{\mathrm{T}}& z_2\end{array}\right]}^{\mathrm{T}}\in \mathrm{s}\mathrm{p}\mathrm{a}\mathrm{n}\left\{U_k\right\},\Vert z{\Vert }_2=1\} maximizes \left|z_2\right| such that {\Vert \frac{z_1}{z_2}\Vert }_2 is minimized.

3 Execution of the IGMBACK algorithm

This section considers how to solve the generalized eigenvalue problem (13), i.e., Pu=\lambda Qu.

In the implementation of Algorithm 1, we focus on solving the generalized eigenvalue problem (13). Based on Theorem 2.2, it is sufficient to compute the minimum generalized eigenvalue of \left(L_m^{\mathrm{T}}{L}_m,G_m^{\mathrm{T}}{G}_m\right) and the corresponding combination of eigenvectors \left\{{\lambda }_{m+1},{\mu }_{m+1}\right\}. Since L_m^{\mathrm{T}}{L}_m-\lambda G_m^{\mathrm{T}}{G}_m is symmetric, we can compute \left\{{\lambda }_{m+1},{\mu }_{m+1}\right\} with the help of the inverse iteration method. However, when the condition numbers of L_m^{\mathrm{T}}{L}_m and G_m^{\mathrm{T}}{G}_m are relatively large, it becomes difficult to compute \left\{{\lambda }_{m+1},{\mu }_{m+1}\right\} by inverse iteration. Therefore, in the GMBACK algorithm [13], the authors solve the corresponding generalized eigenvalue problem by solving the associated singular value decomposition problem, which overcomes the difficulties associated with the increasing number of conditions for L_m^{\mathrm{T}}{L}_m and G_m^{\mathrm{T}}{G}_m. Can Algorithm 1 also calculate \left\{{\lambda }_{m+1},{\mu }_{m+1}\right\} by solving the associated singular value decomposition problem?

Noting that the matrix Q contains V_m^{\mathrm{T}}{V}_m, make the Cholesky decomposition of V_m^{\mathrm{T}}{V}_m: V_m^{\mathrm{T}}{V}_m=X{X}^{\mathrm{T}}, where X is the lower triangular matrix, thus

\begin{array}{rl}Q& =\left[\begin{array}{cc}X& 0\\ x_0^{\mathrm{T}}{V}_m{\left(X^{-1}\right)}^{\mathrm{T}}& \sqrt{{\Vert x_0\Vert }_2^2-{\Vert X^{-1}{V}_m^{\mathrm{T}}{x}_0\Vert }_2^2}\end{array}\right]\left[\begin{array}{cc}{X}^{\mathrm{T}}& X^{-1}{V}_m^{\mathrm{T}}{x}_0\\ 0& \sqrt{{\Vert x_0\Vert }_2^2-{\Vert X^{-1}{V}_m^{\mathrm{T}}{x}_0\Vert }_2^2}\end{array}\right]\\ & =Z{Z}^{\mathrm{T}},\end{array}

where

Z=\left[\begin{array}{cc}X& 0\\ x_0^{\mathrm{T}}{V}_m{\left(X^{-1}\right)}^{\mathrm{T}}& \sqrt{{\Vert x_0\Vert }_2^2-{\Vert X^{-1}{V}_m^{\mathrm{T}}{x}_0\Vert }_2^2}\end{array}\right].

By substituting Q=Z{Z}^{\mathrm{T}} into Pu=\lambda Qu, the original generalized eigenvalue problem Pu=\lambda Qu is transformed into an ordinary eigenvalue problem for symmetric matrices:

Z^{-1}P{\left(Z^{-1}\right)}^{\mathrm{T}}v=\lambda v,

where v=Z^{\mathrm{T}}u. When \sqrt{\lambda } is small, based on P={\left[H_m-\beta e_1\right]}^{\mathrm{T}}\left[H_m-\beta e_1\right], to solve exactly for the matrix Z^{-1}P{\left(Z^{-1}\right)}^{\mathrm{T}} with the minimum eigenvalue {\lambda }_{min} and the corresponding eigenvector v is not an easy task. Therefore, we transform the original problem (16) into a solution for the minimum singular value of another matrix

(16)\begin{array}{l}\begin{array}{l}{={\left[\begin{array}{l}{H}_m{\left(X^{-1}\right)}^{\mathrm{T}}-H_m([V_m^{\mathrm{T}}{V}_m\end{array}\right]}^{-1})}^{\mathrm{T}}{V}_m^{\mathrm{T}}{x}_0\cdot xl-\beta \cdot xl\cdot e_1]\end{array},\end{array}

and the corresponding right singular vector, which solves the original problem very well. In equation (17), xl=\frac{1}{\sqrt{{\Vert x_0\Vert }_2^2-{\Vert X^{-1}{V}_m^T{x}_0\Vert }_2^2}}.

So we get the specific execution process of IGMBACK algorithm.

Algorithm 2　Restart IGMBACK algorithm, \mathrm{I}\mathrm{G}\mathrm{M}\mathrm{B}\mathrm{A}\mathrm{C}\mathrm{K}(m):

① Select initial estimated value x_0 and parameter q of (1), where q satisfies 2⩽q⩽m. Calculate r_0:=b-A{x}_0 and v_1:=\frac{r_0}{\beta }, where \beta :={\Vert r_0\Vert }_2.

② Incomplete orthogonalisation process (q): Performing an incomplete orthogonalisation process in m steps, V_{m+1} and H_m were calculated, for j=1,2,\dots ,m.

i. Calculate {\hat{v}}_{j+1}:=A{v}_j.

ii. For i=max\{1,j-q+1\},\dots ,j, compute h_{ij}=v_i^T{\hat{v}}_{j+1}, {\hat{v}}_{j+1}:={\hat{v}}_{j+1}-h_{ij}{v}_i.

iii. Calculate h_{j+1,j}={\Vert {\hat{v}}_{j+1}\Vert }_2, if h_{j+1,j}\ne 0, calculate v_{j+1}=\frac{{\hat{v}}_{j+1}}{h_{j+1,j}}, otherwise stop.

③ Compute the minimum singular value \sigma of the matrix [H_m-\beta e_1](Z^{-1}{)}^T and the corresponding right singular vector v. Calculate u=(Z^{-1}{)}^{T}v=[{\tilde{y}}_m^T,\eta {]}^T. Standardize the vector u to derive y_m:y_m=\frac{{\tilde{y}}_m}{\eta }. Forming an approximate solution x_m^{IGB}:=x_0+V_m{y}_m.

④ Starting over: Define \Vert \left(H_m{y}_m-\beta e_1\right)\Vert x_m{\Vert {}_2^{-2}{x}_m^{\mathrm{T}}\Vert }_F^2=\sigma. If the termination condition is met, stop. Otherwise, redefine x_0:=x_m^{IGB} and calculate r_0:=b-A{x}_0, v_1:=\frac{r_0}{\beta } where \beta :={\Vert r_0\Vert }_2, and return to step ②.

4 Numerical experiments

It is well known that when the coefficient matrix A is a ill-posed matrix, a small perturbation of A or b causes a large change in the approximate solution of (1). This is further exacerbated [9] when b is parallel to the left singular vector corresponding to the minimum singular value of A. In this section, we illustrate the effectiveness of the ICMBACK(m) algorithm in solving the preceding problem by means of several numerical examples, and the IGMBACK(m) algorithm is often more efficient than the GMBACK(m) and CMRES(m) algorithms in solving many problems. The programming software used for the numerical experiments in this paper is MATLAB 7.0.

We mainly compare the following three algorithms:

1) IGMBACK(m,q)(q is the parameter in the incomplete orthogonalisation process).

2) GMBACK(m).

3) CMRES(m). The different convergence rates of the three algorithms during execution are reflected by comparing the change in the backward perturbation norm {\Vert {\mathrm{\Delta }}_{min}\Vert }_F of the approximate solutions of the three algorithms with the increase in the number of selected generations (restarts). In the GMRES(m) algorithm, let the backward perturbation norm {\Vert {\mathrm{\Delta }}_{min}\Vert }_F=\frac{{\Vert r_m\Vert }_2}{{\Vert x_m\Vert }_2}. The stop criteria is {\Vert {\mathrm{\Delta }}_{min}\Vert }_F⩽{10}^{-7}.

4.1 Numerical experiments with strong informal matrices

This part of the numerical experiment will consider two matrices that are sensitive to small perturbations (called sensitivity matrices).

We notate the singular value decomposition of the matrix A with the following symbols: A=US{V}^T, where

U=[u_1,u_2,\dots ,u_n],V=[v_1,v_2,\dots ,v_n],S=(s_1,s_2,\dots ,s_n).

For i=1,2,\dots ,n-1, s_i⩾s_{i+1} in S, then u_1 and v_1 are the left singular vectors and the right singular vector corresponding to the maximum singular value s_n of matrix A, respectively.

Example 4.1　Consider a Toeplitz matrix of order 100\times 100: A=\mathrm{T}\mathrm{o}\mathrm{e}\mathrm{p}\mathrm{l}\mathrm{i}\mathrm{t}\mathrm{z}([\underset{\_}{1},1,1]). Calculate the maximum singular value s_1=2.9990, and minimum singular value s_n=2.6784\cdot {10}^{-2} of the matrix A. For the first experiment Test 1, set b=u_1, x_0=v_1. Thus, we obtain r_0=u_1-A{v}_1, that is to say, r_0 is parallel to u_1 (noted as r_0\Vert u_1). For the second experiment Test 2, assume b=u_n, x_0=v_n. Thus, we obtain r_0=u_n-A{v}_n, that is to say, r_0\Vert u_n. In these two experiments, let m=20, q=15. The experimental results are shown in Fig.1 (a) and Fig.1 (b), respectively. -10 in the vertical coordinate represents {10}^{-10}, and so on.

Fig.1 (a) shows that the convergence rates of the three algorithms in the first experiment are essentially the same. From Fig.1 (b), it can be seen that the convergence rates of IGMBACK(m) and GMBACK(m) in the second experiment are equivalent (IGMBACK (m) converges slightly faster), and both converge significantly faster than GMRES(m). The coefficient matrices are the same in both experiments (both sensitive matrices), what makes IGMBACK(m) and GMBACK(m) significantly better than GMRES(m) in the second experiment? This is mainly because the second experiment was set up with b\Vert u_n and x_0=v_n making r_0\Vert u_n.

Example 4.2　Consider another matrix of order 100\times 100: A=\mathrm{T}\mathrm{o}\mathrm{e}\mathrm{p}\mathrm{l}\mathrm{i}\mathrm{t}\mathrm{z}([-1,\underset{\_}{1},1,1,1]). Matrix A has sensitive eigenvalues. Calculate the maximum singular value s_1=4.2394, and minimum singular value s_n=9.0205\cdot {10}^{-1} of the matrix A. The first experiment Test 1, solves the equation Ax=b, with a setting of b=u_0, x_0=v_n. Let m=20, q=15. The result is shown in Fig.2 (a).

From Fig.2 (a), it can be seen that the convergence rates of IGMBAC-K(m), GMBACK(m), and GMRES(m) are completely consistent, that is, the setting of vector b and the selection of initial estimates for approximate solutions cannot cause differences in the convergence rates of the three algorithms, as shown in Example 4.1 (even if set at b\Vert u_n). This is mainly owing to the fact that the coefficient matrix A is well conditioned. However, by reducing the minimum singular value of A, the convergence rates of the above three algorithms will become different. Use E=0.89u_n{v}_n^T to perturb the matrix A. The first n-1 singular values of the perturbed coefficient matrix A-E are the same as A, while the minimum singular value S_n decreases to 1.2048\cdot {10}^{-2}. The left and right singular vectors also remain unchanged. The second experiment Test 2 considers solving the system of equations (A-E)x=b using the three algorithms described above, still setting b=u_n, x_0=v_n. Let m=20, q=15. The result is shown in Fig.2 (b).

As depicted in Fig.2 (b), stagnation occurs in GMRES, while IGMBACK(m) and GMBACK(m) are convergent, and \mathrm{I}\mathrm{G}\mathrm{M}\mathrm{B}\mathrm{A}\mathrm{C}\mathrm{K}(m) converges at a rate comparable to GMBACK(m). Obviously, the latter two algorithms are clearly superior to the restarted GMRES(m), which is largely attributable to the the foundation of setting b\Vert u_n. In the second experiment, when the minimum singular value of A is reduced, the condition of the coefficient matrix becomes worse and the coefficient matrix becomes sensitive matrix.

From the previous two examples, we can conclude that the GMRES(m) algorithm, which uses the least-squares problem to derive approximate solutions, is not sufficient to guarantee convergence if the coefficient matrix is a sensitive matrix and the vector b on the right-hand side of the system of equations is parallel to the left singular vector corresponding to the least singular value of the coefficient matrix. However, GMBACK(m) and IGMBACK(m) algorithms can get better convergence effect by minimizing or approximately minimizing the backward perturbation norm. Interestingly, when comparing the TLS (Total Least Sguares) algorithm with the classical least square method, we can come to a similar conclusion [26].

In order to further illustrate that the algorithms IGMBACK(m) can be compared with GMBACK(m), we compare the CPU time taken by each of the two algorithms to run in Example 4.1 and Example 4.2. The comparison results are shown in Tab.1. In the table, CPU and CPU{}_a represent the entire CPU time used and the average CPU time used every reboot, respectively, in seconds. iter represents the number of restarts; ratio is the ratio of CPU time used by algorithm IGMBACK(m) and GMBACK(m); BPN represents the norm of backward perturbation. Note that IGMBACK(m) becomes GMBACK(m) when q=m.

From Tab.1, we can see that the algorithm \mathrm{I}\mathrm{G}\mathrm{M}\mathrm{B}\mathrm{A}\mathrm{C}\mathrm{K}(m) is more efficient than GMBACK(m).

4.2 Numerical experiments related to partial differential equations

In practical applications, a large number of large-scale nonsymmetric linear systems arise from the discretization of partial differential equations.

Example 4.3　Consider the convective diffusion equation defined over the region (0,1) × (0,1)

-\frac{{\mathrm{\partial }}^{2}u}{\mathrm{\partial }{x}^2}-\frac{{\mathrm{\partial }}^{2}u}{\mathrm{\partial }{y}^2}+\gamma \left(x\frac{\mathrm{\partial }u}{\mathrm{\partial }x}+y\frac{\mathrm{\partial }u}{\mathrm{\partial }y}\right)+\beta u=f,

where \gamma =1000, \beta =10. The boundary condition is u(x,y)=0. The above convective diffusion equation is discretized by the central difference method with lattice length h=\frac{1}{32} to obtain an nonsymmetric matrix A of order {31}^2. Take b such that x=[11\cdots 1{]}^{\mathrm{T}} is an exact solution to the system of equations Ax=b. Suppose x_0 is a random matrix of order 961\times 1, with the values of the elements in the interval (0.0,1.0), and let m=15, q=10. The experimental result is shown in Fig.3.

In this example, the GMRES algorithm does not converge in 400 steps when m is set at 15, while Fig.3 shows that the ICMBACK and GMBACK can drop to around {10}^{-8} in 40 steps by the backward perturbation norm. This indicates that the latter two algorithms are significantly better than GMRES algorithm. Tab.2 shows the CPU time spent by each of the three algorithms when m is 15 and 25 respectively in Example 4.3.

Tab.2 illustrates that IGMBACK is more efficient than GMBACK and that both algorithms are significantly better than GMRES algorithm.

Example 4.4　Consider the second order elliptic partial differential equation defined on the region [0,1]\times [0,1]:

-e^y\frac{{\mathrm{\partial }}^{2}u}{\mathrm{\partial }{x}^2}–e^x\frac{{\mathrm{\partial }}^{2}u}{\mathrm{\partial }{y}^2}+(x+y)\frac{\mathrm{\partial }u}{\mathrm{\partial }x}+(x-y)\frac{\mathrm{\partial }u}{\mathrm{\partial }y}+u=-(y^2{e}^y+x^2{e}^x)e^{xy}+(y^2+x^2+1)e^{xy},

the boundary condition is u(0,y)=1, u(1,y)=e^y, u(x,0)=1, u(x,1)=e^x. Differentiating the above second-order elliptic partial differential equations by the central difference method with lattice lengths A and h duals, respectively, yields a nonsymmetric difference linear systems Ax=b of coefficient matrix orders {39}^2 and {69}^2. Let x_0={11\dots 1}^T. The experimental results are shown in Tab.3 and Tab.4, respectively. The symbols in Tab.3 and Tab.4 are used as in Tab.1, where ratio represents the ratio of CPU time used by the algorithm \mathrm{I}\mathrm{G}\mathrm{M}\mathrm{B}\mathrm{A}\mathrm{C}\mathrm{K}(m) or GMRES(m) and GMBACK(m) for the same value of m.

As seen in Tab.3 and Tab.4: overall \mathrm{I}\mathrm{G}\mathrm{M}\mathrm{B}\mathrm{A}\mathrm{C}\mathrm{K}(m) and GMBACK(m) converge faster than CMRES(m). GMBACK(m) is more effective than GMRES(m), which in turn is more effective than \mathrm{I}\mathrm{G}\mathrm{M}\mathrm{B}\mathrm{A}\mathrm{C}\mathrm{K}(m). GMBACK(m) is more effective than GMRES(m), and \mathrm{I}\mathrm{G}\mathrm{M}\mathrm{B}\mathrm{A}\mathrm{C}\mathrm{K}(m) is more effective than GMBACK(m). We have additionally done a number of numerical experiments and found that \mathrm{I}\mathrm{G}\mathrm{M}\mathrm{B}\mathrm{A}\mathrm{C}\mathrm{K}(m) and GMBACK(m) often converge better than GMRES(m), and that \mathrm{I}\mathrm{G}\mathrm{M}\mathrm{B}\mathrm{A}\mathrm{C}\mathrm{K}(m) is generally more efficient than GMBACK(m) and GMRES(m).

5 Summary and expectation of the work

In today's increasingly popular Krylov subspace methods [15], this paper presents another method for solving large sparse linear systems of equations Krylov subspace method: The truncated form of GMBACK algorithm, that is, IGMBACK algorithm. The IGMBACK algorithm uses an incomplete orthogonalization process to generate a set of bases for {\mathcal{K}}_m(A,r_0), which overcomes the disadvantage of using Arnold process in CMBACK algorithm with large amount of computation and storage capacity. The detailed theoretical derivation and some theoretical analysis of the new algorithm are given. Numerical experiments have shown that IGMBACK is generally more effective than CMBACK and CMRES.

It should be noted that the GMRES algorithm has been greatly improved in recent years by the efforts of many authors, and has become a major iterative method for solving large nonsymmetric linear systems. Reference [18] lists the GMRES algorithm including some of the major deformations [11, 17, 19, 23] and theoretical developments [2, 16] in recent years. At present, the international exploration of GMRES algorithm is deepening, and GMRES algorithm itself has become a relatively independent research field. At the same time, this algorithm is also widely used in many fields such as science and engineering calculation. In this case, the advantages of IGMBACK and GMBACK over CMRES in solving many problems illustrate the value of IGMBACK and GMBACK.

The convergence problem of the GMBACK algorithm has not yet been solved. In the truncated version of the GMBACK algorithm, the cost is the loss of some important properties of the original method, such as orthogonality or A-orthogonality, minimization of the backward perturbation parametrization, etc., which are the main basis for the analysis of the convergence of the method, thus making the convergence of the IGMBACK algorithm more complicated. Therefore, one of our major research directions in the future is how to establish the convergence theory of these two algorithms and compare them with GMRES algorithm theoretically, so as to further determine theoretically what kind of equations the GMBACK series algorithm is better than the restarted GMRES for solving. In addition, in many cases and applications, direct use of the iterative method does not converge or converges very slowly. The current general solution is to combine iterative methods with preprocessing techniques. Therefore, another major research direction in the future is how to combine IGMBACK and GMBACK algorithm with preprocessing technology, so as to obtain a new algorithm with better convergence effect and higher accuracy.

References

[1] Arioli M, Duff I, Ruiz D. Stopping criteria for iterative solvers. SIAM J Matrix Anal Appl 1992; 13(1): 138–144

[2] Baker A H, Jessup E R, Kolev T V. A simple strategy for varying the restart parameter in GMRES(m). J Comput Appl Math 2009; 230(2): 751–761

[3] Ben-IsraelAGrevilleT N E. Generalized Inverses: Theory and Applications. New York: Wiley Interscience, 1974

[4] Brown P N. A theoretical comparison of the Arnoldi and GMRES algorithms. SIAM J Sci Stat Comput 1991; 12(1): 58–78

[5] Brown P N, Hindmarsh A C. Reduced storage matrix methods in stiff ODE systems. Appl Math Comput 1989; 31: 40–91

[6] Brown P N, Saad Y. Hybrid Krylov methods for nonlinear systems of equations. SIAM J Sci Stat Comput 1990; 11(3): 450–481

[7] Cao Z H. On a deflation method for the symmetric generalized eigenvalue problem. Linear Algebra Appl 1987; 92: 187–196

[8] Cao Z H. Total generalized minimum backward error algorithm for solving nonsymmetric linear systems. J Comput Math 1998; 16(6): 539–550

[9] GolubG HVan LoanC F. Matrix Computations, 2nd ed. Baltimore, MD: John Hopkins Univ Press, 1990

[10] Higham D J, Higham N J. Backward error and condition of structured linear systems. SIAM J Matrix Anal Appl 1992; 13(1): 162–175

[11] Huhtanen M, Permki A. Orthogonal polynomials of the R-linear generalized minimal residual method. J Approx Theory 2013; 167(3): 220–239

[12] Jia Z X. On IOM(q): the incomplete orthogonalization method for large unsymmetric linear systems. Numer Linear Algebra Appl 1996; 3(6): 491–512

[13] Kasenally E M. GMBACK: a generalized minimum backward error algorithm for nonsymmetric linear systems. SIAM J Sci Comput 1995; 16(3): 698–719

[14] Kasenally E M, Simoncini V. Analysis of a minimum perturbation algorithm for nonsymmetric linear systems. SIAM J Numer Anal 1997; 34(1): 48–66

[15] Li X A, Chen Y H, Zhang Y, Wang X P. Development of the Krylov subspace method for solving large sparse linear systems. science & Technology Review 2013; 11: 68–73

[16] LiesenJTichyP. The field of values bound on ideal GMRES, 2012, arXiv: 1211.5969v1

[17] Liu Y Q, Yin K X, WU E H. Fast GMRES-GPU algorithm for solving large scale sparse linear systems. Journal of Computer-Aided Design & Computer Graphics 2011; 23(4): 553–560

[18] Ma X F. An overview of recent developments and applications of the GMRES method. Pure Math 2013; 3: 181–187

[19] Najafi H S, Zareamonghaddam H. A new computational GMRES method. Appl Math Comput 2008; 199(2): 527–534

[20] Saad Y. Krylov subspace methods for solving large unsymmetric linear systems. Math Comp 1981; 37(155): 105–126

[21] Saad Y. Practical use of some Krylov subspace methods for solving indefinite and nonsymmetric linear systems. SIAM J Sci Stat Comput 1984; 5(1): 203–228

[22] Saad Y, Schultz M H. GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM J Sci Stat Comput 1986; 7(3): 856–869

[23] Sterck H D. Steepest descent preconditioning for nonlinear GMRES optimization. Numer Linear Algebra Appl 2013; 20(3): 453–471

[24] StewartG WSunJ-G. Matrix Perturbation Theory. New York: Academic Press, 1990

[25] Sun L, Wang X H, Guan Y. IMinpert: an incomplete minimum perturbation algorithm for large unsymmetric linear systems. Numer Math J Chinese Univ 2007; 16(4): 300–312

[26] Van HuffelSVandewalleJ. The Total Least Squares Problem: Computational Aspects and Analysis. Frontiers in Applied Mathematics, Vol 9. Philadelphia, PA: SIAM, 1991

RIGHTS & PERMISSIONS

Higher Education Press 2023