krylov子空间算法

Krylov子空间的定义：定义：令NR，由1mA，，，A所生成的子空间称之为由与A所生成的m维Krylov子空间，并记,mKAv。主要思想是为各迭代步递归地造残差向量，即第n步的残差向量nr通过系数矩阵A的某个多项式与第一个残差向量0r相乘得到。即0nrpAr。但要注意，迭代多项式的选取应该使所构造的残差向量在某种内积意义下相互正交，从而保证某种极小性（极小残差性），达到快速收敛的目的。Krylov子空间方法具有两个特征：1.极小残差性，以保证收敛速度快。2.每一迭代的计算量与存储量较少，以保证计算的高效性。投影方法线性方程组的投影方法方程组Axb,A是nn的矩阵。给定初始0x,在m维空间K(右子空间）中寻找x的近似解1x满足残向量1rbAx与m维空间L(左子空间）正交，即1bAxL，此条件称为Petrov-Galerkin条件。当空间K=L时，称相应的投影法为正交投影法，否则称为斜交投影法.投影方法的最优性：１.(误差投影)设A为对称正定矩阵,0x为初始近似解,且K=L,则1x为采用投影方法得到的新近似解的充要条件是01minzxKxz其中,12,zAxzxz2.(残量投影)设A为任意方阵,0x为初始近似解,且LAK,则1x为采用投影方法得到的新近似解的充要条件是01minzxKxz其中122,zbAzbAzbAz矩阵特征值的投影方法对于特征值问题Axx,其中A是n×n的矩阵，斜交投影法是在m维右子空间K中寻找ix和复数i满足iiiAxxL,其中L为m维左子空间.当L=K时，称此投影方法为正交投影法.误差投影型方法：取L=K的正交投影法非对称矩阵的FOM方法（完全正交法）对称矩阵的IOM方法和DIOM方法对称矩阵的Lanczos方法对称正定矩阵的CG方法残量投影型方法：取L=AK时的斜交投影法GMERS方法（广义最小残量法）重启型GMERS方法、QGMERS、DGMERSArnoldi方法标准正交基方法：Arnoldi方法是求解非对称矩阵的一种正交投影方法。Arnoldi算法就是对非对称矩阵A，产生Krylov子空间0,mAr的一组标准正交基的方法。该算法构造0,mAr的一组标准正交基12,,,mvvv和Hessenberg矩阵1112131212223232331,,1000mmmmmmmmmhhhhhhhhhhHhhh，1112131212223232331,,11,000000mmmmmmmmmmmhhhhhhhhhhHhhhh基于Gram-Schmit正交化方法首先，选取一个Euclid范数为1的向量1v,对,mA,通常可取112,，在12,,,jvvv已知的情况下,不妨设12,,,,jjvvvAv线性无关(否则构造完毕),则可求出与每个都正交的向量1212jjjjjjjwAvhvhvhv而不难看出,,1,2,,jiijhAvvij，再记121,,jjjjhww，得到与12,,,jvvv都正交的向量11,jjjjwvh，重复此过程，即可得到一组标准正交基。若期间某个j使得1,0jjh，则说明v的次数是j，且j是A的不变子空间。Arnoldi算法：（1）取向量1v，满足11,1vj（2）按（2）式计算,1,2,,ijhij，再按（1）式计算jw（3）按（3）式计算1,jjh，若1,0jjh，则停止，否则按（4）式计算1jv（4）若jm，则1jj，转（2），否则停止1212jjjjjjjwAvhvhvhv（1）,,1,2,,jiijhAvvij（2）121,,jjjjhww（3）11,jjjjwvh（4）定理：如果记以12,,,mvvv为列构成的矩阵为mV,由ijh定义的(m+1)×m阶上Hessenberg矩阵（假设一个nn阶矩阵A，在1ij时，它的0ija，那么这个矩阵A就叫做Hessenberg矩阵）为mH,删除最后一行得到的矩阵为mH，则：1TmmmmmmmAVVHweVHTmmmVAVH在Arnoldi算法中，可能有较大舍入误差，改写：1111,,,1,2,,jiiijjijhAvhvhvvij修正的Arnoldi算法：（1）取向量1v，满足11,1vj（2）计算jjwAv（3）依次对1,2,,ij，计算,jiijhwv与jjiijwwhv（4）计算1,jjh，若1,0jjh，则停止，否则计算1jv（5）若jm，则1jj，转（2），否则停止FOM（完全正交化）方法非对称矩阵的FOM方法：解方程组的投影法的矩阵表示设nm阶矩阵12,,,mVvvv与12,,,m的列分别构成K与L的一组基。记10,zxxzVy，有00TWrAVy当TWAV非奇异时，有10TTyWAVWr，从而得到迭代公式：1100TTxxVWAVWrFOM算法：（1）计算00rbAx，1200,rr，01rv，置1j（2）计算jjwAv（3）依次对1,2,,ij，计算,jiijhwv与jjiijwwhv（4）计算1,jjh，若1,0jjh，则置mj，转（6）（5）计算1jv，若jm，则置1jj，转（2）（6）按下式计算mx011,mmmmmxxVyyHe不难看出，当采用上述FOM算法时，需要存储所有的iv，(i=1,2,…m)，当m增大时，存储量以Omn量级增大.而FOM计算量是2Omn.可见其代价十分高昂.因此我们考虑重启的FOM算法重启型FOM算法：（1）计算10200001,,,rrbAxrrv（2）生成0,mAr的一组标准正交基，得到mH（3）按下式计算mx，若mx满足精度要求，则停止，否则置0mxx，转（1）。011,mmmmmxxVyyHeIOM方法非对称矩阵的IOM方法所谓不完全正交化方法(IOM),是指在正交化过程中,1jv仅与最近k个11,,,jkjjvvv正交,这样做虽然破坏的正交性,但是降低了计算量.当然k选得越小，对每个j对应的计算量也越小,但可能要选更大的m才能取得满足精度要求的近似解.IOM算法仅仅是把FOM算法的第三步改为max1,1,,ijkj，计算ijh与jw。但采用IOM后，仍然需要存储12,,,mvvv,因为在第(vi)步0mmmxxVy中仍然需要这些向量.解决这个问题可以考虑采用H的LU分解，通过自身分解的迭代更新以减少每一步的存储量DIOM算法：（1）计算00rbAx，1200,rr，010rvr，置1m（2）计算mmwAv（3）对max1,1,,imkm，依次计算,miimhwv与mmiijwwhv（4）计算1211,,,mmmmmmhwwvwhm（5）按（4）式更新mH的LU分解，若0mmu，则停止（6）按（3）式计算m，按（2）式计算mp，其中0i时，0iimup（7）按（1）式计算mx，若符合精度要求，则停止，否则1mm，转（2）1mmmmxxp（1）11mmmimmimpuvup（2）1,110,1{mmmmlm（3）1,1,mkmmkmuh,11,,1,,1imimiiimuhluimkm（4）,1,1,11,1,1,mmmmmmmmmmmmmmhluhluuLanczos方法Lanczos方法是求解大规模稀疏对称矩阵端部特征问题的一种常用的正交投影方法。它在Krylov子空间上对对称矩阵采用Rayleign-Ritz方法，得到某些特征值和特征向量的近似。基本思想是给定一初始向量，逐步地构造出由该初始向量和对称阵生成的Krylov子空间的标准正交基。通过简单的三项递推公式将大规模对称阵投影为小规模对称三对角阵，然后用此三对角阵的特征对来得出原始矩阵的近似特征对。由于三对角阵好的结构和小的维数，在准确运算下，每一步所需存储量和计算量是常量，不会随着子空间维数的增加而改变。Lanczos算法：（1）计算00rbAx，1200,rr，01rv，置1j（2）计算1jjjjwAvv，其中当1j时10jjv（3）计算,jjjwv与jjjjwwv（4）计算121,jjjww，若10j，则置mj转（5），否则置11jjjwv，若jm，则1jj，转（2）（5）置111,,,,,mmiiimTtridiagVvv，计算11mmyTe（6）计算0mmmxxVyLanczos双正交化方法在双正交化过程中，取00001,{,,,}mmKArspanrArAr10000,{,,,}mTTTmLArspanrArAr容易看出A和TA在其中扮演对偶的角色，此方法特别适用于需要求解两个系数矩阵分别是A和TA的方程组.基于Lanczos双正交化过程的双共轭梯度法BiCG算法：（1）计算00rbAx，选取0tr使得00,0trr，置方向向量00sr，00ttsr，并置m=0（2）计算,,mmmmmttrrAss与向量1mmmmxxs，1mmmmrrAs，mmmTttmrrAt，如果1mx满足精度要求，则停止（3）计算参数11,,mmmmmttrrrr，与向量11mmmmsrs，11mmmttmtsrs，置m=m+1,转（2）CG方法CG法用来解对称且正定方程组十分有效，但若是拿来解非对称或是非正定的线性方程组则会发生中断。它是借由迭代的方式产生一序列的方向向量用来更新迭代解以及残向量，虽然产生的序列会越来越大，但是却只需要存储少数的向量。当系数矩阵A相当大而且稀疏，此时CG法几乎就是高斯消去法。CG法理论上虽然保证最多n步能解出线性方程组Axb的解，但是若系数矩阵是病态的，此时累进误差会让CG法在n步后无法求得充分精确的近似解。CG算法：（1）计算00rbAx，选取00sr，置0j（2）计算参数,,jjjjjrrAss，更新向量1jjjjxxs与残向量1jjjjrrAs，若jx满足精度要求，则停止（3）计算111111,