第8章 矩阵特征值计算

第八章 矩阵特征值计算1 特征值性质和估计工程实践中有许多种振动问题，如桥梁或建筑物的振动，机械机件的振动，飞机机翼的颤动等，这些问题的求解常常归纳为求矩阵的特征值问题。另外，一些稳定分析问题及相关问题也可以转化为求矩阵特征值与特征向量的问题。1.1 特征值问题及性质设矩阵 （或 ） ，特征值问题是：nARnC求 和非零向量 ，使Cx(1.1)其中 是矩阵 属于特征值 的特征向量。 的x A全体特征值组成的集合记为 。sp()求 的特征值问题(1.1) 等价于求 的特征方A程(1.2)()det)0pIA的根。因为一般不能通过有限次运算准确求解的根，所以特征值问题的数值方法只()0p能是迭代法。反之，有时为了求多项式 11()nnqaa的零点，可以把 看成矩阵123010 的特征多项式（除 因子不计） 。这是一个()nHessenberg 矩阵，可用 方法求特征值，从QR而求出代数方程 的根。q矩阵特征值和特征向量的计算问题可分为两类：一类是求矩阵 的全部特征值及其对应的A向量；另一类是求部分特征值（一个或几个、按模最大或最小）及其对应的特征向量。本章介绍部分特征值和特征向量的幂法、内积法；求实对称矩阵全部特征值的雅可比法、Given方法和 Householder 方法；求任意矩阵全部特征值的 QR 算法。在第 5 章已给出特征值的一些重要性质，下面再补充一些基本性质。定理 1 设 ，则nRA(1) 设 为 的特征值，则 为 的特AI征值；(2) 设 是 的特征值， 是一多12,n A()px项式，则矩阵 的特征值是()p。特别地， 的特征值是1(),p k。2kkn定理 (1)设 可对角化，即存在非奇nR异矩阵 使P121nA的充分必要条件是 具有几个线性无关的特征向量。(2) 如果 有 个 不同的特征值m()，则对应的特征向量 线性12,m 12,mx无关。定理 设 为对称矩阵，则nRA(1) 的特征值均为实数。(2) 有 个线性无关的特征向量。(3) 存在一个正交矩阵 使P12TnPA且 为 的特征值，而(1,2)in的列向量 为 对应于 的特征u juj向量。定理 设 为对称矩阵（其特征值依nRA次记为 ） ，则12(1) （对任何非零向量 ） 。1(,)nxnRx(2) (1.3)100(,)ma,in()nxRxRA记 ，称为矩阵 的瑞利,()A(Rayleigh)商。证明 只证(1), (2)留作习题。由于 为实对A称矩阵，可将 对应的特征向量12,n正交规范化，则有 。设12,nx (,)ijijx为 中任一向量，则有0R，1ix21/20ni于是2211(,)nniixA从而(1)成立。结论 (1)说明瑞利商必位于 和n之间。112 特征值估计与扰动定义 1 设 。令()ijnaA(1) ；,2niijjr(2) 集合 。,iiiDzrzC称复平面上以 为圆心，以 为半径的所有圆iai盘为 的格什戈林(Gershgorin) 圆盘A定理（格什戈林圆盘定理）(1) 设，则 的每一个特征值必属于下述某()ijn个圆盘之中(1。4)1,2,niiijjarn或者说， 的特征值都在复平面上 个圆盘的A并集中。(2) 如果 有 个圆盘组成一个连通的并集 ，mS且 与其余 个圆盘是严格分离，则 中恰Sn有 的 个特征值，其中重特征根按其重数重复计算。特别地，如果 的一个圆盘 是与其他圆盘AiD分离的（即孤立圆盘） ，则 中只包含 的一A个特征值。证明 只给出 (1)的证明。设 为 的任一特征值， 是相应的特征向量，即 。记xx，考虑 的第 个方程，1makiinAk即，或1kjkjx()kkjjax于是 kkjkjj ja即。kkjjar这说明， 的每一个特征值必位于 的一个圆AA盘中，并且相应的特征值 一定位于第 个圆k盘中（其中 是对应特征向量 绝对值最大的kx分量的下标） 。利用相似矩阵性质，有时可以获得 的特征值进一步的估计，即选取非奇异对角矩阵121nD作相似变换 。适当选取 1ijaA(,i可使某些圆盘半径及连通性发生变化。2,)n例 1估计矩阵 410特征值的范围。解 的个圆盘为A123:4,:,:42.DD由定理，可知 的个特征值位于个圆盘的并集中，由于 是孤立圆盘，所以 内恰好1 1包含 的一个特征值 （为实特征值） ，即。35的其他两个特征值 ， 包含在 ， 的并A2323集中。现选取对角矩阵10.9D做相似变换 14190.A的个圆盘为123:41,:,:41.89EE显然，个圆盘都是孤立圆盘，所以，每一个圆盘都包含 的一个特征值（为实特征值）A且有估计 1235,9.8.下面讨论当 有扰动时产生的特征值扰动，A即 有微小变化时特征值的敏感性。定理(Bauer-Fike 定理) 设 是的一个特征值，且nRE，则有112diag(,nPD(1.5)1()minpAPE其中 为矩阵的 范数，p=,2.证明 只要考虑 。这时 非奇异，()DI设 是 对应于 的特征向量，由xE左乘 可得()0I1P1()()DxEx11I是非零向量。上式两边取范数有1Px。1()()p而对角矩阵 的范数为I，1()(),minpAD所以有 1pPE这就得到(1.5) 式。这时总有 中的一个 取()到 值。m由定理可知 是特征值扰动1condP的放大系数，但将 对角化的相似变换矩阵A不是唯一的，所以取 的下确界P()(1.6)11()infcod()iag,)n称为特征值问题的条件数只要 不很大，A矩阵微小扰动只带来特征值的微小扰动但是难以计算，有时只对一个 ，用 代()APcond()替 特征值问题的条件数和解线性方程组时的矩阵条件数是两个不同的概念，对于一个矩阵 ，A两者可能一大一小，例如矩阵 ，10diag(,)有 ，但解线性方程组的矩阵条件数()1A0cond本章将介绍一些计算机上常用的两类方法，一类是幂法及反幂法（迭代法） ，另一类是正交相似交换的方法（变换法） 2 幂法及反幂法2.1 幂法把矩阵 的按绝对值（模）最大的特征值，A叫做 的主特征值。幂法是一种计算矩阵主特征值及对应特征向量的迭代方法，特别适用于大型稀疏矩阵设 是非亏损矩阵，即 有A个线性无关的特征向量 ，设其对应n12,nx的特征值是 ，而且满足12,n(2.1)3,现讨论求 及 的方法x设 是任一非零向量，则必存在 个不全为0v n零的数 ，使得 （并假定(1,)ian 01njavx） 。按照迭代公式10(2.2)1 (,2)kvA由初始向量 开始计算，就可得到一个向量列0，并且kv2101=()kkknnjjjavAvx(2.3)121knjk jjx由于 ，所以当 充分大时，1 (,3)j由上式可知有(2.4)11kkvAx因为 ()()k所以当 充分大时， 就是 的属于特征值 的k 1特征向量的近似向量。用 表示向量 的按模为最大的分量，max()x容易证明对任何实数 ，总有 。tmax()()tx由(2.4)式，得111max()ax()kkv因此，当 充分大时，有(2.5)1()axkkv用公式(2.5)、 (2.4)计算矩阵 的主特征值及A主特征向量的方法叫幂法。因为它使用了 的幂与初始向量 的乘积。0必须指出，使用公式(2.5) 、 (2.4)计算矩阵的主特征值及对应特征向量时，有一个巨大A的隐患，这就是：当 时， 不等于10kAv零的分量，将随 ，而无限变大。这样就k有可能发生上益；而当 时， 的各分量又k都将随着 而趋于零。为了从根本上消除这个隐患，在实用中，在迭代法的每一个步骤，都要实行规格化措施（这就是实用中的幂法） 。可归结为：定理 7 设 有 个线性无关的特征向量An，其对应的特征值为 ，且满足1,nx 1,n123从任一非零向量 ( )出发，按下列公0vua式构造向量列 及数列 ：kk(2.6)1max() ,2)kkkAvu则有(2.7)11liax()kku证明 设 ，且 ，则得010njv1201120, ,max()max()AAvvuvu22 0020011)ax(,m)ax()kkk kAAvvuu由于 ，所以02(nkjjv11211()max( )()nk kjjkkkjjaxu同理可得 。1ax (kkv算法 11. 输入 ，初始向量 ，误()ijA1(,)nu差限 ，最大迭代次数 ；N2. 置 ；,0k3. 计算 max()vu4. 若 ，输出 ，停机；否则转 5；,u5. 若 ，置 ，转 3；否则，kN1k输出失败信息，停机。评注：(1)若主特征值 是 重根，即1r12 n 矩阵 仍有 个线性无关的向量，则定理 7 中An的结论，仍然成立。事实上，前面的公式变成为 (1) 111()()kk krrrknnxuu 显然，只要 不全为零，当 充分大时，,r k就有 (1)1kkrxu于是最大特征值为 ，特征向量为 。()ikx(1)kx即对这种情况幂法仍然有效。注意：一般来说，用 个不同的初始向量，r可以得到 对应的 个不同的特征向量。1(2) 若所选择的 ，恰使 ，虽然01njavx10a在理论上幂法不收敛，但由于舍人误差的存在，幂法仍能收敛，但收敛的很慢。遇到这种情况，应另选 重新开始计算。0v例 1 用幂法求矩阵 231046A的主特征值及对应的特征向量。解 取 ，由式(2.6)式求得T0(,1)vTT1234(0,1)(2,4)6A可知 ， 。依次14T1,.5,.2u继续迭代，计算结果列表 2.1表 2.1k（归一化向量）Tumax()kkv0 0 0 1.0000 1.00001 0.5000 1.0000 0.2500 4.00002 0.5000 1.0000 0.8611 9.00003 0.5000 1.0000 0.7306 11.44004 0.5000 1.0000 0.7535 10.92245 0.5000 1.0000 0.7493 11.01406 0.5000 1.0000 0.7501 10.99277 0.5000 1.0000 0.7500 11.00048 0.5000 1.0000 0.7500 11.0000于是得主特征值的近似值 ，对应的1=.0特征向量为 。其实，T1(0.5,.,75)x该矩阵的准确特征值为 。它的效率为32。31例 用幂法计算 1.0.52A的主特征值和相应的特征向量计算过程如表2.2.表 2.2 计算结果k(规范化向量)Tumaxkv0 (1, 1, 1)1 (0.9091, 0.8182, 1) 2. 750 000 05 (0.7651, 0.6674, 1) 2. 558 791 810 (0.7494, 0.6508, 1) 2 .538 002 915 (0.7483, 0.6497, 1) 2 .536 625 616 (0.7483, 0.6497, 1) 2. 536 584 017 (0.7482, 0.6497, 1) 2. 536 559 818 (0.7482, 0.6497, 1) 2. 536 545 619 (0.7482, 0.6497, 1) 2. 536 537 420 (0.7482, 0.6497, 1) 2. 536 532 3下述结果是用位浮点数字进行运算得到的，的分量值是舍入值于是得到ku12.536及相应的特征向量 和相应T(0748,9,1)的特征向量的真值（位数字）为 1., T1(.26,6)x2.2 加速方法由前面讨论知道，应用幂法计算 的主特征A值的收敛速度主要由比值 来决定，但当21r接近于时，收敛可能很慢这时，一个补r救的办法是采用加速收敛的方法幂法的加速有许多办法。（一） 原点移位法若 的特征值为 ，则矩阵A12,n的特征值为 ，而B=pI2,npp且 特征向量相同。如果对矩阵 按前面的, B幂法计算，则有 (1)(0)1112211kkk knnxx pppaaa Ixx适当选取 ，使得 且i211 (,3)ip这样，用幂法求 的主特征值 及相应的特B1p征向量的收敛速度要比求 的快。这种加速收A敛的方法称为原点移位法。例 设 有特征值4RA15,2,34jj比值 ，做变换210.9r,1,pBI则 的特征值为 1234,0,应用幂法计算 的主特征值 的收敛速度的比1值为 2211.9p例 计算例中矩阵 的主特征值A做变换 ，取 ，则BI0.75.21对 应用幂法，计算结果如表 2.2.表 2.2 计算结果k(规范化向量)Tumaxkv0 (1, 1, 1)5 (0.7516, 0.6522, 1) 1. 791 401 16 (0.7491, 0.6511, 1) 1. 788 844 37 (0.7488, 0.6501, 1) 1 .787 330 08 (0.7484, 0.6499, 1) 1 .786 915 29 (0.7483, 0.6497, 1) 1. 786 658 710 (0.7482, 0.6497, 1) 1. 786 591 4由此得 的主特征值为 , 的主特B1.8659A征值 为1，10.752.3与例的结果比较，上述结果比例迭代次的还好若迭代次， （相1.78652应的 ） 12.368原点移位法使用简便，不足之处在于 的选p取十分困难，但在一些简单情形， 是可以估计的。如矩阵 的特征值满足A或1230n 1230n时，取 ，则有p21 (,)ipi且（因为 ）21112p因此，用原点移位求 可使收敛速度加快。评注：1）以上 的选择只有理论意义。因p为实践中很难提供 和 的值。2n2）在实践中，通常需要对特征值的分布有一大概的了解，才能粗略地估计 ，并通过计p算，不断进行修改。由圆盘定理可知，矩阵的特征值满足A1 1mini() max()n niijiiijiij jaa 三3）以上原点移位法很难形成一个自动选择的程序，所以计算机上并不实用。然而这种p加速的思想是重要的，它常在其它一些加速收敛的方法中表现出来。（二） 埃特肯(Aitken) 加速设矩阵 的特征值为 ，用幂A12n法迭代，产生序列 ，将 Aitken 加速技巧应k用于序列 ，就可获得新的序列 ，计算k k式是(2.8)212()kk k例 2 用 Aitken 加速法(2.8) 计算例 1 中矩阵的主特征值。A解 由表 2.1 可知， ，于1234,9,.4是利用公式(2.8) 得1(9)4.81.再利用 ，求 ，有234,2(1.9)91.0298.4同理得 23(0)1.4.619.已知矩阵 的准确主特征值为 11，用 AitkenA加速公式计算，迭代 3 次（用到幂法中 的值）3，即可获得较高精度的解，相当于用幂法迭代7 到 8 次。（三）幂指数加速法如果 ，则 。矩阵uA, 1,2ku的迹为 12trace()n而矩阵 的迹为kkkk若 ，则当 充分大有12n 21 111trace()kkk nA故 1trace()kA取 得矩阵序列21,mk 2482,mA 为了求出 对应的特征向量，取任意初始向量1，则0v(2.9)200mkvv就是特征向量 的近似值。1x（四）瑞利商加速法现在讨论用内积法求对称矩阵 的按模最大A特征值和特征向量。由于 对称，可设 有个线性无关单位正交特征向量 ，对n 12,nx应的特征值是 ，而且满足12,n n设 是任一非零向量，则有 （并假定0v01jav） ，于是1a(2.10)01112nnkkkkii jjaAxx利用 ，得(,)ijij211212,kkkk naav于是21212 21112(,)nkk kjk kkjaAO v(2.11)由(2.10)和(2.11)式可知，当 充分大有，1(,)kAv1max()a()kvx并且比幂法的收敛速度快。算法 21. 输入 ，初始单位向量 ，()ija1(,)nu误差限 ，最大迭代次数 ；N2. 置 ；01,k3. 计算 2TAyu4. 若 ，输出 ，停机；否则转0,u5；5. 若 ，置 ，转 3；否则，kN01,k输出失败信息，停机。2.3 反幂法反幂法用是计算矩阵按模最小特征值和特征向量的方法，也是修正特征值、求相应特征向量最有效的方法。设 为 阶非奇异矩阵，则 存在且 的An1A特征值均不为零，设 有特征值为 ，2,n相应的特征向量为 ，且12,nx 0则矩阵 的特征值为 ，相应的特1(,)i征向量也为 ，且2,nx 110因此对矩阵 应用幂法求主特征值 ，就是1An对 求按模最小的特征值。用 代替 作幂法1A计算，称为反幂法。反幂法的求解过程为：(2.12)01max()kkkuvv类似幂法的分析可得到 时，有，max()nku111knknnO用反幂法迭代一次，需要解一个方程组。实际计算时，可以事先把 作 分解，这样，每ALU迭代一次只要解两个三角方程组。例 3 用反幂法求矩阵 289347按模最小的特征值及其对应的特征向量。解 首先对 作 分解，可得ALU12894, 32.50311.0取初始向量 ，用公式(2.12)计算，T(,)uv结果列于表 8-3表 8-3k（归一化向量） max()kv0 1.0000 1.0000 1.0000 1.00001 0.4348 1.0000 -0.4783 4.56522 0.1902 1.0000 -0.8834 0.98773 0.1843 1.0000 -0.9124 0.82454 0.1831 1.0000 -0.9329 0.81345 0.1832 1.0000 -0.9130迭代 5 次，可得 ，对应的特征向10.834量为 。于是矩阵 按T3(0.82,.,9)x A模最小的特征值为 ，对应的1.2.特征向量为 。3当矩阵 有一个近似的特征值为已知时，用A反幂法可以很快地使其准确化。如果矩阵存在，设 是 的特征值 的一个近似1()pIpi值，显