《矩阵计算》课程教学资源（讲义）第七讲 Krylov子空间迭代法

7.2GMRES方法-249.7.2GMRES方法7.2.1算法描述GMRES方法是目前求解非对称线性方程组的最常用算法之一。在该算法中，“最佳近似解”的判别方法为“使得1lrmll2=116-Az(m)l2最小",即r(m) = arg ..mink(7.6)6 - Arl|2rEr(0)+Km下面我们就根据这个最优性条件来推导出GMRES方法设送代初始向量为a(0),则对任意向量Er(0)+Km，可设=r(0)+Vmy,其中yERm.于是有r=b-Ar= b - A(r(0) - Vm3)= ro - AVmy= βu1- Vm+1Hm+1,my= Vm+1(βei-Hm+1,my),这里β=roll2.由于Vm+1列正交，所以IIrl2 = IIVm+1(Bei - Hm+1,my)l2 = IIBe1 - Hm+1,myll2于是最优性条件（7.6）就等价于g(m) = arg in IlBe1 Hm+1,ml2.(7.7)这是一个最小二乘问题。由于Hm+1,m是一个上Hessenberg矩阵，且通常m不是很大，所以我们可以用基于Givens变换的QR分解来求解.下面就是GMRES方法的基本框架算法7.4.GMRES方法基本框架1:选取初值(0),停机标准>0,以及最大选代步数IterMax2: ro = b - Ar(0), β= Iroll23: U1= ro/β4: for j = 1 to IterMax do5:w= Auj6:fori=1 to jdo% Arnoldi过程7:hi,j = (i, w)8:w=w-hi.jvi9:end for10:hj+1,j=Iw/211:if hj+1.j =0 then12:m = j, break13:end if

7.2 GMRES 方法 · 249 · 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

-250.第七讲Krylov子空间选代法14:Uj+1 = w/hj+1,j15:%相对残量relres =ri2/B16:if relres< e then%检测是否收敛17:m= j, breakend if18:19:endfor20:解最小二乘问题(7.7),得到g(m)21: 2(m) = 2(0) + Vmy(m)7.2.2具体实施细节需要解决的问题有：(1）如何计算残量r,≤b－ArG）的范数?(2）如何求解最小二乘问题(7.7)?这两个问题可以同时处理.首先采用QR分解来求解最小二乘问题.设Hm+1m的QR分解为Hm+1,m=Qm+1Rm+1,m其中Qm+1ER(m+1)x(m+1)是正交矩阵，Rm+1,mER(m+1)xm是上三角矩阵.则(7.8)Iβe1 - Hm+1,myll2 = βQm+1e1 - Rm+1,myl2 =βq1其中RmERmxm是非奇异上三角矩阵(这里假定Hm+1.m不可约).所以问题(7.7)的解为y(m) = βRmiqi(1:m),且Irmll2 = [16 - Ar(m)12 =I/βe1 - Hm+1,my(m)|2 = β [q1(m + 1)],其中q1(m+1)表示q1的第m+1个分量.Hm+1,m的QR分解的递推计算方法由于Hm+1.m是上Hessenberg矩阵，因此我们采用Givens变换hul(1)当m=1时,H21，构造Givens变换Gi使得h21R21，即H21=GIR21GiH21-0(2）假定存在G1,G2...,Gm-1,使得(Gm-1..G2G1)Hm,m-1= Rm,m-1,即Hm,m-1 = (Gm-1-.-G2G1)Rm,m-1 ≤ QmRm,m-1

· 250 · 第七讲 Krylov 子空间迭代法 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) 如何计算残量 rj ≜ b − Ax(j) 的范数? (2) 如何求解最小二乘问题 (7.7)? 这两个问题可以同时处理. 首先采用 QR 分解来求解最小二乘问题. 设 Hm+1,m 的 QR 分解为 Hm+1,m = Q T 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 , (7.8) 其中 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 T 1 R21. (2) 假定存在 G1, G2, . . . , Gm−1, 使得 (Gm−1 · · · G2G1)Hm,m−1 = Rm,m−1, 即 Hm,m−1 = (Gm−1 · · · G2G1) TRm,m−1 ≜ Q T mRm,m−1

7.2GMRES方法-251.为了书写方便，这里假定G，的维数自动扩张，以满足矩阵乘积的需要(3）考虑Hm+1.m的QR分解.易知[Hm,m-1hm其中hm=[h1m,h2m,...,hmm]THm+1,m0hm+1,m所以有hQmR.0QmhmhmmHm+1,m000hm+1,m0hm+1,m其中hm-1是Qmhm的前m-1个元素组成的向量，hmm是Qmhm的最后一个元素.构造Givens变换Gm:[Im-100E R(m+1)x(m+1)Gm=0CmSm0-8mCmJ其中im.mhm+1,mhm.mhm + hm+1,mCmSm-hm,mhm,m令Qm0Qm+1= Gm01则RmhmhRmihm.mRm+1,m.Qm+1Hm+1.m=Gm00000hm+1,m所以可得Hm+1,m的QR分解Hm+1,m=Qm+Rm+1,m凸由Hmm-1的QR分解到Hm+1,m的QR分解，我们需要(1）计算Qmhm，即将之前的m一1个Givens变换作用到Hm+1,m的最后一列的前m个元素上，所以我们需要保留所有的Givens变换；(2)残量计算:1rmll2=Iβq1(m+1)/=[BQm+1(m+1,1)l,即GmGm-1-..G2Gi(Bei)的最后一个分量的绝对值.由于在计算rm-1时就已经计算出Gm-1G2Gi(βe1),因此这里只需做一次Givens变换即可

7.2 GMRES 方法 · 251 · b 为了书写方便, 这里假定 Gi 的维数自动扩张, 以满足矩阵乘积的需要. (3) 考虑 Hm+1,m 的 QR 分解. 易知 Hm+1,m =   Hm,m−1 hm 0 hm+1,m   , 其中 hm = [h1m, h2m, . . . , hmm] T . 所以有   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 = q 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ˆm,m 0 hm+1,m      =      Rm−1 h˜m−1 0 h˜m,m 0 0      ≜ Rm+1,m. 所以可得 Hm+1,m 的 QR 分解 Hm+1,m = QT m+1Rm+1,m. b 由 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 变换即可;

-252第七讲Krylov子空间选代法(3）y(m)的计算：当相对残量满足精度要求时,需要计算y(m)=Rmlq1(1:m),而qi即为GmGm-1*.G2Gi(βe1)算法7.5.实用GMRES方法1：选取初值r(0),停机标准e>0,以及最大送代步数IterMax2: ro = b - Ar(0), β= Iroll23: if β/l/bll2 < e then停止计算，输出近似解(0)45: end if6: U1= ro/β7:=βBe1%记录q18: for j = 1 to IterMax do9:w=Auj10:fori=1tojdo%Arnoldi过程11:hi,j =(u,w)12:w=w-hijuiend for13:14:hj+1,j = Ilwll2%选代中断15:if hj+1.j = 0 then16:m = j, breakend if17:18:Uj+1= w/hj+1,519:fori=ltoj-1do%计算Gj-1**.G2GiHi+1,j(1:j,j)hij20:hi+1,3hi+13-SiCiend for21:if [hjil >|hj+1j/ then% 构造 Givens 变换 Gj22:23:T=hij+1,j/hj,cj=1/V1+T2,sj=cjT24:else25:T=hjj/hj+1,j,8j=1/V1+T2,cj=8jTend if26:%计算 G,Hi+1,(1:j,5)27:hj = cjhjj + 8jhi+1,j28:hj+1,j = 0% 计算 G,(BGj-1...G2Gie1)29:Ej+1S2Ci%相对残量30:relres = |5j+1l/β31:if relres<e then

· 252 · 第七讲 Krylov 子空间迭代法 (3) y (m) 的计算: 当相对残量满足精度要求时, 需要计算 y (m) = R−1 m q1(1 : m), 而 q1 即为 GmGm−1 · · · G2G1(βe1). 算法 7.5. 实用 GMRES 方法 1: 选取初值 x (0) , 停机标准 ε > 0, 以及最大迭代步数 IterMax 2: r0 = b − Ax(0) , β = ∥r0∥2 3: if β/∥b∥2 < ε 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|/β % 相对残量 31: if relres < ε then

7.2GMRES方法-253.32:m=j,breakendif33:34:end for35:m=j36: y(m) = H(1 : m,1 : m))s(1 : m)%求最小二乘问题的解，回代求解37: r(m) = r(0) + Vmy(m)38: if relres < e then39:输出近似解及相关信息40: else41:输出算法失败信息42: end if7.2.3GMRES方法的中断在上面的GMRES方法中，当执行到某一步时有hk+1,=0,则算法会中断(breakdown).如果出现这种中断，则我们就找到了精确解定理7.4设 A E Rnxn 非奇异且 ro≠0.若 hi+1, ≠0,i=1,2,..,k-1,则 hk+1,k=0当且仅当r(K)是方程组的精确解.(不考虑舍入误差)(板书)证明.设hk+1,=0,则有AV=-VHk,y() -H-1(Bei),所以[rk/2 = [6 - Ar(k)l2 = [16 - A((0) + Vky(k)]l2= ro-VHky()l2= IlBu1-V(Be1)l2= 0反之，设()是精确解，则0 = b - Ar(k) = ro - Vk+1Hk+1,k3(k) = V+1(Be1 - Hk+1,ky(k)反证法，假设hk+1,≠0,则vk+1≠0.因此Vs+1单位列正交,故列满秩，所以由上式可知Be1 - Hk+1,ky(k) = 0.由于Hk+1,是上Hessenberg矩阵，且hi+1,≠0,i=1,2..,k.通过向后回代求解可得y(k)=0,口于是β=0.这与r0≠0矛盾.所以hk+1,=07.2.4带重启的GMRES方法由于随着选代步数的增加，GMRES方法的每一步所需的运算量和存储量都会越来越大.因此当选代步数很大时，GMRES方法就不太实用通常的解决方法就是重启，即事先设定一个重启选代步数k，如k=20或50等等，当GMRES达到这个选代步数时仍不收敛，则计算出方程组在(0)+K中的最佳近似解(k),然后令z(0)=2(k),并重新开始新的GMRES选代.不断重复该过程，直到收敛为止

7.2 GMRES 方法 · 253 · 32: m = j, break 33: end if 34: end for 35: m = j 36: y (m) = H(1 : m, 1 : m)\ξ(1 : m) % 求最小二乘问题的解, 回代求解 37: x (m) = x (0) + Vmy (m) 38: if relres < ε then 39: 输出近似解 x 及相关信息 40: else 41: 输出算法失败信息 42: end if 7.2.3 GMRES 方法的中断 在上面的 GMRES 方法中, 当执行到某一步时有 hk+1,k = 0, 则算法会中断 (breakdown). 如果 出现这种中断, 则我们就找到了精确解. 定理 7.4 设 A ∈ R n×n 非奇异且 r0 ̸= 0. 若 hi+1,i ̸= 0, i = 1, 2, . . . , k − 1, 则 hk+1,k = 0 当且仅 当 x (k) 是方程组的精确解. (不考虑舍入误差) (板书) 证明. 设 hk+1,k = 0, 则有 AVk = VkHk, y(k) = H −1 k (βe1). 所以 ∥rk∥2 = ∥b − Ax(k) ∥2 = ∥b − A(x (0) + Vky (k) )∥2 = ∥r0 − VkHky (k) ∥2 = ∥βv1 − Vk(βe1)∥2 = 0. 反之, 设 x (k) 是精确解, 则 0 = b − Ax(k) = r0 − Vk+1Hk+1,ky (k) = Vk+1(βe1 − Hk+1,ky (k) ). 反证法, 假设 hk+1,k ̸= 0, 则 vk+1 ̸= 0. 因此 Vk+1 单位列正交, 故列满秩, 所以由上式可知 βe1 − Hk+1,ky (k) = 0. 由于 Hk+1,k 是上 Hessenberg 矩阵, 且 hi+1,i ̸= 0, i = 1, 2, . . . , k. 通过向后回代求解可得 y (k) = 0, 于是 β = 0. 这与 r0 ̸= 0 矛盾. 所以 hk+1,k = 0 □ 7.2.4 带重启的 GMRES 方法 由于随着迭代步数的增加, GMRES 方法的每一步所需的运算量和存储量都会越来越大. 因 此当迭代步数很大时, GMRES 方法就不太实用. 通常的解决方法就是重启, 即事先设定一个重启 迭代步数 k, 如 k = 20 或 50 等等, 当 GMRES 达到这个迭代步数时仍不收敛, 则计算出方程组在 x (0) + Kk 中的最佳近似解 x (k) , 然后令 x (0) = x (k) , 并重新开始新的 GMRES 迭代. 不断重复该过 程, 直到收敛为止