（计算数学专业论文）稀疏近似逆预处理方法研究.pdf

j ， 独创性声明 本人声明所呈交的学位论文是本人在导师指导下进行的研究工 作及取得的研究成果。据我所知，除了文中特别加以标注和致谢的地 方外，论文中不包含其他人已经发表或撰写过的研究成果，也不包含 为获得电子科技大学或其它教育机构的学位或证书而使用过的材料。 与我一同工作的同志对本研究所做的任何贡献均己在论文中作了明 确的说明并表示谢意。 签名：日期：年月日 论文使用授权 本学位论文作者完全了解电子科技大学有关保留、使用学位论文 的规定，有权保留并向国家有关部门或机构送交论文的复印件和磁 盘，允许论文被查阅和借阅。本人授权电子科技大学可以将学位论文 的全部或部分内容编入有关数据库进行检索，可以采用影印、缩印或 扫描等复制手段保存、汇编学位论文。 ( 保密的学位论文在解密后应遵守此规定) 签名：翩繇堑塑 局巧天纵 导师签名：尘：曼 日期：年月日 i q ： 一。ln一一-h片oo萼 k ， l l 启 量 1 矗 0 j，111|1 1 i*。一，。4。； 譬u霉 摘要 摘要 大型稀疏线性方程组的求解是许多科学和工程计算中的重要问题。当前计算 机技术发展飞速，大型科学计算已经进入大规模并行计算时代，基于并行计算环 境研究大型稀疏线性系统的高效并行算法显得尤为重要。稀疏近似逆方法具有良 好的并行性特点，且已经被证实具有较强的健壮性，它也能克服诸如不稳定等问 题，所以对稀疏近似逆方法研究具有较大的理论研究和实际意义。 本文提出了基于目标矩阵正交投影的稀疏近似逆预处理方法。这种方法考虑 的是基于f r o b e n i u s 范数最小化i i t p a l l 。的近似逆技术，其中p s ，s 是尺“”的 某一向量子空间，r 是目标矩阵，那么丁_ 1 p 是预处理子。这时原预处理矩阵 p a ( 或者4 p ) 以f 为新的目标，而不是以单位矩阵，为目标，这就是目标矩阵的 含义。由最小化f r o b e n i u s 范数l i r p a l l 。得到的“最优矩阵”只是在某种意义上的 最优，但是当p s 时是最小的。 本文首先分类介绍了大型稀疏线性方程组的迭代法和预条件迭代法。然后介 绍了稀疏近似逆方法的主要理论和当前较为成功的算法。并结合数值实验对这几 种算法的健壮性、有效性、并行性进行了比较。本文还介绍了基于目标矩阵的稀 疏近似逆方法，并进行了数值实验。最后详细阐述了基于目标矩阵的正交投影的 稀疏近似逆预处理方法。根据这个方法，在目标矩阵特殊化为正交矩阵的情况下， 得到了预处理后的矩阵的一些谱分析结果。 关键词：稀疏矩阵，稀疏近似逆预处理子，正交投影，目标矩阵 童 r l 、 z j a b s t r a c t t h el a r g es p a r s el i n e a rs y s t e m sa r e 强i m p o r t a n tp r o b l e mo fl a r g es c i e n t i f i ca n d p r o j e c tc o m p u t i n g n o w a d a y s ，t h ec o m p u t e rt e c h n o l o g y h a sd e v e l o p e dq u i c k l ya n dt h e s c i e r l t i f i c 觚dp r o j e c tc o m p u t i n gi si nt h ea g eo fp a r a l l e lc o m p u t i n g ，s ot h ee f f i c i e n t p a r a l l e la l g o r i t h mi nt h ep a r a l l e le n v i r o n m e n tb e c o m e sm o r ei m p o r t a n t t h es p a r s e a 1 ，p r o x i m a t ei n v e r s em e t h o dh a sg o o dp a r a l l e l i z a b i l i t y , a n di t i sp r o v e dt oh a v eg o o d r o b u s 协e s s t o g e t h e rw i t hi tc a ns o l v et h es t a b i l i t yp r o b l e m ，s ot h er e s e a r c hf o rs p a r s e a p p r o x i m a t ei n v e r s em a k e s a c a d e m i ca n dp r a c t i c a ls e n s ed e e p l y as p a r s ea p p r o x i m a t ei n v e r s em e t h o d ，n a m e dt h es p a r s ea p p r o x i m a t ei n v e r s e m 如o db a u s e do nt h eo r t h o g o n a lp r o j e c t i o no f t h et a r g e tm a t r i xi sp r o p o s e dm t h et h e s i s t h i sm 甜i o dc o n s i d e r st h ef r o b e n i u sn o r mm i n i m i z a t i o n0 丁一p a l l ，o fa p p r o x i m a t e i n v e 娼et e c h n i q u e s ，w h e r ep s ， si sac e r t a i nv e c t o rs u b s p a c eo f r “4a n dti s 也et a r g e tm a t r i x ，s ot 一1 pi st h en e wp r e c o n d i t i o n e r t h ep r i m a r yp r e c o n d i t i o n e d m a t r i x ，p a ( o ra p ) t a r g e t st h e m a t r i xt ，n o tt h ei d e n t i t ym a t r i xz ，t h a t i st h e m e a n i n go ft 鹕e tm a t r i x t h e r e f o r et h eo p t i m a lm a t r i xw eo b t a i nf r o m 咿一刚i i ，i s o p t i m a li ns o m es e n s e ，b u t 忙一刚忆i s m i n i m a lw h e npi ss u b j e c tt o s t h ei t e r a t i v em e t h o d sa n dt h ep r e c o n d i t i o n i n gi t e r a t i v em e t h o d sf o rs o l u t i o nt h e l a r g es p a r s el i n e a rm a t r i xs y s t e ma r ei n t r o d u c e da tf i r s t s e c o n d l y , t h ep r i m a r yt h e o r yo f s p a r s ea p p r o x i m a t ei n v e r s em e t h o d sa n dt h em e r i t sa n dd r a w b a c k so f t y p i c a ls p a r s e a p p r o x i m a t ei n v e r s em e t h o d sa r ep r o p o s e d a c c o r d i n gt ot h e n u m e r i c a lt e s t s ，t h e s e a l 硎t 衄si na s p e c t so fr o b u s t n e s s ，e f f i c i e n c ya n dp a r a l l e l i s m a r ec o m p a r e da n d a n a l y z e d t h es p a r s ea p p r o x i m a t ei n v e r s em e t h o d s t h a tb a s e do nt h et a r g e tm a t r i c e sa i l d m a ：k e sn u 加甜c a lt e s t sa r e a l s oi n t r o d u c e d f i n a l l y , t h es p a r s ea p p r o x i m a t er e v e r s e m e t h o db a s e do nt h eo r t h o g o n a lp r o j e c t i o no ft h et a r g e tm a t r i xi se l a b o r a t e di nd e t a i l ， t h e n m es p e c t r a la n a l y s i so ft h ep r e c o n d i t i o n e dm a t r i xw h e nt h et a r g e tm a t r i x i st h e o r t h o g o n a lm a t r i x i sp r o p o s e d k e y w o r d s ：l a r g es p a r s em a t r i x ，s p a r s ea p p r o x i m a t ei n v e r s ep r e e o n d i t i o n e r , o r t h o g o n a l p r o j e c t i o n ，t a r g e tm a t r i x i l 啪2 2m 2 04眦7i-y 目录 目录 第一章引言。l 1 1 课题研究背景1 1 2 预备知识4 1 2 - 1 稀疏矩阵的存贮技术。4 1 2 2 迭代法5 1 2 2 1 迭代法概论5 1 2 2 2 经典迭代法6 1 2 2 3k r y l o v 子空间迭代法。8 1 2 3 预条件迭代法8 1 3 本文工作及创新9 第二章稀疏近似逆技术1 1 2 1 稀疏近似逆技术( s a i ) 发展现状1 1 2 2 基于f r o b e n i u s 范数最小化稀疏近似逆方法1 l 2 - 3 基于矩阵分解的近似逆方法1 3 2 4 数值实验1 4 2 5 本章总结17 第三章基于目标矩阵正交投影的稀疏近似逆方法。1 8 3 1 基于目标矩阵的稀疏近似逆方法1 8 3 1 1 目标矩阵1 8 3 1 2 数值实验2 0 3 2 基于目标矩阵正交投影的稀疏近似逆方法2 7 3 2 1 概论：2 7 3 2 2 正交投影：2 8 3 2 3 目标矩阵的正交投影2 9 3 2 4 目标矩阵的特殊化3 0 3 2 5 预处理后的矩阵的谱分析3 1 3 2 5 1 预处理后的矩阵的特征值和奇异值分析3 1 i i i 3 2 5 2 剩余矩阵( m a t r i xr e s i d u a l ) 的f r o b e n i u s 范数分析3 3 3 2 5 3 预处理后的矩阵的谱条件数分析3 3 3 2 6 实例分析3 4 3 3 本章总结3 5 第四章总结和展望3 6 致谢3 7 参考文献3 8 攻硕期间取得的研究成果4 2 i v 第一章引言 1 1 课题研究背景 第一章引言 计算问题是现代社会各个领域普遍存在的共同问题，工业、农业、交通运输、 医疗卫生、文化教育、社会经济等，都有许多的数据需要加工计算，通过数据计 算和分析，以便掌握隐藏在数据背后的规律和信息。现代科学技术需要强有力的 计算能力，人类计算能力的提高包含两个方面：一是计算机性能的提高；二是计 算方法效率的提副1 1 。近二十年来，人们用计算机解决的问题越来越复杂，问题的 规模也在不断增大，有不少例子表明，计算的应用需求超过了计算机性能的提高 速度，现代许多问题是大规模的，这时，传统的数学方法几乎无能为力，这也是 计算数学的一大挑战。大型稀疏线性系统的求解是计算数学中十分重要的研究课 题，因为在许多科学与工程的问题，如设计和计算机分析，电力系统网络，计算 电磁学，排队系统【2 】等，都对大型稀疏线性系统的求解提出了要求。对这个问题的 研究具有较大的理论研究和实际意义。 求解线性方程组a x = b 的方法分为直接法和迭代法。直接法是把系数矩阵彳分 解变成可逆矩阵。它的优点是没有舍入误差，所以在某些方面得到了一定的应用。 但是由于直接法会消耗大量的运算时间和存储空间，所以这种方法只适用于对中 小型的方程组( ，z 1 0 0 0 ) 和某些大型的稀疏线性方程组的求解，当需要求解大型矩 阵时( 也就是本论文考虑的情况) ，这种方法就无能为力了。求解大型矩阵往往需要 使用迭代法求解。迭代法通过某种规则建立迭代格式，从某一初始值x ( o 出发，得 到迭代序列x ( ，石( 扪，x ( 扪，x ，逐次逼近方程组的准确解x 。迭代法以期待收敛 在一个合理的、误差在可接受的范围内为原则，被广泛应用于偏微分方程组和计 算流体学等。迭代法往往只需要较少的存储空间。但是当前还没有一个公认的最 好方法，每个方法都有它的优点和缺点，这也说明这个领域仍然有待研究和发展。 对于迭代法一般要考虑以下几个问题： ( 1 ) 怎样构造迭代序列? ( 2 ) 判断构造的迭代序列是否收敛? 以及在什么条件下收敛? ( 3 ) 若迭代序列收敛，收敛的速度是多少? ( 4 ) 得到的近似解的误差估计。 电子科技大学硕士学位论文 其中，迭代序列的收敛性和收敛速度是一个关键问题。 假设x 是方程组a x = b 的唯一解，将a x = 6 转变为等价的方程组： x = g x + f ， 得到迭代公式： a x = b ，x m = g x + 厂，( 后= 0 ，1 ，2 ，) 。 对给定的初始向量x ) ，按照上面的公式计算得到收敛于方程组近似解的向量 序列 x 。) ，则称此方法为迭代法。若对于任意的初始值x 们，当迭代次数无限增 加是，序列 x ) 都收敛于同一工，即 这时有 l i m x ( 2 ) = 工， 膏 工= g x + + 厂， 这时称迭代法是收敛的，否则称之为发散的，矩阵g 称为迭代矩阵。 当前的迭代法一般分为两类：一类是基于矩阵分裂的定常迭代法，例如 j a c c o b i ，g a u s s s e i d e l ，s o r 和a o r 等；另一类是非定常迭代法，其中应用广泛 的有c g 、g m r e s 、b i c g s t a b 等k r y l o v 子空间迭代法【3 j 。 然而迭代法的健壮性和效率上的不足影响了迭代法在处理线性方程组问题上 的应用。为了弥补迭代法的这些不足之处，可以在求解方程组之前对方程组进行 预处理，即采用预处理技术。预处理技术是构造预处理子尸使得： p a x = p b ( a p y = b ，x = e y ) ， 这时称p 为左( 右) 预处理子。选择合理的预处理子可以加快迭代法的收敛速度，甚 至可以使原来不收敛的问题变得收敛。对于k r y l o v 子空间技术，迭代速度与系数 矩阵彳的特征值的分布紧密相关，即特征值的分布越集中在1 附近，收敛速度就 越快，反之，收敛速度越慢甚至是发散的。对预处理子的构造要考虑以下几点因 素 1 】： ( 1 ) 线性方程组的背景和来源。 ( 2 ) 系数矩阵的结构特点。如是否是稀疏的，对称的，是否是特殊矩阵。 ( 3 ) 预处理部分计算量较小。 2 第一章引言 ( 4 ) 预处理后的矩阵，即p 。a 与单位矩阵，的近似程度。 般来说，较好的预处理予( f i gc p u 运行时间少，收敛地较快) 应该满足以下 几点： ( 1 ) 在某种意义上，预处理子尸的逆是系数矩阵彳的一个较好的近似。事实上， 预处理的主要目的就是使得剐为单位矩阵的一个近似。 ( 2 ) 构造预处理子p 的代价是可以接受的( 即代价不大) 。 ( 3 ) 新系统比原系统更易于求解。 与迭代法相对应，预条件迭代法也可以分为以下两大类： ( 1 ) 预条件定常迭代法：主要是基于矩阵的分裂，构造预处理子，再应用于经 典迭代法，如g a u s s s e i d e l ，s o r ，a o r t 4 1 。在此类预处理子中，修正预条件方法 较多。此类方法的特点：理论支持强、计算代价小，但是效果相对而言要差些。 修正预条件方法在某种程度上可以看成是g a u s s 消元法的扩展。 ( 2 ) 预条件非定常迭代法：主要为预处理与k r y l o v 子空间方法( 如c g ， m i n r e s ，g m r e s ，b i c g s t a b ) 相结合。基于纯代数技术构造预处理子的方法已经 被广泛的应用到科学与工程计算领域，包括不完全分解( i l u ) 口】，稀疏近似逆( s a i ) 和代数多水平等。 随着并行和分布式处理器的发展，构造可并行化的预处理子成为判断大型稀 疏线性系统技术是否高效的关键。为了构造可并行化的预处理子，在最近几年有 很多学者把目光投向稀疏近似逆方法的研究。基于稀疏近似逆方法的预处理技术 已经被证实具有较强的健壮性，它能克服诸如不稳定等问题。所有的近似逆技术 的假设是，对某一稀疏矩阵么，有可能找到一个在某种意义上好的稀疏矩阵尸， 使得p 近似的等于彳。但是，要构造这样的预处理子尸是比较困难的，因为稀疏 矩阵的逆通常是密集的，更确切地说，已经证明不可约稀疏矩阵的逆是密集的。 所以，如何保持p 的稀疏性仍是需要研究的问题。 近二十年，稀疏近似逆方法由于便于并行计算，得到了快速发展。目前，有 两种方法求得稀疏近似逆预处理子，它们的区别之处在于稀疏近似逆预处理子能 否表示成一个矩阵。也就是说方法一在于稀疏近似逆预处理子是由一个矩阵求得， 如f r o b e n i u s 范数最小化稀疏近似逆方法。方法二得到的稀疏近似逆预处理子表示 成几个矩阵的乘积，这种方法叫做因子化稀疏近似逆，例如a i n v 、s a i n v 方法( 参 见 9 1 4 ) 。s a p i 方法在 1 5 ，1 6 中有介绍，b e n z i 和t u m e i n 通过数值实验比较了 a i n v 、s a i n v 两种方法【7 j 。有关向量和平行结构稀疏近似逆技术见文献 5 8 。 基于f r o b e n i u s 范数最小化稀疏近似逆技术是近似逆方法中应用最广泛、并行 电子科技大学硕士学位论文 最好的一类方法。f r o b c n i u s 范数最小化的主要思想是把原方程组血= b 转化为 小值问题： m i n i - a p i i ，= 喜曾卜和川， 中s 是一非零元模式，印e 分别代表，、尸的第列。这样，得到的p 就是彳的 个近似逆。而问题就变成了求解n 个独立的线性方程问题，所以它能很好的用于 行计算。 1 2 预备知识 1 2 1 稀疏矩阵的存贮技术 如果一个矩阵中只有少数元素不为零，则称这样的矩阵为稀疏矩阵，反之是 稠密矩阵。大型稀疏矩阵由于其元素较多，如何节省存贮空间是必须要考虑的问 题。在实际求解问题中稀疏矩阵的阶数往往达到几十万阶甚至几亿阶，如果仍然 采用稠密矩阵的存贮技术，将对存贮空间和计算能力提出巨大的挑战 2 3 1 。如果用 稠密矩阵方法求解稀疏线性方程组，存贮量和计算量都将十分大。例如，对一个咒 阶稀疏线性方程组，如果采用稠密g a u s s 消元法求解，其时间复杂度为o ( n 3 ) ，存 贮量为o ( n 2 ) 。 稀疏矩阵中非零元素往往较少，特别是从有限元方法和偏微分方程中得到的 稀疏矩阵，常常每行每列只有几个非零元，如何利用这种稀疏性，尽量减少存储 量和计算量是值得研究和考虑的问题。下面介绍几种主要的稀疏存贮技术。 ( 1 ) 对角线存贮法【2 剐 设矩阵a r ”“中非零元为 a “+ l ：l j 珑，1 f ，z ) ，则可以用b = ( ) r “ 与整型向量y = ( ，l ，乞，m ) r 存贮彳，其中 分。 当1 i + t ，l 其他 采用这种方法，对于对称矩阵，就可以只存贮稀疏矩阵的上三角或下三角部 当l if 2 ，乙为连续整数时，这种方法称为等宽带存贮法。此时，不需要存贮 4 + ， n v 口 ，l = 第一章引言 整型向量y ，只需要记下y 的最大分量和最小分量。 已经被证明的是，带状矩阵的l u 分解会保持原矩阵的带状结构，不会在带状 结构之外增加填入元，所以，等宽带存贮法的优点十分明显。 ( 2 ) 对称矩阵的变宽带存贮法1 2 列 设彳是一个n 阶的对称矩阵，现在只存贮它的下三角部分。设m ；为第f 行的局 部宽带，则可以对4 的下三角部分的元素进行逐行存贮，且各行的元素从最左边 的非零元开始到对角元，这样可以将这些元素存贮到一个向量x 中。为了访问么元 素，需要在引入一个, 1 维整型数组y ，其中，第f 个元素指明么的第f 个对角元存贮 在x 中的位置。 对对称矩阵进行l d l r 分解时，将不会在其包络之外的位置进入填入元，所以 采用变宽带存贮法相当合适。 ( 3 ) 坐标存贮法【2 3 】 对于非零元结构不规则的稀疏矩阵来说，有多种可供选择的存贮法，其中坐 标法就是很好的一种。坐标法存贮一个稀疏矩阵彳时，需要利用一个数据类型与彳 相同的向量x ，以及两个整型向量x 1 、工2 。其中，x 用于存放彳中的非零元，x 1 、x 2 分别记录下这些元素的行和列。 这种存贮法的缺点是不能实现对矩阵元素的有效访问，在最坏情况下，访问 矩阵的某元素需要遍历整个向量x 。 1 2 2 迭代法 1 1 2 2 1 迭代法概论 众所周知，对于中小型的方程组( 力1 0 0 0 ) 和某些大型的稀疏线性方程组时使 用直接法求解，但是当系统是大型矩阵时( 也就是本论文考虑的情况) ，直接法就无 能为力了。这时只能选择迭代法求解。迭代法是通过某种规则建立迭代格式，从 某初始值x o 出发，得到迭代序列石( ，x ( 孙，x ( ”，x ( ，逐次逼近方程组的准确解 ，。迭代法以期待收敛在个合理的范围、误差在可接受的范围内为原则，被广 泛应用于偏微分方程组求解和计算流体学等【2 】。迭代法只需要较少的存储空间，但 是当前还没有一个公认最好的迭代法，每种方法都有它的优点和缺点，这也说明 这个领域仍然有待研究和发展。每一种迭代方法都有一定的应用价值和应用范围， 对于同一具体的问题，有的迭代法收敛，有的不收敛。并且，迭代速度也有差别， 所以对迭代法的研究具有广阔的空间【2 3 1 。 电子科技大学硕士学位论文 常见的迭代法有j a c o b i 迭代法、g a u s s s e i d e l 迭代法。为了改善迭代法的收敛 和收敛速度，h a d j i d i o m s 2 4 】于19 7 8 年提出了一种快速超松弛迭代法( a o r 迭代 ) ，它含有两个参数纸y ，当0 9 = y 时就是逐次超松弛迭代法( 也就是s o r 迭代 ) 。另外。也发展了一些新的迭代法，如s s o r 迭代法、s a o r 迭代法、p s d 迭 法等，目的都在于改善迭代矩阵的收敛性和收敛速度。当前有许多求解线性系 的迭代法的文献，见 4 3 ，4 5 ，5 7 。 1 2 2 2 经典迭代法 本文讨论的是线性方程组 a x = 6 ，a r n x n9 工，b r ”“， 其中，彳= ( ) l 驯；。是系数矩阵，x r 蒯，b r n x l ，x 是未知变量。那么残差变量 是b a x 。可以把系数矩阵彳写成如下形式： a = d e f ， 其中，d 是a 的对角线元素，即d = d i a g ( a i ，口肼) ，一e 是严格对角线以下元素， 一，是严格对角线以上元素。即 或者 d = a 1 1 0 0 4 2 2 00 f = 0 一口1 2 o0 oo 一口i 一a 2 n 0 o o 一a n 2 把系数矩阵彳写成这种形式后，就可以方便地对经典迭代进行介绍。 ( 1 ) j a c o b i 迭代法 在j a c o b i 迭代法中，其迭代公式： x 。+ 1 = d 一1 ( e + ，) x + d b ， 6 第一章引言 z “1 = ( 一弓”+ 6 f ) ， ( 其中f = 1 ，咒) u - i j 亭i 称毋= i - d 。1 b 为雅可比迭代矩阵。 ( 2 ) g a u s s s e i d e l 迭代法 g a u s s s e i d e l 迭代法的迭代公式： x ( k + 1 ) = ( d e ) 一1 f x + ( d - e ) b ， 1i - l n + 1 = ( 一拶“一勖矽+ 包) ， ( 其中f _ 1 ，甩) t - 矗j = lj = i + l 称b o s = ( d e ) - 1 ，为雅可比迭代矩阵。 g a u s s s e i d e l 迭代法可以被看做是对j a e o b i 迭代法的一个改进。在每一步迭代 中，根据算出的新近似解的一个分量z “”来构造下一个分量罐n ，用新分量墨“ 来代替老分量z “”进行计算，这样带来的好处是只需要n 个单位元来存贮近似分 量，而且，这样的新的近似值通常更加接近精确值，这就是g a u s s s e i d e l 迭代法较 j a e o b i 迭代法而言的优越性。 ( 3 ) s o r 迭代法 s o r 迭代法的迭代公式是： 其中 x k + l = l 吆+ c o ( d - t o e ) b ， l w = ( d - w l ) 。1 ( ( 1 一w ) d + w u ) ， 称匕为s o r 迭代矩阵，称称为超松弛因子( 从严格意义上来讲，当c o 1 时才称 为超松弛，当c o 1 时称为低松弛，在这里不3 dx e 另j ) 。下面将对经典迭代法的收敛 速度比较【4 1 。 表1 1 经典迭代法的收敛速度比较 jg ss o r 迭代法收敛速度 1 2 , r 2 h 2 刀2 h 2 2 r c h ，兰2 ( 1 - x h ) 在表1 - 1 中，其中h = l l ( n + 1 ) ，h 是差分网格的边长。从表1 - 1 中可以得到， g m r e s 方法是一个k r y l o v 子空间方法。其中k = 以，l = a k 。，如是第，z 个 k r y l o v 子空间。 算法1 2 231 2 4 3 1 ( g m r e s 算法) 1 计算= b - a x o ，p ：- i i , j l l ：，u ：= p 2 定义疗。= ) 。搿s 。+ l ，班。，疗，是( 肌+ 1 ) m 矩阵。令曰m = o 3 f o rj = 1 ，2 ，md o 4 计算_ ：= a v j 5 f o ri = 1 ，2 ，d o 6 7 i ，：= ( w j ，v ) 7 w j ：= 一u 8 e n d d o 9 乃“，- - w i l l ：i fh j + l ，j = o ，令，z ：= ，g 。t o1 2 1 0 v 川= w j 乃+ f 1 1 e n d d o 1 2 计算l i q 一f = y l i ：的最小值y 。和= + 虼 1 2 3 预条件迭代法 求解大型稀疏线性方程组是工程计算中的重要课题。迭代法健壮性和效率上 的不足影响了迭代法在处理线性方程组的问题上的应用，通过预处理技术就可以 改善这些缺点。选择合理的预处理子可以加快迭代法的收敛速度，甚至可以使原 来不收敛的问题变得收敛。目前，多数采用的是预条件k r y l o v 子空间方法。现在 已经有许多k r y l o v 子空间方法用于求解一般的线性系统，例如g m r e s ， b i c g s t a b 及q m r ( 参见 2 ) 。为了增加有效性，通常对k r y l o v 子空间方法进行 第一章引言 预处理。普遍认为选择一个好的预处理子比选择一个k r y l o v 子空间方法更为重要。 对预处理子的有效选择始终是需要讨论的课题。对预条件迭代法的详细介绍参见 2 5 ，其他的预条件迭代法见 3 7 ，3 8 ，4 7 ，5 2 ，5 3 。 预处理技术通过预处理后得到与原问题等价的方程组，但是使问题变得更简 单。通常使用左或右预处理子p 使原线性方程组变成易于求解的形式 4 毋= 6 ，x = b y 或p a x = p b ， ( 1 1 ) 然后用迭代法求解公式( 1 1 ) 式。构造预处理予的方法有很多，如预条件k r y l o v 子 空间方法、不完全分解方法( i l u ) 【5 4 】、s s o r 方法、多项式方法、稀疏近似逆方法 ( s a i ) 等。 对于预条件k r y l o v 子空间方法，迭代速度与系数矩阵彳的特征值分布情况密 切相关。特征值分布越集中在1 附近，收敛速度就越快，否则就越慢，或者根本 不收敛。因此，预处理子的选择应该考虑预处理后的系统剐有更紧凑的特征值分 布【3 1 。 另一种构造预处理子的方法是基于稀疏近似逆的预处理子技术，这项技术在 近年来得到了重视和发展。已经得到证实的是这类预处理子可解决些i l u 预处 理子【2 】难以解决的问题，并且可以代替传统i l u 预处理子。基于稀疏近似逆的预处 理技术的共同思想是：求得一个稀疏矩阵p a ，使原问题a x = b 的解就变成 工p b 。基于f r o b e n i u s 范数最小化的稀疏近似逆方法的主要思想是转化为求解最 小值问题： m t n l l z - a p ，= 喜曾卜和川f ， 这样预处理过程可以很容易地并行执行。 1 3 本文工作及创新 ( 1 ) 分类讨论了当前的稀疏近似逆技术，并结合数值实验，比较和分析现有的 几种较成功的稀疏近似逆方法的优缺点( 强壮性、稳定性、并行性等) 及其适用的范 围。 ( 2 ) 介绍基于目标矩阵的稀疏近似逆预处理技术，通过数值实验比较了选择不 同的矩阵范数和目标矩阵的结果。 9 电子科技大学硕士学位论文 ( 3 ) 提出了基于目标矩阵的正交投影的稀疏近似逆预处理技术，当目标矩阵的 殊化后对预处理后的矩阵的进行了谱分析。得到的结论是当目标矩阵t 是正交 阵时，预处理后的矩阵的特征值、奇异值、“剩余矩阵”( 即丁一p a ) 的范数和谱 件数与系数矩阵的特征值有关。更确切的，当i 丸i = 1 ( 或者吒= 1 ) 时，其他的特 值、奇异值和谱条件数的模也等于1 ，且l l r - q l ，等于0 。 1 0 第二章稀疏近似逆技术 第二章稀疏近似逆技术 2 1 稀疏近似逆技术( s a i ) 发展现状 大型稀疏线性系统的求解是许多科学和工程计算中重要且基本的问题之一， 随着科学技术的不断创新和发展，具体问题的复杂性不断提高，在高维微分方程 数值解、有限元分析、核物理与流体力学计算、电力系统的优化设计、地震石油 数据处理及数值天气预报等许多领域的大规模科学工程计算和数值处理中都会遇 到对称或非对称大型稀疏线性系统的求解问题。另一方面，当前计算机技术发展 飞速，大型科学计算技术已经进入大规模并行计算时代，基于并行计算环境研究 大型稀疏线性系统的高效并行算法显得尤为重要。当前国内外有很多学者专注与 此项技术的研究4 1 ，并且得到了相当好的成果。现在存在着一些比较优秀的算法， 在一些问题上可以与传统的不完全i l u 方法相媲美，在某些问题中，甚至优于i l u 方法。在接下来的文章中也对它们适用情况和效果做了对比。 当前，稀疏近似逆技术主要分为两类：基于f r o b e n i u s 范数最小化和基于矩阵 分解的近似逆方法。这两类中都有一些代表性的算法，例如，基于f r o b e n i u s 范数 最小化方法中，有m r 算法、s p a i 算法、s e l f - p r e c o n d i t i o nm r 等；第二类中有a i n v 算法、f s a i 算法及逆i l u 算法【2 j 。 2 2 基于f r o b e n i u s 范数最小化稀疏近似逆方法 基于f r o b e n i u s 范数最小化稀疏近似逆方法是近似逆方法中应用最广泛、并行 性最好的一类方法。f r o b e n i u s 范数最小化的主要思想是转化为求解的最小值问题： 曾卜胛| | ，= 喜刊旷却， 其中s 是一非零元模式，巳、乃分别代表，、p 的第，列。这样，得到的矩阵p 就 是彳的一个近似逆。问题就变成了求解1 1 个独立的线性方程问题，它能很好的用于 并行计算。 确定s 的稀疏模式就可以直接计算尸了，这是关键的问题。对于s 的静态选 电子科技大学硕士学位论文 可以采用彳或a - 1 的稀疏模式，但是这样得到的结果往往是不理想的。已经被 的是，即使么是一个稀疏矩阵，它的逆4 - 1 往往是稠密的，而且，不可约矩阵 肯定是稠密的【7 8 】。但是有一点值得庆幸，彳- 1 往往含有许多绝对值很小的元素， 省略可以化简为稀疏的。t h o m a sh u c k l e t 9 提出了预测彳一的稀疏模式的方法。 ，对于一般的矩阵预测它的逆的稀疏模式较为困难，有许多学者把目光投向 适应的动态稀疏模式s 的研究。这种方法通常猜测出一个初始稀疏模式s ，然 过迭代达到要求。其中，g o u l d 和h u c k l e 27 】提出的s p a i 算法较为成功，现在 到了很多的应用。算法具体如下： 算法2 2 1 【2 7 1 ( s p a i 算法) 然后删除指标k ，除了使虞最小的m a x 个的k 。 6 决定新的指标集，并对a ( 1w i ，jw j ) 进行q r 分解，然后求解新的最小平方 问题，计算新的剩余向量，= e j a m ，1 w i ，j u - 厂。 这个算法需要设定一个初始的稀疏模式和其他几个参数，包括了最大残量的 阈值占及每次迭代的新的非零元的最大个数m a x 。s p a i 的具体实现参见 1 0 。 下面介绍其他的基于f r o b e n i u s 范数最小化稀疏近似逆方法m i n i m a lr e s i d u a l 算法。 算法2 2 2 【2 1 ( m i n i m a lr e s i d u a l 算法) 令初始值m = m o f o rj = 1 ，n ，d o f o rf = 1 ，z ，d o 1 2 第二章稀疏近似逆技术 m ，= m e ， j 乃= e j a m j 口，= ( 0 ，么0 ) ( a r j ，彳，：，) m j 。m i q i r j 对m ，进行数值筛选，当元素小于某给定值孝时，改为0 e n d d o e n d d 0 2 。3 基于矩阵分解的近似逆方法 基于矩阵分解的稀疏近似逆方法有许多优点： ( 1 ) 它构造的预处理子可以保持系数矩阵彳的对称正定性，显然，这是它的一 大优点。而且，如果么是对称正定矩阵，那么预处理后的线性系统可以用c g 法求 解。 ( 2 ) 与其他的近似逆方法相比，它具有计算量小、定义的参数少的优点，这样 减少了系统的不确定性。 ( 3 ) 它对重排( r e c o r d 舐n g ) 很敏感，这就意味着，可以通过适当的重排技术使 得非零元填充现象减少，从而提高收敛效率【3 0 1 。 但是，基于矩阵分解的稀疏近似逆方法也有它的缺点，它的并行性不如基于 矩阵分解的方法，a i n v 法已经被证实难以并行实现。 下面将介绍a i n v 算法。在a i n v 算法中，构造预处理子需要计算两个向量 组 乞) ：l 、 w 。t 1 1 ，它们是a 一双共轭的，即嵋a z ，= 0 ，当且仅当f j 。假设a r “ 是非奇异矩阵，得到彳。1 和 2 i ) ：。、 w 。n ：，的关系式如下： 那么 z = z 1z 2 ，z 。 ，w = w l ，w 2 ，】， 矽7 a z = d = d i a g ( p i ，p 2 ，p 。) ， 其中p i = w t a z i 0 。显然得到w 、z 非奇异且有： 1 ，w = l - 7 ，这里 。如果用a ，r 表示 算法2 3 1 【2 1 ( a i n v 算法) 令z = e i ，p ：= 口l l f o ri = 2 ，nd o z o ) = e t f o rj = 1 ，f - 1d o 曩产 = 口j 考卜d i 门= z 卜一( 卜 p y - 1 ) ) 矽。1 e n d d o 卜1 ) - 口；乏卜” e n d d o 为了保证预处理子的稀疏性，矿和z 在迭代过程中消去了一部分的非零元， 这些元素的绝对值小于一给定正值( 丢弃误差) 。这样导致不完全分解： 一一一一一一1 一r z z w w ，且d d 。因子分解结果为m = z dw 2 4 数值实验 近年来，提出了许多新的稀疏近似逆方法 5 】，基于稀疏近似逆技术的预处理技 术已经被证实具有较强的健壮性，它能克服诸如不稳定等问题。m i c h e l eb e n z i 和 m i r o s i a vt u m e 3 】通过数值实验对其中一些方法进行了比较。下面通过数值实验比较 s p a i 方法与不完全因式分解i l u ( 0 ) t ”】两种典型的方法，对它们在各个方面的优劣 都进行了比较。该实验在c r a yc 9 8 单处理器环境下进行，采用了2 0 个大型矩阵， 它们主要来自于h a r w e l l b o e i n g 收集矩阵【1 1 】。这些矩阵最初出现在石油储存库仿 真、电路设计、半导体设备模型的问题中。 下面先介绍i l u ( 0 ) 算法。 1 4 第二章稀疏近似逆技术 算法2 4 1 t 2 1 ( i l u ( 0 ) 算法) f o ri = 2 ，zd o f o rk = 1 ，i 一1 ，a n df o r ( f ，七) n z ( a ) ，d o 计算= a i k a 珏； f o rj = k + 1 ，m ，a n df o r ( f ，) n z ( a ) ，d o 令口 ，2 口 ，一口擅口材。 e n d d o e n d d o e n d d o 表2 1 实验矩阵及m u ( 0 ) 算法在b i c g s t a b 预处理条件下的运行结果 实验矩阵 nn zi t sp t i m ei t t i m e 3 d c d8 0 0 0 5 3 6 0 01 50 0 5 8o 5 5 1 a l e 3 d1 5 9 0 4 5 0 9 0 3 80 2 1 80 2 9 3 o r s r e g l2 2 0 51 4 1 3 32 9o 0 1 6o 3 0 1 o r s i r r l1 0 3 06 8 5 82 70 0 0 8o 1 3 4 o r s i r r 28 8 65 9 7 02 80 0 0 70 1 2 1 s a y l r 3