（计算数学专业论文）迭代反位移变换的arnoldi算法的一种变形.pdf

迭代反位移的a r n o l d i 方法的一种变形 摘要 近年来，直接投影法成为求解大规模二次特征值问题的一种常用方法该方法 将大规模二次特征值投影到适当选取的低维子空问，从而达到降阶和保持原问题结 构的目的迭代反位移的a r n o l d i 方法足一种新的直接投影法，它结合了反位移变 换，并通过正交投影，利于r a y l e i g h r i t z 过程产生的r i t z 值和r i t z 向量分别作为 原二次特征值问题的近似特征值和近似特征向量然而进步的理论分析表明该算 法具有收敛性态的不规则性：r i t z 值收敛，但r i t z 向量可能收敛得非常慢甚至发 散。为克服这种内在隐患，基于残量范数极小原则，本文提出了利用精化向量来实现 迭代反位移的a r n o l d i 方法的变形的构想，给出了新算法的实现方式，并在理论与实 际算法上体现本文所做的修改对于原来算法的改进作用 本文分一下四个部分：第章主要介绍相关的问题背景，解决这类问题的基本方 法以及与论文相关的研究方向及发展动态；第二章简要的描述了迭代反位移的a r n o l d i 方法；第三章在引入精化向量的基础上，具体给出了迭代反位移的a r n o l d i 方法的变 形的主要思想，并在理论上证明了该算法的优越性；最后一章是数值试验，对于不同 类型的问题进行测试，体现了改进后算法的有效性 关键词：a r n o l d i 过程，r i t z 值，r i t z 向量，迭代反位移的a r n o l d i 方法，精化向 量，精化近似特征向量 迭代反位移的a r n o l d i 方法的一种变形 a b s tr a c t l v t h ed i r e c tp r o j e c t i o nm e t h o di st h em o s tp o p u l a rm e t h o df o rs o l v i n gl a r g es c a l e q u a d r a t i ce i g e n v a l u ep r o b l e m s ( q e p ) t h i sk i n do ft h em e t h o dp r o j e c t st h el a r g eq e p t oaw e l lc h o s e nl o w d i m e n s i o ns u b s p a c ei no r d e rt op r e s e r v et h es t r u c t u r eo ft h eo r i g - i n a lq e p t h ei t e r a t e ds h i f t a n d - i n v e r ta r n o l d ia l g o r i t h mi san e wp r o j e c t i o nm e t h o d ， i tc o m b i n e sw i t has h i f t - a n d - i n v e r tt r a n s f o r m a t i o na n de m p l o y st h er a y l e i g h - - r i t zp r o - c e d u r et of o r mt h ea p p r o x i m a t ee i g e n p a i r s b u tt h ef u r t h e rt h e o r ya n a l y s i sp r o v et h a t t h em e t h o dm a yc o n v e r g ei r r e g u l a r l y ：t h er i t zv e c t o r sa r em o r ed i f f i c u l tt oc o n v e r g e c o m p a r e dw i t ht h ec o r r e s p o n d i n gr i t zv a l u e s b a s e do nar e s i d u a ln o r mm i n i m i z i n g s c h e m e av a r i a n to ft h ei t e r a t e ds h i f t a n d i n v e r ta r n o l d im e t h o di sp r o p o s e di nt h i s p a p e r t h er i t zv a l u ei su s e d 嬲t h ea p p r o x i m a t ee i g e n v a l u e w h i l et h ea p p r o x i m a t e e i g e n v e c t o ri sd e r i v e db ys a t i s f y i n gac e r t a i no p t i m a lp r o p e r t y , a n di tc a nb ec o m - p u t e db yas m a l ls i z e ds i n g u l a rv a l u ep r o b l e m c o m p a r i s o n sa r ed o n eb yn u m e r i c a l e x p e r i m e n t sa n di ts h o w st h a tt h en e wm e t h o dh a sb e t t e rp e r f o r m a n c e t h i sp a p e ri n c l u d e sf o u rp a r t s i nt h ef i r s tp a r t ，r e l a t e dp r o b l e m sa n db a c k g r o u n d i si n t r o d u c e da n dt h eb a s i cm e t h o dt os o l v et h el a r g es c a l eq e pi sa l s oi n c l u d e d i n t h es e c o n dp a r t ，w eb r i e f l yd e s c r i b et h ei t e r a t e ds h i f t - a n d i n v e r ta r n o l d im e t h o d i n t h et h i r dp a r t ，b a s e do i ld e s c r i b i n gt h er e f i n e dv e c t o r ，t h em a i ni d e ao fav a r i a n to ft h e i t e r a t e ds h i k a n d i n v e r ta r n o l d im e t h o di sp r e s e n t e d t h el a s tp a r ti st h en u m e r i c a l e x p e r i m e n t w et e s taf e wp r o b l e m sa n dt h er e s u l t ss h o wt h a tt h en e wv a r i a n th a s b e t t e rp e r f o r m a n c et h a nt h eo r i g i n a la l g o r i t h m k e yw o r d ：a r n o l d ip r o c e s s ，r i t zv a l u e ，r i t zv e c t o r ，i t e r a t e ds h i f t a n d i n v e r ta r n o l d i a l g o r i t h m ，r e f i n e dv e c t o r ，r e f i n d ea p p r o x i m a t ee i g e n v e c t o r 厦门大学学位论文原创性声明 兹呈交的学位论文，是本人在导师指导下独立完成的研 究成果。本人在论文写作中参考的其它个人或集体的研究成 果，均在文中以明确方式标明。本人依法享有和承担由此论 文而产生的权利和责任。 责任人( 签名) ：丢刍为 加睁r 月；j 日 厦门大学学位论文著作权使用声明 本人完全了解厦门大学有关保留、使用学位论文的规定。 厦门大学有权保留并向国家主管部门或其指定机构送交论文 的纸质版和电子版，有权将学位论文用于非赢利目的的少量 复制并允许论文进入学校图书馆被查阅，有权将学位论文的 内容编入有关数据库进行检索，有权将学位论文的标题和摘 要汇编出版。保密的学位论文在解密后适用本规定。 本学位论文属于 2 、不保密( 、) 篡隔端嬲 第一章引言 第一章引言 1 1 问题的来源 在大量的应用科学和工程计算中，如流体力学，电路模拟与设计，结 构力学，计算机分析，动态系统网络和声学的动态分析等科学工程领域， 许多问题都归结为大型的二次特征值问题【1 8 】 l ( a ) x ：= ( a a 2 + b a + c ) z = 0 ，( 1 1 ) 的数值求解，其中a ，b ，c 冗n 煳。若a 可逆，则( 1 1 ) 数学等价于求解 l l ( 入) z ：= ( j 入2 + a 一1 b a + a 一1 c ) z = 0 ，( 1 2 ) 对于在实际问题离散后得到一些阶数为几千阶有时为几万阶甚至几十 万阶的矩阵需要的往往是二次特征值问题的依实部最大或最小，依虚部最 大或最小，依模最大或最小的若干个特征值或某区域中的特征值及其对应 的特征向量 正是由于二次特征值问题在许多学科中的广泛应用，因此二次特征值 问题求解的理论研究，算法的开发和软件的研制等是当今计算数学和科学 与工程计算领域的重大课题，寻求和研究大型二次特征值问题部分特征对 的数值方法成为国际学术界的研究热点之一 1 2 求解大规模二次特征值问题的基本方法 线性化方法是求解大规模二次特征值f 可题的最基本方法之一它是将 ( 1 1 ) 式化为“线性”特征值问题 ( _ 芦) ( ：) = a ( 言；) ( ：) c - 3 ， 则( 1 1 ) 和( 1 3 ) 有相同的特征值，还有其他一些缵眭化形式【2 3 ，2 4 1 ，那么 原本求解( 1 1 ) 的大规模的二次特征值i ；7 题就转化为求解( 1 3 ) 的一般 特征值问题。经过多年的研究，对于( 1 3 ) ，已有成功开发的k r y l o v 子空 第一章引言 2 间【2 】，位移求逆的a r n o l d i 2 3 ，2 5 ，有理k r y l o v 2 6 ，j a c o b i d a v i s o n 【2 7 】等标 准化方法可以使用但线性化方法也存在着这样的弊端：由于线性化，问 题的规模变为了原来的两倍，这使得存储量和计算量显著增加；另外出现 了线l 生化后失去了原问题的矩阵结构特点的现象，因而不能有效的发挥原 问题的矩阵结构特征。 近年来，人们在求解大规模二次特征值问题时往往倾向于运用直接投 影法，其主要的原因就在于该类方法能充分保持原问题矩阵结构特性，它 的主要思想是将问题( 1 1 ) 向适当选取的低维子空间进行特定的投影，将 原大规模问题约化为中小规模的二次特征值问题，用标准化方法求解出新 的中小规模二次特征值问题的特征对，用其作为原大规模问题的部分特征 对的近似。在众多的直接投影法中，比较受欢迎的有：残量迭代法，j a c o b i - d a v i s o n 方法，s o a r 方法等【1 3 ，1 6 ，1 7 ，1 8 ，1 9 ，2 0 ，2 l ，2 8 】 文献1 也提出了一种利用a b 产生的k r y l o v 子空间的直接投影 法一一迭代反位移的a r n o l d i 方法 【1 】中指出，如果直接将原问题正交 投影到a _ b 产生的k r y l o v 子空间，这样产生的投影子空间就不包含矩 阵c 的信息，而【1 】中的残量分析也表明由此得出的近似解的残量范数与 矩阵c 之间存在着密切的联系，因此此时可能无法在投影子空间上得到 较好的近似解。为减小矩阵c 对残量范数的影响，该算法与反位移变换 相结合，将a 以b 产生的k r y l o v 子空间作为投影子空间，并利用正交投 影过程产生的r i t z 值与r i t z 向量作为原问题的近似解，同时【1 】在理论 和实际算法中验证了该算法在消除矩阵c 对残量范数影响力方面的有效 性然而对上述算法进行适当变形，我们从中发现其实它是一种特殊的线 性化后并进行正交投影的方法，贾仲孝在【4 ，51 中指出在一般的特征值问 题中，对于非对称矩阵a ，当我们用正交投影方法求解它的部分特征对 时，其残量范数的收敛形态可能不规则，即近似特征值收敛而近似特征向 量收敛的非常慢甚至不收敛，而这也成为本文中新算法的出发点本文提 出的新算法与迭代反位移的a r n o l d i 方法相比，新算法保留了r i t z 值作为 特征值的近似，而在特征向量选取方面，用在投影子空间满足最优条件的 精化近似特征向量作为所求特征向量的近似。本文中的理论分析和数值结 第一章引言 果都表明这种新的方法更有效。 3 第二章迭代反位移的a r n o l d i 方法 4 第二章迭代反位移的a r n o l d i 方法 2 1a r n o l d i 方法 考虑广义的特征值问题 ( a a + b ) z = 0 ， 其中a ，b 冗n x n 且a 可逆。则( 2 1 ) 数值等价于求解 ( ，a + a b ) x = 0 ( 2 1 ) ( 2 2 ) 给定k 维k r y l o v 子空间k ( 4 b ，9 1 ) ，由a r n o l d i 过程产生它的一组标 准正交基q 1 ，吼，其矩阵表达式为 ( a b ) q 沪q h k + h k + l , k q k + l e 虿， 其中q k = ( 口1 7 一，钒) ，风= q t ( a - 1 b ) q k 是上h e s s e n b e r g 阵 对( 2 2 ) 利用正交投影法， u入+za-1【b。)x。a_kic。bk(，a吼-，1 b q 1 ) 于是我们就有了正交投影后的中小规模特征值问题： ( ，a + 风) z = 0 如果( ，乱) 是一风的特征对，那么( ，q k u ) 就被称为r i t z 值和r i t z 向量，同时也作为a - - 口的近似特征对上述求解一般特征值问题的方法 被称为a r n o l d i 方法近似特征对( 0 0 ，q 七乱) 的残量住为 氟= ( 0 0 i + a 。b ) q = o o q 七u + q 仇乱+ h k + l , k q k + 1 e t u = o o q k u o o q k u + h k + 1 k 吼+ l u k h k + l ，k q k + l u k ， 第二章迭代反位移的a r n o l d i 方法 5 其中“= ( u l ，u ) t 外部特征值u 的最后个分量札t 与 川e 都会随着迭代步数k 的增 大，出现相应减少的情况【2 9 】，因此a r n o l d i 方法在不考虑存储量和计算 量的条件下算法最终将会收敛。 对于大规模二次特征值问题( 1 1 ) ，若a 可逆，我们同样使用上述的 a r n o l d i 方法，在产生矩阵a _ 1 b 在子空间也( a b ，q 。) 的标准正交基饥 下的限制矩阵风的同时，我们也得到矩阵a - 1 c 在子空间( a b ，g 。) 的标准正交基q t 下的限制矩阵g ： g k = q t ( a c ) q k 对( 1 2 ) 进行正交投影， u a 2 + a b z a + 瓦a 茫。- a i c ，) b x ，上级i c ，k ( a 。b 9 1 ) 于是形成投影后的中小规模的特征值问题： l k ( a ) x = ( ，a 2 + h a + g k ) x = 0 ( 2 3 ) 如果( o o ，u ) 是( 2 3 ) 的特征对，那么( o o ，饥t 1 ) 仍被称为r i t z 值和r i t z 向量，并作为( 1 1 ) 的近似特征对近似特征对( 如，仇u ) 的残量【1 】为 r k = 如h k + l ，老u 七4 魄+ l + a k u ， 其中a k = c q k a q k g k ，u = ( u l ，u ) t 同a r n o l d i 方法所做的分析一样，随着k 的增大，外部特征值1 1 的最 后个分量u t 与h m ，七也会出现相应减少的情况【2 9 】，因此残量的前半部 分o o h k + l , k 1 1 k a q 知+ 。将随着k 的增大而减小，但另一方面对于残量后半部分 而言，除了矩阵c 是小矩阵的少数情况外，t 都无法被忽略，所以当残 量前半部分足够小的时候，我们有 n i i l l 七u l i l 第二章迭代反位移的a r n o l d i 方法 6 因此，如果直接应用a r n o l d i 方法于大规模二次特征值问题，则会出 现残量范数停滞在i l t u i | 水平上的现象，针对这一缺陷，叶强 1 】在利用 a r n o l d i 方法的基础上结合了反位移变换的思想，由此得出了迭代反位移 的a r n o d l i 方法。 2 2 迭代反位移的a r n o l d i 方法 在实际应用中，我们往往需要的是某区域中的部分特征对假定欲求 出目标位盯附近的特征值，设入。= 盯为初始值，进行反位移变换，6 ：= l , x = 1 a a o ，则l ( a ) 可变为 三( 入)= a ( a + a o ) 2 + b ( 入+ 入o ) + c = a a 7 2 + ( 2 a o a + b ) 入+ l ( a o ) ， 则 l ( 0 ) = 臼2 l ( a ) = 0 2 a + 口台+ d ， 其中a = l ( x o ) ，宫= 2 a o a + b ，e = a 如果d e t ( f i ) = 0 ，那么显然知是( 1 1 ) 的特征值反之，若a 是非奇 异阵，则l ( o ) 可转化为二次特征值问题 l ( o ) z = a l ( o ) x = ( p 2 ，+ o a 一1 雪+ a 一1 d ) z = 0 ，( 2 4 ) 易知( a ，妒) 是( 1 1 ) 的特征对数值等价于( 8 ，垆) 是( 2 4 ) 的特征对，其 中0 = i a 一知因此，计算o r 附近的特征值数值等价于计算( 2 4 ) 的模最 大特征值。 关于a 一- 雪与初始值v 。的a r n o l d i 过程的矩阵表达形式为： a 一1 雪= 日。+ h m + l , m v m + 1 e 。t ，( 2 5 ) 其中+ 1 = ( v m ，u r n + 1 ) = 【u l ，t ，m + l 】是n ( m + 1 ) 阶矩阵，它 的列构成了j i c m + 1 ( a 一1 台，v 1 ) 的一组正交基，是m m 上h e s s e n b e r g e 阵 第二章迭代反位移的a r n o l d i 方法7 对( 2 4 ) 利用正交投影法，即用满足 p 枷分篙黠幢川】) 【 p 瓦m ( a _ 1 雪，u 1 ) 、 的( 百，驴) 称为( 2 4 ) 的近似特征对 因为j i c m ( a 一1 台，u 1 ) ，所以 驴= 9 ，( 2 6 ) 由关系式( 2 5 ) ? ( 2 6 ) ，( 1 ) 可转化为下面的小规模二次特征值问题 ( 俨，+ 莎，mq - g 。) 9 = 0 ，( 2 7 ) 如果( 百，g ) 是( 2 7 ) 的模最大特征值所对应的特征对，那么( 百9 ) 就 是( 2 4 ) 的模最大特征值所对应的特征对。令天= 盯+ i 0 ，则迭代反位移 a r n o l d i 方法用( 天，乒) 来作为( 1 1 ) 的近似特征对 定义近似特征对( 天，p ) 的相应残量为： r ：= ( a a 2 + b a + c ) 9 ， 则该算法的残量范数的上界由下面定理给出 定理1 ：如果每是。( p ) 模最大的特征值，矿= 队y ，y 2 ，y m 】? 为相应 的特征向量，令9 = 雪，及天= a o + 1 否，则有 懈埘删刚料h , , , + l , m y m l l ( a o ) q m + l i i + 错 其中。= b 一l ( a o ) v m o m ，0 m = 昭五一1 b 证明：见” 随着a r n o l d i 算法迭代步数k 的增大，残量范数界的前半部分将逐渐 减少 2 9 】，后半部分虽然含有| f m 训l ，但它与i a ，一a o l 2 同阶，随着位移 a 。的不断更新，a ，与知接近的程度将不断增大，残量后半部分的影响 力也将逐步减弱。在算法的具体实现方式上，将重复a r n o l d i 过程的迭代 直到残量前半部分在某种程度小于后半部分( 也就是说更多的a r n o l d i 过 第二章迭代反位移的a r n o l d i 方法 8 程迭代对残量范数的影响不大了) ，然后利用满足i a 一口i 叩的新生成的 特征值a 作为新的位移进行重新启动 算法1 迭代反位移a r n o l d i 方法 1 输入：初始位移盯，位移更新准则7 7 ，初始单位向量v 。，精度t o l ， a r n o l d i 过程最大迭代步数k ，令a o = 口 2 计算反位移变换后的二次特征值问题的新矩阵a ，宜，0 3 迭代：f o rj = l ：k ( 3 1 ) 利用a r n o l d i 过程确定h e s s e n b e r g 阵马以及b ( _ - 1 雪， 1 ) 子空间的一组标准正交基巧，并计算g 3 = q ( a 一1 c ) q j ； ( 3 2 ) 计算托2 + 屿p + g j 的模最大特征值对应的特征对( p o , u ) ； ( 3 3 ) 计算原问题的近似特征对入，= 知+ 厮1 ； ( 3 4 ) 计算残量范数界前半部分7 ，= 乌等h 。+ 1 ，m i a l ( x 。) 口m + 。i i ( 3 5 ) 计算残量范数界后半部分7 2 = 与样ii ( b v , n l ( a o ) o 。) y l l ； ( 3 6 ) 满足，y 1 o 1 7 2 ，跳出循环； 4 计算近似特征向量：z = v i a ； 5 位移更新：若l 入l 一仃i 7 7 ，则a o = a l ； 6 计算残量r = ( 增a + a 1 b + c ) z ，如果阡丽矸煳岛司而 t o l 则停 机，否则转步骤2 上述算法的具体实现细节可参考文献 1 】然而对于算法1 ，我们进 行适当的变形，发现了该算法的个内在的性质。 性质l ：在精确计算的前提下，迭代反位移的a r n o l d i 方法等价于线性 化后并利用印。n t ( 苫曼) ，作为投影子空间的正交投影化方法 第二章迭代反位移的a r n o l d i 方法 证明：定义 耻( 一对 p = ( 一0 1 台一a 0 _ l 。) ， 及。= ( 曼) ，男s 么有 b 。= k t m p 。 容易看出( 2 7 ) 的特征对( 百，9 ) 满足 b 。( ? ) = 舀( ? ) ( 2 8 ) 所以迭代反位移的a r n o l d i 方法可看作是一种特殊的线性化后并结合 正交投影的一种投影化方法。 9 第三章迭代反位移的a r n o l d i 方法的一种变形 1 0 第三章迭代反位移的a r n o l d i 方法的一种变形 3 1 精化投影方法 对于非对称矩阵a ，当用传统的正交投影方法求解它的部分特征对 的时候，r i t z 向量收敛可能比r i t z 值收敛要困难的多，特别是矩阵为非 h e r m i t i a n 阵时，为了解决这个问题，贾【4 , 5 】提出了对于确定的近似特征 值，用在投影子空间上满足最优条件的向量( 称为精化向量) 来代替r i t z 向量，其思想如下： 对于任意的求解子空间配，以及近似特征值天( 可以是r i t z 值或其他 近似特征值) ，精化投影类方法是用满足 ( a 一入j r ) 蟊i | = m i n l i ( a a r ) uj l ， u 瓦 i i 钱= 1 的豆作为与近似特征值天相对应的近似特征向量，并称缸为a 在子空间 芘上与天相对应的精化向量。对指定的求解子空间疋，精化向量可以通 过计算一个小规模的特征值问题或小规模的奇异值分解得到 文献【4 , 6 】中，对，对般的投影空间瓦，建立了精化向量豇的先验 误差界和其它一些重要结果目前，已开发成功的精化型算法有隐式重新 开始的精化a r n o l d i 方法【7 】 隐式重新开始的精化调和a r n o l d i 方法【8 1 精 化子空间迭代法【9 】和精化的位移求逆的a r n o l d i 方法【1 0 】等 3 2 迭代反位移的a r n o l d i 方法的一种变形法 由第二章的性质1 及传统的正交投影法的收敛不规则形态，我们可以 在理论上得出r i t z 向量的不收敛l 生可能出现在迭代反位移的a r n o l d i 算法 中 为了改善r i t z 向量的收敛性能，此时根据精化投影方法的原则，应该 在子空间厄。( a 。台，u 。) 当中找个单位长度的向量u 满足下面的最优条 第三章迭代反位移的a r n o l d i 方法的一种变形 件： l i ( 否2 十每a 一1 台+ a 一1 0 ) u 1 1 = m i n i l ( 莎2 ，+ 舀a 一1 唐+ a 0 ) y t l ， y 瓦仇( a 一1 台，u 1 ) l l y l l = 1 ( 3 1 ) 此条件等价于找一个三使得l , = 三，也就是下式成立 i l ( 驴+ 百a 一1 台+ a 一1 e ) 引i = m i n l ( 驴+ 百a 一1 台+ a 一1 d ) z ， z c m i = 1 ( 3 2 ) 对于( 3 2 ) ，可通过小规模矩阵萨+ 否a 一宜+ a 一- e 的奇异值 分解来得到，它的最小奇异值所对应的右奇异向量三即满足最优条件，用 所得的向量三求出u = 互，那么( 天，u ) 其中天= 口+ l 百可作为( 1 1 ) 的 特征对的近似，那么u 是关于莎在c 。( a - 1 雪，u 。) 上最佳的近似向量，也 被称为精化近似向量【2 2 】，那么它们所对应的的残量范数为： 恫i = i i ( i 2 a + 支b + c ) 钍l l ， 日 l l ( 天2 a + 5 , 3 + c ) 引l a 咕+ 耕+ b ( 吾+ 糊+ c 】划 ( a 嘉+ 2 a p o a + a ：a + 万b + b a o + c ) 三i i a 嘉+ 万1 ( 2 沁a + b ) + 增a + b a o + q 三 l i a ( 天一盯) 2 + 亩( 天一盯) - fc 】引i = l i a ( 5 , 一仃) 2 ( 百2 j + 舀a 一1 台+ a 一1 0 ) 引l 1 5 , 一口1 2 l | a i ( 百2 ，+ 百a 一1 台+ a 一1 e ) 引l( 3 3 ) l 天一盯1 2 j i a l j t n ( 舀2 - 4 - 舀a 1 啻+ a 一1 d ) 从上述讨论中我们得到了新的算法。 第三章迭代反位移的a r n o l d i 方法的一种变形1 2 算法2 迭代反位移a r n o l d i 方法的一种变形法 1 输入：初始位移口，位移更新准则7 7 ，初始单位向量u 。，精度t o l ， a r n o l d i 过程最大迭代步数k ，令a 。= 盯 2 计算反位移变换后二次特征值问题的新矩阵a ，台，e 3 迭代：f o rj = l ：k ( 3 1 ) 计算7 = a 一1 b v j ； ( 3 2 ) 正交化：r = 7 一v j h j ，其中h j = 吁7 ； ( 3 3 ) v 3 + l = r h j + l j ，其中+ l j = l | r i i ； ( 3 4 ) 计算q = 曙a 一1 e k ； ( 3 5 ) 计算口2 j4 - 口+ q 的模最大特征值痧； ( 3 6 ) 近似特征值天= a o - f 百一1 ； ( 3 7 ) 由关系式( 3 1 ) 及( 3 2 ) 计算近似特征向量； ( 3 8 ) 计算残量r = ( 增a + a b + c ) x ，如果残量范数满足停机准则， 跳出循环； 4 位移更新：若h d r i 叩，则a o = a 1 ； 5 重新启动：_ z ，转步骤2 注：( 1 ) 矩阵与向量的乘积运算y = a _ 1 雪z 可分两步进行，首先计算童= 台z ，接着利用y = u 1 ( l 。仝) ，其中l ，u 由矩阵a 的l u 分解求出 ( 2 ) 求解( 3 2 ) ，我们把算法中( 3 1 ) ，( 3 4 ) 每一步计算的a - 1 雪，a 。d k 存储起来，这样矩阵k = 舀z v m4 - 私一 h v m + a 。d 就不需要额外的矩阵 乘向量的运算量，接着( 3 2 ) 可通过n m 阶矩阵k 的奇异值分解来求解 出数值解，这需要o ( n m 2 ) 的计算量。虽然这与迭代反位移a r n o l d i 方法 相比，增加了奇异值分解这部分的相应计算量，但却加快了收敛速度，关 于这情况可以见数值例子 第三章迭代反位移的a r n o l d i 方法的一种变形 ( 3 ) 下面的数值试验中，我们把停机准则设定为邴丽丽揣岛_ 可而 3 3 算法的收敛性分析 在以下的算法收敛性分析中，我们设定( p ：妒) 为( 1 1 ) 的精确解，( 莎，u ) 为算法2 求解出的近似解下面的结论揭示出在投影子空间所包含的信息 足够完整的情况下，即妒与的夹角正弦值s i n g ( 妒，) 接近0 时，算法2 的良好的收敛性态。令1 m 。为仇列向量张成的空间，丌m 为j i c m ( a 一- 雪，u 。) 上的正交投影算子，砚。为k 上的正交投影算子。 定理2 ：假设e = i i ( z 丌m ) p i i ，则 i i ( 萨，+ 百a 一，房+ a 一- o ) 让i i 堕二二旦1 l 二盟止坌二垒l 之翼鱼圣 上业兰二堡旦上迹 、l 一 证明：令雪= 号尚为在瓦。( a 一1 豆，杉z ) 上的单位投影向量且 e = ( j 一7 r m ) 妒， l ( 8 一) = 萨j + 舀a 一1 雪+ a 一1 e ， 我们有 嘲= 辫 l ( 曰) ( 妒一e ) 一弋厅乏f ( 占2 一目2 ) 妒+ ( 百一e ) a 一1 房妒一( 百) e 圻压i 则 i i l ( o ) # l l 生刿絮譬业继， 因为 l i l ( # ) u l i = m i n l i l ( 0 ) ：i i i ， 圣厄m ( a _ b ，u 1 ) 忪| | = 1 1 3 第三章迭代反位移的a r n o l d i 方法的一种变形 所以 幡渤训堕刿幽螋等譬盟逊世业 定舢定义川咱。“孑n 则有一。的充分必要条件 是，y 一0 证明：由于丌m = 瑶，则有 砘。= ( 0 曼) ( 罾o ) 丌2 m = ii l 。i o 吆 = ( 警意) 7 = i i ( ，一嵋j 一曼v ：) ( 苫) l i i i ( 苫二焉斗 = 而眦一螵) 妒 = 以f 订 于是，我们得到- ，0 的充要条件是，y 0 定理4 ：定义= i i a i i i i p l i ，若一0 ，则近似特征对( 天，u ) 满足 咿忪小叫咝必辱磬乒登型，掣 其中 m 钟枷+ 1 1 4 一跳= i i p i i + 伊剐+ 唏， 1 4 第三章迭代反位移的a r n o l d i 方法的一种变形 如果百是单特征值，则 忪忪小刮喾删伺， ( 3 5 ) i i f | | 引天一盯l ! 去三，二r + o ( 2 ) ， ( ) v1 一巴。 其中s 是舀的谱条件数。 证明：矩阵p 2 。的特征值百满足 百l l i p 2 。| l = i l k 二尸k m i i i i 尸i i ， 根据【4 】，我们有 i 百一e l 4 1 1 p 1 ( 2 + 7 、研) 1 1 2 ”( 1 、研) 1 2 m 令 由定理2 和( 3 3 ) 得到 i t 乒l t l 爻一a l i t a i l l l ( 萨r + 否a 一1 雪+ a 一1 c ) u l l i i - o i i i a l l 坐幽螋世耻掣型塑业堡 枷+ l | 一b 1 1 ) 驯邮i i + 帮矧 利h i 咝必譬磬乒弛盟 其中 m - - 4 枷+ i i 枷1 ) ，- | | p l i + i i a 。刖+ 帮， 如果舀是单特征值，则根据【5 】中的定理3 7 ，我们有 i 百- e l _ s 俐南知( 书， 因此可以得到下面的不等式： 1 5 第三章迭代反位移的a r n o l d i 方法的一种变形 i f i i l 天一盯l ! ! 喾+ 。( y z ) i f i | l 天一盯l ! 去兰 + o ( y 2 ) v1 一e 定理4 中的( 3 4 ) 和( 3 5 ) 表明在投影子空间具备充足的信息的条件 下，算法2 的良好的全局收敛性态。 1 6 第四章数值实验 第四章数值实验 这一部分给出了m a t l a b 7 0 里得到的数值实验结果。所有的数值测试 是在一个操作系统为w i n d o w s x p ，c p u 为i n t e l ( r ) c e l e r o n ( r ) 2 6 6 g h z ，主存 为2 5 6 hr i b ，机器精度为e p s = 2 2 2 1 0 6 的机器上进行的。我们将算法 1 同算法2 进行比较来看改善效果。 所有的例子中，o r 为初始选择的位移，停机准则设定为l e - 8 ，位移更 新准则7 7 设定为o 1 ，我们用“r v s i a m ”代表算法2 ， “v s i a m ” 代表算法1 ，在下面的图像中，虚线代表的是“i s i a m ”，实线代表的 是“r v s i a m ”，横坐标代表算法的迭代步数，纵坐标代表算法的残量范 数。 例1 我们考虑阶数为1 1 8 2 的二次特征值问题，其中a = j r ，b = 2 0 0 木t r i d i a g ( 一1 ，2 0 ，一1 ) ，c 是来自于m a t r i xm a r k e tc o l l e c t i o n 3 1 中的矩阵 c a v i t y 0 5 ，我们分别选取两个初始位移盯= 7 ，o r = 9 ，计算该二次特征值问 题在位移口附近的特征值从图1 中我们可以看出“r v s i a m ”收敛所需 的迭代步数明显少于“v s i a m ”，在实验中我们也发现虽然“r v s i a m ” 每一步迭代所耗的时间多于“v s i a 皿”，但在给定精度l e - 8 下达到收敛 时“r v s i a 1 2 1 ”比“v s i a m ”耗时要少，例如当伊= 7 r 时，“r v s i a i l l ” 收敛所需时间为1 4 7 2 秒，“v s i a i n ”则需1 6 3 7 秒 第四章数值实验1 8 图1 c a v i t y0 5 的数值实验，t o l = l e 一8 ，上：矿= z 下：a = 9 例2 我们考虑阶数为1 1 0 7 的二次特征值问题，其中a = j ，b = 1 0 0 车t r i d i a g ( - 1 ，9 ，- 1 ) ，c 是来自于m a t r i xm a r k e tc o l l e c t i o n 3 】中的矩 阵g r e l l 0 7 ，我们同样选取了两个初始位移d r = 7 ，盯= 9 ，计算该二次特 征值问题在位移t o 附近的特征值。从图2 中我们可以看出与例1 相同， r v s i a m ”收敛所需的迭代步数仍少于“v s i a m ”，在实验中我们也发 eloc一毋丁口历 e l o lj帚：dl西皂 第四章数值实验1 9 现，在给定精度l e - 8 下达到收敛时“r v s i a m ”比“v s i a m ”耗时要 少，例如当盯= 9 时，“r v s i a m ”收敛所需时间为1 5 6 2 秒，“v s i a 肌， 则需1 7 8 7 秒。 t h en u m b e ro fi l e r a t i o n 图2 矩阵g r e1 1 0 7 的数值实验，d kl e 一8 ，上：仃= 7 ，下：c r - - 9 eloc面r口鬲e 巨oc一再，里2 参考文献 参考文献 1 q y e ，a ni t e r a t e ds h i f t - a n d i n v e r ta r n o l d ia l g o r i t h m 加7 q u a d r a t i cm a t r i x e i g e n v a l u ep r o b l e m ，a p p l e d m a t h c o m p u t ，1 7 2 ( 2 0 0 6 ) ，8 1 8 8 2 7 2 jz b a i ，j d e m m e l ，j d o n g a r r a ，a r u h e ，a n dh v a nd e rv o r s t ，t e m p l a t e s o r t h es o l u t i o no ja l g e b r a i ce i g e n v a l u ep r o b l e m s ：ap r a c t i c a lg u i d e s i a m p h i l a d e l p h i a 【3 】z b a i ，r b a r r e t ，d d a y ，j d e m m e la n dj d o n g a r r a ，t e s tm a t r i xc o l - l e c t i o nf o rn o n h e r m i t i a ne i g e n v a l u ep r o b l e m s ，t e c h n i c a lr e p o r t ，d e p a r t - m e n to fm a t h e m a t i c s ，u n i v e r s i t yo fk e n t u c k y , u s ，1 9 9 5 a v a i l a b l ea t h t t p ：m a t h n i s t g o v m a t r i x m a r k e t 【4 】z j i aa n dg w s t e w a r t ，a na n a l y s i s 巧t h er a y l e i g h r i t zm e t h o d o r a p p r o x i m a t i n ge i g e n s p a c e s ，m a t h c o m p u t ，7 0 ( 2 0 0 i ) ，6 3 7 - 6 4 7 【5 】5z j i a ，t h ec o n v e r g e n c e 巧g e n e r a l i z e dl a n c z o sm e t h o d sy o rl a r g eu n s y m m e t t i ce i g e n p r o b l e m ，s i a mj m a t r i xa n a l a p p l ，1 6 ( 1 9 9 5 ) ，8 4 3 - 8 6 2 【6 】z j i a ，r e f i n e di t e r a t i v ea l g o r i t h m sb a s e do na r n o l d i sp r o c e s s o rl a r g eu n s y m m e t r i ce i g e n p r o b l e m s ，l i n e a ra l g e b r aa p p l ，19 9 7 ，25 9 ：1 2 3 7 1z j i a ，五p o l y n o m i a lc h a r a c t e r i z a t i o n so ft h ea p p r o x i m a t ee i g e n v e c t o r s b yt h er e f i n