第七讲 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 的算法