.216.第七讲子空间选代方法7.2GMRES方法7.2.1算法描述GMRES方法是目前求解非对称线性方程组的最常用算法之一.在该算法中，“最佳近似”解的判别方法为"使得rml2=6-A（m)12最小"即(m) = arg min r6-Arll2(7.6)rEr(0)+K下面我们就根据这个最优性条件来导出GMRES方法设选代初始向量为r(0),则对任意向量E(0)+Km，可设=(0）+Vmy,其中yERm.于是有r=b- Ar= b - A(r(0) - Vmy)=ro-AVmy=Bu1-Vm+1Hm+1.my= Vm+1(βei-Hm+1my),这里β=llroll2由于Vm+1列正交，所以Irll2=IIVm+1(Be1-Hm+1.my)ll2=llBe1-Hm+1,myll2于是最优性条件（7.6）就等价于y(m) = arg min I/Be1 - Hm+1,m/ll2(7.7)VE这是一个最小二乘问题由于Hm+1.m是一个上Hessenberg矩阵，且通常m不是很大，所以我们可以用基于Givens变换的QR分解来求解.下面就是GMRES方法的基本框架算法7.4.GMRES方法基本框架1:选取初值r(0),停机标准ε>0,以及最大选代步数IterMax2: ro = b- Ar(0), β = [roll23: U1=To/B4: for j =1 to IterMax do5:w= Auj6:fori=lto j do%Arnoldi过程7:hij = (ui,w)8:w=w-hijviend for9:10:hij+1.,j = [wll211:if hj+1j= 0 then12:m=j,breakend if13:
· 216 · 第七讲 子空间迭代方法 7.2 GMRES 方法 7.2.1 算法描述 GMRES 方法是目前求解非对称线性方程组的最常用算法之一. 在该算法中, “最佳近似” 解的判别 方法为 “使得 ∥rm∥2 = ∥b − Ax(m)∥2 最小”, 即 x (m) = arg min x∈x(0)+Km ∥b − Ax∥2. (7.6) 下面我们就根据这个最优性条件来导出 GMRES 方法. 设迭代初始向量为 x (0) , 则对任意向量 x ∈ x (0) + Km, 可设 x = x (0) + Vmy, 其中 y ∈ R m. 于是有 r = b − Ax = b − A(x (0) − Vmy) = r0 − AVmy = βv1 − Vm+1Hm+1,my = Vm+1 (βe1 − Hm+1,my), 这里 β = ∥r0∥2. 由于 Vm+1 列正交, 所以 ∥r∥2 = ∥Vm+1(βe1 − Hm+1,my)∥2 = ∥βe1 − Hm+1,my∥2. 于是最优性条件 (7.6) 就等价于 y (m) = arg min y∈Rm ∥βe1 − Hm+1,my∥2. (7.7) 这是一个最小二乘问题. 由于 Hm+1,m 是一个上 Hessenberg 矩阵, 且通常 m 不是很大, 所以我们可以用 基于 Givens 变换的 QR 分解来求解. 下面就是 GMRES 方法的基本框架. 算法 7.4. GMRES 方法基本框架 1: 选取初值 x (0) , 停机标准 ε > 0, 以及最大迭代步数 IterMax 2: r0 = b − Ax(0) , β = ∥r0∥2 3: v1 = r0/β 4: for j = 1 to IterMax do 5: w = Avj 6: for i = 1 to j do % Arnoldi 过程 7: hi,j = (vi , w) 8: w = w − hi,jvi 9: end for 10: hj+1,j = ∥w∥2 11: if hj+1,j = 0 then 12: m = j, break 13: end if
7.2GMRES方法.217.14:Uj+1 = w/hj+1.j15:%相对残量relres = r;ll2/β16:if relres< e then%检测是否收敛17:m= j, breakendif18:19:end for20:解最小二乘问题(7.7),得到g/(m)21: 2(m) = r(0) + Vmg/(m)7.2.2具体实施细节需要解决的问题有：(1)如何计算残量rm会b-Ar（m）的范数?(2）如何求解最小二乘问题(7.7)?这两个问题可以同时处理.首先采用QR分解来求解最小二乘问题.设Hm+1m的QR分解为Hm+1m=Qm+1Rm+1,m,其中Qm+1ER(m+1)x(m+1)是正交矩阵，Rm+1,mER(m+1)xm是上三角矩阵.则IBe1 - Hm+1,myll2 = IIBQm+1e1 - Rm+1,m3/l2 =Bq1其中RmERmxm是非奇异上三角矩阵(这里假定Hm+1,m不可约).所以问题(7.7)的解为y(m) = βRm'q(1:m),且I/rm2 = b - Ar(m)2 = I|Be1 - Hm+1,my(m)l/2 = β - [q1(m + 1)/其中gi(m+1)表示gi的第m+1个分量Hm+1,m的QR分解的递推计算方法由于Hm+1.m是上Hessenberg矩阵，因此我们采用Givens变换[h11，构造Givens变换Gi使得(1)当m=1时,H21 :*=R21,即H21=GIR21:GiH210(2)假定存在G1,G2...Gm-1，使得(Gm-1...G2G1)Hm.m-1=Rm.m-1,即Hmm-1=(Gm-1...G2G)TRm,m-1≤QmRm.m-1
7.2 GMRES 方法 · 217 · 14: vj+1 = w/hj+1,j 15: relres = ∥rj∥2/β % 相对残量 16: if relres< ε then % 检测是否收敛 17: m = j, break 18: end if 19: end for 20: 解最小二乘问题 (7.7), 得到 y (m) 21: x (m) = x (0) + Vmy (m) 7.2.2 具体实施细节 需要解决的问题有: (1) 如何计算残量 rm ≜ b − Ax(m) 的范数? (2) 如何求解最小二乘问题 (7.7)? 这两个问题可以同时处理. 首先采用 QR 分解来求解最小二乘问题. 设 Hm+1,m 的 QR 分解为 Hm+1,m = Q ⊺ m+1Rm+1,m, 其中 Qm+1 ∈ R (m+1)×(m+1) 是正交矩阵, Rm+1,m ∈ R (m+1)×m 是上三角矩阵. 则 ∥βe1 − Hm+1,my∥2 = ∥βQm+1e1 − Rm+1,my∥2 = βq1 − [ Rm 0 ] y 2 , 其中 Rm ∈ R m×m 是非奇异上三角矩阵 (这里假定 Hm+1,m 不可约). 所以问题 (7.7) 的解为 y (m) = βR−1 m q1(1:m), 且 ∥rm∥2 = ∥b − Ax(m) ∥2 = ∥βe1 − Hm+1,my (m) ∥2 = β · |q1(m + 1)|, 其中 q1(m + 1) 表示 q1 的第 m + 1 个分量. Hm+1,m 的 QR 分解的递推计算方法 由于 Hm+1,m 是上 Hessenberg 矩阵, 因此我们采用 Givens 变换. (1) 当 m = 1 时, H21 = [ h11 h21] , 构造 Givens 变换 G1 使得 G1H21 = [ ∗ 0 ] = R21, 即 H21 = G ⊺ 1R21. (2) 假定存在 G1, G2, . . . , Gm−1, 使得 (Gm−1 · · · G2G1)Hm,m−1 = Rm,m−1, 即 Hm,m−1 = (Gm−1 · · · G2G1) ⊺Rm,m−1 ≜ Q ⊺ mRm,m−1
第七讲子空间选代方法.218.十为了书写方便，这里假定G，的维数自动扩张，以满足矩阵乘积的需要(3）考虑Hm+1.m的QR分解.易知[Hm,m-1hm其中hm=[him,h2m...,hmm]THm+1,m0hm+1,m所以有Rm-0RmQmQmhmm-1h.H.0241.0hm+1m00hm+1,m其中hm-1是Qmhm的前m一1个元素组成的向量，hmm是Qmhm的最后一个元素.构造Givens变换Gm：00Im-1E R(m+1)x(m+1)Gm=0CmS10-smCm其中hm.mhm+1.mhm.m/hm,m+hm+1.mhmmhm,m令[Qm0Qm+1= Gm0-则hm-1Rm-1Rm-1hmhjjhj.j00≤Rm+1,mQm+1Hm+1.m=Gm000hm+1,m所以可得Hm+1.m的QR分解Hm+1,m=Qm+1Rm+1m+由 Hm,m-1 的 QR 分解到 Hm+1,m的 QR 分解,我们需要(1）计算Qmhm，即将之前的m=1个Givens变换作用到Hm+1m的最后一列的前m个元素上，所以我们需要保留所有的Givens变换；(2)残量计算:llrmll2=[Bq1(m+1)/=[BQm+1(m+1,1)/,即GmGm-1...G2G1(Be1)的最后一个分量的绝对值.由于在计算rm-1时就已经计算出Gm-1.：G2G1(Be1),因此这里只需做一次Givens变换即可；(3)y(m）的计算：当相对残量满足精度要求时,需要计算 y(m）=Rmqi(1：m),而 q1即为GmGm-1...G2Gi(βe1)
· 218 · 第七讲 子空间迭代方法 † 为了书写方便, 这里假定 Gi 的维数自动扩张, 以满足矩阵乘积的需要. (3) 考虑 Hm+1,m 的 QR 分解. 易知 Hm+1,m = [ Hm,m−1 hm 0 hm+1,m] , 其中 hm = [h1m, h2m, . . . , hmm] ⊺ . 所以有 [ Qm 0 0 1] Hm+1,m = [ Rm,m−1 Qmhm 0 hm+1,m] = Rm−1 h˜m−1 0 hˆmm 0 hm+1,m , 其中 h˜m−1 是 Qmhm 的前 m−1 个元素组成的向量, hˆmm 是 Qmhm 的最后一个元素. 构造 Givens 变换 Gm: Gm = Im−1 0 0 0 cm sm 0 −sm cm ∈ R (m+1)×(m+1) , 其中 cm = hˆm,m h˜m,m , sm = hm+1,m h˜m,m , h˜m,m = √ hˆ2 m,m + h 2 m+1,m. 令 Qm+1 = Gm [ Qm 0 0 1] , 则 Qm+1Hm+1,m = Gm Rm−1 h˜m−1 0 hˆ j,j 0 hm+1,m = Rm−1 h˜m−1 0 h˜ j,j 0 0 ≜ Rm+1,m. 所以可得 Hm+1,m 的 QR 分解 Hm+1,m = Q ⊺ m+1Rm+1,m. † 由 Hm,m−1 的 QR 分解到 Hm+1,m 的 QR 分解, 我们需要 (1) 计算 Qmhm, 即将之前的 m − 1 个 Givens 变换作用到 Hm+1,m 的最后一列的前 m 个元素 上, 所以我们需要保留所有的 Givens 变换; (2) 残量计算: ∥rm∥2 = |βq1(m + 1)| = |βQm+1(m + 1, 1)|, 即 GmGm−1 · · · G2G1(βe1) 的最后一个分量的绝对值. 由于在计算 rm−1 时就已经计算出 Gm−1 · · · G2G1(βe1), 因此 这里只需做一次 Givens 变换即可; (3) y (m) 的计算: 当相对残量满足精度要求时, 需要计算 y (m) = R−1 m q1(1 : m), 而 q1 即为 GmGm−1 · · · G2G1(βe1)
7.2GMRES方法.219.算法7.5.实用GMRES方法1：选取初值(0)，停机标准>0,以及最大选代步数IterMax2: ro = b - Ar(o), β=Ilroll23: if β<Ethen停止计算，输出近似解(0)4中5: endif6: U1 = ro/β7:=e1%记录qi8:for j =1 to IterMax do9:w=Auj10:for i=ltojdo%Arnoldi过程11:hi.j = (vi,w)12:w=w-hi.juiend for13:hj+1,j = I/wl214:%选代中断if hij+1j=0 then15:m=j,break16:endif17:18:Uj+1=w/hij+1.jfori=1 to j-1do%计算 Gj-1-..G2GiHj+1,;(1:j,3)19:hijhijCSi20:hi+1j]hi+1.-Siendfor21:if|hjil>[hi+1.jlthen%构造Givens变换G22:23:T=hj+1,j/hj,Cj= 1/V1 +T2,8j = Cj724:elseT= hjs/hj+1j,sj =1/V1 +T2,cj =sjt25:endif26:27:%计算G,Hi+15(1:j.i)hjj =c,hjj + sjhj+1.328:hj+1j= 0SiECSi29:%计算G,(BGj-1-.G2Giei)10Sj+1-SiCi30:relres =I/j+1ll2/β%相对残量if relres< e then31:m=j,break32:endif33:34:endfor35:m=j
7.2 GMRES 方法 · 219 · 算法 7.5. 实用 GMRES 方法 1: 选取初值 x (0) , 停机标准 ε > 0, 以及最大迭代步数 IterMax 2: r0 = b − Ax(0) , β = ∥r0∥2 3: if β < ε then 4: 停止计算, 输出近似解 x (0) 5: end if 6: v1 = r0/β 7: ξ = βe1 % 记录 q1 8: for j = 1 to IterMax do 9: w = Avj 10: for i = 1 to j do % Arnoldi 过程 11: hi,j = (vi , w) 12: w = w − hi,jvi 13: end for 14: hj+1,j = ∥w∥2 15: if hj+1,j = 0 then % 迭代中断 16: m = j, break 17: end if 18: vj+1 = w/hj+1,j 19: for i = 1 to j − 1 do % 计算 Gj−1 · · · G2G1Hj+1,j (1: j, j) 20: [ hij hi+1,j] = [ ci si −si ci ] [ hij hi+1,j] 21: end for 22: if |hjj | > |hj+1,j | then % 构造 Givens 变换 Gj 23: τ = hj+1,j/hjj , cj = 1/√ 1 + τ 2, sj = cj τ 24: else 25: τ = hjj/hj+1,j , sj = 1/√ 1 + τ 2, cj = sj τ 26: end if 27: hjj = cjhjj + sjhj+1,j % 计算 GjHj+1,j (1: j, j) 28: hj+1,j = 0 29: [ ξj ξj+1] = [ cj sj −sj cj ] [ξj 0 ] % 计算 Gj (βGj−1 · · · G2G1e1) 30: relres = ∥ξj+1∥2/β % 相对残量 31: if relres< ε then 32: m = j, break 33: end if 34: end for 35: m = j