第五章矩阵的特征值与特征向量问题.doc

第四章 矩阵的特征值与特征向量问题

物理、力学和工程技术中的许多问 题在数学上都归结为求矩阵的特征值和特征向量问题.计算方阵A的特征值，就是求特征方程

即

的根。求出特征值后，再求相应的齐次线性方程组

的非零解，即是对应于的特征向量。这对于阶数较小的矩阵是可以的，但对于阶数较大的矩阵来说,求解是十分困难,所以用这种方法求矩阵的特征值是不切实际的.
我们知道，如果矩阵A与B相似，则A与B有相同的特征值.因此人们就希望在相似变换下,把A化为最简单的形式。一般矩阵的最简单的形式是约当标准形。由于在一般情况下，用相似变换把矩阵A化为约当标准形是很困难的,于是人们就设法对矩阵A依次进行相似变换,使其逐步趋向于一个约当标准形，从而求出A的特征值.
本章介绍求部分特征值和特征向量的幂法，反幂法;求实对称矩阵全部特征值和特征向量的雅可比方法；求特征值的多项式方法；求任意矩阵全部特征值的QR方法.

第一节幂法与反幂法
一 幂法
幂法是一种求任意矩阵A的按模最大特征值及其对应特征向量的迭代算法.该方法最大的优点是计算简单，容易在计算机上实现，对稀疏矩阵较为合适,但有时收敛速度很慢.
为了讨论简单，我们假设
（1)n阶方阵A的特征值按模的大小排列为
(1）
(2）是对应于特征值的特征向量；
(3) 线性无关.
任取一个非零的初始向量，由矩阵A构造一个向量序列
（2）
称为迭代向量。由于线性无关,构成n维向量空间的一组基，所以,初始向量 可唯一表示成
(3）
于是
（4)
因为比值所以
(5）
当k充分大时有
(6)
从而
（7)
这说明当k充分大时,两个相邻迭代向量与近似地相差一个倍数,这个倍数便是矩阵A的按模最大的特征值.若用表示向量的第个分量，则
（8)
也就是说两个相邻迭代向量对应分量的比值近似地作为矩阵A的按模最大的特征值.
因为，又,所以有,因此向量可近似地作为对应于的特征向量。
这种由已知的非零向量和矩阵A的乘幂构造向量序列以计算矩阵A的按模最大特征值及其相应特征向量的方法称为幂法。
由（4）式知,幂法的收敛速度取决于比值的大小.比值越小，收敛越快，但当比值接近于1时，收敛十分缓慢。
用幂法进行计算时，如果,则迭代向量的各个不为零的分量将随着k无限增大而趋于无穷。反之,如果，则的各分量将趋于零.这样在有限字长的计算机上计算时就可能溢出停机。为了避免这一点，在计算过程中，常采用把每步迭代的向量进行规范化，即用 乘以一个常数,使得其分量的模最大为1。这样,迭代公式变为
（9)
其中是模最大的第一个分量。相应地取
(10）
例1 设

用幂法求其模为最大的特征值及其相应的特征向量（精确到小数点后三位）。
解 取，计算结果如表4—1所示.
表4-1

1
1
0
1
1
1
0
1
2
2
—2
2
2
1
—1
1
3
3
—4
3
—4
-0.75
1
-0.75
4
—2.5
3.5
-2.5
3.5
—0.714
1
—0.714
5
—2。428
3。428
—2。428
3.428
—0。708
1
—0。708
6
—2.416
3.416
-2。416
3。416
-0.707
1
—0。707
7
-2。414
3.414
-2。414
3.414
-0。707
1
-0。707

当k=7时, 已经稳定，于是得到

及其相应的特征向量为

应用幂法时,应注意以下两点：
（1）应用幂法时，困难在于事先不知道特征值是否满足（1)式，以及方阵A是否有n个线性无关的特征向量。克服上述困难的方法是：先用幂法进行计算,在计算过程中检查是否出现了预期的结果。如果出现了预期的结果，就得到特征值及其相应特征向量的近似值；否则,只能用其它方法来求特征值及其相应的特征向量.
(2）如果初始向量选择不当，将导致公式(3)中的系数等于零。但是,由于舍入误差的影响，经若干步迭代后，。按照基向量展开时，的系数可能不等于零。把这一向量看作初始向量,用幂法继续求向量序列,仍然会得出预期的结果,不过收敛速度较慢.如果收敛很慢，可改换初始向量。

二 原点平移法

由前面讨论知道,幂法的收敛速度取决于比值的大小.当比值接近于1时，收敛可能很慢。这时，一个补救的方法是采用原点平移法.
设矩阵
(11)
其中p为要选择的常数。
我们知道与除了对角线元素外,其它元素都相同,而A的特征值与的特征之间有关系,并且相应的特征向量相同.这样,要计算的按模最大的特征值,就是适当选择参数，使得仍然是的按模最大的特征值，且使

对应用幂法,使得在计算的按模最大的特征值的过程中得到加速，这种方法称为原点平移法.
例2 设4阶方阵A有特征值

比值,令作变换

则的特征值为

应用幂法计算的按模最大的特征值时，确定收敛速度的比值为

所以对B应用幂法时,可使幂法得到加速。虽然选择适当的p值，可以使得幂法得到加速，但由于矩阵的特征值的分布情况事先并不知道，所以在计算时,用原点平移法有一定的困难。
下面考虑当的特征值为实数时,如何选择参数，以使得用幂法计算时得到加速的方法.
设的特征值满足

则对于任意实数, 的按模最大的特征值或。
如果需要计算及时，应选择使

且确定的收敛速度的比值

当,即时，为最小。这时用幂法计算及时得到加速.
如果需要计算及时，应选择使

且确定收敛速度的比值

当即时,为最小.这时用幂法计算及时得到加速.
原点平移的加速方法,是一种矩阵变换方法.这种变换容易计算，又不破坏A的稀疏性，但参数p的选择依赖于对A的特征值的分布有大致了解.
三 反幂法

反幂法用于求矩阵A的按模最小的特征值和对应的特征向量,及其求对应于一个给定的近似特征值的特征向量。
设n阶方阵A的特征值按模的大小排列为

相应的特征向量为.则的特征值为

对应的特征向量仍然为。因此，计算矩阵A的按模最小的特征值，就是计算 的按模最大的特征值.这种把幂法用到上,就是反幂法的基本思想.
任取一个非零的初始向量，由矩阵构造向量序列
(12)
用(12)式计算向量序列时，首先要计算逆矩阵。由于计算时，一方面计算麻烦，另一方面当A为稀疏阵时,不一定是稀疏阵，所以利用进行计算会造成困难。在实际计算时,常采用解线性方程组的方法求.（12)式等价于
（13）
为了防止溢出，计算公式为
(14）
相应地取
（15）
（13）式中方程组有相同的系数矩阵A,为了节省工作量,可先对矩阵A进行三角分解
(16)
再解三角形方程组
(17）
当A是三对角方阵，或是非零元素较少且分布规律的方阵时,无论存储或计算都比较便.根据幂法的讨论，我们知道,在一定条件下,可求得的按模最大的特征值和相应的特征向量,从而得到A的按模最小的特征值和对应的特征向量，称这种方法为反幂法。反幂法也是一种迭代算法，每一步都要解一个系数矩阵相同的线性方程组.
设p为任一实数，如果矩阵可逆，则的特征值为

对应的特征向量仍为。
如果p是矩阵A的特征值的一个近似值,且

则 是矩阵的按模最大的特征值.因此，当给出特征值的一个近似值p时，可对矩阵应用反幂法，求出对应于的特征向量.反幂法迭代公式中的通过方程组

求得。
例3 用反幂法求矩阵

的对应于特征值的特征向量.
解 取解方程组

得

再解方程组

得

与的对应分量大体上成比例，所以对应于的特征向量为

第二节 雅可比方法

雅可比方法是用来计算实对称矩阵A的全部特征值及其相应特征向量的一种变换方法.在介绍雅可比方法之前，先介绍方法中需要用到的线性代数知识与平面上的旋转变换。

一 预备知识
（1） 如果n阶方阵满足

则称为正交阵。
（2) 设是阶实对称矩阵，则的特征值都是实数，并且有互相正交的个特征向量.
（3) 相似矩阵具有相同的特征值.
(4） 设是阶实对称矩阵，为阶正交阵,则也是对称矩阵.
(5) 阶正交矩阵的乘积是正交矩阵.
(6) 设是阶实对称矩阵，则必有正交矩阵，使
(1）
其中的对角线元素的是的个特征值，正交阵的第列是的对应于特征值的特征向量。
由（6）可知，对于任意的阶实对称矩阵,只要能求得一个正交阵，使(为对角阵）,则可得到A的全部特征值及其相应的特征向量,这就是雅可比方法的理论基础。

二 旋转变换
设

为二阶实对称矩阵，即。因为实对称矩阵与二次型是一一对应的,设对应的二次型为
(2)
由解析几何知识知道，方程表示在平面上的一条二次曲线.如果将坐标轴旋转一个角度，使得旋转后的坐标轴与该二次曲线的主轴重合,如图4—1所示,

则在新的坐标系中，二次曲线的方程就化成
(3）
这个变换就是
（4）
变换(4）把坐标轴进行旋转，所以称为旋转变换。其中
(5)
称为平面旋转矩阵.显然有，所以是正交矩阵。上面的变换过程即。由于
所以只要选择满足

即
（6)
(当时，可选取） 就成对角阵，这时的特征值为

相应的特征向量为

三 雅可比方法
雅可比方法的基本思想是通过一系列的由平面旋转矩阵构成的正交变换将实对称矩阵逐步化为对角阵,从而得到的全部特征值及其相应的特征向量。首先引进中的平面旋转变换.变换
（7）
记为,其中
（8）

则称为中平面内的一个平面旋转变换，称为平面内的平面旋转矩阵.容易证明具有如下简单性质：
①为正交矩阵。
②的主对角线元素中除第个与第个元素为外,其它元素均为1；非对角线元素中除第行第列元素为，第行第列元素为外,其它元素均为零.
③只改变的第行与第行元素,只改变的第列与第列元素，所以只改变的第行、第行、第列、第列元素。
设为阶实对称矩阵,为一对非对角线元素.令

则为实对称矩阵,且与有相同的特征值.通过直接计算知
(9)
当取满足关系式
(10）
时，，且
(11)
由于在正交相似变换下，矩阵元素的平方和不变，所以若用表示矩阵的对角线元素平方和，用表示的非对角线元素平方和,则由（11)式得
(12)
这说明用对作正交相似变换化为后，的对角线元素平方和比的对角线元素平方和增加了，的非对角线元素平方和比的非对角线元素平方和减少了，且将事先选定的非对角线元素消去了(即）。因此,只要我们逐次地用这种变换,就可以使得矩阵的非对角线元素平方和趋于零,也即使得矩阵逐步化为对角阵.
这里需要说明一点：并不是对矩阵的每一对非对角线非零元素进行一次这样的变换就能得到对角阵.因为在用变换消去的时候，只有第行、第行、第列、第列元素在变化，如果或为零,经变换后又往往不是零了.
雅可比方法就是逐步对矩阵进行正交相似变换，消去非对角线上的非零元素,直到将的非对角线元素化为接近于零为止,从而求得的全部特征值，把逐次的正交相似变换矩阵乘起来,便是所要求的特征向量。
雅可比方法的计算步骤归纳如下:
第一步 在矩阵的非对角线元素中选取一个非零元素 。一般说来，取绝对值最大的非对角线元素；
第二步 由公式求出，从而得平面旋转矩阵；
第三步 ,的元素由公式（9)计算.
第四步 以代替，重复第一、二、三步求出及，继续重复这一过程，直到的非对角线元素全化为充分小(即小于允许误差）时为止。
第五步 的对角线元素为的全部特征值的近似值,的第j列为对应于特征值（为的对角线上第j个元素）的特征向量.
例1 用雅可比方法求矩阵

的特征值与特征向量.
解 首先取，由于,故取,所以

=
再取由

得

所以

=
继续做下去，直到非对角线元素趋于零，进行九次变换后,得

的对角线元素就是A的特征值，即

相应的特征向量为

相应的特征值的精确值

相应的特征向量为

由此可见，雅可比方法变换九次的结果已经相当精确了。
用雅可比方法求得的结果精度都比较高,特别是求得的特征向量正交性很好,所以雅可比方法是求实对称矩阵的全部特征值及其对应特征向量的一个较好的方法.但由于上面介绍的雅可比方法,每次迭代都选取绝对值最大的非对角线元素作为消去对象，花费很多机器时间。另外当矩阵是稀疏矩阵时，进行正交相似变换后并不能保证其稀疏的性质，所以对阶数较高的矩阵不宜采用这种方法.
目前常采用一种过关雅可比方法.这种方法是选取一个单调减小而趋于零的数列（即,且)作为限值，这些限值称为”关” ，对矩阵的非对角线元素规定一个顺序（例如先行后列、自左至右的顺序).首先对限值按规定的顺序逐个检查矩阵的非对角线元素,碰到绝对值小于的元素就跳过去，否则就作变换将其化为零.重复上述过程,直到所有的非对角元素的绝对值都小于为止.再取类似处理,直到所有的非对角线元素的绝对值都小于时,迭代停止。这时的应小于给定的误差限。
实际运算中常用如下的办法取限值：对于矩阵，计算的非对角线元素平方和,任取,取

4．3多项式方法求特征值问题
4。3.1 F-L方法求多项式系数
我们知道，求n阶方阵A的特征值就是求代数方程
（4。3。1）
的根。称为A的特征多项式。上式展开为
(4。3.2）
其中为多项式的系数。
从理论上讲,求A的特征值可分为两步：
第一步 直接展开行列式｜|求出多项式；
第二步 求代数方程的根，即特征值。
对于低阶矩阵,这种方法是可行的。但对于高阶矩阵，计算量则很大,这种方法是不适用的。这里我们介绍用F—L（Faddeev—Leverrier）方法求特征方程（4。3.2）中多项式的系数。由于代数方程求根问题在第2章中已经介绍，所以本节中解决特征值问题的关键是确定矩阵A的特征多项式，所以称这种方法为多项式方法求特征值问题.
记矩阵A=的对角线元素之和为
(4。3.3）
利用递归的概念定义以下n个矩阵
(4.3。4)
可以证明,（4.3.4)式中即是所求A的特征多项式的各系数.用（4。3.4）式求矩阵的特征多项式系数的方法称为F-L方法.相应特征方程为：

（4。3。5）
而且可证矩阵A的逆矩阵可表示为
（4。3。6)
例1 求矩阵

的特征值与。
解 用F—L方法求得

所以A的特征方程为

此方程的根，即特征值为

从例1中的计算结果可知Faddeev曾经证明: 对n阶矩阵A,按(4.3。4）式计算出的总有
（4.3。7)
4。3.2 特征向量求法
当矩阵A的特征向量确定以后,将这些特征值逐个代入齐次线性程组(）x=0中，由于系数矩阵的秩小于矩阵的阶数n,因此虽然有n个方程n个未知数,但实际上是解有n个未知数的相互独立的r个方程(r<n）. 当矩阵A的所有特征值互不相同时,这样的问题中要解的齐次方程组中有n—1个独立方程,其中含有n个特征向量分量，因此特征向量分量中至少有一个需要任意假设其值,才能求出其他特征分量。
在计算机中解这样的齐次线性程组，可用高斯-若当消去法，以便把一组n个方程简化为等价的一组n—1个方程的方程组.然而,用高斯—若当消去法简化一个齐次线性程组时，方程之间不都是独立的,在消去过程中系数为零的情况较多.必需交换方程中未知数的次序,以避免主元素位置上为零的情况.因此，为了提高精度和避免零元素的可能性，我们总是用主元素措施把绝对值最大的系数放于主元素位置.
例如,假设矩阵A为

其特征方程为
=0
展开后为

故特征值分别为

下面求特征向量,将代入方程组中，得
(4.3.8）
以-5为主元素，交换上式第一与第二个方程得
（4。3.9）
用高斯—若当消去法消去-5所在列中的,并把主元素所在行调到最后，得
（4.3。10）
再以16/5为主元素,消去它所在列中的,并把主元素所在的行调到最后,得
(4.3。11）
这就是用高斯若当消去法实现把一组三个方程简化为等价的一组两个独立方程的情形.因为这个等价的方程组包含两个独立的方程，而有三个未知数，所以只要假定其中一个值，则其它两个值就可以通过两个独立方程解出.比如，令，则得到矩阵A的对应于的一个特征向量为

对另外两个特征值的对应特征向量求法与上述对的推导过程相同。
计算机中实现求解这样的齐次线性方程组的消去步骤是,用第3章讨论过的高斯-若当消去法的公式，方程组（4。3。9)的系数矩阵经过第一次消去后的矩阵B为
（4.3.12）
以矩阵为方程组（4。3.10）的系数矩阵，其中省略了有0和1元素的第一列.
在进行第二次消元之前，要应用完全主元素措施对前两行进行最大主元素选择，然后再进行必要的行或列交换.每完成一次消元过程，总省略只有0和1元素的第一列,并且计算机仅寻找矩阵的前n-k行中的最大主元素，其中k是消元过程应用的次数.对(4。3。12）式再进行一次消元过程，则得到列矩阵
(4。3。13）
此矩阵是对应于方程组(4。3。11）的系数矩阵,不过省略了含0和1元素的前两列.一般来说，最后矩阵列的数目等于矩阵的阶数和秩的差值。
由于方程组（4。3.8)有三个未知数，两个独立方程，所以计算机必须任意给定一个未知数的值,以便可以从其他两个独立方程中解出另外两个未知数.为方便,在计算机决定特征向量时,要恰当地设定任意选取的未知数的值。例如，令,由方程组（4.3.11）知道,其他两个分量的值正好能从含的非零系数项得出。为此,从计算机所存储的最终矩阵中，令最上面的0元素为-1，并把它顺次调到最下面第三行的位置上,就得到所求的特征向量。
在工程问题中，从特征方程所求出的特征值，少数情形也有相同的。一般地,当一个特征方程有k重根时,矩阵的秩可能比其阶数少1，或2,或3,…,或k，当然对应于的线性无关的特征向量的个数也就是1，或2,或3,…，或k，下面通过一个特征值对应两个线性无关特征向量的例子进一步说明计算机求特征向量的方法。
设矩阵A为

其特征方程为

展开后得

所以特征值为

为了决定的特征向量，将代入方程组（)x=0,得
（4。3.14）
应用一次高斯-若当消去法,得

（4.3.15）
写成矩阵形式，(4。3。15）式的系数矩阵为
(4.3.16)
因为方程组（4.3。15)的系数矩阵的秩为1,它比矩阵阶数少2，因此对应于有两个线性无关的特征向量,必须给两个未知数任意规定值，才能确定这两个线性无关的特征向量,由(4.3。15）式可看出,一般总是选择求一个特征向量;选择求另一个特征向量；这样有两个线性无关的特征向量
，
计算机中求两个线性无关的特征向量的办法是，在(4.3.16）式的B中,把第一列中第一个0元素用-1代替,第二列中第二个0元素也用-1代替,然后把第一、第二行顺次调到最下面一行的位置上，第三行自然就成了第一行，如此调换后矩阵的第一列和第二列就是所求的两个线性无关的特征向量。对应于的全部特征向量为

其中与是任意常数，且不同时为零。
为了说明列交换的必要性，避免主元素为零，再举一个例子,设矩阵A为

其特征方程为

特征值为

对应于的特征向量可由解下列方程组而求得
（4。3。17)
用一次高斯-若当消去法，得
（4.3。18)
若不进行列交换，则下一个消元过程只能在第一行的第二个元素与第二行的第二个元素中找最大主元素，而它们都是零，我们不得不对(4.3.17）式进行列交换,即交换未知数之间的次序，之后再进行消去过程.
对(4.3.17)式进行列交换，即把绝对值最大系数放在主元素位置，显然是第一列与第三列的交换,交换后成为
(4。3.19)
其中未知数列矩阵中与也进行了交换，这样才能保证（4.3。17)式与(4。3。19）式等价，对(4。3。19）式进行一次高斯—若当消去法，得
(4。3.20)
再进行一次消去过程，得
（4.3.21)
在计算机中计算,剩下一个最终的列矩阵
（4。3.22）
将（4.3.22)式中的列矩阵B中第一个0元素用—1代替，并随即调到最下面一行，便得到
（4。3.23)
这就是对应于方程组(4。3。19）的解,在计算机程序中应把原来进行列交换的列号次序记住，重新把(4.3。23)式中各分量排列一下，即交换第一行和第三行的元素，就得到对应于的特征向量

对应于的全部的特征向量为
k
其中k为不等于零的任意常数.

4。4 QR算 法
QR算法也是一种迭代算法，是目前计算任意实的非奇异矩阵全部特征值问题的最有效的方法之一.该方法的基础是构造矩阵序列,并对它进行QR分解。
由线性代数知识知道,若A为非奇异方阵,则A可以分解为正交矩阵Q与上三角形矩阵R的乘积，即A=QR,而且当R的对角线元素符号取定时,分解式是唯一的.
若A为奇异方阵，则零为A的特征值.任取一数p不是A的特征值，则A-pI为非奇异方阵。只要求出A-pI的特征值,就很容易求出A的特征值,所以假设A为非奇异方阵,并不妨碍讨论的一般性.
设A为非奇异方阵，令，对进行QR分解，即把分解为正交矩阵与上三角形矩阵的乘积
=
做矩阵
继续对进行QR分解
并定义
一般地，递推公式为

QR算法就是利用矩阵的QR分解，按上述递推公式构造矩阵序列。只要A为非奇异方阵，则由QR算法就完全确定。这个矩阵序列具有下列性质。
性质1 所有都相似，它们具有相同的特征值.
证明 因为

若令,则为正交阵,且有

因此与A相似,它们具有相同的特征值.
性质2 的QR分解式为

其中
证明 用归纳法.显然当k=1时,有

假设有分解式
于是

因为,所以

因为都是正交阵,所以也是正交阵,同样也是上三角形阵，从而的QR分解式为

由前面的讨论知.这说明QR算法的收敛性有正交矩阵序列的性质决定.
定理1 如果收敛于非奇异矩阵为上三角形矩阵,则存在并且是上三角形矩阵.
证明 因为收敛,故下面极限存在

由于为上三角形矩阵,所以为上三角形矩阵。又因为

所以存在，并且是上三角形矩阵。
定理2 （QR算法的收敛性)设A为n 阶实矩阵，且
1) A的特征值满足：
2) ,其中且设有三角分解式=LU（L为单位下三角阵，U为上三角阵），则由QR算法得到的矩阵序列本质上收敛于上三角形矩阵.即满足
当
当
的极限不一定存在
证明 因为,矩阵决定的收敛性。又我们利用求,然后讨论的收敛性.
由定理条件得

令
其中的(i,j）元素为

于是

由假设，当i>j时，故

设方阵X的QR分解式为

由

由知，对充分大的非奇异，它应有唯一的QR分解式，并且

于是

但上三角阵的对角线元素不一定大于零。为此，引入对角矩阵

以便保证(）的对角线元素都是正数,从而得到的QR分解式

由的QR分解式的唯一性得到

从而

由于,所以

从而
其中

于是

因为为上三角阵，为对角阵,且元素为1或-1,所以
当
当
的极限不一定存在
例 用QR算法求矩阵

的特征值.A的特征值为-1,4,1+2i，1—2i.
解 令,用施密特正交化过程将分解为

将与逆序相乘，求出

用代替A重复上面过程,计算11次得

由不难看出，矩阵A的一个特征值是4,另一个特征值是—1，其他两个特征值是方程

的根.求得为