第五章 求矩阵特征值和特征向量

n阶方阵A的n个特征值就是其特征方程

det(AI)0

的n个根，方程A属于特征值的特征向量x是线性方程组

Axx

的非零解。本章讨论求方阵A的特征值和特征向量的两个常用的数值方法。以及求实对称矩阵特征值的对分法。

5.1 幂 法

在实际问题中，矩阵的按模最大特征根起着重要的作用。例如矩阵的谱半径即矩阵的按模最大特征根的值，它决定了迭代矩阵是否收敛。本节先讨论求实方阵的按模最大特征根的常用迭代法：幂法。

5.1.1幂法的基本思想

幂法是求实方阵A按模最大特征值及其特征向量的一种迭代方法。它的基本思想是：先任取非零

初始向量x0，然后作迭代序列

xk1Axk，k0,1, （5。1）

再根据k增大时，xk各分量的变化规律：按模最大的特征向量会愈来愈突出，从而可求出方阵A的按模最大特征值及其特征向量。

先看一个计算实例。 例1 设矩阵

1A2用特征方程容易求得A的两个特征值为

2 111，23

下面用幂法来计算，取初始向量x01,0，计算向量序列 xk1Axk，k0,1, 具体结果如表5.1所示.

表5.1 幂法计算结果

Tk

0

x1(k)

1

(k) x20

1 2 3 4 5 6 7

考察两个相邻向量对应分量之比：

1 5 13 41 121 365 1093

2 4 14 40 122 3 1094

x1(2)x1(3)x1(4)x1(5)x1(6)x1(7)5,(2)2.6,(3)3.1,(4)2.951,(5)3.016，(6)2.994 (1)x1x1x1x1x1x1(2)(3)(4)(5)(6)(7)x2x2x2x2x2x22,(2)3.5,(3)2.857,(4)3.05,(5)2.983, (6)3.005 (1)x2x2x2x2x2x2由上面计算看出，两相邻向量对应分量之比值，随k的增大而趋向于一个固定值3，而且这个值恰

好就是矩阵A的按模最大的特征值。这一现象是否有普通性？下面进行具体分析。

5.1.2 幂法的计算公式

为简便起见，设矩阵A的几个特征值按模的大小排列如下：

12n

其相应特征向量为u1,u2,un，并且是线性无关的，因此可作为n维向量空间的一组基。

任取初始向量x0x1,x2,(0)(0)(0)xn0，首先将x0表示为

Tx0a1u1a2u2作迭代序列

anun

xk1Axk, k0,1,

则

x1Ax0a11u1a22u2annun

…… ……

kkxkAxk1a11ku1a22u2annun

于是

kkkn2xk1a1u1a2u2anun 11为了得出计算1和u1的公式，下面分三种情况讨论。 1.1为实根，且

12

2

当a10，k充分大时，则有

xka11ku1 xk1a11k1u11xk

所以

xi(k1) 1(k) ， u1xk （5.2）

xi2. 1为实根，且12，

23

当a1,a2不为0，k充分大时，则有

xkka1k1u1(1)ka21u2 xk111k1a11u1(1)ka2k1u2

xk2k2k2a11u1(1)k2a21u221xk于是得

x2k2A2xk1xk

(A221I)xk0

所以

1(k2)1xi2x(k),u1xk11xki 1x(k2)i22x(k),u2xk12xki3．1uiv,2uiv,23

当k充分大时，则有

xkka11uk1a22u2

xk1k1k1a11u1a22u2 xk2k2a21k1u1a22u2

于是得

3

5.3）

（xk2(12)xk112xk[若令

k21(12)k11k112]a1u1[k22(12)k12k112]a2u20

12p,12q

得

xk2pxk1qxk0 （5.4）

式（5.4）是以p,q为变量，以xk2,xk1,xk的几个分量为系数的矛盾方程组。

用最小二乘法解矛盾方程组（5.4），求出p,q，然后再解一元二次方程

x2pxq0

得到的两个根便是1、再由

2的近似值。

[A2(12)A12]xk0

可得

(A1I)(A2I)xk(A1I)(xk12xk)0

和

A2ΙA1ΙxkA2Ιxk11xk0

综上所述，可得

1p1P24q 22（5.5）

u1xk12xk

2p122p24q

u2xk11xk

在实际应用幂法时，可根据迭代向量各分量的变化情况来判定属于哪种情况。

若迭代向各分量单调变化，且有关系式xk1cxk，则属于第1种情况;若迭代向量分量变化不单调，但有关系式xk2cxk，则属于第2种情况；若迭代向量各分量变化不规则，但有关系式

xk2pxk1qxk0，则属于第3种情况。

5.1.3 幂法的实际计算公式

4

当k时，若

11，则x的分量会趋于无穷大；若11，则x的分量会趋于零。因此会使

计算机出现上溢或下溢现象。为了防止溢出，可采用如下迭代公式

mkmax(xk)xk yk

xkxkmkxk k0,1,

（5.6）

xk1Αyk

y0(1,1,n,1)T

式（5.6）的更详细的带计算过程的计算公式为：

xi(k)aijy(jk1) i1,2,j1,n

（5.7）

yi(k)xi(k)mk mk||xin(k)||maxxi(k)x(pk) k1,2,i

s(myki1(k)ixi(k))

y0(1,1,,1)T

(p)(p)注：当xk0时，按模最大特征值为正，故计算s时取xi(k),当xk0时，取xi(k)。

由式（5.6）知

xk1Ayk2xk2A2yk3A2xk3Axk2xk3xk2xk3… …

在式（5.8）两端同乘以A，得

xk2Ak1x0xk3x0 （5.8）

 Axk1因为

xk2Akx0xk3x0 （5.9）

xkAyk1Axk1xk1Axk1xk1

将式（5.8）、（5.9）两端分别取范数后，代入上式得

5

xkkkkn21a1u1a2u2anun11所以

k112na1u1a2u2anun11kk

当k时，mkxk因此当k充分大时，mk就是按模最大的特征值

利用式（5.9）可得

1 （5.10）

1的近似值。

ykxkAyk1Axk1xkxkxkxk1

另一方面，有

Akx0xkxk1x0 （5.11）

Akx0Aky0x0k1Aky0x0x0

AAy0x0Ak1x1

Ak1y1x1x0 （5.12）

… …

xk将式（5.12）代入式（5.11），有

xk1x0k

ykAkx0Ax0k2k1a1u1a21k1a1u1an21u2u2kknanun1knanun1

从而

当k时，yku1 （5.13） u1这说明归一化向量序列yk收敛于按模最大的特征值所对应的特征向量。因此，当k充分大时，yk就是特征向量u1的近似值。

6

5.1.4 幂法的计算步骤

为节省篇幅，这里仅介绍1231．输入

矩阵的阶n，系数aij(i1,n。

,n;j1,,n)，初始向量y0=(1,1,,1)T，允许误差ep

2．输出

特征值和特征向量yi（i1,2,,n）

3．计算步骤

1)给出迭代先导条件 计算公式 s2ep

对应语句 s=2*ep;

2)用幂法求按模最大特征值及特征向量 计算公式 式(5.7) 对应语句

while(s>ep)

{ for(j=1;j<=n;j++)

x[i]=x[i]+a[i][j]*y[j]; m=fabs(x[1]); p=1;

for(i=2;i<=n;i++) if(fabs(x[i])>m) { m=fabs(x[i]); P=i; } S=0;

for(i=1;i<=n;i++) s=s+(x[i]-m*y[i]); for(i=1;i<=n;i++) y[i]=x[i]/m;

}

5.1.5 幂法的计算实例

例2 用幂法求矩阵

21A1201的按模最大特征值和相应的特征向量(ep104)。

解 取xT0(1,1,1)，用幂法迭代公式

7

012

 mkmax(xk)xk yk

xkxkmkxk k0,1,

 xk1Αyk 计算结果如表5.2所示。

表5.2 幂法迭代公式计算结果表 k 0 1 2 3 4 5 2.416 7 -3.416 7 2.416 7 3.416 7 0.707 3 -1.000 0 0.707 3 6 2.414 6 -3.414 6 2.414 6 3.414 6 0.707 1 -1.000 0 0.707 1 1.000 0 2.000 0 3.000 0 2.500 0 2.428 6 xi(k) 1.000 0 -2.000 0 -4.000 0 -3.500 0 -3.128 6 1.000 0 2.000 0 3.000 0 2.500 0 -2.428 6 mk 1.000 0 2.000 0 4.000 0 3.500 0 3.428 570 1.000 0 1.000 0 0.750 0 0.714 3 0.708 3 -1.000 0 yik 0.000 0 -1.000 0 -1.000 0 -1.000 0 1.000 0 1.000 0 0.750 0 0.714 3 0.708 3 T所以1m63.414634,u1y6(0.707143,。事实上，矩阵A的最大特征值1,0.707143)为122，其对应的特征向量u1(22T,1,)。 22 5.2 逆幂法

逆幂法的功能是求A按模最小特征值和和特征向量。

设A为非奇异方阵，12n为A的特征值，u1,,un为其相应的特征向量，则A的特征值为

111121n其相应的特征向量仍为u1,,un。A按模最大特征值的倒数则为阵A按

1模最小特征值。

5.2.1 逆幂法的计算公式

取x0，计算向量序列

xk1A1xk,k0,1, （5.14）

它等价于

Axk1xk,k0,1, （5.15）

8

这样一来，我们可以通过反迭代过程，即通过解方程组（5.14）求得xk1。

当

n1n,an0,k充分大时，则有

xi(k)n(k1),unxk1 (5.15)

xi在实际计算中，为了减少运算量，可先将A作三角分解

ALU

然后再解两个三角形方程组就可以了，再考虑到防止数据溢出，即得逆幂法的实际迭代步骤如下： （1）对A作LU分解 （2）解方程组

LUx(k1)y(k)

即解

Lz(k1)y(k)(k1)z(k1)Uxx(k)(k)y||x(k)||(0)T(0)y(1,1,,1)或y(1,0,令

,0)Tmaxxi(k)x(pk)imkmaxxi(k)x(pk)ix(pk)0x(k)p0

min迭代条件为

1 mk||x(k)mky(k)||

5.2.2 逆幂法的计算步骤

逆幂法的迭代过程与幂法一样，前面已经介绍过LU分解法，这里不再详细叙述逆幂法的计算步

骤。

例1 用逆幂法求矩阵

210A121

012按模最小特征值及相应特征向量。（精度0.05）

解 先对A作三角分解ALU

9

11L20取x0(1,0,0)T,用逆幂法迭代公式

00102132103U01

24003yk得计算结果如表5.3所示。

xkxk；Lzkyk；Uxk1zk k0,1,

表 5.3 逆幂法迭代公式计算结果 k 0 1.000 000 1 0.250 000 0.500 000 0.750 000 0.750 000 0.333 333 0.666 667 1.000 000 0.333 333 0.833 333 1.555 551 2 0.500 000 1.333 333 1.166 667 1.333 333 0.375 000 1.000 000 0.875 000 0.375 000 1.875 000 1.666 667 3 0.937 500 1.500 000 1.250 000 1.500 000 0.625 000 1.000 000 0.833 333 0.625 000 1.312 500 1.708 333 4 1.177 083 1.729 167 1.281 250 1.729 160 0.680 722 1.000 000 0.749 0 0.680 722 1.340 361 1.634 538 5 1.180 722 1.710 843 1.225 904 1.710 843 xi(k) 0.000 000 0.000 000 mk 1.000 000 1.000 000 yi(k) 0.000 000 0.000 000 1.000 000 zi(k) 0.500 000 0.333 333 从上表可看出，

min10.5845作为按模最小的特征值的近似值，m5x(5)(1.1807,1.7108,1.2590)T作为相应的特征向量的近似值，而A的按模最小的特征值的理论值为220.5859.

~5.2.4 用逆幂法求在附近的特征值的计算公式

设与最接近的特征值为i，即有

~|i||j|,ji,j1,作矩阵AI，它的特征值及相应特征向量为

,n

jj和uj,j1,,n

10

若用逆幂法求AI，则有

(AI)xk1xk,k0,1,

则可求得AI按模最小特征值和相应特征向量为

xi(k)ii(k1) xiuxk1i于是得矩阵A在附近特征值和相应特征向量为

~xi(k)ii(k1) xiuxk1i

~5.2.5 用逆幂法求在附近的特征值的计算实例

例2 用逆幂法求矩阵

220A021

012在2.93附近的特征值及相应特征向量(精度ep105)。

解 对矩阵A2.93I作三角分解

100.93A2.93I00.931010.93100取x0(0,0,1),由迭代公式

T

000.93101000.931111000.930.930.93xkykxk得计算结果如表5.4所示。

；Lzkyk；Uxk1zk k0,1,

表5.4逆幂法的计算

11

k 0 0.000 00 1 7.959 06 6.883 79 7.959 06 1.000 00 -0.930 00 0.8 90 1.000 00 -0.930 00 1.8 90 2 12.692 31 -12.803 85 12.837 58 12.837 58 0.988 68 -0.997 37 1.000 00 0.988 68 -0.997 37 2.072 44 3 14.278 43 -14.267 63 14.266 27 14.266 27 1.000 00 -0.999 24 0.999 15 1.000 00 -0.999 24 2.073 60 1m3

xi(k) 1.0000 00 -7.401 92 1.000 00 mk 1.0000 00 0.000 00 yi(k) 0.000 00 1.000 00 0.000 00 z(k)i 0.000 00 1.000 00 1由表5.4知，(A2.93I)的按模最大的特征值为m3，即A2.93I的按模最小的特征值为

，

所以矩阵A的特征值为

1m3

+2.93≈3.000 36，相应的特征向量为(1,-0.999 24,0.999 15)

T

5.3实对称阵特征值的对分法

首先讨论三对角对称矩阵的情形，再讨论一般对称矩阵的情形。

5.3.1 求实三对角阵特征值的对分法

1．实对称三对角的Sturm序列 设实对称三对角阵

c1b1Cb1c2b2b2c3b3bn2cn1bn1 bn1cn其中bi0(i1,,n)，用p()表示CI的i阶主子式，并规定p0()1,b00,则

pn()det(CI)为矩阵C的特征多项式，且容易验证

p0()1 p1()c1

p2()(c2)p1()b12p0() (5.17)

… …

12

pi()(ci)pi1()bi21pi2()

称多项式序列{p0(),p1(),,pn()}为矩阵C的Sturm序列。 Sturm序列具有以下性质：

性质1 pi()0(i1,,n)仅有实根。

性质2 相邻两个多项式pi1()和pi()无公共零点。 性质3 设是pi()0的根，则 pi1()pi1()0。

性质4 pi()0(i1,开来。

性质5 设1in1,且是pi()0的根，则当正数足够小时，pi1()和pi()在区间

并且 pi()0的根把pi1()0的根严格地隔离,n)的根全是单根，

~~(,)内同号, pi()和pi1()在区间(,)内同号。

2．序列在某点的连号数

序列Sturm在{p0(),p1(),,pn()}在时的连号数pi()，由以下规则确定

（1） 若pi1()和pi()同号，则从pi1()到pi()有1个连号数；若符号相反，则无连号数。 （2）pi()0，则以pi1()的符号作为它的符号。 （3）p()={p0(),p1(),,pn()}的连号数。 如：p(2)1,0,1,0,1的连号数=2。 3．Gerschgorin定理 定义1 设 A(aij)Rnn~~~~~~~，定义

Di{zzaiiri}ni1,,n （5.18）

为第i个Gerschgorin盘（或称圆盘），其中riaj1jiij为Di的半径。

定理1（Gerschgorin定理或称圆盘定理）设A(aij)Rnn，则A的全部特征值都在区域

DD1D2Dn内。

,n，有zaiiri。所以矩阵AzI是严格对角占优的，而严格

证明 当zD内，即对i1,对角占优矩阵是非奇异的，即det(AzI)0，故z不是A的特征值。换句话说，方阵A的全部特征

13

值都属于D。

由Gerschgorin定理立即可得 推论1 矩阵A的特征值满足

minimin(aiiaij)iin

j1ji （5.19）

nmaxiimax(iaiiaij)jj1i推论2 设C为实对称三对角方阵，令

mmin1jncjbjbj1,Mmax1jncjbjbj1则C的所有特征值都属于区间m,M。

4.求实对称三对角阵的对分法

定理2 矩阵C在区间[,)内特征值的个数等于(~)。 证明略。

例1 求三对角矩阵

2100C1210 01210012在2,0内特征值的个数。

解 首先写出C的Sturm序列

P01， P12，

P22P1P0 P32P2P1 P42P3P2

再分别计算Sturm序列在-2和0两点的连号数

(2)1,0,1,0,1}2(0)1,2,3,4,5}0

所以在2,0上有202个特征值。

14

5.20）

（利用定理2和Gerschgorin定理，由对分法就可将矩阵C的特征值进行隔离，具体步骤如下： 1. 求矩阵C的Sturm序列；

2. 根据Gerschgorin定理求出C的全部特征值的上界M和下界m； 3. 将区间m,M对分，取中点1mM，若1k，则矩阵C在区间1,M有k个特征2值，在区间m,1内有nk个特征值，再分别计算m,1及1,M的中点2，分别计算2，3，

3，继续下去，总可将m,M分成若干个小区间，使得矩阵C在每个小区间上至多有一个特征值，

这样就可以将C的n个特征值隔离开了。

4. 继续对有根区间使用对分法，就可求出满足给定精度的特征值的近似值，这就是求实对称三对角方阵特征值的对分法。

例2 用对方法将三对角方阵

12C002120021200 21的特征值进行隔离，并求出其最大特征值，使它至少有2位有效数字。

解 设C的四个特征值的次序为1234，则im,M，由圆盘定理可得特征值的上界，下界分别为

Mmaxcjbj1bj5

1j4mmincjbj1bj3

1j4即矩阵C的特征值属于区间[3,5]。

C的特征多项式序列为

P01;P11

P1P i2,3,4 ii14Pi2（1）分别计算各点的连号数，具体结果见表5.5

表5.5各点连号数的计算 iPi 0 1 2 3 4 pi(3)1 4 12 32 80 Pi(1) Pi(1) Pi(3) Pi(5) 连号数 4 3 2 1 0 1 2 +0 -8 -16 1 +0 -4 -0 16 1 -2 -0 8 -16 1 -4 12 -32 80 因此，区间3,1,1,1,1,3,3,5各有C的一个特征值，这样就将C的特征值进行了隔离。 （2）最大特征值13,5，计算1列表如下。

15

表5.6 最大特征值所在区间划分

Pi i 0 1 2 3 4 连号数 1所在区间 14,5 Pi4 Pi4.5 Pi4.25 Pi4.22 1 -3 5 -3 -11 1 -3.5 8.25 -14.875 19.062 1 3.25 6.5625 -8.328 0.816 1 -3.22 6.3684 -7.626 -0.917 1 0 0 1 14,4.5 14,4.25 14.22,4.25 由于1区间4.22,4.25,取该区间内任意一点作为最大特征值的近似值都具有2位有效数字。

5.3.2实对称阵的三对角化

为了用对分法求任意实对称阵A的特征值，首先需要正交变换将其化为三对角矩阵，此过程可以通过Househoulder矩阵实现.

n定义2 设vR*Rn\\{0}，n阶Househoulder矩阵HH(v)定义为：

HI其中I为n阶单位阵。

显然，对任意xR，有

n2Tvv （5.21） vTvvTxHxx2Tv

vvnnT因此向量Hx,x以及v共面。特别地，如果xR，且vx0，有Hxx，于是由在R上垂直于

v的所有向量x组成的超平面H在映射xHx下是不变的。最后，对任意xRn

vTHxvTx

因此，如果x和v之间的夹角用表示，则v和Hx的夹角等于。易知：Hx是x关于超平面H的像。由于此原因，映射xHx又称为Househoulder反射，或者称为镜面反射(见图5-1)。

x v

u H Hx

16

图5-1 镜面反射图

引理1 Househoulder阵是对称且正交的。

证明 因为ITI, (vTv)TvTv,(vvT)TvvT且vTv是一个正数，H的对称性显然。又

HTHHHTH2I44Tvv(vvT)(vvT) TT2vv(vv)(vvT)(vvT)v(vTv)vT(vTv)(vvT)

所以

HTHI

因此H为正交阵。

引理2 设1kn, Hk为k阶Househoulder阵，则下列分块形式的矩阵

HInk 0 （5.22） T0 Hk仍是Househoulder阵，其中Ink为nk阶单位阵，0为nkk阶零阵。

n定理3 设b,uR* , bu, bn证明 取vbuR*，取

2u2，则存在Househoulder阵H，使得Hbu。

HI且

2Tvv TvvTHbb b2bubuTbubuT2bub bubTbbTuuTbuTuTnb

因为 b2u2， 所以bbuu。又因为bubiuiuTb，故

TTi1Hbb2bTbuTb2bbubTTbu

bbuu

17

Hbusgn(br1)b22br1bn

00b1注 ur1的符号由下列原则决定：“使ur1与br1异号”，否则，在bu中第r1个分量的计算中可能出现两个相近的数作减法。

Househoulder阵H的形状讨论如下： 令

000rvbubr1ur1记

Pbr2nrbn显然有v22Pnr22，且

0rTT20PrnrTP2vvHITInr2vvPnr20r,nr0r,r2T0PPnr,rnrnrI2Pnr2Ir,r0nr,r0r,nr2PnrPnTrPnr22

r1n，定理4 设已知向量b的nr1个分量不全为量，则必有一个Householder阵H，

使得Hb的前r个分量与向量b的前r个分量对应相等，而Hb的后nr1个分量都等于零。

证明 设

18

b1u1bu22  bbr， uur

bur1r1  bunn令

bi i1,,r2uisignbr1br21bn0 ir2,,nn满足ubR*且u2ir1

b2。

T由定理3知：对于向量b及与其等长的不同向量u，必定存在一个Househoulder阵

HI2bubuTbubu，使得对称阵A经过Househoulder变换后具有以下性质：

性质1 A与HAH相似。

性质2 HAH为对称阵。

性质3 HA的前r行与A的前r行相同。 证明 将H、A分别写成分块形式：

Ir,rH0nr,r则

0r,nr2PnrPnTrPnr22I 0Ar;记A0 DCB An-rAHArDCB DAnr性质4 HAH的前r列与HA的前r列相同。 证明

ArHAHDC BIr DAn-r00 D

ArDCnn

DAnrD,Hk kn2,使由递推公式

B定理5 设ARsym，则存在Househoulder阵H1,H2,

19

A1A

Ai1HiAiHi i1,得到的对称矩阵Ak1是三对角阵。

证明略。

例3 用Househoulder变换将对称阵

,k （5.23）

62A3128349118 17化为三对角阵。

解 n4，则至多通过n2422次Househoulder变换可将对称阵A化为三对角阵。

第一步 令A1A,记

62b1

31取

6223212u100则

614 000214, vv1b1u11312242.9666

0001T0-0.5345 -0.8018-0.2673v1v1 H1I22v120 -0.8018-0.5811 -0.13960-0.2673-0.1396 0.9535所以

63.7413.741713.875H1A1H101.8070-3.139

70061.8079 -3.1394A2

94.2912-6.00794-6.00792.851220

第二步 记

3.74113.875b21.8073.139取

76 943.741713.87563.6228 03.741713.8756u2221.80793.13940则

00, vv2b2u225.43073.13941T02v2v2H2I2v2200所以

2239.3483

000100 0-0.49900.866600.86660.49906 -3.741700 -3.741713.8576 -3.62270A H2A2H230 -3.62278.4058-3.638700-3.6387-1.2634显然，A3为所求对称三对角阵。

用Househoulder变换将对称阵A化为三对角阵算法所涉及程序设计知识较多，这里不再介绍。

21