第九章_矩阵特征值问题的数值方法

1、第9章 矩阵特征值问题的数值方法 9.1 特征值与特征向量9.1 特征值与特征向量 设A是n阶矩阵，x是非零列向量. 如果有数存在，满足 ， (1)那么，称x是矩阵A关于特征值的特征向量. Axx如果把(1)式右端写为 ，那么(1)式又可写为:x()0IA x| 0IA即1110( ) |.nnnfIAaaa记它是关于参数的n次多项式，称为矩阵A的特征多项式， 其中a0=(-1)nA. (2) 显然，当是A的一个特征值时，它必然是 的根. 反之，如果是 的根，那么齐次方程组(2)有非零解向量x，使(1)式成立. 从而，是A的一个特征值. A的特征值也称为A的特征根. ( )0f( )0f矩阵特

2、征值和特征向量有如下主要性质： 定理阶矩阵A是降秩矩阵的充分必要 条件是A有零特征值. 定理设矩阵A与矩阵B相似，那么它们 有相同的特征值. 定理阶矩阵A与AT有相同的特征值. 定理设ij是n阶矩阵A的两个互异特 征值，x、y分别是其相应的右特征向 量和左特征向量，那么，xTy=0 . 9.2 Hermite矩阵特征值问题 设A为n阶矩阵，其共轭转置矩阵记为AH. 如果A=AH，那么，A称为Hermite矩阵.Hermite矩阵的有关性质 设 是Hermite矩阵A的n个特征值. 有以下性质: 全是实数.12,.,n12,.,n 有相应的n个线性无关的特征向量，它们可以化为一组标准酉交的特征向

3、量组 ,即 12,.,n 12,.,nu uuHiju uij 是酉空间中的一组标准酉交基. 12,.,nu uu 记U=( ),它是一个酉阵，即UHU=UUH=I，那么即A与以 为对角元的对角阵相似.12,.,nu uu1HnUAUD12,.,n A为正定矩阵的充分必要条件是 全为正数. 12,.,n定理设 是Hermite矩阵A的n个特征值，那么证:12,.,n21maxiinA 21niFiA2222()()( ( )HAA AAA由21maxiinA 因此2221()()nHiFiAtr A Atr A又由21niFiA得 设x是一个非零向量，A是Hermite矩阵，称 为矩阵A关于向

4、量x的Rayleigh商，记为R(x). HHxAxxx定理如果A的n个特征值为 其相应的标准酉交的特征向量为 那么有 12.n12,.,nu uu1( )nR x定理设A是Hermite矩阵 ，那么100min( )min( )kn kkkx Cxx CxR xR x 且且或极值定理 定理极值定理) 设Hermite矩阵的n个特征值为 ，其相应的标准酉交特征向量为 . 用Ck表示酉空间Cn中任意的k维子空间，那么12.n12,.,nu uu1100max min( )maxmin( )kkn kn kkCx CxkCx CxR xR x 且且或Hermite矩阵特征值问题的性态 矩阵特征值问

5、题与求解线性方程组问题一样，都存在当矩阵A的原始 数据有小变化(小扰动)时，引起特征值问题的变化有大有小的问题，如果引起的变化小，称 该特征值问题是良态的. 反之，称为病态的. 矩阵特征值问题的性态是很复杂的，通常分别就单个特征值或整体特征值给出状态数进行分 析. 对于Hermite矩阵，由于其特征值问题的特殊性质，其特征值都是良态的.下面先证明Hermite矩阵特征值的扰动定理. 定理设矩阵A，E，A+E都是n阶Hermite矩阵，其特征值分别为 那么,证 设矩阵A关于特征值1，2，n 的标准酉交特征向量为u1，u2，un， 是由ui，ui+1，un生成的n-i+1维子空间. 对 中任意非零

6、向量x,由极值定理,有12n12n12n1inii1n iC 1n iC 111000()maxmaxmaxn in in iHiHx CxHHHHx Cxx CxxAE xx xx Axx Exx xx x 且且且由定理，又由定理，对任意x0，有从而有另一方面, A=(A+E)-E. 记 为矩阵-E的特征值，那么，重复上面的过程,可得从而有10maxn iHiHx Cxx Axx x 且110maxn iHnHx CxxExxx 且1ii12n1in i 1iiiin定理通常又称为Hermite矩阵特征值的扰动定理 定理设矩阵A和A=A+E都是n阶Hermite矩 阵，其特征值分别为 和 ，

7、那么 这个定理表明，扰动矩阵E使A的特征值的变化不会超过 E2. 一般E2小，因此，Hermite矩阵特征值是良态的. 12n12n222iiEE9.3 Jacobi方法理论上，实对称矩阵A正交相似于以A的特征值为对角元 的 对角阵. 问题是如何构造这样的正交矩阵呢? Jacobi方法就是通过构造特殊的正交矩阵 序列，通过相似变换使A的非对角线元素逐次零化来实现对角化的. 平面旋转矩阵与相似约化先看一个简单的例子. 设 是二阶实对称矩阵，即a21=a12，其特征值为1，2. 令 使得 记 容易验证BT=B, 且11122122aaAaacossinsincosR11TR AR11122122T

8、bbBR ARbb22111112222212212211122222111222cos2sincossin()sincos(cossin)sin2sincoscosbaaabbaaabaaa解之得:当 时当 时 并规定1122aa12112221arctan2aaa1122aa21arctan2pqppqqaaa4经典的Jacobi方法 设A是实对称矩阵，记A1=A.Jacobi方法的基本思想是用迭代格式 Ak+1=QTkAkQk , k=1,2, 构造一个相似矩阵序列，使Ak收敛于一个对角阵. 其中 Qk为平面旋转矩阵，其旋转角k由使Ak的绝对值 最大元a(k)pq=a(k)qp=0 或按

9、列依次使A的非对角元 零化来确定. 定理设A是n阶实对称矩阵，那么由Jacobi方法产生的相似矩阵序列Ak的非对角元收敛于0. 也就是说，Ak收敛于以A的特征值为对角元的对角阵. 记 其中Ek是Ak除主对角元外的矩阵.由平面旋转矩阵的性质 中,对于 ,有因此,( )()kkiikAdiag aE1TkpqkpqAR A R,ip q(1)2(1)2( )2( )2()()()()kkkkipiqipiqaaaa22( )212()kkkpqFFEEa又由假设,因此,这样,便有从而,当22( )( )( )( )1,1max(1)nkkkkpqijpqpqi j niji jijaaan na且

10、且且22( )(1)kkpqFEn na222112222(1)(1)kkkFFFEEEnnnn1,0kFkE 实用的Jacobi方法 循环Jacobi方法必须一次又一次扫描，才能使Ak收敛于对角阵 ，计算量很大. 在实际计算中，往往用一些特殊方法来控制扫描次数，减少计算量. 下面介 绍一种应用最为广泛的特殊循环Jacobi方法阈Jacobi方法. 阈Jacobi方法首先确定一个阈值，在对非对角元零化的一次扫描中，只对其中绝对值 超过阈值的非对角元进行零化. 当所有非对角元素的绝对值都不超过阈值后，将阈值减少， 再重复下一轮扫描，直至阈值充分小为止. 减少阈值的方法通常是先固定一个正整数Mn，

11、扫描一次后，让 . 而阈值的下界是根据实际问题的精度要求选定的. M用Jacobi方法计算特征向量 假定经过k次迭代得到Ak+1=RTkRT1AR1Rk，(15) 这时Ak+1是满足精度要求的一个近似的对角阵. 如果记Qk=R1R2Rk=Qk-1Rk，(16) 那么，Qk是一个正交矩阵，且(15)式又可表示为Ak+1=QTkAQk.当Ak+1的非对角元素充分小，Qk的第 j列qj可以看成是近似特征值a(k+1)jj相应的特征向量了. 在实际计算中，可以按(16)式在迭代过程中形成Qk，把Qk看成是Qk-1右乘一个平面旋转矩阵得到. 不妨记 Q0=I，Qk的元素按下式计算：( )(1)(1)(

12、)(1)(1)( )(1)cossin,sincos, ,1,2,kkkipipkipkkkkiqipkiqkkkijijqqqqqqqqjp qin 9.4 对分法 理论上，一个实对称矩阵正交相似于一个以其特征值为对角元的 对角阵. 但是，经典的结果告诉我们，一个大于4次的多项式方程不可能用有限次四则运算 求根. 因此，我们不可能期望只用有限次相似变换将一个实对称矩阵约化为一个对角阵.下面先介绍将一个实对称矩阵相似约化为实对称三对角矩阵的方法，再讨论求其特征值的对 分法. 相似约化为实对称三对角矩阵 将一个实对称矩阵正交相似约化为一个实对称三对角矩阵的算法，可归纳如下： 记A(1)=A,对k

13、=1,2,n-2 按(4)式、(5)式和(8)式计算 ； 按(9)(12)式，计算A(k+1). ,kkku和( )1,( )2,( )(4)kkkkkkkkknkaaua ( )( )21,1()() (5)nkkkkkiki ksignaa ( )1,()(8)kkkkkka1( )1(1)( ),(9),(10)1,(11)2(),(12)kkkkTkkkkkkkkkkTTkkkkgA uu gyguAAu yy u序列的性质 设实对称三对角矩阵为 其中i0 (i=1,2,，n-1) 其特征矩阵为T-I. 记T-I的第i阶主子式为 111221nnaTaa111211( )iiiiaap

14、a 这是关于的i次多项式，当i=n时， pn()=T-I是矩阵T的特征多项式. 令p0()1，则有p1()=1-，pi()=(i-)pi-1()-2i-1pi-2()，i=2,3，n.(15) 多项式序列pi() (i=0,1,，n)称为Sturm序列 111211( )iiiiaapa定理pi() (i=1,2,，n)的根都是 实根. 证 由(14)式，pi()是i阶实对称矩阵的特征多项式，因此，pi() (i=1,2,，n)的根全是实根. 定理定理设是pi()的一个根，那么 pi-1()pi+1()0，即相邻的两个多项式无公共根； pi-1()pi+1()0，即pi-1()与pi+1( )

15、反号. lim( )0,iplim( )0,lim( )0,iipp当i为奇数当i为偶数定理pi()的根都是单根，并且将pi+1()的根严格隔离. 同号数和它的应用 定义1 设p0()1，pi()(i=1,2,，n) 是一个Sturm序列，称相邻的两个数中符号一致的数目为同号数，记为ai(). 若某个pi()=0，规定与pi-1()反号. 定理设两个实数x3n，可以用乘幂法计算2及其相应的 特征向量. 在计算1和v后，按(15)式形成n-1阶矩阵B的计算过程称为收缩方法. 9.6 反幂法 反幂法可以求一个非奇异矩阵A的逆矩阵A-1的按模最小的特征值及相应的特征向量，又可以求A的一个近似特征值相

16、应的特征向量. 求按模最小特征值及相应特征向量的反幂法,又称为反迭代法.1(1)1,0,1,kkkkkTkkkTkkxzxLUxzz xkz z求近似特征值的特征向量的反幂法 先对矩阵 进行LU分解，记 那么， (7) 下面介绍一种选取特殊的初始向量x0的反幂法半迭代法. lIAlIALU1,0,1,kkkkkkkxzxLyzUxy k 假设 ,选取初始向量x0满足x0=1，这时z0=x0.对照(7)式中的第二个式子.可把z0看成满足Le=z0.(8) 这里，e=(1,1，1)T，而z0的各个分量的取值多少是无关重要的.这样，在第一个迭代步的计算中，只需求解(7)式中的上三角方程组Ux1=e.

17、 “半迭代法”的命名也由此而得.lIA9.7 QR方法 定理设A是n阶矩阵，其n个 特征值为 .那么存在一个酉矩阵U，使UAU是以为 对角元的上三角矩阵. 12,n 12,n 两个基本定理定理设A是n阶实矩 阵，那么，存在一个正交矩阵Q，使QAQ为一个准上三角矩阵，它的对角元是A的一个特征值，对角元上的二阶块矩阵的两个特征值是A的一对共轭复特征值. 相似约化为上Hessenberg矩阵 对一般n阶矩阵，QR算法的每一个迭代步需要O(n)次乘法运算.如果矩 阵阶数稍大，这个算法几乎没有实际的应用价值.通常采用的方法是先将矩阵相似约化为上Hessenberg形式的矩阵，在此基础上应用QR迭代.这时

18、，一个QR迭代步的乘法运算次数只需O(n）次. 所谓上Hessenberg矩阵是指一个n阶矩阵A，如果当ij+1时，aij=0，称A为上Hessenberg矩阵.例如：一个5阶的上Hessenberg矩阵具有如下的形式： *0*00*000*A下面介绍QR方法时，都假设矩阵A是一个上Hessenberg矩阵. 算法 设A是n阶矩阵且有QR分解AQR，(2) 这里，Q是酉矩阵，R是上三角矩阵.如果A是满秩并规定R有正对角元，这个分解是惟一的. 一、QR算法的基本思想 记AA且有AQ 1R1.将等号右边两个矩阵因子的次序交换，得ARQ，且 ，(3) 即AA.不难证明:即Ak+1AkA，矩阵序列Ak

19、有相同的特征值. 记12111AQ AQ111111kkkkkkkAQ A QQQ A QQ12kkQQQQ12kkRR RR容易得到 是Ak的一个QR分解kkkAQ R 如果A是一个满秩的上Hessenberg矩阵，可以证明，经过一个QR迭代步得到的A2QA1Q仍然是上Hessenberg矩阵. 因为上Hessenberg矩阵次对角线以下的元素全为0，因此，只要证明，当k时，由迭 代格式(4)产生的矩阵Ak的次对角元趋向于零就可以了. 二、 QR算法的收敛性 定理设n阶矩阵A的n个特征值满足|n|0，其相应的n个线性无关特征向量为x1，x2，xn. 记X(x1,x2,xn), Y= X.如果Y存在LU分解，那么，由(4) 式产生的矩阵Ak基本收敛于上三角矩阵R.这里，基本收敛的含义指Ak的元素中除对角线以下的元素趋于零外，可以不收敛于R的元素. 三、 QR算