（计算数学专业论文）二阶双正交过程和二次特征值问题求解.pdf

复旦大学硕士学位论文 摘要 在本文中，我们提出了二阶双正交过程的算法，通过这过程可以得n - 阶 k r y l o v 子空间的一组双正交基。有别于标准的k r y l o v 子空间，二阶k r y l o v 子空 间的构造需要基于一对矩阵和初试向量。给定一组双正交基之后，我们采用斜投 影技术来求解二次特征值问题。由于在二阶双正交过程中只涉及到矩阵向量乘法 运算，同时通过揭示其中的矩阵关系，我们对存储量进行了优化；又因为在算法 中隐式的构造了线性化问题的k r y l o v 子空间，使得算法不仅在存储空间上而且在 运算时间上达到较优的效果，所以这个方法适用于大规模和稀疏问题的求解。我 们给出的数值例子也证实，通过这一算法得到的r i t z 值和对应的r i t z 向量可以快 速和准确的收敛到真实解，并且通常是那些具有较大模的特征值。 关键词二阶k r y l o v 子空间，二次特征值问题，双正交化，斜投影法 复旦大学硕士学位论文 i i a b s t r a c t i nt h i sp a p e r ，w eg i v eap r o c e d u r e ，i e ，t h es e c o n d o r d e rb i o r t h o g o n a l i z a t i o np r o c e - d u r e ，t og e n e r a t et h eb i o r t h o n o r n l a lb a s i so ft h es e c o n d - o r d e rr i g h ta n d t e f tk r y l o v s u b s p a e e sw h i c hb a s e d o na p a i ro fm a t r i c e sa n ds t a r t i n gv e c t o r s t h ea p p l i c a t i o no f t h i sp r o c e d u r ei st os o l v et h el a r g e s c a l eq u a d r a t i ce i g e n v a l u ep r o b l e m sv i ao b l i q u e p r o j e c t i o nt e c h n i q u e t h i sm e t h o dc a nt a k ef u l la d v a n t a g eo ft h es p a r s e n e s sf o r l a r g e - s c a l es y s t e ma sw e l la 8t h es u p e r i o rc o n v e r g e n c eb e h a v i o ro fk r y l o vs u b s p a c e b a s e dm e t h o db yi m p l i c i t l yl i n e a r i z a t i o n ，w h i c hm a k e st h es o l u t i o na c c e p t a b l ei n t e r m so fb o t hc o s ta n dt i m e ，o u rn u m e r i c a le x p e r i m e n t ss h o wt h a tt h er i t zv a l u e s ， o f t e nt h o s ew i t ht h el a r g e s tm a g n i t u d e ，a n dt h ec o r r e s p o n d i n gr i t zv e c t o r sa l w a y s c o n v e r g et ot h ee x a c to n ef a s ta n da c c u r a t e l y k e y w o r d sb i o r t h o g o n a i i z a t i o n ，o b l i q u ep r o j e c t i o nm e t h o d ，q u a d r a t i ce i g m w a l u e p r o b l e m s ，s e c o n d - o r d e rk r y l o vs u b s p a c e 复旦大学硕士学位论文 1 前言 二次特征值问题( q e p ) 是工程应用中一类十分重要的问题，在声学系统的动 态分析 1 ，2 、流体力学的稳定性分析 3 以及信号处理等应用中经常涉及到求解 二次系统的特征值问题。 在本文中，为了表述的方便，记n 阶二次a 一矩阵 4 ，5 】为 q 。( a ) = a 2 m + a d + k( 1 ) 其中的系数矩阵m ，d ，k 为给定的nx 礼的实方阵( 可以简单的将本文中的内容推 广到复数的情况) 。a c 是q 。( a ) 的一个特征值并且非零向量。是对应的右特征 向量，如果成立 q 。( a ) 。= ( a 2 m + a d + k ) z = 0 ( 2 ) 类似的，非零向量y 是对应的左特征向量，如果成立 y hq 。( a ) = y 日( a 2 m + a d + k ) = 0 ( 3 ) 满足以上要求的三元组( a ，z ，y ) 称为二次a 一矩阵的特征组 4 。所谓的二次特征值 问题就是要求解瓯( a ) 的特征组。 二次特征值问题与一般的矩阵特征值问题以及广义特征值问题相比，最大的 区别在于，q 。( a ) 有2 n 个特征值( 包括无穷大的情况) ，以及至多2 n 个右( 左) 特 征向量。由于这些向量都是n 维空间中的向量，所以当特征向量的个数超过礼 时，它们之间必然存在某种线性相关性 6 ，7 。另一个区别在于，对于二次特征值 问题，尚不存在一个简单的类似于( 广义) s c h u r 分解的标准型 8 ，9 。这两点区别 增加了求解这种非线性问题的难度。 求解二次特征值问题的通常方法为，首先将问题( 2 ) 和( 3 ) 线性化为等价的广 义特征值问题 ( a g h ) z = 0 和w h ( a g h ) = 0 ，( 4 ) g = ，h = - 警 复旦大学硕士学位论文 并且 。= ： ，叫h = g h ( a m + d ) hj 其他常用的线性化形式可以参考文献 4 ，7 1 。如果阶数n 不是很大，我们可以使 用q z 算法1 0 1 或者其他直接法来求解问题( 4 ) 所有的特征值和对应的特征向量。 对于大规模系统可以采用基于k r y i o v 子空间的方法，例如：a r n o l d 或者l a n c z o s 方法，来近似的求解其中一部分需要的特征值和对应的特征向量。上面的线性化 方法直接也易于实现，但是不免存在着一些问题。一方面，在作线性化时，我们 将原始的礼阶问题转化为等价的2 n 阶问题求解，当n 很大时，这样作可能会增加 存储量并带来计算上的困难。另一方面，原始问题所具有的特殊性质，尤其是隐 含的谱性质( 见表1 ) ，经过线性化后不一定会继续保留，这也会给问题的求解带 来困难同时影响结果的准确程度。对于大规模二次特征值问题的求解在最近几年 得到广泛的关注，一些迭代的方法被提出。例如，j a c o b i d a v i d s o n 方法1 1 ，1 2 1 ， 不过这个方法一次只关注于一个特征值，这样虽然可以在局部快速收敛，但全局 的收敛速度慢。文献1 3 1 中提出的降阶技巧可以用于求出其他的特征对。 矩阵性质特征值性质特征向量性质 1m ，d ，k 为实矩阵特征值为实数或者共。为a 的右特征向量，那 轭对么牙为j 的右特征向量 2m ，为实对称正定矩阵且特征值为纯虚数 d = 一d ? 3m ，k 为复共轭矩阵，m 正特征值为纯虚数或者z 为a 的右特征向量，那 定且d = 一d h毗( a ，一i ) 对的形式么为一i 的左特征向量 4 m 复共轭正定，d ，复共& ( 入) 0 轭半正定 表1 ：二次特征值问题所具有的一些谱性质 7 。 最近，文献 1 4 】中提出了一种基于r a y l e i g h r i t z 投影技术的方法来近似求解 ( 1 ) 的一部分特征值及其对应的右特征向量。其中包含的主要思想为，使用二阶 a r n o l d 过程构造二阶k r y l o v 子空问瓯。( a ，b ；饥) 的一组正交基。然后使用这组正 交基将原始问题投影到一个m 阶的二次特征值问题求解。由于通常来说投影后问 题的阶数m 并不会很大，可以使用上面提到的直接法来求解。这个方法最显著的 2 复旦大学硕士学位论文 特征是，由于投影的过程是直接应用于原始问题而并非某一个线性化的形式，因 此矩阵m ，d ，k 所具有的性质将显式的保留在投影后的系统中。这也就意味着， 原始问题中所具有的谱性质也将肯定保留在投影后的系统中。而且，在构造二阶 k r y l o v 子空间正交基的过程中隐式的构造了线性化问题的k r y l o v 子空间，所以此 方法也具有基于k r y l o v 子空间方法的快速收敛性。 由于文献 1 4 中的方法只构造子空间g 。( a ，b ；u ) 的一组正交基，并利用这组 基作投影，因此可以认为这个方法是正交投影法或者单边法。在本文中，我们将 这个方法推广到双边的情况，也就是说，需要构造二阶k r y l o v 子空间g 。( a ，b ；u ) 和g 。( a t ，b r ； ) 的一组双正交基。这个构造过程，我们称之为二阶双正交过程， 通过采用双边的g r a m s c h m i d t 过程将二阶k r y l o v 序列双正交化。除了不具有 短递推形式外，二阶双正交过程以及由此得到的矩阵关系类似于标准的非对称 l a n c z 0 8 过程f 1 5 ，1 6 1 。有了这样的一组双正交基之后，我们可以利用斜投影的技 术对原始问题作投影。通过求解投影问题的真实特征值和特征向量，可以得到原 始二次入一矩阵q 。( ) 的近似特征组，或者称为r i t z 组。关于将双边l a n c z o s 方法 应用于线性化之后的广义特征值问题( 4 ) ，来求解二次特征值问题的方法可以参考 文献 1 7 ，1 8 ，1 9 】。与正交投影的方法相比，本文中提出的斜投影算法的优势主要 体现在两个方面。其一，采用斜投影算法求得的r i t z 值可以更快的逼近真实解， 其收敛的速度比正交投影的方法快几乎一倍，关于这一现象可以见数值例子中的 图示。其二，用正交投影法无法直接逼近二次特征值问题的左特征向量，尤其是 对于非对称的问题。而斜投影的方法正好可以克服这一缺点，我们不仅可以得到 r i t z 值，而且可以同时得到对应的左右r i t z 向量。 在本文中，我们将采用标准的线性代数记号。也就是说，矩阵和向量将分别 用大写和小写字母表示，标量用小写希腊字母表示。单位阵记为，零向量或者 矩阵记为0 ，单位矩阵的第j 列记为e f 。向量和矩阵的维数取决于所处的上下文 环境，并与公式中出现的其他矩阵和向量的维数保持一致。矩阵a 在( i ，j ) 位置 的元素记为a 咖符号腿和c 分别表示全体实数和虚数集合。相应的，用符号r n 表示礼维实空间中向量全体，符号r “表示m n 实矩阵全体。矩阵的1 范 数、2 一范数和f r o b e n i u s 一范数分别用记号忆| 1 2 和”怯表示。矩阵的转置和 共轭装置用- t 和h 表示。我们用s p a n q l ，q 2 ，) 和s p a j l q 表示由向量序 3 复旦大学硕士学位论文 列 ) 和矩阵q 的列所张成的子空间，d i a g ( p iq ) 表示以矩阵只q 为对角元的块 对角矩阵，s i g n ( u 1 表示标量u 的符号。 本文是这样安排的。在2 中，我们给出二阶k r y l o v 子空间的定义以及这个 定义提出的初衷。接下来，我们给出本文的一个主要算法，也就是二阶双正交 过程。通过揭示其中的矩阵关系，我们可以对算法作迸一步的改进，以减少存 储量。在3 中，我们给出推广的二阶双正交过程，并利用这个过程给出求解二 次特征值问题的斜投影算法，简单的向后误差分析电将在这一节中涉及。在4 中，我们给出几个数值例子来验证算法的可行性，同时和二阶a r n o l d i 算法得到 的结果进行比较。我们在5 总结全文内容。 4 复旦大学硕士学位论文 2 二阶双正交过程 在本节中，我们将首先给出二阶k r y l o v 子空间的定义，以及提出这个定义 的初衷。接着我们给出本篇论文的一个主要算法，这个算法可以用于构造二阶 k r y l o v 子空间的一组双正交基。 2 1 二阶k r y l o v 子空间 求解大规模矩阵问题包括线性方程组 2 0 j 和特征值问题 2 l j 等是当今科学计 算研究中的重大课题。从上世纪八十年代开始，这方面的研究工作取得了许多重 大进展，其中尤以k r y l o v 子空问方法的影响最大。在求解这类问题时，我们通常 先要使用a r n o l d i 过程构造k r y l o v 子空间 咒m ( a ；“) = s p a n u ，a u ，a m 一1 ) 的一组f 交基，或者使用l a n c z o s 过程构造子空间瓦。( a ，u ) 和赶，。( a t ，u ) 的一组 双正交基，然后利用子空间的基来求得近似解。因为在构造基的过程中只涉及到 矩阵向量乘法运算，所以这类方法特别适用于大规模稀疏问题。除了矩阵问题 外，在模型降阶问题中，也会涉及到k r y l o v 子空间的应用【2 2 。但是对于二次问 题，例如：二次特征值问题和二次系统的模型降阶问题，标准k r y l o v 子空问由于 只涉及到一个矩阵a ，所以无法充分反映原始问题的非线性特征，也就无法利用 这些予空间来直接求解类似的二次问题。 在文献f 1 4 1 中，作者将标准的k r y l o v 子空间思想推广n - ：阶情形，定义了 二阶k r y l o v 子空间毋。( a ，b ；u ) ，这个子空间的构造需要涉及到两个给定方阵 a ，b 和初试向量u 。当矩阵a ，b 非对称时，类似于标准的左右k r y l o v 子空间 厄。( a ，u ) 和坛。( a r ， ) 的定义，我们也可以提出二阶k r y l o v 子空间岛。( a ，b ；u ) 和9 。( a t ，b t ； ) 的定义。 定义1 ：设a ，b 是n 阶实方阵，非零向量u ，v 的长度为礼。那么满足 r 0 r 1 q ( 5 ) 5 十 日2 一 当2 一 q b+ l一 0 u a a 复旦大学硕士学位论文 辛口 80= s i = a 丁3 0( 6 ) s j = a t 8 卜1 + b t s 卜2 当j 2 时 的向量序列 ，r 小，r r a - - 1 ) 和( s o ，8 1 ，8 m - l ，称为二阶k r y l o v 序列。由这 两个向量序列张成的子空间，分别记为 鲰( a ，b ；u ) = s p a n ( r 0 ，1 1 ，r m - 1 和 ( a t ，b r ； ) = s p a n s o m s ，8 r n - 1 ) 称为二阶k r y l o v 子空间。 ( 8 ) 当b = 0 时，上述定义中n - g # k r y l o v 子空间可以看作是对一般k r y l o v 子 空间的个简单推广。只不过，依照( 5 ) 和( 6 ) 生成的向量满2 = - - 阶的线性递推关 系。也就是说，向量r j 和8 j 可以通过关系 _ = 功( a ，b ) u 和5 = p a a 了，b t ) u 求得。这里的功( o ，卢) 是关于a 和的多项式，并且满足递推关系 p o ( a ，p ) = 1 p l ( a ，卢) = a p j ( a ，卢) = a p j 一1 ( o ，卢) + 卢功一2 ( n ，p )当j 2 时。 下面，我们给出之所以提出二阶k r y l o v 予空间定义的初衷。为简单起见，当 k = n 时，我们用r ( ) 表示 q 。( a ) = a 2 i a a b ， ( 9 ) 或者当南= 2 n 时，表示q 。( a ) 的线性化形式 6 m 何下 复旦大学硕士学位论文 正如文献 7 ，2 3 】中所指出的，采用投影法求解p k ( a ) 的特征值问题，需要构 造两个m ( s 七) 维子空间。和7 k 来逼近死( a ) 的左右特征向量所张成的子空 间。假设矩阵只。和q 。的列构成。和7 已。的一组双正交基，那么通过斜投影技 术可以把原始的阶系统投影到如下的m 阶系统 只。( a ) = 焉r ( a ) q 。 已知只。( a ) 的特征值和左右特征向量( p ，f ，g ) ，可以通过下面的关系 ( 0 ，z ，y ) = ( p ，q 。f ，p m g ) 得到巩( a ) 的r i t z 值和对应的r i t z 向量，也即原问题的近似特征组。同时，这个 方法确保右r i t z 向量所对应的残向量r = r ( 日) z 正交于子空间c 。，而左r i t z 向 量对应的残向量，= y 打r ( p ) 正交于子空间冗。 当r ( a ) = a g h 时，可以选择k r y l o v 子空涮尼。( 日，砬) 和。( 日t ，o ) 作为 冗。和c 。，即 7 k 。= 尼。( 日，n ) 和c 。= j c 。( 日t ，0 ) 如果记初试向量 也t = u ro 和。t = ”to ， 其中u ，口础。那么，可以很容易的推出，长度为n 的二阶k r y l o v 序列 q 和 ( s ， 与长度为2 n 的k r y l o v 序列( 伊町和 ( 日r ) 之间存在关系 日血= r g ， 和c 日t r 。= b 之一。 c - 。， 如果重新假设矩阵q 。和p 竹。的列构成二阶k r y l o v 子空间( 7 ) 和( 8 ) 的组双正交 基，那么由关系( 1 0 ) 可知，予空间冗。可以嵌入由矩阵d i a g ( q 。，q 。) 的列张成的 空间，也就是说 咒m c 日，也，s p a n ( ： ) 同样的，子空间c 。可以嵌入由矩阵d i a g ( p 。，b 丁p m ) 的列张成的空间，也就是洗 酽雌伽b 圳 7 复旦大学硕士学位论文 因此，我们可以得出结论由矩阵d i a g ( q 。，q 。) 和d i a g ( p m ，b 7 p m ) 的列所张成 的子空间比冗。和c 。蕴涵了更多的信息。那么同样的，如果我们使用矩阵 d i a g ( q 。，q ，。) 和d i a g ( p m ，b t 尸竹。) 对特征值问题a g h 作投影，得到的系统为 警之一( a g - h ) 一曼 = a ：只：b o 。 一 p 只三：a b q q ”m 瑶b 。q m c ， 一个重要的事实是，系统( 9 ) 可以线性化为 a 一a 玎 其中矩阵n 可以是任意非奇的n 阶方阵f 7 】。这意味着，如果矩阵磁b q 。非奇， 我们可阱将式( 1 1 ) 等价的复原为如下的二次形式 2 ，一a 只：a q 。一p t b q 。= 暖r ( a ) q 。 所以说，矩阵q 。和p m ，或者说二阶k r y l o v 子空间( 7 ) 和( 8 ) ，己经提供了充分 多的信息，我们完全可以使用它们来对原始系统进行投影，而不需要多余的线性 化步骤。这一点正是我们提出定义1 的初衷。整个过程中唯一的难点在于如何构 2 2 二阶双正交过程 这罩，我们给出本文的一个主要算法。这个算法可以构造二阶k r y l o v 子空阳j 的一组双正交基，也就是说，通过算法可以得到向量序列 0 成立 q ，= q ；屿+ + 1 ，j + 1 e 歹 因为我们要求初始向量q i = 0 ，那么 q = iq ；必 、。：。 如果二阶双正交过程不产生异常中断，那么上三角阵玩的对角元非零。也就是 说，矩阵玩的逆阵露1 存在。这样，x t j 于j = 1 ，2 ，m 成立 彰+ z = 西靠嘭+ 。 e j = q j 露1 e ， ( 2 2 ) 上面的等式揭示了嘭+ ，和q 2 _ f n n 关系，同时也提供了一个有效的计算嘭+ 。的 方法。将关系( 2 2 ) 应用到算法i 第5 行，可以得到 r = a + 8 q j 一1 【年二勺一l 同样的想法适用于式( 1 6 ) ，我们可以推出 b r 弓= j 癌瞒 v 2 1 。 v 2 j y 。 k 如果算法不产生中断，那么对于j = 1 ，2 ，m 成立 巧+ ，i 勺= b t 弓e q 勺 将上面的关系应用到算法i 第6 行，可以得到 s = a t p j + b t b 一1 踞牛1 通过式( 2 3 ) 和( 2 4 ) ，我们得到如下的精简存储的二阶双正交过程。 ( 2 3 ) ( 2 4 ) 1 4 复旦大学硕士学位论文 算法2 ：精简存储的二阶双正交过程 给定矩阵ab 础“，初试向量u ， 黔满足v t , i 0 ，以及m 1 。本算 法给出二阶k r y l o v 子空间舀。( a ，b ；“) 和g 。( a t ，b t ；u ) 的一组双正交基。 “，= u i “ q 1 = u 洞，p 1 _ s i g n ( 0 3 ) v 洞 ，= 0 ，g = 0 f o rj = 1 ，2 ，md o r = a 劬+ b f s = a t p j + b 丁g f o ri 一1 ，2 ，jd o u i j = p ，r = r q | l u i j v 玎。8 t q i ，s = 8 一p t v i j e n df o r u ：，r i f w = 0t h e nb r e a k d o w n = 洞，圳= s i g n ( w ) u j 扎j q j + 1 = = r 码+ 1 ，j ，p j + 1 = 8 v j + 1 , j f = q i o i l e j 。g = p 3 苛i l e 3 e n df o r 可以看到在算法2 中，我们不再需要显式的计算和保留辅助向量序列 “) 和 或) 。这样做带来的唯一开销是我们需要通过式( 2 3 ) 和( 2 4 ) 来隐式的求解西和 或。注意到矩阵和k 都是j 阶的上三角阵，我们可以使用向后回代法求解问 题彳1 勺和盯1 e ，。向后回代法的工作量为o ( j 2 ) ，因为通常j 不会很大，所以这 些额外开销相比较其他运算而言可以忽略不计。由此可知，相比算法1 ，算法2 在不增加额外开销的情况下，节省了一半的存储量，这一点对于求解大规模问题 而言是及其重要的。 1 5 l 互 置 出 曩 n z & 玑 n 挖 坞 m 垢 ” 复旦大学硕士学位论文 3 斜投影法求解二次特征值问题 在文献 1 4 中，作者使用二阶a r n o l d i 过程来构造二阶k r y l o v 子空间的一组 正交基，利用这组基可以采用正交投影的技术来求解二次特征值问题。在本文 中，我们将这个思想推广到采用斜投影的技术。 3 1 推广的二阶双正交过程 在给出求解二次特征值问题的斜投影方法之前，我们先绘出二阶双正交过程 的一个推广形式。在推广的算法中我们将构造子空间 9 。( 一m 一1 d ，一m 一1 k ；) 和毋。( 一m 一丁d t ，一m t k 丁；口)( 2 5 ) 的组m 一双正交基。也就是说，通过算法可以得到肉量序列 q j ) 和 辫 ，使得 s p a n q t ，q 2 ， = g 。( 一m d ，一m - 1 ；u ) ， s p a n p 1 ，p 2 ，p m ) = ( 一m t d 下，一m 一丁k t ； ) ， 并且 p t m q j = 屯对于i ，j = 1 ，2 ，m 算法3 ：推广的二阶双正交过程( g s o b ) 给定矩阵m ，d ，k r “，初试向量n ， 郾满足u t m u 0 ，以及 m 1 。本算法给出子空间( 2 5 ) 的一组m 一双正交基。 1 u = 口t m u 2 q l = 札厕，p 1 = s i g n ( w ) v v 阿l 3 ，= 0 ，g = 0 4 f o rj = l ，2 ，md o 5 r = 一m 一1 ( d q j + g f ) 6 8 = 一m t ( d r p j + k r g ) 7 f o ri = 1 ，2 ，jd o 8 珏订= p 朋一，r ；r 一哦“哲 9 v o = 8 t m q l s = s p i v o 1 6 复旦大学硕士学位论文 1 7 e n df o r u = s 7 m r i f u = 0t h e nb r e a k d o w n u z + l ，j = i u i ，+ 1 ，j = s i g n ( w ) u i + l ，j + 1 ：2r ，乱，+ 1 ，p j + 12 。s v j + 1 ，j f = q 3 0 i l e j 。g = p i 审i l e l e n df o r 关于算法3 的一些说明 在推广的二阶双正交过程的5 - 6 行，对于任意给定的右端项b ，我们需要求 解。满足 m x = b 或者m 7 0 = b ( 2 6 ) 当矩阵m 的阶数很大时，求解上面的线性代数问题将花费大量的时间。所 以，一个合理的处理方式是在开始执行算法之前。给出矩阵m 的某种分解 形式，例如：l u 分解，或者当m 对称正定时，c h o l e s k y 分解。这个分解 过程只需要进行一次，接着在算法运行过程中我们可以采用向前及向后回代 法求解问题( 2 6 ) 。由于求解三角方程组的工作量至多为o ( 舻) ，所以这样作 可以提高计算的效率。 当矩阵m ，o ，对称并且初始向量相同时，那么算法将得到q r r 。= f k 。也 就是说，对于对称问题，在整个循环过程中我们只需要使用一半的工作量就 可以得到需要的结果。即便如此，在对称的情形下算法仍然不具有短递推的 形式。 如同算法1 一样，推广的二阶双正交过程也会产生中断( 算法第1 2 行) 。只 不过，这时产生中断的判断条件为s t m r = 0 。在实际计算中，我们通常需 要预先设定一个闽值，当1 8 t m r l - c 时，即认为算法发生中断。 我们在表2 中列出了算法中涉及的各部分运算及其工作量。其中r l ( l u ) 表 示计算矩阵m 的l u 分解所涉及的浮点运算次数，肛( a ) 表示计算矩阵向量 乘法运算a z 或者a r x 的浮点运算次数，而( a ) 表示求解问题4 茁：b 或者 m n 他n m m 复旦大学硕士学位论文1 8 a t 。= b 的浮点运算次数。可以发现，对于大规模问题，算法的主要工作量 将集中在计算矩阵向量乘法( 尤其对于卢( m ) 的计算) 和求解关于l 和u 的 线性代数问题上。 r i ( l u ) 1 p ( m )( m + 1 ) 2 p ( d ) 十p ( k ) 2 m d l ) + v ( v ) 2 m ( 4 n 一1 ) ( m + 1 ) 2 + 2 n m ( m - t - 1 ) + m ( m + 1 ) ( 2 m + 1 ) 3 其他向量运算 = 6 n m 2 + l o n m + o ( m 3 1 表2 ：推广的二阶双正交过程涉及的运算量。 在算法中，我们采用了2 3 中提到的精简存储技巧。如果不考虑系数矩阵 m ，d ，k 以及l ，矿的存储，这个算法需要大约 的空间来保存必要的中间变量和最终得到的双正交基。 3 2 斜投影法求解二次特征值问题 由算法3 得到的矩阵q 。和p 仃。满足爿t 。m q 。= 。利用这个结果，我们可以 采用2 1 中提及的斜投影技术将原始的n 阶问题投影到m ( 礼) 阶问题来求解。 也就是说，我们需要求解口及非零向量f ，g 满足 q 。( 日) f = ( 口2 j + 目d 。+ k m ) ，= 0( 2 7 ) 和 g 日q m ( 口) = g h ( 日2 ，+ 日d m + k k ) = 0 ， 其中= 瑶d q 。和= 职k q 。分别为m xm 的实方阵。 如果已知( 目，f ，g ) 是m 阶二次a 一矩阵q 。( p ) 的真实特征组 系式 ( 日，z ，y ) = ( p ，q 。f ，只。g ) ( 2 8 ) 那么可以通过关 ( 2 9 ) 复旦大学硕士学位论文1 9 得到原来的n 阶二次a 一矩阵( 1 ) 的近似特征组( p ，x ，可) ，或者说是r i t z 组。这样 得到的r i t z 组对真实特征组的逼近程度可以通过残向量 r = ( 口2 m + o d + k ) x 和 8 h = y h ( 矿m + o d + k ) 来评价。 通过关系式( 2 9 ) ，易得右r i t z 向量对应的残向量满足r j _ s p a n p t 。 ，或者说 对于任意向量z s p a n r ，) ，；n t 等式 z r ( 目2 m + o d + k ) z = 0 成立。类似的，左r i t z 向量对应的残向量满足s l s p a n q 。) ，或者说对于任意向 量。s p a n q 。) ，如下等式 y 日( 日2 m + o d + k ) z = 0 成立。更进一步，我们可以得到斜投影法的简单向后误差估计。采用l a n c z o s 方 法求解非对称特征值问题的详细误差分析可以参考文献 2 8 ，高次特征值问题的 误差分析可以参考文献 4 】。由关系式( 2 9 ) ，可以推出残向量r 8 和r i t z 向量。，y 之间存在关系 日r = 0 和8 h x = 0 因此，如果令向后误差矩阵e 为 e = 蘸+ 淼 那么，易得等式 ( p 2 m + o d 十k e ) x = 0( 3 0 ) 和 y h ( p 2 m + o d + k e ) = 0 ( 3 1 ) 成立。当误差矩阵e 的f r o b e n i u s 一范数 ；= 牒+ 雌i l y l l i ( 3 2 ) 复旦大学硕士学位论文 充分小时，通过等式( 3 0 ) 和( 3 1 ) ，我们可以知道r i t z 组( 0 ，z ，y ) 是原先的二次 ，矩阵q ( a ) 加上一个小扰动矩阵之后的真实特征组。同样的，我们也可以通过 e 的范数来评估r i t z 组的准确性。不过需要注意的是，这里的误差矩阵依赖于特 定的r i t z 值和对应的左右r i t z 向量。 下面，我们给出采用斜投影技术求解q ( a ) 的特征值问题的算法框架。 算法4 ：斜投影法求解二次特征值问题 l ，给定矩阵m ，k ，d ，初始向量，v 通常可以随机生成) 和正整数m ，通过执 行推广的二阶双正交过程得到子空间f 2 5 ) 的一组肛双正交基。 2 将原始问题作斜投影，即计算矩阵口。= 碟d q 。和k 。= 碟k q 。 3 ，隶解投影后的饥阶二次入一矩阵的特征值阊题( 2 7 ) 和( 2 8 ) 。通过求得的特征 组( p ，g ) ，由关系式( 2 9 ) 求出原来问题的r i t z 组( 0 ，9 5 ，) ，并规范化 x = q 。f l l q 。刑2 和y = f m 9 1 1 p m 圳。 4 通过相对残量范数 瑞箫箫鼎和揣器薷群( 3 3 0 1 1 1 d t lo i ) i d i ， 2 m 1 1 1 十l，+ 1 1 k 1 1 - ” 1 p 俐 彳1 + 11 + 1 i k l l l 7 来评估r i t z 组的准确性。 关于算法4 的一些说明： 正如我们在31 中所指出的，当矩阵m ，d ，k 对称并且初始向量相同时， 那么算法3 将得到q 。= r 。这时，作斜投影后的矩阵d 。，k 。也是对称 的。这意味着，q 。( a ) 所具有的特殊结构以及隐含的谱性质将显式的在投影 后的m 阶系统中保留。 为了求解投影后的m 阶二次a 一矩阵的特征值问题，我们采用51 中提到的直 接法。也就是说，首先将问题线性化为一个一般矩阵的特征值阀题，然后采 用成熟的方法，例如：q r 算法，来求解投影后系统的所有特征值和对应的 特征向量。通常来说，投影后的阶数m 要远小于原始问题的阶数n ，因此通 过求解投影问题来逼近原始问题可以极大的减少计算的工作量。 2 0 复旦大学硕士学位论文 2 1 对于二次特征值问题，我们必须通过求出相对残量范数来判断r i t z 组是否 已经收敛到真实特征组。这一点不同于求解一般矩阵特征值问题的a r n o l d i 方法，可以通过上h e s s e n b e r g 矩阵的次对角元来判断结果是否已经收敛。 我们知道在使用k r y l o v 子空间方法求解标准的特征值问题时，正交性的保 持对求解结果的好坏会带来很大的影响，对于二次问题同样如此。在使用推 广的二阶双正交算法求m 一双正交基的过程中，可能会出现正交性丧失的情 况。为了得到高精度的解，我们必须增加重正交化的过程。当然我们也可以 在投影步骤中，显式的产生矩阵a = 碟m q 。，而不是使用单位阵j 。为 了克服由于正交性丧失所带来的不良影响，在执行了固定步长的迭代之后对 算法进行重启 2 9 l 是一个有效的手段。但是，如何在二次特征值问题的求解 过程中采用重启技术f 3 0 1 仍是一个重要的研究课题。 在本文的算法中，我们构造二阶k r o l o v 子空间( 2 5 ) 的一组m 一双正交基来 对原问题进行投影。当然，我们也可以采用其他子空间的基来作投影。比如 说，当矩阵m 对称正定时，有c h o l e s k y 分解m = l 工t 。这时，可以使用 算法2 来构造二阶k r y l o v 子空间 ( a ，b ；) 和( a 丁，b r ； ) 的一组双正交基，其中 a = 一l 一1 d l t b = 一l 一1 k l t 然后，我们使用这组双正交基对原问题进行投影，这样得到的投影系统也将 显式的保留原问题所具有的谱性质。只不过，这个方法得到的r i t z 组和真 实特征组之间存在关系 ( d ，z ，) _ ( 0 ，l - t q 。，l - t p m 9 ) 采用不同的子空间进行投影，需要的存储量、浮点运算次数和得到的结果都 不会完全相同，在这里我们不作迸一步的讨论。 复旦大学硕士学位论文 2 2 4 数值例子 在本节中，我们给出几个数值例子来说明采用斜投影的方法，也就是算法 4 ，可以近似的求得二次特征值问题的解。同时，我们将计算得到的结果( 在图中 用g s o b 表示) 同问题的真实解( 使用m a t l a b 函数p o l y e i g 得到的结果，在图 中用e x a c t 表示) ，以及使用文献f 1 4 1 中提出的二阶a r n o l d i 方法求得的结果( 在 图中用s o a r 表示) 进行比较。在下面的算例中，如果不加特别说明，都采用各 个分量为1 的向量作为初试向量，用1 0 。o 作为中断判断的闽值。 算例1 ( 随机非对称系统) ：在这个算例中，m ，d ，都是随机生成的2 0 0 阶非对称 实矩阵，矩阵的元素服从均值为0 ，方差和标准方差为l 的正态分布。图1 显示 了系统的真实特征值，以及当m = 2 0 时，由g s o b 和s o a r 计算得到的近似特 征值。在图2 中显示了由式( 3 3 ) 计算得到的左右相对残量范数( 分别用g s o b l e f t 和r i g h t 表示) ，同时作出了由式( 3 2 ) 计算得到的向后误差矩阵的范数( 用e r r o r m a t r i x 表示1 。由二阶a r n o l d i 方法得到的r i t z 对的相对误差也在图中给出。从图 中我们可以看出，左右r i t z 向量的收敛情况几乎一致，并且r i t z 组是通过对原始 问题作一个小扰动之后得到的真实解。 1 5 f ； i 卜，一毒辫争一 。卜”嚣瀚瓤妒一 一5 铲甲 1 。 佃r 匿 r e a ip a r t 图l ：特征值图2 ：误差 为了说明通过算法4 可以求得近似特征值，并且通常是具有最大模的那些特 征值，我们在图3 和4 中分别作出了模最大的特征值( a = 一3 2 2 - t - 1 5 5 i ) 和模次 大的特征值( a = 0 3 0 74 - 9 4 3 i ) 的收敛情况，图中的横坐标表示m 的取值从l 到 em己m妄口璺t一 复旦大学硕士学位论文 5 0 ，纵坐标表示特征值的相对误差，即l a e l lj ) , l 。由s o a r 得到的收敛情况同 时显示在图中。我们发现，对于这个例子而言，采用算法4 可以更快的收敛到真 实特征值，并且对于模大的特征值，其收敛的速度也越快。如果收敛到相同的精 度，相比s o a r 算法，斜投影法只需要一半的迭代次数。 图3 ：模最大特征值 01 02 03 04 05 0 r e d u c e do r d e r 图4 ：模次大特征值 算例2 ( 随机g y r o s c o p i c 系统) ：这个算例采用的是g y r o s c o p i c 系统，这是一类十 分重要的系统，并且已经得到广泛研究，可以参考文献5 ，7 1 。在该系统中系数 矩阵具有很好的性质，矩阵m 对称正定，d 斜对称，k 对称负定。由表1 易 知，g y r o s c o p i c 系统的特征值或者为实数或者关于实轴和虚轴对称。在这里我们 随机生成2 0 0 阶实矩阵作为算例。图5 显示了真实特征值，以及当m = 2 0 时，由 g s o b 和s o a r 计算得到的近似特征值。在图6 中显示了由式( 3 3 ) 计算得到的相 对残量范数，以及由式( 3 2 ) 计算得到的向后误差矩阵的范数。同算例1 ，我们可 以得出结论，左右r i t z 向量的收敛情况几乎一致，并且r i t z 组是通过对原始问题 作一个小扰动之后得到的真实解。 虽然这里的矩阵d 是斜对称的，但是如果取相同的初试向量，那么通过算法 3 仍然可以得到q 。= f k 。因此，通过投影得到的系统将保持原先的矩阵特性， 更进一步，隐含的特征值属性也将显式的保留。正如在图5 中显示的一样，采用 投影算法求得的近似特征值或者为实数或者关于实轴和虚轴对称。 复旦大学硕士学位论文 a d p r o x i m a t ee i g e n v a l u e s 审 v ，一一扣州一v ；一v ： j e x a c ti xg s o bi 】审s o a ri v：v r e a id a n 图5 ：特征值 i n d e x 图6 ：误差 算例3 ( n a s t r a n r b 的算例) ：这个算例来源于n a s t r a n 。n a s t r a n 【3 1 】是一个有限元 数值包，其中也涉及到需要求解结构力学中碰到的二次特征值问题。在这个算例 中系数矩阵m ，d ，k 是3 6 0 0 阶的稀疏对称阵。关于矩阵的一些性质在表3 中列 出，其中表格的最后一列为使用m a t l a b 函数c o n d e s t 估算的矩阵1 一范数条件 数的一个下界。 非零元个数1 一范数条件数 m5 5 2 13 6 0 0i n f d1 9 5 7 01 0 2 5i n f k5 9 0 6 22 1 9 1 0 1 28 4 2x1 0 6 表3 ：算例3 的矩阵性质。 在这个算例中，我们使用了求解特征值问题常用的带位移的反幂技巧 6 ，1 4 。也就是说，我们需要求解如下的新的二次a 一矩阵的特征值问题 矿砑+ 肛西+ 露，( 3 4 ) 这里 m = 口2 m + a d + k d = d 2 a m 詹：m eoc$；日一2 elog 复旦大学硕士学位论文 并且 1 p 2 j i 石 在这个算例中我们设定位移盯= 1 0 4 。通过求出问题( 3 4 ) 的特征值，我们可以得 到原始问题的特征值为 a = 口+ 二 肛 这样，只要的模越大，求出的特征值就越接近位移一。 由于这个问题的规模较大，我们无法使用p o l y e i g 函数求得其真实的特征 值，所以也就无法用图示比较近似特征值的逼近程度。在图7 中，显示了当 m = 5 0 时，由g s o b 和s o a r 得到的相对残量范数。可以看到，虽然这个问题 的规模较大并且系数矩阵是坏条件的，但是两种算法都可以得到较好的近似结 果。 i n d e x 图7 ：误差 2 5 复旦大学硕士学位论文 2 6 5 小结 在本文的开始，我们提出了二阶k r y l o v 子空间的定义。这子空间的构造需要 基于两个矩阵和初试向量，可以看作是对标准k r y l o v 子空间的一个推广。接着， 我们给出了二阶双正交过程的算法，通过这一过程可以得n - - 阶k r y l o v 子空间的 一组双正交基。除了不具有短递推的性质外，二阶双正交过程的输入、输出以及 从中得到的矩阵关系都类似于标准的非对称l a n c z o s