数值模拟导论 -第六讲 Krylov子空间矩阵求解方法 雅克比 ·怀特 感谢 Deepak Ramaswamy, Michal Rewienski,Karen Veroy and Jacob White 常规的子空间极小化算法 ——回顾学过的正交化和投射定理 GCR算法 － krylov-子空间 －对称矩阵的简化 －收敛条件 回顾特征值和特征向量 －范数和谱半径 －谱映射定理 概述 任意的子空间方法 SMA-HPC ?2003 MIT 可以近似为向量 的和 任意的子空间方法 近似的解 Mx＝ b 选择一个 k维的子空间 1 0 k k ii i x wα ? = ?= ∑ null {} 10 ,, ? ??? k ww nullnull { } 01 ,, k ww ? ??? null null ·残向量定义为 如果 残向量最小化的思路为：取 ，最小化式： 1 0 k k ii i xwα ? = =? ∑ null kk rbMx=? 1 0 k ii i bMwα ? = =? ∑ null 任意子空间方法 SMA-HPC ?2003 MIT 最小残向量 11 2 2 ... NN x MxM xM b+++ = nullnull null kk Mxbr ?= i sα ′ i sα ′ ()() 11 2 2 00 T kk T kkk ii ii rrrb Mwb Mwαα ?? == ???? ≡=? ? ???? ???? ∑∑ null null kk rbMx=? 如果 或者 正交于 ，那么我 们将很容易计算出 的最小值。 如果 或者 于 非正交 那么要建立一组正交的向量 使向量 空间 ＝向量空间 并且 ， ( ) ( ) 0 ij Mw Mw = nullnull ( ) i Mw null ( ) j Mw null 任意子空间方法 SMA-HPC ?2003 MIT 最小残向量 算法 2 2 2 2 k ii rb Mwα=? ∑ null { } 01 ,, k pp ? ??? null null {} 01 ,, k ww ? ??? null null () ( ) 0 T ij Mp Mp = nullnull i j≠ ( ) i Mw null ( ) j Mw null ()( ) 0 T ij Mp Mp = nullnull 给定 M， b并且一组搜索方向 1）通过将 正交化生成 for j＝ 0 to k－ 1 2）解 计算 r 的最小值 {} 01 ,, k ww ? ??? nullnull j Mws ′ j p s ′ null 任意子空间方法 SMA-HPC ?2003 MIT 最小残向量 运算步骤 ( ) () ()() 1 0 T k ji jj i T i ii Mw Mp p wp Mp Mp ? = =? ∑ ( ) () ()() ( ) () ()() 0 11 00 TT i kk ii k ii ii ii rMp rMp x pp Mp Mp Mp Mp ?? == ∑∑ k x 1）正交化 2）解计算 r 的最小值 任意子空间方法 SMA-HPC ?2003 MIT 最小残向量 图示计算步骤 j Mws ′ for j＝ 0 to k－ 1 for i＝ 0 to j-1 正交搜索方向 标准化向量 更新结果 更新残向量 00 rbAx=? 任意子空间方法 SMA-HPC ?2003 MIT 最小化算法 j j pw= ( ) () T j jjii p pMpMpp←? ()() 1 j j T jj p p Mp Mp ← ( ) ( ) 1 T jjj jj xxrMpp + =+ ( ) ( ) 1 T jjj j j rrrMpMp + =+ 选择 的标准 对所有在空间 中的 任意 ， 的值都很小 当 时 在向量空间 一种选择，单位向量， 向量空间 如果 k＝ N进行 QR分解 如果 k<N情况会很糟糕 01 {,, } k ww ? ??? null null i sα ′ 1 0 k k ii i bMxb Mwα ? = = ?? ∑ null 任意子空间方法 SMA-HPC ?2003 MIT 子空间的选择 标准 kN?? 1 A b ? ≈ { } 1 , k k xee∈ ??? null null 01 ,, k ww ? ??? null null 01 {,, } k ww ? ??? nullnull 求 的最小值，其中假设 （对称 矩阵）并且 推导出 x最小化 f. 取向量空间 这便是 f的最快下降方向，但 f并不是残余值。 这种方法不能用于非对成矩阵，和不满足 的情况 () 1 2 TT f xxMxxb=? 任意子空间方法 SMA-HPC ?2003 MIT 子空间的选择 传统方法 T M M= 0 T xMx> ( ) 1 x fx Mxb x Mb ? ?=??= {} ( ) ( ) { } 01 01 ,, ,, k kx x w w fx fx ? ? ??? = ? ???? null null 0 T xMx> 注意： 向量空间 ＝向量空间 如果： 向量空间 ＝向量空间 那么 并且向量空间 ＝向量空间 krylov子空间 01 { ,, } k rr ? ??? ( ) ( ) { } 01 ,, k xx fx fx ? ???? { } 01 ,, k ww ? ??? null null 任意子空间方法 SMA-HPC ?2003 MIT 子空间的选择 krylov子空间 1 0 0 k ki i i rr Mrα ? = =? ∑ { } 00 10 ,,, k rMr M r ? ??? 01 { ,, } k rr ? ??? 01 { ,, } k rr ? ??? 求解第 k步搜索方向的步长 更新结果 更新残向量 计算新的正交搜索方向 krylov方法 SMA-HPC ?2003 MIT GCR算法 GCP的第 k步 ( ) () ()() T k k k T kk rMp MpMp α = ( ) () ()() T k k k T kk rMp MpMp α = 1kk kk rr Mpα + =? ( ) ( ) ()() 1 11 0 T k k j kk j T j jj Mr Mp p rp Mp Mp + ++ = =? ∑ 向量内积， O(n)矩阵向量的内积 O(n)如果是稀疏矩阵需要 O(n)步运算 向量加， O(n) O(K)内积，总共需要 O(nk)步运算 如果 M是稀疏矩阵，当用 k步来逼近 n时，总共需 要运算步＝ O(n)+O(2n)+…+O(kn)=O(n 3 )具有很好 的收敛性。 ( ) () ()() T k k k T kk rMp MpMp α = 1kk kk xx pα + =+ 1kk kk rr Mpα + =? krylov 方法 SMA-HPC ?2003 MIT GCR算法 k步逼近的计算耗时 1kk kk rr Mpα + =? 我们会发现下面的情况： if M=MT 那么 j<k 一步完成正交 如果 ，那么矩阵对称，稀疏， GCR需要运 算步 更好的收敛速度 1kj rMp + ⊥ ( ) ( ) ()() ( ) () ()() 11 1 11 0 TT kk k jk k kjkk j kk jj Mr Mp Mr Mp p rppp Mp Mp Mp Mp ++ + ++ = =? ?? ∑ krylov 方法 SMA-HPC ?2003 MIT GCR算法 对称矩阵的情况 2 ()On knnull 吸热 krylov 方法 SMA-HPC ?2003 MIT “ 不与外界进行热交 换 ”的例子 绝缘棒和矩阵 近端温度 远端温度 将棒离散化 节点平衡方程 krylov 方法 SMA-HPC ?2003 MIT “不与外界进行热交 换 ”的例子 电路和矩阵 krylov 方法 SMA-HPC ?2003 MIT “与外界进行热交换 ” 的例子 导体棒和矩阵 近端温度 远端温度 离散化 节点平衡方程 krylov 方法 SMA-HPC ?2003 MIT “与外界进行热交换 ” 的例子 电路和矩阵 节点平衡方程 残余误差 迭代次数 反复迭代后的 log（残余误差）对比图 GCR性能 SMA-HPC ?2003 MIT （随机的 Rhs） GCR性能 SMA-HPC ?2003 MIT ( Rhs＝－ 1，＋ 1，－ 1，＋ 1….） 反复迭代后的log（残余误差）对比图 Krylov 子空间方法 SMA-HPC ?2003 MIT 收敛性分析 多项式逼近方法 {} { } 00 0 0 ,.... , ,.... k k ww rMrMr= () 100 0 k ki ik i k xMrMrαξ + = == ∑ nullnullnullnullnull 次多项式 (()) 10 10 0 0 k ki ik i rr MrIMMrαξ ++ = =? =? ∑ 0 0α ≠ { } { } 00 0 0 ,.... , ,.... k k ww rMrMr= 如果向量空间 注：对任意 的 存在向量空间 。 Krylov 方法 SMA-HPC ?2003 MIT 收敛性分析 基本知识 0 j α ≠ jk≤ { } { } 00 0 01 , ,.... , ,.... k k pp p rMr Mr= 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ξξ ++ + =? = ? = ? ≡? 0 1 () k M r + ? 2 1 2 k r + 1 (0) 1 k+ ? = 如果 ，在 GCR中对所有的 1）向量空间 2），是k次多项式，这个多项式可以最小化 这里 是（ k＋ 1）次多项式，要想最小化 ，必须使 3） ，有 Krylov 方法 SMA-HPC ?2003 MIT 收敛性分析 GCR多项式的最优 化处理 22 10 1 () k k rMr + + ≤? null 1k+ ? null 1 (0) 1 k+ ? = 2 1 2 k r + GCR最优化的性质 这里 是任意的 k次多项式，并且有 因此任意的满足零约束的多项式都可以获得上述同 的联系。 是矩阵 M的一个特征值 如果 为 特征向量 或者：如果满足 特征值和特征向量 简介 SMA-HPC ?2003 MIT 基本定义 iii Muuλ= null null i λ i u null i M Iλ? i λ i u null 矩阵的特征值和特征向量满足 其中 为特征值， 是奇异矩阵， 是矩阵 M的一个特征向量 ()0 i MIλ? = 特征值和特征向量 简介 SMA-HPC ?2003 MIT 基本定义 例子 特征值？ 特征向量？ 如果是一个下三角矩阵会怎样？ 特征值和特征向量 简介 SMA-HPC ?2003 MIT 注意问题 几乎所有的 N*N矩阵都有 N个线性独立的特征向量。 矩阵 M的一组特征值被称之为 M的型谱 i λ 特征值和特征向量 简介 SMA-HPC ?2003 MIT 注意问题 几乎所有的 N*N矩阵都有 N个线性独立的特征向量。 并不意味着独立的特征值， 可以等于 并不意味着 M是非奇异矩阵。 i λ i λ j λ 特征值和特征向量 简介 SMA-HPC ?2003 MIT 谱半径 i λ 矩阵 M的谱半径为：圆心位于原点，可以包含 所有矩阵 M的特征值的圆的最小半径。 特征值和特征向量 简介 SMA-HPC ?2003 MIT 热流例子 i λ 棒长 吸入热 特征值和特征向量 简介 SMA-HPC ?2003 MIT 热流例子（续） i λ 特征值和特征向量 简介 SMA-HPC ?2003 MIT 热流例子（续） i λ 四个特征向量－选择哪一个？ 特征参数 SMA-HPC ?2003 MIT 谱半径定理 给定一个多项式 将多项式应用到矩阵 那么： ( ) 01 ... p p f xxxαα α=+ ++ ( ) 01 ... p p f MMMαα α=+ ++ ( )() ( )( ) fM f M=谱半径 谱半径 特征参数 SMA-HPC ?2003 MIT 谱半径定理证明 注意矩阵下面的性质 应用到多项式矩阵的形式 分解因式 1121 1pp MMUUUU UU MUU λλ λ λ ? ?? ? == ?= () 11 1 01 ... p p f MUUUU UUααλ αλ ? ?? =+ ++ () ( ) () () 1 01 01 ... ... p p p p f MUI U fMU U I ααλ αλ ααλ αλ ? =+++ =+++ 11 2 2 ... NN xu u uα αα=+++ nullnull null 特征参数 SMA-HPC ?2003 MIT 谱分解 用特征向量的成分分解任意的 x 通过解 来计算 x 用矩阵 M来代替 x，得： 11 2 2 111 2 22 ( ... ) ... NN NNN Mx M u u u uu u α αα αλ αλ α λ =+++ =+++ null nullnull null nullnull Krylov 方法 SMA-HPC ?2003 MIT 收敛性分析 重要结论 1） GCR运算法则在最多 n步内收敛于真实解。 2）如果 M只有 q个特征值，那么 GCR运算法则最多以需 要 q步就可收敛 。 ( ) ( )( ) ( ) ( ) () 12 00 .... , nni n xx x x M Mr r λλ λ λλ?=? ? ? ∈ ?? = = null null 令其中 ，因此 证明： 总结 SMA-HPC ?2003 MIT ·任意子空间运算法则 －正交化搜索方向 ·GCR算法 － Krylov子空间 －对称矩阵的简化 －泄漏和绝缘例子 ·回顾特征值和特征向量 －谱半径理论 ·GCR极限情况 － Q步收敛法