第五章矩阵的特征值与特征向量问题.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 设 用幂法求其模为最大的特征值及其相应的特征向量（精确到小数点后三位）。解 取，计算结果如表41所示.表4-1 1 1 0 1 1 1 0 12 2 2 22 1 1 13343 4-0.75 1-0.7542.53.5-2.5 3.50.714 10.71452。4283。4282。428 3.4280。708 10。70862.4163.416-2。4163。416-0.707 10。7077-2。4143.414-2。4143.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)由解析几何知识知道，方程表示在平面上的一条二次曲线.如果将坐标轴旋转一个角度，使得旋转后的坐标轴与该二次曲线的主轴重合,如图41所示,则在新的坐标系中，二次曲线的方程就化成 (3）这个变换就是 （4）变换(4）把坐标轴进行旋转，所以称为旋转变换。其中 (5)称为平面旋转矩阵.显然有，所以是正交矩阵。上面的变换过程即。由于所以只要选择满足 即 （6)(当时，可选取） 就成对角阵，这时的特征值为 相应的特征向量为 三 雅可比方法雅可比方法的基本思想是通过一系列的由平面旋转矩阵构成的正交变换将实对称矩阵逐步化为对角阵,从而得到的全部特征值及其相应的特征向量。首先引进中的平面旋转变换.变换 （7）记为,其中 （8） 则称为中平面内的一个平面旋转变换，称为平面内的平面旋转矩阵.容易证明具有如下简单性质：为正交矩阵。的主对角线元素中除第个与第个元素为外,其它元素均为1；非对角线元素中除第行第列元素为，第行第列元素为外,其它元素均为零.只改变的第行与第行元素,只改变的第列与第列元素，所以只改变的第行、第行、第列、第列元素。设为阶实对称矩阵,为一对非对角线元素.令 则为实对称矩阵,且与有相同的特征值.通过直接计算知 (9)当取满足关系式 (10）时，且 (11)由于在正交相似变换下，矩阵元素的平方和不变，所以若用表示矩阵的对角线元素平方和，用表示的非对角线元素平方和,则由（11)式得 (12) 这说明用对作正交相似变换化为后，的对角线元素平方和比的对角线元素平方和增加了，的非对角线元素平方和比的非对角线元素平方和减少了，且将事先选定的非对角线元素消去了(即）。因此,只要我们逐次地用这种变换,就可以使得矩阵的非对角线元素平方和趋于零,也即使得矩阵逐步化为对角阵.这里需要说明一点：并不是对矩阵的每一对非对角线非零元素进行一次这样的变换就能得到对角阵.因为在用变换消去的时候，只有第行、第行、第列、第列元素在变化，如果或为零,经变换后又往往不是零了.雅可比方法就是逐步对矩阵进行正交相似变换，消去非对角线上的非零元素,直到将的非对角线元素化为接近于零为止,从而求得的全部特征值，把逐次的正交相似变换矩阵乘起来,便是所要求的特征向量。雅可比方法的计算步骤归纳如下:第一步 在矩阵的非对角线元素中选取一个非零元素 。一般说来，取绝对值最大的非对角线元素；第二步 由公式求出，从而得平面旋转矩阵；第三步 ,的元素由公式（9)计算.第四步 以代替，重复第一、二、三步求出及，继续重复这一过程，直到的非对角线元素全化为充分小(即小于允许误差）时为止。第五步 的对角线元素为的全部特征值的近似值,的第j列为对应于特征值（为的对角线上第j个元素）的特征向量.例1 用雅可比方法求矩阵的特征值与特征向量.解 首先取，由于,故取,所以=再取由得所以=继续做下去，直到非对角线元素趋于零，进行九次变换后,得的对角线元素就是A的特征值，即相应的特征向量为相应的特征值的精确值 相应的特征向量为由此可见，雅可比方法变换九次的结果已经相当精确了。用雅可比方法求得的结果精度都比较高,特别是求得的特征向量正交性很好,所以雅可比方法是求实对称矩阵的全部特征值及其对应特征向量的一个较好的方法.但由于上面介绍的雅可比方法,每次迭代都选取绝对值最大的非对角线元素作为消去对象，花费很多机器时间。另外当矩阵是稀疏矩阵时，进行正交相似变换后并不能保证其稀疏的性质，所以对阶数较高的矩阵不宜采用这种方法.目前常采用一种过关雅可比方法.这种方法是选取一个单调减小而趋于零的数列（即,且)作为限值，这些限值称为”关” ，对矩阵的非对角线元素规定一个顺序（例如先行后列、自左至右的顺序).首先对限值按规定的顺序逐个检查矩阵的非对角线元素,碰到绝对值小于的元素就跳过去，否则就作变换将其化为零.重复上述过程,直到所有的非对角元素的绝对值都小于为止.再取类似处理,直到所有的非对角线元素的绝对值都小于时,迭代停止。这时的应小于给定的误差限。实际运算中常用如下的办法取限值：对于矩阵，计算的非对角线元素平方和,任取,取43多项式方法求特征值问题4。3.1 F-L方法求多项式系数 我们知道，求n阶方阵A的特征值就是求代数方程 （4。3。1）的根。称为A的特征多项式。上式展开为 (4。3.2）其中为多项式的系数。 从理论上讲,求A的特征值可分为两步： 第一步 直接展开行列式|求出多项式； 第二步 求代数方程的根，即特征值。 对于低阶矩阵,这种方法是可行的。但对于高阶矩阵，计算量则很大,这种方法是不适用的。这里我们介绍用FL（FaddeevLeverrier）方法求特征方程（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 求矩阵 的特征值与。 解 用FL方法求得 所以A的特征方程为 此方程的根，即特征值为 从例1中的计算结果可知Faddeev曾经证明: 对n阶矩阵A,按(4.3。4）式计算出的总有 （4.3。7)4。3.2 特征向量求法 当矩阵A的特征向量确定以后,将这些特征值逐个代入齐次线性程组(）x=0中，由于系数矩阵的秩小于矩阵的阶数n,因此虽然有n个方程n个未知数,但实际上是解有n个未知数的相互独立的r个方程(rj时，故 设方阵X的QR分解式为 由 由知，对充分大的非奇异，它应有唯一的QR分解式，并且 于是 但上三角阵的对角线元素不一定大于零。为此，引入对角矩阵 以便保证(）的对角线元素都是正数,从而得到的QR分解式 由的QR分解式的唯一性得到 从而 由于,所以 从而 其中 于是 因为为上三角阵，为对角阵,且元素为1或-1,所以 当 当 的极限不一定存在 例 用QR算法求矩阵 的特征值.A的特征值为-1,4,1+2i，12i. 解 令,用施密特正交化过程将分解为 将与逆序相乘，求出 用代替A重复上面过程,计算11次得 由不难看出，矩阵A的一个特征值是4,另一个特征值是1，其他两个特征值是方程 的根.求得为