第七讲子空间迭代方法子空间迭代算法的基本思想是在一个维数较低的子空间中寻找解析解的一个“最佳”近似，子空间选代方法的主要过程可以分解为下面三步：(1）寻找合适的子空间；(2）在该子空间中求“最佳近似”解；(3）若这个近似解满足精度要求，则停止计算；否则，重新构造一个新的子空间，并返回第(2)步这里主要涉及到的两个关键问题是：（1）如果选择和更新子空间(2）如何在给定的子空间中寻找“最佳近似”解关于第一个问题，目前较成功的解决方案就是使用Krylov子空间7.1Krylov子空间7.1.1Arnoldi过程与Lanczos过程设AERnxn，rER"，则由A和r生成的Krylov子空间为Km(A,r) =span[r,Ar,A?r,...,Am-r),m≤n.通常简记为Km.设解析解在Km中的“最佳近似”为z（m).我们假定r,Ar,A?r，..,Am-1r线性无关，则dimKmm.令u1,U2....,Um是Km的一组基，则Km中的任意向量a均可表示为=YiUi+y2U2+...+YmUm=Vmy,其中y=[31,92,….，9mJT为线性表出系数,Vm=[31,U2,….,Um].于是,寻找“最佳近似"(m)就转化为（1）寻找一组合适的基1,V2,...，Um；(2)求出r(m）在这组基下面的线性表出系数y(m)Arnoldi过程首先考虑基的选取：假定r,Ar,A?r,...,Am-1r线性无关，因此它们就自然地组成Km的一组基但为了确保算法的稳定性，一般来说，我们通常希望选取一组标准正交基，这并不困难，只需对向量组{r，Ar，A?r，...，Am-1r进行单位正交化即可.对这个过程略加修改，就就得到下面的Arnoldi过程211
第七讲 子空间迭代方法 子空间迭代算法的基本思想是在一个维数较低的子空间中寻找解析解的一个 “最佳” 近似. 子空间 迭代方法的主要过程可以分解为下面三步: (1) 寻找合适的子空间; (2) 在该子空间中求 “最佳近似” 解; (3) 若这个近似解满足精度要求, 则停止计算; 否则, 重新构造一个新的子空间, 并返回第 (2) 步. 这里主要涉及到的两个关键问题是: (1) 如果选择和更新子空间; (2) 如何在给定的子空间中寻找 “最佳近似” 解. 关于第一个问题, 目前较成功的解决方案就是使用 Krylov 子空间. 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] ⊺ 为线性表出系数, Vm = [v1, v2, . . . , vm]. 于是, 寻找 “最佳近似” x (m) 就转化为 (1) 寻找一组合适的基 v1, v2, . . . , vm; (2) 求出 x (m) 在这组基下面的线性表出系数 y (m) . Arnoldi 过程 首先考虑基的选取. 假定 r, Ar, A2 r, . . . , Am−1 r 线性无关, 因此它们就自然地组成 Km 的一组基. 但为了确保算法的稳定性, 一般来说, 我们通常希望选取一组标准正交基. 这并不困难, 只需对向量组 {r, Ar, A2 r, . . . , Am−1 r} 进行单位正交化即可. 对这个过程略加修改, 就就得到下面的 Arnoldi 过程. 211
第七讲子空间选代方法.212.算法7.1.Arnoldi过程(MGS)1: U = r/llrll22: for j=1 to m -1do3:z=Auj4:fori=lto jdo% MGS5:hij = (vi,2)6:z=z-hijviend for7:8:hi+1,3 = /zll29:if hj+1.j= 0 thenbreak10:end if11:12:Uj+1 = 2/hj+1.j13: end for可以证明，Arnoldi过程生成的向量u1，2，..，Um构成Km的一组标准正交基.引理7.1如果Arnoldi过程不中断，则Km=span[u1,U2,...,Um]证明只需证明vkEKm即可.事实上，我们可以通过归纳法证明viEK，k=1,2，..，m.具体证明过口程留作练习.记Vm=[v1,02,...,Um]h1,1h1,2h1,3h1,mh2,mh2,1h2,2h2,3h3,2h3,3h3.mER(m+1)xmHm+1,m=:.hm,mhm,m-1hm+1,m则由Arnoldi过程可知hi+1,juj+1=Auj-h1,ju1-h2ju2-..-hjjuj
· 212 · 第七讲 子空间迭代方法 算法 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}. 证明. 只需证明 vk ∈ Km 即可. 事实上, 我们可以通过归纳法证明 vk ∈ Kk, k = 1, 2, . . . , m. 具体证明过 程留作练习. □ 记 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−1 hm,m hm+1,m ∈ R (m+1)×m, 则由 Arnoldi 过程可知 hj+1,jvj+1 = Avj − h1,jv1 − h2,jv2 − · · · − hj,jvj
7.1Krylov子空间.213-即h1.j:3+1hj+1.jAvi ==Vm+1Hm+1,m(:,j),hi.jui=Vm+110i=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.由于V是列正交矩阵，上式两边同乘VT可得(7.2)VIAVm=Hm等式(7.1)和(7.2)是Arnoldi过程的两个重要性质.这两个性质也可以通过下图来表示.HmVmVm+1VmAHm+1,m十由于Hm+1,m是上Hessenberg矩阵，因此Arnoldi过程也称为部分上Hessenberg化过程需要指出的是，如果r,Ar,A2r,..,Am-1r线性相关，则Arnoldi过程就会提前中断。此时，我们会得到一个不变子空间定理7.1如果Arnoldi过程在第k步时中断，即hk+1.k=0,其中k<m.则有AV=ViHk，即K是A的一个不变子空间口证明.直接由等式（7.1)可知结论成立Lanczos过程如果A是对称矩阵，则Hm为对称三对角矩阵，此时将其记为Tm，即β1[a1B1Tm =(7.3)βm-.βm-1m与Arnoldi过程类似，我们有下面的性质AVm=VmTm+BmUm+1em(7.4)
7.1 Krylov 子空间 · 213 · 即 Avj = ∑ 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 ⊺ m, (7.1) 其中 Hm 是由 Hm+1,m 的前 m 行组成的矩阵, 即 Hm = Hm+1,m(1 : m, 1 : m), em = [0, . . . , 0, 1]⊺ ∈ R m. 由于 Vm 是列正交矩阵, 上式两边同乘 V ⊺ m 可得 V ⊺ mAVm = Hm. (7.2) 等式 (7.1) 和 (7.2) 是 Arnoldi 过程的两个重要性质. 这两个性质也可以通过下图来表示. A Vm = Vm+1 = Hm+1,m Vm Hm + † 由于 Hm+1,m 是上 Hessenberg 矩阵, 因此 Arnoldi 过程也称为部分上 Hessenberg 化过程. 需要指出的是, 如果 r, Ar, A2 r, . . . , Am−1 r 线性相关, 则 Arnoldi 过程就会提前中断. 此时, 我们会 得到一个不变子空间. 定理 7.1 如果 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 ⊺ m, (7.4)
-214.第七讲子空间选代方法(7.5)VIAVm=Tm考察关系式（7.4)两边的第i列可知β,uj+1= Auj-Qjuj -βj-1uj-1,j =1,2,.这里我们令vo=0和Bo=0.根据这个三项递推公式，Arnoldi过程可简化为下面的Lanczos过程算法7.2.Lanczos过程1:Set Vo= 0 and Bo= 02: U1 = r/rll23:forj=1tom-1do4:Z= Auj5:αj =(uj,2)6:z = z - αjVj - βj-1Uj-17:β, = 212if β, = 0 then8:9:breakend if10:11:Uj+1= z/Bj12:end for可以证明，由Lanczos过程得到的向量组{u1,U2...Um】是单位正交的定理7.2设{u1,v2....,Um］是由Lanczos过程得到的向量组，则i=j,1.(ui,Ug) = ji.j=1,2....,m.0，ij口证明.留作练习7.1.2Krylov子空间算法一般格式Krylov子空间送代算法的一般过程为：(1) 令m= 1(2）定义Krylov子空间Km(3）找出仿射空间(0）+Km中的“最佳近似”解；(4）如果这个近似解满足精度要求，则选代结束；否则令m←m+1,即将Krylov子空间的维数增加一维，并返回第(3)步
· 214 · 第七讲 子空间迭代方法 V ⊺ 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.2 设 {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) 步
7.1Krylov子空间.215.算法7.3.Krylov子空间选代算法1：选取初始向量(0)2: 计算 ro =b- Ar(o),u=To/lroll23:寻找“最佳近似"解：z(1)Ea(0)+Ki=r(0)+span[ui)4:ifr(1）满足精度要求then终止送代5:6: end if7: for m = 2 to n do调用Arnoldi或Lanczos过程计算向量um8:寻找“最佳近似"解：(m)E(0)+Km=(0)+span[u1,..,Um)9:ifr(m）满足精度要求then10:11:终止选代endif12:13:endfor子空间的选取和更新问题可以通过Krylov子空间来解决下面需要考虑如何寻找方程组在仿射空间(0）+Km中的“最佳近似”解（m).首先，我们必须给出“最佳”的定义，即（m）满足什么条件时才是“最佳”的、不同的定义会导出不同的算法，我们很自然地想到用近似解与真解之间的距离来衡量“最佳，即使得（m）一工，在某种意义下最小，如（m)工2达到最小.但是由于不知道，因此这种方式往往是不实用实用的“最佳”方式有：(1)）使得Irm2=16-Az(m)l2最小,即极小化残量rm=b-Az(m).这种方式是可行的，当A对称时，相应的算法即为MINRES方法，当A不对称时，相应的算法即为GMRES方法(2）若A对称正定，我们也可以极小化残量的能量范数rmA-1,相应的算法即为CG方法.事实上，我们有IIrmllA-1 α (rMA-1rm) = (b - Ar(m)TA-1(b - Ar(m)= (A-Ib- r(m)TA(A-Ib - 2(m)= (z, - 2(m)TA(r, -2(m)= I* - 2(m)llA.因此CG方法也可以看作是极小化z(m）一工。的能量范数。一(m)A的算法
7.1 Krylov 子空间 · 215 · 算法 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 ⊺ mA −1 rm ) 1 2 = ( (b − Ax(m) ) ⊺A −1 (b − Ax(m) ) ) 1 2 = ( (A −1 b − x (m) ) ⊺A(A −1 b − x (m) ) ) 1 2 = ( (x∗ − x (m) ) ⊺A(x∗ − x (m) ) ) 1 2 = ∥x∗ − x (m) ∥A. 因此 CG 方法也可以看作是极小化 x (m) − x∗ 的能量范数 ∥x∗ − x (m)∥A 的算法
.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 β[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 β |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
第七讲子空间选代方法.220.36: y(m) = H(1 : m,1 : m)/s(1 : m)%求最小二乘问题的解，回代求解37: r(m) = r(0) + Vmy(m)38:if relres<then39:输出近似解工及相关信息40:else41输出算法失败信息42:end if7.2.3GMRES方法的中断在上面的GMRES方法中，当执行到某一步时有hj+1,j=0,则算法会中断(breakdown)。如果出现这种中断，则我们就可以找到精确解定理7.3设AeRnxn非奇异且ro≠0.若 hi+1,i≠0,i=1,2,...,k-1,则 hk+1,k=0当且仅当(K)是方程的精确解.(不考虑舍入误差)证明.必要性.设hk+1,=0,则有AV=V,Hk, y()=H-(Bei)所以[rk/2 = [b - Ar(k)12= 1b - A((0) + Vky(k)2= [Iro - VHky()ll2 = Bu1- V(Bei)ll2 = 0充分性.设(k)是精确解，则0 = b - Ar() = ro - Vk+1Hk+1,ky(k) = V+1(Be1 - Hk+1,ky(k))反证法，假设hk+1.≠0,则Uk+1≠0.因此Vk+1单位列正交，故列满秩，所以由上式可知Be1 - Hk+1,xy(k) = 0.由于Hk+1,是上Hessenberg矩阵，且hi+1,≠0,i=1,2,,k.通过向后回代求解可得y(k)=0,于是口β=0.这与ro≠0矛盾.所以hk+1,k=07.2.4带重启的GMRES方法由于随着迭代步数的增加，GMRES方法的每一步所需的运算量和存储量都会越来越大因此当选代步数很大时，GMRES方法就不太实用.通常的解决方法就是重启，即事先设定一个重启选代步数k如k=20或50等等，当GMRES达到这个选代步数时仍不收敛，则计算出方程组在(0)+Kk中的最佳近似解r(),然后令(0)=(k),并重新开始新的GMRES选代.不断重复该过程,直到收敛为止
· 220 · 第七讲 子空间迭代方法 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 方法中, 当执行到某一步时有 hj+1,j = 0, 则算法会中断 (breakdown). 如果出现 这种中断, 则我们就可以找到精确解. 定理 7.3 设 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 迭代. 不断重复该过程, 直到收敛为止