第3章矩阵特征值与特征向量的计算

1、 引言引言 在科学技术的应用领域中，许多问题都归为求解一个在科学技术的应用领域中，许多问题都归为求解一个特特征系统。如动力学系统和结构系统中的振动问题，求系统的征系统。如动力学系统和结构系统中的振动问题，求系统的频频率与振型；物理学中的某些临界值的确定等等。率与振型；物理学中的某些临界值的确定等等。3.1 3.1 乘幂法及其变体乘幂法及其变体3.1.1 乘幂法乘幂法定理定理 设设A Rn n有完全特征向量系有完全特征向量系，若若 1, 2, n为为A的的n个特征值且满足个特征值且满足n21 对任取初始向量对任取初始向量x(0) Rn，对乘幂公式对乘幂公式)1()(kkAxx确定的迭代序列确定的

2、迭代序列xk，有下述结论有下述结论： （1）当当 时时，对对i = 1, 2, , n211)()1(lim kikikxx收敛速度取决于收敛速度取决于 的程度的程度，r 越小收敛越快越小收敛越快，r 1收敛慢收敛慢，112r且且x(k)（当当k充分大时充分大时）为相应于为相应于 1的特征向量的近似值的特征向量的近似值。321（2）当当 时时a）若若 1 = 2，则主特征值则主特征值 1及相应特征向量的求法同（及相应特征向量的求法同（1）；）；(2)21( )limkikkixx收敛速度取决于收敛速度取决于 的程度的程度。向量向量 、113r)(1)1(kkxxc）若）若 ，则连续迭代两次则连

3、续迭代两次，计算出计算出x(k+1)，x(k+2)，210)()1()2( kjkjkjqxpxx)(1)1(kkxx分别为主特征值分别为主特征值 1、 2相应的特征向量的近似值相应的特征向量的近似值。然后对然后对j = 1, 2, , n 解方程解方程b）若若 1 = - 2，对，对i = 1, 2, , n求出求出 、 后后，由公式由公式pq2122 pqip 2222 pqip 解出主特征值解出主特征值 1、 2。此时收敛速度取决于此时收敛速度取决于 的程度的程度。113r)(2)1(kkxx)(1)1(kkxx向量向量 、 分别为相应于分别为相应于 1， 2的特征向量的近似值的特征向量

4、的近似值。令令max(x)表示向量表示向量x分量中绝对值最大者。即如果有某分量中绝对值最大者。即如果有某i0，使，使iniixx 1max0则则 max (x) = xi对任取初始向量对任取初始向量x(0)，记，记)max()0()0()0(xxy 则则)0()1(Ayx 一般地，若已知一般地，若已知x(k)，称公式，称公式 ), 1, 0()max()()1()()()(kAyxxxykkkkk定理定理 设设A Rn n具有完全特征向量系具有完全特征向量系， 1, 2, , n为为An21则对任初始向量则对任初始向量x(0)，由规范化的乘幂法公式确定的由规范化的乘幂法公式确定的1)()max

5、(lim kkx(1)（2）y(k)为相应于主特征值为相应于主特征值 1的特征向量近似值的特征向量近似值的的n个特征值个特征值，且满足且满足向量序列向量序列y(k)，x(k)满足满足3.1.2 反幂法反幂法若若 A 有有| 1 | | 2 | | n |，则，则 A 1 有有11111 nnA 1 的主特征根的主特征根 A的绝对值最小的特征根的绝对值最小的特征根)()1(kkxAx )(1)1(kkxAx 如何计算如何计算解线性方程组解线性方程组对应同样一组特征向量。对应同样一组特征向量。设设A Rn n可逆，则无零特征值，由可逆，则无零特征值，由)0( xxAx 有有 xxA 11 规范化反

6、幂法公式为规范化反幂法公式为 ), 1, 0()max()()1()()()(kyAxxxykkkkk为节省工作量，可先对为节省工作量，可先对A进行进行LU分解，再解三角形方程分解，再解三角形方程( )(1)(0,1,)kkLxyUxxk3.1.3 乘幂法的加速（乘幂法的加速（原点位移法）原点位移法） njjkjjkkkvAxx111)1()( ，希望，希望 | 2 / 1 | 越小越好。越小越好。取取 0（常数），用矩阵（常数），用矩阵B = A - 0I 来代替来代替A进行乘幂迭代。进行乘幂迭代。0 ii (i = 1, 2, , n)iiiiiivvAvvIABv)()(000 设设 i

7、 (i = 1, 2, , n)为矩阵为矩阵B B 的特征值，则的特征值，则B与与A特征值之间特征值之间应有关系式：应有关系式：关于矩阵关于矩阵B的乘幂公式为的乘幂公式为 )0(0)0()()(xIAxBxkkk 为加快收敛速度，适当选择参数为加快收敛速度，适当选择参数 0，使使00210()maxjj n 达到最小值。达到最小值。 0101 1210()knjkjjjvv jkjnjjkvv12111 当当 i (i = 1, 2, , n)为实数，且为实数，且 1 2 n时，取时，取)(212*0n 则为则为 ( 0) 的极小值点。这时的极小值点。这时122122122*01*022212

8、12121 nnnn若已知矩阵若已知矩阵A的特征值的特征值 i的一个近似值的一个近似值 0 ，要求，要求 i和对和对记记B = A - 0I，对任取初始向量对任取初始向量x(0) Rn， ), 1, 0()max()()1()()()(kyBxxxykkkkk应的特征向量，可考虑利用原点移位加速的反幂法。应的特征向量，可考虑利用原点移位加速的反幂法。3.2 子空间迭代法子空间迭代法斯密特斯密特（Schmidt）正交化过程：正交化过程： 设设 1， 2， 3 为为R3上的三个线性无关的向量，上的三个线性无关的向量，令令 ，则，则 1为单位长度的向量，再令为单位长度的向量，再令2111222211

9、222,),( 可以验证可以验证（ 1, 2）= 0，即，即 1与与 2正交。若令正交。若令22311333),(),( 则则0),(),(2313 2333 即与即与 1, 2正交，将其单位化为正交，将其单位化为于是向量组于是向量组 1, 2, 3构成构成R3上一组标准正交基，且上一组标准正交基，且232322131221321321),(),(),(,QR其中其中Q = 1, 2, 3为正交矩阵，为正交矩阵，R是上三角阵。是上三角阵。对对n维向量空间，设维向量空间，设 1, , n为为Rn上上n个线性无关的向量，个线性无关的向量，类似有类似有11 2111 11222),( 2222 22

10、311333),(),( 2333 jjnnjnn ),(11 2nnn 232322322113122111),(),(),(),(),(),(,nnnnnn 即即QR Q为正交阵为正交阵，R 为上三角阵为上三角阵将将n个线性无关向量变换为个线性无关向量变换为n个两两正交向量的方法称为个两两正交向量的方法称为 斯密特正交化方法。斯密特正交化方法。斯密特正交化过程将可逆阵斯密特正交化过程将可逆阵A分解为正交阵与上三角阵的乘积。分解为正交阵与上三角阵的乘积。3.3 雅克比雅克比 (Jacobi) 旋转法旋转法 1预备知识预备知识1）若若B是上（或下）三角阵或对角阵，是上（或下）三角阵或对角阵，则

11、则B的主对角元素即是的主对角元素即是B的特征值的特征值。2）若矩阵若矩阵P满足满足PTP = I，则称则称P为正交矩阵为正交矩阵。显然显然PT = P-1，且且P1, P2, 是正交阵是正交阵时时，其乘积其乘积P = P1P2Pk仍为正交矩阵仍为正交矩阵。3）称矩阵称矩阵jiPij11cossin11sincos11为旋转矩阵为旋转矩阵 2雅克比方法雅克比方法设矩阵设矩阵A Rn n是对称矩阵是对称矩阵，记记A0 = A，对对A作一系列旋转相似变换作一系列旋转相似变换TCP AP其中其中A与与C的的 元素有下面的关系元素有下面的关系11c o ss i n11s i nc o s11Ti jP

12、Pij列列ij行行P是一个正交阵是一个正交阵，我们称它是我们称它是（i, j）平面上的旋转矩阵平面上的旋转矩阵 PTAP只改变只改变A的第的第i行行、j行行、i列列、j列的元素列的元素；A和和C的元素仅在第的元素仅在第i行（列）和第行（列）和第j行（列）不同行（列）不同，它们之间有如下的关系它们之间有如下的关系：cossin,cossinikkiikjkjkkjjkikccaaki jccaa22(1)22(1)cossinsin2sincossin21sin2(cos22kiiiijjijkjjiijjijijjiiijjijcaaacaaaccaaa( , )lklkcal ki j其它元

13、素不变22,1,1( )( )nnijiji ji jijCaAa引入记号22( )( )( )( )22(3.3.8)ijijCACAac易知我们选取我们选取Pij，使得使得 ，因此需使因此需使 满足满足0ijjicc2tg2(3.3.9)ijiijjaaa将将 限制在下列范围内限制在下列范围内44 如果如果 0iijjaa0ija 4 0ija 4 直接从三角函数关系式计算直接从三角函数关系式计算sin 和和cos ，记，记)1()1()1()1()1(2sinkijkjjkiikjikiiaaaxaay则则yx2tg当当 时，有下面三角恒等式：时，有下面三角恒等式：422222tg112

14、cos1cos2yxy于是于是 2221cos2yxy采用下面公式计算采用下面公式计算 sin 222cos2tgcossin22sinyxx3.3.2 Jacobi旋转法 由于一次正交相似变换AC=PTAP可将A的两个非对角元素化为零。因此可选一系列正交变换矩阵Pk，对A进行正交相似变换，直至将化为近似对角矩阵首先，在首先，在A=A1=(als(1)的非对角元素中选取绝对值最大元素，记为的非对角元素中选取绝对值最大元素，记为11111 1111,1211 1(2)(2)(2)2,=()0ijTlsijj iijPPAAP APAaaa对确定的，用公式(3.3.9)确定出 ，从而产生平面旋转矩

15、阵约化 为且的非对角元素1111(1)(1),(0)maxijijl saa然后，再在然后，再在A2的非对角元素中选取绝对值最大元素的非对角元素中选取绝对值最大元素22222,222223222211 12,12112( ,),( ,),( , ),1,2,kkijTTTkijTTTkkkijPPAijj iAP A PP P APPPPkAPP P A PPP确定出其位置，按公式(3.3.9)确定 角，从而产生平面旋转矩阵并将 中位置的元素化为零 得到如此做下去 可得到一系列平面旋转矩阵使得(1)2( )2(1)2111111(3.3.8)()()2()2()()2()22()()()1()

16、(1)(1)21()0 ()(1)kkkkkijijkijkkkkkAAaaAaAAAAn nn nAkn n 由关系式知并有P78 例 Jacobi 过关法 Jacobi方法中 ,每变换一次都要在非对角元素中扫描选取绝对值最大元素,这难免增加很大的工作量.为减少工作量和提高运行效率,下面介绍一种改进的方法,即Jacobi过关法 具体做法:0112131232(1)111212( ),( ),.nnnnijijrrAAAnaaaaaaaanAn计算矩阵 的非对角元素平方和预取阀值对 的非对角元素进行扫描对的元素进行旋转点还 将化为零 其余元素视为过关 不作变换 当所有的非对角元

17、素绝对值都小于 后,缩小阀值,如取重复前面的步骤.按此方法,阀值不断缩小 直到满足之后停止计算这种方法称为Jacobi过关法.3.4 Household方法 Householder方法是计算实对称矩阵A的部分或全部特征值及其特征向量,计算过程是:先利用正交相似变换将A约化为对称的三对角矩阵C,其次应用对分法计算C的特征值,最后计算特征向量.3.4.1 实对称矩阵的三对角化22(, ).cossin0 (, )1/,sinsjkkjjkkjjkikikkjjkikccki jccaaki jaaacoa 选取 使得某个变为零只要取, , , ,2,3,12,4,12, ,13,4,23,5,23

18、, ,21, ,2, ,1,0(, ),;,;,2,3,1,1,2,.Ti j ki j ki j kjkkjnnnn ni j iPPCPAPccki jPPPPPPPinPjiinAAC 记上述旋转矩阵为则相似变换使元素取或对用一次对 进行正交相似是变换,便可将 化为三对角矩阵Household变换 设 为单位向量,即 nRu1uuT2(3.4.2)THIuu称此为Household矩阵，简称H矩阵。容易验证 H 矩阵是对称的且正交的矩阵，即IHHHHHTT2,由于 uuuuIHuT)2(且对如何与u 直交的向量v, 都有 0Tu v vvuuIHvT2因此，对任意 ，可设 ，则其H变换为

19、nRxvuxvuHvHuHx若将 中所用与向量u正交的方向视为一个镜面，有上述公式看到，H 变换不改变向量在镜面上的投影，并将向量沿法向量的投影改变为反方向等长度的向量，因此H变换也称为镜面反射变换nR由H变换的性质不难知道，对任意非零向量 ，如果 则必存在则必存在H矩阵，使得矩阵，使得 nRyx,22yxyHx 事实上，当取 时，即可验证由（3.4.2）式所定义的矩阵满足 的要求 2yxyxuyHx 121/221211,()(,).(1),(,() ,0,0) ,TnnTrrii rHAaa aanraca aasign assarn 利用变换将 化为三对角矩阵的过程 可逐列 行 的方式进行.设为矩阵的某一列向量如想将此向量的后个分量化为零 即将 变成其中1221122112,(0,0,() ,) ,21()rrnTrrrnrra aaacacHHacwacasign as aas ssign aaw假设不全为零由于且故存在矩阵 满足若令11200() ,) ,Tn rTrrrnHIHIwwwasign as aa 则 矩阵可取为其中P83 例3.5 3.4.2 求对称三对角矩阵特征值的对分法1122210,1,2,1nnniabbabAbbabin 考虑对称的三对角矩阵这里约定00112112( ),( )1.(