矩阵特征值与特征向量的计算yjs课件

第第8章章 矩阵特征值和特征向量的计算矩阵特征值和特征向量的计算 很多工程计算中，会遇到特征值和特征向量的计算，如：机械、结构或电磁振动中的固有值问题；物理学中的各种临界值等。这些特征值的计算往往意义重大。求解线性方程组的迭代法，重要一点是判断迭代法的收敛性；判断方法之一就是看迭代矩阵的特征值的模是否都小于1。1PA()是的高次的多项式，它的求根是很困难的。设法通过数值方法是求它的根。通常对某个特征值，可以用些针对性的方法来求其近似值。若要求所有的特征值，则可以对A做一系列的相似变换，“收敛”到对角阵或上（下）三角阵，从而求得所有特征值的近似。n阶方阵A的特征值是特征方程 PA()=det(A-E)=0 的根.A的特征向量是齐次线性方程组 (A-E)x=0 的非零解.特征根和特征向量的定义（复习）2定理定理1：A R n n，1,n为为A的特征值的特征值，则则（2）A的行列式值等于全体特征值之积，即的行列式值等于全体特征值之积，即（1）A的迹数等于特征值之和，即的迹数等于特征值之和，即特征根和特征向量的基本结论。定理定理2 设为A R n n的特征值且Ax=x，其中x不为0,则（1）c为为cA的特征值（的特征值（c为常数且不为为常数且不为0);（2）-p为为A-pI 的特征值，即的特征值，即(A-pI)x=(-p)x;（3 3）k为为Ak的特征值；的特征值；（4 4）设设A A为非奇异阵，那么为非奇异阵，那么 且且 为为 特征值，即特征值，即3定义定义 设矩阵设矩阵A,B R n n，若有可逆阵若有可逆阵P，使使 则称则称A与与B相似相似。定理定理 若矩阵若矩阵A,B R n n且相似且相似，则则（1）A与与B的特征值完全相同的特征值完全相同；（2）若若x是是B的特征向量的特征向量，则则Px便为便为A的特征向量的特征向量。相似矩阵及定义其性质48.1 幂法和反幂法幂法和反幂法 8.1.1 幂法幂法 幂法是用来求矩阵A按模最大最大的特征值和相应的特征向量的方法.也称为主特征值主特征值和主特征向量主特征向量。设A是单构矩阵,即A有n个线性无关的特征向量.A的n个特征值为|1 2 n 对应的特征向量为1,2,n 线性无关.我们要求1 和1.幂法的基本思想是取初始非零向量x0Rn,作迭代 xk+1=Axk=Ak+1x0,k=0,1,2,产生迭代序列xk.由于1,2,n 线性无关,从而有 x0 =11+22+nn (8.3)5xk =Akx0=11k1+22k2+nnkn设|12n,这时,上式可写成若10,则对充分大的k有 因而有 从而特征向量1 xk.乘幂法的收敛速度取决于|2/1|的大小.故有 8.1.1 幂法幂法6因此，常把每一步计算的迭代向量xk规范化。对非零向量x,用max(x)表示x的按绝对值最大的分量,称向量y=x/max(x)为向量x的规范化向量.例如,设向量x=(2,1,-5,-1)T,则max(x)=-5,y=(-0.4,-0.2,1,0.2)T.可见规范化向量y总满足y=1.幂法的规范化计算公式为:任取初始向量x0=y0 0 0,计算可得实际计算时，考虑到当11时，xk的非零向量趋于无穷；当11时,xk趋于零；导致计算机会出现上溢或下溢。8.1.1 幂法幂法7所以其收敛速度由比值|2/1|来确定.又由于所以因此,当k充分大时可取:1 mk,1 xk.8.1.1 幂法幂法8算法8.1 幂法 程序见p174。(1)输入矩阵A，非零初始向量y0，最大迭代次数N，精度，置k:=0，u=0；(2)计算 mk=max(yk)；(3)计算(4)若|mk-u|n|0对应的特征向量为1,2,n,则有A-1的特征值为对应的特征向量为n,n-1,1.要想求n和n只需对A-1应用乘幂法,任取初始向量x0=y00,作14也可将上式改写成式(8.8)称为反幂法.显然有每一步求xk需要求解线性方程组,可采用LU分解法求解.8.1.3 反幂法反幂法15(8.9)8.1.3 反幂法反幂法16程序见P178，例8.3则每次迭代只需解二个三角形方程组，因此实用的公式为8.1.3 反幂法反幂法17 Jacobi方法是求实对称矩阵全部全部特征值和特征向量的一种矩阵变换方法。8.2 Jacobi 方法方法 实对称矩阵A具有下列性质：(1)A的特征值均为实数;其对应的特征向量线性无关且两两正交。(2)存在正交矩阵Q，使QTAQ=diag(1,2,n),而且 Q的第i个列向量恰为i的特征向量;（3）若记A1=QTAQ,则A1仍为对称矩阵.直接找直接找Q不大可能。我们可以构造一系列特殊形式的正交阵不大可能。我们可以构造一系列特殊形式的正交阵Q1,.,Qn对对A作正交变换，作正交变换，使得使得对角元素比重对角元素比重逐次增加逐次增加，非非对角元变小。对角元变小。当非对角元已经小得无足轻重时，当非对角元已经小得无足轻重时，可以近似认为对角元就是可以近似认为对角元就是A的所有特征值。的所有特征值。Jacobi方法就是这样一类方法。方法就是这样一类方法。18平面解析几何中的平面坐标旋转变换表示平面上坐标轴旋转角的变换.8.2.1 平面旋转矩阵平面旋转矩阵(旋转正交相似变换）旋转正交相似变换）在三维空间直角坐标系中,ox1y1平面绕着oz1轴旋转角的坐标变换为 一般地,在n维向量空间Rn中,沿着xi yj平面旋转角的变换矩阵为19称Rij()为平面旋转矩阵或平面旋转矩阵或Givens变换矩阵变换矩阵.Rij()具有下列性质:i j 8.2.1 平面旋转矩阵平面旋转矩阵(旋转正交相似变换）旋转正交相似变换）208.2.1 平面旋转矩阵平面旋转矩阵(旋转正交相似变换）旋转正交相似变换）21 (2)Rij()为正交矩阵,即Rij-1()=RijT();(3)如果A为对称矩阵,则RijT()ARij()也为对称矩阵,且与A有相同的特征值.(4)RijT()A仅改变A的第i行与第j行元素,ARij()仅改变A的第i列与第j列元素.8.2.1 平面旋转矩阵平面旋转矩阵(旋转正交相似变换）旋转正交相似变换）22 设实对称矩阵A=(apq)nn,记B=RijT()ARij()=(bpq)nn则它们元素之间有如下关系:所以有8.2.1 平面旋转矩阵平面旋转矩阵(旋转正交相似变换）旋转正交相似变换）（8.14）23从而由上面两式可得 如果aij0,适当选取角,使只需角满足8.2.1 平面旋转矩阵平面旋转矩阵(旋转正交相似变换）旋转正交相似变换）24 由式(8.15),令t=tan,则t满足方程t2+2dt-1=0为保证|/4,取绝对值较小的根,有于是且且8.2.1 平面旋转矩阵平面旋转矩阵(旋转正交相似变换）旋转正交相似变换）（8.20）25非对角线元素的平方和，非对角线元素的平方和，由8.2.1 平面旋转矩阵平面旋转矩阵(旋转正交相似变换）旋转正交相似变换）得得26 我们有以下的收敛性定理保证上述计算过程。我们有以下的收敛性定理保证上述计算过程。8.2.2 Jacobi方法方法 27 Jacobi收敛性定理收敛性定理证证则有反复利用上式，即得28 设设k充分大时，有充分大时，有的近似特征值.的对角线元素就是为对角阵,则kAAD 这里需要说明一点这里需要说明一点:并不是对矩阵并不是对矩阵A的每一对非的每一对非对角线非零元素进行一次这样的变换就能得到对对角线非零元素进行一次这样的变换就能得到对角阵角阵.因为在用变换消去因为在用变换消去 的时候的时候,只有第只有第 i 行行,第第 j 行行,第第 i 列列,第第 j 列元素在变化列元素在变化,如果如果 或或 为零为零,经变换后又往往不是零了经变换后又往往不是零了.因此,Qk=RT1RT2RTk 的列向量xj(j=1,2,n)为A的近似特征向量.8.2.2 Jacobi方法方法 29的全部特征值.解解 记 A0=A,取i=1,j=2,aij(0)=a12(0)=2,于是有例例 用Jacobi 方法计算对称矩阵从而有30所以 再取i=2,j=3,aij(1)=a23(1)=2.020190,类似地可得以下依次有例题例题31从而A的特征值可取为 12.125825,28.388761,34.485401 特征向量为R1TR2TRkT例题例题32 为了减少搜索非对角线绝对值最大元素时间,对经典的Jacobi方法可作进一步改进.1.循环循环Jacobi方法方法:按(1,2),(1,3),(1,n),(2,3),(2,4),(2,n),(n-1,n)的顺序,对每个(i,j)的非零元素aij作Jacobi变换,使其零化,逐次重复扫描下去,直至S(A)为止.2.过关过关Jacobi方法方法:取单调下降收敛于零的正数序列k,先以1为关卡值,依照1中顺序,将绝对值超过1的非对角元素零化,待所有非对角元素绝对值均不超过1时,再换下一个关卡值2,直到关卡值小于给定的精度.8.2.2 Jacobi方法方法 具体算法和程序见p183，p184。33 用用Jacobi方法求得的结果精度一般都比较高，特别是求方法求得的结果精度一般都比较高，特别是求得的特征向量正交性很好。所以得的特征向量正交性很好。所以Jacobi方法是求实对称矩阵方法是求实对称矩阵全部特征值和特征向量的一个较好的方法。全部特征值和特征向量的一个较好的方法。它的弱点是计算量大，对原矩阵是稀疏矩阵，旋转变换后它的弱点是计算量大，对原矩阵是稀疏矩阵，旋转变换后不能保持其稀疏的性质。不能保持其稀疏的性质。一般适用于阶数不高的矩阵.8.2.2 Jacobi方法方法 348.3 QR方法方法60年代出现的QR算法是目前计算中小型矩阵的全部特征值与特征向量的最有效方法最有效方法。（实矩阵、非奇异。）理论依据理论依据：任一非奇异实矩阵都可分解成一个正交矩阵Q和一个上三角矩阵R的乘积，而且当R的对角元符号取定时，分解是唯一的。同理可得：Ak相似于A(k=2,3,),故他们有相同特征根。35QR方法方法收敛性收敛性 36QR方法收敛性方法收敛性37 QR方法运算量很大，为了减少运算量，常在使用QR方法之前把矩阵A简化为拟上三角矩阵。或称之为海森伯格矩阵(次对角元以下的元素全为零)。8.3.2 化一般矩阵为拟上三角矩阵化一般矩阵为拟上三角矩阵 形状为形状为可以用镜面反射矩阵将可以用镜面反射矩阵将A化为化为Hessenberg形形,下面介绍。下面介绍。38定义定义8.2为为镜面反射矩阵镜面反射矩阵，或，或Householder变换矩阵变换矩阵。Houholder矩阵矩阵H=H(v)有如下性质：有如下性质：(2)(3)记记S为以为以v为法向量的平面为法向量的平面,则几何上则几何上x与与y=Hx关于平面关于平面S对称对称。因为因为上式表明向量上式表明向量x-y与与v平行，注意到平行，注意到y与与x的长度相等，于是的长度相等，于是x经变换后的象经变换后的象y=Hx是是x关于关于s对称的向量，如下图所示。对称的向量，如下图所示。镜面反射变换镜面反射变换39xvyx-y据前面定义和性质，据前面定义和性质，有下面的定理。有下面的定理。定理定理8.4得得Hx=y。证证镜面反射变换镜面反射变换由此可得由此可得40镜面反射变换镜面反射变换程序见P187定理得证。定理得证。41从而有镜面反射变换镜面反射变换再利用c的值反过来计算42与平面旋转变换不同的是，镜面反变换可成批的消去向量的非零元与平面旋转变换不同的是，镜面反变换可成批的消去向量的非零元.将任意矩阵将任意矩阵A简化为海森伯格矩阵的步骤如下：简化为海森伯格矩阵的步骤如下：镜面反射变换镜面反射变换43镜面反射变换镜面反射变换44镜面反射变换镜面反射变换45镜面反射变换镜面反射变换例例：用Householder变换将矩阵A化为上Hessenberg阵解解：求Housholder矩阵其中46镜面反射变换镜面反射变换47镜面反射变换镜面反射变换 用用Household方法对矩阵方法对矩阵A作正交相似变换作正交相似变换,使使A相似与上相似与上Hessenberg阵，程序见阵，程序见P19048、用、用 GivensGivens变换对上变换对上HessenbergHessenberg阵作阵作QRQR分解分解具体步骤：具体步骤：假设假设（否则进行下一步）（否则进行下一步）取旋转矩阵R(1,2)49、用、用 GivensGivens变换对上变换对上HessenbergHessenberg阵作阵作QRQR分解分解50、用、用 GivensGivens变换对上变换对上HessenbergHessenberg阵作阵作QRQR分解分解51、用、用 GivensGivens变换对上变换对上HessenbergHessenberg阵作阵作QRQR分解分解52、用、用 GivensGivens变换对上变换对上HessenbergHessenberg阵作阵作QRQR分解分解53、用、用 GivensGivens变换对上变换对上HessenbergHessenberg阵作阵作QRQR分解分解因此，最多做n-1旋转变换，即得54、用、用 GivensGivens变换对上变换对上HessenbergHessenberg阵作阵作QRQR分解分解因为均为正交矩阵，故其中Q=R21TR32TRnn-1T仍为正交矩阵。可算出完成这一过程的计算量O(4n2),比一般矩阵的QR分解的计算量O(n3)少一个数量级。需要说明的是：通常用QR方法计算全部特征值，然后用反幂法求其相应的特征向量。55解：首先将A化为上Hessenberg矩阵，取56H即与A相似的上Hessenberg阵。将H进行QR分解。记B1=H5758重复上述过程，迭代11次，得算法和程序见P192，P19359练习题练习题第第194页页 习题习题88.2,8.3-8.8 60