Re-Orthogonalized/Affine GMRES and Orthogonalized Maximal Projection Algorithm for Solving Linear Systems

Abstract

GMRES is one of the most powerful and popular methods to solve linear systems in the Krylov subspace; we examine it from two viewpoints: to maximize the decreasing length of the residual vector, and to maintain the orthogonality of the consecutive residual vector. A stabilization factor, $\eta$, to measure the deviation from the orthogonality of the residual vector is inserted into GMRES to preserve the orthogonality automatically. The re-orthogonalized GMRES (ROGMRES) method guarantees the absolute convergence; even the orthogonality is lost gradually in the GMRES iteration. When $\eta <1/2$, the residuals’ lengths of GMRES and GMRES(m) no longer decrease; hence, $\eta <1/2$ can be adopted as a stop** criterion to terminate the iterations. We prove $\eta =1$ for the ROGMRES method; it automatically keeps the orthogonality, and maintains the maximality for reducing the length of the residual vector. We improve GMRES by seeking the descent vector to minimize the residual in a larger space of the affine Krylov subspace. The resulting orthogonalized maximal projection algorithm (OMPA) is identified as having good performance. We further derive the iterative formulas by extending the GMRES method to the affine Krylov subspace; these equations are slightly different from the equations derived by Saad and Schultz (1986). The affine GMRES method is combined with the orthogonalization technique to generate a powerful affine GMRES (A-GMRES) method with high performance.

1. Introduction

In the paper, we revisit GMRES and GMRES(m) to solve

$$\mathbf{A}\mathbf{x}=\mathbf{b},\mathbf{x},\mathbf{b}\in {\mathbb{R}}^n,\mathbf{A}\in {\mathbb{R}}^{n\times n}.$$

(1)

In terms of an initial guess, ${\mathbf{x}}_0$, and the initial residual ${\mathbf{r}}_0=\mathbf{b}-\mathbf{A}{\mathbf{x}}_0$, Equation (1) is equivalent to

$$\mathbf{A}\mathbf{u}={\mathbf{r}}_0,$$

(2)

where

$$\mathbf{u}=\mathbf{x}-{\mathbf{x}}_0$$

(3)

is the descent vector or called a corrector.

Equation (2) can be used to search a proper descent vector, $\mathbf{u}$, which is inserted into Equation (3) to produce a solution, $\mathbf{x}$, closer to the exact solution. We suppose an m-dimensional Krylov subspace:

$${\mathcal{K}}_m(\mathbf{A},{\mathbf{r}}_0):=\mathrm{Span}\{{\mathbf{r}}_0,\mathbf{A}{\mathbf{r}}_0,\dots ,{\mathbf{A}}^{m-1}{\mathbf{r}}_0\}.$$

(4)

Krylov subspaces are named after A.N. Krylov in a paper [1], who described a method for solving a secular equation to determine the frequencies of small vibrations of material systems.

Many well-developed Krylov subspace methods for solving Equation (1) have been discussed in the text books [2,3,4,5]. Godunov and Gordienko [6] used an optimal representation of vectors in the Krylov space with the help of a variational problem; the extremum of the variational problem is the solution of the Kalman matrix equation, in which the 2-norm minimal solution is used as a characteristic of the Krylov space.

One of the most important classes of numerical methods for iteratively solving Equation (1) is the Krylov subspace methods [7,8], and the preconditioned Krylov subspace methods [9]. Recently, Bouyghf et al. [10] presented a comprehensive framework for unifying the Krylov subspace methods for Equation (1).

The GMRES method in [11] used the Petrov–Galerkin method to search $\mathbf{u}\in {\mathcal{K}}_m$ via a perpendicular property:

$${\mathbf{r}}_0-\mathbf{A}\mathbf{u}\perp {\mathcal{L}}_m=\mathbf{A}{\mathcal{K}}_m.$$

(5)

The descent vector $\mathbf{u}\in {\mathcal{K}}_m$ is achieved by minimizing the length of the residual vector [2]:

$$min\parallel \mathbf{r}\parallel ,$$

(6)

where

$$\mathbf{r}=\mathbf{b}-\mathbf{A}\mathbf{x}={\mathbf{r}}_0-\mathbf{A}\mathbf{u},$$

(7)

and $\parallel \mathbf{r}\parallel$ is the Euclidean norm of the residual vector, $\mathbf{r}$.

In many Krylov subspace methods, the central idea is the projection. The method of GMRES is a special case of the Krylov subspace methods with an oblique projection. Let $\mathbf{V}$ and $\mathbf{W}=\mathbf{A}\mathbf{V}$ be the matrix representations of ${\mathcal{K}}_m$ and ${\mathcal{L}}_m$, respectively. We project the residual Equation (2) from an n-dimensional space to an m-dimensional subspace:

$${\mathbf{W}}^{\mathrm{T}}\mathbf{A}\mathbf{u}={\mathbf{W}}^{\mathrm{T}}{\mathbf{r}}_0.$$

(8)

We seek $\mathbf{u}$ in the subspace ${\mathcal{K}}_m$ by

$$\mathbf{u}=\mathbf{V}\mathbf{y};$$

(9)

hence, Equation (8) becomes

$${\mathbf{W}}^{\mathrm{T}}\mathbf{W}\mathbf{y}={\mathbf{W}}^{\mathrm{T}}{\mathbf{r}}_0.$$

(10)

Solving $\mathbf{y}$ from Equation (10) and inserting it into Equation (9), we can obtain $\mathbf{u}$, and, hence, the next $\mathbf{x}={\mathbf{x}}_0+\mathbf{u}$; $\mathbf{x}$ is a correction of ${\mathbf{x}}_0$, and we hope $\parallel \mathbf{r}\parallel < \parallel {\mathbf{r}}_0\parallel$.

Meanwhile, Equation (7) becomes

$$\mathbf{r}={\mathbf{r}}_0-\mathbf{W}\mathbf{y},$$

(11)

where $\mathbf{W}\mathbf{y}=\mathbf{A}\mathbf{u}$. By means of Equations (5), (7), and (11), we enforce the perpendicular property by

$${\mathbf{r}}_0-\mathbf{W}\mathbf{y}\perp {\mathcal{L}}_m;$$

(12)

it is equivalent to the projection of ${\mathbf{r}}_0$ on ${\mathcal{L}}_m$. $\mathbf{W}\mathbf{y}=\mathbf{A}\mathbf{u}$ is the projection of ${\mathbf{r}}_0$ on ${\mathcal{L}}_m$, and $\mathbf{r}$ is perpendicular to ${\mathcal{L}}_m$, which ensures that the length of the residual vector is decreasing. For seeking a fast convergent iterative algorithm, we must keep the orthogonality, and simultaneously maximize the projection quantity $\parallel \mathbf{W}\mathbf{y}\parallel =\parallel \mathbf{A}\mathbf{u}\parallel$. Later on, we will modify the GMRES method from these two aspects.

Liu [12] derived an effective iterative algorithm based on the maximal projection method to solve Equation (1) in the column space of $\mathbf{A}$. Recently, Liu et al. [13] applied the maximal projection method and the minimal residual norm method to solve the least-squares problems.

There are two reasons to cause the slowdown and even the stagnation of GMRES: one is the loss of orthogonality, and other is that its descent vector, $\mathbf{u}=\mathbf{V}\mathbf{y}$, is not the best one. There are different methods to accelerate and improve GMRES [14,15,16,17]. The restart technique is often applied to the Krylov subspace methods, which, however, generally slows their convergence. Another way to accelerate the convergence is updating the initial value [18,19]. The preconditioning technology is often used to speed up GMRES and other Krylov subspace methods [20,21,22,23,24]. The multipreconditioned GMRES [25] is also proposed by using two or more preconditioners simultaneously. The polynomial preconditioned GMRES and GMRES-DR were also developed in [26]. As an augment of the Krylov subspace, some authors sought better descent vectors in a larger space to accelerate the convergence of GMRES [27,28,29]. Recently, Carson et al. [30] presented mathematically important and practically relevant phenomena that uncover the behavior of GMRES through a discussion of computed examples; they considered the conjugate gradient (CG) and GMRES methods crucially important for further development and practical applications, as well as other Krylov subspace methods.

Contribution and Novelty:

The major contribution and novelty of the present paper are that we propose a quite simple strategy to overcome the slowdown and stagnation behavior of GMRES by inserting a stabilization factor, $\eta ={\mathbf{r}}^{\mathrm{T}}\mathbf{A}\mathbf{u}/{\parallel \mathbf{A}\mathbf{u}\parallel }^2$, into the iterative equation. In doing so, the orthogonality property of the residual vector and the stepwise decreasing property of the residual length are automatically preserved. We extend the space of GMRES to an affine Krylov subspace for seeking a better descent vector. When the maximal projection method and the affine technique are combined into the orthogonalization technique, two very powerful iterative algorithms for solving linear systems with high performance are developed.

Highlight:

• Examine the iterative algorithm for linear systems from the viewpoint of maximizing the decreasing quantity of the residual and maintaining the orthogonality of the consecutive residual vector.
• A stabilization factor to measure the deviation from the orthogonality is inserted into GMRES to preserve the orthogonality automatically.
• The re-orthogonalized GMRES guarantees the preservation of orthogonality, even if the orthogonality is gradually lost in the iterations by means of GMRES.
• Improve GMRES by seeking the descent vector to minimize the length of the residual vector in a larger space of the affine Krylov subspace.
• The new algorithms are all absolute convergence.

Outline:

The Arnoldi process and the conventional GMRES and GMRES(m) are described in Section 2. In Section 3, we propose a new algorithm modified from GMRES by inserting a stabilization factor, $\eta$, which is named a re-orthogonalized GMRES (ROGMRES). We discuss ROGMRES from two aspects: preserving the orthogonality and maximally decreasing the length of the residual vector; we also prove that ROGMRES can automatically maintain $\eta =1$, and preserve the good property of orthogonality and the maximality of reducing the residual. In Section 4, a better descent vector is sought in a larger affine Krylov subspace by using the maximal projection method; upon combining it with the orthogonalization technique, we propose a new algorithm, the orthogonalized maximal projection algorithm (OMPA). The iterative equations in GMRES are extended to that in the affine Krylov subspace by using the minimum residual method in Section 5, which results in an affine GMRES (A-GMRES) method. The numerical tests of some examples are given in Section 6. Section 7 concludes the achievements, novelties, and contributions of this paper.

2. GMRES and GMRES(m)

The Arnoldi process [2] is often used to orthonormalize the Krylov vectors ${\mathbf{A}}^j{\mathbf{r}}_0,j=0,\dots ,m-1$ in Equation (4), such that the resultant vectors ${\mathbf{v}}_i,i=1,\dots ,m$ satisfy ${\mathbf{v}}_i·{\mathbf{v}}_j={\delta }_{ij},i,j=1,\dots ,m$, where ${\delta }_{ij}$ is the Kronecker delta symbol. The full orthonormalization procedure, known as the Algorithm 1, is set up as follows.   

Algorithm 1: Arnoldi process

1:  Select m and give an initial ${\mathbf{r}}_0$ 2:  ${\mathbf{v}}_1=\frac{{\mathbf{r}}_0}{\parallel {\mathbf{r}}_0\parallel }$ 3:  Do $j=1:m$ 4:  ${\mathbf{w}}_j=\mathbf{A}{\mathbf{v}}_j$ 5:     Do $i=1:j$ 6:     $h_{ij}={\mathbf{v}}_i·{\mathbf{w}}_j$ 7:     ${\mathbf{w}}_j={\mathbf{w}}_j-h_{ij}{\mathbf{v}}_i$ 8:     Enddo of i 9:  $h_{j+1,j}=\parallel {\mathbf{w}}_j\parallel$ 10: ${\mathbf{v}}_{j+1}=\frac{{\mathbf{w}}_j}{h_{j+1,j}}$ 11: Enddo of j

A dot between two vectors represents the inner product, like as $\mathbf{a}·\mathbf{b}={\mathbf{a}}^{\mathrm{T}}\mathbf{b}$. $\mathbf{V}$ denotes the Arnoldi matrix whose jth column is ${\mathbf{v}}_j$:

$$\mathbf{V}:=[{\mathbf{v}}_1,\dots ,{\mathbf{v}}_m]\in {\mathbb{R}}^{n\times m}.$$

(13)

After m steps of the Arnoldi process, we can construct an upper Hessenberg matrix, ${\mathbf{H}}_m$, as follows:

$${\mathbf{H}}_m=\left[\begin{array}{cccc}{h}_{11}& h_{12}& \dots & h_{1m}\\ h_{21}& h_{22}& \dots & h_{2m}\\ & h_{32}& \dots & h_{3m}\\ & \ddots & & ⋮\\ & & h_{m,m-1}& h_{mm}\end{array}\right].$$

(14)

Utilizing ${\mathbf{H}}_m$, the Arnoldi process can be expressed as a matrix product form:

$$\mathbf{A}\mathbf{V}=\mathbf{V}{\mathbf{H}}_m+h_{m+1,m}{\mathbf{v}}_{m+1}{\mathbf{e}}_m^{\mathrm{T}},$$

(15)

where ${\mathbf{H}}_m\in {\mathbb{R}}^{m\times m}$, and ${\mathbf{e}}_m={[0,\dots ,0,1]}^{\mathrm{T}}\in {\mathbb{R}}^m$. Now, the augmented Hessenberg upper matrix, ${\overline{\mathbf{H}}}_m\in {\mathbb{R}}^{(m+1)\times m}$, can be formulated as

$${\overline{\mathbf{H}}}_m=\left[\begin{array}{c}{\mathbf{H}}_m\\ h_{m+1,m}{\mathbf{e}}_m^{\mathrm{T}}\end{array}\right]=\left[\begin{array}{cccc}{h}_{11}& h_{12}& \dots & h_{1m}\\ h_{21}& h_{22}& \dots & h_{2m}\\ & h_{32}& \dots & h_{3m}\\ & \ddots & & ⋮\\ & & h_{m,m-1}& h_{mm}\\ 0& \dots & 0& h_{m+1,m}\end{array}\right].$$

(16)

According to [2], the Algorithm 2 method is given as follows.

Algorithm 2: GMRES

1: Give $m_0$, ${\mathbf{x}}_0$, and $\epsilon$ 2: Do $m=m_0,\dots$ ($k=m-m_0$), until $\parallel {\mathbf{r}}_k\parallel <\epsilon$ 3: ${\mathbf{r}}_k=\mathbf{b}-\mathbf{A}{\mathbf{x}}_k$ 4: ${\mathbf{v}}_1^k=\frac{{\mathbf{r}}_k}{\parallel {\mathbf{r}}_k\parallel }$ 5: ${\mathbf{V}}_k=[{\mathbf{v}}_1^k,\dots ,{\mathbf{v}}_m^k]$  (by Arnoldi process) 6: Solve $\left({\overline{\mathbf{H}}}_k^{\mathrm{T}}{\overline{\mathbf{H}}}_k\right){\mathbf{y}}_k=\parallel {\mathbf{r}}_k\parallel {\overline{\mathbf{H}}}_k^{\mathrm{T}}{\mathbf{e}}_1$ 7: ${\mathbf{x}}_{k+1}={\mathbf{x}}_k+{\mathbf{V}}_k{\mathbf{y}}_k$

In the above, ${\mathbf{x}}_k$ is the kth step value of $\mathbf{x}$; ${\mathbf{r}}_k$ is the kth residual vector; ${\mathbf{y}}_k$ is the vector of m expansion coefficients used in ${\mathbf{V}}_k{\mathbf{y}}_k$; ${\mathbf{e}}_1$ is the first column of ${\mathbf{I}}_{m+1}$; ${\mathbf{V}}_k$ denotes the kth step $\mathbf{V}$ in Equation (13); ${\overline{\mathbf{H}}}_k\in {\mathbb{R}}^{(m+1)\times m}$ is the kth step augmented Hessenberg upper matrix in Equation (16); and ${\mathbf{V}}_k{\mathbf{y}}_k$ is the kth step descent vector. Notice that ${\mathbf{V}}_k{\mathbf{y}}_k\in {\mathcal{K}}_m(\mathbf{A},{\mathbf{r}}_k)$.

In addition to Algorithm 3, there is a restarted version (GMRES(m)) [11,31] described as follows.

Algorithm 3: GMRES(m)

1: Give $m_0$, $m_1$, ${\mathbf{x}}_0$, and $\epsilon$ 2: Do $j=1,\dots$ 3: Do $m=m_0,\dots ,m_1$ ($k=m-m_0$) 4: ${\mathbf{r}}_k=\mathbf{b}-\mathbf{A}{\mathbf{x}}_k$ 5: ${\mathbf{v}}_1^k=\frac{{\mathbf{r}}_k}{\parallel {\mathbf{r}}_k\parallel }$ 6: ${\mathbf{V}}_k=[{\mathbf{v}}_1^k,\dots ,{\mathbf{v}}_m^k]$  (by Arnoldi process) 7: Solve $\left({\overline{\mathbf{H}}}_k^{\mathrm{T}}{\overline{\mathbf{H}}}_k\right){\mathbf{y}}_k=\parallel {\mathbf{r}}_k\parallel {\overline{\mathbf{H}}}_k^{\mathrm{T}}{\mathbf{e}}_1$ 8: ${\mathbf{x}}_{k+1}={\mathbf{x}}_k+{\mathbf{V}}_k{\mathbf{y}}_k$ 9: If $\parallel {\mathbf{r}}_{k+1}\parallel <\epsilon$, stop 10: Otherwise, ${\mathbf{x}}_0={\mathbf{x}}_{k+1}$; go to 2

Instead of $m_0=1$ used in the original GMRES and GMRES(m), a suitable value of $m_0\ne 1$ can speed up the convergence; hence, $m=m_1-m_0+1$ is the frequency for restart. Rather than the original GMRES and GMRES(m) [2,11,31], we take ${\mathbf{x}}_k$, not ${\mathbf{x}}_0$, in (7) of Algorithm GMRES and in (8) of Algorithm GMRES(m). Because ${\mathbf{x}}_0$ is not updated in the original GMRES and GMRES(m), they converge slowly. In Algorithm GMRES(m), if we take $m_1=m_0$, then we return to the usual iterative algorithm for GMRES, with a fixed-dimension, $m=m_0$, of the Krylov subspace.

The restart remedies the long-term recurrence and accumulations of the round-off errors of GMRES. In order to improve the restart process, several improvement techniques have been proposed, as mentioned above. Some techniques based on adaptively determining the restart frequency can be found in [32,33,34,35]. For a quicker convergence, we can update the current value of $\mathbf{x}$ at each iterative step. Because ${\mathbf{x}}_k$ is new information for determining the next value of ${\mathbf{x}}_{k+1}$, the result of ${\mathbf{x}}_k$ cannot be wasted for saving computational cost.

3. Maximally Decreasing Residual and Automatically Preserving Orthogonality

In this section, we develop a new iterative algorithm based on GMRES to solve Equation (1). The iterative form is

$${\mathbf{x}}_{k+1}={\mathbf{x}}_k+{\mathbf{u}}_k,$$

(17)

where ${\mathbf{u}}_k$ is the kth step descent vector, which is determined by the iterative algorithm. For instance, ${\mathbf{u}}_k$ in GMRES is ${\mathbf{u}}_k={\mathbf{V}}_k{\mathbf{y}}_k$.

Lemma 1.

A better iterative scheme than that in Equation (17) is

$${\mathbf{x}}_{k+1}={\mathbf{x}}_k+\frac{{\mathbf{r}}_k^T\mathbf{A}{\mathbf{u}}_k}{\parallel \mathbf{A}{\mathbf{u}}_k{\parallel }^2}{\mathbf{u}}_k.$$

(18)

Proof. 

Let

$${\mathbf{x}}_{k+1}={\mathbf{x}}_k+\eta {\mathbf{u}}_k;$$

(19)

it follows that

$${\mathbf{r}}_{k+1}={\mathbf{r}}_k-\eta \mathbf{A}{\mathbf{u}}_k.$$

(20)

Taking the squared norms of both sides yields

$$\parallel {\mathbf{r}}_{k+1}{\parallel }^2=\parallel {\mathbf{r}}_k{\parallel }^2-(2\eta {\mathbf{r}}_k^{\mathrm{T}}\mathbf{A}{\mathbf{u}}_k-{\eta }^2\parallel \mathbf{A}{\mathbf{u}}_k{\parallel }^2)=:\parallel {\mathbf{r}}_k{\parallel }^2-h.$$

(21)

To maximally reduce the residual we consider

$$\underset{\eta }{max}\left\{h=2\eta {\mathbf{r}}_k^{\mathrm{T}}\mathbf{A}{\mathbf{u}}_k-{\eta }^2{\parallel \mathbf{A}{\mathbf{u}}_k\parallel }^2\right\},$$

(22)

which leads to

$${\eta }_k=\frac{{\mathbf{r}}_k^{\mathrm{T}}\mathbf{A}{\mathbf{u}}_k}{\parallel \mathbf{A}{\mathbf{u}}_k{\parallel }^2}.$$

(23)

Inserting ${\eta }_k$ into Equation (19), the proof of Equation (18) is finished.    □

Lemma 2.

In Equation (18) the following orthogonal property holds:

$${\mathbf{r}}_{k+1}·\left(\mathbf{A}{\mathbf{u}}_k\right)=0.$$

(24)

Proof. 

Apply $\mathbf{A}$ to Equation (18):

$$\mathbf{A}{\mathbf{x}}_{k+1}=\mathbf{A}{\mathbf{x}}_k+\frac{{\mathbf{r}}_k^{\mathrm{T}}\mathbf{A}{\mathbf{u}}_k}{\parallel \mathbf{A}{\mathbf{u}}_k{\parallel }^2}\mathbf{A}{\mathbf{u}}_k.$$

(25)

Using $\mathbf{A}{\mathbf{x}}_{k+1}=\mathbf{b}-{\mathbf{r}}_{k+1}$ and $\mathbf{A}{\mathbf{x}}_k=\mathbf{b}-{\mathbf{r}}_k$, Equation (25) changes to

$${\mathbf{r}}_{k+1}={\mathbf{r}}_k-\frac{{\mathbf{r}}_k^{\mathrm{T}}\mathbf{A}{\mathbf{u}}_k}{\parallel \mathbf{A}{\mathbf{u}}_k{\parallel }^2}\mathbf{A}{\mathbf{u}}_k.$$

(26)

Taking the inner product with $\mathbf{A}{\mathbf{u}}_k$ yields

$${\mathbf{r}}_{k+1}·\left(\mathbf{A}{\mathbf{u}}_k\right)={\mathbf{r}}_k^{\mathrm{T}}\mathbf{A}{\mathbf{u}}_k-\frac{{\mathbf{r}}_k^{\mathrm{T}}\mathbf{A}{\mathbf{u}}_k}{\parallel \mathbf{A}{\mathbf{u}}_k{\parallel }^2}{\parallel \mathbf{A}{\mathbf{u}}_k\parallel }^2=0.$$

(27)

Equation (24) is proven.    □

It follows from Equation (5) that $\eta =1$ is implied by GMRES. However, during the GMRES iteration the numerical values of $\eta$ may deviate from 1 to a great extent, even being zero or a negative value. Therefore, we propose a re-orthogonalized version of GMRES, namely the ROGMRES method (Algorithm 4), as follows.

Algorithm 4: ROGMRES

1: Give $m_0$, ${\mathbf{x}}_0$, and $\epsilon$ 2: Do $m=m_0,\dots$ ($k=m-m_0$), until $\parallel {\mathbf{r}}_k\parallel <\epsilon$ 3: ${\mathbf{r}}_k=\mathbf{b}-\mathbf{A}{\mathbf{x}}_k$ 4: ${\mathbf{v}}_1^k=\frac{{\mathbf{r}}_k}{\parallel {\mathbf{r}}_k\parallel }$ 5: ${\mathbf{V}}_k=[{\mathbf{v}}_1^k,\dots ,{\mathbf{v}}_m^k]$  (by Arnoldi process) 6: Solve $\left({\overline{\mathbf{H}}}_k^{\mathrm{T}}{\overline{\mathbf{H}}}_k\right){\mathbf{y}}_k=\parallel {\mathbf{r}}_k\parallel {\overline{\mathbf{H}}}_k^{\mathrm{T}}{\mathbf{e}}_1$ 7: ${\mathbf{u}}_k={\mathbf{V}}_k{\mathbf{y}}_k$ 8: ${\mathbf{x}}_{k+1}={\mathbf{x}}_k+\frac{{\mathbf{r}}_k^{\mathrm{T}}\mathbf{A}{\mathbf{u}}_k}{\parallel \mathbf{A}{\mathbf{u}}_k{\parallel }^2}{\mathbf{u}}_k$

Correspondingly, the ROGMRES(m) Algorithm 5 is given as follows.

Algorithm 5: ROGMRES(m)

1: Give $m_0$, $m_1$, ${\mathbf{x}}_0$, and $\epsilon$ 2: Do $j=1,\dots$ 3: Do $m=m_0,\dots ,m_1$ ($k=m-m_0$) 4: ${\mathbf{r}}_k=\mathbf{b}-\mathbf{A}{\mathbf{x}}_k$ 5: ${\mathbf{v}}_1^k=\frac{{\mathbf{r}}_k}{\parallel {\mathbf{r}}_k\parallel }$ 6: ${\mathbf{V}}_k=[{\mathbf{v}}_1^k,\dots ,{\mathbf{v}}_m^k]$  (by Arnoldi process) 7: Solve $\left({\overline{\mathbf{H}}}_k^{\mathrm{T}}{\overline{\mathbf{H}}}_k\right){\mathbf{y}}_k=\parallel {\mathbf{r}}_k\parallel {\overline{\mathbf{H}}}_k^{\mathrm{T}}{\mathbf{e}}_1$ 8: ${\mathbf{u}}_k={\mathbf{V}}_k{\mathbf{y}}_k$ 9: ${\mathbf{x}}_{k+1}={\mathbf{x}}_k+\frac{{\mathbf{r}}_k^{\mathrm{T}}\mathbf{A}{\mathbf{u}}_k}{\parallel \mathbf{A}{\mathbf{u}}_k{\parallel }^2}{\mathbf{u}}_k$ 10: If $\parallel {\mathbf{r}}_{k+1}\parallel <\epsilon$, stop 11: Otherwise, ${\mathbf{x}}_0={\mathbf{x}}_{k+1}$; go to 2

Upon comparing with GMRES(m) in Section 2, the computational cost of ROGMRES(m) is slightly increased by computing an extra term, ${\mathbf{r}}_k^{\mathrm{T}}\mathbf{A}{\mathbf{u}}_k/{\parallel \mathbf{A}{\mathbf{u}}_k\parallel }^2$, at each iteration; it needs one matrix-vector product and two inner products of vectors.

Theorem 1.

For GMRES, if ${\eta }_k\ne 1$ is happened at the kth step, the orthogonality in Equation (5) cannot be continued after the kth step, i.e.,

$${\mathbf{r}}_{k+1}·\left(\mathbf{A}{\mathbf{u}}_k\right)\ne 0;$$

(28)

moreover, if ${\eta }_k<1/2$, the residual does not decrease, i.e.,

$$\parallel {\mathbf{r}}_{k+1}\parallel >\parallel {\mathbf{r}}_k\parallel .$$

(29)

For ROGMRES in Equation (18), the residual is absolutely decreased:

$$\parallel {\mathbf{r}}_{k+1}\parallel <\parallel {\mathbf{r}}_k\parallel ,$$

(30)

and the orthogonality ${\mathbf{r}}_{k+1}·\left(\mathbf{A}{\mathbf{u}}_k\right)=0$ holds.

Proof. 

Applying $\mathbf{A}$ to ${\mathbf{x}}_{k+1}={\mathbf{x}}_k+{\mathbf{V}}_k{\mathbf{y}}_k$ in Algorithm GMRES yields

$${\mathbf{r}}_{k+1}={\mathbf{r}}_k-\mathbf{A}{\mathbf{u}}_k,$$

(31)

where ${\mathbf{u}}_k={\mathbf{V}}_k{\mathbf{y}}_k$.

It follows from Equation (23) that

$$({\mathbf{r}}_k-{\eta }_k\mathbf{A}{\mathbf{u}}_k)·\left(\mathbf{A}{\mathbf{u}}_k\right)=0,$$

(32)

which can be rearranged to

$$({\mathbf{r}}_k-\mathbf{A}{\mathbf{u}}_k)·\left(\mathbf{A}{\mathbf{u}}_k\right)=({\eta }_k-1){\parallel \mathbf{A}{\mathbf{u}}_k\parallel }^2,$$

(33)

and then by Equation (31),

$${\mathbf{r}}_{k+1}·\left(\mathbf{A}{\mathbf{u}}_k\right)=({\eta }_k-1){\parallel \mathbf{A}{\mathbf{u}}_k\parallel }^2.$$

(34)

If ${\eta }_k\ne 1$, we can prove Equation (28). Taking the squared norm of Equation (31) and using Equation (23) generates

$$\parallel {\mathbf{r}}_{k+1}{\parallel }^2=\parallel {\mathbf{r}}_k{\parallel }^2-2{\mathbf{r}}_k·\left(\mathbf{A}{\mathbf{u}}_k\right)+\parallel \mathbf{A}{\mathbf{u}}_k{\parallel }^2=\parallel {\mathbf{r}}_k{\parallel }^2+(1-2{\eta }_k){\parallel \mathbf{A}{\mathbf{u}}_k\parallel }^2.$$

(35)

If ${\eta }_k<1/2$, $(1-2{\eta }_k){\parallel \mathbf{A}{\mathbf{u}}_k\parallel }^2>0$ and Equation (29) is proven.

For ROGMRES, taking the squared norm of Equation (26), we have

$$\parallel {\mathbf{r}}_{k+1}{\parallel }^2={\parallel {\mathbf{r}}_k\parallel }^2-\frac{{\left({\mathbf{r}}_k^{\mathrm{T}}\mathbf{A}{\mathbf{u}}_k\right)}^2}{\parallel \mathbf{A}{\mathbf{u}}_k{\parallel }^2}.$$

(36)

Because ${\left({\mathbf{r}}_k^{\mathrm{T}}\mathbf{A}{\mathbf{u}}_k\right)}^2/{\parallel \mathbf{A}{\mathbf{u}}_k\parallel }^2>0$, Equation (30) is proven. The orthogonality ${\mathbf{r}}_{k+1}·\left(\mathbf{A}{\mathbf{u}}_k\right)=0$ was proven in Equation (24) of Lemma 2.    □

Corollary 1.

For GMRES, the orthogonality of the consecutive residual vector is equivalent to the maximality of h in Equation (22), where

$${\eta }_k=1\iff {\mathbf{r}}_k^T\mathbf{A}{\mathbf{u}}_k={\parallel \mathbf{A}{\mathbf{u}}_k\parallel }^2.$$

(37)

Proof. 

By means of Equation (34), to satisfy the orthogonality condition of the residual vector we require ${\eta }_k=1$; it implies Equation (37) by the definition in Equation (23). At the same time, Equation (36) changes to

$$\parallel {\mathbf{r}}_{k+1}{\parallel }^2=\parallel {\mathbf{r}}_k{\parallel }^2-{\parallel \mathbf{A}{\mathbf{u}}_k\parallel }^2.$$

(38)

Inserting ${\eta }_k=1$ into Equation (22), we have

$$maxh=2{\mathbf{r}}_k^{\mathrm{T}}\mathbf{A}{\mathbf{u}}_k-\parallel \mathbf{A}{\mathbf{u}}_k{\parallel }^2={\parallel \mathbf{A}{\mathbf{u}}_k\parallel }^2,$$

(39)

of which ${\mathbf{r}}_k^{\mathrm{T}}\mathbf{A}{\mathbf{u}}_k={\parallel \mathbf{A}{\mathbf{u}}_k\parallel }^2$ in Equation (37) was used. For GMRES, obtaining the maximal value of h and preserving the orthogonality are the same.    □

Remark 1.

To maintain ${\eta }_k=1$ in GMRES is a key issue for preserving the orthogonality and for the maximality of reducing the residual. However, for the traditional GMRES, it is not always true that it can maintain ${\eta }_k=1$ during the iteration process.

Theorem 2.

Let

$${\mathbf{w}}_k=\frac{{\mathbf{r}}_k^T\mathbf{A}{\mathbf{u}}_k}{\parallel \mathbf{A}{\mathbf{u}}_k{\parallel }^2}{\mathbf{u}}_k$$

(40)

be the descent vector of ROGMRES in

$${\mathbf{x}}_{k+1}={\mathbf{x}}_k+{\mathbf{w}}_k.$$

(41)

The following identity holds:

$${\eta }_k=\frac{{\mathbf{r}}_k^T\mathbf{A}{\mathbf{w}}_k}{\parallel \mathbf{A}{\mathbf{w}}_k{\parallel }^2}=1;$$

(42)

such that ROGMRES can automatically maintain the orthogonality of the residual vector, and achieve the maximality of the decreasing length of the residual vector.

Proof. 

Let

$$\gamma =\frac{{\mathbf{r}}_k^{\mathrm{T}}\mathbf{A}{\mathbf{u}}_k}{\parallel \mathbf{A}{\mathbf{u}}_k{\parallel }^2};$$

(43)

Equation (40) is written as

$${\mathbf{w}}_k=\gamma {\mathbf{u}}_k.$$

(44)

Inserting it into Equation (42) produces

$${\eta }_k=\frac{1}{\gamma }\frac{{\mathbf{r}}_k^{\mathrm{T}}\mathbf{A}{\mathbf{u}}_k}{\parallel \mathbf{A}{\mathbf{u}}_k{\parallel }^2};$$

(45)

by means of Equation (43), it follows that

$${\eta }_k=\frac{\parallel \mathbf{A}{\mathbf{u}}_k{\parallel }^2}{{\mathbf{r}}_k^{\mathrm{T}}\mathbf{A}{\mathbf{u}}_k}\frac{{\mathbf{r}}_k^{\mathrm{T}}\mathbf{A}{\mathbf{u}}_k}{\parallel \mathbf{A}{\mathbf{u}}_k{\parallel }^2}=1.$$

(46)

Thus, Equation (42) is proven. The proof to satisfy the orthogonality condition of the consecutive residual vector and the maximality of h is the same as that given in Corollary 1.    □

Remark 2.

The properties in ROGMRES are crucial so that it can automatically maintain ${\eta }_k=1$, for preserving its good properties of the orthogonality, and the maximality for reducing the residual. The method in Equation (18) can be viewed as an automatically orthogonality preserving method. Numerical experiments will verify these points.

Theorem 3.

For ROGMRES in Equation (18), the iterative point, $x_k$, is located on an invariant manifold:

$$\parallel {\mathbf{r}}_k{\parallel }^2=\frac{\parallel {\mathbf{r}}_0{\parallel }^2}{Q_k},$$

(47)

where $Q_0=1$; $Q_k>0$ with $Q_{k+1}>Q_k$ is an increasing sequence.

Proof. 

We begin with the following time-invariant manifold:

$${\parallel \mathbf{r}\parallel }^2=\frac{\parallel {\mathbf{r}}_0{\parallel }^2}{Q\left(t\right)},$$

(48)

where $Q\left(0\right)=1$, and $Q\left(t\right)$ is an increasing function of t.

We suppose that the evolution of $\mathbf{x}$ in time is driven by a descent vector, $\mathbf{u}$:

$$\dot{\mathbf{x}}=\gamma \mathbf{u}.$$

(49)

Equation (48) is a time-invariant manifold, such that its time differential is zero, i.e.,

$$\frac{1}{2Q\left(t\right)}\dot{Q}\left(t\right){\parallel \mathbf{r}\parallel }^2-\mathbf{r}·\left(\mathbf{A}\dot{\mathbf{x}}\right)=0,$$

(50)

where $\dot{\mathbf{r}}=-\mathbf{A}\dot{\mathbf{x}}$ was used.

Inserting Equation (49) into Equation (50) yields

$$\gamma =\frac{\dot{Q}\left(t\right)}{2Q\left(t\right)}\frac{{\parallel \mathbf{r}\parallel }^2}{{\mathbf{r}}^{\mathrm{T}}\mathbf{A}\mathbf{u}}.$$

(51)

Inserting $\gamma$ into Equation (49) and using $\mathbf{v}=\mathbf{A}\mathbf{u}$, we can derive

$$\dot{\mathbf{x}}=q\left(t\right)\frac{{\parallel \mathbf{r}\parallel }^2}{{\mathbf{r}}^{\mathrm{T}}\mathbf{v}}\mathbf{u},$$

(52)

where

$$q\left(t\right):=\frac{\dot{Q}\left(t\right)}{2Q\left(t\right)}>0.$$

(53)

By applying the forward Euler scheme to Equation (52), we have

$$\mathbf{x}(t+\Delta t)=\mathbf{x}\left(t\right)+\beta \frac{{\parallel \mathbf{r}\parallel }^2}{{\mathbf{r}}^{\mathrm{T}}\mathbf{v}}\mathbf{u},$$

(54)

where

$$\beta =q\left(t\right)\Delta t>0.$$

(55)

An iterative form of Equation (54) is

$${\mathbf{x}}_{k+1}={\mathbf{x}}_k+{\beta }_k\frac{\parallel {\mathbf{r}}_k{\parallel }^2}{{\mathbf{r}}_k^{\mathrm{T}}{\mathbf{v}}_k}{\mathbf{u}}_k.$$

(56)

It generates the residual form:

$${\mathbf{r}}_{k+1}={\mathbf{r}}_k-{\beta }_k\frac{\parallel {\mathbf{r}}_k{\parallel }^2}{{\mathbf{r}}_k^{\mathrm{T}}{\mathbf{v}}_k}{\mathbf{v}}_k.$$

(57)

Taking the squared norms of both sides yields

$$\parallel {\mathbf{r}}_{k+1}{\parallel }^2=\parallel {\mathbf{r}}_k{\parallel }^2-2{\beta }_k\parallel {\mathbf{r}}_k{\parallel }^2+{\beta }_k^2\frac{\parallel {\mathbf{r}}_k{\parallel }^4}{{\left({\mathbf{r}}_k^{\mathrm{T}}{\mathbf{v}}_k\right)}^2}{\parallel {\mathbf{v}}_k\parallel }^2.$$

(58)

Dividing both sides by $\parallel {\mathbf{r}}_k{\parallel }^2$ renders

$$a_k{\beta }_k^2-2{\beta }_k+1-\frac{\parallel {\mathbf{r}}_{k+1}{\parallel }^2}{\parallel {\mathbf{r}}_k{\parallel }^2}=0,$$

(59)

where

$$a_k:=\frac{\parallel {\mathbf{r}}_k{\parallel }^2{\parallel {\mathbf{v}}_k\parallel }^2}{{\left({\mathbf{r}}_k^{\mathrm{T}}{\mathbf{v}}_k\right)}^2}\ge 1.$$

(60)

By means of Equation (47), Equation (59) changes to

$$a_k{\beta }_k^2-2{\beta }_k+1-\frac{Q_k}{Q_{k+1}}=0.$$

(61)

If we take

$${\beta }_k=\frac{1}{a_k}=\frac{{\left({\mathbf{r}}_k^{\mathrm{T}}{\mathbf{v}}_k\right)}^2}{\parallel {\mathbf{r}}_k{\parallel }^2{\parallel {\mathbf{v}}_k\parallel }^2},\frac{Q_k}{Q_{k+1}}=1-\frac{1}{a_k},$$

(62)

Equation (61) is satisfied. In terms of ${\beta }_k$, Equation (56) is recast to

$${\mathbf{x}}_{k+1}={\mathbf{x}}_k+\frac{{\left({\mathbf{r}}_k^{\mathrm{T}}{\mathbf{v}}_k\right)}^2}{\parallel {\mathbf{r}}_k{\parallel }^2{\parallel {\mathbf{v}}_k\parallel }^2}\frac{\parallel {\mathbf{r}}_k{\parallel }^2}{{\mathbf{r}}_k^{\mathrm{T}}{\mathbf{v}}_k}{\mathbf{u}}_k={\mathbf{x}}_k+\frac{{\mathbf{r}}_k^{\mathrm{T}}{\mathbf{v}}_k}{\parallel {\mathbf{v}}_k{\parallel }^2}{\mathbf{u}}_k.$$

(63)

Noticing ${\mathbf{v}}_k=\mathbf{A}{\mathbf{u}}_k$, Equation (63) is just the iterative Equation (18) for ROGMRES.    □

Remark 3.

For the factor defined in Equation (23), ${\eta }_k$ plays a vital role in Equation (34) to dominate the orthogonal property of GMRES, and also in Equation (35) to control the residual decreasing property of GMRES. When GMRES blows up for some problems, by introducing ${\eta }_k$ in GMRES, we can stabilize the iterations to avoid blowing up. In this viewpoint, ${\eta }_k$ is a stabilization factor for GMRES.

Inspired by Remark 3, two simple methods can be developed for GMRES and GMRES(m): we can compute ${\eta }_k$ at each step, and consider $|{\eta }_k-1|>{10}^{-10}$ or ${\eta }_k<1/2$ to be the stop** criteria of the iterations; they are labeled as GMRES($\eta$) and GMRES(m,$\eta$), respectively.

Remark 4.

h given in Equation (21) is a decreasing quantity of the residual, whose maximal value by means of Equation (36) is

$$maxh=\frac{{\left({\mathbf{r}}_k^T\mathbf{A}{\mathbf{u}}_k\right)}^2}{\parallel \mathbf{A}{\mathbf{u}}_k{\parallel }^2}.$$

(64)

Upon comparing with $a_k$ in Equation (60), we have

$$maxh=\frac{\parallel {\mathbf{r}}_k{\parallel }^2}{a_k}.$$

(65)

Because $a_k\ge 1$, the best value of $maxh$ that can be obtained is $\parallel {\mathbf{r}}_k{\parallel }^2$ when $a_k=1$; in this situation, it will directly lead to the exact solution at the $k+1$ step, owing to $\parallel {\mathbf{r}}_{k+1}\parallel =0$, which is obtained by inserting $\mathbf{A}{\mathbf{u}}_k={\mathbf{r}}_k$ into Equation (38). In this sense, $a_k$ can also be an objective function. $a_0$ correlates to $f={\left({\mathbf{r}}_0^{\mathrm{T}}\mathbf{A}\mathbf{u}\right)}^2/{\parallel \mathbf{A}\mathbf{u}\parallel }^2$, which is the quantity of the projection from ${\mathbf{r}}_0$ on ${\mathcal{L}}_m=\mathbf{A}{\mathcal{K}}_m$, by

$$a_0=\frac{\parallel {\mathbf{r}}_0{\parallel }^2}{f}.$$

(66)

In the next section, we will maximize f, i.e., minimize $a_0$, to find the best descent vector for $\mathbf{u}$.

4. Orthogonalized Maximal Projection Algorithm

As seen from Equation (38), we must make ${\parallel \mathbf{A}\mathbf{u}\parallel }^2$ as large as possible. $\mathbf{A}\mathbf{u}=\mathbf{A}\mathbf{V}\mathbf{y}$, mentioned in Section 1, is the projection of ${\mathbf{r}}_0$ on ${\mathcal{L}}_m=\mathbf{A}{\mathcal{K}}_m$. Therefore we consider

$$\underset{\mathbf{u}}{max}\left\{f=\frac{{({\mathbf{r}}_0·\mathbf{v})}^2}{{\parallel \mathbf{v}\parallel }^2}\right\}$$

(67)

to be the maximal projection of ${\mathbf{r}}_0$ on the direction $\mathbf{v}/\parallel \mathbf{v}\parallel$, where $\mathbf{v}=\mathbf{A}\mathbf{u}$. In view of Equation (66), the maximization of f is equivalent to the minimization of $a_0$.

We can expand $\mathbf{u}\in {\mathbf{r}}_0+{\mathcal{K}}_m$ via

$$\mathbf{u}={\mathbf{r}}_0+\mathbf{V}\mathbf{y},$$

(68)

which is slightly different from the descent vector $\mathbf{u}=\mathbf{V}\mathbf{y}\in {\mathcal{K}}_m$ used in GMRES.

Theorem 4.

For $\mathbf{u}\in {\mathbf{r}}_0+{\mathcal{K}}_m(\mathbf{A},{\mathbf{r}}_0)$, ${\mathbf{r}}_0\ne \mathbf{0}$, the optimal $\mathbf{u}$ approximately satisfying Equation (2) and subjecting to Equation (67) is given in Equation (68), where

$$\begin{array}{ccc}\hfill & & \mathbf{W}=\mathbf{A}\mathbf{V},\mathbf{C}={\mathbf{W}}^T\mathbf{W},\hfill \end{array}$$

(69)

$$\begin{array}{ccc}\hfill & & \mathbf{C}\mathbf{y}={\mathbf{W}}^T({\mathbf{r}}_0-\mathbf{A}{\mathbf{r}}_0).\hfill \end{array}$$

(70)

Proof. 

Due to $\mathbf{v}=\mathbf{A}\mathbf{u}$ and Equation (68),

$$\mathbf{v}={\mathbf{v}}_0+\mathbf{W}\mathbf{y},$$

(71)

where $\mathbf{W}=\mathbf{A}\mathbf{V}$ and

$${\mathbf{v}}_0=\mathbf{A}{\mathbf{r}}_0.$$

(72)

By Equation (71), one has

$$\begin{array}{ccc}\hfill & & {\mathbf{r}}_0·\mathbf{v}={\mathbf{r}}_0·{\mathbf{v}}_0+{\mathbf{r}}_0^{\mathrm{T}}\mathbf{W}\mathbf{y},\hfill \end{array}$$

(73)

$$\begin{array}{ccc}\hfill & & {\parallel \mathbf{v}\parallel }^2={\parallel {\mathbf{v}}_0\parallel }^2+2{\mathbf{v}}_0^{\mathrm{T}}\mathbf{W}\mathbf{y}+{\mathbf{y}}^{\mathrm{T}}\mathbf{C}\mathbf{y},\hfill \end{array}$$

(74)

where

$$\mathbf{C}:={\mathbf{W}}^{\mathrm{T}}\mathbf{W}>\mathbf{0}\in {\mathbb{R}}^{m\times m}.$$

(75)

Resorting to the maximality condition for f in Equation (67), we can obtain

$${\mathbf{r}}_0·\mathbf{v}{\mathbf{y}}_2-{\parallel \mathbf{v}\parallel }^2{\mathbf{y}}_1=\mathbf{0},$$

(76)

where

$$\begin{array}{ccc}\hfill & & {\mathbf{y}}_1:={\nabla }_{\mathbf{y}}({\mathbf{r}}_0·\mathbf{v})={\mathbf{W}}^{\mathrm{T}}{\mathbf{r}}_0,\hfill \end{array}$$

(77)

$$\begin{array}{ccc}\hfill & & {\mathbf{y}}_2:=\frac{1}{2}{\nabla }_{\mathbf{y}}{\parallel \mathbf{v}\parallel }^2={\mathbf{W}}^{\mathrm{T}}{\mathbf{v}}_0+\mathbf{C}\mathbf{y}.\hfill \end{array}$$

(78)

From Equation (76), ${\mathbf{y}}_2$ is equal to ${\mathbf{y}}_1$:

$${\mathbf{y}}_2={\eta }^{-1}{\mathbf{y}}_1\Rightarrow {\mathbf{y}}_2={\mathbf{y}}_1,$$

(79)

because of

$${\eta }^{-1}=\frac{{\parallel \mathbf{v}\parallel }^2}{{\mathbf{r}}_0·\mathbf{v}}=\frac{{\parallel \mathbf{A}\mathbf{u}\parallel }^2}{{\mathbf{r}}_0^{\mathrm{T}}\mathbf{A}\mathbf{u}}=1,$$

(80)

according to Corollary 1.

Then, it follows from Equations (77)–(79) that

$$\mathbf{C}\mathbf{y}+{\mathbf{W}}^{\mathrm{T}}{\mathbf{v}}_0={\mathbf{W}}^{\mathrm{T}}{\mathbf{r}}_0,$$

(81)

which can be arranged to Equation (70).    □

We pull back the minimization problem used in the GMRES method to the following one:

$$min\{\parallel {\mathbf{r}}_0{-\mathbf{v}\parallel }^2\}=min\{\parallel {\mathbf{r}}_0-\mathbf{A}\mathbf{u}{\parallel }^2\},$$

(82)

which is derived from Equations (6) and (7).

Theorem 5.

For $\mathbf{u}\in {\mathbf{r}}_0+{\mathcal{K}}_m(\mathbf{A},{\mathbf{r}}_0)$, derived from the minimization problem in Equation (82), $\mathbf{u}$ in Equation (2) is optimized to be that in Theorem 4.

Proof. 

From

$${\nabla }_{\mathbf{y}}\parallel {\mathbf{r}}_0{-\mathbf{v}\parallel }^2={\nabla }_{\mathbf{y}}{\parallel \mathbf{v}\parallel }^2-2{\nabla }_{\mathbf{y}}({\mathbf{r}}_0·\mathbf{v})=\mathbf{0},$$

(83)

and Equations (78) and (77) it follows that

$${\mathbf{W}}^{\mathrm{T}}{\mathbf{v}}_0+\mathbf{C}\mathbf{y}={\mathbf{W}}^{\mathrm{T}}{\mathbf{r}}_0,$$

(84)

which is just Equation (70).    □

According to Theorems 4 and 5, the maximal projection is equivalent to the minimization of residual. Therefore, the maximal projection algorithm (MPA) is an extension of GMRES to a larger space of the affine Krylov subspace, not merely in the Krylov subspace.

By considering Lemma 1 and Theorem 5, we come to an improvement of ROGMRES as well as the double-improvement of GMRES to the following iterative formulas for the orthogonalized maximal projection algorithm (OMPA):

$$\begin{array}{ccc}\hfill & & \left({\mathbf{W}}^{\mathrm{T}}\mathbf{W}\right)\mathbf{y}={\mathbf{W}}^{\mathrm{T}}({\mathbf{r}}_0-\mathbf{A}{\mathbf{r}}_0),\hfill \end{array}$$

(85)

$$\begin{array}{ccc}\hfill & & \mathbf{u}={\mathbf{r}}_0+\mathbf{V}\mathbf{y},\hfill \end{array}$$

(86)

$$\begin{array}{ccc}\hfill & & \mathbf{x}={\mathbf{x}}_0+\frac{{\mathbf{r}}_0^{\mathrm{T}}\mathbf{A}\mathbf{u}}{{\parallel \mathbf{A}\mathbf{u}\parallel }^2}\mathbf{u}.\hfill \end{array}$$

(87)

Equation (85) can be further simplified in the next section. The algorithms of OMPA (Algorithm 6) and OMPA(m) (Algorithm 7) are given as follows.

Algorithm 6: OMPA

1: Give $m_0$, ${\mathbf{x}}_0$, and $\epsilon$ 2: Do $m=m_0,\dots$ ($k=m-m_0$), until $\parallel {\mathbf{r}}_k\parallel <\epsilon$ 3: ${\mathbf{r}}_k=\mathbf{b}-\mathbf{A}{\mathbf{x}}_k$ 4: ${\mathbf{v}}_1^k=\frac{{\mathbf{r}}_k}{\parallel {\mathbf{r}}_k\parallel }$ 5: ${\mathbf{V}}_k=[{\mathbf{v}}_1^k,\dots ,{\mathbf{v}}_m^k]$  (by Arnoldi process) 6: ${\mathbf{W}}_k=\mathbf{A}{\mathbf{V}}_k$ 7: Solve $\left({\mathbf{W}}_k^{\mathrm{T}}{\mathbf{W}}_k\right){\mathbf{y}}_k={\mathbf{W}}_k^{\mathrm{T}}({\mathbf{r}}_k-\mathbf{A}{\mathbf{r}}_k)$ 8: ${\mathbf{u}}_k={\mathbf{r}}_k+{\mathbf{V}}_k{\mathbf{y}}_k$ 9: ${\mathbf{x}}_{k+1}={\mathbf{x}}_k+\frac{{\mathbf{r}}_k^{\mathrm{T}}\mathbf{A}{\mathbf{u}}_k}{\parallel \mathbf{A}{\mathbf{u}}_k{\parallel }^2}{\mathbf{u}}_k$

Algorithm 7: OMPA(m)

1: Give $m_0$, $m_1$, ${\mathbf{x}}_0$, and $\epsilon$ 2: Do $j=1,\dots$ 3: Do $m=m_0,\dots ,m_1$ ($k=m-m_0$) 4: ${\mathbf{r}}_k=\mathbf{b}-\mathbf{A}{\mathbf{x}}_k$ 5: ${\mathbf{v}}_1^k=\frac{{\mathbf{r}}_k}{\parallel {\mathbf{r}}_k\parallel }$ 6: ${\mathbf{V}}_k=[{\mathbf{v}}_1^k,\dots ,{\mathbf{v}}_m^k]$  (by Arnoldi process) 7: ${\mathbf{W}}_k=\mathbf{A}{\mathbf{V}}_k$ 8: Solve $\left({\mathbf{W}}_k^{\mathrm{T}}{\mathbf{W}}_k\right){\mathbf{y}}_k={\mathbf{W}}_k^{\mathrm{T}}({\mathbf{r}}_k-\mathbf{A}{\mathbf{r}}_k)$ 9: ${\mathbf{u}}_k={\mathbf{r}}_k+{\mathbf{V}}_k{\mathbf{y}}_k$ 10: ${\mathbf{x}}_{k+1}={\mathbf{x}}_k+\frac{{\mathbf{r}}_k^{\mathrm{T}}\mathbf{A}{\mathbf{u}}_k}{\parallel \mathbf{A}{\mathbf{u}}_k{\parallel }^2}{\mathbf{u}}_k$ 11: If $\parallel {\mathbf{r}}_{k+1}\parallel <\epsilon$, stop 12: Otherwise, ${\mathbf{x}}_0={\mathbf{x}}_{k+1}$; go to 2

In Algorithm OMPA(m), if we take $m_1=m_0$, then we return to the usual iterative algorithm with a fixed-dimension $m=m_0$ of the Krylov subspace. In Algorithm OMPA(m), $m=m_1-m_0+1$ is the frequency for restart.

Remark 5.

In view of Equations (67) and (80), f and $\eta =1$ have the following relation:

$$f=\frac{{\mathbf{r}}_0·\mathbf{v}}{\eta }={\mathbf{r}}_0·\mathbf{v}={\parallel \mathbf{A}\mathbf{u}\parallel }^2.$$

(88)

Therefore, the maximization of f is equivalent to the maximization of the length of the projection vector $\mathbf{A}\mathbf{u}$. By means of Equation (64),

$$f=maxh;$$

(89)

hence,

$$\underset{\mathbf{u}}{max}f=\underset{\mathbf{u}}{max}\underset{\eta }{max}\left\{2\eta {\mathbf{r}}_0^T\mathbf{A}\mathbf{u}-{\eta }^2{\parallel \mathbf{A}\mathbf{u}\parallel }^2\right\}$$

(90)

follows from Equation (22). Through a double-maximization, we have derived the maximal projection algorithm.

5. An Affine GMRES Method

In the derivation of GMRES, the following technique was employed [11]:

$$\underset{\mathbf{y}}{min}\parallel r_0{\mathbf{e}}_1-{\overline{\mathbf{H}}}_m\mathbf{y}\parallel ,$$

(91)

where $r_0=\parallel {\mathbf{r}}_0\parallel$.

Instead of Equation (91) used in the GMRES method to develop the iterative algorithm from $\mathbf{u}\in {\mathcal{K}}_m(\mathbf{A},{\mathbf{r}}_0)$, we seek the descent vector, $\mathbf{u}$, in a larger subspace of $\mathbf{u}\in {\mathbf{r}}_0+{\mathcal{K}}_m(\mathbf{A},{\mathbf{r}}_0)$ by

$$\underset{\mathbf{y}}{min}\parallel {\mathbf{r}}_0-\mathbf{A}{\mathbf{r}}_0-\mathbf{W}\mathbf{y}\parallel ,$$

(92)

which is obtained from Equation (82) by inserting Equation (86) for $\mathbf{u}$, and using $\mathbf{W}=\mathbf{A}\mathbf{V}$.

Now, Equation (92) is re-written as

$$\underset{\mathbf{y}}{min}\parallel {\mathbf{r}}_0-\mathbf{A}{\mathbf{r}}_0-\mathbf{A}{\mathbf{V}}_m\mathbf{y}\parallel .$$

(93)

We denote ${\mathbf{V}}_m$ the $m\times m$-dimensional Arnoldi matrix, and ${\mathbf{V}}_{m+1}$ the $(m+1)\times (m+1)$-dimensional Arnoldi matrix. It is known that

$$\mathbf{A}{\mathbf{V}}_m={\mathbf{V}}_{m+1}{\overline{\mathbf{H}}}_m,{\mathbf{V}}_m^{\mathrm{T}}{\mathbf{V}}_m={\mathbf{I}}_m,{\mathbf{V}}_{m+1}^{\mathrm{T}}{\mathbf{V}}_{m+1}={\mathbf{I}}_{m+1}.$$

(94)

Theorem 6.

For $\mathbf{u}\in {\mathbf{r}}_0+{\mathcal{K}}_m(\mathbf{A},{\mathbf{r}}_0)$, derived from Equation (93), the orthogonalized affine GMRES possesses the following iterative formulas:

$$\begin{array}{ccc}\hfill & & \left({\overline{\mathbf{H}}}_m^T{\overline{\mathbf{H}}}_m\right)\mathbf{y}=r_0{\overline{\mathbf{H}}}_m^T{\mathbf{e}}_1-r_0{\left[{\overline{\mathbf{H}}}_m^T{\overline{\mathbf{H}}}_m\right]}_1,\hfill \end{array}$$

(95)

$$\begin{array}{ccc}\hfill & & \mathbf{u}={\mathbf{r}}_0+{\mathbf{V}}_m\mathbf{y},\hfill \end{array}$$

(96)

$$\begin{array}{ccc}\hfill & & \mathbf{x}={\mathbf{x}}_0+\frac{{\mathbf{r}}_0^T\mathbf{A}\mathbf{u}}{{\parallel \mathbf{A}\mathbf{u}\parallel }^2}\mathbf{u},\hfill \end{array}$$

(97)

where $r_0=\parallel {\mathbf{r}}_0\parallel$, and ${\left[{\overline{\mathbf{H}}}_m^T{\overline{\mathbf{H}}}_m\right]}_1$ is the first column of ${\overline{\mathbf{H}}}_m^T{\overline{\mathbf{H}}}_m$.

Proof. 

The objective ${\mathbf{r}}_0-\mathbf{A}{\mathbf{r}}_0-\mathbf{A}{\mathbf{V}}_m\mathbf{y}$ in Equation (93) can be simplified as follows. In terms of the first Krylov vector ${\mathbf{v}}_1={\mathbf{r}}_0/r_0$, where $r_0=\parallel {\mathbf{r}}_0\parallel$, and using the first one in Equation (94), we have

$${\mathbf{r}}_0-\mathbf{A}{\mathbf{r}}_0-\mathbf{A}{\mathbf{V}}_m\mathbf{y}=r_0{\mathbf{v}}_1-r_0\mathbf{A}{\mathbf{v}}_1-{\mathbf{V}}_{m+1}{\overline{\mathbf{H}}}_m\mathbf{y}.$$

(98)

We write

$${\mathbf{v}}_1={\mathbf{V}}_{m+1}{\mathbf{e}}_1,\mathbf{A}{\mathbf{v}}_1=\mathbf{A}{\mathbf{V}}_m{\mathbf{e}}_2={\mathbf{V}}_{m+1}{\overline{\mathbf{H}}}_m{\mathbf{e}}_2,$$

(99)

where ${\mathbf{e}}_2$ is the first column of ${\mathbf{I}}_m$. Then, it follows that

$${\mathbf{r}}_0-\mathbf{A}{\mathbf{r}}_0-\mathbf{A}{\mathbf{V}}_m\mathbf{y}={\mathbf{V}}_{m+1}[r_0{\mathbf{e}}_1-r_0{\overline{\mathbf{H}}}_m{\mathbf{e}}_2-{\overline{\mathbf{H}}}_m\mathbf{y}].$$

(100)

Taking the squared norms of both sides yields

$$\begin{array}{ccc}\hfill & & \parallel {\mathbf{r}}_0-\mathbf{A}{\mathbf{r}}_0-\mathbf{A}{\mathbf{V}}_m{\mathbf{y}\parallel }^2\hfill \\ & & ={[r_0{\mathbf{e}}_1-r_0{\overline{\mathbf{H}}}_m{\mathbf{e}}_2-{\overline{\mathbf{H}}}_m\mathbf{y}]}^{\mathrm{T}}{\mathbf{V}}_{m+1}^{\mathrm{T}}{\mathbf{V}}_{m+1}[r_0{\mathbf{e}}_1-r_0{\overline{\mathbf{H}}}_m{\mathbf{e}}_2-{\overline{\mathbf{H}}}_m\mathbf{y}]\hfill \\ & & ={[r_0{\mathbf{e}}_1-r_0{\overline{\mathbf{H}}}_m{\mathbf{e}}_2-{\overline{\mathbf{H}}}_m\mathbf{y}]}^{\mathrm{T}}[r_0{\mathbf{e}}_1-r_0{\overline{\mathbf{H}}}_m{\mathbf{e}}_2-{\overline{\mathbf{H}}}_m\mathbf{y}]\hfill \\ & & =\parallel r_0{\mathbf{e}}_1-r_0{\overline{\mathbf{H}}}_m{\mathbf{e}}_2-{\overline{\mathbf{H}}}_m{\mathbf{y}\parallel }^2,\hfill \end{array}$$

(101)

where the last one in Equation (94) was used.

The minimization problem in Equation (93) is now simplified to

$$\underset{\mathbf{y}}{min}\parallel r_0{\mathbf{e}}_1-r_0{\overline{\mathbf{H}}}_m{\mathbf{e}}_2-{\overline{\mathbf{H}}}_m\mathbf{y}\parallel .$$

(102)

To satisfy the minimality condition, we can derive Equation (95); using the technique in Lemma 1, Equation (97) can be derived.    □

Remark 6.

Notice that the minimization problem in Equation (102) is slightly different from that in Equation (91), where an extra term, $r_0{\overline{\mathbf{H}}}_m{\mathbf{e}}_2$, appeared in Equation (102). However, the resulting affine GMRES (A-GMRES) (Algorithm 8) in the affine Krylov subspace can significantly improve the performance rather than GMRES, when the orthogonality is taken into account. We remind that (6)–(8) in the following Algorithm A-GMRES are slightly more computational expensive than (6) and (7) in Algorithm GMRES. Also, (7)–(9) in the following Algorithm A-GMRES(m) (Algorithm 9) are slightly more computational expensive than (7) and (8) in Algorithm GMRES(m).

Algorithm 8: A-GMRES

1: Give $m_0$, ${\mathbf{x}}_0$, and $\epsilon$ 2: Do $m=m_0,\dots$ ($k=m-m_0$) 3: ${\mathbf{r}}_k=\mathbf{b}-\mathbf{A}{\mathbf{x}}_k$, $r_k=\parallel {\mathbf{r}}_k\parallel$ 4: ${\mathbf{v}}_1^k=\frac{{\mathbf{r}}_k}{\parallel {\mathbf{r}}_k\parallel }$ 5: ${\mathbf{V}}_k=[{\mathbf{v}}_1^k,\dots ,{\mathbf{v}}_m^k]$  (by Arnoldi process) 6: $\left({\overline{\mathbf{H}}}_k^{\mathrm{T}}{\overline{\mathbf{H}}}_k\right){\mathbf{y}}_k=r_k{\overline{\mathbf{H}}}_k^{\mathrm{T}}{\mathbf{e}}_1-r_k{\left[{\overline{\mathbf{H}}}_k^{\mathrm{T}}{\overline{\mathbf{H}}}_k\right]}_1$ 7: ${\mathbf{u}}_k={\mathbf{r}}_k+{\mathbf{V}}_k{\mathbf{y}}_k$ 8: ${\mathbf{x}}_{k+1}={\mathbf{x}}_k+\frac{{\mathbf{r}}_k^{\mathrm{T}}\mathbf{A}{\mathbf{u}}_k}{\parallel \mathbf{A}{\mathbf{u}}_k{\parallel }^2}{\mathbf{u}}_k$ 9: If $\parallel {\mathbf{r}}_{k+1}\parallel <\epsilon$, stop

Algorithm 9: A-GMRES(m)

1: Give $m_0$, $m_1$, ${\mathbf{x}}_0$, and $\epsilon$ 2: Do $j=1,\dots$ 3: Do $m=m_0,\dots ,m_1$ ($k=m-m_0$) 4: ${\mathbf{r}}_k=\mathbf{b}-\mathbf{A}{\mathbf{x}}_k$, $r_k=\parallel {\mathbf{r}}_k\parallel$ 5: ${\mathbf{v}}_1^k=\frac{{\mathbf{r}}_k}{\parallel {\mathbf{r}}_k\parallel }$ 6: ${\mathbf{V}}_k=[{\mathbf{v}}_1^k,\dots ,{\mathbf{v}}_m^k]$  (by Arnoldi process) 7: $\left({\overline{\mathbf{H}}}_k^{\mathrm{T}}{\overline{\mathbf{H}}}_k\right){\mathbf{y}}_k=r_k{\overline{\mathbf{H}}}_k^{\mathrm{T}}{\mathbf{e}}_1-r_k{\left[{\overline{\mathbf{H}}}_k^{\mathrm{T}}{\overline{\mathbf{H}}}_k\right]}_1$ 8: ${\mathbf{u}}_k={\mathbf{r}}_k+{\mathbf{V}}_k{\mathbf{y}}_k$ 9: ${\mathbf{x}}_{k+1}={\mathbf{x}}_k+\frac{{\mathbf{r}}_k^{\mathrm{T}}\mathbf{A}{\mathbf{u}}_k}{\parallel \mathbf{A}{\mathbf{u}}_k{\parallel }^2}{\mathbf{u}}_k$ 10: If $\parallel {\mathbf{r}}_{k+1}\parallel <\epsilon$, stop 11: Otherwise, ${\mathbf{x}}_0={\mathbf{x}}_{k+1}$; go to 2

6. Examples

In this section, we apply the Algorithms GMRES, GMRES(m), A-GMRES, A-GMRES(m), ROGMRES, ROGMRES(m), OMPA, and MPA(m) to solve some examples endowed with a symmetric cyclic matrix, randomized matrix, ill-posed Cauchy problem, diagonal dominant matrix, and highly ill-conditioned matrix. They are subjected to the convergence criterion with $\parallel {\mathbf{r}}_k\parallel <\epsilon$, where $\epsilon$ is the error tolerance.

6.1. Example 1: Cyclic Matrix Linear System

A cyclic matrix, $\mathbf{A}$, is generated from the first row $(1,\dots ,n)$. The exact solution is assumed to be $x_i=1,i=1,\dots ,n$, and $\mathbf{b}$ can be obtained from Equation (1). First we apply GMRES to this problem with $n=300$ and $\epsilon ={10}^{-10}$. The initial values are $x_i^0=1+1/i,i=1,\dots ,n$. For this problem, we encounter the situation that ${\eta }_k<1$ after 150 steps, not ${\eta }_k=1$ as shown in Figure 1a by a solid line. As shown in Figure 1b, the residuals blow up for GMRES.

This problem shows that the original GMRES method cannot be applied to solve this linear system, since ${\eta }_k<1$ happened after 150, and the iteration of GMRES blows up after around 157 steps. In this situation, by Equation (35),

$$\parallel {\mathbf{r}}_{k+1}\parallel >\parallel {\mathbf{r}}_k\parallel ,$$

the residual obtained by GMRES grows step-by step; hence, the GMRES method is no longer stable after 157 steps.

To overcome this drawback of GMRES, we can either relax the convergence criterion, or employ the version of GMRES$\left(\eta \right)$ before it diverges.

Table 1. For example 1, ME and NS obtained by different methods.

	GMRES(m)	ROGMRES	ROGMRES(m)	GMRES($\mathit{\eta }$)	GMRES(m,$\mathit{\eta }$)	MPA(m)	A-GMRES(m)
ME	$7.39\times {10}^{-14}$	$7.22\times {10}^{-14}$	$6.50\times {10}^{-14}$	$8.22\times {10}^{-14}$	$9.04\times {10}^{-14}$	$5.66\times {10}^{-14}$	$5.88\times {10}^{-14}$
NS	28	44	18	16	16	10	15

compares the maximum error (ME) and number of steps (NS) for different methods under $\epsilon =3\times {10}^{-10}$. $m_0=10$ is used in GMRES and ROGMRES; $m_0=1$ is used in GMRES($\eta$); $m_0=30$ and $m_1=40$ are used in GMRES(m); $m_0=35$ and $m_1=40$ are used in ROGMRES(m) and GMRES(m,$\eta$); $m_0=50$ and $m_1=50$ are used in MPA(m) and A-GMRES(m). As shown in Figure 1a by a dashed line and a dashed–dotted line, both ROGMRES and ROGMRES(m) can keep the value of ${\eta }_k=1$ up to the termination; as shown in Figure 1b, they converge faster than GMRES. ROGMRES(m) converges faster than ROGMRES. This example shows that ROGMRES can stabilize GMRES when it is unstable, as shown in Figure 1b.

It is interesting that the MPA(m), with a fixed value $m=m_0=m_1=50$, can attain a highly accurate solution with ten steps, as shown in Figure 1b; A-GMRES(m) obtains ME = $5.88\times {10}^{-14}$ with 15 steps by using $m=m_0=m_1=50$. If the orthogonality is not adopted in A-GMRES(m), ME = $7.51\times {10}^{-14}$ is obtained with 14 steps; it means that ${\eta }_k=1$ can be kept, even though the orthogonalization technique was not employed. Within ten steps, A-GMRES with $m_0=50$ can obtain ME = $7.06\times {10}^{-14}$. Within nine steps, A-GMRES(m) with $m_0=50$ and $m_0=51$ can obtain ME = $8.02\times {10}^{-14}$.

6.2. Example 2: Random Matrix Linear System

The coefficient matrix $\mathbf{A}$ of size $50\times 50$ is randomly generated with $1<a_{ij}<2$. The exact solution is $x_i=1,i=1,\dots ,50$, and $\mathbf{b}$ can be obtained by Equation (1).

We consider the initial values $x_i^0=0.5$ and take $\epsilon ={10}^{-10}$.

Table 2. For example 2, ME and NS obtained by different methods.

	GMRES(m)	ROGMRES	OMPA	A-GMRES	ROGMRES(m)	GMRES($\mathit{\eta }$)	GMRES(m,$\mathit{\eta }$)
ME	$3.51\times {10}^{-2}$	$9.77\times {10}^{-14}$	$1.88\times {10}^{-13}$	$3.36\times {10}^{-9}$	$3.51\times {10}^{-2}$	$3.65\times {10}^{-2}$	$3.25\times {10}^{-2}$
NS	>550	51	43	52	>550	45	163

compares ME and NS for different methods under $\epsilon ={10}^{-10}$. $m_0=10$ is used in ROGMRES, OMPA, A-GMRES, and GMRES($\eta$); $m_0=30$ and $m_1=40$ are used in GMRES(m), ROGMRES(m) and GMRES(m,$\eta$). ROGMRES is better than OMPA and A-GMRES; other methods are not accurate. For this example, the restart ROGMRES(m) is not good.

6.3. Example 3: Inverse Cauchy Problem and Boundary Value Problem

Consider the inverse Cauchy problem:

$$\begin{array}{ccc}\hfill & & \Delta u=u_{rr}+\frac{1}{r}{u}_r+\frac{1}{r^2}{u}_{\theta \theta }=0,\hfill \end{array}$$

(103)

$$\begin{array}{ccc}\hfill & & u(\rho ,\theta )=h\left(\theta \right),0\le \theta \le \pi ,\hfill \end{array}$$

(104)

$$\begin{array}{ccc}\hfill & & u_n(\rho ,\theta )=g\left(\theta \right),0\le \theta \le \pi .\hfill \end{array}$$

(105)

The method of fundamental solutions is taken as

$$u\left(\mathbf{x}\right)=\sum _{j=1}^n{c}_{j}U(\mathbf{x},{\mathbf{s}}_j),{\mathbf{s}}_j\in {\Omega }^c.$$

(106)

where

$$U(\mathbf{x},{\mathbf{s}}_j)=ln{r}_j,r_j=\parallel \mathbf{x}-{\mathbf{s}}_j\parallel .$$

(107)

We consider

$$\begin{array}{ccc}\hfill & & u(x,y)=cosxcoshy+sinxsinhy,\hfill \end{array}$$

(108)

$$\begin{array}{ccc}\hfill & & \rho \left(\theta \right)=\exp \left(\sin \theta \right){\sin }^2\left(2\theta \right)+\exp \left(\cos \theta \right){\cos }^2\left(2\theta \right),\hfill \end{array}$$

(109)

where $u(x,y)$ is an exact solution of the Laplace equation, and $\rho \left(\theta \right),0\le \theta \le 2\pi$ describes the boundary contour in the polar coordinates.

By GMRES with $n=60$, through 60 iterations the residual is shown in Figure 2a, which does not satisfy the convergence criterion with $\epsilon ={10}^{-2}$. If we run it to 70 steps, GMRES will blow up. In contrast, ROGMRES converges with 72 steps. As compared in Figure 2b, ME obtained by ROGMRES is smaller than 0.27, as shown in Figure 2c. However, GMRES obtains ME = 0.5 with a worse result.

Because the inverse Cauchy problem is a highly ill-posed problem, GMRES does not converge within a loose convergence criterion $\epsilon ={10}^{-2}$, and after 70 steps it will blow up. It is apparent that GMRES cannot solve the Cauchy problem adequately. The ROGMRES method is stable for the Cauchy problem, and can improve the accuracy to 0.27; however, it leaves room to improve the methods based on the Krylov subspace for solving the Cauchy problem.

Next, we consider the mixed boundary value problem. Under $\epsilon ={10}^{-3}$, for GMRES(m) with $m_0=3$ and $m_1=10$, as shown in Figure 3a, ${\eta }_k$ after 12 steps is irregular; it causes the residuals to decrease non-monotonically, as shown in Figure 3b. The GMRES(m) does not converge within 200 steps, and obtains an incorrect solution with ME = 0.92. If ${\eta }_k<1/2$ is used as a stop** criterion, through 14 steps, ME = $6.64\times {10}^{-3}$ is obtained. When the ROGMRES method is applied, we obtain ME = $4.16\times {10}^{-4}$ with 187 steps. As shown in Figure 3a by a dashed line, ${\eta }_k=1$ is kept well, and the residuals monotonically decrease, as shown in Figure 3b by a dashed line for ROGMRES. ME = $4.87\times {10}^{-4}$ is obtained through 167 steps by OMPA with $m_0=4$, which is convergent faster than ROGMRES. ME = $4.09\times {10}^{-4}$ is obtained through 109 steps by A-GMRES with $m_0=4$, which is convergent faster than ROGMRES and OMPA.

6.4. Example 4: Diagonal Dominant Linear System

Consider an example borrowed from Morgan [27]. The matrix $\mathbf{A}$ consists of main diagonal elements from 1 to 1000; the super diagonal elements are 0.1, and other elements are zero. We take $\mathbf{b}=\mathbf{1}$.

Table 3. For example 4, NS obtained by different methods.

	Original GMRES(m)	ROGMRES(m)	OMPA	A-GMRES	[27]	Current GMRES(m)
NS	>300	44	32	31	300	44

compares the number of steps (NS) for different methods under $\epsilon ={10}^{-10}$. $m_0=1$ is used in OMPA and A-GMRES; $m_0=1,m_1=25$ are used in original GMRES(m), ROGMRES(m) and current GMRES(m). The NS of ROGMRES(m), and current GMRES(m) are the same, because of ${\eta }_k=1$ for the current GMRES(m) before convergence. Figure 4 compares the residuals.

The method in Equation (18) is an automatically orthogonality preserving method. When the OMPA and A-GMRES methods do not consider the automatically orthogonality preserving method, the results are drastically different, as shown in Figure 5. Because $\eta =1$ is not preserved, both OMPA and A-GMRES diverge very quickly after 40 steps.

6.5. Example 5: Densely Clustered Non-Symmetric Linear System

Consider an example with $a_{ij}=(i+j/2)/n,i,j=1,\dots ,n$, where $n=2000$. Suppose that $x_i=1,i=1,\dots ,2000$ are exact solutions, and the initial values are $x_i^0=0,i=1,\dots ,2000$. The values of matrix elements are densely clustered in a narrow range with $a_{ij}\in [0.00075,1.5]$. This problem is a highly ill-conditioned non-symmetric linear system with the condition number over the order ${10}^{15}$.

Table 4. For example 5, ME and NS obtained by different methods.

Method	OMPA	A-GMRES	Current GMRES(m)	ROGMRES(m)
NS	4	3	9	2
ME	$2.10\times {10}^{-9}$	$1.46\times {10}^{-11}$	$1.18\times {10}^{-13}$	$5.98\times {10}^{-13}$

compares the maximal errors (MEs) and number of steps (NS) for different methods under $\epsilon ={10}^{-10}$. $m_0=1$ is used in OMPA and A-GMRES; $m_0=1,m_1=3$ are used in current GMRES(m), and $m_0=2,m_1=3$ are used in ROGMRES(m). ROGMRES(m) is convergent faster than other algorithms. Even for this highly ill-conditioned problem, the proposed novel methods are highly efficient in finding accurate solutions.

7. Conclusions

In this paper the GMRES method was re-examined from the two aspects of preserving the orthogonality and maximizing the decreasing length of the residual vector, both of which can significantly enhance the stability and accelerate the convergence speed of GMRES. If ${\eta }_k={\mathbf{r}}_k^{\mathrm{T}}\mathbf{A}{\mathbf{u}}_k/{\parallel \mathbf{A}{\mathbf{u}}_k\parallel }^2$ in GMRES is kept to ${\eta }_k=1$ during the iterations, it is stable; otherwise, the GMRES is unstable. For any iterative form of ${\mathbf{x}}_{k+1}={\mathbf{x}}_k+{\mathbf{u}}_k$, we can improve it to ${\mathbf{x}}_{k+1}={\mathbf{x}}_k+{\mathbf{u}}_k\left({\mathbf{r}}_k^{\mathrm{T}}\mathbf{A}{\mathbf{u}}_k\right)/{\parallel \mathbf{A}{\mathbf{u}}_k\parallel }^2$, such that the orthogonality is preserved automatically, and simultaneously the length of the residual vector is reduced maximally. It brings the GMRES method to a re-orthogonalized GMRES (ROGMRES) method, preserving the orthogonality automatically and having the property of absolute convergence. We derived the new iterative form ${\mathbf{x}}_{k+1}={\mathbf{x}}_k+{\mathbf{u}}_k\left({\mathbf{r}}_k^{\mathrm{T}}\mathbf{A}{\mathbf{u}}_k\right)/{\parallel \mathbf{A}{\mathbf{u}}_k\parallel }^2$ also from the invariant-manifold theory.

In order to find a better descent vector, we solved a maximal projection problem in a larger affine Krylov subspace to derive the descent vector. We proved that the maximal projection is equivalent to the minimal residual length in the m-dimensional affine Krylov subspace. After taking the automatically orthogonality preserving method into account, a new orthogonalized maximal projection algorithm (OMPA) method was developed. We derived formulas similar to GMRES in the affine Krylov subspace, namely the affine GMRES (A-GMRES) method, which had superior performance for the problems tested in the paper; it is even better than OMPA for some problems. The algorithm A-GMRES possessed two advantages: solving a simpler least squares problem in the affine Krylov subspace, and preserving the orthogonality automatically.

Through the studies conducted in the paper, the novelty points and main contributions are summarized as follows.

• The GMRES was examined from the viewpoints of maximizing the decreasing quantity of the residual and to maintain the orthogonality of the consecutive residual vector.
• A stabilization factor, $\eta$, to maintain orthogonality was inserted into GMRES to preserve the orthogonality automatically.
• GMRES was improved in a larger space of the affine Krylov subspace.
• A new orthogonalized maximal projection algorithm, OMPA, was derived; a new affine A-GMRES was derived.
• The new algorithms ROGMRES, OMPA, and A-GMRES guarantee the absolute convergence.
• Through examples, we showed that the orthogonalization techniques are useful to stabilize the methods of GMRES, MPA, and A-GMRES.
• Numerical testings for five different problems confirmed that the methods of ROGMRES, ROGMRES(m), OMPA, OMPA(m), A-GMRES, and A-GMRES(m) can significantly accelerate the convergence speed compared with the original methods of GMRES and GMRES(m).