（计算数学专业论文）大型线性方程组不完全分解预条件方法的研究.pdf

摘要 摘要 对于稀疏矩阵彳来说，完全分解所产生的预条件子一般不能保证具有和矩阵 彳一样的稀疏性，往往稠密了很多。因此，为了使预条件子的稀疏结构不那么稠 密并且预条件效果也不受很大的影响，不完全分解方法被提了出来。 不完全c h o l e s k y 分解所产生的预条件子对于许多大型线性方程组来说是一类 非常有效的预条件子，文中提出的不完全c h o l e s k y 分解算法3 _ 4 所产生的预条件 子三中每列允许保留的非零元素个数介于仇与n 。+ p 之间，每列具体保留多少个， 通过最优参数f 来确定，由此预条件子的存储空间也就得到了控制。通过数值实验 可以看出，最优参数f 的值选择为该矩阵的f r o b e n i u s 范数的数量级的倒数值时似 乎为最佳，文中提出的新算法的p c g 迭代步数与每列保留以+ p 个元素的l i n 和 m o r 6 的算法3 2 的p c g 迭代步数相同，但存储空间比算法3 2 的要小。 从稀疏模式s 产生的时刻与不完全分解的整个过程的关系来说，稀疏模式可 以分为两种情况。第一种稀疏模式是在分解过程之前事先确定好了的静态稀疏模 式，第二种是在分解过程中产生的不可预测的动态稀疏模式。对于第二种不可预 测的动态稀疏模式，文中在被分解矩阵么是对称m 一矩阵的情况下给出了它的稳 定性证明。同时，从证明过程中也可以看出，在相同的参数p 和相同的p c g 迭代 步数下，文中提出的算法3 4 具有至少和l i n 和m o r 6 的算法3 2 一样的稳定性。 关键词：预条件子，不完全c h o l e s k y 分解，存储空间，舍弃策略，稳定性 a b s t r a b c t a 。b s t r a c t f o rt h es p a r s em a t r i xa ，t h ep r e c o n d i t i o n e rw h i c hg e n e r a t e sf r o mc o m p l e t e d e c o m p o s i t i o nu s u a l l yd o e sn o th a v et h es a m es p a r s ea st h em a t r i xa ，g e n e r a l l ym u c h d e n s e s o ，i no r d e rt og a i nt h et w oe n d st h a tp r e c o n d i t i o n e rk e e ps p a r s ea n dh a v e s u f f i c i e n te f f e c t , i n c o m p l e t ed e c o m p o s i t i o na p p r o a c hh a sb e e np r o p o s e d t h ep r e c o n d i t i o n e rt h a tc o m e sf r o mt h ei n c o m p l e t ec h o l e s k yf a c t o r i z a t i o nh a s b e e ns h o w nt ob eav e r ye f f e c t i v ep r e c o n d i t i o n e rf o raw i d ev a r i e t yo fl a r g es y s t e m so f l i n e a re q u a t i o n s t h ep r e c o n d i t i o n e rlw h i c hi sg e n e r a t e df o r mi n c o m p l e t ec h o l e s k y f a c t o r i z a t i o na l g o r i t h m3 - 4i nt h i sp a p e rr e t a i n st h en u m b e ro fn o n z e r oo fe a c hc o l u m n b e t w e e n n k a n dn k + p t h es p e c i f i cn u m b e ro fn o n z e r oe l e m e n t sw h i c hi s r e t a i n e d i ne a c hc o l u m ni sd e f i n e db yt h eo p t i m a lp a r a m e t e rf s o ，t h em e m o r yo f p r e c o n d i t i o n e ri sc o n s t r a i n e d f r o mt h en u m e r i c a le x p e r i m e n t , i tc a l lc o n c l u d et h a ti t s e e m st h a tt h eo p t i m a lp a r a m e t e rfi st h eb e s tw h e ni tc h o o s e st h ec o u n t d o w no ft h e o r d e ro fm a g n i t u d eo ft h ef r o b e n i u sn o r mo ft h em a t r i xa n dt h a tt h en u m b e ro fp c g i t e r a t i o no ft h en e wa l g o r i t h mi nt h i sp a p e ri ss a m ea st h en u m b e ro fp c g i t e r a t i o no f a l g o r i t h m3 - 2o f l i na n dm o r 6 ，b u tt h em e m o r yi sl e s st h a na l g o r i t h m3 - 2 f r o mt h er e l a t i o nb e t w e e nt h em o m e n to fs p a r s ep a t t e r nsg e n e r a t e da n dt h e w h o l ep r o c e s so fi n c o m p l e t ef a c t o r i z a t i o n ，s p a r s ep a t t e r nsc a l lb ed i v i d e di n t ot w o c a s e s t h ef i r s tc a s ei st h es t a t i cs p a r s ep a t t e r nw h i c hi sd e t e r m i n e dp r i o rt ot h ep r o c e s s o ff a c t o r i z a t i o n ，t h es e c o n dc a s ei st h ed y n a m i cs p a r s ep a t t e r nw h i c hc a nn o tb e p r e d e t e r m i n e di sd e t e r m i n e dd u r i n gt h ef a c t o r i z a t i o n f o rt h es e c o n ds p a r s ep a t t e r n , t h i s p a p e rp r o v e dt h es t a b i l i t yw h e nt h em a t r i x aw h i c hi st ob eb r o k e nd o w ni sa s y m m e t r i c a lm m a t r i x a tt h es a m et i m e ，f r o mt h ep r o c e s so fp r o o f w ec a i lf i n d 仃l 如 a tt h es a m ep a r a m e t e rpa n dt h es a m en u m b e ro fp c gi t e r a t i o n ，a l g o r i t h m3 - 4i sa t l e a s ta ss t a b l ea st h ea l g o r i t h m3 - 2o fl i na n dm o r 6 k e yw o r d s ：p r e c o n d t i o n e r , i n c o m p l e t ec h o l e s k yf a c t o r i z a t i o n ，m e m o r y , d r o ps t r a t e g i e s ， s t a b i l i t y i i i 独创性声明 本人声明所呈交的学位论文是本人在导师指导下进行的研究工 作及取得的研究成果。据我所知，除了文中特别加以标注和致谢的地 方外，论文中不包含其他人已经发表或撰写过的研究成果，也不包含 为获得电子科技大学或其它教育机构的学位或证书而使用过的材料。 与我一同工作的同志对本研究所做的任何贡献均已在论文中作了明 确的说明并表示谢意。 签名： ! 蜀立! 璺日期：沙。墨年r 月“日 关于论文使用授权的说明 本学位论文作者完全了解电子科技大学有关保留、使用学位论文 的规定，有权保留并向国家有关部门或机构送交论文的复印件和磁 盘，允许论文被查阅和借阅。本人授权电子科技大学可以将学位论文 的全部或部分内容编入有关数据库进行检索，可以采用影印、缩印或 扫描等复制手段保存、汇编学位论文。 ( 保密的学位论文在解密后应遵守此规定) 签名：困立! 璺 第一章弓| 畜 1 1 研究背景与意义 第一章引言 自从1 9 4 6 年第一台电子计算机问世以来，经过半个多世纪的发展，科学与工 程计算已经成为本世纪最重要的科学进步之一。科学计算己与理论研究及科学试 验并列成为当今世界科学活动的三耪主要方式之一。在许多科学与工程领域逡如 果没有计算就不可能有第一流的研究成果。为众多的科学与工程问题提供计算方 法，提高计算的可靠性、有效性和精确性，便是科学与工程计算这一领域的主要 研究内容。 数值线性代数又称矩阵计算，它是科学与工程计算的核心。数僵线性代数研 究的主要内容是，如何针对各类科学与工程问题所提出的矩阵计算问题的特点， 设计出相应的快速可靠的算法。可以毫不夸张地讲，大部分科学与工程问题最终 都要癌维力一个矩阵计算阅遂，其中具有挑战性的润题是大规模矩阵计算阕题。 事实上，许多大规模矩阵是稀疏的而不是稠密的，因此稀疏线性代数方程组的离 效求解是计算数学中十分重要的课题之一。稀疏矩阵的定义如下：如果矩阵 a r ”被称l 乍是稀疏的，那么它的任何一行中的非零元素的个数。，n 或者任 何一列中的非零元素的个数气。，m 。 科学与工程计算的许多重要领域，如计算流体力学、电磁场计算、材料模拟 与设计、电力系统的优化设计、石漓地震数据处理、数值天气预报及核爆炸数值 模拟都离不开微分方程的数值求解，焉这通常可通过有限体积、有限元或有限差 分法等方法进行离散，化为墨仁线性方程组或稀疏线性方程组，非线性方程组的求 解又可通过n e w t o n 迭代或简单迭代等方法进行线性化，最后也化为稀疏线性代数 方程组的求解闻题。大量经验表骧，在以上物理问题的数值模拟中，解线性代数 方程组所需要的时间在整个问题中的总计算时闻中往往占很大比重，有时甚至达 到百分之八十，成为整个计算中的瓶颈，从而探讨其高效解法十分重要。 线性代数方程组的解法可以分为蛊接解法( d i r e c tm 鼓h o d ) 帮迭代解法 ( i t e r a t i v em e t h o d ) 两大类。赢接解法具有在不计舍入误差的情况下麓褥到准确解 的优点，它还具有可靠性的优点。由于直接解法的这些优点，对于中、小型方程 电子科技大学硕士学位论文 组( 其阶数小于一千) 及其某些特定类型的大型的稀疏方程缝，常用此方法求解。 但当系数矩阵的条件数很大时，由于舍入误差的影响，所求出的解常与准确解相 差甚远，而且直接法的内存需求与计算时闻均很大。随着计算技术的发展，计算 机的存储量日益增大，计算速度也迅速提高，直接法在计算机上求解的线性方程 组的规模也越来越大，但直接法大多数均嚣对系数矩阵么进行分解，因而一般不 能保持么的稀疏性，因此求解起来比较困难。而在实际应用中，常常遇到大型稀 疏线性方程组( 其阶数是几百万阶甚至上亿阶) 的求解问题，这时一般采用迭代 法求解线性方程组。用迭代法求解线性方程组的优点是方法简单、易于控制，便 于编程，而且存储需求与每步迭代的计算量都很少，尤其是迭代法还具有能充分 利焉稀疏性的特性。儇必须选取合适的迭代格式及初始向量，以使迭代过程尽快 地收敛。与迭代法只知道系数矩阵与向量乘积的计算法则，而不必知道具体系数 矩阵就能求解相应线性方程组的能力来说，直接法更是无法比拟的。 迭代法种类繁多，但其收敛速度最终与矩阵的谱分布有关。近来研究最多的 迭代法是k r y l o v 子空间方法 1 】，该方法的收敛速度依赖于系数矩阵的谱分布。谱 的分布越集中，收敛速度越快，及之，收敛速度越慢，甚至不收敛 2 - 4 。为了改进 k r y l o v 子空间方法的可执行性和稳定性，便引入了预处理技术。后来，人们都普 遍意识到预处瑷技术是解决科学计算孛具有挑战性闯题的最有力的方法。 确实，在最近的不多几年中，人们找到了许多高效的预条件子。而且与崴接 方法和i , x y l o v 予空间方法相比，预条件技术现在已经成为一个比较活跃的研究方 向。但一个最好的、具有一般性的预条件予看起来是不太可能存在的，而且这种 状况在将来也是不会发生改变的。 1 。2 研究现状 在十九世纪期间，g a u s s 、j a c o b i 、s e i d e l 和n e k r a s o v 等人第一次用迭代的方 法来求解线性方程组。在二十毽纪前半叶，迭代法得到了重要的发展。僵是系统 的研究用来解大型线性方程组的迭代法仅仅开始于数字电子计算机发展之后。从 二十世纪开始五十年代开始，我们可以把迭代法的发展分为嚣个主要的历史时期。 第一个时期是以稳定迭代法为主的时期，大致是从上世纪五十年代到上世纪七十 年代前期。例如j a c o b i 方法、g a u s s - s e i d e l 方法、s o r 方法以及这些方法的改进与 加速形式。但随著科学技术的发展，需要求解的线性方程组规模越来越大，矩阵 2 筹一牵弓l 言 的类型越来越广，此时这些方法收敛速度一般太慢，甚至不收敛。其中最有效的 s o r 方法依赖于参数的选取，但一般无法得到最优参数，从而限制了该方法。第 二个时期是以k 坷l o v 子空阆迭代法为主的时期，这个时期开始予上世纪七十年代 中期。k r y l o v 子空间方法县有存储量小，计算量小且易于并行等优点，非常适合 于并行求解大型稀疏线性方程组。 1 2 。 k r y io v 子空间 假定n 阶线性代数方程组 a x = b ， 其中彳是大型稀疏非奇异方阵，b 是n 维列向量。设初值，定义r 0 = 6 一魄，迭 代m 步产生子空间 邑( 么，r o ) - s p a n ( t o ，么，a ”q r o ) 和近似解x o + k ，并且相应的残向量与另外某个一维子空间岛正交，即 毛= 刍一瓴土k ， 则称疋为k r y l o v 子空间，且上述方法称为k r y l o v 子空间方法。 依据匕不同的选取，可将k r y l o v 子空间方法分成三类【1 】【5 】【6 】。 第一类取乙= 疋，此时称为正交投影方法，h e s t e n e s 和s t i e f e l 与1 9 5 2 年提 出的c g 法为其中之一，此时要求a 对称正定。当a 不为对称正定矩阵时，f o m 1 j 方法也是正交投影方法。 第二类取乙= a 玩，这类方法具有残量极小性，郎名的e u c l i d 范数在子空闻 + 趸。中极小。g m r e s 7 1 【s 】为此类方法中的典型代表，该方法具有良好的数值稳 定性，但每步迭代的计算基与存储量都线性增长，同时难以有效并行化，从而提 戡了截断g m r e s 和重扇g m r e s 等方法作为其近似进行计算。 第三类取三埘= k ，卅tr o ) ，称为双正交化方法，此类方法用于a 为非对称矩阵 的情形。1 9 5 2 年由l a n c z o s 提出并于1 9 7 4 年由f l e t c h e r 改进的b i c g 法【l 】为此类 方法的经典算法，该算法是基于l a n c z o s 双正交化算法。此外，c g s 、b i c g s t a b 、 q m r 、t f q m r 等方法都是此类方法的演化与改进。b i c g 具有短迭代计算公式， 每步计算量和存储量不变，具有计算实现的高效性，健b i c g 方法要用到彳r 的运 电子科技大学硕士学位论文 算，著在许多情形下出现运算中断现象，同时收敛不稳定。1 9 8 9 年s o n n e v e l d 提 出的c g s 方法【9 】避免了使用么r 的运算，但该方法基于残向量多项式的平方，当 b i c g 运算中断时，c g s 也中断。为了克服c g s 方法的缺点，1 9 9 2 年v o r s t 提出 了b i c g s t a b 方法【l 们，采用了另种方法来构造残向量，其收敛性比c g s 有所改 善，同时不需要计算转置，但仍未克服运算中断的缺陷。1 9 9 1 年f r e u n d 和n a c h t i g a l 提出的q m r 方法【l l 】，具有局部残量极小，收敛性态较好，且不会中断，僵仍需使 用进行计算。为解决此问题，1 9 9 3 年f r e u n d 又提出了t f q m r 方法【l z 1 。 1 2 2 预条件 k r y l o v 子空间迭代法的收敛速度依赖予系数矩阵特征值的分布，对于很多问 题，比如特征值很分散时，直接使用k r y l o v 子空间迭代法的收敛速度特别慢，或者 根本不收敛。爵此使用预条件技术改变其收敛性，使中断问题可解，并且加速其 收敛速度。其基本思想是将原线性方程组转化为性质更好的线性方程组。为此， 需要对原线性方程组进行变换。设矩阵掰是矩阵蠢的非奇异近似逆矩阵( 在菜种意义 下的近似1 ，即 掰名。 那么线性方程组a x = b 可有三种基本变换形式。其一为左预条件法，即 m a x = m b ； 英二为右预条件法，即 a m y = b ，其中x = 坳； 除此之外，还有分裂型的预条件法 掰l 膨2 y = 掰l b ，其d ? x = m 2 y ，m = 掰l m 2 选择哪种类型的预条件子主要依赖于所选用的迭代方法，以及该问题的具体特殊 性等等吲。 一般来说，经过预处理的迭代将会比没经过预处理的迭代的计算时间和迭代 步数少的多。一个姆的顸条件子嬲应该满足下面的要求： 一、经过预处理以后的线性方程组系统应该容易求解，也就是说经过预处理 4 第一章弓| 言 后的迭代应该很快的收敛。 二、预条 孛子旋该容易构造和应用，也就是说保证每一次迭代都不会花费较 多的时间。 预条俘技术的关键问题是矩阵掰即预条件子的构造，基本的预条件有矩阵分 裂型预条件、多项式预条件、截断级数预条件、不完全分解预条件、稀疏近似逆 预条件等。 矩阵分裂型预条件【l 】是描对原矩阵进行分裂，分成两个矩阵和或差的形式，再 取其中一个比较接近a 的矩阵之逆作为预条件予。设爿的对角元素组成的对角矩 阵为d ，保留其严格上、下三焦部分中元素不变磊取其它元素全为零，分别记失 与【，则基于j a c o b i 的预条件子为肘= d ，该预条件又称对角预条件。类似地， 对应予g a u s s s e i d e l 迭代的预条件子为m = ( d 一三) ，对应于s o r 的预条件子为 m = ( d 一珊五) ，对应于s s o r 的预条件子为m = ( d o j u ) d ( d c o l ) 一。为了使 预条件更为有效，常把原系数矩阵进行适当地分块，然艨再进行预处理。 多项式预条件方法【5 】的基本思想是选择一个多项式s ( 力作为预条件子m ，使 得s ( x ) a 在某种意义下接近单位矩阵。早年r u t i s h a u s e r 、s t i e f e l 等人就曾用过多项 式方法，而近年来用多项式预条件是受并行机发展的影响，新的机器能高效地实 现矩阵运算。最有多项式s ( x ) 的选取，就是希望矩阵s ( 并m 能够尽可能的接近单位 阵。侧如，使s ( x ) a 的谱接近单位阵的谱，记么的谱为艿( 么) ，次数不超过k 的多项 式空间记为峨，则s ( x ) 饥求，使m a x 1 一以( 五) l 最小。其中e 是一个包含8 ( a ) 的 连续集合，实际上对么的谱只知道个大概，所以需要对的么谱有一个粗略的信计。 对于对称正定的矩阵4 ，可以取e 为一个包含爿的所有特征值的区间 口，】。 截断级数预条佟溺的思想是把系数矩阵a 的逆矩阵么焉级数表示，然后截断 级数的前m 项构成预条件阵。具体的做法是：假设矩阵a 分裂成两部分a p q ， 其中d e t ( p ) 雾0 ，记h = p - 1 q ，萁谱半径夕( 嚣) l ，则a _ 可雳级数( 置。) 芦表示， 取截断级数的前m 项构成预条件子m = p ( h 。) 一。 不完全分解预条件 1 4 - 1 9 】麓基本思想是把稀疏矩阵么分解戒一个上三焦矩阵汐 和一个下三角矩阵五，使得a = l 7 u ，分别取出它们中的某些元素形成新的上三 熊矩阵u 和下三角矩阵，使褥a l u ，那么矩阵l u 就是预条件子必。当稀疏 矩阵进行分解的时候，经常会发生元素填充( f i l l ，i n ) ，也就是说分解后产生的矩阵 和u 比原矩阵么的稀疏性大大降低了，相应地就会增加计算时间和存储空间。 电子科技火学硕士学位论文 毒予已有大量实验证明，不完全分解是串行计算是最有效的方法之一，所以一直 有大量研究在进行。按原矩阵的稀疏结构且不做任何填充的不完全分解技术效果 往往不佳，因此人们更多的是研究允许有更多填充的不完全分解技术。对于这种 技术，其核心是选择何种舍弃策略。基本策略有对某些特殊的位置进行元素填充 的，也有基于对分解因予中元素相对大小的舍弃策略，一般来说是大于某个阙值 的留下，而小予这个阂值的舍去。 稀疏近似逆预条件t 2 0 埘】是又一类近来研究最多的预条件方法，其基本思想是 构造系数矩阵逆矩阵的一个近似逆。讨论最多的构造方法是基于f r o b e n i u s 范数极 小化，即按极小化，一肼的f r o b e n i u s 范数来计算稀疏近似逆矩阵m 。如果直接计 算近似逆的话，英计算量缀大，故常常篱化为按鬈一a m , 的彝量2 范数极小化来求膨 的每个列，此处的岛是一个第i 个分量为l 、其它分壁为0 的列向量。该方法的核心 问题是如何确定膨中每列的稀疏结构。一种方法是事先静态地确定，但对一般的 稀疏矩阵，难以确定有效的结构。另一种也是研究最多的方法是动态地确定，该 方法中首选一初始稀疏结构，常选对角阵结构，再在计算中通过判断动态地变更。 动态方法中，常用方法是对每个么磁= 利餍迭代法求磁的近似值，丽在求近似值 的过程中采用舍弃策略。用上面的方法计算稀疏近似逆时，预条件予m 的形成时 闻褶对于不完全分解面言一般较长，为此爨现了分解型近似逆。构造方法主要有 两种，其中一种首先由k o l o t i l i n a 等提出，直接从矩阵分解的思想出发计算得到。 另一种构造方法是从不完全双共轭过程出发计算得到。该方法无需知道稀疏结构 且可应用于一般稀疏矩阵，其基本思想是如能计算如a 共轭的一对单位正交基底 乏：i = 1 ，2 ，n 与 心：f = l ，2 ，n ，且么非奇异，则知形r a z = d 为对角矩阵， 其中z 与渺分别是由列向量乏与m 组成的矩阵，从而a = z d q w f ，所以可望在 计算过程中采用适合的舍弃策略得到有效的稀疏近似逆。 1 3 本文的研究工作与论文结构 第二章主要介绍j a c o b i 、g a u s s + s e i d e l 、s o r 、a o r 、2 p p j 迭代法、文中将要用 劐的g c p 算法和不完全分解等预备知识。 第三章研究具有有限存储空间的带阈值的不完全c h o l e s k y 分解。众所周知，在 不完全分解中允许更多的元素填充能有效她减少迭代步数，但是允许更多的元素 填充进来就意味着构造的预条件子m 需要更多的存储空间，这对于上亿阶的大型 6 第一章引言 的稀疏矩阵来说可以说是一个挑战。在这里提出了一种可控制存储空间的算法3 4 ， 并做了数值实验，实验结果与l i n 矛1 m o r 6 r 7 】的算法3 2 的实验结果进行了比较，结 果表明在相同的参数和相同的p c g 迭代步数下该算法所需的存储空间更少，而且 该算法的不完全c h o l e s k y 分解更稳定。 第四章主要内容是研究了不完全c h o l e s k y 分解的存在性与稳定性。从稀疏模式 s 产生的时刻与不完全分解的整个过程的关系来说，稀疏模式可以分为两种情况， 第一种稀疏模式是在分解过程之前事先确定好了的静态稀疏模式，第二种是在分 解过程中产生的不可预测的动态稀疏模式。第四章则证明了在被分解矩阵么是对 称m 一矩阵的情况下第二种动态稀疏模式的稳定性，由该定理也可以看出新算法 的不完全c h o l e s k y 分解至少和l i n 和m o r 6 的算法的不完全c h o l e s k y 分解一样稳定。 7 电子科技大学硕士学位论文 2 。 迭代法 第二章预备知识 弟一早 耿宙刘状 迭代法一般可以表述为： x = 纸( x 一n ，z 。一) ， k = l ，z + l ， ( 2 1 ) 其中饵是从r ”r ”r “到r “上的算予，称为迭代算子，x 们，x 卜1 为迭代初 筐，可以事先给定，通常称形如( 2 1 ) 的迭代法为z 步迭代法，当z = | 时，就是 单步迭代法，如果迭代算子纭与k 无关，即纯兰妒，则称迭代法( 2 一1 ) 为定常迭 代，否则称为不定常迭代。这里将着重讨论单步定常线性迭代法： 工( 。) = 溉( + c k = l ，2 ， 其中b r “”为迭代矩阵，x f 。称为初值，e 为常向量。 对于线性方程组 a x = b 来说，把系数矩阵分解成 彳= m 一 其中坛是非奇具的，然蜃代入( 2 3 ) 式可褥 m x = n x + b ， 对上式两边同乘膨。可得 x ：m 一1 n x 十肘一1 及 给定一个初始向量可得迭代格式 x + 1 = 黜2 ) + e k = 0 ，l ，2 ，3 ， ( 2 2 ) ( 2 _ 3 ) ( 2 4 ) 第二章预备知识 其中b 称为迭代矩阵。迭代法( 2 - 4 ) 收敛的充分必要条件是 l i m b o = 0 t _ 羞p ( b ) l ，则称该迭代矩阵收敛，否则，该迭代矩阵发散。若产生的迭代序 列 x ) 收敛到一个确定的向量7 ，那么x 就是原方程组的解。 令a - d - l - u ，其中d = d i a g ( a l l ，8 2 2 ，a n n ) ， l = o a 2 l 0 一a 3 i a 3 2 0 - a n l一a n 2一a n 3 ，u = 0 一a 1 2 一a 1 3 一a l 0 一a 2 3 一a 2 * - a 一1 月 o 则可得j a c o b i 迭代格式 x 村) = d 一1 ( 三+ u ) x 。+ d b k = 0 ，1 ，2 ，3 其中鹭= 移一) 称为j a c 曲i 迭代矩降。 在计算霹“”时，如果用“n ，搿 代替霹柚，础，则可能会得到更好的 收敛效果此时的迭代公式为 解得 x ( + 1 = d 一1 ( 6 + 厶卅十乙k ) ， 菇+ = ( d 一三) 一u x 。+ ( d 一) bk = o ，1 ，2 ，3 此迭代格式称为g a u s s _ s e i d e l 迭代，简称g s j , 塞代，b c = ( d 一三) 一u 称为g s 迭代矩阵。 在g s 迭代中，为了得到更好的迭代效果，可在磐正项前乘以一个参数国，于 是就得到所谓的逐次超松弛迭法，简称s o r 迭代法，其中国称为松弛因子。此时 解褥 。+ 1 = 并+ c o d 一( 办+ 三。+ u x 瓣一d x 辑) ， x + 1 。( d c o l ) 一1 ( 1 一甜) d + c o u x + c o ( d o j l ) 一b ， 9 电子科技大学硕士学位论文 b s = ( d 一缈) 一1 ( 1 一c o ) d + u g n s o r s k 代9 l i 阵- t 2 5 1 。 s o r 迭代法只引入了一个参数缈，为了得到更好的效果，于是人们又引入了一 个参数，得到了么的a o r 迭代矩阵 0 ，。= ( d y 三) - 1 【( 1 - ( o ) d + ( 缈一7 ) 三+ 缈u 】， 可以看出，当y = 缈时，a o r 迭代法即为s o r 迭代法；当y = f _ o = 1 时，a o r 迭代法 即为g s 迭代法；当y = 0 ，c o = 1 时，a o r 迭代法即为j a c o b i 迭代法。 1 9 8 3 年，m i s s i r l i s 提出了并行j a c o b i 型方法，并讨论了它的收敛性，其中假 设了j a c o b i 矩阵，的特征值均为实数，所得的收敛条件还很复杂。我国学者胡家 赣将这个方法推广到了两个参数的情形，称为两参数并行j a c o b i 型方法，记为2 p p j 方法【2 6 l 。 对方程组( 2 3 ) 的迭代格式 石( “) = x 似j t o m 。1 ( a x ( 一6 ) ， k = o ，1 ， ( 2 5 ) 一般要取m 尽量接近于彳以取得较好的收敛性。仍设d = d i a g a ，c = d a ， j = d c 。若p ( j ) 1 ，则有 a 一= ( ，一玎d1 = ( 7 ) d - l ， r = 0 k 因此，可取m 一为( ，7 ) d 一，此处尼为一整数。现若取后= l 并引进另一参数，亦 r = 0 即取 m = ( ，+ r j ) d 。 代入到( 2 5 ) ，即得2 p p j 格式 z 似) = 石似一c o ( i + r j ) d 。1 ( 止似一6 ) ，k = 0 ，l 迭代矩阵为 i c o ( + r j ) d a = i 一缈( ，+ ) ( ，一，) = ( 1 - 0 ) ) 1 + 缈正， 其中 1 n 第二章预备知识 当缈= f 时，它就是1 p p j 格式。 2 2p c g 算法 i = ( 1 一r ) i + r j j 在讨论p c g 算法之前【2 7 1 ，我们先了解一下c g 算法。给了初始迭代值z ( 们，可 得r o = f a x 们，取p 1 = ，o 即最速下降方向，然后，对k = 1 ，2 ，3 ，可顺次得 出 口：芷型竺：丝! 丝! ：( r ( k - 1 ) , r ( k - 1 ) ) “( ，却) ( p ，和。) ( p ，a p ) 石( ) = x ( 一1 ) + 口 p ( ， 厂= 厂一厶= r 一吼却， ( 2 6 ) 如= 芳筹= 芳为， p 川= + 反+ l p n ( 2 - 6 ) 称为c g 算法的迭代公式。 c g 算法的收敛性 x ( k ) - - x * 卅p m 。虬2 ( 糕r 其中口和为彳的最小和最大的特征值，也可表示为 ，。忆州，虬2 ( 镑r c a = p a 此处白为的彳谱条件数。由上式可知，c 爿愈小，即愈接近1 ，c g 法收敛愈快。实 际上上式左端为第七次迭代的误差和初始按彳+ 莫l l x l l 爿= ( 么x ，x ) v 2 的比值。令上式右 边小于占，易知若要使这个比值小于占，即误差压- n 到g 倍，约需k = i n t ( 寺巳1 1 1 兰) 1 r 次迭代，此处i n t 为整数部分。 电子科技大学硕士学位论文 因为迭代法的收敛速度都与迭代矩阵的条件数( 系数矩阵的最大特征值与最 小特征值之比) 有关，条件数越大，收敛速度越慢。虽然共轭梯度法在理论上最 多经过n 次迭代就可以达到精度解，但实际上由于计算机舍入误差的存在和矩阵彳 的一些病态特性， 们，p c q 的a 正交性以及 ，1 ，似 的正交性随着 k 的增加而变差，故一般来说降低了收敛的速度而且x 不是精确解。对于大型和 超大型的线性方程组，即使n 次迭代收敛，也是实际计算中不能接受的。于是出现 了预处理共轭梯度法( p c g ) ，它是通过通过引入预处理矩阵m ，使系数矩阵的特 征值分布更为集中，降低矩阵条件数，改善矩阵病态特性，以达到提高收敛速度 的目的。 考虑方程( 2 3 ) ，其中为彳正定对称矩阵，又设m = 厦为么的一个近似分解 m = 口兰a 。这样，可用方程组 代替( 2 3 ) ，或者表示成 此处 事实上，粗略估计可知 ( f 1 肚r ) ( f x ) = 1 7 1 f f y = g ， f = l - 1 么f r ，y = r x ，g = 1 2 1 f ， f 兰l - 1 ( l l r ) l - r = ( l - 1 l ) ( l r l - 7 ) = j ( 2 7 ) 从而当口的确很近似于彳的完全分解时，f 近似于，。f 的条件数近似于条件数 的最小值l ，这里用的是，的谱条件数，即f 的最大特征值和最小特征值之比，它 是等于或大于1 的数。然后将c g 法公式( 2 - 6 ) 用之于方程组( 2 7 ) ，即可得i c c g 方法的迭代公式。 上面的i c c g 方法是在彳为对称正定的假设下来推导的。在实际问题中彳往往 不是正定甚至是不对称的。怎样解这样的方程组自然会引起许多学者的兴趣，从 而也提出了一些算法，比如i l u c g 方法和t c g 方法。 i c c g 方法用的是c h o l e s l ( 、r 分解来降低系数矩阵么的条件，对于彳是不对称 正定的话，就用l u 分解来降低a 的条件数。爿不对称( 更不为正定) ，则彳有近 1 2 第二章预备知识 似分解a 兰l u ，我们取f = 1 7 1 a u 一，则当l u 近似于么时，f 近似于j ，其条件 数就近似于l 。我们也有方程组 f y = (2ey g 2 8 ) = ， l 。芍) 此时，f = 1 7 a u ，y = u x ，g = f 1 厂。然后将c g 法公式用之于方程组( 2 8 ) ， 可得i l u c g 方法。 对于彳是不对称正定的情形，一般认为可采用i c c g 方法解( 2 3 ) 的法方程 a r a x = a r b ，但由于a 丁a 的条件数增大后而稀疏性减弱，故不宜采用。但实际上 并不需要实际求出a r a 后再用i c c g 方法。另一方面，i l u c g 方法是将( 2 3 ) 预条件 后得至l j ( 2 8 ) 。然后将c g 法公式用之于方程组( 2 8 ) 。我们也可对( 2 8 ) 的法方程 f r f y = f r g ( 2 - 9 ) 用c g 方法。若f = a ( l u ) ，则f 7 f = ( l u ) 可a r a ( l u ) = u 可1 2 r a r a u 一1 7 1 ， y = ( l u ) x ，g = b 。同样，我们也可取f = ( l u ) 叫a 或f = 1 2 1 a u 一得到方程组( 2 9 ) ， 然后用c g 法公式( 2 6 ) 解方程组( 2 9 ) ，这些迭代法称之为t c g 迭代法。 2 3 不完全分解 不完全分解预条件技术来源于对系数矩阵么的完全分解。19 7 7 年m e i j e r i n k 和 v o r s t 首先提出了不完全l u 即i l u 技术【1 4 1 ，由于m = u ) 一，所以此时在迭代法 中形如z = m r 的计算可通过求解上下三角线性方程组得以进行。 假定一个聆阶系数矩阵么，它的元素为( f ，j = l ，z ) ，l * nu 分别为稀疏下 三角矩阵和上三角矩阵，如果 a = l u ， ( 2 1 0 ) 那么就称( 2 1 0 ) 式为4 的完全l u 分解。如果彳= 三矽+ 尺，那么它就为a 的不完 全l u 分解。一般来说，不完全l u 分解过程计算出了一个稀疏下三角矩阵三和一 个上三角稀疏矩阵痧，并且误差矩阵尺= 三痧一彳满足一定的条件限制，例如某些 特定位置的元素值为零。 对任意一个非零模式集s 电子科技大学硕士学位论文 s c 肛j ) l i 硝f ，z ， 不完全l u 分解都可以表示成下面的算法。根据指标f ，j ，k 不同的循环顺序，可产 生31 = 6 种不同的算法，这六种可能性分别是彬，咖，匆，驰，声，危，下面我们只介绍 协，i l 萝，j k 三种算法，剩下的三种算法类似。 算法2 1 协型i l u 图2 1 对于算法2 1 来说，第k 步计算需要用到从第k + 1 行到第n 行和第k + 1 列到第 刀列的元素。如图2 1 所示，在第k 步，黑色部分是己计算好了的元素，阴影部分 是未计算好的元素。左下角的黑色部分是不完全l u 分解中的下三角矩阵中的一 部分元素，右上角的黑色部分是不完全l u 分解中的上三角矩阵u 中的一部分元 素。在实际计算中，算法2 1 执行起来比较繁琐，因为在第k 步，从第k + 1 行到第 n 行的元素都要变动。一般使用下面的两种算法，这两种算法分别是按行和列进行 计算的。 算法2 - 2i k j 型i l u 图2 2 算法2 2 是按行进行计算的，在第i 步，计算的就是下三角矩阵和上三角矩 1 4 第二章预餐知识 阵u 的第i 行。如图2 - 2 所示，阴影部分怒已计算好的元素，白色部分是未计算好 的元素，左下麓的阴影部分是下三角矩阵五，右上角的阴影部分是上三角矩阵u 。 算法2 3 是按列进行计算的，在第，步，计算的就是下三霜矩阵五和上三角矩阵秽 的第，列。与图2 2 一样，图复3 的阴影部分是已计算好的元素，白色部分是寒计 算好的元素，左下角的阴影部分是下三角矩阵，右上角的阴影部分是上三角矩阵 秽。 算法2 - 3 霸型i l u l ： | 隧2 ，3 按原矩阵a 的稀疏结构且不做任何填充的不完全l u 分解称为i l u ( 0 ) ，若爿是 对称躲，则称为无填充不完全c h o l e s k y 分解( 1 c ( 0 ) ) 。i l u ( 0 ) 和l c ( ) 由于不允 簪 任德填充瓶可能与完全l u 和完全c h o l e s k y 分解相差较大磊效果不健，所以许多 学者就希望通过允许更多元素的填充来改善i l 叭o ) 和i c ( o ) 的效果。 对予允许更多填充的不完全分解技术，其核心是选取何种舍弃策略。基本策 酪有两种： 一种是基于对矩阵结构的分析，如对五点差分形式的矩阵，如纂对称正定， 已有大量分析与实验表明，取( p ，g ) 为( 1 ，1 ) ，( 1 ，2 ) ，0 ，3 ) ，( 2 ，4 ) ，( 3 ，5 ) 时，i c ( p ，譬) 预条件效果穰好 2 8 - 3 0 l 。此外，在不完全l u 分解得到的预条件中，对比较规则的稀 疏矩阵，有类舍弃技术是基予对分解因子枣元素的层次值进行判断的( 经典的 有i l u ( p ) 算法) ，g u s t a f s s o n 萼f 入了i l u ( p ) 算法f 3 l 】，当填充元索的层次值大于某 个阂值p 时就禽弃掉，见算法2 ，4 。i l u ( p ) 的精确定义见定义2 。l 【l j 。 1 5 电子科技大学硕士学位论文 算法2 4i l u ( p ) 算法p 5i l u t ( 弘 定义2 。p 1 矩阵名中元素吩的裙始层次值为 魄= 罡 如果魄o 或者；= 矗 其它 1 6 第二章预备知识 每当不完全叫分解进行一次时，元素的层次值就会发生变化，层次值由下式确定 l e v u = r a i n l e v 口，l e v f , + 伽材+ 1 ) ， 从i l u ( p ) 的定义我们可以看出，当p = 0 时，不完全分解的非零模式集就是矩阵彳 的稀疏结构，i l u ( p ) 算法就是m u ( o ) 算法。 另一种舍弃策略是基于对分解因子中元素相对大小的判断。如在不完全l u 分 解中，将分解因子中的元素的绝对值与该元素对应的原系数矩阵中相应行的2 范 数与以参数f 的乘积值相比，当前者小于后者时就将该元素舍去。为限制存储量， 之后还可以只保留分解因子每行中p 个绝对值较大的元素，该方法即为经典的 i l u t ( p ，f ) 方法【1 9 】，见算法2 5 。 1 7 电子科技大学硕士学位论文 第三章具有有限存储空间的带阈值不完全c h o le s k y 分解 当矩阵么为对称矩阵时，l u 分解就可以看作是c h o l e s k