数值分析原理第八章

1、第八章 矩阵特征值和特征向量的计算对于nn阶的实矩阵，线性代数理论中是通过求解特征多项式的零点而得到特征值l，然后通过求得齐次线性方程组的非零向量而得到矩阵的相应于特征值l的特征向量当矩阵阶数较高时，这种方法计算量很大，故常用数值方法近似求解特征值与特征向量目前常用的数值方法有迭代法(幂法)和变换法(Jacobi方法等)两类8.1 乘幂法与反幂法一、乘幂法乘幂法是求矩阵按模最大的特征值(主特征值)和相应的特征向量的一种迭代法设，初始向量，令(8.1)生成迭代向量序列由递推公式(8.1)，知(8.2)这表明等于用矩阵的k次幂左乘，故称此方法为乘幂法下面分析当k时，向量序列的变化规律设，为矩阵的n

2、个特征值，且满足(8.3)相应于特征值，的n个线性无关的特征向量构成向量空间上的一组基任取非零的初始向量，则可由这组特征向量线性表出(8.4)其中为线性组合系数将式(8.4)代入式(8.2)，得(8.5)由和式(8.5)，得(8.6)当时，由式(8.3)知，特征值下面针对进行讨论由式(8.6)有由于，故若，当k充分大时，此时有(8.7)上式表明，与只近似相差一个常数因子，故可取作为矩阵的相应于主特征值的近似特征向量当k充分大时，若，则有(8.8)这表明主特征值可由式(8.8)近似求得如果矩阵的特征值满足则根据式(8.6)有(8.9)则当k充分大时，由于，故有 (8.10)由于都是矩阵的特征值对

3、应的特征向量，故也是矩阵的特征值对应的特征向量由式(8.10)知，k较大时，就是与主特征值的对应的近似特征向量类似于式(8.8)，可求得主特征值的近似由于此时的特征向量子空间不是一维的，故由式(8.10)得到的近似特征向量只是该子空间的一个特征向量，对于不同的初始向量可能得到的特征向量空间中线性无关的近似特征向量对于矩阵的其它主特征值情形，如，等，同样可以用乘幂法求解，具体过程可参阅文献6以上讨论说明了乘幂法的基本原理通过上述对乘幂法过程的分析可知，乘幂法是一种迭代法，公式计算简单，便于上机实践，可以方便地用于近似求矩阵按模最大的一个(或几个)特征值及相应的特征向量需要注意的是：(1) 从理论

4、上讲，对于任意给定的初始向量，有可能使式(8.4)中的，但因舍入误差的存在，随着迭代过程的进行，等效于从的出发进行迭代(2) 在用乘幂法(8.1)进行迭代计算时，迭代向量的分量的绝对值可能会出现非常大(当)或者非常小(当)的现象，甚至出现溢出为此，实用中每进行m步就需要对迭代向量进行一次规范化，即用(其中表示向量的按模最大的分量)代替继续迭代由于特征向量允许相差一个常数因子，故按前面乘幂法的理论依然得到正确的近似特征向量在每次规范化后，用乘幂计算前后两个向量的分量的比值作为主特征值的近似，这种规范化并不影响主特征值的近似计算。规范化的乘幂法避免了溢出的可能性，至于m取多少，取决于实际情形，可以

5、取m=5或m=1等等下面给出乘幂法的具体算法算法8.1 (乘幂法)输入 矩阵、误差限、任意初始向量、m；输出 主特征值的近似值及相应的近似特征向量；Step 1 利用式(8.1)求，并由式(8.8)求，置k=1；Step 2 利用式(8.1)求，并利用式(8.8)求；Step 3 判断是否满足？若满足，则取，结束；否则，转向Step4；Step 4 置k=k+1，判断mod(k1,m)=0？若满足，取；转向Step2例8.1 用乘幂法计算矩阵的主特征值及相应特征向量，要求解 取初始向量，用乘幂法公式进行计算，且取，每迭代5步进行一次规范化，计算结果见表8.1表8.1 乘幂法的计算结果k0123

6、453637(1.0000000 1.0000000 1.0000000)(-6.0000000 2.0000000 8.0000000)(102.0000000 -32.0000000 34.0000000)(-1218.0000000 206.0000000 608.0000000)(17058.0000000 -4664.0000000 190.0000000)(-218118.0000000 46130.0000000 61832.0000000)(-13.2201800 3.1081369 2.2688627)(174.7731588 -41.0901285 -28.9947753)

7、-6.0000000-17.0000000-11.9411765-14.0049261-12.7868449-13.2201800-13.2201800取-13.2201800作为主特征值的近似值，相应于的特征向量的近似值取为(174.7731588, -41.0901285, -28.9947753)T二、原点平移法乘幂法的收敛速度取决于或的程度，r1时收敛速度较快，r1时收敛速度较慢对收敛较慢的乘幂法，可以采用原点平移法加快其收敛速度设矩阵，这里p是可以选取的参数当的特征值为时，的特征值，且与有相同的特征向量若的主特征值为，则要选择适当的参数p，使其满足是的主特征值，即 ()，且对矩阵应用

8、乘幂法，可以使得计算的主特征值的过程得到加速这种方法通常称为原点平移法对于参数p的选择依赖于对矩阵的特征值分布的大致了解通常可以用Gerschgorin圆盘定理得到矩阵的特征值分布情况定理8.1(Gerschgorin圆盘定理)19 设为n阶实矩阵，则(1) 的每一个特征值必定属于下述n个闭圆盘(称为Gerschgorin圆) (8.11)的并集；(2) 由矩阵的所有Gerschgorin圆组成的连通部分中任取一个，如果它是由k个Gerschgorin圆构成，则在这个连通部分中有且仅有的k个特征值(Gerschgorin圆相重时重复计数，特征值相同时也重复计算)求得矩阵的主特征值后，可得矩阵的

9、主特征值，同时得到对应的特征向量的近似例8.2 对例8.1的矩阵，取p=4.6，用原点平移法求其主特征值及相应的特征向量解 对应用乘幂法算法8.1，每迭代1步进行一次规范化，计算结果见表8.2由表8.2可知，的主特征值= -17.8201809+4.6= -13.2201809，相应的特征向量的近似值取为(1.0000000, -0.2351055, -0.1716216)T表8.2 原点平移法的计算结果k01231112(1.0000000 1.0000000 1.0000000)(-10.6000000 -2.6000000 3.4000000)(-16.8264160 2.7584906

10、 1.7396227)(-17.4019737 3.7969499 3.0797486)(-17.8201809 4.1896219 3.0583200)(-17.8201809 4.1896216 3.0583203)(1.0000000 1.0000000 1.0000000)(1.0000000 0.2452830 -0.3207547)(1.0000000 -0.1639281 -0.1033864)(1.0000000 -0.2181908 -0.1769770)(1.0000000 -0.2351055 -0.1716212)(1.0000000 -0.2351055 -0.171

11、6216)-10.6000000-16.8264160-17.4019737-17.8201809-17.8201809事实上，如果对于矩阵的特征值能够分离的很清楚，可以利用原点平移法求得矩阵的所有特征值及其相应的特征向量但需要说明的是，虽然常常能够选择合适的p值使乘幂法得到加速，但设计一个自动选择适当参数的过程非常困难原点平移法的价值不在于直接使用它使迭代过程加速，而在于把原点平移法与其它方法(如反幂法等)结合使用以获得更好的效果三、反幂法反幂法又称逆幂法，它是近似求矩阵按模最小的特征值及其相应特征向量的一种方法设，且无零特征值，则若l为矩阵的特征值，则必为矩阵的特征值，且与具有相同的特征向

12、量如果的n个特征值满足则的n个特征值满足因此，若用做乘幂矩阵，由乘幂迭代格式 (8.12)便可求出的按模最大特征值，进而得到矩阵的按模最小的特征值因此，对任取初始向量，称式(8.12)为求矩阵按模最小特征值的反幂法在应用式(8.12)计算时，高效的算法并不是先计算的逆矩阵，而是通过求解如下线性方程组求得 (8.13)由于在迭代过程中求解的线性方程组(8.13)的系数矩阵不变，故可以利用矩阵的三角分解，则每次迭代只需解两个三角形方程组 (8.14)下面给出反幂法的具体算法算法8.2 (反幂法)输入 输入矩阵、误差限、任意初始向量；输出 按模最小的特征值的近似值及相应的近似特征向量；Step 1

13、利用式(8.14)求，并由式(8.8)求，置k=1；Step 2 利用式(8.14)求，并利用式(8.8)求；Step 3 判断是否满足？若满足，则取，结束；否则，置k=k+1，转向Step2根据反幂法与乘幂法的上述关系，若p是的相对分离较好的近似值()，可以结合原点平移法与反幂法来求得的更精确的特征值与相应的特征向量的近似此时只需利用反幂法算法8.2求解的按摸最小的特征值及其相应的特征向量即可8.2 Jacobi方法Jacobi方法是近似求实对称矩阵全部特征值及相应特征向量的一种变换方法其基本思想是将对称矩阵经一系列正交相似变换约化为一个近似对角阵，从而该对角阵的对角元就是的近似特征值，由各

14、个正交变换阵的乘积可得对应的特征向量在本节中，将用到下列线性代数知识：(1) 若矩阵与相似，即存在可逆阵使，且与具有完全相同的特征值(2) 若矩阵满足，则称为正交矩阵显然，且当，是正交矩阵时，其乘积仍为正交矩阵(3) 实对称矩阵的特征值均为实数(4) 对任何实对称矩阵，总存在正交阵，使得=，其中为的特征值，的各列为相应的特征向量(5) 设为实对称矩阵，为正交阵，则(6) 称矩阵为n阶的Givens旋转矩阵，它是在单位阵的p行、q行和p列、q列的四个交叉位置上分别置上cosq，sinq，-sinq和cosq而成的容易验证旋转矩阵是正交矩阵，即，用它作相似变换时十分方便Jacobi方法就是用这种旋

15、转矩阵对实对称阵作一系列的旋转相似变换，从而将约化为近似对角阵一、Jacobi方法首先考虑的二阶实对称矩阵的对角化取二阶矩阵对矩阵做旋转变换，即其中 (8.15)由式(8.15)的最后一式知，要使的相似矩阵成为对角阵，只须适当选取q，使即时取 (8.16)当时，取确定q后，旋转矩阵则随之确定此时得到矩阵的特征值为，它们对应的特征向量分别是的各列，即对应于的特征向量分别是可以看出，对二阶实对称矩阵，用适当的正交相似变换，一次即可把化为对角阵当为n阶实对称矩阵时，需使用n阶Givens旋转矩阵容易验证，对于n阶Givens旋转矩阵来说，只改变的第p行与第q行元素，只改变的第p列与第q列元素，只改变

16、的第p行、第q行、第p列与第q列元素Jacobi方法是通过一系列旋转相似变换逐渐将实对称矩阵化为对角阵的过程，即 (8.17)恰当地选取每个旋转矩阵，就可以使趋于对角化设，其中p和q分别指向矩阵的严格上三角阵中绝对值最大的元素的行号和列号变换过程中产生的也是实对称阵与的差别仅在于第p、q行与第p、q列的元素由矩阵乘法可得 (8.18) (8.19) (8.20)由式(8.18)(8.20)知，再结合正交矩阵的性质，可得若，选取q使，只需q满足 (8.21)或 (8.22)此时则有 (8.23)引入记号 (8.24)式中表示矩阵的对角元的平方和，表示矩阵的非对角元的平方和由于，故有 (8.25)

17、这说明，只要，则按上述方法构造的Givens旋转矩阵对变换后就会使对角线元素的平方和增加，而非对角线元素的平方和减少由于的对称性，实际上只要计算的上三角元素，而下三角元素由对称性获得，这样即节省了计算量，又避免了舍入误差对对称性的影响需要指出的是，若在某一步已有，则可能又变为非零元素因此，并不能保证通过有限次旋转相似变换将矩阵化为对角阵Jacobi方法中每一步都有两个主要步骤(1) Givens旋转矩阵的确定选取Givens旋转矩阵，使得，则需使q满足式(8.21)或式(8.22)考虑到舍入误差的影响，当时，取 (8.26)当时，记 (8.27)由式(8.21)及式(8.27)有 (8.28)

18、 (8.29)为避免相近数相减(当非常大时)，取 (8.30)进而有 (8.31)(2) 特征向量的计算由式(8.17)知，若记，则有下述递推关系 (8.32)如果趋于对角阵，则的各列就是近似特征向量的计算可由式(8.32)得到，具体地其元素间的关系为 (8.33)下面给出Jacobi方法求矩阵特征问题的具体算法算法8.3 (Jacobi算法)输入 输入矩阵、误差限；输出 矩阵的所有特征值及其相应的特征向量；Step 1 置，k=0；Step 2 选主元，即确定，使；Step 3 按式(8.27)或式(8.31)计算；Step 4 按式(8.33)计算新正交阵的元素；Step 5 按式(8.1

19、8)(8.20)计算Ak+1的元素；Step 6 计算；若，则停止计算，否则，置k=k+1，转向Step2下面给出Jacobi方法收敛性的结论定理8.2 设为实对称矩阵，则由出发，Jacobi方法所产生的矩阵序列收敛于对角阵，且就是实对称矩阵的全部特征值证明 由式(8.26)可知由知从而有进而有递推得到当时，显然有，即非对角元的平方和趋于零，趋于对角阵，Jacobi方法收敛例8.3 用Jacobi方法求例8.1中实对称矩阵的全部特征值与特征向量，取解 按Jacobi算法8.3编程上机计算，结果见表8.3从表8-3可以看出，的特征值分别为，对应的特征向量分别为表8.3 Jacobi方法的计算结果k0AI1267，二、Jacobi法的变形使用Jacobi方法进行计算时，每次先要在非对角元中循环扫描以挑选主元，这要花费较多的机时实用中，常常对Jacobi方法进行修正在开始扫描之前，首先确定一个阀值，在实对称矩阵的严格