第七讲 Krylov子空间选代法子空间选代法的基本思想在一个维数较低的子空间中寻找解析解的一个“最佳”近似.子空间选代法的主要过程可以分解为下面三步(1）寻找合适的子空间；(2）在该子空间中求“最佳近似解”；(3）若这个近似解满足精度要求，则停止计算否则，重新构造一个新的子空间，并返回第（2）步主要涉及到的两个关键问题(1）如果选择和更新子空间；(2）如何在给定的子空间中寻找“最佳近似解”关于第一个问题，目前较成功的解决方案就是使用Krylov子空间关于Krylov子空间选代法的相关资料[1] Y. Saad, Iterative Methods for Sparse Linear Systems, 2nd edition, 2003 [113][2] J. Demmel, Applied Numerical Linear Algebra, 1997 [30][3] A.Greenbaum, Iterative Methods for Solving Linear Systems, 1997 [59][4] G. H. Golub and C. F Van Loan, Matrix Computations, 4th, 2013 [57][5] L.N. Trefethen and D. Bau, Numerical Linear Algebra, 1997 [125]7.1Krylov子空间7.1.1Arnoldi过程与Lanczos过程设AeRnxn,rERn，则由A和r生成的Krylov子空间定义为Km(A,r) = span(r,Ar,A?r,..., Am-1r), m ≤n,通常简记为Km.设解析解在Km中的“最佳近似解”为r(m)，我们假定r,Ar,A?r,...,Am-1r线性无关，则dimKm=m.令U1,u2，,Um是Km的一组基，则Km中的任意向量a均可表示为=yiui+y2u2+...+YmUm=Vmy,其中y=[31,J2，9mJT为线性表出系数,Vm=[u1,U2Um]于是，寻找“最佳近似解”r(m)244
第七讲 Krylov 子空间迭代法 子空间迭代法的基本思想 在一个维数较低的子空间中寻找解析解的一个 “最佳” 近似. 子空间迭代法的主要过程可以分 解为下面三步: (1) 寻找合适的子空间; (2) 在该子空间中求 “最佳近似解”; (3) 若这个近似解满足精度要求, 则停止计算; 否则, 重新构造一个新的子空间, 并返回第 (2) 步. 主要涉及到的两个关键问题 (1) 如果选择和更新子空间; (2) 如何在给定的子空间中寻找 “最佳近似解”. 关于第一个问题, 目前较成功的解决方案就是使用 Krylov 子空间. 关于 Krylov 子空间迭代法的相关资料 [1] Y. Saad, Iterative Methods for Sparse Linear Systems, 2nd edition, 2003 [113] [2] J. Demmel, Applied Numerical Linear Algebra, 1997 [30] [3] A. Greenbaum, Iterative Methods for Solving Linear Systems, 1997 [59] [4] G. H. Golub and C. F. Van Loan, Matrix Computations, 4th, 2013 [57] [5] L.N. Trefethen and D. Bau, Numerical Linear Algebra, 1997 [125] 7.1 Krylov 子空间 7.1.1 Arnoldi 过程与 Lanczos 过程 设 A ∈ R n×n , r ∈ R n , 则由 A 和 r 生成的 Krylov 子空间定义为 Km(A, r) = span{r, Ar, A2 r, . . . , Am−1 r}, m ≤ n, 通常简记为 Km. 设解析解在 Km 中的 “最佳近似解” 为 x (m) . 我们假定 r, Ar, A2 r, . . . , Am−1 r 线性无关, 则 dim Km = m. 令 v1, v2, . . . , vm 是 Km 的一组基, 则 Km 中的任意向量 x 均可表示为 x = y1v1 + y2v2 + · · · + ymvm = Vmy, 其中 y = [y1, y2, . . . , ym] T 为线性表出系数, Vm = [v1, v2, . . . , vm]. 于是, 寻找 “最佳近似解” x (m) 244
7.1Krylov子空间-245.就转化为下面两个子问题：(1）寻找一组合适的基U1,U2.,m(2)求出2m)在这组基下面的线性表出系数0m=[m]TArnoldi过程首先考虑基的选取.假定r.Ar，A2r.....Am-1r线性无关，因此它们就自然地组成Km的一组基.但为了确保算法的稳定性，一般来说，我们通常希望选取一组标准正交基.这并不困难，只需对向量组{r,Ar,A?r,...，Am-1r)进行单位正交化即可.对这个过程略加修改，就就得到下面的Arnoldi过程算法7.1.Arnoldi过程(MGS)1: U1 =r/llrll22: for j = 1 to m - 1 do3:Z=Auj4:fori=ltojdo% MGS5:hij =(us,z)G.z=z-hijviend for7:8:hj+1,j=Il/29:if hj+1.j = 0 then10:breakend if11:12:Uj+1= z/hj+1.j13: end for可以证明，Arnoldi过程生成的向量u1,u2，...,m构成Km的一组标准正交基引理7.1如果Arnoldi过程不中断，则Km=span[ui,U2,...,Um](留作练习，归纳法记Vm=[u1,2,...,m],11h1,2hi,3hi,m...h2,1h2,2h2.3h2,m..h3,2h3,mh3.3...ER(m+1)xmHm+1,m =:.hm,m.hm+1,m
7.1 Krylov 子空间 · 245 · 就转化为下面两个子问题: (1) 寻找一组合适的基 v1, v2, . . . , vm; (2) 求出 x (m) 在这组基下面的线性表出系数 y (m) = h y (m) 1 , y (m) 2 , . . . , y (m) m iT . Arnoldi 过程 首先考虑基的选取. 假定 r, Ar, A2 r, . . . , Am−1 r 线性无关, 因此它们就自然地组成 Km 的一 组基. 但为了确保算法的稳定性, 一般来说, 我们通常希望选取一组标准正交基. 这并不困难, 只需 对向量组 {r, Ar, A2 r, . . . , Am−1 r} 进行单位正交化即可. 对这个过程略加修改, 就就得到下面的 Arnoldi 过程. 算法 7.1. Arnoldi 过程 (MGS) 1: v1 = r/∥r∥2 2: for j = 1 to m − 1 do 3: z = Avj 4: for i = 1 to j do % MGS 5: hi,j = (vi , z) 6: z = z − hi,jvi 7: end for 8: hj+1,j = ∥z∥2 9: if hj+1,j = 0 then 10: break 11: end if 12: vj+1 = z/hj+1,j 13: end for 可以证明, Arnoldi 过程生成的向量 v1, v2, . . . , vm 构成 Km 的一组标准正交基. 引理 7.1 如果 Arnoldi 过程不中断, 则 Km = span{v1, v2, . . . , vm}. (留作练习, 归纳法) 记 Vm = [v1, v2, . . . , vm], Hm+1,m = h1,1 h1,2 h1,3 · · · h1,m h2,1 h2,2 h2,3 · · · h2,m h3,2 h3,3 · · · h3,m . . . . . . . . . . . . hm,m hm+1,m ∈ R (m+1)×m
-246第七讲Krylov子空间选代法则由Arnoldi过程可知hj+1,juj+1=Auj -h1,j1-h2,j2 -hijuj即h1.j:j+1hj+1.jAui==Vm+1Hm+1m(:,j).hijui=Vm+1>0i=1..0所以有AVm=Vm+1Hm+1.m=VmHm+hm+1.mUm+1em,(7.1)其中Hm是由Hm+1.m的前m行组成的矩阵，即Hm=Hm+1,m(1：m,1:m),em=[0..,0,1]TeRm.由于Vm是列正交矩阵，上式两边同乘VT可得VAVm=Hm.(7.2)等式（7.1）和(7.2)是Arnoldi过程的两个重要性质.这两个性质也可以通过下图来表示HmVmVmVm+1AHm+1,m需要指出的是，如果rAr.A?r.....Am-1r线性相关，则Arnoldi过程就会提前中断.此时，我们会得到一个不变子空间。定理7.2如果Arnoldi过程在第k步时中断，即hk+1,=0,其中k<m.则有AVk=VHk,即Kk是A的一个不变子空间（留作课外自习，直接利用等式（7.1））Lanczos过程如果A是对称矩阵，则H㎡为对称三对角矩阵，此时将其记为Tm，即β1α1：β1Tm=(7.3)βm-1βm-1Qm与Arnoldi过程类似，我们有下面的性质(7.4)AVm=VmTm+βmUm+1em
· 246 · 第七讲 Krylov 子空间迭代法 则由 Arnoldi 过程可知 hj+1,jvj+1 = Avj − h1,jv1 − h2,jv2 − · · · − hj,jvj , 即 Avj = X j+1 i=1 hi,jvi = Vm+1 h1,j . . . hj+1,j 0 . . . 0 = Vm+1Hm+1,m(:, j). 所以有 AVm = Vm+1Hm+1,m = VmHm + hm+1,mvm+1e T m, (7.1) 其中 Hm 是由Hm+1,m 的前m行组成的矩阵, 即 Hm = Hm+1,m(1 : m, 1 : m), em = [0, . . . , 0, 1]T ∈ R m. 由于 Vm 是列正交矩阵, 上式两边同乘 V T m 可得 V T mAVm = Hm. (7.2) 等式 (7.1) 和 (7.2) 是 Arnoldi 过程的两个重要性质. 这两个性质也可以通过下图来表示. A Vm = Vm+1 = Hm+1,m Vm Hm + 需要指出的是, 如果 r, Ar, A2 r, . . . , Am−1 r 线性相关, 则 Arnoldi 过程就会提前中断. 此时, 我 们会得到一个不变子空间. 定理 7.2 如果 Arnoldi 过程在第 k 步时中断, 即 hk+1,k = 0, 其中 k < m. 则有 AVk = VkHk, 即 Kk 是 A 的一个不变子空间. (留作课外自习, 直接利用等式 (7.1)) Lanczos 过程 如果 A 是对称矩阵, 则 Hm 为对称三对角矩阵, 此时将其记为 Tm, 即 Tm = α1 β1 β1 . . . . . . . . . . . . βm−1 βm−1 αm . (7.3) 与 Arnoldi 过程类似, 我们有下面的性质 AVm = VmTm + βmvm+1e T m, (7.4)
7.1Krylov子空间-247.VMAVm = Tm.(7.5)考察关系式(7.4)两边的第i列可知βjuj+1= Avj -Qjuj -Bj-1uj-1,j =1,2,...这里我们令uo=0和Bo=0.根据这个三项递推公式，Arnoldi过程可简化为下面的Lanczos过程算法7.2.Lanczos过程1: Set vo = 0 and βo = 02: U1 =r/llrl23: for j = 1 to m - 1 do4:z= Auj5:Qj =(uj,2)6:z=z-ajuj-βj-iuj-17:β, = I/z/128:if βj=0 then9:breakend if10:11:Uj+1 = 2/βj12: end for可以证明，由Lanczos过程得到的向量组[u1,U2，..，Um}是单位正交的定理7.3设[u1,v2,..,Um} 是由 Lanczos 过程得到的向量组,则i=j,(i,Uj)=Oiji,j=1,2,...,m.0,ij,（留作练习）7.1.2Krylov子空间方法-般格式Krylov子空间选代法的一般过程(1)令m=1(2）定义Krylov子空间Km;(3）找出仿射空间(0）+K㎡中的“最佳近似”解；(4）如果这个近似解满足精度要求，则迭代结束；否则令m←m+1，即将Krylov子空间的维数增加一维，并返回第(3)步
7.1 Krylov 子空间 · 247 · V T mAVm = Tm. (7.5) 考察关系式 (7.4) 两边的第 j 列可知 βjvj+1 = Avj − αjvj − βj−1vj−1, j = 1, 2, . . . , 这里我们令 v0 = 0 和 β0 = 0. 根据这个三项递推公式, Arnoldi 过程可简化为下面的 Lanczos 过程. 算法 7.2. Lanczos 过程 1: Set v0 = 0 and β0 = 0 2: v1 = r/∥r∥2 3: for j = 1 to m − 1 do 4: z = Avj 5: αj = (vj , z) 6: z = z − αjvj − βj−1vj−1 7: βj = ∥z∥2 8: if βj = 0 then 9: break 10: end if 11: vj+1 = z/βj 12: end for 可以证明, 由 Lanczos 过程得到的向量组 {v1, v2, . . . , vm} 是单位正交的. 定理 7.3 设 {v1, v2, . . . , vm} 是由 Lanczos 过程得到的向量组, 则 (vi , vj ) = δij = 1, i = j, 0, i ̸= j, i, j = 1, 2, . . . , m. (留作练习) 7.1.2 Krylov 子空间方法一般格式 Krylov 子空间迭代法的一般过程 (1) 令 m = 1 (2) 定义 Krylov 子空间 Km; (3) 找出仿射空间 x (0) + Km 中的 “最佳近似” 解; (4) 如果这个近似解满足精度要求, 则迭代结束; 否则令 m ← m + 1, 即将 Krylov 子空间的 维数增加一维, 并返回第 (3) 步
-248第七讲Krylov子空间选代法在实际计算中，为了充分利用选代初值中所包含的有用信息，我们通常是在仿射空间(0)+Km中寻找“最佳近似解”算法7.3.Krylov子空间选代算法1:选取初始向量z(0)2: 计算 ro = b - Ar(0), u1 = ro/llroll23:寻找“最佳近似解"：r(1)Er(0)+Ki=z(0)+span[ui)4:ifr(1)满足精度要求then5:终止送代6: end if7: for m = 2 to n do8:调用Arnoldi或Lanczos过程计算向量Um寻找“最佳近似解":a(m)E2(0)+Km=a(0)+span[ui,..,Um)9:ifa(m）满足精度要求then10:11:终止送代12:end if13: end for子空间的选取和更新问题可以通过Krylov子空间来解决.下面需要考虑如何寻找方程组在仿射空间r(0)+Km中的“最佳近似”解r(m).首先，我们必须给出“最佳”的定义，即r(m）满足什么条件时才是“最佳”的，不同的定义会导出不同的算法。我们很自然地想到用近似解与真解之间的距离来衡量“最佳”，即使得（m）一*在某种意义下最小，如[（m）一+2达到最小.但是由于不知道，因此这种方式往往是不实用.实用的“最佳”判别方式有：(1）使得lrmll2=I16-Aaz(m)l2最小，即极小化残量rm=b-Az(m).这种方式是可行的，当A对称时，相应的算法即为MINRES方法，当A不对称时，相应的算法即为GMRES方法(2）若A对称正定，我们也可以极小化残量的能量范数lrmllLA-1，相应的算法即为CG方法.事实上，我们有IIrmllA-1 ≤ (A-1rm)* = (b - Ar(m)TA-1(b - Ar(m)= (A-Ib - 2(m)TA(A-1b - 2(m)= ( -α(m)TA(r - (m)= C- (m)|/A.因此CG方法也可以看作是极小化能量范数一r(m)lA的算法
· 248 · 第七讲 Krylov 子空间迭代法 b 在实际计算中, 为了充分利用迭代初值中所包含的有用信息，我们通常是在仿射空间 x (0)+ Km 中寻找 “最佳近似解”. 算法 7.3. Krylov 子空间迭代算法 1: 选取初始向量 x (0) 2: 计算 r0 = b − Ax(0) , v1 = r0/∥r0∥2 3: 寻找 “最佳近似解”: x (1) ∈ x (0) + K1 = x (0) + span{v1} 4: if x (1) 满足精度要求 then 5: 终止迭代 6: end if 7: for m = 2 to n do 8: 调用 Arnoldi 或 Lanczos 过程计算向量 vm 9: 寻找 “最佳近似解”: x (m) ∈ x (0) + Km = x (0) + span{v1, . . . , vm} 10: if x (m) 满足精度要求 then 11: 终止迭代 12: end if 13: end for 子空间的选取和更新问题可以通过 Krylov 子空间来解决. 下面需要考虑如何寻找方程组在 仿射空间 x (0) + Km 中的 “最佳近似” 解 x (m) . 首先, 我们必须给出 “最佳” 的定义, 即 x (m) 满足什 么条件时才是 “最佳” 的, 不同的定义会导出不同的算法. 我们很自然地想到用近似解与真解之间的距离来衡量 “最佳”, 即使得 x (m) − x∗ 在某种意义 下最小, 如 ∥x (m) − x∗∥2 达到最小. 但是由于 x∗ 不知道, 因此这种方式往往是不实用. 实用的 “最佳” 判别方式有: (1) 使得 ∥rm∥2 = ∥b − Ax(m)∥2 最小, 即极小化残量 rm = b − Ax(m) . 这种方式是可行的, 当 A 对称时, 相应的算法即为 MINRES 方法, 当 A 不对称时, 相应的算法即为 GMRES 方法; (2) 若 A 对称正定, 我们也可以极小化残量的能量范数 ∥rm∥A−1 , 相应的算法即为 CG 方法. 事 实上, 我们有 ∥rm∥A−1 ≜ r T mA −1 rm 1 2 = (b − Ax(m) ) TA −1 (b − Ax(m) ) 1 2 = (A −1 b − x (m) ) TA(A −1 b − x (m) ) 1 2 = (x∗ − x (m) ) TA(x∗ − x (m) ) 1 2 = ∥x∗ − x (m) ∥A. 因此 CG 方法也可以看作是极小化能量范数 ∥x∗ − x (m)∥A 的算法
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 |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 |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 迭代. 不断重复该过 程, 直到收敛为止