QR方法求矩阵全部特征值

1、值分析QR方法求矩阵全部特征值用 QR 算法求矩阵特征值：621( i ) (A)2 3 1111要求：ii2345644567H 0367800289000101）根据QR 算法原理编制求（i ）与（ ii计算结果（要求误差10 5）（ 2）直接用现有的数学软件求（i ） ， （ ii ）的全部特征值，并与（1）的结果比较。QR 方法是目前公认为计算矩阵全部特征值的最有效的方法。它适用于任一种实矩。QR方法的原理是利用矩阵的正交分解产生一个与矩阵A 相似的矩阵迭代序列， 这个序列将收敛于一个上三角阵或拟上三角阵，从而求得原矩阵A的全部特征值。QR是一个迭代算法，每一步迭代都要进行QR分解，再

2、作逆序的矩阵乘法。要是对矩阵A 直接用QR方法，计算量太大，效率不高。因此，一个实际的QR方法计算过程包括两个步骤，首先要对矩阵A 用初等相似变换约化为上Hessenberg矩阵，再用QR方法求上Hessenberg矩阵的全部特征值。示意如下：用 Householder阵作正交相似变换第一步A 不需迭代，对A的每一列计算一次上 Hessenberg阵 Bxx用 QR变换产生迭代序列，迭代计算第二步 BBkQkRkBkRkQk对 B 矩阵的约化只需将每列次对角线上的元素约化为0。因此常常用平面旋转阵( Givens 变换阵)来进行约化。一、QR方法原理及收敛性设 A Rn n已实现了QR分解，记

3、AA1Q1R1其中Q1是正交阵，R1 是上三角阵。因为Q1TQ1 1，用Q1对A1作正交相似变换有Q1T A1Q1A2可改写为A2Q1 1 A1Q1 Q1 1Q1R1Q1R1Q1显然A2 只是A1 的QR分解因子阵的逆序相乘，而且A2 与原矩阵A1 有相同的特征值。对矩阵A2再重复以上过程并继续下去，可以得到一个与原矩阵A有相同特征值的矩阵序列。其过程如下：记A1A对A1 作正交分解A1Q1R1作矩阵A2 R1Q1，A2 A1 A对 A2 作正交分解A2 Q2 R2作矩阵A3 R2 Q2 ， A3 A2 A1 ， A3 A1 A重复以上过程可得一般的形式为对 Ak 作正交分解Ak Qk Rk(

4、A1 A)构成矩阵序列Ak 1 RkQk ( k=1,2 )Ak 1 A从矩阵 A开始得到一个矩阵序列A1 ,A2, A3, Ak ,这个矩阵序列中每一个矩阵都与原矩阵A( A1) 相似，即都有与A相同的特征值。这个矩阵序列Ak 在实质上收敛于依次以1, 2, , n为对角元的上三角阵。具体可以表示为1x x2lim Akkxn其中aii (k)i (k )i j时 aij (k)0 (k )i j时 aij ( k)的极限不一定存在。二、用正交相似变换约化矩阵为上Hessenberg 阵用 Householder 变换可以将一个向量指定的某个分量以下的各分量变为0。 我们只要求消掉A的次对角

5、线以下的元素，即将A约化为上Hessenberg阵。为了使变换前后矩阵的特征值不变，需要用 Householder 矩阵对A作相似变换，即用正交阵同时左乘和右乘A 时，原来已变为0 的元素不再改变。若设H 1 是Householder 矩阵，用它对A的第一列元素的变换示意如下：xxxxxxxxH 1 AH 10 x0x依次对 A 的各列进行类似的变换，一共要进行n 2 次变换，最终可以得到一个与原矩阵A有相同特征值的上Hessenberg 阵。xxxH n 2H n 3 H 2H 1AH 1 H 2 H n 3H n 2xxHouseholder 矩阵来实现。A为上 Hessenberg 阵的

6、过程可以用一系列H kCkk e1其中，Rk计算公式为Ck(ak 1,k(k),ank(k)TRn k,CkCk/Ckn(k)(k)2 1/2k sgn(ak 1,k )( (aik / Ck ) ) ,其中，对于每一个H k 有i k 1kukCkke1 ,k k(ak 1,k(k) k),RkI1ukukT.经过 n 2步约化就可得到一个上Hessenberg阵（ A的第 n 1 列不需要约化）An 1Hn 2 H2H1AH1H2 Hn 2(记为)Ba11(1)(2)a222n2(n 1)(n 1)(n 1)(n 1) n(n 1)(n 1)nn三、 Hessenberg 阵的QR算法设矩

7、阵 A Rn n ，其特征值都是实数。若已用Hessenberg 阵xx xPAPT x x记为 BxxHouseholder 变换约化为上对已得到的上Hessenberg 阵可用QR变换，经过迭代过程约化为上三角形矩阵以求出A的特征值。只要A的特征值是实数，将B约化为上三角形矩阵总是可以做到的。由于 B 矩阵结构上的特点，对B 矩阵的约化只需将每列次对角线上的元素约化为 0. 可用平面旋转阵（Givens 变换阵）来进行约化。n 阶方阵 R（i, j）为平面旋转阵11cossin1R(i , j)1sincos(j i)Rij非奇异，显然有det(Rij ) 1。 Rij 还是一个非对称的正

8、交阵，有RijRij(T)I ， Rij(T)也是一个平面旋转阵。Rij 有以下几种作用：1 、Rij左乘向量x (x1,x2,xn)T只改变 x 的第 i 个和第 j 个分量。现构造Rij 对 x 作变换使Rij x的结果将x 的第 j 个分量约化为0。yixi cos xj sin令 y Rijx，有yjxi sin xj cosykxk(k 1,2, ,n;k i, j)调整 角可使yj 0。若记 C cos ,S sin ，按下式选取2222C xi / xi x j xi /r , S xj / xi xj xj / r于是yi Cxi Sxj rxix jy jSxi Cx j 0

9、有 Rij x(x1 , xi 1 , r, xi 1 , x j 1 ,0, xj 1 , xn ) 。2 、对非零的n 维向量 x连续左乘Ri,i 1, Ri,(i 2), Rin ，可将 x的第 i+1 到第 n 个分量都约化为零；即RinRi(n 1) Ri(i 2) Ri(i 1)x(x1,xi 1,ri,0,0)T其中ri xi2 xi 1 2xn 23 、用Rij 对矩阵A作变换得到的结论是RijT左乘A只改变A的第i,j行；RijT右乘A只改变A的第i,j列；用Rij对 A作正交相似变换Rij T ARij 改变了A的i 行和 j 行以及 i 列和 j 列。用一系列R(i, j

10、)连续左乘矩阵A，可以将矩阵A化为上三角阵。数据结构描述主要数据成员说明double Ann存放矩阵Adouble Qnn,存放QR分解式的正交矩阵Qdouble Rnn存放QR分解式的上三角阵Rdouble pnnGivens 矩阵pdouble InnN 阶单位阵double Vnn存放Q矩阵的转置double Tnn初等反射阵Tdouble eps精度double max最大值double det存放行列式的值int count存放迭代次数主要函数成员说明double Det(double Lnn)用高斯列主元方法求行列式int Non_singularMatrix(double Lnn

11、)判断是否是非奇异矩阵void Disp(double Hnn)输出矩阵int IsZero(double a,int j)判断数组是否全为0int sgn(double y)符号函数void Hessenberg(double Ann)将矩阵化为上Hessenberg 矩阵int IsHessenberg(double Enn)判断是否是上Hessenberg 矩阵void QRAlgorithm(double Ann)QR算法求特征值void SeekEigenvalue(double Ann)判断是否满足QR算法条件，满足则进行QR方法求特征值算法的描述（流程图）源程序C程序为：#inc

12、lude#include#define n 3#define eps 1E-5double Det(double Lnn)/ 用高斯列主元方法求行列式double det=1,t;for(int k=0;kn-1;k+)double max=Lkk;int Ik=k;for(int j=k+1;jn;j+)if(maxLjk)max=Ljk;Ik=j;if(max=0) return 0;if(Ik!=k)for(int j=k;jn;j+)t=LIkj;LIkj=Lkj;Lkj=t;det*=-1;for(int i=k+1;in;i+)t=Lik/Lkk;Lik=t;for(int j=k

13、+1;jn;j+)Lij=Lij-t*Lkj;det*=Lkk;if(Ln-1n-1=0) return 0;elseprintf( 矩阵 A 的行列式为:%.5fn,det*Ln-1n-1);return det*Ln-1n-1;int Non_singularMatrix(double Lnn)/判断是否是非奇异矩阵printf( 判断矩阵A 的行列式是否为0?n);if(Det(L)!=0) return 1;return 0;void Disp(double Hnn)/ 输出矩阵for(int i=0;in;i+)for(int j=0;jn;j+)printf(%.6f ,Hij);

14、printf(n);int IsZero(double a,int j)/判断数组是否全为0for(int i=0;i0) return 1;else if(y=0) return 0;else return -1;void Hessenberg(double Ann)/ 将矩阵化为上Hessenberg 矩阵printf( 原矩阵为:n);Disp(A);/ 输出原矩阵doubleTnn,Bnn,Cnn,cn-1,vn-1,un-1,Rn-1n-1,In-1n-1, t,w,s;int i,j,k,m;for(k=0;kn-2;k+)for(i=0;in-k-1;i+)for(j=0;jn-

15、k-1;j+)if(i=j) Iij=1;else Iij=0;/定义单位阵Idouble max=fabs(Ak+1k);/ 求最大值for(i=0;in-k-1;i+)if(maxfabs(Ai+k+1k)max=fabs(Ai+k+1k);for(i=0;in-k-1;i+)ci=Ai+k+1k/max;/ 标准化数组if(IsZero(c,n-k-1)/ 判断数组是否全为0continue;/ 数组为0，则这一步不需要约化for(i=0,t=0.0;in-k-1;i+)t+=ci*ci;vk=sgn(Ak+1k)*sqrt(t);/u0=c0+vk;for(j=1;jn-k-1;j+)

16、uj=cj;/ 求矩阵 uw=vk*(c0+vk);for(i=0;in-k-1;i+)for(j=0;jn-k-1;j+)Rij=Iij-ui*uj/w;for(i=0;in;i+)for(j=0;jn;j+)if(i=j) Tij=1;else Tij=0;/定义矩阵T 为单位阵for(i=0;in-k-1;i+)for(j=0;jn-k-1 ;j+)Ti+k+10+k+1=Rij;/初等反射阵 Tfor(i=0;in;i+) / 矩阵 T左乘矩阵 Afor(j=0;jn;j+)for(m=0,s=0.0;mn;m+)s+=Tim*AmD;BiU=s;)for(i=0;in;i+) / 矩

17、阵 T右乘矩阵 Afor(j=0;jn;j+)for(m=0,s=0.0;mn;m+)s+=Bim*Tmj;CiU=s；)for(i=0;in;i+)for(j=0;jn;j+)AiU=Cij;)printf(原矩阵化成上 Hessenberg 矩阵后为:n);Disp(A);)int lsHessenberg(double Enn)/判断是否是上 Hessenberg 矩阵for(int i=2;in;i+)for(int j=0;j+1i;j+)if(Eij!=0)return 0;return 1;int count=1;void QRAlgorithm(double Ann)/QR 算

18、法求特征值int i,j,k,m;double pnn,Qnn,Rnn,Fnn,Vnn,c,s,t,y,max;for(i=0;in;i+)/ 赋Q为单位阵for(int j=0;jn;j+)if(i!=j) Qij=0;else Qij=1;for(m=0;mn-1;m+)/ 对矩阵 A 进行QR分解for(i=0;in;i+)/ 赋 p 为单位阵for(j=0;jn;j+)if(i!=j) pij=0;else pij=1;c=Amm/sqrt(Amm*Amm+Am+1m*Am+1m);s=Am+1m/sqrt(Amm*Amm+Am+1m*Am+1m);为 Givenspmm=c;pmm+

19、1=s;pm+1m=-s;pm+1m+1=c;/p矩阵for(i=0;in;i+)/p 矩阵左乘W矩阵,p 矩阵左乘Q矩阵for(j=0;jn;j+)t=0;y=0;for(k=0;kn;k+)t+=pik*Akj;y+=pik*Qkj;Rij=t;Fij=y;for(i=0;in;i+)/A,R赋新值for(j=0;jn;j+)Aij=Rij;Qij=Fij;for(i=0;in;i+)/Q 转置for(j=0;jn;j+)Vij=Qji;for(i=0;in;i+)for(j=0;jn;j+)Qij=Vij;for(i=0;in;i+)/ 矩阵 A为矩阵R和矩阵Q的乘积for(j=0;jn

20、;j+)t=0;t+=Rik*Qkj;Aij=t;count+;max=fabs(A10);/ 求A矩阵下三角的绝对值最大值for(i=1;in;i+)for(j=0;jmax)max=fabs(Aij);if(maxeps)/ 满足精度条件，输出Q矩阵，R矩阵以及特征值printf( 在精度%.1e下 ,QR算法迭代次数为:%d次 n 输出矩阵A%d:n,eps,count-1,count);Disp(A);/输出A矩阵printf(输出矩阵Q:n);Disp(Q);printf(输出矩阵R:n);Disp(R);printf( 输出 A 矩阵的全部特征值:);for(i=0;in;i+)if(i%3=0) printf(n);printf(a%d=%.6ft,i+1,Aii);printf(n);return;QRAlgorithm(A);判断是否满足QR算法条件，满足则void SeekEigenvalue(double Ann)/进行QR方法求特征值double Znn;for(int i=0;in;i+)for(int j=0;jn;j+)Zij=Aij;printf( 判断矩阵A 是否是非奇异矩阵?n);if(Non_singularMatrix(Z)=0)printf( 矩阵 A 不是非奇异矩阵! 不符合分解条件!n);return;elsepr