数值分析第8章——矩阵特征值问题计算

1,第八章矩阵特征值问题计算,对n阶方阵A求数和非零向量x,使其满足Ax=x这样的值称为矩阵A的特征值，非零向量x称为矩阵A的与特征值相对应的一个特征向量。,2,定义1设矩阵A,BRnn，若有可逆阵P，使则称A与B相似。,定理1若矩阵A,BRnn且相似，则（1）A与B的特征值完全相同；（2）若x是B的特征向量，则Px便为A的特征向量。,8.1预备知识,3,定理2：设ARnn具有完全的特征向量系，即存在n个线性无关,其中i为A的特征值，P的各列为相应于i的特征向量。,的特征向量构成Rn的一组基底，则经相似变换可化A为,对角阵，即有可逆阵P，使,4,定理3：ARnn，1,n为A的特征值，则,（2）A的行列式值等于全体特征值之积，即,（1）A的迹数等于特征值之和，即,5,定理4,6,定理5设ARnn为对称矩阵，其特征值12n，则,（1）对任意ARn，x0，,（2）,（3）,7,定理6（Gerschgorin圆盘定理）设ARnn，则,表示以aii为中心，以半径为的复平面上的n个圆盘。,（2）如果矩阵A的m个圆盘组成的并集S（连通的）与其余,（1）A的每一个特征值必属于下述某个圆盘之中，,nm个圆盘不连接，则S内恰包含m个A的特征值。,8,9,10,11,12,13,定理7,14,15,一幂法,1基本思想,任取非零向量v0,则可唯一表示为,设n阶矩阵A的特征值,满足,且其对应有n个线性无关的特征向量x1,x2,xn，即,8.2幂法和反幂法,16,则,17,其中,由假设条件,所以当k充分大时，有,从而,所以,18,即为矩阵A的对应特征值1的近似特征向量。,且,19,则对任意非零初始向量,，下面的式子成立,定理,20,迭代公式(1)实质上是由矩阵A的乘幂Ak与非零向量v0相乘来构造向量序列xk,从而计算主特征值及其对应的主特征向量,故称这种方法为幂法。,21,则对任意非零初始向量,定理,按照下述方法构造的向量序列,则有,22,2.幂法实用计算公式,23,解取v0=(0,0,1)T,则,24,直到k=8时的计算结果见下表,从而,25,二、幂法的加速,1、原点平移法,如果是矩阵A的特征值,则对任意的实数p,矩阵A-pE的特征值为-p,且A与A-pE的特征向量相同.据此,如果要计算A的主特征值1,只要选择合适的数p,使1-p为矩阵A-pE的主特征值,且,那么,对矩阵A-pE应用幂法求其主特征值1-p,收敛速度将会加快.这种通过求A-pE的主特征值和特征向量,进而得到A的主特征值和特征向量的方法叫原点平移法。,26,且使,显然,当2-p=-(n-p),即P=(2+n)2时,上式取最小值;如果希望计算n,类似的讨论可知应选取p=(1+n-1)2。,则对任意实数p,矩阵A-pE的主特征值为1-p或n-p,要求1,则选p使,27,例2用原点平移加速法求例1中矩阵A的主特征值与其对应的特征向量。解取p=-2.5,做平移变换B=A-pE，则,对B应用幂法，仍取x0=(0,0,1)T,则,28,迭代5步的计算结果见下表,0.5,1,0.7500,13.5000,6.7500,13.5000,10.1250,5,0.5,1,0.7500,13.5007,6.7503,13.5007,10.1256,4,0.5,1,0.7507,13.5179,6.76,13.5179,10.1406,3,0.5,1,0.7545,14,7,14,10.5625,2,0.5,1,0.875,4,2,4,3.5,1,k,可得到B的主特征值113.5000特征向量v1(0.5,1.0,0.7500)T因此,A的主特征值为1=1+p11.0000,特征向量仍为v1=(0.5,1,0.7500)T。,29,30,31,设A为n阶实对称矩阵，称,为向量x的瑞利商，其中(x,x)=xTx为内积。不难证明，对实对称矩阵A，如果其特征值满足,2、瑞利商加速,由幂法公式生成的xk的瑞利商满足,由此可见，R(xk)比mk更快的收敛于1。,32,幂法的瑞利商加速迭代公式为,其中A为n阶实对称矩阵。对给定的误差限，当|mkmk-1|时，取,33,三、反幂法,反幂法是用于求非奇异矩阵A的按模最小的特征值和对应特征向量的方法.而结合原点平移法的反幂法则可以求矩阵A的任何一个具有先验了解的特征值和对应的特征向量。,设矩阵A非奇异,其特征值i(i=1,2,n),满足,其相应的特征向量x1,x2,xn线性无关,则A-1的特征值为1/i，对应的特征向量仍为xi(i=1,2,n).,34,此时,A-1的特征值满足,因此,对A-1应用幂法,可求出其主特征值1/n和特征向量xnuk,从而求得A的按模最小特征值n1/和对应的特征向量xnuk，这种方法称为反幂法。,35,为了避免求A-1,可通过解线性方程组Avk=uk-1得到yk，采用LU分解,即先对A进行LU分解A=LU,此时反幂法的迭代公式为,36,对给定的误差，当|j时,故,设方阵X的QR的分解式为,由,49,由知,对充分大的k,非奇异,它应有唯一的QR分解式,并且,于是,但上三角阵的对角线元素不一定大于零.,为此,引入对角阵,以便保证的对交线元素都是正数,50,从而得到的QR分解式,由的QR分解式的唯一性得到,从而,由于所以,51,从而,其中,于是,因为为上三角阵,为对角阵,且元素为1或-1,所以,52,的极限不一定存在,53,例1用QR算法求矩阵特征值.A的特征值为-1,4,1+2i.,54,解令用施密特正交化过程将分解为,55,将与逆序相乘,求出,用代替A重复上面过程,计算11次得,56,由不难看出,矩阵A的一个特征值是4,另一个特征值是-1,其它两个特征值是方程,的根.求得为,57,上Hessenberg化,58,59,60,61,62,63,64,用正交变换化对称矩阵为对称三对角阵,65,带原点位移的QR算法,66,67,用单步QR方法计算上Hessenberg特征值,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,隐式QR算法,88,89,90,91,92,93,定理6,设A为n阶实对称矩阵，则必存在正交矩阵,P使得,其中,为A的n个特征值。,证,设A的互不相同的特征值为,根据性质1和性质3知，,对应于特征值,94,把它们标准正交化，,特征向量组，,特征向量共有n个，,并有,其中,为A的n个特征值。,这样的,特征值的特征向量是正交的，,向量两两正交，以它们为列向量构成正交矩阵P，,又由性质2知，,A的属于不同,故这n个单位特征,95,由定理6可知，实对称矩阵的对角化问题，实质上是求正交矩阵P的问题，计算P的步骤如下：,（1）,（2）,求出齐次线性方,标准正交的特征向量，,求出实对称矩阵A的全部特征值,对于各个不同的特征值,对基础解系,这个向量组所含向量,96,（3）,（4）,则P为正交矩阵且使得,为对角阵，对角线上的元素为,相应特征向量的特征值。,取,97,例13,解,得特征值,98,99,7.2对称QR方法,对称矩阵的三对角化,100,带原点位移的QR迭代,101,102,隐式QR迭代,103,104,定义Jacobi方法是用来计算实对称矩阵A的全部特征值及其相应特征向量的一种变换方法.,Jacobi方法的基本思想是通过一系列的平面旋转矩阵构成的正交变换将对称矩阵逐步化为对角阵,从而得到A的全部特征值及其相应的特征向量.,定理(1)如果n阶方阵A满足则称A为正交阵.,7.3Jacobi方法,105,(2)设A是n阶实对称矩阵,则A的特征值都是实数,并且有互相正交的n个特征向量.,(3)相似矩阵具有相同的特征值.,(4)设A是n阶实对称矩阵,P为n阶正交阵,则B=PAP也是对称矩阵.,(5)n阶正交矩阵的乘积是正交矩阵。,设A是n阶实对称矩阵,则必有正交矩阵P，使,106,由(6)可知,对于任意的阶实对称矩阵A,只要能求得一个正交阵P,使,则可得到A的全部特征值及其相应的特征向量,这就是雅克比方法的理论基础.,其中的对角线元素是A的n个特征值，正交矩阵P的第i列是A的对应于特征值的特征向量。,下面我们详细介绍雅可比方法。首先引进中的平面旋转变换.变换,107,记为,其中,108,X=,Y=,称为平面内的平面旋转矩阵.容易得到如下性质:,则称x=为中平面内的一个旋转变换,(2)的主对角线元素中除第i个与第j个元素为外,其它元素均为1;非对角元素中除第i行第j列元素为,第j行第i列元素为外,其它元素均为零.,(1)为正交矩阵,109,(3)只改变A的第i行第j行元素,AP只改变A的第i列第j列元素,所以只改变A的第i行,第j行,第i列,第j列元素.,设A=为n阶实对称矩阵,为一非对角线元素.令则为实对称矩阵,且与A有同的特征值.通过直接计算知,110,111,当取满足关系式时,且,由于在正交相似变换下,矩阵元素的平方和不变,所以若用D(A)表示矩阵A的对角线元素的平方和,用S(A)表示A的非对角元素平方和,则由(11)式得,112,的非对角元素平方和A的非对角元素平方,且将事先选定的非对角元素消去了,这说明用对A作正交相似变换化为后,的对角线元素平方和比A的对角元素平方和增加了,和减少了,113,因此,只要我们逐次地用这种变换,就可以使得矩阵A的非对角线元素平方和趋于零,也即使得矩阵A逐步化为对角阵.,这里需要说明一点:并不是对矩阵A的每一对非对角线非零元素进行一次这样的变换就能得到对角阵.因为在用变换消去的时候,只有第i行,第j行,第i列,第j列元素在变化,如果或为零,经变换后又往往不是零了.,114,雅克比方法就是逐步对矩阵A进行正交相似变换,消去非对角线上的非零元素,直到将A的非对角线元素化为接近与零为止,从而求得A的去全特征值,把逐次的正交相似变换矩阵乘起来,便是所要求的特征向量.,Jacobi的计算步骤归纳如下:,第一步在矩阵A的非对角元素中选取一个非零元素.一般来说,取绝对值最大的非对角线元素.,115,第二步由公式求出,从而得平面旋转矩阵,第三步,的元素由公式(9)计算,第四步以代替A,重复第一,二,三步求出及继续重复这一过程,直到的非对交线元素全化为充分小(即小于允许误差)时为止.,第五步的对角化元素为A的全部特征值的近似值,的第j列为对应于特征值(为的对角线上第j个元素)的特征向量.,116,例1用雅克比方法求矩阵的特征值与特征向量,解首先i=1,j=2,由于故取所以,117,再取i=1,j=3,由,得,所以,118,继续做下去,直到非对角线元素趋于零,进行九次变换后,得,的对角线元素就是A的特征值,即,119,继续做下去,直到非对角线元素趋于零,进行九次变换后,得,的对角线元素就是A的特征值,即,120,用雅可比方法求得的结果都比较高,特别是求得的特征向量正交性很好,所以雅可比方法是求实对称矩阵的全部特征值及其对应特征向量的一个较好的方法.,但由于上面介绍的雅可比方法每次迭代都选取绝对值最大的非对