A Convergent Restarted GMRES Method For Large Linear Systems Minghua Xu2,Jinzi Zhaol,Jiancheng Wu2 and Hongjun Fan 1.State Key Laboratory for Novel Software Technology,Nanjing University,Nanjing 210093,P.R.China 2.Department of Information Science,Jiangsu Polytechnic University,Changzhou,213016,P.R.China October 3,2003 Abstract The GMRES method is popular for solving nonsymmetric linear equations.It is generally used with restarting to reduce storage and orthogonal- ization costs.However,it is possible to show that the restarted GMRES method may not converge,i.e.,it may be stationary.To remedy this difficulty,a new convergent restarted GMRES method is discussed in this paper. Key words GMRES,Krylov subspace,iterative methods,nonsymmetric sys- tems. AMS(MOS)subject classifications 65F10 1.Introduction The restarted GMRES algorithm GMRES(m)1]proposed by Saad and Schultz is one of the most popular iterative methods for solving large linear systems of equations Ax=b,A∈Rnxm,x,b∈R, (1.1) with a sparse,nonsymmetric,and nonsingular matrix A.It is known that when A is positive real,the restarted GMRES method will produce a sequence of ap- proximates tk that converge to the exact solution.However,when A is not positive real,this method often slows down convergence and stagnates.The analysis and implementation of the restarted GMRES algorithm continue to re- ceive considerable attention245.6.7..For example,Y.Saad suggested a flexible inner-outer preconditioned GMRES method FGMRES(m).R.B.Morgan gave a restarted GMRES method augmented with eigenvectors ],and Cao Zhihao et. al.presented a convergent restarted GMRES algorithm based on the algorithm FGMRES(m)14l.We will now briefly review the algorithm GMRES in this sec- tion.A new restarted GMRES method and its analysis will be given in section 2,section 3 gives the examples and comparisons,and conclusions are given in section 4.The restarted GMRES can be briefly described as follows. Algorithm 1:GMRES(m)for systems (1.1)
A Convergent Restarted GMRES Method For Large Linear Systems M inghua Xu1,2 ,J inxi Zhao1 ,J iancheng W u2 and Hongjun F an1 1.State Key Laboratory for Novel Software Technology ,Nanjing University, Nanjing 210093,P.R.China 2.Department of Information Science, Jiangsu Polytechnic University,Changzhou,213016,P.R.China October 3, 2003 Abstract The GMRES method is popular for solving nonsymmetric linear equations. It is generally used with restarting to reduce storage and orthogonalization costs. However, it is possible to show that the restarted GMRES method may not converge, i.e., it may be stationary. To remedy this difficulty, a new convergent restarted GMRES method is discussed in this paper. Key words GMRES, Krylov subspace, iterative methods, nonsymmetric systems. AMS(MOS) subject classifications 65F10 1. Introduction The restarted GMRES algorithm GMRES(m)[1] proposed by Saad and Schultz is one of the most popular iterative methods for solving large linear systems of equations Ax = b, A ∈ R n×n , x, b ∈ R n , (1.1) with a sparse, nonsymmetric, and nonsingular matrix A. It is known that when A is positive real, the restarted GMRES method will produce a sequence of approximates xk that converge to the exact solution. However, when A is not positive real, this method often slows down convergence and stagnates. The analysis and implementation of the restarted GMRES algorithm continue to receive considerable attention [2,3,4,5,6,7,8]. For example, Y.Saad suggested a flexible inner-outer preconditioned GMRES method FGMRES(m)[2]. R.B.Morgan gave a restarted GMRES method augmented with eigenvectors [3] ,and Cao Zhihao et. al. presented a convergent restarted GMRES algorithm based on the algorithm FGMRES(m)[4]. We will now briefly review the algorithm GMRES in this section. A new restarted GMRES method and its analysis will be given in section 2, section 3 gives the examples and comparisons, and conclusions are given in section 4. The restarted GMRES can be briefly described as follows. Algorithm 1: GMRES(m) for systems (1.1) 1
1.Start:Choose xo and compute ro =b-Axo and B=lroll,v1 ro/B. 2.Iterate:For j=1,2,...,m do: h,i=(（Avj,v),i=1,2,…,j +1=4妈-含4 合1 hj+1,=l⑦j+1ll: Uj+1=Di+1/h+1 3.Form the approximate solution: Em=o+Vmym;where ym minimizes lBe1-Hmyll,yE Rm.Here Hm is the (m+1)by m matrix whose only nonzero entries are the elements hij defined in step 2.Vm [v1,v2,...,Um]and the vector el is the first column of the (m +1)x(m +1)identity matrix. 4.Restart: Compute rm =b-Azm,if satisfied then stop else compute zo : Im,ro :=rm,B :=lroll,v1 :ro/B and go to 2. If A is not positive real,then ro span{Aro,A2ro,.,Amro}may happen. In this situation the restarted GMRES method is stationary.To avoid this disad- vantage,we introduce and analyze a new convergent restarted GMRES method. Conveniently,we use the term CGMRES(m)to denote the method. 2.CGMRES(m) The linear systems associated with (1.1)can be taken as the following form [-a[]=[ (2.1) where I E Rnxn is the identity matrix,while u*E Rn is a given vector and f u*+b,g =-ATu*E R".Since A is nonsingular,thus the system (2.1) u* has an unique solution z*= x* Let zo is the initial approximate solution 2B-[]%=] B2o.Solving the systems (2.1)with GMRES(m),where m >2,we have the following results: Proposition 2.1 Denoting by Bi,i=1,2,...,n,the eigenvalues of ATA and supposing 3≥32≥…≥3m21/4, (2.2) then we have that the eigenvalues of matrix B have positive real part. Proof We have .- (2.3) 2
1. Start: Choose x0 and compute r0 = b − Ax0 and β = kr0k, v1 = r0/β. 2. Iterate: For j = 1, 2, · · · , m do: hi,j = (Avj , vi), i = 1, 2, · · · , j, v¯j+1 = Avj − P j i=1 hi,jvi , hj+1,j = kv¯j+1k, vj+1 = ¯vj+1/hj+1,j . 3. Form the approximate solution: xm = x0 + Vmym, where ym minimizes kβe1 − Hmyk, y ∈ Rm. Here Hm is the (m+ 1) by m matrix whose only nonzero entries are the elements hi,j defined in step 2. Vm = [v1, v2, · · · , vm] and the vector e1 is the first column of the (m + 1) × (m + 1) identity matrix. 4. Restart: Compute rm = b − Axm, if satisfied then stop else compute x0 := xm, r0 := rm, β := kr0k, v1 := r0/β and go to 2. If A is not positive real, then r0 ⊥ span{Ar0, A2 r0, · · · , Amr0} may happen. In this situation the restarted GMRES method is stationary. To avoid this disadvantage, we introduce and analyze a new convergent restarted GMRES method. Conveniently , we use the term CGMRES(m) to denote the method. 2. CGMRES(m) The linear systems associated with (1.1) can be taken as the following form " I A −AT 0 # " u ∗ x # = " f g # , (2.1) where I ∈ Rn×n is the identity matrix , while u ∗ ∈ Rn is a given vector and f = u ∗ + b, g = −AT u ∗ ∈ Rn . Since A is nonsingular, thus the system (2.1) has an unique solution z ∗ = " u ∗ x ∗ # . Let z0 is the initial approximate solution of (2.1), B = " I A −AT 0 # , r¯0 = " f g # − Bz0. Solving the systems (2.1) with GMRES(m),where m ≥ 2 , we have the following results: Proposition 2.1 Denoting by βi , i = 1, 2, · · · , n, the eigenvalues of AT A and supposing β1 ≥ β2 ≥ · · · ≥ βn ≥ 1/4, (2.2) then we have that the eigenvalues of matrix B have positive real part. Proof We have |λI1 − B| = ¯ ¯ ¯ ¯ ¯ (λ − 1)I −A AT λI ¯ ¯ ¯ ¯ ¯ , (2.3) 2
where I1 E R2nx2n is a identity matrix.If A=1 is the eigenvalue of matrix B,by (2.3),we get 0. This is a contradiction with that matrix A is nonsingular.Hence,we have A1. According to(2.3)and入≠l,we have |λI1-B|= (A-1)1 -A AT λI (入-1)1 -A 0 +(-1)-1ATA =I(A-1)λI+ATA If (ijli =1,2,...,n;j =1,2}denote the eigenvalues of B,then Ai.j,j=1,2, can be given by solving the following equation (入-1)入+62=0，i=1,2,…,n,j=1,2. (2.4) Solving (2.4)we obtain 1=1+-4g 2 (2.5) and 2=1--4g (2.6) 2 i=1,2,...,n.Using (2.2)yields the desired result. According to [1]and Proposition 2.1,we know that using the restarted GMRES method to solve(2.1)will produce a sequence of approximations which converges to the exact solution of (2.1)when Bi>1/4,i =1,2,...,n.If the conditions 2≥1/4，i=1,2,·,n do not hold,we can use Proposition 2.2 Assume that yk,k =1,2,...,m minimizes lBe1-Hkyll,yE Rk,Hk is the (k+1)x k matrix whose nonzero entries are the elements hij defined by GMRES(m)for (2.1),z=z0+Vyk is the approximate solution of (2.1),where V=[v1,v2,...,vx],vi is the Arnoldi vector generated by GMRES(m) for(2.1),i=1,2,…,k,k=1,2,·,m(m≥2).Suppose that zk= uk R2n,uk,kRn and residual Bak ,then the following results hold: (1)lFmll2 <lroll2 and zk tents to the exact solution z* x* of(2.1) (2)x*is the exact solution of(1.1). Proof of (1)The residual vector of the approximate solution zk can be written as a=[g]-a-[r]-[a]=a1m 3
where I1 ∈ R2n×2n is a identity matrix.If λ = 1 is the eigenvalue of matrix B,by (2.3),we get ¯ ¯ ¯ ¯ ¯ 0 −A AT I ¯ ¯ ¯ ¯ ¯ = 0. This is a contradiction with that matrix A is nonsingular.Hence,we have λ 6= 1. According to (2.3) and λ 6= 1,we have |λI1 − B| = ¯ ¯ ¯ ¯ ¯ (λ − 1)I −A AT λI ¯ ¯ ¯ ¯ ¯ = ¯ ¯ ¯ ¯ ¯ (λ − 1)I −A 0 λI + (λ − 1)−1AT A ¯ ¯ ¯ ¯ ¯ = |(λ − 1)λI + AT A|. If {λi,j |i = 1, 2, · · · , n; j = 1, 2} denote the eigenvalues of B ,then λi,j , j = 1, 2 , can be given by solving the following equation (λi,j − 1)λi,j + βi = 0, i = 1, 2, · · · , n, j = 1, 2. (2.4) Solving (2.4) we obtain λi,1 = 1 + √ 1 − 4βi 2 (2.5) and λi,2 = 1 − √ 1 − 4βi 2 , (2.6) i = 1, 2, · · · , n. Using (2.2) yields the desired result. According to [1] and Proposition 2.1,we know that using the restarted GMRES method to solve (2.1) will produce a sequence of approximations which converges to the exact solution of (2.1) when βi ≥ 1/4, i = 1, 2, · · · , n. If the conditions βi ≥ 1/4, i = 1, 2, · · · , n do not hold, we can use Proposition 2.2 Assume that yk, k = 1, 2, · · · , m minimizes kβe1−Hkyk, y ∈ Rk , Hk is the (k + 1) × k matrix whose nonzero entries are the elements hi,j defined by GMRES(m) for (2.1),zk = z0 + Vkyk is the approximate solution of (2.1),where Vk = [v1, v2, · · · , vk], vi is the Arnoldi vector generated by GMRES(m) for (2.1),i = 1, 2, · · · , k, k = 1, 2, · · · , m(m ≥ 2). Suppose that zk = " uk xk # ∈ R2n , uk, xk ∈ Rn and residual ¯rk = " f g # − Bzk ,then the following results hold: (1) kr¯mk2 < kr¯0k2 and zk tents to the exact solution z ∗ = " u ∗ x ∗ # of (2.1). (2) x ∗ is the exact solution of (1.1). Proof of (1) The residual vector of the approximate solution zk can be written as r¯k = " f g # − Bzk = " u ∗ + b − uk − Axk g + AT uk # = " r¯k1 r¯k2 # , k = 0, 1, · · · , m. 3
Suppose further that=o/lo2= 1,1 ≠o.We find 1,2 B01=101,1+1A01,2-2AT⑦，1=101,1之0 (2.7) and if 1.1 =0,we have =-4r.[82=4a0 (2.8) Let Km span{ro;Bro;...,Bm-17o.We have m2=mi典 [5 -B[zo+zlll2 2EKm 9 mi迎l7o-Bz2 (2.9) 2年1m min 601 -BVmyll2, y∈Rn where z=Vmy. Using BVm Vm+1Hm;we get 热301-Bm2=i典3e1-Hn2, (2.10) where el is the first column of the (m+1)x (m +1)identity matrix.By (2.9) and (2.10),we obtain Ilk=映历-B2≤- 3BB12: (2.11) B Let BDy,R1 fo-cI B. C1=B13 According to .10 and (2.7),we have C>0,RfB1=0 and lfoll2 llfo -c1B01+c1Bll2 VB+clB13>‖R2. (2.12) = By (2.11)and (2.12),we can get lIFmll2≤lRll2<oll2: (2.13) In similar way,if⑦，1=0then⑦1,2≠0 we have BB20B2元l2, IFm2=mlo-B2≤6-B度 (2.14) 2E气m 4
Suppose further that ¯v1 = ¯r0/kr¯0k2 = " v¯1,1 v¯1,2 # 6= 0. We find v¯ T 1 Bv¯1 = ¯v T 1,1 v¯1,1 + ¯v T 1,1Av¯1,2 − v¯ T 1,2A T v¯1,1 = ¯v T 1,1 v¯1,1 ≥ 0 (2.7) and if ¯v1,1 = 0,we have v¯ T 1 B 2 v¯1 = h −v¯ T 1,2A T , 0 i " Av¯1,2 0 # = −v¯ T 1,2A T Av¯1,2 ≤ 0 (2.8) Let Km = span{r¯0, Br¯0, · · · , Bm−1 r¯0}. We have kr¯mk2 = min z∈Km k " f g # − B[z0 + z]k2 = min z∈Km kr¯0 − Bzk2 = min y∈Rm kβv¯1 − BVmyk2, (2.9) where z = Vmy. Using BVm = Vm+1Hm, we get min y∈Rm kβv¯1 − BVmyk2 = min y∈Rm kβe1 − Hmyk2, (2.10) where e1 is the first column of the (m + 1) × (m + 1) identity matrix. By (2.9) and (2.10),we obtain kr¯mk2 = min z∈Km kr¯0 − Bzk2 ≤ kr¯0 − βv¯ T 1 Bv¯1 kBv¯1k 2 2 Bv¯1k2. (2.11) Let c1 = βv¯ T 1 Bv¯1 kBv¯1k 2 2 , R1 = ¯r0 − c1Bv¯1. According to ¯v1,1 6= 0 and (2.7),we have c1 > 0, RT 1 Bv¯1 = 0 and kr¯0k2 = kr¯0 − c1Bv¯1 + c1Bv¯1k2 = q kR1k 2 2 + c 2 1 kBv¯1k 2 2 > kR1k2. (2.12) By (2.11) and (2.12), we can get kr¯mk2 ≤ kR1k2 < kr¯0k2. (2.13) In similar way, if ¯v1,1 = 0 then ¯v1,2 6= 0 we have kr¯mk2 = min z∈Km kr¯0 − Bzk2 ≤ kr¯0 − βv¯ T 1 B2 v¯1 kB2v¯1k 2 2 B 2 v¯1k2, (2.14) 4
where m≥2. Let c227 B20,R2=0-cB21. B2ō1修 According to 01.20 and (2.8),we have c2<0,B2⑦1=0 and ol2=I0-c2B2⑦1+c2B22 VR23+引B2⑦3>2l2 (2.15) Thus,we can get lrm2≤R22<Iol2 (2.16) Applying (2.13)and (2.16),we know if ro0 then lmll2<oll2,the result (1) holds. Proof of(2)Since z*= 2* satisfies (2.1),we can find -ATu*g (2.11) and Ax*=(f-u*)=b (2.12) By (2.12)we have thus obtained the result (2). According to Propositions 2.1 and 2.2,using algorithm GMRES(m)to solve system(2.1),we can obtain and approximate solution of(1.1).Now the CGM- RES(m)algorithm can be briefly described as follows: Algorithm 2:CGMRES(m)for systems(1.1) 1.Start:Choose z0 and compute ro f 9 Bzo and B=lroll,v1=ro/B. 2.Iterate:For j=1,2,...,m do: h=(Buj,,i=1,2,…,j, 0+1=BU5- Shijvi, h+1=⑦j+1, Uj+1=j+1/hj+1 3.Form the approximate solution of(2.1): 2m =20+Vmym:where Vm [v1,v2,..,Um]and ym minimizes lBe1-Hmyll,yE Rm.Here Hm is the (m +1)by m matrix whose only nonzero entries are the elements hi.i defined in step 2,and el is the first column of the (m +1)x(m+1)identity matrix. 5
where m ≥ 2. Let c2 = βv¯ T 1 B2 v¯1 kB2v¯1k 2 2 , R2 = ¯r0 − c2B 2 v¯1. According to ¯v1,2 6= 0 and (2.8),we have c2 < 0, RT 2 B 2 v¯1 = 0 and kr¯0k2 = kr¯0 − c2B2 v¯1 + c2B2 v¯1k2 = q kR2k 2 2 + c 2 2 kB2v¯1k 2 2 > kR2k2. (2.15) Thus, we can get kr¯mk2 ≤ kR2k2 < kr¯0k2. (2.16) Applying (2.13) and (2.16),we know if ¯r0 6= 0 then kr¯mk2 < kr¯0k2, the result (1) holds. Proof of (2) Since z ∗ = " u ∗ x ∗ # satisfies (2.1),we can find −A T u ∗ = g (2.11) and Ax∗ = (f − u ∗ ) = b (2.12) By (2.12) we have thus obtained the result (2) . According to Propositions 2.1 and 2.2, using algorithm GMRES(m) to solve system (2.1), we can obtain and approximate solution of (1.1). Now the CGMRES(m) algorithm can be briefly described as follows: Algorithm 2: CGMRES(m) for systems (1.1) 1. Start: Choose z0 and compute r0 = " f g # − Bz0 and β = kr0k, v1 = r0/β. 2. Iterate: For j = 1, 2, · · · , m do: hi,j = (Bvj , vi), i = 1, 2, · · · , j, v¯j+1 = Bvj − P j i=1 hi,jvi , hj+1,j = kv¯j+1k, vj+1 = ¯vj+1/hj+1,j . 3. Form the approximate solution of (2.1): zm = z0 + Vmym, where Vm = [v1, v2, · · · , vm] and ym minimizes kβe1 − Hmyk, y ∈ Rm. Here Hm is the (m + 1) by m matrix whose only nonzero entries are the elements hi,j defined in step 2, and e1 is the first column of the (m + 1) × (m + 1) identity matrix. 5