大规模矩阵问题的Krylov 子空间方法综述科学工程培训讲义.ppt

大规模矩阵问题的krylov子空间方法综述 krylov子空间方法综述 背景介绍投影方法krylov子空间及其标准正交基krylov子空间方法求解线性方程组krylov子空间方法求解矩阵的特征值研究热点和尚未解决的问题 背景 大规模线性方程组的求解很多科学工程计算问题都转化为求解方程组ax b 如偏微分方程组的差分格式 有限元方法离散得到刚度矩阵 大规模矩阵特征值和特征向量的计算工程计算领域十分常见 如量子物理中的kohn sham方程求解化为哈密顿矩阵某些关键特征值对的计算 投影方法 线性方程组的投影方法方程组ax b a是n n的矩阵 给定初始x 0 在m维空间k 右子空间 中寻找x的近似解x 1 满足残向量r b ax 1 与m维空间l 左子空间 正交 即b ax 1 l此条件称为petrov galerkin条件 当空间k l时 称相应的投影法为正交投影法 否则称为斜交投影法 k l x 0 x 1 k l x 0 x 1 解方程组的投影法的矩阵表示设n m阶矩阵v v 1 v 2 v m 与w w 1 w 2 w m 的列分别构成k与l的一组基 记z x 1 x 0 z vy 有当非奇异时 通常情况成立 见定理1 1 有从而得到迭代公式 投影方法的最优性 误差投影 设a为对称正定矩阵 x 0 为初始近似解 且k l 则x 1 为采用投影方法得到的新近似解的充要条件是其中 且x为 1 1 的精确解 残量投影 设a为任意方阵 x 0 为初始近似解 且l ak 则x 1 为采用投影方法得到的新近似解的充要条件是其中 矩阵特征值的投影方法对于特征值问题ax x 其中a是n n的矩阵 斜交投影法是在m维右子空间k中寻找xi和复数 i满足axi ixi l 其中l为m维左子空间 当l k时 称此投影方法为正交投影法 kyrlov子空间 定义为m维krylov子空间 为v的次数 即使得q a v 的非零首一多项式的最低次数 krylov子空间的标准正交基 arnoldi方法基于gram schmit正交化方法首先 选取一个euclid范数为1的向量v 1 对 通常可取在已知v 1 v 2 v j 的情况下 不妨设v 1 v 2 v j av j 线性无关 否则构造完毕 则可求出与每个都正交的向量 而不难看出 再记 得到与v 1 v 2 v j 都正交的向量重复此过程 即可得到一组标准正交基 若期间某个j使得hj 1 j 0 则说明v的次数是j 且kj是a的不变子空间 定理如果记以v 1 v 2 v m 为列构成的矩阵为vm 由hij定义的 m 1 m阶上hessenberg矩阵为 删除最后一行得到的矩阵为hm 则在arnoldi算法中 可能有较大舍入误差 改写 krylov子空间方法解线性方程组 误差投影型方法取l k时的正交投影法 非对称矩阵的fom方法 解方程组的投影法的矩阵表示设n m阶矩阵v v 1 v 2 v m 与w w 1 w 2 w m 的列分别构成k与l的一组基 记z x 1 x 0 z vy 有当非奇异时 通常情况成立 见定理1 1 有从而得到迭代公式 回顾 不难看出 当采用上述fom算法时 需要存储所有的vi i 1 2 m 当m增大时 存储量以o mn 量级增大 而fom计算量是o m2n 可见其代价十分高昂 因此我们考虑重启的fom算法 krylov子空间方法解线性方程组 误差投影型方法取l k时的正交投影法 非对称矩阵的fom方法2 非对称矩阵的iom方法和diom方法 所谓不完全正交化方法 iom 是指在正交化过程中 v j 1 仅与最近k个v j k 1 v j 1 v j 正交 这样做虽然破坏的正交性 但是降低了计算量 当然k选得越小 对每个j对应的计算量也越小 但可能要选更大的m才能取得满足精度要求的近似解 iom算法仅仅是把fom算法的第三步改为 iii 对i max 1 j k 1 j 计算hij与w j 但采用iom后 仍然需要存储v 1 v 2 v m 因为在第 vi 步中仍然需要这些向量 解决这个问题可以考虑采用h的lu分解 通过自身分解的迭代更新以减少每一步的存储量 krylov子空间方法解线性方程组 误差投影型方法取l k时的正交投影法 非对称矩阵的fom方法2 非对称矩阵的iom方法和diom方法 对称矩阵lanczos算法 回顾 定理如果记以v 1 v 2 v m 为列构成的矩阵为vm 由hij定义的 m 1 m阶上hessenberg矩阵为 删除最后一行得到的矩阵为hm 则在arnoldi算法中 可能有较大舍入误差 改写 a是非对称阵时 h是hessenberg阵 则当a是对称阵时 h也为对称阵 故应为三对角阵 记为对它采用修正arnoldi正交化过程 就得到lanczos方法 krylov子空间方法解线性方程组 误差投影型方法取l k时的正交投影法 非对称矩阵的fom方法2 非对称矩阵的iom方法和diom方法 对称矩阵lanczos算法4 正定对称阵的cg法 cg法是解正定对称线性方程组最有效的方法之一该方法充分利用了矩阵a的稀疏性 每次迭代的主要计算量是向量乘法 而实际上 cg方法也是一种krylov子空间方法后面将给出一个数值例子 krylov子空间方法解线性方程组 残量投影型方法取l ak时的斜交投影法 gmers方法 gmers方法把方程组的求解化为一个极小问题的求解 回顾之前介绍投影类方法 回顾 残量投影 设a为任意方阵 x 0 为初始近似解 且l ak 则x 1 为采用投影方法得到的新近似解的充要条件是其中 gmers方法把方程组的求解化为一个极小问题的求解 回顾之前介绍投影类方法不去直接求x 1 转而去解此极小问题 这是gmers方法的独到之处 再由之前的定理 回顾 定理如果记以v 1 v 2 v m 为列构成的矩阵为vm 由hij定义的 m 1 m阶上hessenberg矩阵为 删除最后一行得到的矩阵为hm 则在arnoldi算法中 可能有较大舍入误差 改写 gmers方法在arnoldi正交化过程之后把方程组的求解化为一个极小问题的求解 回顾之前介绍投影类方法不去直接求x 1 转而去解此极小问题 这是gmers方法的独到之处 再由之前的定理 krylov子空间方法解线性方程组 残量投影型方法取l ak时的斜交投影法 gmers方法2 重启型gmers方法 qgmers dgmers gmers算法具有很好的收敛性 在不考虑舍入误差的情况下 该算法能在在有限步得到精确解 但是 和fom算法类似 它需要保存大量正交基 因此我们可以采取类似重启型fom的算法实现重启型gmers算法 同样可以采取不完全正交的方法 在正交化过程中 v j 1 仅与最近k个v j k 1 v j 1 v j 正交就能得到qgmers方法 为了避免存储v 1 v 2 v m k 可以对进行qr分解 这种算法被称为dqgmers krylov子空间方法解矩阵特征值问题 arnoldi方法求解非对称矩阵特征值lanczos方法求解对称矩阵特征值 arnoldi方法再次回顾前面的定理而所以有以下结论 1 若 则是a的不变子空间 从而hm的所有特征值是a的特征值子集2 若 设hm的特征值对是 即 由上述定理可得 以上讨论说明 可以用hm的特征值 又称ritz值 来近似a的特征值 用向量 又称ritz向量 来近似相应的特征向量 事实上 hm为a在kyrlov子空间上的正交投影 由于hm是上hessenberg阵 它的特征值求解难度远小于原问题 例如可以采取分解法来求解 arnoldi算法构造标准正交基1需要存储所有的基向量 当m很大时 存储量大2理论上为了保证收敛速度 m越大越好矛盾 lanczos方法a是对称阵时 hm是三对角阵 仍然采用hm的特征值来近似a的特征值 相应的ritz向量来近似a的特征向量 问题转化为三对角阵特征值的求解问题 而它的求解 复杂度大大减小 有很多成熟的办法 如分而治之法 qr法等 可以在并行机上达到to logn 的量级 对称lanczos算法 在没有舍入误差的情形下 得到的是标准正交基相互正交的 并且至多n步必然终止 但是在出现舍入误差的情况下 计算得到的将失去正交性 因此 长期以来被人们认为是不稳定的 直到1971年 paige通过仔细的舍入误差分析 发现失去正交性恰与近似特征值的精度提高有关 之后 经过大量的理论分析和数值实验 人们充分认识到lanczos方法是求解大型稀疏对称矩阵特征值问题的最有效方法之一 目前 lanczos方法是求大型稀疏对称矩阵部分极端特征值的最有效方法 改进技术 预条件法krylov子空间方法的收敛性一般都和矩阵的特征值分布有关 特征值分布越集中 收敛速度越快 若特征值分布较分散 则收敛的很慢 此时可以采用预条件子来改善矩阵的特征值分布 选取矩阵使得a的特征值比a集中 并且m 1容易求得 于是将原问题化为问题m 1ax m 1b 这是左预条件法 还可采用右预条件法和分裂预条件法 lanczos双正交化方法对于非对称线性方程组 也可采用类似于lanczos算法将非对称阵化为三对角阵 它的好处在于可以可以形成短迭代 而fom iom gmers的计算量和存储量都随着m的增大而增大在双正交化过程中 取容易看出a和在其中扮演对偶的角色 此方法特别适用于需要求解两个系数矩阵分别是a和at的方程组 多项式加速技术在求解矩阵特征值问题时 可以利用chebeshev多项式来加快收敛所谓多项式加速 就是采用作为迭代初始向量 其中为多项式 如果 其中pi为a的特征值对应的特征向量则将特征值按实部递减顺序排列 其中为r个 想要 的特征值 将 不想要 的特征值点集用s表示 如果多项式 能使得s尽可能小 则相当于让求解过程产生加速 而chebyshev多项式由于它的特殊性质 恰好能满足这种要求 不过 至今 仍无有效方法确定chebyshev多项式的某些参数 块方法 调和投影法当矩阵的特征值分布较密集或有重特征值时 常规的arnoldi方法或lanczos方法就必须取m很大时才能收敛 此时可以采用块方法取可以使得无须很大的m就可让迭代收敛 并且此方法适合改为并行算法 当特征值密集时 还可采用调和投影法 这也是最近研究的一个热点 当左右子空间的基满足wm avm时 称为调和投影 它除了可以通过选取位移点改善原矩阵的特征值分布外 还可以将内部特征值问题化为外部特征值问题 尚待解决的问题 在arnodli过程和lanczos双正交化过程中 会出现崩溃 怎样的矩阵容易出现崩溃 目前已有一些处理此崩溃的办法如前瞻技术 但并不十分成熟 在投影类方法中krylov方法是否是最优的 k是否有更好的选择 怎样将预条件技术运用到求解矩阵特征值问题 cg法解方程组ax b 输入 cg a b tol tol为误差限 输出 初始向量x0 随机产生 cputime 方程的解x 并绘制残差随迭代步变化的曲线 functionx cg a b tol t cputime n numel b x0 rand n 1 r b a x0 s b a x0 x x0 j 1 while norm r tol count j j 1 res j