（计算数学专业论文）解非对称线性系统的两类方法.pdf

摘要 摘要 直接法和迭代法是解非对称线性系统的两类重要的方法与直接法相比， 缺少稳定性是迭代法一个公认的缺点，而预处理技术不仅可以提高迭代法的 稳定性也可以提高其收敛效率本文将侧重讨论k r y l o v 子空间迭代法的预处理 子构造和重启型左共轭方向法 第二章首先综述了b r u 提出的i s m ( a i s m ) 【1 构造预处理子的方法分析 表明这个过程不仅可得到非对称矩阵a 的三角分解还可得到其逆的三角分解， 借助于该过程，我们给出了不完全三角分解预处理子数值试验说明了新预处 理子的有效性 第三章讨论了重启左共轭方向法针对l c d ( m ) 中选择一个“合适”的重启 参数的困难，我们采用了选择变化的重启参数的技巧，这种技巧首次由b a k e r 等 在 2 提出并成功地应用至i j g m r e s ( m ) 方法在新算法中仅需要给出最大重启 次数和最小重启次数，每次迭代算法自适应地在已设定的范围内选取合适的 重启参数文章给出的数值结果表明新的算法较传统的l c d ( m ) 方法有效且稳 定 关键词预处理矩阵，非对称线性系统，广义极小残量法，左共轭方向法，重启 参数 a b s t r a c t a b s t r a c t d i r e c tm e t h o da n di t e r a t i v em e t h o da r et w oi m p o r t a n tm e t h o d sf o rs o l v i n gn o n s y m m e t r i cl i n e a rs y s t e m s c o m p a r e dt od i r e c ts o l v e r s ，l a c ko fr o b u s t n e s si sw i d e l y r e c o g n i z e da st h ew e a k n e s so fi t e r a t i v es o l v e r s b o t ht h ee f f i c i e n c ya n dr o b u s t n e s so f i t e r a t i v et e c h n i q u ec a nb ei m p r o v e db yp r e c o n d i t i o n i n g t h ec o n c e n t r a t i o n so ft h i s p a p e rw i l lb et h ec o n s t r u c t i o n so fp r e c o n d i t i o n e r sf o rk r y l o vi t e r a t i v em e t h o da n dt h e r e s t a r tl e f tc o n j u g a t eg r a d i e n tm e t h o d i nc h a p t e rt w o ，as i m p l er e v i e wo ft h ei s m ( a i s m ) m e t h o d p r o p o s e db yb r ue t a 1 j 】，f o rc o n s t r u c t i n gp r e c o n d i t i o n e ri sf i r s tg i v e n t h e o r e t i c a la n a l y s i ss h o w st h a t i nt h i sp r o c e s s ，b o t hl o w e ra n d u p p e rt r i a n g u l a rf a c t o r so f t h el ，d e c o m p o s i t i o no fa a n dt h ei n v e r s e so ft h e mc a nb eo b t a i n e d a sar e s u l t a ni n c o m p l e t et r i a n g u l a rf a c t o r - i z a t i o np r e c o n d i t i o n e ri sa c h i e v e d s o m ep r e l i m i n a r yn u m e r i c a lr e s u l t sd e m o n s t r a t e t h a tt h en e w p r o p o s e dp r e c o n d i t i o n e ri se f f i c i e n t t h er e s t a r tl e f tc o n j u g a t eg r a d i e n tm e t h o d ( l c d ( m ) ) i sd i s c u s s e di n c h a p t e r t h r e e 。t oo v e r c o m et h ed i f f i c u l t yo fc h o o s i n g s u i t a b l e f i x e dr e s t a r tp a r a m e t e r s ，a t e c h n i q u ef o rc h o o s i n gav a r y i n gr e s t a r tp a r a m e t e r si sa d o p t e d ，w h i c hw a sf i r s tp u t f o r w a r di n 2 】a n d a p p l i e d t og m r e s ( m ) s u c c e s s f u l l y i ti so n l yn e e d e dt op r o v i d et h e m a x i m u mr e s t a r tp a r a m e t e ra n dt h em i n i m u mo n e i nt h eg i v e ni n t e r v a l ，t h es u i t a b l e r e s t a r tp a r a m e t e r sc a nb ep r o d u c e da d a p t i v e l yf o re a c hi t e r a t i o n w er e p o r ts o m e n u m e r i c a lr e s u l t s ，w h i c hi n d i c a t et h a tt h ea d a p t i v el c d ( m ) o u t p e r f o r m st h el c d ( m ) b o t hi ne f f i c i e n c ya n di ns t a b i l i t y k e y w o r d sp r e c o n d i t i o n e dm a t r i x ，n o n s y m m e t r i cl i n e a rs y s t e m s ，g m r e s ，l e f tc o n j u g a t ed i r e c t i o n ，r e s t a r tp a r a m e t e r i v 第1 章前言 第1 章前言 数学，物理，流体力学和经济学等中的许多问题最终都可以归结为求解如 下线性系统： a x = b ， ( 1 0 1 ) 其中a r 似n 为给定的大型稀疏非对称矩阵，b r 几为给定的向量本文将考 虑解线性系统( 1 0 1 ) 的两类迭代方法 对于迭代法的研究一直是国内外研究的热点问题，早在l9 2 3 年， g a u s s 给c e r l i n g 的信中提及迭代法的巧算技巧 3 迭代法一般分为两类：一 类是基于矩阵分裂的定常迭代法，如j a c o b i ，g a u s s s e i d e l ，s o r 和a o r 等，这类 迭代法形式简单易于计算机实现，因此受到了工程人员的青睐然而随着社会 与科技的进步所求线性系统的规模越来越大，定常迭代法往往非常的低效，已 经很少单独使用了另一类是非定常迭代法，其中最广泛的便是k r y l o v 子空间 方法，最早可以追溯至5 0 年代初l a n c z o s 4 1 和a m o l d i 5 的上作我们知道子空 间方法的稳定性及收敛速率与系数矩阵的特征值分布有着密切的关系，特征 值分布越集中，算法就越稳定且收敛速度较快；反之则越慢借助于预处理技 术，k r y l o v 子空间方法的稳定性和收敛速度可以得到显著的提高近几年来预 处理子空间方法已成为解决大规模的问题的首要选择，如【6 一1 9 】 所谓预处理技术就是将原来的线性系统转化成与其等价的线性系统，而 这个等价系统的系数矩阵的特征值要比原来集中，从而当用子空间方法去解 时，收敛速度会有明显的提高s a a d 在 8 中指出，预处理子m 通常是矩阵4 在某 种意义下的近似，例如使得1 1 1 一a m 怯最小时的m 矩阵m 还应满足：以其为 系数矩阵的线性系统m z = 妙易于求解，这是因为在预处理的过程中每次都需 要解线性系统： 或 m a x = m b a m u = b ( 1 0 2 ) ( 1 0 3 ) 第1 章前言 许多情况下，矩阵m 形如分解形式m = l ，则原线性系统( 1 0 1 ) 也可以转化成 如下等价形式： l - 1 a u 一1 秒= l 一1 b ，z = u 一1 耖( 1 0 4 ) 目前，预处理子可分为两类：显式预处理子和隐式预处理子一个预处理 子是隐式的，是指在应用中每步都要解一个或多个线性系统，现有的方法有 矩阵的不完全分解法，例如i l u ，儿q ，q r ，s s o r 及a d i ( 见 8 ，2 0 等) 而一 个预处理子是显式的是指在每步的迭代中只需计算一次矩阵向量乘积，方法 包括f r o b e n i u s 范数极小化( 该方法i 由b e n s o n 在【1 5 中首次提出) 以及基于矩阵逆 的三角分解的方法 9 1 4 ，1 6 ，2 1 由于显式预处理子的高效性及可并行计算 的特点，越来越多的专家开始研究它若矩阵1 存在l d u 分解，即4 = l d u ， 那么基于逆不完全三角分解的方法就是要计算形如m = z b w 丁a _ 1 的 稀疏近似逆，其中z u 一，w u - t 为上三角矩阵，且d d 例如a i n v 方 法【1 2 ，1 3 通过1 正交化一组单位向量得到这个分解事实上，这种分解也可 以通过对矩阵a 进行秩一递减公式得到，见 2 2 1 r a f i e i ；f i t o u t o u n i a n 在 2 3 】分析 中说明a i n v 方法不仅能计算出系数矩阵a 的三角分解还能得到其逆的三角分 解b r u 在 1 中基于逆s h e r m a n m o r r i s o n 公式得到了矩阵的近似逆，在2 0 0 8 年他 进一步推导出，该方法还可以得到矩阵逆的三角分解，另外对于对称正定矩阵 它可同时计算出矩阵的三角因子以及逆三角因子，那么该方法对于非对称矩 阵是否能得到同样的结论，这就是本文第二部分要讨论的问题 本文的第三章将讨论解非对称线性系统的另一类方法：左共轭方向 法( l c d ) y u a n ，g o l u b ，p l e m m o n s $ 1 c e c i l i o 在 2 4 定义了针对实正定线性方 程组的半共轭方向法，并提出了l c d 方法d a i ，y u a n 2 5 将该方法推广为适 用非奇异，非斜对称线性系统情形，证明了实左共轭梯度( l c o ) 方向的存在 性，提出了克服算法中断的技巧以及节省空间的有限内存l c g 算法w a n g ， d a i 2 6 给出l c g 的复形式算法，证明了其本质性质，并给出克服算法中断的 块技巧c a t a b r i g a 2 7 ，2 8 用于解决对流扩散问题，与经典的子空间方法b i c g s t a b 2 9 ，q m r 3 0 及g m r e s 3 1 算法相比，数值结果表明l c g 是求解非 对称线性方程组的有效方法 类似于g m r e s 及其它子空间方法，随着迭代步数的增加，左共轭方向法 所需存储空间以及计算量都会大大的增加，使得求解时间变得很长为了 克服这样的困难，s a a d 和s c h u l t z 在 3l 】对g r m e s 方法采用了一种重启技巧，记 一2 一 第1 章前言 为g r m e s ( m ) 其中m 为其最大迭代步数，即最大近 以k r y l o v 子空间维数在最 大迭代次数m 之内，每次循环算法达到最大的迭代步数后，把当前所计算的近 似解当作下一个循环的的初始值重新开始这种方法可以固定算法的存储量， 从而减少运算时间但重启型算法丢掉了前面的一些信息，因此较非重启型算 法迭代步数会有所增加这种重启技巧在l c d 方法中也有应用如 2 7 】 g r m e s ( m ) 及l c d ( m ) 可以克服快速增长的存储量以及工作量，但重启参 数m 的选择直接影响其效率经典g m r e s ( m ) 着i l c d ( m ) 方法中参数m 的值是固 定的，但在实际中这个固定参数的值是很难选择的，若m 的值太小，则收敛速度 会很慢，从而导致问题的解决时间变得很长：若所选参数仇的值太大，重启策 略就失去了意义因此，很难选出一个合适的重启参数值m 对于g r m e s ( m ) 方 法很多专家学者开始研究变化的重启值m ! t l 2 ，3 2 3 5 但对于l c d ( m ) ，由于残 量范数非单调递减的，很少有专家学者提出改进的重启策略其中 2 7 ，2 8 中的 重启格式仍是经典的重启型方法，即固定重启参数值在文章的第三部分通过 分析了l c d 方法残量的特点，提出了一种自适应的选取参数的策略 本文其它部分安排如下：第_ 二章我f f j 首先回顾了b r u 等在 1 中介绍的 逆s h e r m a n m o r r i s o n 的分解过程，进而重点研究了算法当中所产生矩阵之间的 一些关系，及算法的中断条件，最后利用所找的关系构造了一种不完全三角分 解预处理子，并通过数值试验验证了其有效性第三章研究左共轭方向法所产 生向量之间的关系，从而利用这些关系来提出一种自适应的l c d ( m ) 方法第 四部分是结论分析和展望 第2 章 基于s h e r m a n m o r r i s o n 公式的不完全三角 分解预处理子 2 1 i s m 分解及其性质 本小节首先描述基于精确的逆s h e n l l a n m o r r i s o n 公式得到的矩阵分解过 程，然后揭示所得矩阵之间的一些重要关系 设b r n n 是非奇异矩阵，z ，y r n 是任意向量，若r = l + r b 一1 x 0 ， 则口的秩一校正b + x y r 非奇异，其逆可由s h e r n l a n m o r r i s o n 公式给出，如下： ( b + x y 丁) 一1 = b 一1 一r - 1 b 一1 x y 丁b ，( 2 1 1 ) 令 z 南) 怎1 ， 批，1 为r n 中的两个向量序列，4 0 为一个非奇异矩阵它的逆 己知或容易计算，且有如下的等式成立： 记a 七：= a o + 叁1z 秒于，后= 1 ，他，贝 a k + 12a k + z 后+ 1 y k + 1 t ， ( 2 1 2 ) 且a n = a 若z 七，y k 满足 d k = 1 + y k 丁a 芒1 x k 0 ，则由( 2 1 1 ) 式，可得矩阵a 南的 逆： a - ；1 = a 芒1 一五1l 七- 一1 1 z 惫) ( 吾肖乏1 ) ，后= 1 ，孔 由于a _ 1 = 4 二1 ，故下式成立 一= 4 0 - 喜去( 锚训( 棚， ( 2 ) 令机= a 芒1 x k ，讥= a 晶y k ，k = 1 ，n 贝l j ( 2 1 3 ) 式可以写成如下形 一4 一 丁南 扩七 z n 随 + o 4 f f a 第2 章丛于s h c r m a n m o r f i s o n 公式的4 i 完伞二角分解预处理予 式： 其中： a = a 0 1 c d 一1 妒t ， = 渺1 ，2 ，咖竹】，妒= 【矽1 ，妒2 ，妒n 】， d = 2 ( 2 1 4 ) 从上述的过程中，注意到妣，讥及d 七的计算依赖于过程矩阵a k 一1 ，k = 1 ，2 ，n 事实上，并不需要计算矩阵a k 一1 ，k = 1 ，2 ，竹引入r n 中的两组 向量序列 u 凫) ：l ， v 。n ：1 ，b r u 等在【1 中给出了如下的结论： 定理2 1 1 令a ，a o 为两个非奇异矩阵， z 七 ：1 ，( 挑 ：1 为瓞n 中的两个向量 序列，且满足( 2 1 2 ) 假定对于七= 1 ，2 ，竹，如= 1 + 甄t n 七- 一1 1 x k 0 ，其 中a 南一1 = a o + 洁k - 1 1z i 订，则 及 u 七：z 扩曼v t a 。o 。l x k 地， u 七= z 凫一 - 地， l = l 讥：_ ，t 。- 。1 丝忱， 讥= 瓤一 。兰忱， 0 = ：l 有定义，并且有如下关系成立： ( 2 1 5 ) ( 2 1 6 ) a 芒1 x k = a i l u 七， 3 ，七t a 庇一1 = v t a o ， ( 2 1 7 ) d 七= 1 + y t a 0 1 u k = 1 + v t a 0 1 x k ( 2 1 8 ) 令矩阵u ，y 分别是以( 2 1 5 ) 和( 2 1 6 ) 定义的向量u ，为列向量所构成的 矩阵由定y 2 ：2 1 1 ，可以看出乱忌，v k 分别是对z 七，玑做一个双共轭化的过程得到 的因此矩阵u ，y 可能分别与矩阵x ，y 具有相似的性质或结构，其中x ，y 是以 向量z t ，y i 为列向量所构成的矩阵下面的引理便揭示了它们之间的关系 引理2 1 1 1 假设定理2 1 1 成立，存在唯一的上三角矩阵取，砖使得 x = u t x ，y = v r r( 2 1 9 ) 其中 t x = t r = 1 t 1 2 0 1 ； 00 00 1t 1 2 01 ； 00 00 = 学 i 0 ，z 七= e 凫，弧= ( o 七一口3 ) t ， ( 2 1 1 0 ) 一6 一 n 佗 竹 1 。 n n 1 “1 第2 章丛于s h e r m a n m o r r i s o n 公式的不完伞三角分解预处理子 及 a o = s 厶，s 0 ，z 七= a k 一8 e k ，y k = e k ，( 2 1 11 ) 其中n 七表示矩阵a 的第七列，a 七，n 6 分别表示矩阵a 及a o 的第尼行，e 惫表示单位矩 阵的第后列在上述的选取下，分解形式( 2 1 4 ) 可以简化为： a 一1 = 8 - 1 ，一s 2 以w 1 哆， ( 2 1 1 2 ) 令( z ) 危表示向量z 的第后个元素，则( 2 1 5 ) g q ( 2 1 6 ) 可分别写成如下形式： 及 z擎-；e、,uk x ku i ， = 一 k 一1 个 v k2y k 一善警 2 一乙焉f 忱， 向一1 u k2 x k 一萎等u i ， 2 一乙 万， 口k = y k 一喜等y i ， r _ 、i “ij b 口2 一厶i f ， ( 2 1 1 3 ) ( 2 1 1 4 ) ( 2 1 1 5 ) ( 2 1 1 6 ) 在这些假设下，b r u 等在 1 中已经给出算法1 来计算a i l 一a 一1 的分解( 详见 参考文献【1 ) 为了更好的了解i s m 分解，b r u 等【l 】中引入了一个辅助单位上三角矩阵w ， 并给出了下面的一个重要引理( 行形式) 引理2 1 2 1 】假定当s 0 时，u s ，k 和d 。= d i a g ( d i ，蠼，d s n ) 是由精确i s m 分 解得到的矩阵当s = 1 ，i g u ，v f i ；i d = d i a g ( d l ，以，以) 是由精确i s m 分解得 一 一 第2 章丛- i 二s h e r m a n m o r r i s o n 公式的4 i 完伞，i 角分解预处删了 到的矩阵则 其中w 的第凫列是 u8 = u k = v 一( 8 1 ) 彬 d s = 8 - 1 d ， 一一三k - 1 警t 毗 其中让七为u 的第七列，x k ，y k 是由( 2 1 1 0 ) 令8 = l 所得向量 ( 2 1 1 7 ) ( 2 1 1 8 ) ( 2 1 1 9 ) ( 2 1 2 0 ) 下面，我们揭示矩阵彤与矩阵野之间的重要关系，基于此提出了一个新的 预处理子 定理2 1 2 假设引理2 1 1 ，2 1 2 及条件( 2 1 1 0 ) 成立，则 u = 巧1 ，w = 丐1 ( 2 1 2 1 ) 证明： 由( 2 1 1 0 ) 可得x k = e 七，即引理2 1 1 中的矩阵x 为单位矩阵，因 此c ，= 砭1 为非奇异单位上三角矩阵欲证= t v l ，注意到w 和r 都是单 位上三角矩阵，所以是可逆的，由引理2 1 2 及仉= u ，可得 露磁= ( 城) t 啦， i 0 时精确的i s m 分 解( 2 1 1 2 ) 存在则 且 a 一1 ：u d 一1 w t 曙= d u 一s w t 结合定理2 1 2 ，便得到矩阵4 的l d u 分解 ( 2 1 2 7 ) ( 2 1 2 8 ) 定理2 1 5 假设引理2 1 1 ，定理2 1 2 以及定理2 1 4 成立，则有如下等式： = 碍d t x ：(：2129ad2 9 )= 彤 ： ( ： 一1 0 一 第2 章丛于s h e r m a n m o r r i s o n 公式的f i 完伞二角分解预处理了 证明：由( 2 1 2 1 ) 和( 2 1 2 7 ) 成立，即得上述结论 若选取条件( 2 1 1 1 ) ，则有如下类似的结论 定理2 1 6 假定条件( 2 1 1 1 ) 和引理2 1 3 y g 立，当s 0 时，精确的i s m 分 解( 2 1 1 2 ) 存在则 且 a 一1 ：形d 一1 9 r 玩= 矿刁d s f v ( 2 1 3 0 ) ( 2 1 3 1 ) 证明：类似于参考文献 7 中定理2 1 的证明 定理2 1 7 假设引理2 1 1 ，定理2 1 3 以及定理2 1 6 成立，则如下等式成立： a = 酃d 豉 ( 2 1 3 2 ) 比较( 2 1 2 9 ) 和( 2 1 3 2 ) ，它们都是矩阵a 的l d u 分解，因此取= 豉，t v = 元r d = d ( 该结论可由矩阵a 的l d u 分解的唯一性得到) 所以无论选 择( 2 1 1 0 ) 还是( 2 1 11 ) 都得到一个令人惊奇的结论，臣 j i s m 分解不仅可得到矩 阵w ( w 一，矿) ，还可得到它们的逆耳，取 2 2 基于s h e r m a n m o r s s i o n 公式的不完全三角分解算法 由于行形式及列形式具有相同的性质及结论，因此在下面的讨论中仅讨 论行形式( 2 1 1 0 ) 并且在 1 】中，b r u 等通过数值试验已经说明行形式要优于列 形式 本小节中我们首先给出基于( 2 1 2 9 ) 构造隐式预处理子的算法过程，然后 将所得新预处理子与近似逆预处理子进行一个简单的比较e h ( 2 1 2 9 ) 矢n ，当采 取一个舍弃策略时，便得到矩阵a 的一个不完全三角分解算法如下： 算法1 ( 计算不完仝三角分解( 2 2 9 ) ) 1 a o = s i ，x 七= e 七，y k = ( a 惫一8 e 忌) t ，k = 1 ，礼 2 f o rk = 1 ，nd o ： u k2x k v k2 凫 f o ri = l ，佗d o ： 计算c _ s d i 雨l t e m p = 毕 i f i t e m p i t o l 令( t x h k = t e m p 计算u k = u k 一( t x h k u t e n d l f 计算e m p = 牮 i f l t e m p i t o l 令( r ) t 南= t e m p 计算v k = v k 一( r ) 伽忱 e n d l f e n d d 0 舍弃向量u 七及讥中小于t o l 的元素 计算蜈1 + 学 若d 七很小，算法终止 e n d d o 3 返回t x ，t y 和仇= d i a g ( d ，d 8 ，蛲) 注2 2 1 在文 【v l ，y 2 ，u 礼】， 章 1 中，算法1 返回的矩阵是u = 【? 2 1 ，u 2 ，u n 】，v = 及d = d i a g ( d 1 ，d 2 ，如) 从而，定义了如下的一个预处理 一1 2 一 第2 章丛十s h e r m a n m o r r i s o n 公式的f i 完伞i 角分解预处理了 矩阵： m 1 = 8 - 1 厶一8 - 2 u d s v t ( 2 2 1 ) 而我们的算法返回的是矩阵取，乃及d 。= d i a g ( d ，a i ，嚷) 从而我t f - i 以 得到如下预处理矩阵： 鸲= $ - 1 贩1 d i l 丐丁 ( 2 2 2 ) 当t o l = 0 时，预处理矩阵尬和m 2 是等价的当t o l 0 时，从算法1 我 们容易看出计算两个预处理矩阵的运算是相同的，但实际上预处理矩阵m 和 如是不同的，因为蚴是稀疏近似逆预处理矩阵( 显式预处理子) ，而m 2 是基 于矩阵a 的不完仝分解( 隐式预处理子) 为了保持稀疏性，b r u 等在【1 中通过数值试验指出对矩阵u 采用绝对舍 弃策略，对矩阵y 采用相对舍弃策略是最优的，最后，考虑参数s 取值，b r u 等 在 1 中也指出选取s = 1 5 i | a i i o 。可以使得迭代次数与存储空间达到一个平衡 因此在下面的数值试验我们也选择这样的舍弃策略及参数s 的值 2 3i s m 分解及a i s m 分解的存在性 本小节我们主要讨论的是何时i s m 分解及a i s m 分解可能出现失败的情 形 若i s m 分解过程可以不中断地执行，则所得矩阵r 及取为单位上三角的， 且满足( 2 1 2 9 ) ，令4 = l b o 是矩阵4 的l d u 分解，其中为单位下三角矩阵， d 为单位上三角阵，d 为对角矩阵且元素是递减排列的由于这种分解的唯一 性，我们得到j r s m 分解所得矩阵矸，取及d 满足： t r = l t ，t x = d ，d = d ， 从而由算法所计算的讲o = 1 ，2 ，n ) 为矩阵a 的主元如果令i 表示矩 阵a 的第i 阶( 1 i 几) 主子式的值并且令o = 1 ，则等式( 2 1 2 9 ) 隐含下列事 实： d a = # ，i = 1 ，2 ，n 籀2 章丛1 j s h c n n a n m o r r i s o n 公式的不完伞i 角分解预处删! 了 从上面的分析可得到如下事实： 定理2 3 1 假设定理2 1 5 成立，精确的i s m 分解可以执行当且仅当矩阵a 的所 有主子式值均不为0 类似于l u 分解，对于任意的非奇异矩阵a ，总存在置换矩阵尸( 或q ) 使 得i s m 过程作用于矩阵p a ( 或a q ) 不中断在 1 中，b m 等已经证明若月是非奇 异m 矩阵，则所得到的主元均是正的，并且证明了不完全i s m 分解所得主元 不小于由精确i s m 所计算出的主元，从而说明了a i s m 对于m - 矩阵不中断对 于矩阵有这样的性质吗? 我们首先来回顾一下日- 矩阵的定义，若a = 随j 】 为一个m 矩阵，则4 = o t l 为日矩阵，其中 铲譬“翼： 由上述定义我们可以得到所有对角占优矩阵为日一矩阵 下面，我们将说明若a 为日矩阵，由精确的i s m 分解所计算的主元也不小 于与a 的伴随m 矩阵的主元文，并且由近似i s m 分解所得的主元也不小于矩 阵a 的伴随m 矩阵的主元以 定理2 3 2 假设a 为日矩阵，a 为其伴随m 一矩阵若d j f l l d i 分别表示由精 确i s m 分解作用于4 和石所计算得到的主元，则哦哦且若盔表示由近 f n i s m ( a i s m ) 作用于矩阵a 所得到的主元，则以d i 证明：为了分析方便，记( 痧，伊) 是i s m 分解过程作用于矩阵a 所得到的矩阵， 并且记( 岔惫，纨) 是对矩阵a 选择行形式所得到的向量 以下采用数学归纳法证明矩阵u ，v ，杪及矿的列向量满足 及 也七i u 七l20 ，( 2 3 2 ) ( 讥) t - i ( v k ) , js0 ， f o r i k ( 2 3 3 ) 如也 0( 2 3 4 ) 第2 章丛于s h e r m a n m o r r i s o n 公式的不完伞二角分解预处理了 f l t u l = x l = 圣1 = c 1 o 得矗1 | u 1i 0 另一方面，y l = ( a 1 和雪1 = ( a 1 8 e 1 ) t 除去第一个元素外分别与a 和a 的第一行相同， 秒1 及0 1 = 雪1 ，可得 则 一l ( u 1 ) t l = ( 0 1 ) i ，f o ri 1 d l = l + 1 8 秒；1 u l = 1 + 兰8 并色1 o 现假定( 2 3 2 ) - , - ( 2 3 4 ) 对于任意的i k l 都成立对于k ，有 七一1 0 l u k = 愀一 另外，对于i k ，有： j = 1 七一l ( v j ) k 。 i 3 蚓+ i ( v j ) k l l s 勺| _ 1 i j = 1 蛐一蒌1 警奶 未知一警奶 j 2 。 。 a k 七一l 0 - i ( ) 一= 一一 j = 1 亟v s r j j ) 七一1 - i ( y k ) t l 一 - i ( y k ) t i 1 5 一 j = 1 ( s 可1 y 吾) t s e l ) 丁 由u 1 = 卜肛咐丁恕 y 一 勺 s h 触 由于一i ( 玑) j i = ( 雪) j ，j 七及( 2 3 3 ) ，则 由( 2 3 3 ) 一( 2 3 5 ) ，有 弧a t 一i 吾l i u j i 歹 l + 考虑定理的第二部分，类似可得证毕 由定理2 3 2 ，当近似i s m 分解作用于矩阵时，其主元不小于精确i s m 分 一1 6 一 缸 u 岛 倒 ，二b 触 凫 u 岛 岸 第2 章坫于s h e r m a n m o r r i s o n 公式的不完伞三角分解预处理了 解作用于其伴随m 矩阵时所得到的主元，但是它们不一定大于精确i s m 分解 作用于日矩阵a 例如考虑h 矩阵 41 一l 4 一ll 一三 1 8 一1一三 1 8 一l l 41 l4 近似i s m 分解算法对于阀值不小于t o l = o 0 5 时返回的主元画小于由精 确i s m 产生的主元d 4 如矩阵a 不是日矩阵，近似i s m 分解过程可能中断例 如，将阀值为t o l = o 0 3 的近似i s m 分解过程作用于s p d 矩阵 有五= o ( 算法中断) 2 0 0 0 4 0 0 1 0 0 0 5 0 4 0 1 0 8 2 0 0 0 o l 0 1 0 2 0 0 3 9 6 o 0 2 0 0 5 0 0 1 0 0 2 2 0 0 2 4 数值实验 本小节给出一些数值例子，验证提出的不完仝三角分解预处理子结合子 空间方法g m r e s 的有效性，有关该方法的详细介绍见参考文献 8 ，3 i 所有数值试验都是在m a t l a b 7 0 上实现，方法为g m r e s 左预处理所用 测试矩阵选自矩阵市场 3 6 及f l o r i d a 大学提供的稀疏矩阵库 3 7 表2 1 中列出 了测试矩阵的一些特征，包括其名字，维数，非零元个数( n n z ) ，矩阵的来源领域 以及在没有预处理子的情况下g m r e s 方法的计算时间( 伽) 在所有的数值 试验中令右端向量b = a e ，e = 【1 ，1 ，1 n 初始解向量为x o = 【0 ，0 ，o n 算法的终止条件是i h | | 1 0 _ 6 | | r o i i ，r k = b x k ，或者迭代次数超过1 0 0 0 步， 本章的第_ 二小节已经指出参数s = 1 5 1 1 a l i 。是最优的选取，因此在下面 的数值试验中只选用这个参数值针对矩阵h o r l 3 l 并t t o r s i r r l ，表2 2 和2 3 列 出在不同的阀值下矩阵m 1 和m 2 的非零元素个数，迭代次数和迭代时间，这 里n n z ( m 1 ) = n n z ( u ) + n n z ( v ) ，n n z ( m 2 ) = n n z ( t y ) + n n z ( t x ) 从表2 2 5 f 1 2 3 可 以看出，两种预处理子尬和m 2 作用于g m r e s 方法时，迭代步数几乎是相同 笫2 章坫1 ：s h e r m a n m o r r i s o n 公式的不完全i 角分解预处州了 表2 1 测试矩阵的特征 m a t r i xn7 0 ，凡2 d e s c r i p t i o n i t e r c p u h o r l 3 14 3 44 7 1 0n e t w o r kf l o w4 0 7o 7 5 f s 5 4 1 45 4 l4 2 8 5c h e m i c a lk i n e t i c s1 8 30 2 6 o r s i r r l1 0 3 06 8 5 8r e s e r v o i rs i m u l a t i o n4 3 8 1 3 4 o r s i r r 2 8 8 6 5 9 7 0 r e s e r v o i rs i m u l a t i o n 3 3 5o 7 2 s a y l r 31 0 0 03 7 5 0r e s e r v o i rs i m u l a t i o n2 6 70 5 2 s h e r m a n l1 0 0 03 7 5 0r e s e r v o i rs i m u l a t i o n2 6 7o 5 2 o r s i 辽g 12 2 0 51 4 1 3 3r e s e r v o i rs i m u l a t i o n1 4 5o 3 1 3 b f w 7 8 2 a7 8 27 5 1 4 e l e c t r o m a g n e t i c sp r o b l e m 2 0 30 31 p d e 9 0 09 0 04 3 8 0p a r t i a ld i f f e r e n t i a le q u a t i o n1 0 30 0 7 8 o l m l 0 0 01 0 0 03 9 9 6 c o m p u t a t i o n a lf l u i dd y n a m i c s 4 5 71 3 9 表2 - 2 矩阵h o r l 3 1 的试验结果 m i 丁d f n n zi t e r t c 础 n 1 3 , 2i t e r 印u 0 】4 1 4 56 20 0 4 717 5 36 20 0 4 7 0 0 l 1 8 4 1 3 4 2o 0 3 1 3 5 7 9 4 20 0 3 0 9 0 0 0 52 9 1 5 34 10 0 3 r4 5 8 64 l0 0 3 l 0 0 0 15 5 4 0 92 00 0 1 67 6 6 52 0o 0 1 5 0 0 0 0 56 8 1 0 31 6o 0 1 59 l 6 31 6 0 0 1 5 0 0 0 0 l8 9 3 1 380 0 1 61 2 4 6 880 0 1 5 0 0 0 0 0 11 0 8 5 9 640 0 0 81 7 5 7 540 0 0 8 的，但用矩阵m 2 做预处理子所需的迭代时间较短且矩阵m 2 所需存储空间 比m 1 要小我们也可以通过矩阵舰a 与m 2 a 的特征值分布图来验证上述结果 由表2 2 和2 3 中的计算结果可知，当阀值t o l = 0 1 或t o l = o 0 1 时，迭代次 数与存储空间达到一个平衡，因此在下面的计算中仅选取这两个阀值表2 4 z j l j 出了计算结果，该表验证了算法的有效性同时画出了矩阵s h e r m a n l 的特征 图以及用矩阵m 1 a 年d m 2 a 预处理后的特征分布图，说明了新的预处理子方法 可以解决一些具有坏条件数的矩阵 第2 章坫于s h e r m a n m o r r i s o n 公式的不完令二角分解预处理了 表2 3 矩阵o r s i r r l 的试验结果 m lm 2 t d 2 n n zi t e r c p u 礼n z i t e r t c p “ o 14 4 3 46 00 2 1 82 9 9 66 0o 1 0 9 0 0 11 0 6 4 5 2 9 0 0 9 3 4 4 5 l 2 90 0 4 7 0 0 0 51 1 6 2 72 80 0 9 3 4 6 3 3 2 80 0 4 7 o 0 0 11 5 2 5 32 80 0 9 85 2 6 12 80 0 6 3 0 0 0 0 52 5 5 0 2180 0 9 46 8 7 21 80 0 3 2 0 0 0 0 l5 0 2 9 01 20 0 7 99 6 7 91 2o 0 3 1 0 0 0 0 0 12 0 2 0 9 450 0 3 1 1 7 2 6 0 5 0 0 1 8 图2 1 当丁d f = 0 1 时，矩阵m 1 a ( a 为图2 - 2 当t o l = o 1 时，矩阵m 2 a ( a 为 矩阵h o r l 3 1 ) 的特征分布图矩阵h o r l 3 1 ) 的特征分和图 2 5 小结 本章在文献 1 】及【7 基础之上进一步研究了逆s h e r m a n m o r r i s o n 分解的 性质