（应用数学专业论文）saoral预条件共轭梯度法.pdf

太原理工大学硕士学位论文 摘要 1 5 1 8 1 9 0 本文研究了解决大型稀疏对称正定线性方程组的类预条件共轭梯 度法。全文共分为四章。 第一章是对目前国内外研究现状的一个描述。 第二章提出了一种新的预条件共轭梯度法s a o r a l p c g 。这种 方法基于s a o r 迭代法以及交错法，构造了预条件子m ，然后利用共 轭梯度法来求解预条件方程m a x = m b 。 第三章对这种预条件共轭梯度法进行了分析推导出了它的条件数 比原来系数矩阵的条件数要低。 第四章用实例证明了这种预条件共轭梯度法的收敛速度比古典的迭 代法( 如j a c o 觇g s s o r ) 和传统的c g 以及s s o r p c g 要快一些 关键词： s a o r 迭代法，交错法，预条件共轭梯度法，条件数，线 性代数方程组，对称正定矩阵 2 a b s t r a c t 太原理工大学硕士学位论文 t h i sp a p e ri si n v e s t i g a t e dap r e c o n d i t i o n e dc o n j u g a t eg r a d i e n tm e t h o d i ns o l v i n gal i n e a ra l g e b r a i cs y s t e mo fl a r g es p a r s es y m m e t r i ca n dp o s i t i v e e q u a t i o n s t h e r ea r ef o u rc h a p t e r si nt h i sp a p e ra l t o g e t h e r i nc h a p t e rl ，t h eo v e r v i e wi sg i v e na b o u tt h es t u d ya tp r e s e n ti nt h e w o r l d i nc h a p t e r2 ，ip r o p o s ean e wt y p eo fp r e c o n d i t i o n e dc o n j u g a t eg r a d i e n t m e t h o d ，w h i c hi sc a l l e ds a o r a l p c gi nb r i e f t h ep r e c o n d i t i o n e rm i sd e r i v e do nt h eb a s eo fs a o ra n dt h ea l t e r n a t i n gm e t h o d t h e nw eu s e c o n j u g a t eg r a d i e n tm e t h o d t os o l v et h ep r e c o n d i t i o n e ds y s t e mm a x = m b i nc h a p t e r3 ，ia n a l y z et h i sa l g o r i t h ma n dd r a wac o n c l u s i o ni t sc o n d i t i o n n u m b e ri sl o w e rt h a nc g a l g o r i t h m s i nc h a p t e r4 ，s o m ee x a m p l e sa r eg i v e nt oi l l u s t r a t et h a tt h ec o n v e r g e n c e o fs a o r a 三p c gm e t h o di sb e t t e rt h a nt h ec l a s s i c a li t e r a t i v em e t h o d ( s u c h a sj a c o n ，g a u s s s e i d e l ，s o r ) a n dt h et r a d i t i o n a lc gm e t h o da sw e l la 8 s s o r p e g k e y w o r d s ：s a o ri t e r a t i v em e t h o d ，a l t e r n a t i n gm e t h o d s ，p r e c o n d i t i o n e dc o n j u g a t eg r a d i e n tm e t h o d ，c o n d i t i o nn u m b e r ，l i n e a ra l g e b r a i c s y s t e m ，s y m m e t r i ca n dp o s i t i v em a t r i x 4 太原理工大学硕士学位论文 符号说明 r m nm 行n 列实元素矩阵的全体 r n所有实”维列向量的全体，即r ” a r矩阵a 的转置 a h矩阵a 的共轭转置 a 一1矩阵a 的逆矩阵 一2 ( a )矩阵4 的谱条件数( n 2 ( a ) = l t a _ 1 忆h a h 2 ) l l 刮，a 的2 模，其值为( a ( a h a ) ) 1 2 a ：= ba 定义为b y ) z 和y 二向量的内积 a ( a ) 或a ：( 4 ) a 的某个特征值或第i 个特征值 i l l a x ；或m i n ；对所有可能的i 求最大( 或最小) a = ( o 玎)第i 行第j 列元素为n 巧的矩阵( 可为n n 阵、亦可 为m 阵等) a = d 一三一矿，其中d 为a 的对角阵，工是a 负的严格下三角 部分，u 是a 负的严格上三角部分，即 ，o l l、 。：f 嘞 l ，一l l l a n n 一u = 习 o a 2 1 0 a 3 【a 3 2 0 o n 【n n 2 0 k 爱 0 0 2 们0 o ，lil-ii、 太原理工大学硕士学位论文 a 为块矩阵，即 翔 5 2 2 2 i 2 n a a a 2 n a a a ，-illi、 ； a 第一章前言 在这一章给出我要研究的问题的基本情况及目前国内外研究状况 1 1 问题的基本情况 大型线性代数方程组的求解是计算数学中的重要课题之一科学和 工程的许多重要领域，如偏微分方程数值解，核物理与流体力学计算，结 构与非结构问题的有限元分析，电力系统的优化设计，石油地震数据处 理以及数值天气预报等，最终都要化成求解a x = b 这样的线性方程组 看似简单的求解问题，要找到消耗小且很逼近精确解的方法，却比较困 难所以国内外的众多学者为此作出了不懈的努力，并且得出了一些比 较好的算法正是因为这些原因，使得这门学科得到了迅速的发展，并 且得到了大家的重视 伴随着计算技术的发展，计算机的存储量日益增大，计算速度也迅 速提高，在计算机上可以求解的线性方程组的规模也越来越大通常求 解大型稀疏线性方程组a x = b 有两种途径：一种是直接法( 如g a u s s 消 去法，c h o l e s k y 分解法) 。它是通过对系数矩阵进行变换将原方程组化 为三角或三对角等容易求解的形式，然后通过回代或追赶等方法得到方 程组的解。 另一种是迭代法这类方法通过产生近似解序列 z ( 女) 来逼近精确 解矩阵a 只有在矩阵向量乘法时才用到，故计算曩相对小一些从而 迭代法的关键问题就是收敛性与收敛速度。不收敛的格式自然不能用， 那些虽然收敛，但收敛很慢的格式自然也不能用所以必须寻找收敛速 度较快的迭代格式。一般来说，一种迭代格式只能对某些类型的系数矩 阵有效，并不能适应任何种类系数矩阵的线性方程组。 研究某一些或某一类特殊系数矩阵的大型稀疏线性方程组的高效的 迭代方法，使之能达到比较理想的收敛效果已成为当前数值计算研究的 6 太原理工大学硕士学位论文 主要对象。 1 2国内外研究现状 7 对于一般大型稀疏线性方程组直接解法的研究，目前主要集中于三 对角、块三对角、带状线性方程组的并行求解迟利华 7 】对此进行了 比较全面的总结，并且基于并行环境提出了改进后的新算法如d p p l 及 p p d 2 等等，同时也通过实验说明了这些新算法比原来的效果要好 直接法在求解过程中容易带来大量的非零元素，增加了计算量，同 时不易并行，故而不能满足求解大规模问题的需要，因此通常使用迭代 法。迭代法包括j a c o b i ，g a u s s s e i d e l ( 简记为g s ) ，s o r ，s s o r 等方 法，y o n g 9 和v a r g e 【1 对此进行了细致的描述与探讨他们的著作已 成为迭代法方面的经典著作但随着科学技术的发展，有越来越多的特 殊的、大型的、稀疏的问题出现，仍用这些方法已显得收敛太慢( 甚至 不收敛) 。而且s o r ，s s o r 等方法往往依赖于最优参数的选取，而一般 方程组无法直接得到最优参数，从而限制了这些方法的高效使用这些 方法现已很少直接用于求解大型稀疏线性方程组 目前众多学者的注意力集中于如下两个方面：其一是对迭代法的改 进。如文【1 3 1 和文 1 5 】都是用单位矩阵与原矩阵a 的负的上三角矩阵之 和( 即+ u ) 乘以原方程组a x = b ，然后使用j a c o b i 法，或g s 法，或 s o r 法进行迭代，就得到了加速的迭代方法，并且用实例说明加速后的 方法比原来的方法收敛快文 2 0 提出用c h e 6 s s h e v 法和共轭梯度法来 加速s a o r 方法文 1 4 中分析了预条件a o r 方法，仍是用j + s 乘以 原方程组a x = b ，然后使用a o r 方法进行迭代其中s 为如下形式： s ：f0 0 _ 0 1 1 双向并行分裂算法 2 严格对角占优三对角系统的一种并行算法 或 站( i ! - - a l n o ： o 太原理工大学硕士学位论文 而文【3 3 】分析的修正的g s 法也用了类似的手段，只是s 形式不同，这 里 s = o l 1 2 0 0 - c 2 a 2 3 0 o 0 0 。0 1 i f o n 一1 a n 一】，ni o ) 所有这些方法实质上都是预条件技术，即先将原来的方程组加以变化然 后使用相应的迭代法通常使用左或右预条件子m 使方程组a x = b 变 成易于求解的形式 或 a m y = b x = m 一1 y( 1 1 ) f l a z = m 一1 6 ( 1 2 ) 文f 2 9 1 基于预条件技术对迭代法的多种加速技巧进行了总结与讨论，并 进一步提出了适合向量化和并行化的预条件处理法。其二是将迭代法作 为预条件子与其他方法结合使用。最为常见的就是和共轭梯度法结合使 用 共轭梯度法( c o n j u g a t eg r a d i e n t 简记为c g 法) 是5 0 年代初期由 h e s t e n s e 和s t i e f e l 首先提出的，它在理论上是一种直接方法，具有n 步 终止性，且具有短迭代计算公式，保证了计算实现的高效性。但实际上 由于计算中的舍入误差，导致其数值稳定性不好故而直到1 9 7 1 年r e i d 太原理工大学硕士学位论文 9 首次将c g 法作为一种迭代法使用后，这一方法才得到广泛使用而c g 法的收敛性又依赖于系数矩阵的条件数 4 】，为了提高收敛速度就要设 法改善c g 算法f 2 6 j 或矩阵的条件数。同迭代法一样，改善条件数的最 重要的方法也是预处理技术，即构造预条件子m ，使a m - 1 或m _ 1 a 尽可能接近于单位阵。到目前为止，对大型稀疏对称正定方程组人们所 研究的预条件技术可分为以下三类： 1 不完全c h o l e s k y 技术 最重要的预处理方式之一涉及计算a 的不完全c h o l e s k y 分解 其思想是计算一下三角阵l ，它具有可利用的稀疏结构，且“接近 于”c h o l e s k y 因子g ，即a = g g t = z 三t + r 预条件子m 则变为 m ；三z 了1 通过求解三t z = p 得z = m 一1 p 亦即，首先从孔= p 解出u ，然后从z 丁z = u 解出2 注意到这些求解的步骤是简单的 回代。因此，不需要显式得计算m 或其因子的逆。 依据原矩阵a 的稀疏模式且不做任何填充( 即不增加非零元的 个数) 的i l u 分解称为i l u ( o ) 若a 是对称的则称为无填充的不完 全c h o l e s k y ( i c ( o ) ) 分解，即i c 是i l u 的特殊情况若还是正定的 则用，g ( o ) 预条件g 即得i c c g ( o ) 27 j 。由于i c ( o ) 不允许任何填 充而与完全c h o l e s k y 分解可能相差较大而效果不佳。故许多作者， 对i l u 和j g 分解进行了大量的研究和改进，如文 2 8 ， 3 0 等，提 出了许多更有效的i l u 和，c 预条件子，如对比较规则的稀疏矩阵 的多水平填充的i l u ( k ) ，使预条件子m 的行和与原矩阵a 的行和 相等的修正的i l u ，块i l u ( o 口b i l u ) 以及基于多层技术的儿矿m ， b l u m 等等。 2 基于矩阵分裂的多项式预条件技术 这一技术由j o h n o n ，m i c h e l l e 和p a u l 1 6 于1 9 8 3 年提出。其思 想基础是既然预条件子凹实质上是系数矩阵4 的某种近似，因此 我们可将方程组m z = r 的解。看作是a z = r 的近似解。求a z = r - 近似解的一种有效的方法就是利用矩阵分裂，记g = m 1 n z ，那么 1 0 就取 z = 知= ( j + g + + g p - t ) m i 一1 r 太原理工大学硕士学位论文 ( 1 3 ) 作为其近似解，从而就可取 m 1 = ( j + g + + g p 一1 ) m l 一1( 1 4 ) 作为预条件子。即用如下公式迭代p 步： m 1 z 女+ 1 = n i z k + 7 ，z o = 0( 1 5 ) 其中假定4 = m 。一1 是a 的一种较好的分裂因为它是矩阵g 的 多项式，因此称之为多项式预条件这类预条件子从向量并行的 观点看很有吸引力，因而引起了广泛的重视。如文 1 7 1 8 1 5 具体给 出了利用前面提到的一些迭代法作为多项式预条件子，如j a c o b i 、 s g s ( 对称的g s 法) 、s s o r ，实质上是选取这些迭代法的迭代矩 阵作为上面提到的g ，然后形成预条件子m 当然这里的m 应当 是对称正定的，这就限制了公式( 1 4 ) 中m 1 ，l 和p 的选取通常 计算中选取的p 值都比较小( 如p = l 一3 ，否则计算量很可观) 于 是许多作者的注意力就转向了p = 1 的特殊情况。如文 4 】就a 的多 种分裂给出了选取m 的比较全面的总结基于块分裂文【3 1 1 提出了 比块s s o r 预条件共轭梯度法( 简记为b s s o t t p c g ) 更优的不完 全分解的块s s o r 预条件共轭梯度法( 即i b s s o r p c g ) 3 结合实际问题用区域分解作预条件子 6 j 区域分解法是并行计算驱动的，但实践证明它也可用于构造整体 预条件子，这通常是对由离散偏微分方程组( 简记为p d e 产生的线 性代数方程组来进行的其基本思想是：将一个给定的区域分解为子 区域，并在每个子区域上计算解的一个近似。c h a n 和g o o v a e r t s 【3 2 】 证明了区域分解法可实质性地改进收敛速度，至少当子区域个数不 太多时是这样。区域分解常结合迭代法使用，且通常内边界条件( 即 对对角块的校正) 基于相邻予区域上近似锯的信息，它一般由前一个 迭代步得到。 太原理工大学硕士学位论文 r o d i c o t i t 和r o b e r tf 3 4 1 对一个矩形网格利用区域分解法，其中 每个子区域内部执行一步s s o r ，为了进一步改进并行性能，其中 表达式中的逆用低阶截断n e u m a n n 级数来近似代替。在一个给定 矩阵a 的交叉对角内计算t l u 因子。当应用预条件子于一个向量v 时，重叠区域上的值取做用重叠i l u 因子计算的两个值的平均。 由于预条件共轭梯度法在解决实际问题中的高效性，将预条件共轭 梯度法稍做变化后用于处理非对称线性方程组的情形也引起了众多学者 的兴趣，有不少的文章对此进行了讨论，如 4 】中分析的i l u c g 和t c g 算法建立在这两种算法的基础上，f l o 对用迭代法作预条件子进行了 总结， 1 9 1 用修正的不完全三u 分解作预条件， 2 l ，2 2 ，2 4 ，2 3 】通过构造 各种稀疏近似逆作为预条件子。 总之，预条件共轭梯度法仍是我们研究的一个热点问题。 1 3本文的研究内容 在工程领域，如控制，网络等，许多实际问题都归结为解一个系数 矩阵4 为对称正定矩阵的大型稀疏线性方程组 a 茹= b ( 1 6 ) 目前，与各种预条件子结合的共轭梯度法是解决这类问题的主要方法。 本文提出了一种新的预条件共轭梯度法主要讨论如下两个问题： 1 以s a o r 迭代法和交错法为基础，提出了一种新的预条件子m ，并 且结合一般的共轭梯度法，导出了s a o r 一4 三一p c g 方法在第四章 通过实例计算说明了s a o r a 三一p c g 是一种收敛性比较好的方法。 2 对s a o ra l p c g 算法进行了分析。主要是预估系数矩阵的条件 数最后通过实例计算说明s a o ra 三一p c g 方法比c g 方法， s s o rp c g 方法收敛性都好。 第二章s a o r - a l 预条件共轭梯度法 本章介绍求解大型稀疏对称正定方程组a x = b 的预条件共轭梯度 法( p r e c o n d i t i 。n e dc o n j u g a t eg r a d i e n tm e t h o d s ，简记为p c g 方法) 的 一些基本知识及其算法，然后对如何构造预条件子s a o r a l 做一详细 的说明 2 1 几个重要概念 本节将对后面要涉及的几个重要概念做一下介绍，如s a o r 迭代 法、交错法( a l t e r n a t i n gm e t h o d ) 。 2 1 1s a o r 迭代法 在解线性代数方程组( 1 6 ) 时 和， a = m 一 且用迭代 常将系数矩阵4 分裂成两个矩阵m ( 2 1 ) m x + l = n x a + b ( 2 2 ) ，来解线性代数方程组( 1 6 ) ，这叫做解线性代数方程组的分裂法许多 迭代法都可以写成这样的形式。 设“，y r ，0 ，把系数矩阵a 写成其分解形式a = d 一三一u ， 则快速超松弛( a c c e l e r a t e do v e r r e l a x a t i o n ，简记为a o r ) 迭代方法 1 1 】 可写为 ( d 一7 l ) x k + l = ( ( 1 一u ) d + ( u 一7 ) l + w u ) z 女+ w b( 2 3 ) 这里和以后为避免繁琐起见，我们省略了“k = 0 ，1 ， ” 太原理工大学硕士学位论文 此时分裂法中的m ，分别为 相应的迭代矩阵为 1 3 ( 2 4 ) 妒= m 一1 n = ( d 一7 l ) 一1 ( ( 1 一w ) d + ( u 一1 ) l + w u ) ( 2 5 ) 若将迭代法中的l 和u 互换，则a o r 迭代方法为 ( d 一7 u ) z k + 1 = ( ( 1 一w ) d + ( u 一7 ) u + u l ) 。女+ w b ( 2 6 ) 此时分裂法中的m ，分别为 ( 2 7 ) 相应的迭代矩阵为 = m 一1 n = ( d 一7 u ) 一1 ( ( 1 一u ) d + ( u 一7 ) u + w l )( 2 8 ) 这两次a o r 迭代合起来就构成一次对称a o r 迭代，记为s a o r 相 应的迭代格式为 ( d - 7 l ) x k + j = ( ( 1 u ) d + ( u 一，y ) l + 叫矿) 。七+ 6 ( 2 9 ) ( d 一7 u ) z k + l = ( ( 1 一w ) d + ( u 一1 ) 矿+ w l ) x + ；+ w b 此时的迭代矩阵为 s = 西妒 显然，s 是s s o r 迭代矩阵的推广 矿u + 己 7 一 u + ， d 江 西 一 一 1u，u 1 | 1 【 m l +矿 7 一w+ 忉 加 一 一 淝 lu一u = | | m 1 4太原理工大学硕士学位论文 2 1 2 交错法 c h u a n l o n g 一$ v a n 9 和t i n g z h u h u a n 9 1 2 分析了如下的交错法。 交错法：对任意给定的初始向量z o ， + 2 。+ m 。b f 2 1 0 1 x k + 1 = 尸“q 。十；+ p 一1 b 这里a m n p q 是系数矩阵a 的两个分裂。设 r = p 一2 ( q m 一1 + ，1 ( 2 1 1 ) t = p “q 。n 那么交错法可以用如下简单的迭代形式等价的描述： x k + l = t x k + 硒( 2 1 2 ) 当a 是日矩阵或单调矩阵时，他们证明了交错法的收敛性且提出了单 调收敛理论，并且提出如下定理： 定理2 1 1 设a _ 1 兰0 ，a = m n p q 是正规分裂如果 ( a + 十q ) 一1 至0 ，那么由丁导出的唯一分裂a = b g 仍是正规分 裂。 定理中提及的口，e 分别为 b = p ( 肘+ 尸a ) 一。m ，g = n ( a + n + q ) 一1 q( 2 1 3 ) 2 2 共轭梯度法 自1 9 7 1 年r e i d 提出将g g 看成一种迭代法使用后，近3 0 年来，有 关的研究得到了前所未有的发展但由于c g 的收敛性依赖于系数矩阵 的条件数，为了提高收敛速度，就要设法改善矩阵的条件数，于是t t 预 条件”技术应用而生。随着各种预条件不断提出和完善，当前预条件+ 共轭梯度法已经成为求解大型稀疏对称正定方程组的最有效办法之一 太原理工大学硕士学位论文 1 5 p c g 方法是将方程组预条件后用c g 法求解，故对c g 法还应徽一些介 绍。在文 2 中对导出c g 做了详尽的阐述，这照我仅列出其具体算法。 葵法2 2 。1 拱轭檄度法p g 川 如果a r n “对称正定，b r ”，。o r “为初始近似俊口z o 6 ，则下列算汝可计算盘x r “，使得a x = b k 一0 r 。b4 z o w h i l er k 0 = + 1 缮= 1 p l = r 0 e l s e 咎k = r k _ l t r 一t 弦t _ 2 t 7 t 一2 p k = r k 一1 + p p k l e n d = = r 一l ，7 女】，p k t a p k x k 。z 女一1 + a k p k = 一1 一a k a p k e n d g 然x k 这个算法每步迭 弋只用到系数矩蓐a 徽一次短阵廊量莱法运算，赦 计算是糨娄篱罄熬 佟为本节的结尾，下蘧讨论c g 产生的迭代点列 茁女) 的收敛性。这 里给出两个结果矧，它们都说明了不管媳在低秩扰动的意义下逐是在 范数的意义下，当4 接近于单位阵时，c g 法非常有效。 定理2 2 。1 + 如巢a = i + b 为n 盼对称正定晦，鼠r a n k ( b ) 一r ，则算 法r 2 g ，1 ) 基多r 七1 旁收敛， 证明见袁亚湘译矩阵计算 ( 2 0 0 1 ，p 6 1 4 ) 从这一结果可得如下重要的结论： * 如果a 接近于单位阵的秩r 校正，则算法( 2 2 1 ) 经避k + l 步后差 不多就可收敛。 另一种形式的更常用的误差界是糟a 范数表示的定义a 范数： - 、兀了石。 1 6太原理工大学硕士学位论文 定理2 2 2 假设a r “对称正定，b r “，若 x k 为算法俾2 纠 产生的迭代序列，k = 一2 ( a ) ，则 忙咱”忪一忆( 筹) 证明见胡家赣著线性代数方程组的迭代解法 ( 1 9 9 1 p 1 7 8 1 8 1 ) 。 z k 的逼近程度经常比定理估计要好得多。然而，上述定理有一个 更实用的形式： * 若n 2 ( a ) z1 ，则在a 范数意义下，c g 法收敛很快 总之，当矩阵a 为良态或有少数几个不同的特征值时，c g 法用起 来很有效 2 3s a o r 一4 l 预条件共轭梯度法 这一节先介绍如何导出p c g 算法，然后是具体的预条件子的选取 2 3 1 推导 p c g 法的思想是对于对称正定线性方程组( 1 6 ) ，把正规的c g 法 运用到变换了的方程组： a 茅= b ( 2 1 4 ) 其中彳= c 一1 a c ，茁= c t z ，i = c - 1 b ，月为对称正定阵。如果把 算法( 2 2 1 ) 用到方程组( 2 1 4 ) ，则得到： 算法2 3 1 k = 0 z o 给定r a5 0 = 6 j 而= b a i o w h i l e 敢0 k = k + 1 太原理工大学硕士学位论文 玎k = 1 声l = 而 e l s e 凤= 墨，瓦一椎。死一z p k = r k l + 凤诹一l e n d o k = i o l 魂一1 磁e a c 。诹 i 女= i k 一1 + o k 诹 r k = r k 一1 一a k c 一1 a c 一7 诹 e n d i = i 这里瓢为苗的近似，氟是变换了坐标系的残量，即 张= b a 苗 当然，一旦得到窑的值，则通过方程x = cr l z 就可得到。的值。 我们也可以直接用r k 来计算解实际上 1 7 氏= i 一五巩= c 一1 ( 6 一a x 女) = c 一1 r k 若令p = c 。氩，则p l = c 一了而= c 。c t o ，再一次运用算法( 2 21 ) 则得到下面算法： 算法2 3 2 k = 0 给定初始值x o ( a x o b ) 7 0 = b a x o w h i l ec 一1 r k 0 k = k + 1 i f k = 1 1 8太原理工大学硕士学位论文 c p l = c 一”o e l s e 敝= ( g 一1 n 一1 ) 7 1 ( g 一1 “一1 ) ( g 一1 r 七一2 ) 丁( e r k 一2 ) p 女= c 一7 c 一1 “一l + 卢m l e n d n = ( c 一1 r 七一1 ) 丁( c 一1 7 一1 ) p k 了1 a p x k2x k 一1 十c o k p k r k 2r k 一1 一a k a p k e n d z = x k 本算法可避免直接利用矩阵c 但如果定义预条件子m = c 2 ( 也 是正定的) ，且令张为方程m z k = r k 的解，则上述算法可简写成： 算法2 3 3 疗员条件共轭梯度法f p c g 给定对称正定矩阵a r n “， b r “及对称正定的预条件子m ，和初值x o f a x o ”，则本算法可求 出方程组a x = b 的解。 k = 0 r o = b a x o w h i l er k 0 解方程m z k = “ k = k + 1 玎k = l pz = z o e | s e 8 k = r k _ 1 t z k q r k _ 2 t z k 一2 p k = z k l + 卢k p k 一1 太原理工大学硕士学位论文 e n d o k = r k _ 1 t z k 一1 p 丁a p k x k = x k 一1 + n p r k2r k 一1 一a k a p k e n d 1 9 z = z k 关于这一算法，应注意到虽然变换阵c 在算法的导出上起了很重要 的作用，但它只有通过预条件子m = c 2 起作用。为了使算法( 2 3 3 ) 在 处理稀疏矩阵时有效，线性方程组m z = r 必须很容易求解，且收敛速 度还要快。于是，设法选择一个好的预条件子就成为问题的关键所在 2 3 2s a o r a l 预条件子 文4 1 中证明了经过s s o r 预条件，方程组( 2 1 4 ) 的系数矩阵a = c a c o 的条件数等于原方程组系数矩阵4 的条件数的平方根，也就 是说用s s o r 预条件确实可以改善c g 的收敛速度因而这里考虑使用 更优的s a o r 预条件但由于s a o r 迭代矩阵s 的中间两项的乘积一 般不能交换，导致其不象s s o r 法那样容易确定对称分裂矩阵，故而将 s a o r 迭代法与交错法结合使用来寻找预条件子m 根据交错法，对a m np qb c 进行如下方式的分裂 选取：m ，由( 2 4 ) 给出，只q 分别对应于( 2 7 ) 中的m 和，即 而b e 就由( 2 1 3 ) 给出 7 ) l + w u ) ( 2 1 5 ) 7 ) u + w l ) u u + 十 q 归 归 叫 叫 一 一 一 一 d 姐 d n 一，一u，一u一u | | l | = = m p q 太原理工大学硕士学位论文 故而对于大型稀疏对称正定方程组( 1 6 ) 可选取m = b 作为预条件 子( 此时对于分裂a d l u 有u = 五了1 ) 其具体的形式如下： ，= 口= 尸( 肘4 - p a ) 一1 m = l d 一7 上t ) ( ( 2 一“) d 一( 7 一u ) ( l + 五丁) ) 一1 ( d 一7 l ) ( 2 1 6 ) 、 ：而1 ( d - 7 l t ) ( 。一7 2 型w w ( l + l t ) ) ( d - r l ) 2 3 3 一些实用的细节 前面给出的算法( 2 2 1 ) 和( 2 3 3 ) 的终止准则是不现实的虽然在理 论上已经证明了c g 法至多1 2 步就能得到方程组的精确解，由于舍入误 差的影响，残量间的正交性有所损失，从而在数学上不能保证有限步终 止。此外，在实际应用算法时，n 通常很大，使得o ( n ) 步迭代工作量仍 是不可接受的因此，用最大迭代次数。以及残量的范数来给出终止 准则即算法( 2 2 1 ) 和( 2 3 3 ) 中的r k 0 由 ( 、，佤 e 1 2 ) a ( k f 删f 2 ) a k m n 。) 解方程m z = r 的2 七= 七+ l 玎k = 1 p = o p 2 r t z 太原理工大学硕士学位论文 2 l e l s e 8 = p l 再 p = z + 8 p e n d “= a p n = p p 丁“ z = z + n p ? = r n “ 声= p p 2 r t z e n d 此外，对算法( 2 3 3 ) 中的求解m z = r 问题，也有许多学者予以了 关注，如文 2 5 提议仍用迭代法来解决问题在这儿我仍是采用经典的 做法直接进行求解。 第三章算法分析 这一部分主要是对系数矩阵a 为块三对角矩阵时，用s a o ra l 预 条件后系数矩阵的条件数比原系数矩阵的条件数要小进行了证明对非 三对角矩阵的情况今后还需从理论上进行深入研究 3 1 块三对角矩阵中预条件子m 的形式 设大型稀疏对称正定线性代数方程组的系数矩阵a 为块三对角矩 阵，即a = d l l 了1 ，其矩阵形式如下： a = d 1l 2 。 l 2d 2l 3 t l 3d 3 l n 一1d n 一1 l n 7 1 l nd n 其中d ：为m 阶方阵 令t = d 一；岩( 三+ l t ) ，将其表示成矩阵形式为 t= d 1 群三2 z 。：旦。l 2 下 d ， 措l 3 壬老如丁 d 3 考l 4 7 搿k 一1d 壬岩l 。丁 ：z z r - “! l n d n 显然丁为三对角形式。设j ；为m 阶单位阵，由p e 方法【4 】可得 太原理工大学硕士学位论文 此处 t - r = r = r r l n _ 五n l n & ) 品 ，刀一j s i l 珂、 i 厶j 。罐。工：l i ni r 当a 正定时& 均为正定，故可取 b = ( ；：? ：s 。) f s _ 5 墨 ， 簖； 2 3 s a o r a 三预条件子m k r f r ，。一，。 、， 碍又 一晶 l o 卜 dld 岩 + 矗 d d = = & & 整理可得 太原理工大学硕士学位论文 m = 硒b ( d - 7 l 7 ) ( 。一爿( l 彬) ) ( d - t l ) = 币b ( 。- 7 l t ) t 。( d 一倒 = 杀面( 。- v l r ) ( b b t ) 。( d 一7 工) = i 南( ( _ d 一，y 上t ) b t ) ( ( d 一7 l t ) b t ) 7 1 3 2 条件数分析 下面具体分析s a o r - a l 预条件后矩阵方程条件数的变化情况。为 计算方便起见，这里考虑去掉系数1 ( w ( 2 一) ) 后的预条件子m ，即 m = ( d v l t ) b 一丁) ( 日一1 ( d 一- y l ) ) 记 c = ( d 一7 l r ) b 一丁 c 丁= b 一1 f d 一7 l ) 故而所找的预条件子m = c c t 由( 2 1 4 ) 式，这里要估计彳= e 此处c 和c 丁由( 3 2 ) 给出 由 ( 3 1 ) ( 32 ) 1 a c r 的条件数畅( 或写为c o n d ( 五) ) k ：幽：( c - 1 a c - t x , x ) 篮( a c i 羔c ：二a y 型一a 们 ( 3 s 一( 一丁z ，一7 z )( ，f ) ( ，9 ) 一 i r 一2 ( c t y , c t y ) 2 赫 太原理二大学颐士学位论文 2 5 此处y = c 一了1 z 由( 2 1 6 ) 和( 3 1 ) 整理得到s a o r a l 预条件子 肚( 。- 7 l t ) ( 。措阶矿) ) ( d - t l ) = u ( 2 一) a + ( ( 1 一u ) d + ( 叫一1 ) 三+ w l 丁) ( d - 暑( 工+ ) ( ( 1 刊d 蛳刊三t - 4 - o j l ) = u ( 2 一u ) a + ( ( 1 一u ) d + ( u v ) l + w l t j b t 1 x ( b 1 ( 1 一w ) d + ( u 一7 ) l r 十u 三j ) 令r = ( ( 1 一u ) d + ( 。一7 ) + w l t ) b r 从而 m 2 ( 2 一w ) a + f f 丁 f 3 4 1 故 ( 玎笋) = “) ( 2 一) ( a p ，p ) + ( r 丁y ，r 丁可) “j ( 2 w ) ( a y ，y ) 当o “ 一 锄 叫。 +一十 舄 + 0 啪 一一 一一 太原理工大学硕士学位论文 由 得 由 得 安= - 一面1 + c 击一z 一志，p 瞄c ，刊2 = 击( 篙) 蕊o s = 舄+ 1 7 9 8 ( 7 1 - w ) 2 一品) p蕊2 万j 矛十一研p 螂c 7 叫2 = 南辔 由( 3 8 ) 和( 3 9 ) 可得到如下简化的结果 1 ：u 2 3 + 4 ( 2 - w ) 3 ( 1 一u ) = 去击 将( 3 1 0 ) ，( 3 1 1 ) 代入( 3 7 ) 式整理得s 的最小值 为 s m i 。= 1 + ( 2 w 一3 ) ( u 1 ) ( i p ) ( 38 ) ( 3 9 ) ( 3 1 0 ) ( 3 1 1 ) ( 3 1 2 ) 下面先求方程( 3 1 1 ) 的根令p ( 9 8 ( 1 一卢) ) = m ，则原方程( 3 儿) 变 4 7 u 3 + 1 8 w 2 2 0 w + 8 一m = 0 ( 3 1 3 ) 南 暑。 嘉 品。 太原理工大学硕士学位论文 引入参数t ，整理( 3 1 3 ) 式为 ( u 。一i 7 w + ；) 2 一( ( 萼- 1 8 + t ) u 2 + ( 一i 2 t + 。) u + 学- ( 8 - m ) ) = o ( 3 1 4 ) 要使方程( 3 t 4 ) 大括号内的二次三项式成为一个完全平方式，当且仅当 ( 一互7 t + z 。) 2 一a ( 挈一- s + c ) ( 学一c s m ，) = 。 c s s ， 即 t 3 1 8 t + ( 1 0 8 4 m ) t 一( 2 1 6 + 2 3 m ) = 0 ( 3 1 6 ) 令t = y + 6 ，则( 3 1 6 ) 式变为 3 + 4 m y + r n = 0 ( 3 1 7 ) 由t a r t a n 8 公式，可知方程( 3 1 7 ) 的任意一个解y o 可表示为： y o =刘一詈+ 譬+ ( 警) 3 + 一筹一警+ ( 等) 3 im ( - 1 + 葡1 6 2 + 洱习葡 ( 3 1 8 ) ( 仁学) _ ( h ) 扎( 弓t + y o + 2 6 ) u 慨 一y 0 2 + _ 1 2 y o f 1 + m ) ：o 太原理工大学硕士学位论文 中 2 9 胪搿铡= 黼等 l z + 2 2 搿( a d - i 2 z d - 1 2 z ) ，( d 一1 1 2 a d 一1 2 # z ) 2 1 i z n l n o 1 i 丁一 = 1 a 。( d 一1 2 a d 一1 2 ) = 1 a 。( d 。a ) = l a m 。( ，一d “c ) = 1 ( 1 一m a x ( j ) ) 此处c = d a ，j = d c 。对模型问题( 见第四章例1 ) m a x a 。( j ) = p ( j ) = 1 o ( h 2 ) 故 p = o ( h 一2 ) 所以 8 j ( 1 a ) 8 m i n ( 3 2 0 ) 8 。关于p 的表达式很复杂，故这里选取近似值对条件数进行估计 由于p 取值一般相对比较大，故而p ( 1 一p ) 一1 代入( 3 1 9 ) 式可 得u 值( 求解的情况还可参见图3 一1 ) ： u i 1 2 或w 2 1 3 1 。( 3 2 1 ) 由式( 3 1 0 ) 可确定相应的，y 值， 7 i = 1 8 4 ：或7 2 = 1 7 9 将( 3 2 1 ) 代入( 3 1 2 ) ，可得 s m l n o 8 8 + 0 1 2 * ，( 3 2 2 ) 3 0 由( 3 2 0 ) 式，得 太原理工大学硕士学位论文 ( 3 2 3 ) 从求最小特征值开始，若将所有方程中的系数9 8 都变换为1 6 ，重 复上面的过程，则w 的值相应为( 求解的情况还可参见图32 ) ： u l 1 0 8 或u 2 1 5 将( 3 2 4 ) 代入( 3 1 2 ) ，可得相应的s 值 s m z n 0 9 4 2 8 + 0 0 6 7 2 # ，或s ，n 伽l ( 3 2 4 ) ( 3 2 5 ) 由( 3 2 0 ) 式，得 6 j ( 1 n ) ( o 9 4 2 8 + 0 0 6 7 2 p ) ，或“j o a ) ( 3 2 6 ) 另一方面 一一=粼=p=o(hinlniai 一2 ) 【aj 故吩的确可以比的条件数要小。 ；li i l ! ! ； 痧 - - 卜 - r ，。，。 ”够” - t r 7 t - 一、 、 、；力 图3 1 系数为9 8 时近似解的情况图3 - 2 系数为1 6 时近似解的情况 第四章实例 a = ( 三27 f ：j ；i ? n d ；nj ， 11 4 l l 4 1 1 4 d 。：i j一，4 ， 一，4 l一1 4 1 a = ( 窨誊主参乏 3 1 3 2 d n 2 04 42 04 42 0一4 42 0 b n = 太原理工大学硕士学位论文 4一l 14一l 1 4 一l ： 这