（计算数学专业论文）求反对称三对角矩阵特征值问题的分治算法研究.pdf

大连理工大学硕士学位论文 摘要 本文研究了求反对称三对角矩阵特征值问题的分治算法把反对称三对角矩阵的特 征值问题转化为对称三对角矩阵的特征值问题，避免了复数运算，减少了计算量并结 合对称三对角矩阵的结构，给出基于不同分治策略的分治算法文章主要包括三部分： 第一部分介绍了矩阵的不同分割策略、l a g u e r r e 迭代的性质以及如何利用l a g t l e l t e 迭代处理病态特征值束问题，最后介绍了迭代的停止标准 第二部分在廉庆荣教授等给出的求反对称三对角矩阵特征值理论的基础上，给出了 一种基于秩1 扰动的分割策略和l a g u e r r e 迭代的分而治之算法该算法较适合求按模较 大的那些特征值，而对于按模较小的特征值可能失去有效数字 第三部分讨论了反对称三对角矩阵与其相伴矩阵( 对角元为零的特殊对称三对角矩 阵) 特征值之间的关系给出了一种基于秩2 扰动的分割策略和l a g u e r r e 迭代的分割胶合 算法。这样做的好处是保持了矩阵的特殊结构和特征值的正负成对出现，实际计算中只 需计算其中的非负特征值即可，减少了计算量 分治算法具有速度快和能根据不同需求而灵活实现的良好性质它既能求矩阵的全 部特征值，又能求指定的若干个特征值或给定区域内的特征值值得指出的是分治算法 具有自然的并行性，适用于求大规模矩阵特征值问题，是一种富有前景的求矩阵特征值 的实用方法 关键词：反对称三对角矩阵；分治策略；分而治之算法；分割胶合算法；相伴矩阵；特 征值；l a g u e rr e 迭代 王文超：求反对称三对角矩阵特征值问题的分治算法研究 t h er e s e a r c ho fd i v i d ea n dc o n q u e r a l g o r i t h m sf o rs k e w s y m m e t r i c t r i d i a g o n a le i g e n v a l u ep r o b l e m s a b s t r a c t t h i st h e s i si n v e s t g a t e sd i v i d e a n d - c o n q u e rm e t h o df o re i g e n v a l u ep r o b l e m so fs k e w - s y m m e t r i ct r i d i a g o n a lm a t r i c e s ，w h i c hr a m st h ee i g e n v a l u ep r o b l e m so fs k e w - s y m m e t r i c t r i d i a g o n a lm a t r i xt ot h ee i g e n v a l u ep r o b l e m so f s y m m e t r i c 心i d i a g o n a lm a t r i x , w h i c ha v o i d st h e c o m p l e xo p e r a t o na n dr e d u c e st h ec o m p u e c a t i o na m o u n t c o n s i d e r i n gt h ed i v e r s es t n l c 缸r e so f t h es y m m e n c 城d i a g o n a lm a t r i c e s ，o n ep r e s e n t sv a r i o u sd i v i d e a n d c o n q u e ra l g o r i t h m sb a s e d o nd i f f e r e n td i v i d e a n d - c o n q u e rs t r a t e g i e s t h ea r t i c l ei n c l u d e st h r e em a i n p a r t s ： t h ef i r s tp a r ti n t r o d u c e sv a r i o u sd i v i d e - a n d c o n q u e rs t r a t e g i e s ，p r o p e r t i e so fl a g u e r r e i t e r a t i o na sw e l la sh o wt od e a lw i t ht h ec l u s t e r so fp a t h o l o g i c a ln g e n v a l u e sb yl a g u e r r e i t e r a t i o n 。a tl a s t , o d eg i v e st h es t o p p i n ge r i t e t a b a s e do nt h et h e o r e mf o u n d a t i o no fs k e w - s y m m e t r i ct r i d i a g o n a lm a t r i xw h i c hg i v e nb y p r o f e s s o rl i a nq i n g r o n ge ta l ，t h es e c o n dp a r tp r e s e n t sad i v i d e a n d - c o n q u e ra l g o r i t h mt h a ti s o nt h eb a s i so ft h er a n k - o n ep e r t u r b a 矗0 1 2o fm a t r i xa n dl t g u e r r ei t e m f i o mn l ea l g o r i t h mi s s u i t a b l ef o re a l c u l a t et h o s el a g e re i g e n v a l u e sb ym o d u l o b u ti tm a yl o s et h es i g n i f i c a n tf i g u r e s i nc a s eo f t h o s es m a l l e re i g e n v a l u e sb ym o d u l o n l et h i r dp a r td i s c u s s e st h ee i g e n v a l u e sr e l a f i o n s u pb e t w e e ns k e w - s y m m e t r i c 趟d i a g o a a l m a t r i xa n di t sc o n c o m i t a n tm a t r i x ( s p e c i a ls y m m e t r i cd i a g o n a lm a t r i xw i 也d i a g o n a le l e m e n t s a r ez e r o ) o n eg i v e sa s p l i t - m e r g ea l g o r i t h mf o rt h ec o n c o m i t a n tm a t r i xb a s e do nt h er a q k - t w o p 嘶b a f i o no fm a t r i xa n dl a g u e r r ei t e r a t i o n i t sa d v a m a g e sa r et ok e e pt h em a t r i c e ss t r u c t u r e a n dp o s i t i v e - n e g a t i v ee i g e n v a l u e si np a i r a n dp r a c t i c a lc o m p u t a t i o np r o c e s sn e e d so n l y c o m p u t e rn o n n e g a t i v ee i g e n v a l u e s ，a n dt h a tc u t sd o w nt h ec o m p u t a t i o na m o u l l t d i v i d e - a n d - c o n q u e rm e t h o dp o s s e s s e st h eg o o dc h a r a c t e r i s t i cw h i c hi ss p e e d f a s ta n d f l e x i b l eu r i n g i n go u t , w h i c hc a nn o to n l ys e a r c hf o ra l lt h ee i g e n v a l u e s b u ta l s og e ts o m ef i x e d e i g e n v a l u e so re i g a n v d u e si n 垂v e ni n t e r v a l k sw o r t hn o t i c i n gt h a to u ta l g o r i t h md o e sw e l li n i t sn a t u r a lp a r a l l e l i s mt h a ti ss u i t a b l ef o rs o l v i n gt h ee i g e n v a l u ep r o b l e m so f l a r g es c a l em a t r i c e s a b o v ea l l ，t h em e t h o di sa ne i g e n v a l u e s e a r c h i n gm e t h o dw i t hw i d e s p r e a dp r o s p e c t k e yw o r d s ：s k e w - s y m m e t r l ct r i d i a g o n a lm a t r i x ；d i v i d e - a n d c o n q u e rs t r a t e g y ； d i v i d e - a n d - c o n q u e ra l g o r i t h m ；s p l i t - m e r g ea l g o r i t h m ；c o n c o m i t a n tm a t r i x ；e i g e n v a l u e ； l a g u e r r ei t e r a t i o n 独创性说明 作者郑重声明：本硕士学位论文是我个人在导师指导下进行的研究工作 及取得研究成果尽我所知，除了文中特别加以标注和致谢的地方外，论文 中不包含其他人已经发表或撰写的研究成果，也不包含为获得大连理工大学 或者其他单位的学位或证书所使用过的材料与我一同工作的同志对本研究 所做的贡献均已在论文中做了明确的说明并表示了谢意 vo1 ，一” 作者签名：至翌日期：! ! ! ! = 兰 大连理工大学硕士研究生学位论文 大连理工大学学位论文版权使用授权书 本学位论文作者及指导教师完全了解“大连理工大学硕士、博士学位论文版权使用 规定”，同意大连理工大学保留并向国家有关部门或机构送交学位论文的复印件和电子 版，允许论文被查阅和借阅。本人授权大连理工大学可以将本学位论文的全部或部分内 容编入有关数据库进行检索，也可采用影印、缩印或扫描等复制手段保存和汇编学位论 文。 作者答勿，茗签名：形生冬竺 靴龇狲蔟 导师签名：z 皇蔓2 1 1 大连理工大学硕士学位论文 1 绪论 1 1 问题的背景及来源 数值线性代数是一个非常活跃的研究领域除了其自身的许多问题时刻在促进着它 的繁荣和发展，许多依靠数值线性代数的科学计算也在推动着这一领域的繁荣和发展不 仅许多传统的科学计算，如物理模型、工程问题，依赖于数值线性代数的计算，而且许 多现代应用学科，如信息恢复、图像还原，也都得益于数值线性代数代数特征值是数 值线性代数的一个重要研究领域之一，对它的研究具有重要的理论意义和实用价值许 多科学和工程问题，如结构力学中的固有频率分析问题以及控制系统中的稳定性问题， 最终都转化为特征值问题 自从1 9 4 6 年第一台计算机e n i a c 问世，经过半个多世纪的发展，使得科学工程计 算已成为当今世界最重要的科学进步之一计算物理学家、诺贝尔奖获得者w i l s o n 在二 十世纪八十年代指出：科学计算、理论研究、科学实验并列为当今世界科学活动的三种 主要方式许多科学和工程领域如果没有科学计算不可能有一流的研究成果 矩阵计算是科学与工程计算的核心，可以毫不夸张地讲，大部分科学与工程问题都 要归结为矩阵计算的问题，其中具有挑战性的问题是大规模矩阵的计算问题 矩阵计算的基本问题 ( 1 ) 求解线性方程组的问题即给定一个n 阶非奇异矩阵a 和n 维向量6 ，求一个抖 维向量z ，使得 a x = b( 1 1 1 ) ( 2 ) 线性最小二乘问题，即给定一个m n 阶矩阵a 和m 维向量b ，求一个”维向量 x ，使得 j 恤一b l l = m m l l a y - b l l ，y e 彤) ( 1 1 2 ) ( 3 ) 矩阵的特征问题，即给定一个聆阶实( 复) 矩阵a ，求它的部分或全部特征值 以及对应的特征向量，也就是求解方程 a x = a x( 1 1 3 ) 一对解，z ) ，其中z r ( c ) ，x r ”( ) 即五为矩阵a 的特征值，工为矩阵爿的属于特 征值五的特征向量 在工程上，矩阵的特征值具有广泛的应用，如大型桥梁或建筑物的振动问题：机械 和机件的振动问题；飞机机翼的颤振问题；无线电电子学及光学系统的电磁振动问题； 王文超：求反对称三对角矩阵特征值问题的分治算法研究 调节系统的自振问题以及声学和超声学系统的振动问题又如天文、地震、信息系统、 经济学中的一些问题都与矩阵的特征值问题密切相关 在科学上，计算流体力学、统计计算、量子力学、化学工程和网络排队的马尔可夫 链模拟等实际问题，最后也都要归结为矩阵的特征值问题 由于特征值问题在许多科学和工程领域中具有广泛的应用，因此对矩阵的特征值问 题的求解理论研究、算法的开发、软件的制作等是当今计算数学和科学与工程计算研究 领域的重大课题，国际上这方面的研究工作十分活跃 1 2 矩阵的特征值问题研究现状及算法概述 对一个 玎阶实( 复) 矩阵一，它的特征值问题，即求方程( 1 1 3 ) 式的非平凡解， 是数值线性代数的一个中心问题这一问题的内在非线性给计算特征值带来许多计算问 题为了求( 1 1 3 ) 式中的五，一个简单的想法就是显式地求解特征方程 d e t ( a - 2 1 1 = 0 ( 1 2 1 ) 除非对于个别的特殊矩阵，由于特征方程的系数不能够用稳定的数值方法由行列式 的计算来求得，既使能精确计算出特征方程的系数，在有限精度下，其特征多项式 ，( 五) = d e t ( a 一五，) 的根可能对多项式的系数非常敏感因此，这个方法只能在理论上是 有意义的，实际计算中对一般矩阵是不可行的首先，若矩阵a 的阶数较大，则行列式 d e t ( a 一五，) 的计算量将非常大；其次，根据g a l o i s 理论，对于次数大于四的多项式求根 不存在一种通用的方法，基于上述原因，人们只能寻求其它途径因此，如何有效地、 精确地求解矩阵特征值问题，就成为数值线性代数领域的一个中心问题 目前，求解矩阵特征值问题的方法有两大类：一类称为变换方法，另一类称为向量 迭代方法变换方法是直接对原矩阵进行处理，通过一系列相似变换，使之变换成一个 易于求解特征值的形式，如j a e o b i 算法。】，g i v e n s 算法吣”1 2 7 】，q 幔算法m 5 1 等变 换方法由于要存储矩阵元素，因而它只适用于求解中小型矩阵，它一般和向量迭代方法 结合起来使用向量迭代方法是通过一系列矩阵向量乘积而求得特征值和特征向量的由 于向量迭代方法可采用压缩存储技术，因而它适合于求大规模矩阵的特征值问题，尤其 是大型稀疏矩阵的部分特征值和特征向量问题，如l a n c z o s 方法口，6 】，a r n o l d i 方法【4 】， d a v i d s o n 方法川等，现在这类问题仍是比较热的研究课题 到1 9 7 0 s 后期，随着中型稠密及大型矩阵特征值的研究，尤其是对称矩阵，人们迫 切需要更加有效地、更加精确地求矩阵特征值的计算方法在1 9 8 1 年c u p p e n ”3 给出求 解对称三对角矩阵特征问题的分而治之算法，并由d o n g a r r a 和s o r e n s e n 一1 于1 9 8 7 年成 大连理工大学硕士学位论文 功地把这一思想用于对称三对角矩阵特征值的并行计算随后，越来越多学者关注这一 算法的研究，并取得了不斐的成果。4 。” 1 3 本文的工作及有关符号 分而治之方法的提出，打破了以往人们认为q 幔算法是求中小型矩阵特征值的最有 效的方法的思想d e r r n n e l 在文 1 8 y 9 提到，计算对称三对角矩阵所有特征值和特征向量， 当矩阵的阶数大于2 5 时，分治算法比q r 算法快值得指出的是这种方法也适用于求大 规模矩霹特征值问题 本文在第二章介绍了矩阵的分割策略和l a g u e r r e 迭代的性质，第三章介绍了改进的 求反对称三对角矩阵特征值模平方的分而治之算法；第四章介绍了反对称三对角矩阵的 特征值与其相伴矩阵的特征值之间的关系，并给出求反对称三对角矩阵特征值的另一分 治算法本文的主要工作分为以下两个部分： ( 1 ) 基于已有理论基础，对反对称三对角矩阵，的叉积7 的压缩矩阵采用分治方 法，并采用矩阵的秩l 扰动，使得特征区间只包含个特征值，在每个包含区间寻求适 当的初始点，使迭代快速收敛于j 的特征值的虚部的平方次幂，然后对其开平方求出相 应的共轭特征值的虚部 ( 2 ) 证明了反对称三对角矩阵与特殊的对称三对角矩阵( 称之为前者的相伴矩阵) 的特征多项式、特征值以及特征向量之间的关系，把求反对称三对角矩阵的特征值问题 转化为其相伴矩阵的特征值问题，采厢秩2 扰动对其相伴矩阵作扰动( 为了保持相伴矩 阵特征值正负成对的性质) ，给出了求反对称三对角矩阵特征值的分割胶合算法 对于上述两种分治算法，由于其内在的并行性，如何设计出更加高效的并行计算是 有待于进一步解决的问题 在本文中用r “和c ”分别表示n 维实的和复的向量空间：彳表示n x n 阶实和复矩阵； 五表示对称三对角矩阵；表示k 阶单位矩阵；上标7 。表示矩阵或向量的转置：气表示 第k 个单位坐标向量：a e t ( a ) 表示矩阵一的行列式；r ( c “”) 为实( 复) 数域上的h m 阶矩阵全体：如无特殊说明| | o | | i o 表示矩阵或向量的e u c l i d e a n 范数 王文超：求反对称三对角矩阵特征值问题的分治算法研究 2 分治方法的基础及理论研究 2 1 分治方法的概述 考虑对称三对角矩阵五的特征值问题 其中 写= l z = l x ， 届 a 口2 孱。 ( 2 1 1 ) ( 2 1 2 ) 1 9 8 1 年c u p p e n t 8 】提出一种求上述对称三对角矩阵所有特征值和特征向量的分而 治之方法( d i v i d e a n d - c o n q u e rm e t h o d ) 其基本思想是先将对称三对角矩阵l 分割为两个 分别为| j 七阶和印一_ ) 一后) 阶低阶对称三对角子矩阵丁( o 和丁( ”r ( o 和，( ”可以用 同样的方法也分别分割为两个更低阶的子矩阵，递归的采用这种分割技术可以把矩阵分 割为一些能嘉接求出特征值的足够小的子矩阵( 比如2 阶或1 阶矩阵) ，或者按照某种 标准分割到适当阶数( 如小于等于2 5 阶) 后，结合其它求矩阵特征值的方法，如q 尺算 法，求出其特征值在求出低阶矩阵特征值的基础上，开始胶合过程在胶合阶段，分 割前的矩阵r 的特征值的求出( 所谓的“治之”) 是建立在其两个子矩阵丁( 0 7 和t c w 的 特征值的基础上的，其中丁( 0 7 和r ”是在分割阶段由r 分割出的低阶子矩阵随后的数 值分析表明，c u p p e n l 8 1 的方法存在着数值不稳定的危险，特别是当存在特征值束时，计 算出的特征向量可能不正交g u 和e i s e n s t a t “1 对c u p p e n ”的方法作了改进，极大地 降低了数值不稳定的危险性 c u p p e n 8 1 的方法在计算的特征值的同时也需要计算对应的特征向量，并且是在 r ( 0 和r ( 1 的特征值和特征向量的基础上进行计算的根据文 1 5 ，当用残量忆z 一 肖0 和正交性x 7 x l | | 作为检验准确性的标准时，c u p p e n t s l 的方法比二分法或多分法精确 i i 的多但根据文【1 3 】，如果用恤一丑| l 作为衡量特征值准确性的标准时，二分法或多分 o 。 法精确些此外c u p p e n 的分而治之方法要求矩阵乘积，存储量为o ( n 2 ) ，而二分法或多 分法的存储量仅为o ( n ) ，比前者少应此，当只需计算特征值时，通常选用后者1 9 8 7 大连理工大学硕士学位论文 年d o n g a r r a 和s o r e n s e n 剀把分治思想应用到求对称三对角矩阵特征值的并行计算，取得 了不错的效果，再次引起了人们对分治方法的极大关注 分割胶合方法1 1 3 1 ( s p l i t - m e r g em e t h o d ) 是后来提出用分治方式求对称三对角矩阵特 征值c 的方法，不同于分而治之方法在分割过程中采用矩阵的秩l 扰动，它采用矩阵的 秩2 扰动与二分法和多分法相似，在分割胶合方法中，特征值的计算独立于特征向量 的计算如果需要计算特征向量，可采用反幂法 i , 2 71 6 】由于分割胶合方法计算特征值时 采用具有三次收敛的l a g u e r r e 迭代，数值试验表明，其计算速度和精度都明显优于二 分法和多分法并且文 1 3 1 中给出的用于计算l a g u e r r e 迭代的非线性三项递归式可以避 免上溢和下溢问题 2 1 1 分割策略 分治方法的第一步就是把原来高阶的对称三对角矩阵特征值问题转化为两个低阶对 称三对角矩阵特征值问题，即所谓的分割阶段设对称三对角矩阵z 表示如下 乃= r 届 以 。 屈一 尼一1 ( 2 1 1 1 ) 不妨假设岛o ( ，= 1 ，2 ，n 一1 ) ，即称不为不可约矩阵否则，若存在某些局= o ，则瓦 就可以约化为若干个低阶主子矩阵特征值问题，l 的特征值就由这若干个低阶主子矩阵 的特征值构成当l 为不可约对称三对角矩阵时，对其作如下分割 = r 丁 + e ， ( 2 1 1 2 ) 记疋= e + e ，其中t o 和丁1 分别为k x k 阶和一后) 一_ j ) ( 或0 一k 一1 ) ( n k 一1 ) ) 阶对称三对角矩阵，通常k = n 2 】 ( 1 ) 当e pp e o e p o z 称此为秩1 扰动矩阵此时有 0 = v 7 时，其中v = ( o ，1 ，臼，0 ) 且p 口= 屈 王文超：求反对称三对角矩阵特征值问题的分治算法研究 丁( o ) = ( 2 ) 当e = 矩阵此时有 r ( o ) 丁( o ) ( 3 ) 当e = q届 届 层 屈口2 0 o op p0 0 屁一。 a k p ，r ( 1 ) = 嘶+ l p 日2 鼠+ 。 。 屈一 孱一。 = p e 女e k + 1 7 + 尸“e k 7 时，其中p = 屏，称此为秩2 扰动 q届 届口2。 屈一 p k da i 8 t 6 吒+ l 屈+ 层+ l o 尺。“，1 ) ，丁( 1 ) = 吼+ 1 屈+ 1 屈+ l c + 2 p n 。 时，称此为秩3 扰动矩阵此时有 屈+ ： + 3 。 屈一， 屈一。 人们关注的问题是：对于上述三种分治策略与的特征值之间的关系利用 h o f f i n a n - w i e l a n d t 定理m 1 的结论，我们可以得到如下定理： 定理2 1 设彳和口是两个n 阶h e n n i t e 矩阵，它们的特征值分别是 五乃和 “鸬以，贝0 1 【o 乃- , u j ) 2 】_ - - l ，则有 下面关系成立 大连理工大学硕士学位论文 k 一五占 ( 2 114 ) 因此，只要( 2 1 1 4 ) 式右端越小，五，就越接近互，即t 与露的特征值就越接近对于 秩1 和秩2 扰动，正与的特征值之间关系更详细的描述，将在后两章节中给出，对于 秩3 扰动的详细的介绍可参见文 2 “ 2 1 2 胶合 在矩阵分割后，先求出矩阵的特征值，即求出，( o 和丁( 1 的特征值剩下的工 作就是如何由的特征值出发求出l 的特征值，这就是胶合阶段主要任务在这个阶段 可以采用不同的迭代方法，如割线法、n e w t o n 迭代、l a g u e r r e 迭代以及路径跟踪等， 以的某个特征值或的特征值构成的区间内的点为初始点经过若干步迭代，最后收敛 到瓦的某个特征值本文中采用l a g u e r r e 迭代从特征区间提取特征值 2 2 l a g u e r r e 迭代 2 2 1 l a g u e r r e 迭代及其计算 根据文【2 7 ，对于不可约对称三对角矩阵c ，它的特征值或特征多项式，) = d e t ( 一五l ) 的根全为互不相同的实数因此，适合求多项式的根都是实单根的具有三次收 敛的l a g u e r r e 迭代州 三f x ) ：x + 1 ；：；：兰二：；：；：一 ( 2 2 1 1 ) “。 ( 一器小_ 1 ) 邮( 一器) 2 _ n ( 舞) 】 非常适合用来从特征区间提取出瓦的特征值这个方法早在1 3 世纪就已经提出1 2 7 。“， 近年来结合分治策略又重新焕发出生机文 1 3 首次对用l a g u e r r e 迭代求不可约对称三 对角矩阵的特征值的实用性进行了研究，而后又被用于求矩阵的奇异值问题【1 ” 为了利用l a g u e r r e 迭代求z 的特征值，我们需要计算的特征多项式 f ( 2 ) = d e t ( t 一五厶) 以及其一阶导数厂研) 、二阶导数厂”( 五) 由文【2 7 】知，计算这些值 的有效的方法是三项递归式众所周知，当矩阵阶数比较大时，三项递归式可能出现上 溢或下溢问题文 1 3 给出了一种改进的非线性的三项递归式，数值试验表明，除极其个 别的情形，改进的三项递归式可以避免上溢或下溢问题 对称三对角矩阵正的特征多项式及其一、二阶倒数的计算公式 王文超：求反对称三对角矩阵特征值问题的分治算法研究 设对称三对角矩阵如( 2 2 1 1 ) 式所示，瓦表示其k 阶顺序主子式，表示如下 正= 崩 届 ( 2 2 1 2 ) 五( 且) = d e t ( t k 一4 i d 为霉的特征多项式，则特征多项式及其导数计算公式如下( 三项递 归式) ： j f o ( 五) = 1 ，石( 五) 2 喁，( 2 2 1 3 ) 【a ( j t ) = ( q 一五) 丘，( i t ) 一属：。正：( ) ，i = 2 ，3 ，- - ，胛 石) = o ，z + ( 兄) = 一，( 2 2 1 4 ) l f ( 五) = ( 一兄) z - ( 五) 一五。( 五) 一羼，f - ：( 五) ，i = 2 ，3 ，- ，” j 石( 兄) = o ，z ( z ) - l ，( 2 n 5 ) l f ( 五) = ( q 一五) ，：。( 丑) 一2 正。( 五) 一属! ，f ：( a ) ，i = 2 ，3 ，一，行 及 ，( 旯) = z ( 兄) ，( 五) = ( 五) ，f ”( 五) = f x 五) 为了避免上述( 2 2 1 _ 3 ) 式、( 2 2 1 4 ) 式和( 2 2 1 5 ) 式在计算中可能出现的上溢或下 溢问题，文 1 3 】给出如下的改进的非线性三项递归式： 令 黝) = 器，h ，2 ，m ( 2 2 1 3 ) 式两边都除以矗，( 五) 得 f 卣) = q 一五 地卅一彘，瑚，。， 2 1 6 令 纵栌一勰和弼) = f 删( 4 ) ，圳，l ，m ( 2 2 1 4 ) 式和( 2 2 1 5 ) 式两边除以z 一。( 五) 并注意下述关系 z ( 五) - - h 彭( 五) ， 大连理工大学硕士学位论文 瞄篱小。瓣 ，f 吃玢一m 泣m卜卜去卜以，弛c m 卜c 蒜，c 五， 乏a ，陀 知( 丑) = o ，g l ( - ) = 0 删= 赤卜佻卅2 钆阶c 急卜z 以“2 2 1 8 分别对( 2 2 1 7 ) 式和( 2 _ 2 1 8 ) 式进行简单的推导，我们可以得到下面两个关系式成立 一器= 删，锵吲饥 我们称( 2 2 1 6 ) 式为改进的三项递归式 ( 2 2 1 7 ) 式和( 2 2 1 8 ) 式分别为改进的计 算特征多项式一阶导数和二阶导数的计算公式，这些是我们在调用l a g u e r r e 迭代时所需 要计算的根据文 1 3 1 知，除了个别极其特殊的情形，它们可以避免计算过程中出现的上 溢或下溢现象但是，值得注意的是，如果对于某个i o i n ) 使得喜( 五) = 0 时，上述改 进的各式就会出现中断现象为了保证排除中断的发生，在计算中作如下处理 。如果专( 五) = 0 ，贝0 令毒( 五) = i 尼g i = 1 ，2 ，- ，”一l ； o 如果磊( 五) = 0 ，则令己( 五) = l 屏一s 1 这里占为机器精度上述处理相当于对磁或引入一个小扰动僻占i 或i 展一，f f 注2 1 ：我们由三项递归式( 2 2 1 3 ) 式经过改进得到( 2 2 1 6 ) 式，在计算改进的非线 性三项递归式( 2 2 1 6 ) 式时，同时可以得到在某个给定值 处，不大于该固定值五的特 征值的个数，即毒( 五) 0 ( 对所有1 i 兰 ) 的个数 根据上述介绍，我们给出计算特征多项式及其导数的改进三项递归式的算法描述： 算法1 改进的三项递归式( m o d i f y d e t b v a l ) 输入：t = 【屈。q ，屈】，旯 姚仉一器舻锵( 这默栌d e t ( t - - i ) ) 硼= n e g 一删脚 ( 后) 表示r 的小于五的特征值的个数) b e g i na l g o r i t h m1 曼望墨j 堕塑堑墅堑堂笙篁蛔整箜坌亟薹麴塞 卣= q a ，如果磊= o ，贝0 令轰= 群占2 。0 ，编= ，氏= 0 ，白= 0 ，n e g c o u n t = 0 如果轰 0 ，贝9 h 曙一c o u n t = 1 f o ri = 2 ：刀 耙叩：( 肇) 专= q 一五一把，印如果缶= o ，则令毒：( 譬2 = l ) s ： 如果眚 o 贝q n e g c o u n t + = 1 r , = ( q 一五) 卡珥一l + 1 一t e m p + 矾一2 白= 【( q 一兄) + 一l + 2 r l , 一i t e m p + 每一2 】 e n df o r e n da l g o r i t h m1 2 2 2 l a g u e r r e 迭代的性质 对于不可约对称三对角矩阵( 如( 2 1 1 1 ) 式) ，其特征多项式 f ( 2 ) zd e t ( t 。一五l ) ， ( 2 2 2 1 ) 的所有根为互不相等的实根 2 7 】设厂( 五) 的根为 五 以 并令五= 一o 。和乃。= 佃对x ( 磊，矗+ 。) ，上掣钟8 迭代t ( 功和t ( 砷分别定义如下 “功爿+ 磊啄菰每需( 2 2 2 2 ) 上一z2工+：弓ii；_：ii二iiii：iniii；i需。 k z , z z j ) c 一锷，一肛小”舒2 叫笔。 并且有下面的命题成立 命题2 1 2 7 。4 对x ( 厶，气+ 1 ) ，州= o ，i ，n ，则有下述结果成立： ( 1 ) 乙 t ( z ) x t ( x ) 以+ l ； 大连理工大学硕士学位论文 ( 2 ) 存在常数t 和c _ 使得 当x 趋于九+ ，时，有阪( x ) 一九+ 。l 0 l x - l 。卜 当z 趋于厶时，有| t ( 石) 一丸i c p 一丸r 根据命题2 1 ，如果令 ，j o 拶j _ t ( x ) = ( t ( t ( 工) ) ) ，! 【一 世1 = 芷( x ) = t ( t ( t ( x ) ) ) ， 则有 厶卜世 趔 互 毫1 1 2 斗九+ 1 ( 2 ，2 ，2 4 ) ( 2 2 2 4 ) 式说明当x ( 丸，丸+ ，) 时，l ( x ) 从z 的左侧趋向于以，t ( x ) 从石的右侧趋 向于厶，文【2 7 】中指出只要x 非常接近九( 或非常接近九+ 。) ，一般调用( 2 2 2 3 ) 式 ( 或( 2 2 2 2 ) 式) 两、三次便可收敛于丸( 或丸。) 2 2 3 病态接近特征值及改进的l a g u e r r e 迭代 根据文献 2 7 】，当多项式有重根时三q 鲥e 迭代只是线性收敛虽然不可约对称三 对角矩阵的特征值为两两互不相同的实数，但是，它们中的某些特征值也可能非常接近 以致于在数值上无法区分，例如w i l k i n s o n 矩阵【2 7 1 当矩阵有一个由r 个非常接近的特征 值构成的束时，在计算中被看作为一个r 重特征值，在这种情形下，l a g u e r r e 迭代的收 敛率仅为线性的 为了在上述情形下加速l a g u e r r e 迭代收敛，对l a g u e r r e 迭代( 2 2 2 2 ) 式和( 2 2 2 3 ) 式作重新定义 2 7 , 1 3 1 ：对正整数，( 0 ，h ) ，定义l a g u e r r e 迭代+ ( x ) 和上，一( z ) 如下 k 。卜h 焉硪意蘸两两 2 3 。 k 卜h 焉孺弱n 孺零两 q 2 3 2 我们称( 2 2 3 1 ) 式和( 2 2 3 2 ) 式为改进的l a g u e r r e 迭代我们需要注意的是：当，= 1 时，( 2 2 3 1 ) 式和( 2 2 3 2 ) 式就分别是( 2 2 2 2 ) 式和( 2 2 2 3 ) 式因此，我们在 调用l a g u e r r e 迭代时，通常是先令，= 1 ，然后计算，来调用( 2 2 3 1 ) 式或( 2 2 ，3 2 ) 式 王文超：求反对称三对角矩阵特征值问题的分治算法研究 由上面的两个式子，可以得到两个序列 x ，( + k = 茸+ ( z ) ；+ ( + ( + ( 功) ) ， ，j i 一 砖? = 茸一( 工) = 一( t 一( 一( 功) ) 当工为，( 五) 的，重根时，如果z z 并且f ( a ) 在区间( x ，石) 内无根，则序列 x ，( + k ，j = 1 ，2 ，以三次收敛速度递增地收敛于x ；类似地有，如果z z 并j | - f ( 2 ) 在区 间( x + ，曲内无根，则序列! ，k = l ，2 ，9 - - 次收敛速度递减地收敛于x 数值试验表 明，即使在病态根束的情形下， ( 2 2 3 1 ) 式和( 2 2 3 2 ) 式也能达到三次收敛 对上述改进的l a g u e r r e 迭代+ ( x ) 和t 一( x ) ，有下面定理邮】： 定理2 2 设 五 以为多项式，( a ) = d e t ( t - j t y o ) 的根，则对x ( 厶，忍+ 。) 和任 一正整数， o ，则厶 x t ( z ) 厶+ ( x ) t + ( x ) 丸+ ，_ ( 2 ) 如果一婴 o ，则厶+ 。 一( 。) 厶一o ) t ( 砷 x 乃1 ， ji x ) 为了方便起见，记凡= 一o 。和五+ 。= + m 注2 2 ：设x ( 以，厶+ 。) 为l a g u e r r e 迭代的初始点并用它逼近特征值九+ l ，在某些青形下， 在离矗+ ，右边非常近的地方可能存在一组特征值，即五。 厶。 屯。，它们相对于 x 离厶+ 。的距离，非常靠近丸。如图2 1 所示 在这种情形下，在达到九+ 。的能够具有三次收敛速度的邻域前，实际计算需要花费 。1 。“。一九 九。x九m + i九r g l + f 图2 1 九+ 1 为病态特征值 f i g u r e2 1 乙+ jb ep a t h o l o g i c a l l yc l o s ee i g e n v a l u e 许多次调用l a g u e r r e 迭代这时，我们可以选择其它初始点来重新开始l a g u e r r e 迭代或 者调用改进的昭e r r b 迭代+ ( x ) 来加速收敛过程 大连理工大学硕士学位论文 2 2 4l a g u e r r e 迭代的调用及停止标准 i l a g u e r r e 迭代的调用 设x ( 丑， + ，) ，如果我们调用工醒8 p 迭代t ( 功来逼近 ，只有当拶进入五的 某个邻域，回时，单调递减序列世“，拶“，才会三次收敛于五简单的观察( 图 2 2 ) 表明， 图2 2 当位于区i 司( ， + 1 ) 的x 接近于 时一- ，( x ) ，( 石) 总是小十零的 f i g u r e 2 2 一厂( x ) ， ) a l w a y sa s s u m e sa n e g a t i v e v a l u e w h e nxi s n e a r 丑i n ( 五，以+ 1 ) 当我们调用l a g u e r r e 迭代丘( x ) 来逼近五。时，有类似的结果：当位于区间( 丑， + 一) 的z 接近矾1 时_ 锵总是大于零的- 综上所述，我们给出调用如下三昭“p h 它迭代标准：设xe ( 五， + ，) 当我们用昭 p r 阳迭代l ( x ) 来逼近 ，如果在x = 拶处有一号等 o ，则调 用t ( 工) ，即：茸( z ) ：e 西：商石) ) ) ，k = 1 ，2 ，；否则，选择其它点为迭代初始 王文超：求反对称- s 2 寸_ 角矩阵特征值问题的分治算法研究 点，如用兰二争丛代替# 1 ，其中假设已经知道 + 。为区间( z ，互+ ，) 中的唯一特征值，并对 新的迭代初始点进行前面的检验并递归到满足条件为止 如注2 所提到的，如果瓦的某些特征值非常接近 ，即在五左侧有病态特征束 4 一。 一。 ，与初始点工= 拶相比较，它们离五较近这时即使条件一婴 o k x ) 在x = 世处能够满足，也需要很多步迭代才能收敛于丑为了避免这个问题，我们可以 选择其它迭代初始点或者用五一( z ) 代替t ( x ) 进行加速另一种情形可以作类似的调整 i i l a g u e r r e 迭代的停止标准 根据文 1 3 ，基于改进三项递归式( 2 2 1 6 ) 一( 2 2 1 8 ) 式l a g u e r r e 迭代计算的 特征值精度为 i ( 丑) 一乃i 2 5 s 。m a x ( 膨+ i 局+ t i ) + l 五i 占， ( 2 2 ，4 ，1 ) 其中( 。) 表示浮点运算，。代表+ ，一，x 和+ 四贝运算，且满足 ( x 。y ) = ( x o y ) o + 占) ， 占为机器精度 本文中l a g u e r r e 迭代的停止标准采用文 2 0 】的停止标准 j x 仕“一x 肚j 2 z 伸一x 阽。1 - 1 x 耻“一z 忙1 ) r ， ( 2 2 4 2 ) 这里f 为容许误差这个标准是基于下述考虑的 令 j 石( “) 一篁( t ) j q * 2 l 币l 此时，当k 寸o o ， z ”弱收敛于五，靠 l 且敛单调递减，即】 o ，= 1 , 2 ，玎一1 若 存在，= o ，l - j