矩阵特征值与特征向量计算.ppt

物理，力学和工程技术中很多问题在数学上归结为求矩阵特征值向量。即下面的数学问题：,第七章 矩阵特征值与特征向量计算,1 引言,1.已知： 求代数方程 的根。 称为A的特征多项式。上式展开即为: 的n个解称为A的特征值。,2.设 为A的特征值，要求相应的齐次方程组 的非零解。即求 的非零解。 此非零解x称为矩阵A的对应于 的特征向量。,回顾几个基本代数结论: 定理1 若 为A的特征值，则有: (1) (2),定理3. (Gerschgorins定理)设 则A的每一个特征值，必 属于下述某个圆盘之中: 且如果一个特征向量的第i个分量绝对值最大,则对应的特征值一定属 于第i个圆盘中.,定义1:设A为n阶实对称矩阵，对于任一非零向量x称 为对应于向量x的Rayleigh商.,关于矩阵A的特征值计算问题，当n=2,3时，可以按行列式展开的方法求 的根。但当n较大时，再按行列式展开法先求 的系数，再求 的根，工作量是非常大的，在实际应用中这种方法是不切实际的。,(1) (2) (3),定理4.设 为对称阵，特征值 ,对应的 特征向量 ,组成规范化正交组，即: 则:,在一些工程实际问题中，通常只要求出矩阵按模最大的特征值(称之为A的主要特征值)和相应的特征向量，对于这种特征值问题应用幂法是较合适的。 幂法是一种计算实矩阵 A 的主特征值的迭代法，最大优点是方法简单，对稀疏矩阵较合适，但有时收敛速度很慢。,设实矩阵 有一个完全的特征向量组,其特征值 为 ， 相应的特征向量为,已知A的主特征值是实数,且满足: 幂法的基本思想是：任取一个初始向量 , 由矩阵A构造一向量序列（称之为迭代向量序列。） : (2),2 幂法与反幂法,一 幂法,这表明： 越来越接近A的对应于 的特征向量 （数量 不计）或者k充分大时， 即迭代向量 为 的特征向量的近似向量(除一个因子外)。,而由假设, 可表示为特征向量 的线性组合：,下面考查主特征值 的计算。用 表示 的第 个分量, 则：,这表明：相邻两次迭代向量的对应分量的比值收敛到主特征值。,这种利用已知非零向量 及矩阵A的乘幂 构造向量序列 以计算A的主特征值 （由(7)式）及相应的特征向量（利用(5)式）的方法称为幂法。,由(6)及 表达式可见： 的收敛速度由比值 来确定，r 越小收敛越快。而 时收敛可能很慢。将上述讨论总结一下有结论:,定理5： 有n个线性无关特征向量。主特征值 满足 ,则对任何非零初始向量v ,有(4)、(7)成立。,若A的主特征值为实的重根，即 ，且 又设A有n个线性无关的特征向量, 对应r个线性无关特征向量为：,从而 (设 )相应的（6）（7）式成立。,为克服这个缺点，我们将迭代向量规范化。 设 ，其规范化向量为 ，其中max(v)表示向量v的绝对值最大的分量。,在定理5的条件下幂法可这样进行：任取一初始向量 构造向量序列：（ 为向量） ,（下面 为向量。） 由(3)式:,这表明：规范化向量序列收敛到主特征值对应的特征向量。,同理可证:,则 收敛速度取决于比值,(9) 有: .,解： 取初始向量 ,按(9)迭代5次得到数据如下 表: k (规范化向量) 0 1 1 1 1 1 1 1 0.2143 0.4821 1 12.00 27.00 56.00 2 0.1875 0.4483 1 8.357 19.98 44.57 3 0.1860 0.4463 1 8.168 19.60 43.92 4 0.1895 0.4460 1 8.157 19.57 43.88 5 0.1859 0.4460 1 8.156 19.57 43.88,例1：用幂法计算下面矩阵的主特征值及对应的特征向量。,二 加速方法,1 原点平移法.,应用幂法计算A的主特征值的收敛速度主要由比值 来决定。,2 Rayleigh商加速,三 反幂法,计算结果如下表：,于是求实对称阵A的特征向量问题就等于寻找正交矩阵P使 为对角阵，其主要困难在于寻找正交阵。,3 Jacobi方法,一 引言,Jacobi方法是用来求实对称矩阵的全部特征值及对应特征向量的一种变换方法。其基本思想是：通过一组平面旋转变换（正交相似变换）将对称阵A化为对角阵，得其全部特征值。,先考察平面上正交变换： ，其中,为平面旋转矩阵。则 ，其中：,先看2阶对称阵 能否找到对称阵P使A经过正交相似 变换化为对角阵？,即 （当 时，可选取 ）从而可使,将上面方法推广，首先在 Rn 中引进平面旋转变换：,二 Jacobi方法,7.4 QR算法,7.4.1 化矩阵为Hessenberg形,定理7.9 （实Schur分解定理）,为（初等）镜面反射矩阵，或Householder变换矩阵。,Houholder矩阵H=H（w）有如下性质：,(1),(2),(3) 记S为与w垂直的平面,则几何上x与y=Hx关于平面S对称。事实上,由,对应于性质（2），有下面的定理。,证,由此可得,定理得证。,（7.4.2）,（7.4.3）,例7.4,解,证,变换，有,如此类推，经n-2步对称正交相似变换，得到Hessenberg形矩阵。,7.4.2 QR算法及其收敛性,上式左边为正交阵，即,这个式子左边是下三角阵，则右边是上三角阵，所以只能是对角阵。设,定理得证。,一般按平面旋转变换或镜面反射变换作出的分解A=QR，R的对角元不,（7.4.5）,证 容易证(1)从它递推得,例7.5 用QR方法求下列矩阵的全部特征值。,该矩阵A非对称，从计算结果看，收敛于上三角阵。,（2）的计算结果为,从计算结果来看，迭代收敛于Schur分块上三角形，对角块分别是1阶和2阶子,一般在实际使用QR方法之前，先用镜面反射变换将A化为Hessenberg形矩阵H，然后对H作QR迭代，这样可以大大节省运算工作量。因为上 Hessenberg阵H的 次对角线以下元素均为零，所以用平面旋转变换作QR分解较为方便。,对i = 1,2,.n - 1,依次用平面旋转矩阵J（i , i+1）左乘H，使J（i , i+1）H的第i +1行第i列元素为零。左乘J(i,i+1)后，矩阵H的第i行与第i+1行零元素位置上仍为零，其他行不变。这样，共n-1次左乘正交矩阵后得到上三角阵R。即 =R， =J（n-1,n）J(n-2,n-1)J(1,2)。可以验证 是一个下Hessenberg阵，即U是一个上Hessenberg阵。这样，得到H的QR分解H=UR。在作QR迭代时，下一步计算RU，容易验证RU是一个上Hessenberg阵。以上说明了QR算法保持了H的上Hessenberg结构形式。,解,重复上面的过程，计算11次得,至此，不难看出，一个特征值是4，另一个特征值是-1，其他两个特征值是方程,上述用QR方法求得的特征值是该特征方程的准确解。,7.4.3 带原点位移的QR算法 前面我们介绍了在反幂法中应用原点位移的策略，这种思想方法也可用于QR算法。一般我们针对上Hessenberg矩阵讨论QR算法，并且假设每次QR迭代中产生的 都是不可约的，否则，可以将问题分解为较小型的问题。这样，带原点位移的QR算法可以描述为：,根据QR算法的收敛性质，位移量有下列两种取法：,例7.7 用带原点位移的QR算法求下列矩阵的特征值：,解 先用镜面反射变换把A化为上Hessenberg矩阵。按（7.4.3）式有,该问题如果不用带原点位移的QR算法，而是用基本QR算法，则收敛速度很慢，计算结果为,7.5 特征值问题的若干Matlab函数文件 7.5.1 幂法的Matlab函数文件 functionz,m=power_m(A,max_it,tol) n,nn=size(A);z=ones(n,1); it=0;error=100; disp(it m z(1) z(3) z(4) z(5)%该文件用于不超过5阶的矩阵。 while it tol w=A*z;ww=abs(w); k,kk=max(ww);%kk为ww最大元素的标号。 m=w(kk); %特征值的近似值。 z=w/w(kk); %特征向量的近似值。 out=it+1 m ;disp(out) Error=norm(A*z-m*z); it=it+1; End error,7.5.2 反幂法的Matlab函数文件 functionz,m=invpower(A,max_it,tol) n,nn=size(A);z=ones(n,1);it=0;error=100; L,U=lu_factor(A); while ittol w=lu_solve(L,U,z);ww=abs(w);k,kk=max(ww); m=( *z)/( *w);z=w/w(kk); out=it+1 m ;disp(out) error=norm(A*z-m*z); it=it+1; end,7.5.3 QR算法的Matlab函数文件 functionQ,R=qr_factor(A) % 矩阵A的QR分解。 n,nn=size(A);R=A;Q=eye(n); for k=1:n-1 x=zeros(n,1);x(k:n,1)=R(k:n,k); g=norm(x);v=x;v(k)=x(k)+g; s=norm(v);w=v/s;u=2* *w; R=R-w* ;Q=Q-2*Q*w* ; end function e=qr_eig(A,max) % 用QR算法求矩阵A的特征值。 for i=1:max Q,R=qr_factor(A); A=R*Q; end e=diag(A),评注 本章介绍了矩阵特征问题的幂法,Jacobi方法和QR算法，它们是求矩阵特征值和特征向量的常用数值方法。本章用到较多的线性代数知识和方法，其中一些是一般线性代数教科书上没有提到的。圆盘定理给出了特征值的大致估计。平面旋转变换和镜面反射变换是两种有力的正交相似变换工具，可以化简矩阵和作QR分解等，用于构造和分析数值方法。 幂法用于求矩阵的主特征值和主特征向量，特别适用于大型稀疏矩阵。它计算简单，但收敛速度往往不能令人满意。可以用反幂法结合位移技巧加速收敛，或求某一指定的特征值。幂法以及它的变形显然适用于对称矩阵，因为这种矩阵的特征值都是实数并且可对角化，结合矩阵的收缩方法可以求出它的全部特征值和特征向量。 Jacobi方法是古典的方法，用于求对称矩阵的全部特征值和特征向量。一般情况下，它和对称的QR方法相比，已经没有多大优越性。但是在求几乎接近对