第 9章 矩阵特征值问题的数值方法
9.1 特征值与特征向量
9.2 Hermite矩阵特征值问题
9.3 Jacobi方法
9.4 对分法
9.5 乘幂法
9.6 反幂法
9.7 QR方法引言工程实践中有多种 振动问题,如桥梁或建筑物的振动,机械机件、飞机机翼的振动,工程实践中有多种振动问题,如桥梁或建筑物的振动,机械机件、飞机机翼的振动,及一些 稳定性分析 和相关分析可转化为求矩阵特征值与特征向量的问题。
London,England,Millennium ('Wobbly') Bridge
(1998-2002,Norman Foster and Partners and Arup
Associates)
I decide that I have to write something today,otherwise I
would not know how to speak English here,
This is a very quick story about a bridge,
London launched three major construction projects to
celebrate the arrival of the Millennium,After all,
Greenwich (pronounced green-ich) is supposed to be
(supposed to be?!) where the prime meridian lies,and
the place where the Millennium officially starts in the
world,The three projects are the Millennium Dome in
North Greenwich,so far the largest single roofed
structure in the world,London Eye right across
Westminster,which becomes so far the largest
observation wheel in the world,and the Millennium
Bridge that links Southeast London with St,Paul’s
Cathedral,which is currently…well...not swinging any
more,it is said,
The bridge was designed by Imperial College,a college of my
former university,On the very first day that the bridge was open to
public,there were simply so many people going there to walk from
the south bank to St,Paul’s that the weight completely exceeded
the architect’s expectation,
The slender steel truss bridge began to vibrate with a million
people on there,The opening ceremony ended up in an
embarrassing vertigo,
Millennium left Londoners a happy adage about swinging bridge,
meaning fancy technology that looks good but functions in a funny
fashion,
Am I using too many F’s here? Or is it simply because my tongue
starts to swing in the same direction when I am writing about this
wobbly bridge?
Next time you visit London,I strongly recommend this place,
After all,with a little swing,this is a shortcut to dash into St,
Paul’s directly from the southeast!
xxG?T
1ex?T
G,Google Matrix,
“the world’s largest matrix computation”,
4,300,000,000
x,PageRank(网页级别) vector
“The $25,000,000,000 Eigenvector”
搜索引擎
9.1 特征值与特征向量设 A是 n阶矩阵,x是非零列向量,如果有数 λ存在,满足,(1)
那么,称 x是矩阵 A关于特征值 λ的特征向量,
如果把 (1)式右端写为,那么 (1)式又可写为,
Ix?
( ) 0I A x
| | 0IA即
1
1 1 0( ) | |,,,
nn
nf I A a a a

记它是关于参数 λ的 n次多项式,称为矩阵 A的特征多项式,其中 a0=(-1)n| A|,
(2)
显然,当 λ是 A的一个特征值时,它必然是 的根,反之,如果 λ是 的根,
那么齐次方程组 (2)有非零解向量 x,使 (1)式成立,从而,λ是 A的一个特征值,
A的特征值也称为 A的特征根,
( ) 0f( ) 0f
矩阵特征值和特征向量有如下主要性质:
定理 9.1.1 n阶矩阵 A是降秩矩阵的充分必要条件是 A有零特征值,
定理 9.1.2 设矩阵 A与矩阵 B相似,那么它们有相同的特征值,
定理 9.1.3 n阶矩阵 A与 AT有相同的特征值,
定理 9.1.4 设 λ i≠λ j是 n阶矩阵 A的两个互异特征值,x,y分别是其相应的右特征向量和左特征向量,那么,xTy=0,
9.2 Hermite矩阵特征值问题设 A为 n阶矩阵,其共轭转臵矩阵记为 AH,如果 A=AH,那么,A称为 Hermite矩阵,
9.2.1 Hermite矩阵的有关性质设 是 Hermite矩阵 A的 n个特征值,有以下性质,
全是实数,
12,,.,,,n
12,,.,,,n
有相应的 n个线性无关的特征向量,它们可以化为一组标准酉交的特征向量组,即
12,,.,,,n
12,,.,,,nu u u Hiju u ij
是酉空间中的一组标准酉交基,
12,,.,,,nu u u
记 U=( ),它是一个酉阵,即
UHU=UUH=I,那么即 A与以 为对角元的对角阵相似,
12,,.,,,nu u u
1
H
n
U A U D




12,,.,,,n
A为正定矩阵的充分必要条件是全为正数,
12,,.,,,n
定理 9.2.1 设 是 Hermite矩阵 A的 n
个特征值,那么证,
12,,.,,,n
2 1m a x iinA
2
1
n
iF
i
A?

2 22
2 ( ) ( ) ( ( ) )
HA A A A A由
2 1m a x iinA因 此
2 22
1
( ) ( )
n
H
iF
i
A tr A A tr A?
又 由
2
1
n
iF
i
A?
得设 x是一个非零向量,A是 Hermite矩阵,
称 为矩阵 A关于向量 x的 Rayleigh商,
记为 R(x),
H
H
x Ax
xx
定理 9.2.2 如果 A的 n个特征值为其相应的标准酉交的特征向量为那么有
12,.,n
12,,.,,,nu u u
1 () nRx
定理 9.2.3 设 A是 Hermite矩阵,那么
100
m i n ( ) m i n ( )
k n k
kkx C x x C xR x R x


且 且或
9.2.2 极值定理定理 9.2.4(极值定理 ) 设 Hermite矩阵的 n个特征值为,其相应的标准酉交特征向量为,用 Ck表示酉空间 Cn中任意的 k维子空间,那么
12,.,n
12,,.,,,nu u u
1 1
0
0
m a x m in ( )
m in m a x ( )
kk
nk nk
k
x C xC
k
C x C x
Rx
Rx



且且或
9.2.3 Hermite矩阵特征值问题的性态矩阵特征值问题与求解线性方程组问题一样,都存在当矩阵 A的原始 数据有小变化 (小扰动 )时,引起特征值问题的变化有大有小的问题,
如果引起的变化小,称 该特征值问题是 良态的,
反之,称为 病态的,
矩阵特征值问题的性态是很复杂的,通常分别就单个特征值或整体特征值给出条件数进行分 析,对于 Hermite矩阵,由于其特征值问题的特殊性质,其特征值都是良态的,下面先证明
Hermite矩阵特征值的扰动定理,
定理 9.2.5 设矩阵 A,E,A+E都是 n阶 Hermite
矩阵,其特征值分别为那么,
证 设矩阵 A关于特征值 λ1,λ2,…,λn 的标准酉交特征向量为 u1,u2,…,un,
是由 ui,ui+1,…,un生成的 n-i+1维子空间,
对 中任意非零向量 x,由极值定理,有
12 n 12 n
12 n 1i n i i
1niC
1niC
1
11
0
00
()
m a x
m a x m a x
ni
n i n i
H
i H
x C x
HH
HH
x C x x C x
x A E x
xx
x A x x E x
x x x x





且且 且由定理 9.2.3,
又由定理 9.2.2,对任意 x≠0,有从而有另一方面,A=(A+E)-E,记 为矩阵 -
E的特征值,那么,
重复上面的过程,可得从而有
1 0
m a x
ni
H
iHx C x
x A x
xx


1
1 0m a x
ni
H
nHx C x
x E x
xx




1ii
12 n
1i n i
1ii
i i n
定理 9.2.5通常又称为 Hermite矩阵特征值的 扰动定理定理 9.2.6 设矩阵 A和 A′=A+E都是 n阶 Hermite矩阵,其特征值分别为和,那么 12 n

12 n 222iiEE
这个定理表明,扰动矩阵 E使 A的特征值的变化不会超过 ‖ E‖ 2,一般 ‖ E‖ 2小,因此,
Hermite矩阵特征值是良态的,
9.3 Jacobi方法
理论上,实对称矩阵 A正交相似于以 A的特征值为对角元 的 对角阵,问题是如何构造这样的正交矩阵呢? Jacobi方法就是通过构造特殊的正交矩阵 序列,通过相似变换使 A的非对角线元素逐次零化来实现对角化的,
9.3.1 平面旋转矩阵与相似约化先看一个简单的例子,
设 是二阶实对称矩阵,即 a21=a12,
其特征值为 λ1,λ2,令使得 记容易验证 BT=B,且
1 1 1 2
2 1 2 2
aaA
aa


c o s s in
s in c o sR





1
1
TR A R?

1 1 1 2
2 1 2 2
T bbB R A R
bb


22
11 11 12 22
22
12 21 22 11 12
22
22 11 12 22
c os 2 sin c os sin
( ) sin c os ( c os sin )
sin 2 sin c os c os
b a a a
b b a a a
b a a a






解之得,
当 时当 时可选取
1 1 2 2aa? 12
11 22
21 a r c ta n
2
a
aa
1 1 2 2aa?
4

为使 RTAR为对角阵,要求 b12=b21=0
一般的 n阶平面旋转矩阵
21
a r c ta n
2
pq
pp qq
a
aa

4

9.3.2 经典的 Jacobi方法设 A是实对称矩阵,记 A1=A.Jacobi方法的基本思想是用迭代格式
Ak+1=QTkAkQk,k=1,2,…
构造一个相似矩阵序列,使{ Ak}收敛于一个对角阵,其中 Qk为平面旋转矩阵,
其旋转角 θk由使 Ak的绝对值 最大元
a(k)pq=a(k)qp=0 或按列依次使 A的非对角元零化来确定,
定理 9.3.1 设 A是 n阶实对称矩阵,那么由 Jacobi方法产生的相似矩阵序列{ Ak}
的非对角元收敛于 0,也就是说,{ Ak}收敛于以 A的特征值为对角元的对角阵,
记 其中 Ek是 Ak除主对角元外的矩阵,由平面旋转矩阵的性质中,对于,有因此,
()() kk ii kA d ia g a E
1
T
k p q k p qA R A R,i p q?
( 1 ) 2 ( 1 ) 2 ( ) 2 ( ) 2( ) ( ) ( ) ( )k k k ki p i q i p i qa a a a
22 ( ) 2
1 2 ( )
k
k k p qFFE E a
又由假设,
因此,
这样,便有从而,当
22( ) ( ) ( ) ( )
1,,1
m a x ( 1 )
n
k k k k
p q i j p q p qi j n i j
i j i j
a a a n n a


且 且且
22 ()( 1 ) k
k p qFE n n a
2 2 2
11 22
22( 1 ) ( 1 ) k
kkF F FE E En n n n
1,0k FkE
9.3.3 实用的 Jacobi方法
循环 Jacobi方法必须一次又一次扫描,才能使{ Ak}
收敛于对角阵,计算量很大,在实际计算中,往往用一些特殊方法来控制扫描次数,减少计算量,
下面介 绍一种应用最为广泛的特殊循环 Jacobi方法 ——阈 Jacobi方法,阈 Jacobi方法首先确定一个阈值 δ,在对非对角元零化的一次扫描中,只对其中绝对值 超过阈值的非对角元进行零化,当所有非对角元素的绝对值都不超过阈值后,将阈值减少,
再重复下一轮扫描,直至阈值充分小为止,减少阈值的方法通常是先固定一个正整数 M≥n,扫描一次后,让,而阈值的下界是根据实际问题的精度要求选定的,M

9.3.4 用 Jacobi方法计算特征向量
假定经过 k次迭代得到 Ak+1=RTk…R T1AR1…R k,(15)
这时 Ak+1是满足精度要求的一个近似的对角阵,如果记 Qk=R1R2…R k=Qk-1Rk,(16)
那么,Qk是一个正交矩阵,且 (15)式又可表示为 Ak+1=QTkAQk.当 Ak+1的非对角元素充分小,Qk的第 j列 qj可以看成是近似特征值
a(k+1)jj相应的特征向量了,
在实际计算中,可以按 (16)式在迭代过程中形成 Qk,把 Qk看成是 Qk-1右乘一个平面旋转矩阵得到,不妨记 Q0=I,Qk的元素按下式计算,( ) ( 1 ) ( 1 )
( ) ( 1 ) ( 1 )
( ) ( 1 )
c o s sin,
sin c o s,
,,,
1,2,
k k k
ip ip k ip k
k k k
iq ip k iq k
kk
ij ij
q q q
q q q
q q j p q
in







9.4 对分法理论上,一个实对称矩阵正交相似于一个以其特征值为对角元的 对角阵,但是,
经典的结果告诉我们,一个大于 4次的多项式方程不可能用有限次四则运算 求根,
因此,我们不可能期望只用有限次相似变换将一个实对称矩阵约化为一个对角阵,下面先介绍将一个实对称矩阵相似约化为实对称三对角矩阵的方法,再讨论求其特征值的对 分法,
9.4.1 相似约化为实对称三对角矩阵将一个实对称矩阵正交相似约化为一个实对称三对角矩阵的算法,可归纳如下:
记 A(1)=A,对 k=1,2,…,n -2
① 按 (4)式,(5)式和 (8)式计算 ;
②按 (9)~ (12)式,计算 A(k+1),
,k k ku 和
()
1,
()
2,
()
( 4)
k
k k k
k
kk
k
k
nk
a
a
u
a








( ) ( ) 2
1,
1
( ) ( ) ( 5 )
n
kk
k k k i k
ik
s ig n a a


() 1,( ) ( 8 )kk k k k ka
1 ( )
1
( 1 ) ( )
,( 9)
,( 10 )
1
,( 11 )
2
( ),( 12 )
k
k k k
T
k k k k
k k k k
k k T T
k k k k
g A u
ug
y g u
A A u y y u



9.4.2 Sturm序列的性质
设实对称三对角矩阵为其中 βi≠0 (i=1,2,…,n-1)
其特征矩阵为 T-λI,记 T-λI的第 i阶主子式为
11
1 2 2
1nn
a
Ta
a









11
12
1
1
()
i
i
ii
a
a
p
a



这是关于 λ的 i次多项式,当 i=n时,pn(λ)=|
T-λI|是矩阵 T的特征多项式,令 p0(λ)≡1,
则有 p1(λ)=α1-λ,pi(λ)=(αi-λ)pi-1(λ)-β2i-1pi-2(λ),
i=2,3,…,n.(15)
多项式序列{ pi(λ)} (i=0,1,…,n)称为
Sturm序列
11
12
1
1
()
i
i
ii
a
a
p
a



定理 9.4.1{ pi(λ)} (i=1,2,…,n)的根都是 实根,
证 由 (14)式,pi(λ)是 i阶实对称矩阵的特征多项式,因此,{ pi(λ)} (i=1,2,…,n)的根全是实根,
定理 9.4.2
定理 9.4.2 设 α是 pi(λ)的一个根,那么
① pi-1(α )pi+1(α )≠0,即相邻的两个多项式无公共根;
② pi-1(α )pi+1(α )<0,即 pi-1(α )与 pi+1( α )反号,
l im ( ) 0,ip
l i m ( ) 0,l i m ( ) 0,iipp当 i 为 奇 数 当 i 为 偶 数定理 9.4.4 pi(λ)的根都是单根,并且将 pi+1(λ)的根严格隔离,
9.4.3 同号数和它的应用
定义 1 设 p0(λ)≡1,{ pi(λ)} (i=1,2,…,n)
是一个 Sturm序列,称相邻的两个数中符号一致的数目为 同号数,记为 ai(λ),若某个 pi(λ)=0,规定与 pi-1(λ)反号,
定理 9.4.5 设两个实数 x<y,那么,形如
(13)式的实对称 三对角矩阵 T的特征多项式在区间 (x,y]上根的数目为 a(x)-a(y),
9.4.4 求 Hermite矩阵特征值的对分法
对分法的计算可归纳为以下 4个部分
①确定 (13)式的矩阵 T的全部特征值的分布区间,
② 在区间[ a,b]中,用区间对分的方法找出只含 T的一个特征值的子区间,
③ 在只含一个特征值的子区间上的对分法,
④ 同号数的计算,
9.5 乘幂法
设 A是 n阶矩阵,其 n个特征值按模从大到小排序为又假设关于 λ1,λ2,…,λn的特征向量 v1,
v2,…,vn线性无关,
1 2 3 n
2
0 1 1 1 2 2
11
[ ( ) ( ) ]k k k kn nA x a v a v v


乘幂法是适用于求一般矩阵按模最大特征值及相应特征向量的算法,
xk→λ k1a1v1 (k→∞).
因此,xk可看成是关于特征值 λ1的近似特征向量,
迭代格式为
1,0,1,
k
k
k
kk
x
z
x
x A z k

按模最大特征值 λ1及其相应的特征向量 v1的乘幂法的计算公式:
1
1
1
,0,1,
k
k
k
kk
T
k kk
T
kk
x
z
x
x Az
z Az
k
zz

9.5.2 收缩方法
设矩阵 A的 n个特征值按模从大到小排序为,其相应的 n个线性无关特征向量为 v1,v 2,…,vn,在计算 A
的最大特征值 λ1及相应特征向量 v1后,可以通过收缩方法,继续用乘幂法计算 λ2及其相应的特征向量 v2,
1 2 3 n
定义 n阶矩阵
把去掉 A1的第 1行和第 1列的 n-1阶矩阵记为
2 1 2 1 1 1 2 2 2 1 1 2 2 2 1 1
11
1 1 1 1 2 1 1 2 1 1
0 0 0
nnT
n n n n n n n n
a v a a v a a v a
A A v
a v a a v a a v a






2 2 2 1 1 2 2 3 2 1 1 3 2 2 1 1
3 2 3 1 1 2 3 3 3 1 1 3 3 3 1 1
2 1 1 2 3 1 1 3 1 1
nn
nn
n n n n n n n n
a v a a v a a v a
a v a a v a a v a
B
a v a a v a a v a





那么,B有与 A1除 λ1外的相同的 n-1个特征值 | λ2| >| λ3| ≥…≥ | λn|,可以用乘幂法计算 λ2及其相应的 特征向量,在计算 λ1和 v
1 后,按 (15)式形成 n-1阶矩阵 B的计算过程称为 收缩方法,
9.6 反幂法
反幂法可以求一个非奇异矩阵 A的逆矩阵 A-1
的按模最小的特征值及相应的特征向量,又可以求 A的一个近似特征值相应的特征向量,
9.6.1 求按模最小特征值及相应特征向量的反幂法,又称为 反迭代法,
1
( 1 ) 1
,0,1,
k
k
k
kk
T
k kk
T
kk
x
z
x
L U x z
zx
k
zz


9.6.2 求近似特征值的特征向量的反幂法
先对矩阵 进行 LU分解,记那么,
(7)
下面介绍一种选取特殊的初始向量 x0的反幂法 ——半迭代法,
l IA l
I A L U
1
,0,1,
k
k
k
kk
kk
x
z
x
L y z
U x y k

假设,选取初始向量 x0满足 ‖ x0‖ ∞=1,
这时 z0=x0.对照 (7)式中的第二个式子,可把 z0
看成满足 Le=z0.(8)
这里,e=(1,1,…,1)T,而 z0的各个分量的取值多少是无关重要的,这样,在第一个迭代步的计算中,只需求解 (7)式中的上三角方程组 Ux1=e.,半迭代法”的命名也由此而得,
lIA
9.7 QR方法定理 9.7.1 设 A是 n阶矩阵,其 n个 特征值为,那么存在一个酉矩阵 U,使
UH AU是以为 对角元的上三角矩阵,
12,,,n
12,,,n
9.7.1 两个基本定理定理 9.7.2 设 A是 n阶实矩 阵,那么,存在一个正交矩阵 Q,使 QT AQ为一个准上三角矩阵,
它的对角元是 A的一个特征值,对角元上的二阶块矩阵的两个特征值是 A的一对共轭复特征值,
9.7.2 相似约化为上 Hessenberg矩阵对一般 n阶矩阵,QR算法的每一个迭代步需要
O(n3 )次乘法运算,如果矩 阵阶数稍大,这个算法几乎没有实际的应用价值,
通常采用的方法是先将矩阵相似约化为上
Hessenberg形式的矩阵,在此基础上应用 QR
迭代,这时,一个 QR迭代步的乘法运算次数只需 O(n2 )次,
所谓上 Hessenberg矩阵是指一个 n阶矩阵 A,如果当 i>j+1时,aij=0,称 A为上 Hessenberg矩阵,例如:一个 5阶的上 Hessenberg矩阵具有如下的形式:
* * * * *
* * * * *
0 * * * *
0 0 * * *
0 0 0 * *
A








下面介绍 QR方法时,都假设矩阵 A是一个上 Hessenberg矩阵,