krylov子空间算法

1、Krylov 子空间的定义：定义：令 RN，由 ，A，L ，Am1 所生成的子空间称之为由 与A所 生成的 m维Krylov 子空间，并记 K m A,v 。主要思想是为各迭代步递归地造残差向量，即第 n 步的残差向量 rn 通过系数矩阵 A的某个多项式与第一个残差向量 r 0相乘得到。即 n0r p A r 。但要注意，迭代多项式的选取应该使所构造的残差向量在某种内 积意义下相互正交，从而保证某种极小性（极小残差性） ，达到快速 收敛的目的。Krylov 子空间方法具有两个特征： 1. 极小残差性，以保证收敛速 度快。 2. 每一迭代的计算量与存储量较少，以保证计算的高效性。投影方法线性方程

2、组的投影方法方程组Ax b,A是n n的矩阵。给定初始 x0,在m维空间 K（右子空 间）中寻找 x的近似解 x1满足残向量 r b Ax1与 m维空间 L（左子空 间）正交，即 b Ax 1 L ，此条件称为 Petrov-Galerkin 条件。当空间 K=L 时，称相应的投影法为正交投影法， 否则称为斜交投 影法.投影方法的最优性：. （误差投影 ）设 A 为对称正定矩阵 ,x0 为初始近似解 ,且 K=L,则 x1 为 采用投影方法得到的新近似解的充要条件是x 1 m0in zz x 0 K其中, z A x z ,x z 122.（残量投影）设 A 为任意方阵 ,x 0为初始近似解

3、,且L AK,则x1 为采用 投影方法得到的新近似解的充要条件是x 1 m0in zz x 0 K其中 z b Az 2 b Az,b Az矩阵特征值的投影方法对于特征值问题 Ax x,其中 A 是 nn 的矩阵，斜交投影法是在 m 维右子空间 K 中寻找 xi和复数 i满足 Axi ixi L,其中 L 为 m维左子空 间.当 L=K 时，称此投影方法为正交投影法 .误差投影型方法：取 L=K 的正交投影法非对称矩阵的 FOM 方法（完全正交法）对称矩阵的 IOM 方法和 DIOM 方法对称矩阵的 Lanczos 方法对称正定矩阵的 CG 方法残量投影型方法：取 L=AK 时的斜交投影法GM

4、ERS 方法（广义最小残量法）重启型 GMERS 方法、 QGMERS、DGMERSArnoldi 方法标准正交基方法：Arnoldi 方法是求解非对称矩阵的一种正交投影方法。 Arnoldi 算 法就是对非对称矩阵 A ，产生 Krylov 子空间 m A,r 0 的一组标准正交基的 方法。该 算法 构造 m A,rHessenberg矩阵h11h12h13Lh1mh21h22h23Lh2mHm0h32h33LMMOOOhm 1, m0L0hm, m 1hmm的一组标准正交基v1, v2,K , vm 和h11h12h13Lh1mh21h22h23Lh2m0h32h33LMHmMOOOh m

5、 1, m0L0hm,m1 mm0L00h m 1, m基于 Gram-Schmit正交化方法首先，选取一个 Euclid 范数为 1的向量 v 1 , 对 m A, , 通常可取112 ，在v 1 ,v 2 ,L ,vj 已知的情况下 ,不妨设 v1,v2,L ,vj ,Avj 线性无 关 （ 否 则 构 造 完 毕 ）, 则 可 求 出 与 每 个 都 正 交 的 向 量 wj Avj h1jv1 h2jv2 L hjjv j而不难 看出 hij Av j ,vi ,i 1,2,L ,j ，再 记 hj1,j wj ,wj ，得 到与 jv1,v2,L ,v j 都正交的向量 v j 1

6、w ，重复此过程，即可得到一组标准 hj 1, j正交基。若期间某个 j 使得hj 1,j 0，则说明 v 的次数是 j，且 j是 A 的不变子空间。Arnoldi 算法：（1）取向量 v1 ，满足 v1 1,j 1（2）按（ 2）式计算 hij,i 1,2,L , j ，再按（ 1）式计算 w j（3）按（3）式计算 hj 1,j ，若hj 1, j 0 ，则停止，否则按（4）式计算 vj 1（4）若 j m，则 j j 1，转（ 2），否则停止w j Avj h1jv1 h2jv2 L hjjv j（1）hijAv j , v i ,i 1, 2,L , jh j 1, j w j , w

7、 j 1 22）3）（4）定理：如果记以 v1,v2,L ,vm 为列构成的矩阵为 Vm ,由hij定义的（m+1） m 阶上 Hessenberg矩阵（假设一个 n n阶矩阵 A，在i j 1时，它的 aij 0，VmT AVmHm那么这个矩阵 A 就叫做 Hessenberg矩阵）为 Hm ,删除最后一行得到的 矩阵为 Hm ，则： AVm VmHm wm emVm 1Hm在 Arnoldi 算法中，可能有较大舍入误差，改写：j 1 i 1 ihij Av j h1jv1 L hi 1,jvi 1,vi ,i 1,2,L , j修正的 Arnoldi 算法：（1）取向量 v1 ，满足 v1

8、 1,j 1（ 2）计算 w j Av j（ 3）依次对 i 1,2,L , j ，计算 hij w j ,v i 与 w j w j hij v i（4）计算hj 1,j ，若hj 1, j 0 ，则停止，否则计算 vj 1（5）若 j m，则 j j 1，转（ 2），否则停止FOM （完全正交化）方法 非对称矩阵的 FOM 方法： 解方程组的投影法的矩阵表示设n m阶矩阵 V v1,v2 ,L ,vm 与W w1,w2 ,L ,wm 的列分别构成 K 与L 的一组基。记 z x1 x0,z Vy，有WT r0 AVy 0当WT AV非奇异时，有 y WTAV WTr 0 ，从而得到迭代公式

9、：10xxVTW AV1 T 0 WrFOM算法：（1）计算 r 0b Ax 000， r 0 ,r 012 ， 1 r ， v r0 ，置 j 1（2）计算 w jAv j（3）依次对 i1,2,L , j，计算 hijw j ,vi 与w j w j hijv（4）计算 hj 1,j，若 hj1,j 0 ，则置m j ，转（6）（5）计算 v j 1，若 jm ，则置 jj 1 ，转（2）（6）按下式计m 算xm 0 m m 1 1 x x Vmy ,y H m e不难看出，当采用上述 FOM 算法时，需要存储所有的 vi ， （i=1,2, m），当 m 增大时，存储量以 O mn 量级

10、增大 .而 FOM 计算量是 O m2n 可见其代价十分高昂 .因此我们考虑重启的 FOM 算法 重启型 FOM 算法：（1）计算 r 0 b Ax0 , r 0 ,r 0 2,v1 r 0（ 2）生成 m A,r 0 的一组标准正交基，得到 Hm（3）按下式计算 xm ，若 xm 满足精度要求， 则停止，否则置 x0 x m 转（1）。 x m x0 Vmy m ,y m Hm1 e1IOM 方法非对称矩阵的 IOM 方法所谓不完全正交化方法 （IOM）, 是指在正交化过程中 ,v j 1 仅与最近 k 个v j k 1,L ,v j 1,v j 正交,这样做虽然破坏的正交性 ,但是降低了计

11、算量当然 k 选得越小，对每个 j 对应的计算量也越小 ,但可能要选更大的 m 才能取得满足精度要求的近似解 .IOM 算法仅仅是把 FOM 算法的第三步改为 i max 1, j k 1,L ,j， 计算 hij 与 w j 。但采用 IOM 后，仍然需要存储 v1,v2,L,vm ,因为在第 (vi)步 xm x0 Vmym 中仍然需要这些向量 .解决这个问题可以考虑采用 H 的 LU 分解，通过自身分解的迭代 更新以减少每一步的存储量DIOM 算法：1)计算 r 0b Ax0 0 0 2 1 0， r ,r ， v r r 0 ，置 m 12)计算 w mmAv3)对 i max1,m

12、k1 ,L ,m ，依次计算 himm i m m i w ,v 与 w whijv4)计算 hm 1,mm,w1 2 m 1 m,v w hm5)4)式更新的 LU 分解，若 umm0，则停止6)3)式计算，按( 2)式计算 p，其中 i 0时， uim pi 07)1)式计算，若符合精度要求，则停止，否则 m m 1，1)转(2)mmum1umm10mlm,m 1m1iuimpm 1,m 1k 1,m m k 1,m2)3)xxuim him li,i 1ui 1,m, i m k 1,L ,m 14)l m,m 1hm, m 1, ummhmm l m,m 1 um 1,mu m 1,m

13、 1Lanczos 方法每一步所需存储量和计算量是常量，不会随着子空间维数的增加而改Lanczos 方法是求解大规模稀疏对称矩阵端部特征问题的一种常 用的 正 交投 影方 法 。它 在 Krylov 子 空间 上对对 称 矩阵采 用 Rayleign-Ritz 方法，得到某些特征值和特征向量的近似。基本思想是 给定一初始向量， 逐步地构造出由该初始向量和对称阵生成的 Krylov 子空间的标准正交基。 通过简单的三项递推公式将大规模对称阵投影 为小规模对称三对角阵， 然后用此三对角阵的特征对来得出原始矩阵 的近似特征对。由于三对角阵好的结构和小的维数，在准确运算下，变。Lanczos算法：1)

14、计算r0b Ax 000 r , r2)计算 wAv jj1jv，其中当1 时 jv j 1 03)计算 jw j ,v j与wjjv4)j 12 ，若j0，则置mj 转 (5)， 否则 置j1 vm ，则 j5)w j 1 ，若 j置 Tm tridiag i , i 1 , i 1 , Vm1，转( 2)v1 ,L ,v，计算m 1 1 yTme6)计算 x m x 0 Vmy mLanczos双正交化方法在双正交化过程中，取Km A,r0spanr 0 ,Ar 0 ,L ,Am 1r 0T 0 0 T 0 T m 1 0 Lm A ,r span r ,A r ,L , A r 容易看出

15、 A和 AT 在其中扮演对偶的角色，此方法特别适用于需要 求解两个系数矩阵分别是 A和 AT的方程组 .基于 Lanczos双正交化过程的双共轭梯度法 BiCG算法：00s r ，1）计算 r 0 b Ax0 ，选取 rt0 使得 r 0 ,rt0 0 ，置方向向量rt 0 ，并置 m=02）计算m m m m r,rtAs ,stm 1 m与 向 量 x m 1x mmms，m1 rmmr m As ，m m T m rtrtm ATt，如果 xm 1 满足精度要求，则停止3）计算参数m 1 m 1 m m r,rtrm m 1 m ,rt， 与 向 量 s rm1stm 1 mrtm st

16、 ，置 m=m+1,转（ 2 ）CG 方法CG 法用来解对称且正定方程组十分有效， 但若是拿来解非对称或 是非正定的线性方程组则会发生中断。 它是借由迭代的方式产生一序 列的方向向量用来更新迭代解以及残向量， 虽然产生的序列会越来越 大，但是却只需要存储少数的向量。当系数矩阵 A 相当大而且稀疏， 此时 CG 法几乎就是高斯消去法。 CG 法理论上虽然保证最多 n 步能 解出线性方程组 Ax b 的解，但是若系数矩阵是病态的，此时累进误 差会让 CG 法在 n 步后无法求得充分精确的近似解。CG 算法：1)计算 r 0 b Ax0 ，选取 s0 r0 ，置 j 02)计算参数 j r j ,r

17、 j As j ,s j ，更新向量xj1xjjs j 与残向量rj13)r j j As j ，若 x j 满足精度要求，则停止计算 j 1 r j1,r j 1 r j ,r j ,sj1 r j 1jj 1s ，置jj 1，转( 2)CG 法是解正定对称线性方程组最有效的方法之一，该方法充分利用了矩阵 A 的稀疏性，每次迭代的主要计算量是向量乘法。GMRES 方法GMRE算S法要求系数矩阵 A是非奇异，非对称的 n维方阵。 GMRES算法利用 Arnoldi 变换构造一正交基 Vm v1, v2 ,L ,vm 来表示 Krylov 子空GMRE方S法把方程组的求解化为一个极小问题的求解，

18、 不去直接求x1 ,转而去解此极小问题，这是 GMRE方S法的独到之处。z b Az0mr 0AVm y m1m vAVm y2Vm 1 e1 Hmy m1m eHmyGMRES算法的收敛性完全取决于系数矩阵 A的特征值的性质。GMRES 算法：1)计算 r 00b Ax ，10 0 2 r ,rr 0 ，置 j 12)计算 w jAv j3)依次对 i1,2,L , j ，计算 hij w j,v i 与 w j w j hijv i4)计算 hj 1,j，若hj 1,j0 ，则置 mj ，转( 6)5)计算 v j 1 ，若 j m ，则置 j j1，转( 2)6)求ymargmin e1

19、 H mz z求解最小二乘问题 y m2，计算m 0 m x x Vm yargmin e Hm z 的算法：（1）令 g1 ， i 1（2）计算hii2 hi21,i 12,c hii ,s hi 1.i ，置 hii（ 3）置向量 t hi.i 1:m ，计算 hi,i 1:m ct shi 1,i 1:m, hi 1,i 1:m chi 1,i 1:m st （ hi.i 1:m 表示 矩阵第 i 行从 i+1 列到 m 列的元素构成的列向量）（4）置 t0 gi ，计算 gi ct0, gi 1 st0（5）若 i m，转（ 2）（6）依次对 j m,m 1,L ,1，计算 gj gj

20、 hjmgm L hj,j 1gj 1 h 所得到的g1,g2,L ,gm 即为所求的 y mGMRES算法允许 Krylov 子空间维数增加到 n，而且总是在最大迭 代次数 n内终止运算；另一种重启型 GMRE算S 法严格要求子空间维数 为一个定值 m，在进行了 m次迭代后，以得到的最后迭代结果 xm 作为 初始点重新进行 Arnoldi 变换，当残余向量 r A bx满足 rTAr 0时， 终止计算。综合考虑时间和空间复杂度， 重启型的 GMRE算S 法更适合。 重启型 GMRES 算法：（1）计算 r 0 b Ax0 , r0 ,r0 2,v1 r0（ 2）生成 m A, r 0 的一组

21、标准正交基，得到 Hm（3）求 ym argmin e1 Hmz ，计算xm x0 Vmym ，若xm 满足精度要 求，则停止，否则置 x0 xm ，转（1）同样可以采取不完全正交的方法， 在正交化过程中 ,v j 1 仅与最近 k个 v j k 1,L ,v j 1,v j 正交就能得到 QGMRES 方法。为了避免存储v1,v2 ,L ,vm k ，可以对 Hm进行 QR分解.这种算法被称为 DQGMRES。QGMRES 算法：2）计算 w j Av j3）对 i max 1,j k 1 ,L , j ，计算 hijijj i j j w ,v 与 w whijvi4）计算 hj 1,j

22、，若hj 1, j 0，则置 mj ，转（ 6）5）计算 vj 1 ，若 j m，则置 j j1，转（ 2）6）求 y m arg min zHmz2，计算 x0mxVm y1）计算 r 0 b Ax0 , r0 ,r0 2,v1Krylov 子空间方法解矩阵特征值问题Arnoldi 方法求解非对称矩阵特征值由定理： AVm VmHm wm emVm 1HmVmT AVm Hm而 w hm 1,mv ，则有如下结论：1）若 hm 1,m 0 ，则A,v 是A 的不变子空间，从而 H的所有特征值是A的特征值子集。2）若 hm 1,m 0，则 H，由上述定理TTm 1 mhm 1,mmhm 1,m

23、veyieyi可得： AVm yi iVmyi的特征值对是 i,yi ，即Hmyi i yi以上讨论说明 ,可以用 Hm的特征值 （又称Ritz值）来近似A的特征值 , 用向量 xi Vmyi （又称Ritz向量）来近似相应的特征向量 .事实上 , Hm为A 在Krylov 子空间上的正交投影 .由于 Hm 是上Hessenberg阵，它的特征值求解难度远小于原问题 ,例如 可以采取分解法来求解 .求解矩阵特征值的 Arnoldi 方法：（1）给定 m，取向量 v1 ，满足 v1 1, j 1（ 2）计算 w j Av j（3）依次对 i 1,2,L ,j，计算 hij wj ,vi 与 wj

24、 wj hijvi（4）计算hj 1,j ，若hj 1, j 0 ，则停止，否则计算 vj 1（5）若 j m，则 j j 1，转（ 2）（6）计算 Hm hij 的特征值对 i,yi 及 A 关于 i的 Ritz 向量 xi Vmyi Arnoldi 算法构造标准正交基存在的问题：1 需要存储所有的基向量 , 当 m 很大时 , 存储量大2 理论上为了保证收敛速度 ,m 越大越好Lanczos方法求解对称矩阵特征值A是对称阵时 , Hm 是三对角阵，仍然采用 Hm 的特征值来近似 A的特 征值，相应的 Ritz向量来近似 A 的特征向量 .问题转化为三对角阵特征值的求解问题，而它的求解，复杂度大 大减小，有很多成熟的办法，如分而治之法， QR法等，可以在并行 机上达到 tO（logN）的量级 .求解对称矩阵特征值的 Lanczos方法：（1）计算 r 0 b Ax0 ， r ,r ， v1 r ，置 j 1（2）计算 w Avj v ，其中当 j 1时 jv0（3）计算 j w ,v 与 w w j v（4）计算 j1 wj,wj 12 ，若 j 1 0，则置 m j 转（5），否则置 vj1 wj ，若 j m，则