计算方法第七章(特征值与特征向量)

1、第七章第七章 矩阵的矩阵的特征值与特征向量特征值与特征向量第一节第一节 乘幂法与反幂法乘幂法与反幂法1.1 乘幂法：乘幂法：用于求矩阵的模（绝对值）最大的特征值用于求矩阵的模（绝对值）最大的特征值。记矩阵记矩阵 A 的的特征值为：特征值为：相应的特征向量为：相应的特征向量为：任取非零任取非零向量向量 ，记，记 则有：则有： 这里，这里， 表示向量的表示向量的第第 i 个分量个分量01 122nnzccc0kkzA z12n11()lim()kikkizz( )iz12,n 11 1kkzc具体计算时，对于任意取的初始向量，按以下格式计算：具体计算时，对于任意取的初始向量，按以下格式计算：10(

2、1,2,)max() (),/kkkkkkkkyAzkzmyyzym任意给定中绝对值最大的分量1limkkm11lim/max()kkzxx11kkkkifmmorzzstop则有则有例子：求矩阵的模最大特征值及其特征向量例子：求矩阵的模最大特征值及其特征向量计算结果计算结果程序程序121241116A%用乘幂法计算矩阵模最大的特征值和相应的特征向量A=-1,2,1;2,-4,1;1,1,-6z0=1,1,1;errtel=1e-6; er=1;k=0;while ererrtel k=k+1; yk=A*z0; c,p=max(abs(yk); mk=yk(p) zk=yk/mk; er=n

3、orm(zk-z0); z0=zk;endk,mk,zk1.2 加速技术：加速技术：显然，乘幂法的收敛速度依赖显然，乘幂法的收敛速度依赖 ，如此比值接近，如此比值接近1，则收敛，则收敛速度会很慢。速度会很慢。21用用 A pI 代替代替A，进行进行乘幂法。迭代速度可能会大大加快。乘幂法。迭代速度可能会大大加快。这叫原点移位法。这叫原点移位法。2211pp12(3,4, )ipppin111,ApIAp求得的按模最大特征值的按模最大特征值埃特金加速：埃特金加速：可以证明：乘幂法线性收敛可以证明：乘幂法线性收敛 称为收敛率称为收敛率11211kkmm0112011kikizz21由于由于 线性收敛

4、于线性收敛于 ，于是可以对之进行埃特金加速，于是可以对之进行埃特金加速，以上是以上是计算特征向量计算特征向量的埃特金加速，同样可以得到关于的埃特金加速，同样可以得到关于计算特计算特征值征值的埃特金加速，的埃特金加速，kz1x22112() ()()(1,2, )()2()()kikikiikikikizzzWinzzz221122kkkkkkkm mmMmmm1kM1.3 反幂法反幂法如果如果 A 非奇异，用其逆矩阵代替非奇异，用其逆矩阵代替 A 进行乘幂进行乘幂法，称为法，称为反幂反幂法法。逆矩阵的特征值与逆矩阵的特征值与A 互为倒数。即为：互为倒数。即为：用用 A 的逆进行乘幂法，得到的逆

5、进行乘幂法，得到 迭代格式为：迭代格式为：111/1/1/nn1/n1/limnkkm110(1,2,)max() (),/kkkkkkkkyA zkzmyyzym任意给定中绝对值最大的分量为避免矩阵的求逆运算，通常也采取如下的算法：为避免矩阵的求逆运算，通常也采取如下的算法：每次迭代需要解每次迭代需要解 ，为此，可，为此，可将将 A 进行进行LR分解，分解，则每次迭代只需解两个三角方程组则每次迭代只需解两个三角方程组最后求得：最后求得：1kkAyz1kkLxzRyx1/limnkkmlimmax()nkknxzx10(1,2,)max()/kkkkkkkAyzkzmyzym任意给定反幂法的一

6、个主要应用是已知矩阵的近似特征值后，求其特征反幂法的一个主要应用是已知矩阵的近似特征值后，求其特征向量。向量。如果已求得矩阵某个特征值如果已求得矩阵某个特征值 的近似值的近似值 ，则，则于是，用反幂法可以求出于是，用反幂法可以求出 的按模最小特征值及相应的按模最小特征值及相应的特征向量。此时，迭代为：的特征向量。此时，迭代为：mm0()mmimimmAIlimmax()mkkmxzx1limkkmmm10()(1,2,)max()/mkkkkkkkAyzkzmyzym任意给定通常，初值选为：通常，初值选为：这里，矩阵这里，矩阵 L 为为 分解中的单位下三分解中的单位下三角矩阵。角矩阵。0(1,

7、1,1)zLeemAILR第二节第二节 对称矩阵的雅可比方法对称矩阵的雅可比方法两个重要的基本性质：两个重要的基本性质：（1）如）如 A 为实对称矩阵，则一定存在正交矩阵为实对称矩阵，则一定存在正交矩阵 Q ，使之相似使之相似于一个对角矩阵，而该对角矩阵的对角元正是于一个对角矩阵，而该对角矩阵的对角元正是 A 的特征值。的特征值。（2）一个矩阵左乘一个正交矩阵或右乘一个正交矩阵，其）一个矩阵左乘一个正交矩阵或右乘一个正交矩阵，其E范范数不变。数不变。1,QQQQIQAQ 22211()()nnijEEijAatrace A Atrace A Q QAQA22222()EEEEEAAQ AAQA

8、Q下面的矩阵是下面的矩阵是一个一个 n 阶正交矩阵：阶正交矩阵：1cossin( , , )1sincos1R p q( p )( q )2.1 雅可比算法雅可比算法算法的思想：算法的思想：设设 A 为对称矩阵，选出为对称矩阵，选出 A 的除对角元外的所有元素中绝对值的除对角元外的所有元素中绝对值最大的一个，然后用前一页中的正交矩阵将最大的一个，然后用前一页中的正交矩阵将。如此，产生一个新的阵，然后再重复上面的步骤，直到最后将如此，产生一个新的阵，然后再重复上面的步骤，直到最后将A 化为对角矩阵，则对角元就是所要求的特征值！化为对角矩阵，则对角元就是所要求的特征值！将上述过程数学化，首先，记将

9、上述过程数学化，首先，记 ，则，则0AA110122121,kkkkAR A RAR ARAR AR(1)1(1)11(),maxkkkkijn npqijij nAaaa 选取选取 ，使得，使得( )( )(),0kkkijn npqAaa第第 k 步步迭代矩阵的元素为：迭代矩阵的元素为：为使为使 ，必须，必须( )0kpqa( )(1)(1)( )( )(1)(1)( )( )(1)2(1)(1)2( )(1)2(1)(1)2( )cossinsincos(, )cos2sincossinsin2sincoscos(kkkkpjpjqjjpkkkkqjpjqjjqkkkkpppppqqqk

10、kkkqqpppqqqkpqqaaaaaaaajp qaaaaaaaaaa (1)(1)(1)22( )( )(1)sincos(cossin)( , )kkkkqpppqqpkkijijaaaaai jp q(1)(1)(1)22kpqkkppqqatgaa在这里，我们通常，限制在这里，我们通常，限制 ，如果，如果 ，当当 时，取时，取 ，当，当 时，时，在具体计算第在具体计算第 k 步步迭代矩阵的元素时，需要计算正弦值和余弦迭代矩阵的元素时，需要计算正弦值和余弦值，通常按如下步骤计算：值，通常按如下步骤计算：/4(1)(1)kkppqqaa(1)0kpqa/4(1)0kpqa/4 (1)(

11、1)(1)(1)(1),() 2kkkkkppqqppqqpqyaaxsign aaa2/tgx y22cos2yxy22sin2xxy1 cos2cos2sin2sin2cos实际计算中，一般预先给实际计算中，一般预先给一个计算精度一个计算精度 ，当第，当第 m 步满足步满足停止计算，这时，停止计算，这时，则对角阵的对角元为特征值近似值，矩阵则对角阵的对角元为特征值近似值，矩阵 P 的列向量为特征向的列向量为特征向量近似值。实际计算中，矩阵量近似值。实际计算中，矩阵 P 是按如下步骤计算：是按如下步骤计算：()1nmijija 1111mmmmmmmAR RR ARRRP AP 01,(1,

12、2,)kkkPIPP Rkm( )(1)(1)( )(1)(1)( )(1)cossinsincos(, )1,2,kkkipipiqkkkiqipiqkkijijppapppppjp qin 最后，雅可比方法的计算步骤可以归纳为：最后，雅可比方法的计算步骤可以归纳为：（1）确定非对角绝对值最大元位置（）确定非对角绝对值最大元位置（p，q），），并并计算计算sin和和cos的值；的值；（2）计算迭代矩阵的元素；）计算迭代矩阵的元素；（3）计算特征向量；）计算特征向量；（4）与计算精度进行比较，以决定是否终止计算，并输出特）与计算精度进行比较，以决定是否终止计算，并输出特征值和特征向量。征值和特

13、征向量。第三节第三节 QR 分解方法分解方法3.1 QR 分解分解设设 u 为为n维实单位向量，称下面矩阵为维实单位向量，称下面矩阵为Householder矩阵：矩阵：容易验证容易验证Householder矩阵为正交阵，同时又是对称阵：矩阵为正交阵，同时又是对称阵：对任意的向量对任意的向量 x 以及单位向量以及单位向量 g，存在存在Householder矩阵，使矩阵，使特别，取特别，取 g = e = ( 1 , 0 , , 0)2HIuu,HHIHH2Hxxg22xxguxxg21221(,0,0)(,0,0)niiHxxexx将将矩阵矩阵 A 记为记为于是，可以求得于是，可以求得House

14、holder矩阵，将矩阵，将 A 的第一个列向量化简。的第一个列向量化简。1212(,),(,)njjjnjAa aaaaaa 1111,s()xagign ae 11 11111111 1211 12,s(),2aeuign aaQIu uae 2111 111 111 111 112() ()2()aeaeaea1112(2)(2)1 12(,)(,)nnQ AQ a aae aa (2)(2)11212(2)(2)(2)222121(2)(2)2nnnnaaaaaAaa对矩阵对矩阵 又再重复前面的过程，即求出又再重复前面的过程，即求出Householder矩阵矩阵于是，我们记于是，我们记

15、则则1A2H2210QH(2)(2)(2)22222 12(,)(,0,0)nHaae(2)(3)(3)1121312(3)(3)(2)2232112(3)(3)212133322(3)(3)3*0nnnnnaaaaaaQ AQ Q AaaAaa如此一直下去，最后得到如此一直下去，最后得到记记 ，注意到这是一个正交矩阵，令，注意到这是一个正交矩阵，令121nnQQQ AR121nnQQQQQQAQR3.2 基本基本 QR 方法方法利用矩阵的利用矩阵的 QR 分解，立即可以得到矩阵的一系列相似矩阵分解，立即可以得到矩阵的一系列相似矩阵其中，其中， 为正交矩阵，为正交矩阵， 为上三角矩阵，为上三角

16、矩阵， 称为称为 QR 序列序列111AAQ R21122ARQQ R32233AR QQ R11kkkkkARQQ RkQkRkA11111111kkkkkkkkkkARQQAQRQAQ11111kkkkkkkAQAQAAAAA相似于相似于最后，可以证明，最后，可以证明， 的对角线下面的元素（不包括对角线）的对角线下面的元素（不包括对角线）收敛于零，由相似关系，不难推出其对角线上的元素收敛收敛于零，由相似关系，不难推出其对角线上的元素收敛到矩阵到矩阵 A 的的特征值！特征值！kA最后，要指出的是，当用最后，要指出的是，当用 QR 方法求出特征值后（准确讲，是方法求出特征值后（准确讲，是特征值

17、的近似值），我们可以用特征值的近似值），我们可以用反幂法反幂法求出求出更加精确的特征更加精确的特征值，更为重要的是可以求出相应的特征向量值，更为重要的是可以求出相应的特征向量。3.3 带原点位移的带原点位移的QR 方法方法为加速收敛，每次选取位移为加速收敛，每次选取位移 ，作，作该矩阵序列有如下性质：该矩阵序列有如下性质：（1）（2）如）如 为拟上三角，则为拟上三角，则 也为拟上三角矩阵（拟上三角也为拟上三角矩阵（拟上三角矩阵指的是次对角线下的元素为零的矩阵）矩阵指的是次对角线下的元素为零的矩阵）（3）如取位移）如取位移 为为 ，则，则 最后一行非对角元二阶收最后一行非对角元二阶收敛于零（特别

18、对于对称矩阵，能达到三阶收敛），其余次对敛于零（特别对于对称矩阵，能达到三阶收敛），其余次对角元收敛于零的速度会慢一些。角元收敛于零的速度会慢一些。1(1,2,)kkkkkkkkAs IQ RAR Qs Ikks1kkAA相似于kAkAks( )knnakA加速技术下的算法：加速技术下的算法：（1）确定计算精度）确定计算精度10E-m（2）对矩阵对矩阵 取加速因子取加速因子 进行加速进行加速（3）判断矩阵）判断矩阵 的最后一行非对角元素是否小于要求的精的最后一行非对角元素是否小于要求的精度，如果不小于，继续加速迭代，如已经小于精度，停止度，如果不小于，继续加速迭代，如已经小于精度，停止计算，并划掉矩阵的最后一行和最后一列，产生一个子矩计算，并划掉矩阵的最后一行和最后一列，产生一个子矩阵阵（4）对子矩阵重复进行上面的加速计算）对子矩阵重复进行上面的加速计算kA( )knnakA一个新的计算方案：一个新的计算方案：对于矩阵对于矩阵取变换取变换11121(1)2122211121(1)12212nnnnnnaaaaaaaAAcAaaa111UR(1)111212111(1)1 11221a