Krylov子空间迭代法

Krylov子空间方法February12,2020内容•子空间和Krylov子空间•Arnoldi算法–Arnoldi过程–Gram-SchmidtArnoldi–HouseHolderArnoldi•FOM–IOM–DIOM2子空间•空间–集合，元素都是向量–线性空间（向量空间）•线性空间（交换律，结合律，幺元性，零元性，可逆性，数乘分配律等）•子空间–线性空间的非空子集•包含零元素，并且满足加法和乘法的封闭性–扩张（符合记作span）•包含所有向量的最小子空间3Krylov子空间选定一个向量v，可以得到一个序列𝑣0=𝑣𝑣1=𝐴𝑣0=𝐴𝑣𝑣2=𝐴𝑣1=𝐴2𝑣……𝑣𝑘+1=𝐴𝑣𝑘=𝐴𝑘+1𝑣,…假定A的特征多项式为𝑃𝐴𝜆=det𝜆𝐼−𝐴=𝜆𝑛+𝛼𝑛−1𝜆𝑛−1+...+𝛼1𝜆+𝛼0根据Cayley-Hamilton定理有𝑣𝑛+𝛼𝑛−1𝑣𝑛−1+...+𝛼1𝑣1+𝛼0𝜈0=0即VP=-𝑣𝑛其中𝑉=[𝑣0,𝑣1,...,𝑣𝑛−1，𝑃=𝛼0,𝛼1,...,𝛼𝑛−1𝑇Krylov子空间：𝐾𝑘(𝐴,𝑣)=𝑠𝑝𝑎𝑛{𝑣,𝐴𝑣,...,𝐴𝑘−1𝑣Krylov矩阵：𝐾𝑘(𝐴,𝑣)=[𝑣,𝐴𝑣,...,𝐴𝑘−1𝑣4𝒚=𝑨𝒙Krylov子空间法•Galerkin原理在子空间中找到近似解z(m)使得残向量r0−Az𝑚与Lm中的所有向量正交，即求（r0−Az𝑚，wi）=0，i=1,2,3….，m记𝑉𝑚=(𝑣1,𝑣2，...，𝑣𝑚)∈R𝑛×𝑚𝑊𝑚=(𝑤1,𝑤2，...，𝑤𝑚)∈R𝑛×𝑚由𝑧𝑚=𝑦1𝑣1+𝑦2𝑣2+...+𝑦𝑚𝑣𝑚=𝑉𝑚𝑦𝑚其中𝑦𝑚=𝑦1,𝑦2,...,𝑦𝑚𝑇所以𝑊𝑚𝑇(𝑟0−𝐴𝑉𝑚𝑦𝑚)=0得到𝑧𝑚=𝑉𝑚𝑦𝑚=𝑉𝑚𝑊𝑚𝑇𝐴𝑚𝑉𝑚−1𝑊𝑚𝑇𝑟05Ax=b,|A|!=0给定任意的x(0)，令x=x(0)+zAz=r(0),r(0)=b-Ax(0)X(m)=x(0)+z(m)选定子空间Km和Lm以及他们的基{vi},{wi}得到基于Galerkin原理构成的算法𝐾𝑚=𝑠𝑝𝑎𝑛{𝑟0,𝐴𝑟0,…，A𝑚−1𝑟0•Arnoldi过程设𝑟0,𝐴𝑟0,…，A𝑚−1𝑟0线性无关构造Km的一组标准正交基（1）令𝑣1=𝑟0||𝑟0||（2）取𝑣2=𝐴𝑣1−ℎ11𝑣1即（𝑣2，𝑣1）=0𝐴𝑣1，𝑣1=（ℎ11𝑣1，𝑣1）ℎ11=𝐴𝑣1，𝑣1记ℎ21=||𝑣2||，则𝑣2=𝑣2/||ℎ21||（3）对于k=1，2，…，m-1取𝑣𝑘+1=𝐴𝑣𝑘−ℎ𝑖𝑘𝑣𝑖𝑘𝑖=1有ℎ𝑖𝑘=(𝐴𝑣k，𝑣i),i=1,2,…,k记ℎk+1,𝑘=||𝑣k+1||,𝑣k+1=𝑣k+1/||ℎk+1,k||（4）当k=m时ℎ𝑚+1,𝑚≠0,𝑣m+1=𝑣m+1/||ℎm+1,𝑚||ℎ𝑚+1,𝑚=0，令𝑣m+1=06Arnoldi方法-基本算法Arnoldi方法-基本算法•Arnoldi过程的矩阵表示由𝑣𝑘+1=𝐴𝑣𝑘−ℎ𝑖𝑘𝑣𝑖𝑘𝑖=1𝑣k+1=𝑣k+1/||ℎk+1,k||得𝐴𝑣𝑘=ℎ𝑖𝑘𝑣𝑖𝑘+1𝑖=1，𝑘=1,2,…𝑚即𝐴𝑣1=ℎ11𝑣1+ℎ21𝑣2𝐴𝑣2=ℎ12𝑣1+ℎ22𝑣2+ℎ32𝑣3…𝐴𝑣𝑚=ℎ1𝑚𝑣1+ℎ2𝑚𝑣2+...+ℎ𝑚+1,𝑚𝑣𝑚+17H𝑚=ℎ11...ℎ1,𝑚−1ℎ1𝑚ℎ21...ℎ2,𝑚−1ℎ2𝑚⋱⋮⋮ℎ𝑚,𝑚−1ℎ𝑚𝑚𝐴𝑉𝑚=𝑉𝑚𝐻𝑚+ℎ𝑚+1,𝑚𝑣𝑚+1𝑒𝑚𝑇𝑒𝑚=0,0,...,1𝑇Arnoldi方法-MGS•Arnoldi-ModfiedGram-Schmidt–在Arnoldi算法的第j步，𝑤𝑗和ℎ𝑖𝑗𝑖=1𝑗可以耦合在一起计算8在数学上是等价的𝒉𝒊𝒋和w𝒋耦合计算有助于减少舍入误差的影响Arnoldi方法-HO•HouseholderArnoldi–目的：用Householder正交化子空间K（A，v1，m）(1)设||𝑣1||=1,令𝑧1=𝑣1，构造H矩阵P1，使得ℎ0=𝑃1𝑧1的后n-1个分量为0更新v1,𝑣1=𝑃1𝑒1(2)对Av2正交化，定义𝑧2=𝑃1𝐴𝑣1，构造P2，使得ℎ1=𝑃2𝑧2的后n-2个分量为0更新v2,𝑣2=𝑃1𝑃2𝑒2(3)继续上述过程9Arnoldi方法-FOM•完全正交化方法–在集合𝑥0+𝐾𝑚={𝑥0+𝑥:𝑥∈𝐾𝑚中找𝑥𝑚，使得估计解𝑥𝑚的残差满足𝑏−𝐴𝑥𝑚⊥𝐾𝑚由于Vm的列向量是Km的标准正交基，即𝑉𝑚𝑇(𝑏−𝐴𝑥𝑚)=0又𝑥𝑚可以表示为𝑥𝑚=𝑥0+𝑉𝑚𝑦𝑚所以：𝑉𝑚𝑇(𝑏−𝐴𝑥𝑚)=𝑉𝑚𝑇(𝑟0+𝐴𝑉𝑚𝑦𝑚)=𝛽𝑒1−𝐻𝑚𝑦𝑚𝛽=||𝑟0||可以得到：𝑦𝑚=𝐻𝑚−1(𝛽𝑒1)10Arnoldi方法-FOM•完全正交化方法𝑉𝑚𝑇(𝑏−𝐴𝑥𝑚)=0𝑥𝑚=𝑥0+𝑉𝑚𝑦𝑚𝑉𝑚𝑇(𝑏−𝐴𝑥𝑚)=𝑉𝑚𝑇(𝑟0+𝐴𝑉𝑚𝑦𝑚)=𝛽𝑒1−𝐻𝑚𝑦𝑚𝛽=||𝑟0||𝑦𝑚=𝐻𝑚−1(𝛽𝑒1)11Arnoldi方法-FOM（m）•完全正交化方法𝑏−𝐴𝑥𝑚⊥𝐾𝑚𝑥𝑚=𝑥0+𝑉𝑚𝑦𝑚𝑦𝑚=𝐻𝑚−1(𝛽𝑒1)当xm的残差很大时，将𝑥𝑚作为𝑥0重复算法即可得到FOM(m)12Arnoldi方法-IOM•不完全正交化方法当m=1时，𝑉𝑚=𝑣1=𝑟0𝛽𝐻𝑚=(𝐴𝑣1,𝑣1)𝑦𝑚=𝐻𝑚−1(𝛽𝑒1)=𝛽𝐻𝑚−1𝑥𝑚=𝑥0+𝐻𝑚−1𝑟0当m很大时，完全正交化需要保存所有的vi，内存开销很大在正交化过程中,v(j+1)仅与最近k个v(j-k+1),…v(j-1),v(j)正交,这样做虽然破坏的正交性,但是降低了计算量IOM算法仅仅是把FOM算法的第三步改为对i=max(1,j-k+1),…,j,计算hij与w(j).13Arnoldi方法-DIOM•直接不完全正交化方法采用IOM后，仍然需要存储v(1),v(2),…v(m),因为在第(vi)步中仍然需要这些向量.解决这个问题可以考虑采用H的LU分解，通过自身分解的迭代更新以减少每一步的存储量使xm的更新依赖于xm-1,14Arnoldi方法-DIOM•直接不完全正交化方法-𝑥𝑚的计算15H𝑚=ℎ11...ℎ1,𝑚−1ℎ1𝑚ℎ21...ℎ2,𝑚−1ℎ2𝑚⋱⋮⋮ℎ𝑚,𝑚−1ℎ𝑚𝑚lowerbidiagonalbandeduppertriangularArnoldi方法-DIOM•直接不完全正交化方法-𝑥𝑚的计算𝑥𝑚=𝑥0+𝑉𝑚𝑦𝑚𝑦𝑚=𝐻𝑚−1(𝛽𝑒1)定义因此有根据Lm的形式16Arnoldi方法-DIOM•直接不完全正交化方法-𝑥𝑚的计算𝑥𝑚=𝑥0+𝑉𝑚𝑦𝑚𝑦𝑚=𝐻𝑚−1(𝛽𝑒1)IOM方法和DIOM方法不再是正交投影而是一种斜投影方法。1718Thanksforyourtime!