矩阵特征值求解

矩阵特征值求解的分值算法

12组

1．1 矩阵计算的基本问题

(1)求解线性方程组的问题.即给定一个n 阶非奇异矩阵A 和n 维向量b ,求一个n 维向量x ,使得

b Ax = （1.1.1）

(2)线性最小二乘问题,即给定一个n m ?阶矩阵A 和m 维向量b ,求一个n 维向量

x ,使得

},min{n R y b Ay b Ax ∈-=- （1.1.2）

(3)矩阵的特征问题,即给定一个n 阶实(复)矩阵A ,求它的部分或全部特征值以及对应的特征向量,也就是求解方程

x Ax λ= （1.1.3）

一对解(λ,x ),其中)(),(n n C R x C R ∈∈λ,即λ为矩阵A 的特征值,x 为矩阵

A 的属于特征值λ的特征向量。

在工程上,矩阵的特征值具有广泛的应用,如大型桥梁或建筑物的振动问题:机械和机件的振动问题；飞机机翼的颤振问题;无线电电子学及光学系统的电磁振动问题;调节系统的自振问题以及声学和超声学系统的振动问题.又如天文、地震、信息系统、经济学中的一些问题都与矩阵的特征值问题密切相关。

在科学上,计算流体力学、统计计算、量子力学、化学工程和网络排队的马尔可夫链模拟等实际问题,最后也都要归结为矩阵的特征值问题.由于特征值问题在许多科学和工程领域中具有广泛的应用,因此对矩阵的特征值问题的求解理论研究算法的开发软件的制作等是当今计算数学和科学与工程计算研究领域的重大课题,国际上这方面的研究工作十分活跃。 1.2 矩阵的特征值问题研究现状及算法概述

对一个n n ?阶实(复)矩阵A,它的特征值问题,即求方程(I.1.3)式的非平凡解,是数值线性代数的一个中心问题.这一问题的内在非线性给计算特征值带来许多计算问题.为了求(l.1.3)式中的λ,一个简单的想法就是显式地求解特征方程

0)det(=-I A λ （1.2.1）

除非对于个别的特殊矩阵,由于特征方程的系数不能够用稳定的数值方法由行列式的计算来求得,既使能精确计算出特征方程的系数,在有限精度下,其特征多项式)det()(I A f λλ-=的根可能对多项式的系数非常敏感.因此,这个方法只能在理论上是有意义的,实际计算中对一般矩阵是不可行的.首先,若矩阵A 的阶数较大,则行列式)det(I A λ-的计算量将非常大;其次,根据Galois 理论,对于次数大于四的多项式求根不存在一种通用的方法,基于上述原因,人们只能寻求其它途径.因此,如何有效地!精确地求解矩阵特征值问题,就成为数值线性代数领

域的一个中心问题.

目前,求解矩阵特征值问题的方法有两大类:一类称为变换方法,另一类称为向量迭代方法.变换方法是直接对原矩阵进行处理,通过一系列相似变换,使之变换成一个易于求解特征值的形式,如Jacobi 算法,Givens 算法,QR 算法等。变换方法由于要存储矩阵元素,因而它只适用于求解中小型矩阵,它一般和向量迭代方法结合起来使用.向量迭代方法是通过一系列矩阵向量乘积而求得特征值和特征向量的.由于向量迭代方法可采用压缩存储技术,因而它适合于求大规模矩阵的特征值问题,尤其是大型稀疏矩阵的部分特征值和特征向量问题,如Lanczos 方法,Arnoldi 方法,Davidson 方法等,现在这类问题仍是比较热的研究课题。 2 分治方法的基础及理论研究 2.1 分治方法的概述

考虑对称三对角矩阵n T 的特征值问题

x x T n λ= （2.1.1）

其中

???

?

??

?

?

?=--n n n n T αββαββα1

12

111

（2.1.2） 1981年Cuppen 提出一种求上述对称三对角矩阵n T 所有特征值和特征向量的分而治之方法（divide 一and 一conquer method).其基本思想是先将对称三对角矩阵n T 分割为两个分别为k k ?阶和k)-(n k)-(n ?阶低阶对称三对角子矩阵

)0(T 和)1(T .)0(T 和)1(T )可以用同样的方法也分别分割为两个更低阶的子矩阵,递

归的采用这种分割技术可以把矩阵分割为一些能直接求出特征值的足够小的子

矩阵(比如2阶或1阶矩阵),或者按照某种标准分割到适当阶数(如小于等于25阶)后,结合其它求矩阵特征值的方法,如QR 算法,求出其特征值。在求出低阶矩阵特征值的基础上,开始胶合过程。在胶合阶段,分割前的矩阵1T 的特征值的求

出(所谓的“治之”)是建立在其两个子矩阵'

)0(T 和')1(T 的特征值的基础上的,其中'

)0(T 和')1(T .是在分割阶段由1T 分割出的低阶子矩阵.随后的数值分析表明,Cuppen 的方法存在着数值不稳定的危险,特别是当存在特征值束时,计算出的特征向量可能不正交。Gu 和Eisenstat 对Cuppen 的方法作了改进,极大地降低了数值不稳定的危险性。

Cuppen 的方法在计算n T 的特征值的同时也需要计算对应的特征向量,并且

是在)0(T 和)1(T 的特征值和特征向量的基础上进行计算的.根据文中,当用残量

X X T n Λ-和正交性n T I X X -作为检验准确性的标准时，Cuppen 的方法比二

分法或多分法精确的多。但文[3中，,如果用n T λ

λ-~

作为衡量特征值准确性

的标准时,二分法或多分法精确些.此外Cuppen 的分而治之方法要求矩阵乘积,存储量为)(2n O ,而二分法或多分法的存储量仅为)(n O ,比前者少.应此,当只需计算特征值时,通常选用后者.1987年Dongarra 和Sorensen 把分治思想应用到求对称三对角矩阵特征值的并行计算,取得了不错的效果,再次引起了人们对分治方法的极大关注。

分割胶合方法（split-merge method ）是后来提出用分治方式求对称三对角矩阵特征值n T 的方法,不同于分而治之方法在分割过程中采用矩阵的秩1扰动,它采用矩阵的秩2扰动.与二分法和多分法相似,在分割胶合方法中,特征值的计算独立于特征向量的计算.如果需要计算特征向量,可采用反幂法。由于分割胶合方法计算特征值时采用具有三次收敛的Laguerre 迭代,数值试验表明，其计算速度和精度都明显优于二分法和多分法.并且文中给出的用于计算Laguerre 迭代的非线性三项递归式可以避免上溢和下溢问题。 2.1.1 分割策略

分治方法的第一步就是把原来高阶的对称三对角矩阵特征值问题转化为两个低阶对称三对角矩阵特征值问题,即所谓的分割阶段.设对称三对角矩阵n T 表示如下

???

?

??

?

?

?=--n n n n T αββαββα1

12

111

（2.1.1.1） 不妨假设)1,,2,1(0-=≠n j j β,即称n T 为不可约矩阵。否则,若存在某些0=j β,则n T 就可以约化为若干个低阶主子矩阵特征值问题, n T 的特征值就由这若干个低阶主子矩阵的特征值构成.当n T 为不可约对称三对角矩阵时,对其作如下分割

E T T T n +???

?

?

?=)1()

0( (2.1.1.2)

记E n T T n +=~

,其中)0(T 和)1(T )分别为k k ?阶和)()(k n k n -?-或

)1()1(--?--k n k n 阶对称三对角矩阵,通常]2/[n k =。

（1）当T

vv O O E ρρθρθρθρ=???

???

?

??=2

时，其中），，，，，010( θ=v 且k βρθ=，称此为秩1扰动矩阵.此时有

??

?

??

??

?

?-=--ραββαββαk k k T

11211

1)

（，??????

? ?

?-=--++++n n n k k k k T αββαββρθα1

12

11

21)

1

（

（2）当T

k k T k k e e e e O O

E 110

0+++=??????

?

?

?=ρρρρ时，其中k βρ=，称此为秩2扰动矩阵.此时有

???????

?

?=--k k k T

αββαββα112111)

（，???

???

?

?

?=--++++n n n k k k k T αββαββα1

1211

1)

1

（ （3）当?????

??

?

?

?=+++0011

1

k k k k k E ββαββ时，称此为秩3扰动矩阵.此时有k k k k k R T ?--∈???

????

??=αββαββα1

121

1

1)

（，

)1()1(1

13

22

2)

1--?----++++∈??????

? ?

?=k n k n n n n k k k k R T αββαββα

（

人们关注的问题是:对于上述三种分治策略n T 与n T ~

的特征值之间的关系.利用Hoffman-Wielandt 定理的结论,我们可以得到如下定理:

定理2.1设A 和B 是两个n 阶Hermite 矩阵,它们的特征值分别是n λλλ≥≥≥ 21和n μμμ≥≥≥ 21，则

F

n

j j j B

A -≤-∑=2

112])([μλ (2.1.1.3)

其中F ?为Frobenius 范数。

根据上述定理,设n T 与

n

T ~

的特征值分别为n λλλ≥≥≥ 21和

n λλ

λ~

2

~

1

~≥≥≥ ，则有下面关系成立

F j E j ≤-λλ~

（2.1.1.4）

因此,只要(2.1.1.4)式右端越小,j λ就越接近j λ~

,即n T 与n T ~

的特征值就越接近。对于秩1和秩2扰动, n T 与n T ~

的特征值之间关系更详细的描述,将在后文给出。

2.1.2 胶合

在矩阵n T 分割后,先求出矩阵n T ~

的特征值,即求出)0(T 和)1(T 的特征值。剩下的工作就是如何由n T ~

的特征值出发求出n T 的特征值,这就是胶合阶段主要任务。在这个阶段可以采用不同的迭代方法,如割线法!Newton 迭代、Laguerre 迭代以及路径跟踪等,以n T ~

的某个特征值或n T ~

的特征值构成的区间内的点为初始点经过若千步迭代,最后收敛到n T 的某个特征值.本文中采用Laguerre 迭代从特征区间提取特征值。 2.2 Laguerre 迭代

2.2.1 Laguerre 迭代及其计算

根据文中,对于不可约对称三对角矩阵n T ,它的特征值或特征多项式

)det ()(n n I T f λλ-=的根全为互不相同的实数.因此,适合求多项式的根都是实

单根的具有三次收敛的Laguerre 迭代

)])

()

(())()()(1)[(1())()(()(2x f x f n x f x f n n x f x f n

x x L ''-'---+'-

+

=--

+ (2.2.1.1)

非常适合用来从特征区间提取出n T 的特征值.这个方法早在13世纪就已经提出,近年来结合分治策略又重新焕发出生机.文中首次对用Laguerre 迭代求不可约对称三对角矩阵的特征值的实用性进行了研究,而后又被用于求矩阵的奇异值问

题。

为了利用Laguerre 迭代求n T 的特征值,我们需要计算n T 的特征多项式

)det()(n n I T f λλ-=以及其一阶导数)(λf '、二阶导数)(λf ''.由文中知,计算这些值的有效的方法是三项递归式.众所周知,当矩阵阶数比较大时,三项递归式可能出现上溢或下溢问题.文中给出了一种改进的非线性的三项递归式,数值试验表明,除极其个别的情形,改进的三项递归式可以避免上溢或下溢问题。

对称三对角矩阵n T 的特征多项式及其一、二阶倒数的计算公式

设对称三对角矩阵n T 如（2.2.1.1）式所示, k T 表示其k 阶顺序主子式,表示如下

???

?

??

?

?

?=--k k k k T αββαββα1

12

111

， （2.2.1.2） )det()(k k k I T f λλ-=为k T 的特征多项式,则特征多项式及其导数计算公式如下(三项递归式):

???=--=-==---n i f f f f f i i i i i

,,3,2),()()()()(,1)(22

11110 λβλλαλλ

αλλ （2.2.1.3） ???='-'-='-='='---n

i f f f f f i i i i i ,,3,2),()()()(1

)(,1)(22

1110 λβλλαλλλ （2.2.1.4） ???

=''-''-''-=''=''=''----n i f f f f f f i i i i i i

,,3,2),()(2)()()(0)(,1)(22

11110 λβλλλαλλλ （2.2.1.5） 及

),()(),()(),()(λλλλλλn n n f f f f f f ''='''='=

为了避免上述(2.2.1.3)式、(2.2.1.4)式和(2.2.1.5)式在计算中可能出现的上溢或下溢问题,文中给出如下的改进的非线性三项递归式: 令

,,,2,1,)

()

()(1n i f f i i i ==

-λλλξ (2．2．1.3)式两边都除以得)(1λ-i f 得

??

?

??=-==--n i i i i i ,,3,2,)(--121

11 λξβλαλξλαλξ）（）

（ （2.2.1.6） 令

)()()(λλληi i i f f '-

=和n i f f i i i ,,1,0)

()

()( ='-=，λλλ? (2．2.1.4)式和(2。2.1.5)式两边除以)(1λ-i f 并注意下述关系

)()(1

λξλj j i f ∏==，

可得下面的递归关系

???

????

=-+-===----n

i i i i i i i i ,,3,2)],())((1)()[()(1)()(1)(,0)(212

11110 ληλξβληλαλξληλξληλη

（2.2.1.7）

??

???=-+-===----n

i i i i i i i i ,,3,2)],())((1)()[()(1)(0)(,0)(212

1110 λ?λξβλ?λαλξλ?λ?λ?

（2.2.1.8）

分别对(2.2.1.7)式和(2.2.1.8)式进行简单的推导,我们可以得到下面两个关系

式)()

()

(),()()(-λ?λλληλλn n f f f f =''=' 我们称(2.2.1.6)式为改进的三项递归式.(2.2.1.7)式和(2.2.1.8)式分别为改

进的计算特征多项式一阶导数和二阶导数的计算公式,这些是我们在调用

Laguerre 迭代时所需要计算的.根据文中知,除了个别极其特殊的情形,它们可

以避免计算过程中出现的上溢或下溢现象.但是,值得注意的是,如果对于某个

)1(n i i ≤≤使得0)(=λξi 时,上述改进的各式就会出现中断现象.为了保证排除中

断的发生,在计算中作如下处理

如果0)(=λξi ，则令1,,2,1.)(-==n i i i εβλξ 如果0)(=λξn ，则令εβλξ1)(--=n n

这里ε为机器精度.上述处理相当于对i α或n α引入一个小扰动εβi 或εβ1-n 。 注 2.1:我们由三项递归式(2.2.1.3)式经过改进得到(2.2.1.6)式,在计算改进的非线性三项递归式(2.2.1.6)式时,同时可以得到在某个给定值λ处,不大于该固定值λ的特征值的个数,即0)(<λξi (对所有n i ≤≤1)的个数. 2.2.2 Laguerre 迭代的性质

对于不可约对称三对角矩阵n T (如(2.1.1.1)式),其特征多项式

)det()(n n I T f λλ-=， （2.2.2.1）

的所有根为互不相等的实根。设)(λf 的根为

n λλλ<<< 21

并令-∞=0λ和+∞=+1n λ。对),(10+∈n x λλ，Laguerre 迭代)(x L +和)(x L -分别定义如下

)])

()

(())()()(1)[(1())()(()(2x f x f n x f x f n n x f x f n

x x L ''-'---+'-

+

=+ （2.2.2.2）

)])

()

(())()()(1)[(1(-))()(()(2-x f x f n x f x f n n x f x f n

x x L ''-'---'-

+

= （2.2.2.3）

并且有下面的命题成立.

命题2.1 对),(1+∈m m x λλn m ,,1,0 =，则有下述结果成立: （1）1)()(++-<<<

当x 趋于1+m λ时，有3

11)(++++-≤-m m x c x L λλ； 当x 趋于m λ时，有3

)(m m x c x L λλ-≤---。 根据命题2.1,如果令

)))((()()(

x L L L x L x k

k

k +++++==

)))((()()(

x L L L x L x k

k

k -----==，

则有

1)

2()1()1()2(+++---<<<<<-m m x x x x x λλ （2.2.2.4）

(2..2.2.4)式说明当),(1+∈m m x λλ时, )(x L -从x 的左侧趋向于m λ,)(x L +从x 的右侧趋向于1+m λ.文中指出只要x 非常接近m λ (或非常接近1+m λ),一般调用(2..2.2.3)式(或(2..2.2.2)式)两!三次便可收敛于m λ(或1+m λ)。 2.2.3 病态接近特征值及改进的Laguerre 迭代

根据文献,当多项式有重根时Laguerre 迭代只是线性收敛.虽然不可约对称三对角矩阵的特征值为两两互不相同的实数,但是,它们中的某些特征值也可能

非常接近以致于在数值上无法区分,例如Wilkinson 矩阵。当矩阵有一个由r 个非常接近的特征值构成的束时,在计算中被看作为一个r 重特征值,在这种情形下,Laguerre 迭代的收敛率仅为线性的。

为了在上述情形下加速迭代Laguerre 收敛,对Laguerre 迭代(2..2.22)式和(2..2.23)式作重新定义，对正整数r （n r ≤<0），定义)(x L r +和)(x L r -如下

)])

()

(())()()(1[())()(()(2x f x f n x f x f n r r n x f x f n

x x L r ''-'---+'-+

=+， （2.2.3.1）

)])

()

(())()()(1[())()(()(2x f x f n x f x f n r r n x f x f n

x x L r ''-'----'-

+

=-， （2.2.3.2）

我们称(2.2.3.1)式和(2.2.3.2)式为改进的Laguerre 迭代.我们需要注意的是:当1=r 时,(2.2.3.1)式和(2.2.3.2)式就分别是(2.2.2.2)式和(2.2.2.3)式.因此,我们在调用Laguerre 迭代时,通常是先令1=r ,然后计算r 来调用(2.2。3.1)式或(2.2.3.2)式。

由上面的两个式子,可以得到两个序列

)))((()()(

x L L L x L x k

r r r k

r k r +++++==

)))((()()(

x L L L x L x k

r r r k

r k r -----

==

当*x 为)(λf 的r 重根时,如果*

)

(k r x +, ,2,1=k 以三次收敛速度递增地收敛于*x ;类似地有,如果x x <*并且

)(λf 在区间),(x x *内无根,则序列 ,2,1

,)(=-k x k r ,以三次收敛速度递减地收敛于*x 。数值试验表明,即使在病态根束的情形下,(2.2.3.1)式和(2.2.3.2)式也能

达到三次收敛.对上述改进的Laguerre 迭代)(x L r +和)(x L r -,有下面定理: 定理 2.2设n λλλ<<< 21为多项式)det()(n n I T f λλ-=的根,则对

),(1+∈m m x λλ和任一正整数n r <,我们有 (1)如果,0)

()

(->'x f x f 则r m r m x L x L x L x ++++<<<<<<λλ)()()(2 ； (2)如果,0)

()

(-

<'x f x f 则121)()()(+---+-<<<<<

注2.2：设),(1+∈m m x λλ为Laguerre 迭代的初始点并用它逼近特征值1+m λ,在某些情形下,在离1+m λ右边非常近的地方可能存在一组特征值,即

r m m m +++<<<λλλ 21它们相对于x 离1+m λ的距离,非常靠近1+m λ。

在这种情形下,在达到1+m λ的能够具有三次收敛速度的邻域前,实际计算需要花费许多次调用La g u e rre 迭代.这时,我们可以选择其它初始点来重新开始

Laguerre 迭代或者调用改进的Laguerre 迭代)(x L r +来加速收敛过程。

3 求反对称三对角矩阵特征值的方法及分而治之算法 3.1 反对称三对角矩阵的算法介绍

反对称矩阵特征值问题在许多实际问题中有广泛应用,如防潮陀螺仪系统的模型分析【23】等,对这类具有特殊结构的矩阵,已有不少求特征值的算法【21-25】它们基本上是把反对称矩阵特征值问题经过正交相似变换,如Housheuodl 变换或Gviens 变换,转化为反对称三对角矩阵特征值问题.众所周知,在相似变换下不改变原矩阵的谱集.对一般矩阵,在采用QL 算法之前,通常也是先把原矩阵化为上

Hessnehegr 型,当矩阵为反对称矩阵时,根据其反对称性,其上Hessnehegr 型矩阵

也就是三对角矩阵.因此,我们在下面讨论的反对称矩阵总是已三对角化了的反对称三对角矩阵,形如

n n n n R J ?--∈???

???

? ?

?--=0001

11

1ββββ

（3.1.1）

不失一般性,假设,1,,2,1,0-=≠n j j β式称之为不可约反对称三对角矩阵，简记为LATS,且近一步假设,1,,2,1,0-=>n j j β若存在,10,0-≤≤=n j j β则原问题可以转化为两个低阶矩阵的特征值问题。

3.1.1 求IAST 特征值的QL 算法

考虑IAST （3.1.1）,如果矩阵阶数n 为奇数,由于IAST 的特征值为互为共扼的纯虚数或零,因此可以通过位移为零的QL 迭代把n J 化为

????

?

?=-10n n J J （3.1.1.1） 其中1-n J 为偶数阶矩阵,其特征值为互为共轭的纯虚数。不失一般性,我们以后不妨假设矩阵的阶数n 为偶数。令J J n =,计算j 的特征值的QL 算法描述如下: （1）计算I J 212β+的第n 列

()()

T

n n n n e I J a 1

22

121212,0,,,,0,0----=+=βββββ

（2）计算er Househould 矩阵0P 使得

,0n e a a P =

这里,20T

ww I P -=其中(

)

h w T

n n n /,0,,,,0,0122

121τ

ββββ--=---

()

()(

)[]

(

)

1

22

11

22

12

1221

1

22

222

1

22112222)(2----------=--=--+=-+=n n n n n n n n sign sign h ββτββτττ

ββ

ββββββτ

（3）计算,00~

JP P J T

=即

??????????

?

?

???

-=0*

0**0*00*

0**0

*00011~

ββJ 然后利用一系列旋转变换把~

J 化为新的LAST

具体地,注意~

J 只是比原矩阵J 在(n,n-3)和（n-3,n ）位置上多两个非零元素.用旋转矩阵

??

?

???

?

?

?

?-=100100c 4-n 1c s s I P

把21~

1

2P P J P P T T

的(n,n-4)和（n-4,n ）位置上多两个非零元素，这时在(n,n-5)和（n-5,n ）位置上多两个非零元素。类似地,通过旋转矩阵3-n 21P P P ，，

， ，把~

J 化为 321~

1

23-n -∧

=n T

T

T

P P P J P P P J

此时J

?为一个新的LAST 。此时对J ?进行压缩。在对J ?的压缩后的矩阵执行双位移QL 算法，经过若干步上面的循环后，最终可以收敛于J 的标准型[25]

()'

2

2/'212,,I I diag P n λλ =, 其中

???

?

??-=0110'2

I .

3.1.2 由LAST 的叉积的压缩矩阵求其特征值的方法

1971年，M.H.C.Paardekooper 曾用于对称矩阵的Jacobi 思想推广到反对称矩阵情形，给出了一种求反对称矩阵特征值得算法。但他的算法仅考虑了反对称性，未利用其纯虚数特征值共轭成对的性质，廉庆荣教授等利用这一性质，构造出矩阵阶数为N n 2=的T JJ （其中J 为LAST ）的压缩矩阵（N 阶不可约对称三对角矩阵），简记CIST (condense irreducible symmetric tridiagonal matrix)。下面的定理可以说明：CIST 的特征值的开方是对LAST 的特征值的虚部。有关

算法是在构造出T JJ 的CIST 后，采用对称QL 算法（带wilkinson 位移）求出CIST 的特征值，然后在对其特征值开方求得LAST 特征值的虚部。

不失一般性，仍考虑形如（3.1.1）式的n 阶不可约反对称三对角矩阵J （不

妨设N n 2=为偶数，n 为奇数时，[]2n N =）。令()

*

*==ij T J JJ α,*J 必为特殊的

五对角对称矩阵：

???????=-=-==-=+===*-*

+-*-+-*-**.

,0,

1,,3,2,,1,,3,211,11,12

21212111其他，，，ij j j j j j j j j jj n nn n j n j αββααββαβαβα （3.1.2.1） 定理 3.1设T JJ J =*如（3.1.2.1）所示，若选择排列阵()

n n ij R p P ?∈=满足

0,,,2,11,1-2,2====+ij j N j j j p N j p p ，(其它)，则T T P J P T ~⊕=*,其中T 和T ~

均

为不可约对称三对角矩阵，()

N N ij R T ?∈=α且*=j i ij 2,2αα，()

N N ij

R T ?∈=α~~且()N j i j i j i ij ,,2,1,,1

~12,12 =≤-=*--αα. 注3.1 当12+=N n 时，T ~

为1+N 阶，定理仍成立。

本文仍采用文[22]中的叫法，称T 和T ~

为*J 的压缩矩阵，其元素可以从*J 中

直接得出，以T 为例：

?

??

??=-=-==-=+==--++-.

,0,1,,2,1,,1,,2,112121,,12221222其他，，ij j j j j j j j j jj N NN N j N j αββααββαβα （3.1.2.2） 下面的定理给出定理3.1中T 的特征值与J 的特征值的关系。

定理 3.2[设J 为n 阶LAST 如（3.1.1）式所示，T 为T JJ J =*的N 阶压缩矩阵（如（3.1.2.2）式所示），则j i μ±为J 的共轭对特征值的充分必要条件是2j μ为

T 的特征值（N j ,,2,1 =）.

由定理3.1易知，2j

μ）（N j ,,2,1 =也是T T ~⊕的二重特征值，而且T 和T ~

各有N 个相异非零实特征值。于是T T ~

⊕有相同的特征值2j μ）（N j ,,2,1 =。 该定理表明，只要求出N 阶对称三对角矩阵T 的特征值，即可得n 阶不可约

反对称三对角矩阵的特征值。在求N 阶对称三对角矩阵T 的特征值时，采用带有wilkinson 位移的QL 算法。

3.2 求LAST 的特征值的分而治之算法 3.2.1 秩1分割

由定理3.1和定理3.2知，求LAST 得特征值就等价于求N 阶对称三对角矩阵T 的特征值，其中T 为T JJ J =*的压缩矩阵。文[22]是用带有wilkinson 位移的QL 算法求T 得特征值本节将对T 利用分治策略，给出一种求不可约对称三对角阵T 的特征值的分而治之算法。如（3.1.2.2）式所示，T 形如

??????

? ??=---N N N N T αγγαγγα111111 , (3.2.1.1)

其中()1,,2,101,-=≠?+N j j j j αα.对T 作分割策略，即对T 作如下的秩1扰动

()()E T E T T T +?+???? ?

?=~

10, (3.2.1.2) 其中

()()???? ??==∈????????

?

? ??=∈=???????

?

?=?102

~,010,T T T R v R vv O O

E k N

N N T βρθθρρθρθρθρ且 ()

??

???

??

??-=---ραγγαγγαk k k k T 1111110 ，()??????

? ??-=---+++N N N N k k k T αγγαγγρθα11111211, k 一般取为[]2N .

这样，把T 分割为一个分块矩阵和一个秩1矩阵的和。如果需要，还可以对

()0T 和()1T 分别进行如上的分割，如此下去，就可以将T 分割为l 2块。

定理2.1说明，当F E 充分小时，T ~

的特征值很接近T 的特征值，但是一般

情况下F E 是不可能充分小的。

3.2.2 胶合

假设分割后所得的子秩阵块()0T 和()1T 的特征系统已求得，即设已求得正交

矩阵1Q ,2Q 使得()0T 和()1T 分别为对角阵21D D ，,即

()()T

T Q T Q D Q T Q D 2

2221011==，, (3.2.2.1) 记???

?

??=21

Q Q Q ,则有 T T T T zz D D zz Q T Q QTQ ρρ+???

?

?

?=+=21~

, (3.2.2.2) 其中Qv z =。

由于()0T 和()1T 都是不可约对称三对角矩阵，它们的特征值为两组互不相同

的实数，则有T ~的特征值至多为二重的，并且当T ~

存在二重特征值时，有下面的

定理成立。

定理3.3 设d （d 是()0T 和()1T 的特征值）为T ~

（如（3.2.1.2）式所示）的二

重特征值，则d 是T 的特征值。

证明：设d 为()0T 的特征值时，对应特征向量为()01≠x ，d 为()1T 的特征值时，对应特征向量()02≠x ，即 ()()221110,dx x T dx x T ==，对任意标量0011≠≠v v 和，使

得向量()

()()()

102211,~,T T d i a g T x v x v x T

T T ==为的特征值d 对应特征向量，令

()21,z z z =，则

()

()

z x z v x z v dx x zz T Tx T

T T 222111~++=+=ρρ， （3.2.2.3） 选取标量1v 和2v 使得0222111=+x z v x z v T T ，即（3.2.2.3）式变为：

()

dx x zz T Tx T =+=ρ~

，这说明d 也是T 的特征值，定理得证。

根据文[1,10]的定理和定理3.3我们可以得到下面的定理。

定理 3.4 设()()N d d d diag D D diag D ,,,,2121 ==,即T ~

的特征值为

N d d d ≤≤≤ 21，设T 的特征值为N λλλ<<< 21,则

（1）当N N T d d z z d ≤≤≤≤≤≤-<λλλρρ 21110时有，; （2）当z z d d d d T N N N ρλλρ+≤≤≤≤≤≤< 2110时有，； 对

某

个

()

N j j ≤≤1λ,其对应的特征向量为

()z I D a u j j j 1--=λ()()

z I D a u j j j +-=λ或，这里()+

-I D j λ为矩阵

)I D j λ-（的Moore-penrose 广义逆。

不妨设0>ρ，此时令z z d d T N N ρ+=+1。有定理

3.4知,

[]()N j d d j j j ,,2,1,1 =∈+λ,即T ~

的特征值为T 的特征值提供了很好的包含空间。

我们可以根据定理3.3来判断区间端点是否为T 的特征值，从而包含区间可以简化为()()j j j j j d d d d ≠∈++11,λ，如果j j d d =+1则j j d =λ。令2/1）（j j j d d x +=+，通过调用算法Modify-DetEval 计算在特征多项式在j x 处的符号变号数，以此判断

()j j j x d ,∈λ或()1,+∈j j j d x λ确定j λ的包含空间，为了从j λ的包含区间内提取出j λ，我们选取初始点j x 用Laguerre 迭代从j λ的包含区间内提取出j λ。首

先判断在j x 处是否满足调用Laguerre 迭代标准，如果不满足，继续用二分法精化j λ的包含区间并重复上述过程。对0<ρ时有类似上述的处理过程。

4 数值实验 4.1 数值算例

我们已经在本文第三章和第四张中对求反对称三对角矩阵特征值的方法进行了讨论，为了对上述算法进行可行性检验，我们用C 语言编程，在 INTER PENTINUM(CPU 为1.5GHz,内存为240M ，机器精度为16102.2-?=esp )的计算机上用QL 算法（调用LAPACK 中的子程序dhseqr ）对下述两类矩阵进行数值试验（只给出对称三对角矩阵的次对角元）

类型一 ()()1,....,2,1,-=-=n k k n k k β，具体特征值为：(){}n

i k n 112-+-

类型二 ()1,....,2,1),(-=∈=n k r b b k β，具体特征为：n

i n k b 11cos 2????

??

+π其中

1=b

4.2精度检验

我们采用下述方法进行精度测试：首先我们知道类型一，类型二的特征值是

已知的，从而我们可以直接计算对应特征值的绝对误差，令~j λ为实际特征值~

j λ的近似值（由上述算法计算所得）表4.1（其中n 表示矩阵的阶数，ε表示机器的精度）给出了误差

?

?????-T j j j

λλ~

max

由表4.1 可以看出，在上述四方法中，这些方法的精度具有同样的数量级。相比而言，S_SM(L)算法精度最高，而S_QL算法和S_DL(L)算法的精度较低，这也符合前面的分析：对于很小的特征值，开方后可能失去意义：QL算法和S_QL 算法的精度随着矩阵阶数的增大而逐渐降级，而S_SM(L)算法似乎不受矩阵的阶数增大的影响；对类型一的矩阵，当阶数n=2000时，QL算法的误差是S_SM(L)算法的误差的735倍，S_QL算法的误差是S_SM(L)算法的误差的1197倍。