第七讲子空间迭代方法子空间迭代算法的基本思想是在一个维数较低的子空间中寻找解析解的一个“最佳”近似，子空间选代方法的主要过程可以分解为下面三步：(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 的算法