第5章矩阵特征值计算
与线性方程组的求解问题一样，矩阵特征值与特征向量的计算也是数值线性代数的重要内容。理论上，矩阵的特征值是特征多项式方程的根，因此特征值的计算可转换为单个多项式方程的求解。然而，对于高阶矩阵，这种转换并不能使问题得到简化，而且在实际应用中还会引入严重的数值误差。因此，正如第2章指出的，人们一般将多项式方程求解转换为矩阵特征值计算问题，而不是反过来。
本章介绍有关矩阵特征值计算问题的基本理论和算法。与非线性方程求根问题类似，计算矩阵特征值的算法也是迭代方法 如果用有限次运算能求得一般矩阵的特征值，则多项式方程求根问题也可用有限次运算解决，这与阿贝尔证明的“高于4次的多项式并不都有用初等运算表示的求根公式”的理论矛盾。。
5.1基本概念与特征值分布
本节首先介绍矩阵特征值、特征向量的基本概念和性质，然后讨论对特征值分布范围的简单估计方法。
5.1.1基本概念与性质
定义5.1： 矩阵A=(akj)∈Cn×n，
(1) 称φ(λ)=det(λI-A)=λn+c1λn-1+…+cn-1λ+cn为A的特征多项式(characteristic polynomial)；n次代数方程φ(λ)=0为A的特征方程(characteristic equation)，它的n个根λ1,λ2，…,λn被称为A的特征值（eigenvalue）。此外，常用λ(A)表示A的全体特征值的集合，也称为特征值谱特征值谱(spectrum of eigenvalue)。
(2) 对于矩阵A的一个给定特征值λ，相应的齐次线性方程组(λI-A)x=0(5.1)有非零解（因为系数矩阵奇异），其解向量x称为矩阵A对应于λ的特征向量（eigenvector）。
根据式(5.1)得出矩阵特征值与特征向量的关系，即Ax=λx 。(5.2)第3章的定义3.5就利用式(5.2)对矩阵特征值和特征向量进行了定义，它与定义5.1等价。另外，同一个特征值对应的特征向量一定不唯一，它们构成线性子空间，称为特征子空间特征子空间(eigenspace)。
我们一般讨论实矩阵的特征值问题。注意，实矩阵的特征值和特征向量不一定是实数和实向量，但实特征值一定对应于实特征向量（式(5.1)的解），而一般的复特征值对应的特征向量一定不是实向量。此外，若特征值不是实数, 则其复共轭也一定是特征值（由于特征方程为实系数方程）。定理3.3表明，实对称矩阵A∈Rn×n的特征值均为实数，图51弹簧质点系统
存在n个线性无关且正交的实特征向量，即存在由特征值组成的对角阵Λ和由特征向量组成的正交阵Q，使得A=QΛQT(5.3)例5.1（弹簧质点系统）： 考虑图51的弹簧质点系统，其中包括3个质量分别为m1、m2、m3的物体，由3个弹性系数分别为k1、k2、k3的弹簧相连，3个物体的位置均为时间的函数，这里考查3个物体偏离平衡位置的位移，分别记为y1(t)、y2(t)、y3(t)。因为物体在平衡状态所受的重力已经和弹簧伸长的弹力平衡，所以物体的加速度只和偏离平衡位置引起的弹簧伸长相关。根据牛顿第二定律以及胡克定律（即弹簧的弹力与拉伸长度成正比）可列出如下微分方程组本书第8章将介绍这种常微分方程组的数值求解方法。： My″(t)+Ky(t)=0其中y(t)=\[y1(t)y2(t)y3(t)\]T, M=m100
0m20
00m3,K=k1+k2-k20
-k2k2+k3-k3
0-k3k3 。一般情况下，这个系统会以自然频率ω做谐波振动，而y的通解包含如下分量： yj(t)=xjeiωt,(j=1,2,3),其中，i=-1，根据它可求解振动的频率ω及振幅xj。由这个式子可得出y″j(t)=-ω2xjeiωt,(j=1,2,3),代入微分方程，可得代数方程-ω2Mx+Kx=0,或Ax=λx,其中，A=M-1K，λ=ω2。通过求解矩阵A的特征值可求出这个弹簧质点系统的自然频率（有多个）。再结合初始条件可确定这3个位移函数，它们可能按某个自然频率振动（简正振动），也可能是若干个简正振动的线性叠加。■
例5.2（根据定义计算特征值、特征向量）： 求矩阵A=5-1-1
31-1
4-21的特征值和特征向量。
【解】矩阵A的特征方程为det(λI-A)=λ-511
-3λ-11
-42λ-1=(λ-3)(λ-2)2=0,故A的特征值为λ1=3，λ2=2（二重特征值）。
当λ=λ1=3时，由(λI-A)x=0，得到方程-211
-321
-422x1
x2
x3=0
0
0 ,它有无穷多个解，若假设x1=1, 则求出解为x=\[1,1,1\]T，记为x1，则x1是λ1对应的一个特征向量。
当λ=λ2=2时，由(λI-A)x=0，得到方程-311
-311
-421x1
x2
x3=0
0
0 ,它有无穷多个解，若假设x1=1, 则求出解为x=\[1,1,2\]T，记为x2，则x2是λ2对应的一个特征向量。■
下面概括地介绍有关矩阵特征值、特征向量的一些性质，它们可根据定义5.1以及式(5.2)证明。
定理5.1： 设λj(j=1,2,…,n)为n阶矩阵A的特征值，则
(1) ∑nj=1λj=∑nj=1ajj=tr(A)。
(2) ∏nj=1λj=det(A)。
这里，tr(A)表示矩阵对角线上元素之和，称为矩阵的迹（trace）。
从上述结论(2)也可以看出，非奇异矩阵特征值均不为0, 而0一定是奇异矩阵的特征值。
定理5.2： 矩阵转置不改变特征值，即λ(A)=λ(AT)。
定理5.3： 若矩阵A为对角阵或上(下)三角阵，则其对角线元素即矩阵的特征值。
定理5.4： 若矩阵A为分块对角阵，或分块上(下)三角阵，例如A=A11A12…A1m
A22…A2m
[3]︙
[4]Amm ,其中，每个对角块Ajj均为方阵，则矩阵A的特征值为各对角块矩阵特征值的合并，即λ(A)=∪mj=1λ(Ajj)。
定理5.5： 矩阵的相似变换相似变换(similarity transformation)不改变特征值。设矩阵A和B为相似矩阵，即存在非奇异矩阵X使得B=X-1AX，则
(1) 矩阵A和B的特征值相等，即λ(A)=λ(B) 。
(2) 若y为B的特征向量，则相应地，Xy为A的特征向量。
通过相似变换并不总能把矩阵转换为对角阵，或者说矩阵A并不总是可对角化的(diagonalizable)。下面给出特征值的代数重数代数重数、几何重数代数重数和亏损矩阵亏损矩阵的概念，以及几个定理。
定义5.2: 设矩阵A∈Rn×n有m个（m≤n）不同的特征值λ～1,λ～2,…,λ～m,若λ～j是特征方程的nj重根，则称nj为λ～j的代数重数(algebraic multiplicity)，并称λ～j对应的特征子空间(Cn的子空间)的维数为λ～j的几何重数(geometric multiplicity)。
定理5.6： 设矩阵A∈Rn×n的m个不同的特征值为λ～1,λ～2,…,λ～m，特征值λ～j(j=1,2，…,m)的代数重数为nj，几何重数为kj，则
(1) ∑mj=1nj=n，且任一个特征值的几何重数都不大于代数重数，即j，nj≥kj。
(2) 不同特征值的特征向量线性无关，并且将所有特征子空间的∑mj=1kj个基(特征向量)放在一起，它们构成一组线性无关向量。
(3) 若每个特征值的代数重数等于几何重数，则总共可得n个线性无关的特征向量，它们是全空间Cn的基。
定义5.3： 若矩阵A∈Rn×n的某个代数重数为k的特征值对应的线性无关特征向量数目少于k（即几何重数小于代数重数），则称A为亏损矩阵（defective matrix），否则称其为非亏损矩阵非亏损矩阵（nondefective matrix）。
定理5.7: 矩阵A∈Rn×n可对角化的充要条件是A为非亏损矩阵。若A可对角化，即存在非奇异矩阵X∈Cn×n使得X-1AX=Λ,其中，Λ∈Cn×n为对角矩阵, 则Λ的对角线元素为矩阵A的特征值，而矩阵X的列向量为n个线性无关的特征向量。
定理5.7中方程的等价形式为A=XΛX-1, 它被称为特征值分解，也被称为谱分解(spectrum decomposition)。特征值分解存在的充要条件是A为非亏损矩阵。但现实中还有很多矩阵是亏损矩阵，如例5.2中的矩阵，它的特征值2的代数重数为2，而几何重数仅为1。这种矩阵不能相似变换为对角矩阵，但存在下面的约当分解约当（Jordan）分解(Jordan decomposition)。
定理5.8： 设矩阵A∈Rn×n, 存在非奇异矩阵X∈Cn×n使得A=XJX-1,矩阵J为形如J1

Jp的分块对角矩阵（称为约当标准型约当标准型），其中Js=λs1
λs
[3]1
[4]λs，1≤s≤p，称为约当块约当块，其对角线元素为矩阵A的特征值。设矩阵A有m个不同的特征值为λ～1,λ～2,…,λ～m，特征值λ～j(j=1,2，…,m)的代数重数为nj，几何重数为kj，则p=∑mj=1kj, λ～j对应kj个约当块(对角元均为λ～j), 其阶数之和等于nj。
在约当分解中，如果所有约当块都是1阶的，则J为对角矩阵，这种分解就是特征值分解，相应的矩阵为非亏损矩阵。约当分解是很有用的理论工具，利用它还可证明下面关于矩阵运算结果的特征值的定理。
定理5.9： 设λj(j=1,2,…,n)为n阶矩阵A的特征值，则
(1) 矩阵cA(c为常数)的特征值为cλ1,cλ2,…,cλn。
(2) 矩阵A+cI(c为常数)的特征值为λ1+c,λ2+c,…,λn+c 。
(3) 矩阵Ak(k为正整数)的特征值为λk1,λk2,…,λkn。
(4) 设p(t)为一多项式函数，则矩阵p(A)的特征值为p(λ1),p(λ2)，…,p(λn) 。
(5) 若A为非奇异矩阵，则λj≠0(j=1,2,…,n),且矩阵A-1的特征值为λ-11,λ-12,…,λ-1n。
5.1.2特征值分布范围的估计
估计特征值的分布范围或它们的界，无论在理论上或实际应用上，都有重要意义。例如，本书前面的内容曾涉及两个问题。 
(1) 计算矩阵的2条件数：  cond(A)2=λmax(ATA)λmin(ATA) 。
(2) 考查1阶定常迭代法x(k+1)=Bx(k)+f的收敛性、收敛速度，收敛的判据是谱半径ρ(B)=max1≤j≤n|λj(B)|<1,收敛速度为R=-lgρ(B) 。
其中，都需要了解矩阵特征值的分布范围。
第4章的定理4.4说明谱半径的大小不超过任何一种算子范数，即ρ(A)≤‖A‖,这是关于特征值的上界的一个重要结论。
下面先给出定义5.4，再介绍有关特征值的界的另一个重要结论。
定义5.4： 设A=(akj)∈Cn×n，记rk=∑nj=1j≠k|akj|(k=1,2，…,n),则集合Dk={z:|z-akk|≤rk,z∈C}(k=1,2,…,n)在复平面为以akk为圆心、rk为半图52复坐标平面，以及3×3矩阵A
的格什戈林圆盘格什戈林圆盘
径的圆盘，称为A的Gerschgorin（格什戈林）圆盘Gerschgorin（格什戈林）圆盘。
图52显示了一个3×3复矩阵的格什戈林圆盘。
定理5.10(圆盘定理圆盘定理)： 设A=(akj)∈Cn×n，则 
(1) A的每个特征值必属于A的格什戈林圆盘中，即对任一特征值λ,必定存在k(1≤k≤n)，使得|λ-akk|≤∑nj=1j≠k|akj|(5.4)用集合的关系说明，这意味着λ(A)∪nk=1Dk。
(2) 若A的格什戈林圆盘中有m个圆盘组成一连通并集S，且S与余下的n-m个圆盘分离，则S内恰好包含A的m个特征值（重特征值按重数计）。
对如图52所示的例子，定理5.10的第(2)个结论的含义是： D1中只包含一个特征值，而另外两个特征值在D2、D3的并集中。下面对定理5.10的结论(1)进行证明，结论(2)的证明超出了本书的范围。
【证明】设λ为A的任一特征值，则有Ax=λx，x为非零向量。设x中第k个分量最大，即|xk|=max1≤j≤n|xj|>0,考虑式(5.2)中第k个方程∑nj=1akjxj=λxk,将其中与xk有关的项移到等号左边，其余项移到右边，再两边取模得|λ-akk||xk|=∑nj=1j≠kakjxj≤∑nj=1j≠k|akj||xj|≤|xk|∑nj=1j≠k|akj|(5.5)最后一个不等式的推导利用了“x中第k个分量最大”的假设。将不等式(5.5)除以|xk|，得到式(5.4)，因此证明了定理5.10的结论(1)。上述证明过程还说明，若某个特征向量的第k个分量的模最大，则相应的特征值必属于第k个圆盘。■
根据定理5.2，还可以按照矩阵的每列元素定义n个圆盘，对于它们定理5.10仍然成立。下面的定理是圆盘定理的重要推论，其证明留给感兴趣的读者。
定理5.11： 设A∈Rn×n，且A的对角元均大于0，则
(1) 若A严格对角占优，则A的特征值的实部都大于0。
(2) 若A为对角占优的对称矩阵，则A一定是对称半正定矩阵。若同时A非奇异，则A为对称正定矩阵。
例5.3(圆盘定理的应用)： 试估计矩阵A=410
10-1
11-4的特征值范围。
【解】直接应用圆盘定理，该矩阵的3个圆盘如下:D1: |λ-4|≤1,D2: |λ|≤2,D3: |λ+4|≤2 。D1与其他圆盘分离，则它仅含一个特征值，且必定为实数（若为虚数，则其共轭也是特征值，这与D1仅含一个特征值矛盾）。所以，对矩阵特征值的范围的估计是3≤λ1≤5,λ2,λ3∈D2∪D3 。再对矩阵AT应用圆盘定理，则可以进一步优化上述结果。矩阵AT对应的3个圆盘为D′1: |λ-4|≤2,D′2: |λ|≤2,D′3: |λ+4|≤1 。这说明D′3中存在一个特征值，且为实数，它属于区间［-5, -3］，经过综合分析可知3个特征值均为实数，它们的范围是λ1∈\[3,5\],λ2∈\[-2,2\],λ3∈\[-5,-3\] 。事实上，使用MATLAB的eigeig命令可求出矩阵A的特征值为4.2030,-0.4429,-3.7601。■
根据定理5.5，还可以对矩阵A做简单的相似变换，如取X为对角阵，然后再应用圆盘定理估计特征值的范围。
例5.4(特征值范围的估计)： 选取适当的矩阵X，应用定理5.5和定理5.10估计例5.3中矩阵的特征值范围。
【解】取X-1=100
010
000.9则A1=X-1AX=410
10-10/9
0.90.9-4的特征值与A的相同。对A1应用圆盘定理，得到3个分离的圆盘，它们分别包含一个实特征值，由此得到特征值的范围估计λ1∈\[3,5\],λ2∈-199,199,λ3∈\[-5.8,-2.2\]此外，还可进一步估计ρ(A)的范围，即3≤ρ(A)≤5.8 。■
例5.4表明，综合运用圆盘定理和矩阵特征值的性质（如定理5.2、定理5.5），可对特征值的范围进行一定的估计。对具体例子，可适当设置相似变换矩阵，尽可能让圆盘相互分离，从而提高估计的有效性。
5.2幂法与反幂法幂法与反幂法反幂法
幂法是一种计算矩阵最大的特征值及其对应特征向量的方法。本节介绍幂法、反幂法以及加快幂法迭代收敛的技术。
5.2.1幂法
定义5.5： 在矩阵A的特征值中，模最大的特征值称为主特征值主特征值，也称为“第一特征值”，它对应的特征向量称为主特征向量主特征向量。
应注意的是，主特征值有可能不唯一，因为模相同的复数可以有很多。例如，模为5的特征值可能是5,-5,3+4i,3-4i, 等等。另外，请注意谱半径和主特征值的区别。
如果矩阵A有唯一的主特征值，则通过幂法能方便地计算出主特征值及其对应的特征向量。对于实矩阵，这个唯一的主特征值显然是实数，但不排除它是重特征值的情况。幂法(power iteration)的计算过程是： 首先随机选一非零向量�瘙
經0∈Rn，然后进行迭代计算�瘙
經k=A�瘙
經k-1,(k=1,2,…)得到向量序列{�瘙
經k}，根据它即可求出主特征与特征向量。下面用定理说明。
定理5.12: 设A∈ n×n，其主特征值唯一，记为λ1，随机选一非零向量�瘙
經0∈ n按迭代公式进行计算：�瘙
經k=A�瘙
經k-1,（k=1,2,…),则有
(1) 当k→ 时，�瘙
經k趋近于λ1的特征向量。
(2)limk→ �瘙
經k+1j�瘙
經kj=λ1,(5.6)
其中�瘙
經kj表示向量�瘙
經k的第j个分量，且j为�瘙
經k、�瘙
經k+1绝对值最大元素的下标。
【证明】这里只考虑矩阵A为非亏损矩阵的情况。设A的主特征值为p-1重特征值，全部n个特征值按模从大到小排列为: λ1=…=λp-1>λp≥…≥λn，它们对应于一组线性无关的单位特征向量x^1,x^2，…,x^n。向量�瘙
經0可写成这些特征向量的线性组合：�瘙
經0=α1x^1+…+αnx^n,一般地，α1≠0，则�瘙
經k=A�瘙
經k-1=Ak�瘙
經0=α1λk1x^1+α2λk2x^2+…+αnλknx^n
=λk1∑p-1j=1αjx^j+∑nj=pαjλjλ1kx^j
=λk1∑p-1j=1αjx^j+εk,(5.7)其中εk=∑nj=pαjλjλ1kx^j。由于λjλ1<1,(j=p,p+1,…,n), 则limk→ εk=0limk→ �瘙
經kλk1=∑p-1j=1αjx^j  。式（5.7）中，∑p-1j=1αjx^j是主特征值λ1的一个特征向量。由于特征向量放大、缩小任意倍数后仍是特征向量，则随k的增大，�瘙
經k越来越趋近于λ1的特征向量。
更进一步，特征向量为非零向量，其绝对值最大元素一定不为0，所以当k足够大时，�瘙
經kj≠0, �瘙
經k+1j≠0。因此，根据式(5.7),有�瘙
經k+1j�瘙
經kj=λ1∑p-1j=1αjx^j+εk+1j∑p-1j=1αjx^j+εkj由于limk→ εk=0，随k的增大，上式等号右边趋于一个常数λ1，这就证明了定理的结论。■
若矩阵A为亏损矩阵，则可利用矩阵的约当分解证明式(5.6)，这里略去。在这种情况下，式(5.6)的收敛速度可能很慢。
关于定理5.12，再说明下列两点。
(1) 理论上讲，定理5.12成立的前提条件应包括α1≠0。特别地，若初始向量�瘙
經0为某个非主特征值对应的特征向量，则式(5.6)求出的就是那个特征值。然而，在实践中，若随机选取�瘙
經0，且由于计算的舍入误差，就不会出现此情况。
(2) 式(5.6)的含义是相邻迭代向量分量的比值收敛到主特征值。其中j在实际计算时取�瘙
經k中最大的分量对应的j。
直接使用幂法，还存在如下两方面问题。 
(1) 溢出： 由于�瘙
經k≈λk1x1，则
|λ1|>1时，实际计算�瘙
經k会出现上溢出（当k很大时）。
|λ1|<1时，实际计算�瘙
經k会出现下溢出（当k很大时）。
(2) 可能收敛速度很慢。由于εk=∑nj=pαjλjλ1kxj，因此εk→0的速度取决于求和式中衰减最慢的因子λpλ1，也就是说，当绝对值第2大的特征值其绝对值接近主特征值的绝对值时，幂法的收敛速度就会很慢。
下面采用规格化向量规格化向量的技术防止溢出，导出实用的幂法。关于加速收敛技术的讨论，见5.2.2节。
定义5.6： 记max(�瘙
經)为向量�瘙
經∈Rn的绝对值最大的分量，即max(�瘙
經)=vj，其中j满足|vj|=max1≤k≤n|vk|, 若j的值不唯一，则取最小的那个，并且称u=�瘙
經/max(�瘙
經)为向量�瘙
經的规格化向量(normalized vector)。
例5.5（规格化向量）： 设�瘙
經=\[3,-5,0\]T，max(�瘙
經)=-5，对应的规格化向量为u=-35,1,0T。
根据定义5.6，容易得出规格化向量的两条性质。
定理5.13: 定义5.6中的规格化向量满足如下两条性质。
(1) 若u为规格化向量，则‖u‖∞=1，并且max(u)=1。
(2) 设向量�瘙
經1和�瘙
經2的规格化向量分别为u1和u2，若�瘙
經1=α�瘙
經2,实数α≠0，则u1=u2。
在幂法的每一步增加向量规格化的操作可解决溢出问题。先看第一步，�瘙
經1=A�瘙
經0，此时计算�瘙
經1的规格化向量u1=�瘙
經1max(�瘙
經1)=A�瘙
經0max(A�瘙
經0)然后使用规格化向量计算�瘙
經2�瘙
經2=Au1=A2�瘙
經0max(A�瘙
經0)(5.8)再进行向量规格化操作，u2=�瘙
經2max(�瘙
經2)=A2�瘙
經0max(A2�瘙
經0)(5.9)式(5.9)的推导利用了式(5.8)和定理5.13的结论(2)。依此类推，得到�瘙
經k=Auk-1=Ak�瘙
經0max(Ak-1�瘙
經0)
uk=�瘙
經kmax(�瘙
經k)=Ak�瘙
經0max(Ak�瘙
經0),(k=1,2,…)(5.10)根据定理5.12的证明过程(为了表述简单，这里不妨设主特征值λ1不是重特征值)，Ak�瘙
經0=λk1α1x^1+∑nj=2αjλjλ1kx^j
uk=Ak�瘙
經0max(Ak�瘙
經0)=α1x^1+∑nj=2αjλjλ1kx^jmaxα1x^1+∑nj=2αjλjλ1kx^jk→∞x1max(x1),即uk逐渐逼近规格化的主特征向量。同理，�瘙
經k=Auk-1=Ak�瘙
經0max(Ak-1�瘙
經0)=λk1α1x^1+∑nj=2αjλjλ1kx^jmaxλk-11α1x^1+∑nj=2αjλjλ1k-1x^j
=λ1α1x^1+∑nj=2αjλjλ1kx^jmaxα1x^1+∑nj=2αjλjλ1k-1x^j,因此，根据定理5.13的结论(1)，有limk→∞�瘙
經k=λ1x1max(x1)limk→∞max(�瘙
經k)=λ1 。基于上述推导，得到如下定理，以及如算法5.1描述的实用幂法实用幂法。应注意的是，由于x1≠0，则max(x1)≠0。
定理5.14: 设A∈Rn×n，其主特征值唯一，记为λ1，随机取一非零初始向量�瘙
經0，设u0=�瘙
經0，按迭代式(5.10)计算，则limk→∞uk=x1max(x1)(5.11)
limk→∞max(�瘙
經k)=λ1(5.12)其中，x1为主特征向量。算法5.1： 计算主特征值λ1和主特征向量x1的实用幂法算法5.1： 计算主特征值λ1和主特征向量x1的实用幂法
输入: A; 输出: x1,λ1.
u∶=随机向量;
While 不满足判停准则 do
�瘙
經∶=Au;
λ1∶=max(�瘙
經);{主特征值近似值}
u∶=�瘙
經/λ1; {规格化}
End
x1∶=u{规格化的主特征向量}
在算法5.1中，可根据相邻两步迭代得到的主特征值近似值之差判断是否停止迭代。每个迭代步的主要计算是算一次矩阵与向量乘法，若A为稀疏矩阵，则可利用它的稀疏性提高计算效率。实用的幂法保证了向量序列{�瘙
經k}、{uk}不溢出，并且向量�瘙
經k的最大分量的极限就是主特征值。
最后针对幂法的适用范围再说明两点。
(1) 若实矩阵A对称半正定或对称半负定，则其主特征值必唯一（而且是非亏损矩阵）。有时也可以估计特征值的分布范围，从而说明主特征值的唯一性。只有满足此条件，才能保证幂法的收敛性。
(2) 对一般的矩阵，幂法的迭代过程有可能不收敛，此时序列{uk}有可能包括多个收敛于不同向量的子序列，它趋向于成为多个特征向量的线性组合。但是，一旦幂法的迭代过程