【孟】北交网络课程讲解特征值

引例：选择旅游地问题

有人打算外出旅游，他在选择旅游地时，主要考虑的因素有景色、费用、居住、饮食和旅途五个因素，若这些因素依次用符号C1,C2，C3,C4，C5表示，则他选择旅游地的重要性可以用如下所谓成对比较矩阵表示，这里,A的任一元素Aij表示Ci与Cj对旅游地重要性之比，其值按Saaty等人提出的1—9尺度确定（见文献［15］，310页）,试求出这些因素对他选择旅游地的权重值。这个问题可以由层次分析法解决，它的层次结构为：

设W1、W2、…、W5依次为景色、费用、居住、饮食及旅游的权重,令向量，由Saaty的研究结果，可知权重W是：成对比较矩阵A的绝对值最大的特征值λmax所对应的归一化特征向量，即W满足：

，于是本问题归结为求矩阵A的绝对值最大的特征值及其对应的特征向量问题对求特征值机特征向量问题，线性代数已经给出解法，但它不适于计算机处理，且当矩阵阶较大时更是求解困难，本章主要介绍怎样在计算机上解决此问题。

问题的描述与基本概念求矩阵的特征值及特征向量的问题在实际问题中也经常遇到,如在工程技术中的振动问题和稳定性问题等，在这些问题中解出特征值或特征向量的计算机解法也称为代数特征问题的计算方法。下面介绍一下有关特征值问题的概念。定义1

设矩阵A∈Rn×n，若存在某个实数或复数λ及非零向量X∈Rn满足AX=λX，则称λ是A的一个特征值，而X称为λ对应的一个特征向量。定义2

称关于变量λ的行列式

为矩阵A的特征多项式，而FA(λ)=0称为特征方程.特征多项式FA（λ)是关于λ的一个n次多项式，线性代数中指出：矩阵A的特征值就是其特征多项式FA（λ）的零点，因此，n阶矩阵A共有n个特征值。求解矩阵A的特征值和特征向量的过程在线性代数中描述为1）求出特征方程FA(λ)=0的根λ1,λ2，…，λn2）对每个特征值λi，求出齐次线性方程组

的基础解系做为λi对应的特征向量，

i=1,2，…，n上述揭发理论很严密，但由于将特征多项式FA(λ）化为一个n次多项式很复杂且特征方程对舍入误差很敏感，特别当n较大时,这些问题更突出，由于这些原因，现在用计算机求解代数特征值问题不用线性代数的方法而用迭代加变换的处理方法，它们具有编程简单，对舍入误差不敏感等优点，本章将介绍具有代表性的这类问题的计算机解法：幂法和反幂法,旋转法及QR方法。

QR方法QR方法是求任意矩阵的全部特征值的一种有效方法,它是JACOBI方法的推广。基本思想利用矩阵的QR分解，通过逆序相乘产生对原矩阵的一系列正交相似变换,使其变化为一个近似的上三角矩阵来求全部特征值.这里QR分解是指将矩阵化为一个正交矩阵Q和一个上三角矩阵左乘的形式。构造原理实对称矩阵可用正交相似变换将其化为对角形矩阵，但对非对称矩阵,一般用正交相似变换化不成对角矩阵,但SCHUR分解定理给我们一个有关这方面的结果。定理3。（实SCHUR分解定理）设矩阵A∈Rn*n，则存在一个正交矩阵Q∈Rn＊n，使QTAQ=其中每个Bii是1*1或2*2的小矩阵,若Bii为1＊1的，其元素就是A的实特征值，否则Bii的特征值是A一对共轭复特征值.此定理的证明可参阅文献[3］。定理3指出了求矩阵A的全部特征值也可用正交相似变换的方法来做，正交相似变换的结果虽然不是对角矩阵，而是分块三角形矩阵，但它同样能很方便地求出全部特征值,有关一般矩阵的正交相似变换，我们不加证明地给出一个结论。定理4。设非奇异矩阵A∈Rn*n，且有n个不同的特征值，记A=A（1)。如果对整数k，有矩阵A(k)的QR分解为A(k）=QkRk,则令A(k+1)=QTkA(k）Qk，当k→∞时有A(k)本质上收敛于分块上三角形矩阵，这里“本质上收敛”指A（k)的主对角线上的元素或子块有确定的极限，其它元素或子块不管是否有极限。此定理给出了求解一般矩阵全部特征值的方法.由定理3，A（k+1)=（Q1Q2.。。.Qk）TA（Q1Q2..。.Qk)，令,则Qk也是正交矩阵，A(k+1）=QTkA（k）Qk说明A(k+1）也是原矩阵A的正交相似变换，从而A（k+1）与A有相同的特征值，n任意，此外，由A（k）=QkRk，则有QTkA（k）=

QTkA（k）Rk=Rk,故有A(k+1)=QkRk（应该是RkQk)，这说明A(k+1)可直接交换Qk与Rk的乘积顺序得到，于是可的如下QR算法.①对A(k)作QR分解A(k)=QkRk.②逆序相乘A(k）的分解矩阵，A(k）=RkQk。③判别A(k+1）是否为主对角线为1＊1或2＊2的子块形式的分块上三角形矩阵，若是对角线上各子块的特征值为所求特征值，终止，否则k+1k,转①.分析从QR算法的构造过程可以看到算法的主要计算量出现在QR分解上，如果直接对矩阵A用QR方法求全部特征值,那麽涉及的计算量是很大的，因此应该先对A作预处理.应用中常先对做正交相似变换将其化为上Hessenberg矩阵H，然后再对H采用QR方法，可以大大减少计算量，这里Hessenberg矩阵也称为拟三角矩阵，它的非零元素比三角矩阵多了一条次对角线,其形式为：

上Hessenberg矩阵

下Hessenberg矩阵实际上，Hessenberg矩阵虽然不是三角矩阵，但它很接近三角矩阵,由于其每列只比三角矩阵多一个非零元，故选用旋转变换做QR分解更简单些，因为对上Hessenberg矩阵H的第1列到第n—1列,依次做旋转变换使H的主对角线下的元素都变为零,则H化为上三角矩阵R了，用矩阵描述就是记，则J为正交矩阵，解出H，可得H=J-1R，因为J—1也是正交矩阵,于是得H的QR分解容易验证按这个方法做对H做QR分解，然后使用QR算法则构造的迭代序列都是上Hessenberg矩阵中进行,于是整个QR算法都在上Hessenberg矩阵中进行,这当然使QR算法的计算量大量减少。下面我们来具体讨论一下一般矩阵相似约化到Hessenberg矩阵的方法，为说明此问题，引入镜面反射阵概念.定义4。设非零向量V=（V1,V2，.。。。。,Vn）T∈Rn,则称矩阵P=I－β-1VVT为Hessenberg矩阵，式中.易验证Householder矩阵P是对称,正交和对合的.从定义可知Householder矩阵主要由一个非零向量V确定，若将V看某一过原点的超平面π的法向量，则有任给一个非零向量α,经Householder矩阵P作用后，记为Pα，则α与Pα是关与超平面π对称，因此也称Householder矩阵是镜面反射阵。Householder矩阵可改变任一向量的方向，这可从下面定理得出。定理5.任取非零向量X=（X1,X2,。.。..,Xn)T∈Rn，可以选择一个Householder矩阵P，使Px=－бｅ1式中e1=(1,0，。...。.,0）T是Rn的单位向量,证：作u=x+бｅ1，用u做一个Householder矩阵P=I－β—1uuT，因为而证毕定理中β=б(б+x1)，为避免б+x1出现两个相似数相减引起有效数字损失,实用中常选这样可得Householder矩阵的计算公式为计算计算u1=x1+б计算β=бu1计算P=I-β-1xxT有了Householder矩阵，现在可以考虑怎样化Hessemberg矩阵了。设矩阵A=（aij）∈Rn＊n，取A的第一列中后n-1个元素利用构造一个Householder矩阵，使，,再取对称的正交变换对A做正交相似变换。类似地取A（1）的第2列后n-2个元素做上述处理,可将A(1)化为：继续下去，做到第n—1次，则可将A化为上Hessenberg阵了。利用Householder矩阵同样可对矩阵做QR分解,有关这方面的内容可参见文献[2]。

例题设矩阵用某种数值方法求A的所有特征和特征向量.解:因为A是实对称矩阵，故可用JACOBI法求A的所有特征和特征向量，由JACOBI算法取第一个旋转矩阵，（因为对）则有同理取得继续下去，最后有故A的所有特征值为：此外于是对应的3个特征向量为:A的准确特征值为:可见JACOBI方法的精度很好设，用Householder相似变换将A化为上Hessenberg矩阵。解：因为A的第一列元素为（1，8,4）T，故取а1=（8,4)T做一个2阶Householder矩阵,取对称正交变换矩阵有矩阵B就是所求.用QR方法求矩阵的全部特征值。解：记A（0)=A，采用旋转变换将A(0）做QR分解，得A（0)

=Q0R0做逆序相乘继续做下去，经过11次计算后，有于是A有两个实特征值为及一对共轭复特征根，它来自2＊2矩阵块的特征根.

简评本章但绍的幂法,Jacobi方法及QR方法是求矩阵特征值和特征向量的常用数值方法，它们都是造构造迭代产生的矩阵序列来达到目的的。幂法计算简单，特别适用于高阶稀疏矩阵，但其收敛速度不能令人满意,要想加快幂法的收敛速度可采用反幂法及位移技术,在关主方面的内容可参看文献。Jacobi方法是古典方法，它收敛快，精度高,便于并行计算且算法稳定。用Jacobi方法求出的特征向量在较好的正交性,不过它的计算量较大，当阶数n增大时收敛速度减慢，因此Jacobi方法适用于求低阶的对称矩阵的全部特征值和特征向量。QR方法是60年代发展起来的，被人们称为数值数学，最值得注意的算法之一，它是目前求任意矩阵全部特征值和特征向量最有效的方法。

Jacobi方法Jacobi方法也称旋转法,是求实对称矩阵的全部特征值和特征向量的方法基本思想将实对称矩阵进行一系列的相似正交变换使其约化成一个近似对角矩阵,然后利用相似正交变换的关系来求全部特征值和特征向量。由于这里的正交相似变换主要采用旋转变换，这就是称Jacobi方法为旋转法的原因。1

构造原理线性代数理论指出:若矩阵A是实对称矩阵，则

一定正交相似于对角形，即存在一个正交矩阵

和对角矩阵

，满足

。因相似矩阵具有相同的特征值，而对角矩阵的特征值就是对角线上的

个元素，由此可得原对称矩阵

的个特征值，此外由相似关系，可得正交矩阵

的

个列向量就是A的n个线性无关特征向量。于是求实对称矩阵A的全部特征值和特征向量问题可转化为求正交相似矩阵P及相似对角形λ的问题。为构造相应的算法，注意到实对称矩阵总与二次型一一对应的，当一个二次型中没有交叉项时，其对应的实对称矩阵就是对角形，于是将一个二次型通过正交变换化为无交叉项的二次型指出了寻找正交矩阵的方法。为找出一种具体的方法，先考虑两个自变量的完全二次型的变化情况：设两个变元的二次型为将写成矩阵形式就是取二阶旋转变换，

则有这里要使二次型无交叉变换,只要适当选择,使即可,由此得应取为使，这样二次型变为无交叉项了。若令A

P

λ则有PTAP=P-1AP=λ。

P也称旋转矩阵，它是正交阵。这说明用旋转矩阵做正交相似变换可将实对称矩阵化为对角形。为得到n阶旋转矩阵，将2阶旋转阵作如下推广定义3.

n阶矩阵称为

n阶旋转矩阵或Givens矩阵.易验证J（i,j，

)—1=

J（i，j,

)T，故J（i,j，

)是正交矩阵。若将用J（i，j,

)对实对称矩阵A=（aij）做正交相似变换，则有J（i,j,

)TA

J（i,j，

)=A(1）=(aij```

)直接计算会发现A只变化了第i,j两行和i，j两列，即选取,则有，于是经一次旋转正交变换的实对称矩阵A变为A（1),它比A更靠近对角形一步(因为有两个非对角线元素变为零了），但有与A相同的特征值.如果每次都选不同的旋转矩阵J（i，j，

）做正交相似变换,就可期望变换若干步后,最终获得所要的对角形矩阵,这个期望可由如下定理给予证实。定理2。

设实对称矩阵是n阶旋转矩阵序列,若记这里A（k)与A(k—1)相比是A（k)的,则有其中λ1,λ2,。。。λn是A的n个特征值。定理1的证明可参看文献[2].利用定理的结果，若记P=J1J2。。.Jn...，则P的列向量就是A的特征向量。实用中不可能做无穷次正交相似变换，一般用非对角元素的平方和是否小于某个小数ξ来终止计算,具体算法为（设P=I,ξ为给定精度)①先选矩阵A（k)的非对角元素的绝对值最大者得到行标P

，列标q②计算旋转角得到旋转矩阵,并做③计算④计算⑤判别,若成立，得A（k+1）的主对角线元素为所求的全部特征值，P的n个列为相应的n个特征向量，否则令,转①通常把按上述算法求实对称矩阵全部特征值及其特征向量的方法称为Jacobi方法.2

分析Jacobi方法每次做由A(k）到A(k+1)的正交相似变换时，只改变了A（k)的ik,jk行和ik,jk列，由计算特点,可以将每次计算的A(k）都存放在同一个矩阵A中，这样可以节约存储空间并减少计算量。此外,Jacobi方法在每次计算前，先要选择非对角矩阵元素的绝对值最大者来确定旋转变换矩阵，这在n很大时是很费时的，再者,前面旋转相似变换矩阵已经变为零的非对角元素在后面的相似变换中可能还会变为非零元素。因此应在这些方面对Jacobi方法做些改进是之能更实用些，这里介绍一个实用的方法称为过关Jacobi方法,此方法是先取阀值，然后按行的顺序a12,a13,.。。a1n，a21，a23.。。an-1n用与之相比，若，则不做变换，若，则做变换将aij化为零，如此多次直至所有非对角元素的绝对值都小于后，再选，重复上述过程，继续下去选阀值，直至，则终止计算，这里，当然也可以选其他的值，只要单调减小即可。过关Jacobi方法使原来的法确定旋转矩阵的过程大大缩短，因此可以加快求解过程。幂法与反幂法一﹑问题的提出:工程技术的许多实际问题，例如振动问题，稳定问题的求解,有时会归结成求矩阵的特征值λ和对应的特征向量χ。学过线性代数后，我们已知求矩阵A的特征值λ和特征向量χ的解法，即先求出A的特征多项式:令ｆ﹙ｘ﹚﹦0。通过求解上述高次多项式方程，所得根λ即为矩阵A的特征值,然后求解方程组﹙Ａ﹣λＩ﹚Ｘ﹦0，就可得出特征值λ对应的特征向量X。上面用定义阐述了如何求解矩阵A的特征值λ和特征向量X。但众所周知，高次多项式求根是相当困难的，而且重根的计算精度较低。同时,矩阵A求特征多项式系数的过程对舍入误差十分敏感,这对最后计算结果影响很大。因此,从数值计算角度来看，上述方法缺乏实用价值。二﹑问题的解决：目前，求矩阵特征值问题实际采用的是迭代法和变换法.这里将介绍通过求矩阵特征向量求出特征值的一种迭代法-——-幂法，而后再介绍一些反幂法的内容。（1）幂法定理：设矩阵A的特征值为并设A有完全的特征向量系(它们线性无关），则对任意一个非零向量所构造的向量序列有其中表示向量的第j个分量.证

仅就为实数的情况来证明。假定于是，由矩阵特征值定义知,得同理可得假定，因为,故得证毕从上述证明过程可得出计算矩阵A的按模最大特征值的方法,具体步骤如下：任取一非零向量，一般可取当k足够大时，即可得到：若按上述计算过程，有一严重缺点，当(或时）中不为零的分量将随K的增大而无限增大，计算机就可能出现上溢（或随K的增大而很快出现下溢)，因此,在实际计算时，须按规范法计算，每步先对向量进行“规范化"，即取中绝对值最大的一个向量记作，用遍除的所有向量，得到规范化向量。幂法的规范化算法：说明

给定非零向量,求n阶矩阵A的按模最大特征值输入

维数n,矩阵A，向量容许误差输出

按模最大特征值的近似值和对应的特征向量步骤：计算3ifthenstop4

,goto2为说明上述算法的正确性，我们试证明下述定理定理二：在定理一的条件下,规范化向量序列收敛于矩阵A按模最大的特征值对应的特征向量，而向量序列的绝对值最大的分量收敛于，即与证

一般故又故现在举一例说明幂法的使用方法:1．用幂法求矩阵按模最大特征值和对应的特征向量解取初始向量，计算出和迭代7次的结果列于下表k0123456711127495—18444.4327714。84322—29。6426244.9233314.97623-29.9504844。9957214。99865-29.9972244。9995914.99988-29.9997444。9995314。99983—29。9996844.9995314.99983—29.9996811110.34672-0。6715310。33413-0。6672710.33337-0。6667010。33334-0.6666710。33333-0.6666710.33333—0.6666710。33333—0。66667由上可见经过7次迭代,的值已稳定到小数后5位,故所求的按模最大特征值和对应的特征向量可取作：本题矩阵的三个特征值和对应的特征向量不难求出为与将规范化即得。幂法的几点解释:(1）若矩阵的按模最大特征值是r重根时，上述两定理仍成立。证：由可得(2)定理一的证明中，所取初始向量的这个假定是不可少的，否则就可能得不到按模的最大特征值。（3）幂法的收敛速度,虽然与初始向量的选取有关,但主要取决于比值的大小，当时,幂法收敛很慢,甚至失去使用价值，若越小，则收敛速度越快，平移原点法就是基于上述基本思想来加快收敛速度的.具体说明：若是矩阵A的特征值，其相应特征向量为

即A=则有：-p为矩阵A-pI的特征值（其中p为一数值）事实上：（A—pI)=（—p）

如果满足：

而且还满足：如果能选择这样的p，使成立,则用同一初始向量用幂法求A—pI及A的特征值,收敛速度要快得多。因此对满足上述条件的p,可先求A-pI的按模最大特征值及特征向量，从而