# Residual-Based Simpler Block GMRES for Nonsymmetric Linear Systems with Multiple Right-Hand Sides

Authors: Qinghua Wu; Liang Bao; Yiqin Lin
Journal: Advances in Mathematical Physics
 (2018)

---

## Abstract

We propose in this paper a residual-based simpler block GMRES method for solving a system of linear algebraic equations with multiple right-hand sides. We show that this method is mathematically equivalent to the block GMRES method and thus equivalent to the simpler block GMRES method. Moreover, it is shown that the residual-based method is numerically more stable than the simpler block GMRES method. Based on the deflation strategy proposed by Calandra et al. (2013), we derive a deflation strategy to detect the possible linear dependence of the residuals and a near rank deficiency occurring in the block Arnoldi procedure. Numerical experiments are conducted to illustrate the performance of the new method.


---

## Body

## 1. Introduction

In this paper, we consider iterative methods for solving a system of linear algebraic equations:(1)AX=B,where A is a nonsingular matrix of order n and X=[x1,…,xs] and B=[b1,…,bs] are rectangular matrices of dimension n×s with s≤n. For solving such systems, the block GMRES [1] and its variants are very popular. Block GMRES is based on the block Arnoldi process and is formally fully analogous to the ordinary GMRES algorithm by Saad and Schultz [2].The following notation is used throughout the paper. Subscripts denote the iteration index and superscripts distinguish between individual columns in a block. We denote by· the Euclidean vector norm and the induced matrix norm and by AF the Frobenius norm. Moreover, for A∈RN×nN>n of rank n, κ(A)=σ1(A)/σn(A) is the spectral condition number, where σ1(A)>σn(A)>0 are the extremal singular values of A.Given an initial approximationX0∈Rn×s to the solution of (1), let(2)R0=B-AX0,and then in analogy to the unblocked case, we build a sequence of iterates Xm∈X0+Km(A,R0) such that(3)RmF=B-AXmF=minX~∈X0+KmA,R0⁡R0-AX~F,where Km(A,R0)=span{R0,AR0,…,Am-1R0}. Equation (3) is equivalent to minimizing every column of Rm, that is,(4)Rmi2=min⁡subject toXmi-X0i∈KmA,R0,and also to the orthogonality condition(5)Rmi⊥AKmA,R0,i=1,…,s,where ⊥ is the orthogonality relation induced by the Euclidean inner product. Assume that Xm is of the form Xm=X0+Ymtm, where Ym=[Y1,Y2,…,Ym] is a basis of Km(A,R0). Then, we obtain the mth residual matrix(6)Rm=R0-AYmtm.Central to the usual implementations of block GMRES is the block Arnoldi process [1], which can be used to construct the orthonormal basis of Km(A,R0). In practice, the possible linear dependence of the residuals of the s systems requires an explicit reduction of the number of right-hand sides. In [3], this was called deflation. If the block residual is nearly rank deficient, block GMRES should be implemented with deflation and there are various sophisticated rank-revealing QR factorizations. For details, see [3] and the references therein. We can write the nondeflated block Arnoldi process as shown in Algorithm 1.Algorithm 1:Nondeflated block Arnoldi process.
( 1 ) Given R0=B-AX0 of full rank. ( 2 ) Compute the QR factorization of R0:R0=Y1ρ0. ( 3 ) For j=1,…,m (3.1) Compute W=AYj. (3.2) For i=1,…,j Hi,j=YiTW. W=W-YiHi,j. End For (3.3) Compute the QR factorization of W:W=YjHj+1,j. End ForFrom Algorithm1, we obtain formally the ordinary Arnoldi relation(7)AYm=Ym+1Hm,where the (m+1)s×ms matrix is(8)Hm≔H1,1H1,2⋯H1,mH2,1H2,2⋯H2,mH3,2⋱⋮⋱Hm,mHm+1,m.In the block Arnoldi algorithm, YiTYi=I holds due to the QR factorizations, and YiTYj=0 when i≠j, where I is a unit matrix and 0 is a zero matrix of order s. This indicates that the whole process is equivalent to the one in which the block vectors are generated column by column using an ordinary modified Gram-Schmidt process.From (6) and (7), we obtain the fundamental block GMRES relation(9)Rm=Ym+1e1ρ0-Hmtm,where e1 is the first s column of the (m+1)s×(m+1)s unit matrix (the size changes with m), ρ0 is an upper triangular matrix obtained in Arnoldi’s initialization step, and tm is the “block coordinates” of Xm-X0 with respect to the block Arnoldi basis.Using (9), the least squares problem (3) is solved by recursive QR factorization of Hm, updated by applying Givens rotations. Once the norm of the residual is small enough, the triangular system with the computed R-factor is solved, and the approximate solution Xm is computed. The detailed algorithm of block GMRES can be found in [3–5].The block GMRES method with deflation at each iteration was proposed in [6]. And a deflation strategy was investigated to detect when a linear combination of approximate solutions is already known; for details, see [7]. In this paper, we deal with a different approach and compare the situation with deflation and without deflation. Let Zm be a block basis of Km(A,R0). Instead of building a block orthogonal basis of Km(A,R0), we look for a block orthogonal basis Vm=[V1,…,Vm] of AKm(A,R0). As a special case, simpler block GMRES (SBGMRES) was proposed by Liu and Zhong [8], where Zm=[R0/R0F,V1,…,Vm-1]. In this paper, we will consider the basis Zm=[R0/R0F,R1/R1F,…,Rm-1/Rm-1F]. We call this case the residual-based simpler block GMRES (RB-SBGMRES).The paper is organized as follows. In Section2, the RB-SBGMRES and SBGMRES algorithms are described (without deflation). In Section 3, some comparison between RB-SBGMRES and SBGMRES is established. In Section 4, the RB-SBGMRES method with deflation and the corresponding algorithm RB-SBGMRES-D (Table 2) are derived. In Section 5, these algorithms are compared using test matrices taken from the Matrix Market [9]. Conclusions are included in the last section.


## 2. Residual-Based Simpler Block GMRES

Suppose thatZm≡[Z1,…,Zm] is a basis of Km(A,R0). The orthogonal basis Vm of AKm(A,R0) is thus obtained from the QR factorization of AZm; that is,(10)AZm=VmUm,where Um is an upper triangular matrix with order ms.Due to the orthogonality propertyRm⊥AKmA,R0, the mth residual matrix Rm can be computed recursively as(11)Rm=R0-VmVmTR0=Rm-1-VmSm,where Sm=VmTRm-1.DefineDm=[S1T,…,SmT]T. Since the columns of Zm form a basis of Km(A,R0), we can represent Xm in the form(12)Xm=X0+Zmtm,where tm∈Rms×s is the “block coordinates” of Xm-X0 with respect to the block basis Zm. Due to Rm=R0-AZmtm, AZm=VmUm, and Rm⊥AKm(A,R0), it follows that(13)Umtm=VmTR0=Dm.Hence, once the residual norm is small enough, we can solve this upper triangular system (13) and then compute the approximate solution Xm.We now present the RB-SBGMRES method without deflation as shown in Algorithms2 and 3 and Table 1. For comparison, we also present the SBGMRES method proposed in [8].Table 1
Computational cost of a cycle of RB-SBGMRES (or SBGMRES).
Step Computational cost Computation ofR0 ns Computation ofWj n s Block Arnoldi procedure (2.1)-(2.2) O(ns2) QR factorization ofWj 2 n s 2 + n s Computation ofSj,Rj 2 s n Computation oftj O ( j s 3 ) Computation ofXj n j sTable 2
Computational cost of a cycle of RB-SBGMRES-D (Algorithm4).
Step Computational cost Computation ofR0 ns Computation ofG0 14 s 3 Computation ofWj n c j Computation ofF1 14 s 3 Block Arnoldi procedure (5.3) O(nkj2) Computation of[Vj,Pj-1] 2 n s 2 QR factorization ofR0 2 n s 2 + n s Computation of[Zj+1,Qj] 2 n s 2 Computation ofXj n ∑ i = 1 j c i QR factorization ofR0 2 n s 2 + n s Computation of[Z1,Q0] 2 n s 2 QR factorization ofW1 2 n c 1 2 + n c 1 Computation of[V1,P0] 2 n s 2 Computation ofFj 4 s j s 2 + 14 s 3 Computation ofSj,Rj 2 k j n Computation ofGj 4 s j s 2 + 14 s 3 Computation oftj O ( ∑ i = 1 j c i ) 3Algorithm 2:RB-SBGMRES(m).
( 1 ) Given X0, set R0=B-AX0. If R0F≤TOL, accept X0 and exit. ( 2 ) For j=1,…,m (2.1) Compute Wj=ARj-1. (2.2) For i=1,…,j-1 Uij=ViTWj. Wj=Wj-ViUij. End For (2.3) Compute the QR factorization of Wj:Wj=VjUjj. (2.4) Compute Sj=VjTRj-1, Rj=Rj-1-VjSj. (2.5) If RjF<TOL, break. End For ( 3 ) Solve the triangular system Ujtj=Dj for tj. ( 4 ) From the approximate solution Xj=X0+[R0,…,Rj-1]tj. ( 5 ) If RjF<TOL, then accept Xj and exit; otherwise, restart: set X0=Xj and go to Step (1).Algorithm 3:SBGMRES(m).
( 1 ) Given X0, set R0=B-AX0. If R0F≤TOL, accept X0 and exit. ( 2 ) For j=1,…,m (2.1) Compute Wj=AVj-1 if j>1, W1=AR0 for j=1. (2.2) For i=1,…,j-1 Uij=ViTWj. Wj=Wj-ViUij. End For (2.3) Compute the QR factorization of Wj:Wj=VjUjj. (2.4) Compute Sj=VjTRj-1, Rj=Rj-1-VjSj. (2.5) If RjF<TOL, break. End For ( 3 ) Solve the triangular system Ujtj=Dj for tj. ( 4 ) Form the approximate solution Xj=X0+[R0,V1,…,Vj-1]tj. ( 5 ) If RjF<TOL, then accept Xj and exit; otherwise, restart: set X0=Xj and go to Step (1).


## 3. Comparison with SBGMRES

In [8], an equivalence between SBGMRES and classical block GMRES had been established. Algorithm 2 indicates that RB-SBGMRES is equivalent to block GMRES; that is, search a solution Xm∈Km(A,R0), such that Rm⊥AKm(A,R0).On the other hand, for single right-hand side, it has been observed in [10] that the gap between the true residual RTm=B-AXm and the updated residual Rm can be strongly influenced by the conditioning of Zm, which is the basis of Km(A,R0), and the choice of the basis Zm has an effect on the conditioning of the matrix Um [10]. Since we compute the coordinates of the correction Xm-X0 in the basis Zm by (13), the approximate solution Xm becomes inaccurate as the conditioning of Um grows. Simpler GMRES [11] is, in general, less accurate than GMRES and is inherently unstable due to the choice of the basis [R0,V1,…,Vm-1]. It is easy to formulate an analogous conclusion in the block case.There is a theorem about condition number of[R0/R0F,V1,…,Vm-1] in [8] stated as follows.Theorem 1.
ForZm=[R0/R0F,V1,…,Vm-1], one has(14)κZm≤2s1/2κRm-1R0FRm-1F,where s is the number of right-hand sides.The condition number ofZm=[R0/R0F,…,Rm-1/Rm-1F] may have an effect on the conditioning of the matrix Um. In the following, we formulate an analogous theorem on the condition number of Zm.Theorem 2.
Suppose thatm steps of RB-SBGMRES have been taken and Riv=Ri-1v, i=1,…,m-1, if and only if v is a zero vector; then, one has(15)κZm≤msκRm-12+∑i=1m-1Ri-1F2+RiF2minv=1⁡Ri-1v2-Riv2,where Zm=[R0/R0F,…,Rm-1/Rm-1F] and s is the number of right-hand sides.Proof.
LetRm-1=Q¯m-1R¯m-1 be a QR factorization of Rm-1. We have (16)Zm=R0R0F,…,Rm-1Rm-1F=V1,…,Vm-1,Q¯m-1Cm,where(17)Cm=S1R0F0⋯0S2R0FS2R1F⋯⋮⋮⋮⋮Sm-1R0FSm-1R1F⋯0R¯m-1R0FR¯m-1R1F⋯R¯m-1Rm-1F.Moreover, we get Rm-1F=R¯m-1F.
From (11), we obtain(18)RiF2=∑j=im-1SjF2+Rm-1F2.Then, it follows that(19)CmF=m.Since the columns of [V1,…,Vm-1,Q¯m-1] are orthogonal, it follows that(20)κZm=κCm.On the other hand,(21)Cm-1=S1-1R0F00⋯00-S1-1R1FS2-1R1F0⋯000-S2-1R2F⋱⋯00⋮⋮⋱⋱⋮⋮000⋱Sm-1-1Rm-2F0000⋯-Sm-1-1Rm-1FR¯m-1-1Rm-1F.Thus, using the triangle inequality and the fact that the 1-norm is bounded from above by the 2-norm, we have(22)Cm-12≤R¯m-1-12Rm-1F2+∑i=1m-1Ri-1F2+RiF2Si-12.The norm of Si-1 can be expressed using (11) as(23)Si-12=1minv=1⁡Siv2=1minv=1⁡Ri-1v2-Riv2.With Rm-1F≤sRm-1, we have(24)Cm-12≤sκRm-12+∑i=1m-1Ri-1F2+RiF2minv=1⁡Ri-1v2-Riv2.The proof then follows from Cm2≤CmF=m.Theorem1 indicates that the conditioning of [R0/R0F,Vm-1] is inversely proportional to the actual relative norm of the residual. Once residuals become small, this will lead to the ill-conditioning of the matrices Zm and the matrices Um, and SBGMRES will behave unstably after some initial residual reduction. However, from Theorem 2, the conditioning of [R0/R0F,…,Rm-1/Rm-1F] is related to the intermediate decrease of the residual norms, not to the residual decrease with respect to the initial residual. For single right-hand sides, it has been observed that [R0/R0F,…,Rm-1/Rm-1F] remains well conditioned while [R0/R0F,Vm-1] becomes ill-conditioned [10].


## 4. RB-SBGMRES with Deflation

When block Krylov subspace methods are used for the solution of linear systems of equations with multiple right-hand sides, the linear dependence of the residual of thes systems may occur, and this is called deflation. Deflation may be possible at startup or in a later step. Sometimes, we need to incorporate a strategy for detecting when a linear combination of systems has approximately converged. Recently, Calandra et al. derived a deflation strategy to detect a near rank deficiency occurring in the block Arnoldi procedure in [7]. We provide a brief overview of the method.Assume that the QR factorization ofR0 has been performed as(25)R0=Y^1ρ^0,with Y^1∈Rn×s having orthonormal columns and ρ^0∈Rs×s. To circumvent deflation at startup, the subspace decomposition at the beginning of the cycle is derived by finding a unitary matrix F1∈Rs×s, such that(26)Y1,P0=Y^1F1,with Y1∈Rn×k1, P0∈Rn×d1, and k1+d1=s.The unitary matrixF1 is determined by the singular value decomposition (SVD) of ρ^0=UΣWT and set F1=U. Choose a relative deflation threshold ε and select k1 singular values of ρ^0 such that σ(ρ^0)>ε.Defineq0=0 and qj=qj-1+kj, for [Yj,Pj-1]∈Rn×(qj+dj) with orthonormal columns, and assume that the following block Arnoldi relation holds at the beginning of the jth iteration:(27)AYj-1=Yj,Pj-1Hj-1,j>1,with Yj-1∈Rn×qj-1, Yj∈Rn×qj, Pj-1∈Rn×dj, and Hj-1∈R(qj-1+p)×qj-1. The jth iteration of the deflated block Arnoldi procedure produces matrices Y^j+1∈Rn×kj and H^j∈R(qj+p)×qj which satisfy(28)AYj-1,Yj=Yj,Pj-1,Y^j+1H^j,with H^j having the following block structure: (29)H^j=Hj-1Hj0kj×qj-1Hj+1,j.DefiningY^j+1=[Yj,Pj-1,Y^j+1], (28) can be reformulated as(30)AYj=Y^j+1H^j-1.The subspace decomposition is performed by finding a unitary matrix Fj+1 such that(31)Yj+1,Pj=Y^j+1Fj+1.Hence, we obtain (32)AYj=Y^j+1Fj+1Fj+1TH^j-1.Defining Hj∈R(qj+p)×qj as Hj=Fj+1H^j, this leads to (33)AYj=Yj+1,PjHj,which is the block Arnoldi relation required at the beginning of the next iteration.From (31), the unitary matrix Fj+1 has the following matrix structure:(34)Fj+1=Iqj0qj×p0p×qjFj.Defining ρ^j∈R(qj+p)×p as ρ^j=FjTρ^j-10kj×p and R^j=ρ^j-H^jYj with R^0=ρ^0, the unitary matrix Fj+1 is determined by the following steps (for details, see [7]): (1)
SVD ofR^j=UΣWT, with U=Uqj+Uqj-Us+Us-.(2)
QR decomposition ofs×s matrix Us+Us-.In the case of simpler block GMRES method, the relationshipRj=Y^j+1R^j, which is an important ingredient for the block GMRES method, cannot be established. We must find another strategy to perform the subspace decomposition.Assume that the QR factorization ofR0 has been performed as (35)R0=Z^1ρ^0.We computeG0 analogously as F1 in (26), which leads to the following formulas:(36)Z1,Q0=Z^1G0,R0=Z1,Q0ρ0,with ρ0=G0Tρ^0, Z1∈Rn×c1 and Q0∈Rn×l1, c1+l1=s. To build a block orthogonal basis of AKm(A,R0), we compute the QR factorization of AZ1:(37)AZ1=V^1U^1,AZ1=V1,P0U1,with [V1,P0]=V^1F1, V^1∈Rn×c1, V1=V1, U^1∈Rc1×c1, P0∈Rn×d1, and k1+d1=c1, and F1∈Rc1×c1 is computed similarly to G0.Assume that the block Arnoldi relation(38)AZj=Vj,Pj-1Uj,holds at the beginning of the (j+1)th iteration of the deflated block Arnoldi procedure, with Zj∈Rn×cj, Vj∈Rn×kj, Pj-1∈Rn×dj, Uj∈Rqj×qj, Zj=[Z1,…,Zj], and Vj=[V1,…,Vj], where qj, sj, and kj are defined as follows: q1=c1, qj=qj-1+cj, s1=k1, and sj=sj+1=sj+kj+1.SetRj=Rj-1-VjSj, with Vj∈Rn×kj, and define Sj∈Rkj×s as Sj=VjTRj-1. We generate Zj+1 by calculating the QR factorization of Rj, which leads to the following formulas:(39)Rj=Z^j+1ρ^j=Zj+1,Qjρj,with Zj+1∈Rn×cj+1, Qj∈Rn×lj+1, cj+1+lj+1=s, [Zj+1,Qj]=Z^j+1Gj, and ρj=GjTρ^j, and Gj is a unitary matrix, which is determined similarly to G0.ComputeAZj+1 and make it orthogonal to the columns of [Vj,Pj-1] by an ordinary modified Gram-Schmidt process, such that (40)AZj,Zj+1=Vj,Pj-1,V^j+1U^j+1,with V^j+1∈Rn×cj+1, and U^j+1∈Rqj+1×qj+1 have the following block structure:(41)U^j+1=UjUj+10Uj+1,j+1=Uj+Uj+1+Uj-Uj+1-0Uj+1,j+1,where Uj+1∈R(sj+dj)×cj+1, sj+dj=qj,Uj+1,j+1∈Rcj+1×cj+1, Uj+1+∈Rsj×cj+1, and Uj+1-∈Rdj×cj+1. Define Λj+1∈R(dj+cj+1)×qj as(42)Λj+1=Uj-Uj+1-0Uj+1,j+1.We formulate(43)Vj+1,Pj=Pj-1,V^j+1Fj+1,by calculating the SVD of Λj+1=UΣWT, instead of computing the SVD of Rj. Note that U∈R(dj+cj+1)×(dj+cj+1); we set Fj+1=U and select kj+1 singular values of Λj+1 such that σ(Λj+1)>ε. Once a near rank deficiency occurs in the block Arnoldi procedure, it will be reflected by the singular values of Λj+1. We now present the RB-SBGMRES method with deflation as shown in Algorithms 4 and 5.Algorithm 4:RB-SBGMRES-D(m).
( 1 ) Given X0∈Rn×s, set R0=B-AX0. If R0F≤TOL, accept X0 and exit. ( 2 ) Compute the QR factorization of R0:R0=Z^1ρ^0. ( 3 ) Determine deflation unitary matrix G0∈Rs×s and c1,l1 such that c1+l1=s (see Algorithm 5), and set q1=c1. ( 4 ) Define [Z1,Q0]=Z^1G0, with Z1∈Rn×c1,(Q0∈Rn×l1) as the first columns of Z^1G0, and define Z1=Z1. ( 5 ) For j=1,…,m (5.1) Compute Wj=AZj. (5.2) If j=1 (5.2.1) Compute the QR factorization of Wj:Wj=AZ1=V^1U^1. (5.2.2) Determine deflation unitary matrix F1∈Rc1×c1 and k1,d1 such that k1+d1=c1 (see Algorithm 5), and set s1=k1. (5.2.3) Define [V1,P0]=V^1F1, with V1∈Rn×k1,(P0∈Rn×d1) as the first columns of V^1F1, and define V1=V1. (5.3) Else (5.3.1) For i=1,…,j-1 Ui,j=ViTWj. Wj=Wj-ViUi,j. End For (5.3.2)  Up=Pj-1TWj. (5.3.3)  Wj=Wj-Pj-1Up. (5.3.4) Define UjT=[U1,j,…,Uj-1,j,Up]T, compute the QR factorization of Wj:Wj=V^jUj,j, leading to AZj=[Vj-1,Pj-2,V^j]U^j, with U^j∈Rqj×qj defined as (41). (5.4) Determine deflation unitary matrix Fj∈R(cj+dj-1)×(cj+dj-1) by the SVD of Λj, and kj,dj such that kj+dj=cj (see Algorithm 5), and set sj=sj-1+kj. (5.5) Define [Vj,Pj-1]=[Pj-2,V^j]Fj, with Vj∈Rn×kj,(Pj-1∈Rn×dj) as the first columns of [Pj-2,V^j]Fj, and set Uj=[Isj-1×sj-1,FjT]U^j, Vj=[Vj-1,Vj]. End If (5.6) Compute Sj=VjTRj-1, Rj=Rj-1-VjSj. (5.7) If RjF<TOL, break. (5.8) Compute the QR factorization of Rj:Rj=Z^j+1ρ^j. (5.9) Determine deflation unitary matrix Gj∈Rs×s and cj+1,lj+1 such that cj+1+lj+1=s (see Algorithm 5), and set qj+1=qj+cj+1. (5.10) Define [Zj+1,Qj]=Z^j+1Gj, with Zj+1∈Rn×cj+1,(Qj∈Rn×lj+1) as the first columns of Z^j+1Gj, and define Zj+1=[Zj,Zj+1]. End For ( 6 ) Solve the triangular system Ujtj=Dj for tj, with Dj=[S1T,…,SjT,(Pj-1TR0)T]T ( 7 ) From the approximate solution Xj=X0+[Z1,…,Zj]tj. ( 8 ) If RjF<TOL, then accept Xj and exit; otherwise, restart: set X0=Xj and go to Step (1).Algorithm 5:Determination of (kj, dj, Fj) or (cj, lj, Gj).
( 1 ) Choose a relative deflation threshold ε. ( 2 ) Compute the SVD of Λj as Λj=UΣWT, with U∈R(cj+dj-1)×(cj+dj-1), Σ∈R(cj+dj-1)×(cj+dj-1), WT∈R(cj+dj-1)×qj; or the SVD of ρ^j as ρ^j=U1Σ1W1T, with U1∈Rs×s, Σ1∈Rs×s, W1∈Rs×s. ( 3 ) Select kj singular values of Λj, such that σl(Λj)>ε; or cj singular values of ρ^j such that σl(ρ^j)>ε. ( 4 ) Define Fj as Fj=U; or Gj as Gj=U1.Denote bykj,m and cj,m the number of Krylov directions keeping at the jth iteration of the mth cycle; Algorithm 4 indicates that kj+1,m<kj,m and cj+1,m<cj,m.


## 5. Numerical Experiments

In this section, RB-SBGMRES is tested and compared with SBGMRES. The test matrices were taken from the Matrix Market [9]. All computations were carried out using Matlab on a PC with the usual double precision, where the floating point relative accuracy is 2.22×10-16. In the following examples, we take X0 as the zero matrix; thus, the initial residual matrix is R0=B-AX0=B, and we set h=π/N, where N is the order of the matrix A. We plot the relative true norm of residual(44)RTmFR0F=B-AXmFB-AX0F=B-AXmFBFand condition number of Um, respectively, in Figures 1, 2, and 3. The condition number of Um is computed by Matlab internal functioncond.Figure 1
Example1, the relative true residual norms versus the number of iterations.Figure 2
Example1, the condition numbers versus the number of iterations.Figure 3
Example1, the relative true residual norms versus the number of iterations.Example 1.
The matrix is fs1836 (A≈1.2·109, κA≈1.5·1011) which is a matrix of order 183, and we have s=2 right-hand sides which are chosen asB(:,1)=A·e,B(:,2)=A·(sin⁡h,sin⁡2h,…,sin⁡(Nh))T,
whereB(:,m) is the mth column of matrix B, h=π/N, and e≡[1,1,…,1]T.In Figures1 and 2, we show the relative true norm of residual B-AXmF/BF and condition number of Um for the SBGMRES and RB-SBGMRES, respectively. It is clearly seen from Figure 1 that the relative true norm of residual of SBGMRES may stagnate at a significantly higher level than that of RB-SBGMRES. The reason for this is that the condition number of the matrix Um of SBGMRES increases significantly faster than that of RB-SBGMRES as Figure 2 shows. In Figure 3, we show the relative true norm of residual B-AXmF/BF for the RB-SBGMRES-D and the BFGMRES-S proposed in [7], with ε=10-9. It is clearly seen from the figure that RB-SBGMRES-D can compete with BFGMRES-S.Example 2.
The matrix is steam1 (A≈2.2·107, κA≈3·107) which is a matrix of order 240, and we have s=5 right-hand sides which are chosen asB:,1=A·e,B(:,2)=A·(sin⁡h,sin⁡2h,…,sin⁡(Nh))T,B(:,3)=A·(h,2h,…,Nh)T,B(:,4)=A·(cos⁡h,cos⁡2h,…,cos⁡(Nh))T,B(:,5)=A·(cos⁡2h,cos⁡2h,…,cos⁡(Nh))T,
whereB(:,m) is the mth column of matrix B, h=π/N, and N is the order of matrix A.It is clear from Figure4 that the condition number of the matrix Um of SBGMRES increases faster than that of RB-SBGMRES. Since the number of right-hand sides s=5, the maximum iteration number is 48. From Figure 5, we can observe that the numerical performance of RB-SBGMRES is better than that of SBGMRES.Figure 4
Example2, the relative true residual norms versus the number of iterations.Figure 5
Example2, the condition numbers versus the number of iterations.Figure6 shows the relative true norm of residual B-AXmF/BF for the RB-SBGMRES-D and BFGMRES-S, with ε=10-9. It is clearly seen that the performances of RB-SBGMRES-D and BFGMRES-S are almost same.Figure 6
Example2, the relative true residual norms versus the number of iterations.Example 3.
The matrix is add20 (A≈0.72, κA≈1.76·104) which is a matrix of order 2395, and we have s=2 right-hand sides which are chosen as B(:,1)=A·e,B(i,2)=sin⁡(ih),  for i=1,2,…,N.It is clear from Figure7 that the numerical performance of RB-SBGMRES is better than that of SBGMRES, and from Figure 8, we see that the condition number of the matrix Um of SBGMRES increases fast while that of RB-SBGMRES remains at a significantly low level. Figure 9 shows the relative true norm of residual B-AXmF/BF for the RB-SBGMRES-D and BFGMRES-S, with ε=10-15. The performances of RB-SBGMRES-D and BFGMRES-S are almost the same for Example 3.Figure 7
Example3, the relative true residual norms versus the number of iterations.Figure 8
Example3, the condition numbers versus the number of iterations.Figure 9
Example3, the relative true residual norms versus the number of iterations.Example 4.
The matrix is fs7601 (A≈3·108, κA≈8.4·103) which is a matrix of order 760, and we have s=3 right-hand sides which are chosen asB(:,1)=A·e,B(i,2)=sin⁡(ih),B(i,3)=ih,  for i=1,2,…,N.It is also obvious from Figure10 that the RB-SBGMRES method is slightly more accurate than SBGMRES in this example, and Figure 11 also shows that the condition number of the matrix Um of SBGMRES increases faster than that of RB-SBGMRES. It is clear from Figure 12 that the performances of RB-SBGMRES-D and BFGMRES-S are almost the same for Example 4.Figure 10
Example4, the relative true residual norms versus the number of iterations.Figure 11
Example4, the condition numbers versus the number of iterations.Figure 12
Example4, the relative true residual norms versus the number of iterations.In order to further verify that the condition number of the matrixUm of SBGMRES increases significantly faster than that of RB-SBGMRES, we compared a broad selection of different matrices from the Matrix Market and presented a comparison of the overall condition number trend. We select matrices, randomly, and set the same convergence threshold for two methods. We compare the iterations required to converge for two methods. It is easy to see from Table 3 that the condition number of the matrix Um of RB-SBGMRES is significantly smaller than that of SBGMRES and the number of iterations for the SBGMRES is slightly larger than that of RB-SBGMRES for most matrices.Table 3 Matrix RB-SBGMRES SBGMRES cond(Um) Iterations required to converge cond(Um) Iterations required to converge ck400 1.5844e+07 99 1.0848e+13 99 nos3 1.8755 e + 06 123 6.6698e+16 135 dwa512 6.4465 e + 06 128 4.5023e+13 130 dwg961b 5.1749 e + 08 238 6.7631e+17 242 dwt607 8.9814 e + 05 136 3.2494 e + 13 137 fidap001 7.1267 e + 05 53 2.3303 e + 11 53 fs5411 2.1648 e + 04 13 1.5152 e + 17 13 gre115 1.2874 e + 04 28 2.0989 e + 11 28 impcola 2.8276 e + 10 51 4.1130e+12 51 jpwh991 1.2250 e + 03 61 2.3160e+16 65 lns511 1.3654 e + 17 70 1.4647e+18 86 mbeacxc 5.8350 e + 16 112 1.1590e+18 112 olm500 2.6871 e + 17 68 1.3127e+18 68 orsirr2 5.1311 e + 07 221 6.6141 e + 18 222 pde900 3.3154 e + 03 94 8.9738 e + 15 98 pores1 2.0692 e + 08 7 3.3978e+12 7 qc324 9.5403 e + 08 81 1.6195 e + 13 92 qh768 1.5526 e + 19 30 4.3068e+17 30 rdb450 3.7541 e + 11 81 2.6783 e + 17 84 saylr1 1.6505 e + 10 59 1.0718e+17 62 sherman2 3.8156 e + 15 268 1.1791e+19 270 shl400 3.1199 e + 10 165 4.0589e+10 165 steam3 9.9777 e + 07 13 6.3224 e + 19 20 st400 3.6364 e + 16 91 5.5214e+17 93 tols340 7.7204 e + 09 48 1.9115e+20 49 utm300 7.5972 e + 11 73 1.5548 e + 17 73 curtis54 1.0376 e + 04 12 3.6682e+06 12 west0132 3.6951 e + 13 32 2.9278e+14 32 nos5 1.8655 e + 12 117 4.3780 e + 17 119 nos7 3.6414 e + 12 130 4.2740 e + 18 183


## 6. Conclusion

In this paper, we have proposed a minimum residual method mathematically equivalent to the block GMRES method for solving systems of linear equations with multiple right-hand sides. Numerical experiments show that, after some initial reduction, the relative true norms of residual SBGMRES may stagnate at a significantly higher level than that of RB-SBGMRES. This difference is clearly caused by the choice of the basisZm, which has an effect on the condition number of the matrix Um. Numerical experiments indicate that Um of RB-SBGMRES remains better-conditioned than Um of simpler block GMRES, which may become a very ill-conditioned triangular matrix. Since the coordinates of the correction Xm-X0 in the basis Zm are computed from (13), its error starts to diverge as k(Um) grows and Xm will become inaccurate. We see that the choice Zm=[R0/R0F,…,Rm-1/Rm-1F] has a better numerical performance. In comparison with the case of deflation, we consider a deflation strategy to detect the possible linear dependence of the residuals of the s systems and a near rank deficiency occurring in the block Arnoldi procedure for RB-SBGMRES method, which was later called RB-SBGMRES-D. Numerical experiments show that the performances of RB-SBGMRES-D and BFGMRES-S are almost the same.



---

*Source: 1369707-2018-04-16.xml*