数值模拟导论-第七讲 Krylov子空间矩阵解法 雅克比·怀特 感谢Deepak Ramaswamy, Michal Rewienski,Karen Veroy and Jacob White ·回顾GCR －最小残向量解法 －Krylov子空间 －与多项式关系 ·回顾特征值和范数 －诱导范数 －谱半径定理 ·收敛速度的评估 －Chebychev多项式 ·预处理 －对角预处理 －近似LU预处理 概要 GCR算法 SMA-HPC ?2003 MIT {} 10 ,, ? ??? k ww GG for j＝0 to k－1 for i＝0 to j-1 正交搜索方向 标准化 更新结果 更新残向量 00 rbAx=? ( ) () T j jjii p pMpMpp←? ()() 1 T jj p p Mp Mp ← ( ) ( ) 1 T jjj jj xxrMpp + =+ ( ) ( ) 1 T jjj j j rrrMpMp + =? j j p r= 残留值便是下一次的搜索方向 GCR算法 SMA-HPC ?2003 MIT 标准化 图示运算步 1）正交化 2）解计算r的最小值 i Mrs ′ k x ·第一步搜索方向 GCR算法 SMA-HPC ?2003 MIT 11 2 2 ... NN x MxM xM b+++ = GG G kk Mxbr ?= 开始的几步 00 rbMxb= ?= 0 0 0 r p Mr = ( ) ( ) 10 00 T xrMpp= () 1 1.0 0 1 1 1.0 0 rp p M rp β β ? = ? ， ·残向量最小值解法 ·第二步搜索方向 GCR算法 SMA-HPC ?2003 MIT 11 2 2 ... NN x MxM xM b+++ = GG G kk Mxbr ?= 开始的几步（续） ·第三步搜索方向 ·残向量最小值解法 ( ) ( ) 21 1 11 T xx rMpp=+ 220020 2,1 2,0 rbMxr Mr Mrγγ=? = ? ? () 1 2,0 0 2,1 1 2 1 2,0 0 2,1 1 rpp p M rpp ββ ββ ?? = ?? GCR算法 SMA-HPC ?2003 MIT 11 2 2 ... NN x MxM xM b+++ = GG G kk Mxbr ?= GCR的第k步 ()() 1 0 k T kk k jj j pr Mr Mpp ? = =? ∑  k k k p p Mp =   ( ) () T k kk rMpα = 1kk kk xx pα + =+ 1kk kk rr Mpα + =? 正交和标准化搜索方向 第k步运算的最佳步长 更新结果和残余值 如果在GCR中对所有的有，那么 1） 2）是k次多项式，最小值为 3） jk≤ GCR算法 SMA-HPC ?2003 MIT 11 2 2 ... NN x MxM xM b+++ = GG G kk Mxbr ?= 多项式梗概 0 j α ≠ {} {} 00 01 , ,..., , ,..., k k p pp rMrMr=向量空间向量空间 () 10 , k k xMrξ + = k ξ 2 1 2 k r + ( ) ( )( ) ( ) 110 0 0 0 1 kk kkk r b Mx r M Mr I M M r Mrξξ ++ + =? = ? = ? ≡? 这里是（k＋1）次多项式，如果那 么他的最小值为 () 0 1k M r + ? ( ) 1 01 k+ ? = 2 1 2 k r + 如果属于向量空间最小化： 1）是k次多项式，最小化 2） 这里是（k＋1）次多项式，如果 那么他可以最小化 这里多项式作为解题工具，只有一个作用。那就是 最小化残向量。 Krylov方法 SMA-HPC ?2003 MIT 残向量最小值 多项式梗概 1k x + { } 00 , ,..., k rMr Mr 2 1 2 k r + () 10 , k k xMrξ + = ( ) ( )()( ) 110 0 0 0 1 kk kkk r b Mx r M Mr I M M r Mrξξ ++ + =? = ? = ? ≡? () 0 1k M r + ? ( ) 1 01 k+ ? = k ξ 2 1 2 k r + 2 1 2 k r + 吸热 Krylov方法 SMA-HPC ?2003 MIT “与外界物热交换 的例子” 绝缘棒和矩阵 近端温度远端温度 离散化 节点平衡方程 Krylov方法 SMA-HPC ?2003 MIT “与外界有热交换”的 例子 导体棒和矩阵 近端温度远端温度离散化 节点平衡方程 GCR性能 SMA-HPC ?2003 MIT （随机的Rhs） 反复迭代后的残向量对数图 GCR性能 SMA-HPC ?2003 MIT ( Rhs＝－1，＋ 1，－1，＋1….） 反复迭代后的残向量对数图 Krylov方法 SMA-HPC ?2003 MIT 残向量最小值 多项式最优化 残向量最小化的最优工具 是任意的k次顺序多项式，所以有 因此任意的满足零约束的多项式都可以获得上述同 的联系 ( ) ( ) 10 0 11 k kk rMrMr + ++ ≤? ≤?  1k+ ?  ( ) 1 01 k+ ? =  1 0 k r r + 诱导范数 SMA-HPC ?2003 MIT 矩阵的放大 倍数问题 假设，那么y比x 大多少？ 或者y相对于x 扩大了多少倍？ yMx= 诱导范数 SMA-HPC ?2003 MIT 回顾向量的 范数 L 2 范数： L 1 范数： 范数：L ∞ 矩阵范数 SMA-HPC ?2003 MIT 矩阵L阶范数 max max 1 l l Mx x x l l Mx M x ≡== 1 1 max n iij j M M = = ∑ 1 max n j ij i M M ∞ = = ∑ 定义 举例 最大列和 最大行和 =最大列绝对值之和 为什么？令 =最大行绝对值之和 为什么？令 不是很容易计算 通过前面的介绍我们知道，解X的值的变动同一个 因子有关，这个因子是关于A的。这个因子是， 起先他被定义为A的条件数，但是由于这个条件数是随 矩阵的范数变化，取代条件数A的是矩阵A的奇异值的 比率。 奇异值超出本课程的范围，要想知道请咨询Trefethen& Bau. 1 1 max n j ij i M M = = ∑ 1 A A ? 矩阵范数 SMA-HPC ?2003 MIT 矩阵L阶范数 [ ] 1,0,...0 T x = 1 max n jij i M M ∞ = = ∑ [ ] 1, 1,... T x = ±± ± 2 M max min () () () A Cond A A σ σ = 特征值和特征向量 应用 SMA-HPC ?2003 MIT 谱半径理论 i λ ( ) 01 ... p p f xxxαα α=+ ++ ( ) 01 ... p p f MMMαα α=+ ++ ( )( ) ( )( ) spectrum f M f spectrum M= 给定一个多项式 将多项式扩展到矩阵 那么就有 Krylov方法 SMA-HPC ?2003 MIT 收敛性分析 矩阵多项式的标准化 M特征空间的条件数 图中英文为：矩阵M的特征向量 Krylov方法 SMA-HPC ?2003 MIT 收敛性分析 矩阵多项式的标准化 Krylov方法 SMA-HPC ?2003 MIT 收敛性分析 重要发现 1）一个Krylov子空间残向量最小化的运算，最多N步收敛于 真实解。 证明：() ( )( ) ( ) ( ) () () 12 .... max 0 0 0 nni n ini n xx x x M Mr λλ λ λλ λ ?=? ? ? ∈ ?=?? = =   令，这里那么， ，并且 2）如果M只有q个独立的特征值，那么Krylov子空间残向量 最小化的运算，最多Q步收敛于真实解。 证明：() ( )( ) ( ) 12 .... qq xx x xλ λλ? =? ? ?  Krylov方法 SMA-HPC ?2003 MIT 对称矩阵的收敛性 多项式的残余量 如果M 是对称矩阵，那么 1）M有标准正交的特征向量 2）M有实数特征值 如果M正定，那么 ( ) 0Mλ > 导热棒矩阵的多项式残余量图 SMA-HPC ?2003 MIT 无热量散失情况（n=10） 保持尽可能的小； 在多项式的适当位置上故意放一些零。 导热棒矩阵的多项式残余量图 SMA-HPC ?2003 MIT 无热量散失情况（n=10） ( ) ki λ? Krylov方法 SMA-HPC ?2003 MIT 对称矩阵的收敛性 多项式的最值问题 ( ) [ ] min max min ,, 0Mλλλλ∈ > 考虑 那么一个多项式的好坏可以通过解最值问题来考察。 （小的多项式性能比较好）。() k pM 多项式的最大最小值问题已经被Chebyshev多项式解决 Krylov方法 SMA-HPC ?2003 MIT 对称矩阵的收敛性 Chebyshev方法解 最值问题 Chebyshev多项式： SMA-HPC ?2003 MIT Chebyshev多项式的最小化超出了[1，10] Krylov方法 SMA-HPC ?2003 MIT 对称矩阵的收敛性 Chebyshev的范围 Krylov方法 SMA-HPC ?2003 MIT 对称矩阵的收敛性 Chebyshev的结果 Krylov方法 SMA-HPC ?2003 MIT 前处理 对角矩阵的例子 是什么原因使GCR收敛更加迅速？ Krylov方法 SMA-HPC ?2003 MIT 前处理 对角矩阵的例子 让M=D+M 其中D是对角矩阵 应用GCR到 矩阵的逆在计算机中很容易求出 经常用来提高收敛性 ( )( ) 111 nd D Mx I DM x Db ??? =+ = 导热棒的例子 SMA-HPC ?2003 MIT 系统离散化 图中：一个小的 x? x? SMA-HPC ?2003 MIT 下面哪个收敛曲线是GCR 迭代 导体棒的例子 SMA-HPC ?2003 MIT 前处理矩阵 特征值 残余值最小化的Krylov子空间运算法则，可以通过 直接设置多项式零点来去除无关的特征值。 Krylov眼中的世界 SMA-HPC ?2003 MIT 热流动的例子 维数密度稀疏度GCR 在2维中，GCR比带状的GE收敛速度快，三维可 以更快，三维矩阵只有个m 3 非零数。GCR收敛 太慢。 Krylov方法 SMA-HPC ?2003 MIT 前处理 LU近似分解 的前处理 使M=LU L为下三角矩阵，U 为上三角矩阵。 将应用GCR到中。 用一个隐式矩阵来表示。 形式就等于解。 ( )( ) ( ) 1 LU M x LU b ? =   ( ) ( ) 1 y LU M x ? =  LUy Mx=  Krylov方法 SMA-HPC ?2003 MIT 前处理 LU近似分解的 前处理（续） 在真实的LU分解中的非零量 填充这个LU分解很麻烦，这里忽略填充 分解2维珊格矩阵 SMA-HPC ?2003 MIT 如果填充了，反而使矩阵分解难度加大 Krylov方法 SMA-HPC ?2003 MIT 前处理 LU近似分解前 处理（续） 放弃LU分解 放弃所有的填充 只放弃含有小数的填充 放弃填充由其他填充产生的零位 放弃填充由填充其他零位产生的零位， 等等 总结 SMA-HPC ?2003 MIT ?回顾GCR －最小残向量解法 －Krylov子空间 －与多项式关系 ?回顾特征值和范数 －诱导范数 －谱半径定理 ?收敛速度的评估 －Chebychev多项式 ?预处理 －对角预处理 －近似LU预处理