（运筹学与控制论专业论文）子空间迭代法的两种rayleigh商加速技术.pdf

南京航空航天大学硕士学位论文 摘要 本文研究求解大型对称矩阵特征值问题的子空间迭代法。为了加速子空间迭代 法的收敛性，我们应用r a y l e i g h 商最小化技术得到两种新的改进算法。 第一种改进算法是用r a y l e i g h 商加速子空间迭代法。它用每次迭代得到的r i t z 矩阵和将r i t z 反迭代得到的矩阵，二者构造一个带参数矩阵的线性组合，适当选取 参数矩阵，使组合后的矩阵的列向量的r a y l e i g h 商达到最小，从而更接近最小特征 向量。 第二个改进算法是用带位移的r a y l e i g h 商加速子空间迭代法。与第一个改进算 法类似，都是构造了一个带参数矩阵的线性组合，不过它选用的矩阵不同，是用 r i t z 矩阵和将r i t z 矩阵带位移反迭代后得到的矩阵构造的，同样通过选取适当的参 数矩阵，使其r a y l e i g h 商达到最小，从而加速子空f d 】的收敛性。 本文分析了这两个改进算法中参数矩阵的选取及其性质，数值稳定性和算法的 收敛性，并给出了数值实验，将新方法和原始子空间方法进行比较，数值实验表明 新改进的两个算法比子空间方法更优越。 关键词：对称正定矩阵，特征值，特征向量，子空间迭代法， r a y l e i g h 商 反迭代，带位移的反迭代。 子空间迭代的两种r a y l e i g h 商加速技术 a b s t r a c t i nt h i sp a p e r ，w ec o n s i d e rt h es u b s p a c ei t e r a t i o nm e t h o df o rs o l v i n gl a r g es y m m e t r i c e i g e n p r o b l e m s ，i no r d e rt o a c c e l e r a t et h ec o n v e r g e n c er a t e ，w ei m p r o v et h eo r i g i n a l m e t h o dw i t ha c c e l e r a t i o nt e c h n i q u e ，a n dp r e s e n tt w on e wa l g o r i t h m s i nm yt h ef i r s tp r o p o s e da l g o r i t h m ，ac o m b i n a t i o no f t h el a t e s tm a t r i xr e c e i v e db y i n v e r s ei t e r a t i o na n dt h er i t zm a t r i xi sf o r m e di n v o l v i n ga nu n d e t e r m i n e dp a r a m e t e r m a t r i x ，w h i c hi s d e t e r m i n e db ym i n i m i z i n gt h er a y l e i g hq u o t i e n t ，t h e ni tw i l ln e a rt h e m i n i m a l e i g e n v e c t o r i nm yt h es e c o n dp r o p o s e da l g o r i t h m ，w ec r e a t eac o m b i n a t i o na st h es a m ea st h e f i r s to n e ，b u ti nt h es e c o n do n et h ec o m b i n a t i o ni n v o l v i n ga nu n d e t e r m i n e dp a r a m e t e r m a t r i x ，w h i c h i sd e t e r m i n e db ym i n i m i z i n gt h er a y l e i g hq u o t i e n ti sf o r m e db yt h e l a t e s tm a t r i xr e c e i v e d b y as h i f t e di n v e r s ei t e r a t i o n a n dt h er i t z m a t r i x ，t h e n a c c e l e r a t et h ec o n v e r g e n c er a t eo f s u b s p a c e i nt h ep a p e r w ea n a l y s i st h ec h o o s i n gm e t h o do f t h ep a r a m e t e rm a t r i xa n di t ss o m e p r o p e r t y , t h en u m e r i c a ls t a b i l i t ya n dc o n v e r g e n c e o u rn u m e r i c a lr e s u l t s s h o wt h a tt h e t w o p r o p o s e da l g o r i t h m sa r es u p e r i o rt ot h eo r i g i n a ls u b s p a c ei t e r a t i o nm e t h o d k e y w o r d s ：s y m m e t r i cm a t i x ，e i g e n v a l u e ，e i g e n v e c t o r ，s u b s p a c e i m r a t i o nm e t h o d r a y l e i g hq u o t i e n t ，i n v e r s ei t e r a t i o n ，t h e s h i f t e di n v e r s ei t e r a t i o n 。 承诺书 本人郑重声明：所呈交的学位论文，是本人在导师指导下，独立进 行研究工作所取得的成果。尽我所知，除文中已经注明引用的内容外， 本学位论文的研究成果不包含任何他人享有著作权的内容。对本论文所 涉及的研究工作做出贡献的其他个人和集体，均已在文中以明确方式标 明。 本人授权南京航空航天大学可以有权保留送交论文的复印件，允许 论文被查阅和借阅，可以将学位论文的全部或部分内容编入有关数据库 进行检索，可以采用影印、缩印或其他复制手段保存论文。 作者签名逊l 挫 日期：! ! ! ! ：墨! 型 南京航空航天大学硕十学位论文 注释表 表示实数域上维向量空间 表示所有n x m 实矩阵的全体 表示向量x 的转置 表示矩阵x 的转置 表示n n 实对称矩阵 表示a 的第f 个特征值 表示爿的属于兄规范化正交特征向量 表示向量2 范数或矩阵2 范数 表示矩阵a 的f r o b e n i u s 范数 表示阶单位矩阵 表示m o d i f i e dg r a ms c h m i d ti t - 交化过程 表示向量x 关于矩阵a 的r a y l e i g h 商 表示向萤x 和y 正交 x 上u 表示向量x 和u 的各列及各列所张成的于空问正 交。 _ = ( x l ，屯，”，矗) 表示一个n x ；，l 矩阵，其中x ，表示x 的第列。 s p a n ( z 、 表示z 的列所张成的了空间。 ， ： i s 订 o 胪 ， 4 “ 懈 删 山 南京航空航天大学硕士学位论文 第一章绪论 矩阵特征值和特征向量的计算是数值线性代数的基本问题之一。这个问题 的研究最早始于1 9 世纪中期，因为求矩阵特征值和特征向量相当复杂，因此直 到2 0 世纪5 0 年代电子计算机的出现，才促使这一领域飞速发展。矩阵特征值 问题的研究相当活跃，其研究具有重要的理论意义和巨大的应用价值。许多复 杂的实际问题，诸如飞机、空间站、海洋石油钻井平台等大型复杂结构的动力 分析和稳定性分析都可转化为大型稀疏矩阵的特征值问题。这些应用一直在推 动着这个领域在整个2 0 世纪的发展。 s i a mjm a t r i x a n d a p p l ) ) 杂志每年有大 约4 0 的文章与特征值问题的研究有关，而且可能会有越来越多的文章关注这 一领域。 求矩阵的特征值已有一些经过理论分析和实践检验的有效数值方法，j a c o b i 最早于1 8 4 6 年提出了j a c o b i 方法“0 1 ，他通过旋转变换将矩阵变为严格对角占 优阵。2 0 世纪5 0 年代，这个方法受到了更多的重视，并由此产生了许多更有 效的方法，例如j a c o b i d a v i d s o n 方法等。至今人们仍密切关注着这类方法。2 0 世纪6 0 年代，出现了更有效的方法，如g i v e n s 方法，h o u s h o l d e r 方法，它们 在求特征值的过程中较j a c o b i 方法具有更大的优越性。 另一个重要的方法是始于1 9 1 3 年的幂迭代1 7 l ，它同样有着重要的意义，而 且是许多其它算法的基础。幂法用于求绝对值最大的特征值，但在求多个极大 特征值的时候不是很有效。在幂法的基础上，人们发展了q r 迭代，l a n c z o s 方法，a m o l d i 方法等。现在相当多的迭代方法仍是通过一些加速技术来提高幂法 的收敛速度。 另外一种方法是在反迭代中不断更新位移，1 9 世纪7 0 年代，r a y l e i g h 首先 用当前迭代向量的r a y l e i g h 商作为位移，给出了r a y l e i g h 商迭代( r q i ) 。r q i 方法把当前迭代向量作为近似特征向量。对于对称矩阵，已经证明g q i 方法是 立方收敛的。 由于正交化技术具有非常优越的稳定性，人们提出了q r 方法用来求稠密 对称矩阵的所有的特征值和相应的特征向量。幂法与q r 正交化方法结合产生 了同时迭代法。1 9 5 7 年，b a u e r 首先提出求解矩阵特征值问题的同时迭代法。 六十年代后期和七十年代初期，j e n n i n g s l 2 1 1 3 】l “1 和r u t i s h a u s e r 【4 l 【5 1 等人进一步发 展了同时迭代法。七十年代中期以来，人们应用r a y l e i g h r i t z 过程不断完善同 时迭代法，得到了求解矩阵特征值问题的子空间迭代法。 子空问迭代的两种r a y l e i g h 商加速技术 随着科学研究和工程技术的不断发展，求解大型矩阵特征值问题在许多科 学研究、工程技术等实际问题中有着举足轻重的作用。最近2 0 多年，大型矩阵 特征值问题的数值计算已取得了许多重要进展。9 0 年代后期以来。大型矩阵特 征值问题的加速技术已成为数值代数领域的热门研究课题。虽然现在有许多有 效的工具和方法，但仍然有许多问题没有解决。 对求解大型对称矩阵特征值问题目前，最常用的方法有子空间方法， l a n e z o s 方法，d a v i d s o n 方法，共轭梯度法等。人们对这些方法的收敛性和稳 定性作了深入细致的研究，著给出了一些变形，例如，广义d a v i d s o n 方法、 j a c o b i d a v i d s o n 方法、d a v i d s o n l a n z o s 方法等，这些新方法在一定程度上加速 了算法的收敛性，优于原始的方法。 目前，子空间迭代法已成为求解大型、稀疏矩阵特征值问题最有效的方法 之一，它的优点是可靠性高，但当要求的特征值与不要求的特征值之间的间隔 很小时，子空间迭代法的收敛性就会变慢，从而运算量加大，因此有必要对子 空间迭代法作进一步的研究，加速它的收敛性，提高算法的有效性。 计算矩阵特征值的一些方法，如坐标关系法1 5 1 16 j 【”，梯度法2 i ，共轭梯度 法1 2 | 1 3 1 【4 1 【7 1 等都采用一个线性组合来优化r a y l e i g h 商。例如在共轭梯度法中，由 第k 步迭代向量以按照如下方法得到第k + l 步迭代向量x 。： x 。= h + t k q i ，这里t 女是参数，q 。是搜索方向。 适当选取参数f 。使r ( x 。) 达到最小从而加速迭代向量收敛到最小特征向量。 我们从共轭梯度法中得到启示，可以将这个线性组合推广成矩阵的形式： x = x + g a 。 其中是对角参数矩阵，x = ( x ，x “，x 7 “1 ) 。 通过适当选取一个参数矩阵使尺( x ) ( j = 1 2 ，p ) 达到最小，从而加速 收敛性。另一方面，对于对称矩阵，带位移的反迭代具有立方收敛的性质，因 此可利用它来加速子空间迭代法的收敛性。 本文提出了子空间迭代法的r a y l e i g h 商加速技术和带位移的r a y l e i g h 商加 速技术两种方法，并对新方法进行了理论分析和数值实验。 本文的主要工作如下： ( 1 ) 对求解大型对称矩阵特征值问题的子空间迭代法提出了r a y l e i g h 商加速技 术，讨论了参数矩阵的选取及其性质，算法的数值稳定性和收敛性。 ( 2 ) 提出了带位移的r a y l e i g h 商加速技术，同样讨论了参数矩阵的选取及性质， 算法的数值稳定性和收敛性。 ( 3 ) 对上述两种新方法进行了数值实验，将结果与原始子空间迭代法进行了比 南京航空航天大学硕士学位论文 较，数值实验表明对于要求的特征值与不要求的特征值分离很小的情况 两种新方法都优于子空间迭代法。 子空间迭代的两种r a y l e i g h 商加速技术 第二章r a y i e i g h - r i t 2 逼近与子空间迭代法 2 1r a y l e i g h r i t z 逼近的相关理论 在本文中所有算法的执行过程中，都应用r a y l e i g h r i m 过程求大型对称矩 阵近似特征值和近似特征向量，所以我们在本章中先给出r a y l e i g h ，r i t z 逼近的 有关理论。 a 是r “”中的一个实对称矩阵，它将尺”映射成r ”。 定义2 1 1 设s 是胄“中的一个m 维子空间沏s n ) ，a 对s 的限制是算子 甲：s s ，其定义为：对任何x s ，w x = x , a x ，其中万。是对s 的正交投影。 设q 是n x 撒列正交规范矩阵，其列形成s 的一组规范正交基，则 万一= q q l 。 因此，对任何x e s ，有疗。x = x 。故 q x = 疋一工= 7 r a z rx = q q 7 a q q 7 x = q h v ， 其中h = q 7 4 q ，v = q x 。 若q 是n x m 列正交规范矩阵，s = s p a n ( q ) ，则r a y l e i g h 商矩阵 p ( q ) = q 7 爿q = h 是在由q 所确定的s 的基下限制甲的矩阵表示。 定义2 1 2 设a 是阶实对称矩阵，q 是n x m 列正交规范矩阵，令 h = q 1 a q 。 设 为目的特征值，y 为相应的特征向量，称“为a 在子空间s = s p a n ( q ) 中的 r i m 值- q = 勿为相应的r i t z 向量。和g 对爿的特征值和特征向量的逼近， 称为r a y l e i g h - r i t z 逼近。 定理2 1 11 2 0 1 子空间s = s p a n ( q ) 中的r i t z 值”。和r i t z 向量q ，( i = 1 ，掰) 是 甲的特征值和特征向量。 如果s 是彳的一个埘维不变子空问，口是由s 的正交规范基为列所构成的 矩阵，则必有n l m 矩阵日，使 r = “q o h = 0 。 显然矩阵日是a 对s 的限制算子甲的矩阵表示，h 的特征值都是a 的特征 值，甲的特征值和特征向量都是a 的特征值和特征向量。 南京航空航天大学硕士学位论文 若s 不是一的不变子空间，q 仍记为以其正交规范基为列所构成的矩阵， 对任何m 矩阵曰，相应的剩余为 r ( b ) = a q q b a 一般说来，r ( b ) 的范数不再为零。显然，若i i r ( b ) 1 | ：越小，说明s 越接近 于a 的不变子空问，甲的特征对越接近于a 的特征对：反之，若l f r ( b ) f l ：很大， 说明s 远离a 的不变子空间，从而甲的特征对一般不是a 的特征对的较好近似。 关于剩余j ir ( b ) l i ：有如下结论： 定理2 1 2 1 对上述n m 列正交规范矩阵q ，则对任何m 埘矩阵口，有 l i 尺( 日) l h - 1 1r ( b ) i i ： 其中h = q 7 一q 。 定理2 1 2 说明当h = q 7 a q 时，剩余范数| l 尺( e ) l l ：达到最d 、。因此 r a y l e i g h r i t z 逼近是一种最佳逼近。 定理2 1 3 f 2 l 设a 是阶对称矩阵，其特征值为五，如，厶，q 是n 棚列 正交规范矩阵，口是m 阶对称矩阵，其特征值为“，“：，“，如果 r ( b ) = 彳q q b 则可以找到册个不同的自然数，f ：，使得 j 一1 1 女i 1 i r ( b ) j 1 2 定理2 。1 4 1 2 4 1 设工是对称矩阵4 对应于特征值五的特征向量，y 是x 的一 个d 研) 近似特征向量，则y 对a 的r a y l e i g h 商是五的d ( 矿) 逼近。 子空间迭代的两种r a y l e i g l 商加速技术 2 2 子空间迭代法 子空间迭代法可以看作是乘幂法的直接推广，与乘幂法不同的是子空间迭 代法是同时用几个( 例如p 个) 正交规范向量迸行类似于乘幂法的迭代，并在 迭代过程中应用r a y l e i g h r i m 原理进行加速。若将这些迭代向量看作p 维子空 间的正交规范基，则每迭代一次就得到一个新的p 维子空间。 设a 是阶可逆实对称矩阵，其特征值为 ，a z , - - , 如，且 五s 如墨0 以“墨墨氐 相应的标准正交特征向量为一，x ：，h 。假设计算a 的，个最小特征值及相应 特征向量，取p 为正整数，g r p n ，x o 是n x p 列正交规范初始矩阵， 则子空间迭代法可概述如下： 算法2 2 1 子空间迭代法 1 ) 随机选取p ( r p ) 个初始向量，将他们正交规范化后记为 x 坤) - 石f “，x p ，o 翔； 2 ) 对v = 1 , 2 ， 2 1 ) 由a x = x t ”，求出z “= 嘣”，x ：v ，? 】； 2 2 ) 将x ( 7 进行旦唬分解得：z ”= q 。r ， 其中q 是n x p 列正交规范矩阵，r ”是p 阶上三角矩阵； 2 3 ) 构造r a y l e i g h 矩阵b = q ”。4q ( ”； 2 4 ) 计算b ”1 的全部特征值“j ”及相应的标准正交特征向量 y j ”( f _ 1 , 2 ，p ) 令d ”= d i a g ( u ：”，“：”，“? ) ，y “；【y ：“，y ，y ? 】。 2 5 ) 计算r i t z 向量矩阵s “= p 9 y “= 耐”，5 ，s p “1 】； 2 6 ) 判断收敛性： 若j | 川5 ，s 卜时”，s 】d 峙 f ，则结束；否则x “1 】- s ”。 6 南京航空航天大学硕士学位论文 关于算法2 2 1 作如下说明 1 ) 通常随机选取p 个向量，然后用m g s 方法得到n p 列正交规范初始 矩阵x ( o ) 。 2 ) 可用q 吧方法或j a c o b i 方法计算p 阶对称矩阵占p 的全部特征值和相应 的标准正交特征向量。 3 ) 在2 1 ) 中，大型线性方程组丘r = j ”，可用c g 方法来求解。 子空间迭代法的收敛性 定理2 2 1 【2 4 1 若a 是对称正定矩阵，则由算法2 2 1 所形成的子空间序列 s ) 不退化，即s ”是p 维的。 定理2 2 2 t 2 3 1 设对称矩阵a 的特征值满足 丑以 只( _ ) ( ，= l 一印) 。 在一元二次方程( 3 1 1 ) 中， 糊q 一。响r ( m + 挚“器叫器卜 、 ，) ，。 r ( m ，) lr ( 州，) i 其中，= f 1 。易知，当，= 尺( j ，) = r ( m ，) 时，取最小值o 。 又因为 r ( s j ) 尺( ， 所以 0 ，即( 3 1 2 ) 必有两个不同的实根。 设( 3 1 2 ) 的两个根为口，岛。 因为 峭2 r 丽i k i 州( 1 - 等) o ， 6 - k ，( 1 觚p ( s j ) i ) 一口j 。 这与( 3 1 3 ) 矛盾，所以一a 不是a 的特征值，故b o 是非退化的。 因为 b 似= x ( ) r c 】 所以x ( “1 也是非退化的。证毕 定理3 1 3 假设x ( o 与a 的特征向量x 完全不正交，即 x t z ：o o ( j = 1 ，p ) ，则x 也与x 完全不正交。 证明采用归纳法。 由假设知z o 与x 向量完全不正交，即l x t x ：o p0 ( j = l 2 ，p ) 。现假设z 。 与x 完全不正交。很明显，m 。也与x 完全不正交， 即 妒m o ( j = l 2 ，p ) 。 假设 m i n ( x r 6 一i ，卜7 蜒”i ，卜t 咋( t i ) = l x 7 q o i ， 则有 南京航空航天大学硕 ：学位论文 x 7 i = lx t ( 一+ q ) 0 所以b 与x 完全不正交。 又因为b ( ) = x ( r “，所以x 1 也与x 完全不正交。证毕 定理3 1 4 在算法3 1 1 中，对于b 的每一列鳄和s 。的每一列 。? ( ，= 1 , 2 ，p ) ，有 r ( 巧。) r ( ) 。 证明为简便表示，我们将r ( q ) ，尺( j 简记为r ( b j ) ，r ( s ，) 。 因为 r ( q ) ：尝，且o ：叩，+ s j ， r ( q ) 2 亏了尹且o 2 口，脚， 所以 月( 屯) = 兰；觜，其中k j = m f a m j , _ = _ 7 一m ， 则 r c b ，) - r ( s j ) = 盟鬻糟型， ，固 两边同时对口，求导得 l r ( s 舰一_ i a ；n j 两- 1 ( 1 - r ( s j ) _ ) 将( 3 1 9 ) 带入( 3 1 8 ) 得： 由引理3 1 1 知， 由引理3 1 1 知， 所以 渺肌户等掣 l - r ( s j ) h j 0 a i 0 1 + a i h i ( 0 r ( b j ) 一r ( s ，) 0 ，即r ( 屯) 月( 0 ) 。 ( 3 1 9 ) 子空间迭代的两种r a y l e i g h 商加速技术 证毕 定理3 1 5 算法3 1 1 所产生的序列 1 ) ：。( = 1 , 2 ，p ) 是收敛的。 证明由定理3 1 4 知， 所以 r ( x ：) ) ：。是单调下降的。 因为 尺( q ) o ，b n o 。，忽略规范化因子， 有 巧，：芎，一杰碟，n ，弓。：q 一，一兰够棚6 f “， ；i ： 这里 彬。= i 秽f s ，= i l e ， 对所有得i n ，使得 这里 s ，= 矽- z 碟 吲i = x ， s ( 3 ，1 t o ) 定理3 1 6 假设z 与特征向量，x 。完全不正交， 则 l i m j ：”= x ，( = l 2 ，p ) 。 r 证明首先证明 受s 一= 。采用反证法 南京航空航天大学硕士学位论文 假设 骢s p _ ，则s “的任何列都不收敛到_ 。 我们忽略才，彰”，砖。的规范化网子，则 巧。= x l + ，+ t _ ，( = l ，2 ，p ) 这里k j 2 ，吃是与，叱，正交的单位向量。 由 1 4 ，我们知道，存在常数，口。，口。，使得 巧n = s j 一+ + t l ， 这里 t = ( 1 + q ) 0 ，o 口。口，s 。，i t 峰pj t i ， 由( 3 1 1 0 ) ，我们得 s ? “1 = t z + “。，+ t 矗，- z h ( 4 z 。+ “h + d _ ) = ( t 一h 巧) 一+ ( 1 一d ) “q 十一_ ， 令足够大，则 r ，面t 2 m i n l ，蔓1 ， 。 + 口m “+ p 。 对所有的i n ，令1 0i 。= s ，则 链翻= 坠铲 垫半避舁业 坐半岩出 习占，i ， 所以，妙“中而与“的系数比要大于。中x 一与“的系数比，因此，可能随着f 的变化而变化，但序列 s ! ) 有如下结论： i x i t s ：剧x i t 5 _ l ，i n 子空间迭代的两种r a y l e i g h 商加速技术 这与我们的假设矛盾，所以必有s ：。中的一个收敛于上。因为s “1 的列是按 r a y l e i g h 商由小到大的顺序排列的，所以，s p 收敛于置。同理，s ：。收敛于x ， ( = 2 ，3 ，p ) 。证毕 定理3 1 7 由算法3 1 1 所产生的迭代序列 秽) 刍单调收敛于五，且是渐进 二次收敛的。 证明由定理3 1 5 知， m 三是收敛的，又由定理3 1 6 我们得到埘 ：。 收敛于 。所以，对于充分大的k ，我们有 d ，= 五+ ( “，这里s ”是一个充分小的实数， 由极大极小值定理 2 4 1 1 2 7 1 1 2 8 1 知， d ，圳：m i 。m 。羔! 【茎! ：! 掣! ：! ! 坦， j v i ，o | y r x ”mxc k + 1 ) y 其巾矿是r ”中j 维子空间。 于是存在一个向量y = _ y 矿，使得 则 其中， d p “) = ! ；毒；等，令“) y 仕) = z 忙 。 v ( 1 ) 7z “彳( + 1 v ) 哆( k + 1 ) 卅- - ，+ 丝筹筹叫，+ 羔 口芦= z 7 ( 一一d p 1 ) z 由对称矩阵爿一巧，特征向量的扰动理论知，对充分小的k ，一个有界常数声” 使得 一= x ，+ p v 、， 其中x ，；l 毛k e r ( a 一，) 的一个单位向量，即 口一2 j ) x ，= 0 - r ) i x j1 1 = 1 。 s ，是正交于k e r ( a 一五，- 0 的一个单位向量，则 南京航空航天大学硕士学位论文 所以有 所以 故 瑾! ) = z ( ) 7 ( 一一2 j i ) 工( n 一占( ) z 耻) 7 z ( ：一 ( 7 + 0 ( s ( ) 2 ， z 7 z - - - t lx ，畦+ d ( s ) 2 ， 巧“”= 彰虬s ”+ o ( s ) 2 = t + d ( 占似1 ) 2 s = o ( e 、2 。 3 2 数值实验 本节给出了应用子空间迭代法的r a y l e i g h 商加速技术计算对称正定矩阵 a 的若干个最小特征值的数值例子。例子中的初始矩阵都是随机产生的。当两 种方法进行比较时，都取相同的精度要求。以下例子中，r 表示我们要求特征 值的个数，p 表示初始矩阵的列数，s 表示r i t z 值d ：o 与相应r i t z 向量s ，作为 的特征对时其残量范数的误差界。i t n 代表算法的迭代数，c p u 代表所需要的 计算时间，单位为秒。s u b s p 代表子空间迭代法，子空间方法的r a y l e i g h 商 加速技术的m a t l a b 程序名用r m i n 表示。 例3 2 1矩阵a 是1 1 0 阶对称正定矩阵 4一lo 14一l l+ 4一l o一14 r = 3 ，p = 6 ，s = 1 0 e 一5 ，用r a y l e i g h 商加速技术和子空间迭代法求a 的三个最 小特征对，其结果见表1 1 。 子空间迭代的两种r a y l e i g h 商加速技术 表1 1 例3 2 1 的计算结果表 程序名特征值i t nc p u f s l s u b s p2 0 0 0 8 0 0 9 8 5 7 6 5 7 56 7 21 3 5 6 2 0 0 3 2 0 3 3 0 1 6 0 1 3 6 2 0 0 7 2 0 5 0 2 2 9 4 9 3 7 r m i n2 o 0 0 8 0 0 9 8 5 7 6 3 3 81 2 60 5 4 7 2 0 0 3 2 0 3 3 0 1 4 7 6 0 5 2 0 0 7 2 0 5 0 2 2 9 3 6 2 3 例3 ，2 2 矩阵a 是1 1 0 0 阶对角矩阵d i a g ( 1 0 1 ，i + o ，0 1 ，1 1 0 0 0 1 ) ( f = 2 ，1 1 0 0 ) 。r = 3 ，p = 6 ，s = 1 0 e 一5 ，用r a y l e i g h 商加速技术和子空间迭代 法求a 的三个最小特征值，其结果见表1 2 衷1 2 例3 2 ，2 的计算结果表 程序名特征值1 1 n c p u ( s ) s u b s p1 0 l 0 0 0 0 0 0 0 0 0 0 0 01 65 9 8 5 2 0 1 0 0 0 0 0 0 0 0 0 0 0 0 3 o 0 0 0 0 0 0 0 0 0 0 4 0 0 r m i n1 0 10 0 0 0 0 0 0 0 0 0 0 01 33 6 4 0 2 0 1 0 0 0 0 0 0 0 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0 0 0 0 2 1 例3 2 3矩阵a 是1 2 0 0 阶对称正定阵 l0 50 5 0 520 5 0 530 5 o 5 0 5 一lo 5 o 5n r = 3 ，p = 6 ，占= 1 o p 一5 。用r a y l e i g h 商加速技术和子空间迭代法求a 的三个最 小特征对，其结果见表1 3 南京航空航天大学硕士学位论文 表1 3 例3 2 3 的计算结果表 程序名特征值i t nc p u ( s ) s u b s p0 7 7 4 3 9 2 8 7 1 9 4 0 0 31 71 1 7 5 0 1 9 7 6 4 9 8 8 9 3l8 5 9 6 2 9 9 8 9 2 3 8 0 4 9 1 0 7 0 r m i n0 7 7 4 3 9 2 8 7 1 9 4 0 0 31 38 0 4 7 1 9 7 6 4 9 8 8 9 3l8 5 9 6 2 9 9 8 9 2 3 8 0 4 9 1 0 3 8 由以上数值例子可见子空间方法的r a y l e i g h 商加速技术比子空间迭代法 收敛速度快，并且节省了计算量和计算时间，特别当要求的特征值密集时效果 更明显。因此，子空间迭代的r a y l e i g h 商加速技术是求解大型稀疏对称矩阵特 征值问题的有效方法之一。 子空间迭代的两种r a y l e i g h 商加速技术 第四章带位移的子空间r a y l e i g h 商加速技术 4 1 带位移的子空间r a y i e i g h 商加速技术 r a y l e i g h 商迭代法( r q i ) 1 3 7 1 1 ”i 是反迭代m p s l l ”1 的一种重要变形，也可看件 是带有变原点位移的反迭代方法，其原点位移的当前值取为当前迭代向量的 r a y l e i g h 商。r q 法的迭代格式如下： 算法4 i 1 r q i 迭代法 1 ) 取初始向量j ”t ，毛o 是膏o 的r a y l e i g h 商，即o = r ( x o ) ； 2 ) i = 0 ，l ，2 ，执行 2 1 1 求x ( ) 似一1 ) x 。= o x “，这里仃是规范化因子。 2 2 ) 计驯筹 带位移的反迭代方法，它的优点是收敛速度快，尤其是当a 是对称矩阵时， 特征值是按立方收敛1 8 1 1 2 4 3 0 1 的。但它最后迭代所得到的特征值与他所选的初始 位移有关，要想得到最小的特征值，则初始向量必须充分接近最小特征向量。 子空问方法虽然收敛速度慢，但却能收敛到最小的特征对，因此可以将两者结 合起来，当a 是可逆的是对称矩阵，先进行若干步子空间迭代，当迭代向量接 近最小特征向量时，再用r q i 方法来加速，得到如下的算法： 算法4 1 2 带位移的子空间r a y l e i g h 商加速技术 】) 执行m 步算法2 i 】，得列正交规范初始矩阵 _ m = ( x f “z ? ) ( r p j v ) ； 2 ) 对i = 0 , l ， 2 1 ) 计算以，；x 叩船”，求4 的全部特征值d f 彰1 和相应的特 征向w f ”，w ：记d “= ( 打口g ( d ；”，d ) ，缈7 ；( w ：”，j ：) 。 南京航空航天大学硕士学位论文 计算s ( 1 ) = z ( 。) w c ) ； 2 2 ) 检验收敛性 若l l 爿 s f ”，”卜【s f ”，o 】d “l i f s ， 则斫。彰“，s s ! ，o 取为所求的特征值及相应特征 向量：否则执行下一步。 2 3 ) 求解下列线性方程组 4 m “一m 。d 。= s ”，其中d 。= d f ：a g ( a ：“，彳? ) ， 即 ( 爿一穆1 ，) m ；= ( j = l ，p ) 解得m “) _ ( 耐”，m 妒) ； 2 4 ) 构造b ( f ) = 埘嗣a ( + s f o 其中a o = 讲昭( 口 ”，口：) 为p 阶的对角矩阵， 即6 j 。= 蟛。脚罗+ 巧( = l ，p ) ，其中彰。使p ( 6 ；) 达到最小。 2 5 ) 将b ( 。进行q r 分解： b ( 。) = x ( “) 尺“) 其中( “为列正交规范阵，j r 为上三角矩阵。 关于算法4 1 2 作如下说明： ( 】) 这里p 的选取与算法3 1 1 中的p 的选取方法是一样的。在算法4 1 2 的2 1 ) 中求解a ，的特征值和特征向量同样也可用j a c o b i 或q l 方法来 求解。在2 3 ) 中，同样也采用不精确的方法来求解。 ( 2 ) 算法4 1 ，2 与算法3 1 1 主要步骤是相同的，不同之处在于2 3 ) 。前者 采用带位移的反迭代来得到迭代向量，而后者则直接采用反迭代方法。 ( 3 ) 经过大量实验证明，算法4 1 2 中肌通常取为算法2 1 1 迭代步数的 ( 4 ) 关于参数矩阵a 的选取问题 在2 4 ) 中， 蟛o = 口0 1 m + 5 ( ，= 1 , 2 ，p ) ， 在不引起混淆的情况下，省略上标简记为 子空间迭代的两种r a y l e i g h 商加速技术 b j = a 卅j 七5 j ( ，= 1 , 2 ，p ) 。 因为 妙喾= 每鬻 将= l ，一d j i ) m j = 0 代入上式得： 删= 端， 其中 ，= t ， r = l + 矿m y $ ， t 一- - 。t j m i 七d mj t f = 一s i t p ；m r m - 对( 4 1 1 ) 式关于口求一阶导数，并令其等于零，得 f a2 + g a + h = 0 其中 ( 4 1 1 ) ( 4 1 2 ) f = t f r e = ( s j m ) 2 m s m ， g = f 一e = ， h = ，一矿= 1 ， 定理4 1 1 一元二次方7 f 呈( 4 1 2 ) 有两个异号实根，且当口取负根时，( 4 1 1