用QR算法求矩阵的特征值
- 格式:doc
- 大小:82.00 KB
- 文档页数:7
普通矩阵特征值的QR 算法摘 要求矩阵的特征值有多种不同的办法,本文主要介绍用QR 法求矩阵的特征值,QR 法是目前求中等大小矩阵全部特征值的最有效方法之一,使用于求实矩阵或复矩阵的特征值,它和雅可比法类似,也是一种变换迭代法。
关键词:QR 分解 迭代序列 特征值 Matlab一 、QR 方法的理论:对任意一个非奇异矩阵(可逆矩阵)A ,可以把它分解成一个正交阵Q 和一个上三角阵的乘积,称为对矩阵A 的QR 分解,即A=QR 。
如果规定R 的对角元取正实数,这种分解是唯一的。
若A 是奇异的,则A 有零特征值。
任取一个不等于A 的特征值的实数μ,则A-μI 是非奇异的。
只要求出A-μI 的特征值和特征向量就容易求出矩阵A 的特征值和特征向量,所以假设A 是非奇异的,不是一般性。
设A=A 1 ,对A 1 作QR 分解,得A 1 = Q 1R 1,,交换该乘积的次序,得A 2 = R 1Q 1=,由于Q 1正交矩阵,A 1到A 2的变换为正交相似变换,于是A 1和A 2就有相同的特征值。
一般的令A 1=A ,对k=1,2,3,…..⎩⎨⎧==+)()(1迭代定义分解kk k k k k Q R A QR R Q A这样,可得到一个迭代序列{A k },这就是QR 方法的基本过程。
二、QR 方法的实际计算步骤Householder A Hessenberg B -----→=用阵作正交相似变换上第阵一步............*::::***⎛⎫⎪ ⎪ ⎪* ⎪** ⎪ ⎪ ⎪ ⎪*⎝⎭Householder 变换:如果 v 给出为单位向量而 I 是单位矩阵,则描述上述线性变换的是 豪斯霍尔德矩阵 (v * 表示向量 v 的共轭转置)H=I -2VV*1k k kGiven k k k B Q R B B R Q +=⎧-------→→⎨=⎩用变换产生迭代序列第二步12***n λλλ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦Householder A B -----→=用阵作正交相似变换(对称阵)三对角阵*****⎛⎫ ⎪ ⎪⎪*⎪**⎪ ⎪* ⎪ ⎪*⎝⎭三、化一般矩阵为Hessenberg 阵称形如11121112122212323331n n n n n nn nn h h h h h h h h h h h H h h ---⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦的矩阵为上海森堡(Hessenberg )阵。
四.实验代码: function [H,B]=Hessenberg(A) n=length(A);B=eye(n);for k=1:n-2X=zeros(n-k,1);H=eye(n);for i=1:n-kX(i)=A(i+k,k);enda=max(abs(X));if a==0.0breakendX=X/a;c=X(1);b1=sqrt(sum(X.^2));if X(1)>=0b1=-b1;endX(1)=X(1)-b1;b=b1^2-b1*c;H0=eye(n-k)-X*X'/b;for i=1:n-kfor j=1:n-kH(i+k,j+k)=H0(i,j);end endA=H*A*H;B=B*H;endH=A;一.实验题目:QR方法求矩阵的特征和特征向量二.设计目的:学会利用镜面变换进行矩阵的QR分解及利用将幂法求特征值和特征向量,熟悉Matlab编程环境。
三.设计原理:利用镜像变换将A相似变换为Hessenberg B矩阵。
记录变换矩阵。
运用Householder矩阵进行QR分解,QR方法为:B1=BB1=Q1R1B2=R1Q1....Bm=QmRmBm+1=RmQmBm+1与Bm相似,从而特征值相等。
再利用原点位移的反幂法求B(或A)的特征向量。
反幂法用来计算矩阵按模最小的特征值及其特征向量,也可用来计算对应与一个给定近似特征值的特征向量。
设A∈R n×n为非奇异矩阵,A 的特征值依次记为|λ1|≥|λ2|≥|λ3|≥…≥|λn |,相应的特征向量为x1 ,x2,…,x n,则A-1的特征值为|1/λn|≥|1/λn-1|≥…≥|1/λ1 | ,相应的特征向量为x n ,x.所以计算A的按模最小的特征值λn的问题就是计算n-1,…,x1A-1的按模最大的特征值问题。
对于A-1应用幂法迭代(称为反幂法),可求得矩阵A-1的主特征值1/λn,从而求得A的按模最小的特征值λn。
特征值的快速求解技巧论文特征值问题是线性代数中一个重要的问题,涉及到矩阵的特征值和特征向量的求解。
求解特征值问题对于很多应用领域来说都是必不可少的,比如在机器学习、信号处理和图像处理中。
由于特征值问题的复杂性,研究人员一直在寻找快速求解特征值的方法。
本文将介绍一些特征值问题的快速求解技巧。
一、QR方法QR方法是一种常用的求解特征值问题的数值算法,它的基本思想是通过QR分解将一个矩阵转化为上三角矩阵,从而求解特征值。
QR方法可以分为隐式QR方法和显式QR方法。
隐式QR方法是通过不断迭代的方式求解特征值,在每次迭代中都会进行QR分解。
显式QR方法是直接计算出QR分解的结果,然后通过计算得到特征值。
二、幂方法幂方法是一种迭代法,主要用于求解矩阵的最大特征值及其对应的特征向量。
它的基本思想是选择一个初始向量,通过不断迭代矩阵的幂次,最终得到矩阵的最大特征值和对应的特征向量。
幂方法是一种简单而有效的求解特征值的方法,但它只能求解最大特征值及其对应的特征向量。
三、反迭代方法反迭代方法是对幂方法的改进,主要用于求解矩阵的特征值和特征向量。
它的基本思想是通过对幂方法的迭代结果应用逆迭代进行修正,从而得到更精确的特征值和特征向量。
反迭代方法通常比幂方法收敛更快,并且可以求解非对称矩阵的特征值和特征向量。
四、雅可比方法雅可比方法是一种常用的求解对称矩阵特征值的方法,它的基本思想是通过不断迭代地将矩阵中的非对角元素置为0,从而得到对称三对角矩阵。
然后,通过QR分解求解对称三对角矩阵的特征值,最终得到原矩阵的特征值。
雅可比方法可以保证得到的特征值是精确的,但是它的计算复杂度比较高,不适用于大规模矩阵的特征值求解。
五、奇异值分解方法奇异值分解方法是一种数值稳定的特征值求解方法,它可以求解任意矩阵的特征值和特征向量。
奇异值分解方法的基本思想是将矩阵分解为三个矩阵的乘积,其中一个矩阵是对角矩阵,其对角线上的元素即为矩阵的奇异值,另外两个矩阵分别是特征向量矩阵和其转置。
实验一:编程实现以下科学计算算法,并举一例应用之。
QR基本法和位移QR法矩阵特征值求解1.QR基本法算法说明:QR基本算法是求矩阵特征值的最有效和应用最广泛的一种方法方法,其基本依据是以下两个定理:1)设A是n阶矩阵,其n个特征值为λ1、λ2、…λm,那么存在一个酉矩阵U使得U T AU是以λ1、λ2、…λm为对角元的上三角矩阵。
2)设A是n阶实矩阵,那么,存在一个正交矩阵Q,使得Q T AQ为一个准上三角矩阵,它的每一个对角元是A的一个特征值,对角元上的二阶块矩阵的两个特征值是A的一对共轭复特征值。
QR基本算法的过程如下:给定循环步数M,A T=A,k=1,2,…M,计算:A k=Q k R kA k+1=R k Q kQR基本算法有如下的收敛性质:如果A的特征值满足|λ1|>|λ2|>|λ2|≥…≥|λm|,则QR基本算法产生的矩阵序列{A k}基本收敛到上三角矩阵(特别,当A为对称阵时,收敛到对角阵),对角元素收敛到A的特征值。
在MATLAB中变成实现的QR基本算法的函数为:qrtz功能:QR基本算法求矩阵全部特征值。
调用格式:l=qrtz(A,M).其中,A为已知矩阵;M为迭代步数;L为矩阵A的全部特征值。
QR基本算法的流程图:l=qrtz(A,M)QR基本算法的MATLAB程序代码如下:function l=qrtz(A,M)for i=1:M[q,r]=qr(A);A=r*q;l=diag(A);endtask11.mformat longA=[1,5,6;4,7,0;8,11,3]l=qrtz(A,20)disp('¾«È·½â')l=eig(A)运行过程和结果:2.位移QR算法位移QR法是为了加快QR算法的收敛速度,其算法的迭代过程如下:给定循环步数M,A1=Hessenberg(A),k=1,2,…M,选择μk,然后计算:A k-μk I=Q k R kA k+1= Q k Q k+μk I一般μk的选择有以下两种考虑方法:(1) 选μk =μk ,即瑞利商位移;(2) 迭代过程中,如果子矩阵的两个特征值为实数时,选最接近a k n,n 的那个作为μk ,即威尔森位移瑞利商位移QR 法流程图如下:在MATLAB 中编程实现的瑞利商位移的QR 算法的函数为:rqrtz 。
QR分解求矩阵特征值C语言程序QR分解是一种常用的矩阵分解方法,其中Q是一个正交矩阵,R是一个上三角矩阵。
利用QR分解可以求解矩阵的特征值。
以下是一个使用C语言实现QR分解求矩阵特征值的程序示例:#include <stdio.h>#include <math.h>#define N 3 // 定义矩阵的维度void eigenvalues(double A[N][N], double lambda[N]);int maindouble A[N][N] = {{2, -1, 0}, {-1, 2, -1}, {0, -1, 2}}; // 待求解的矩阵double Q[N][N], R[N][N]; // 存储QR分解结果double lambda[N]; // 存储特征值eigenvalues(A, lambda); // 求解特征值printf("矩阵的特征值为:\n");for (int i = 0; i < N; i++)printf("%.3f ", lambda[i]);}printf("\n");return 0;//初始化Q矩阵为单位阵for (int i = 0; i < N; i++)for (int j = 0; j < N; j++)if (i == j)Q[i][j]=1;elseQ[i][j]=0;}}//进行QR迭代for (int k = 0; k < N - 1; k++) //计算Ak矩阵的第k列的范数double norm = 0;for (int i = k; i < N; i++) norm += A[i][k] * A[i][k];}norm = sqrt(norm);//计算v向量double v[N] = {0};v[k] = A[k][k] + norm;for (int i = k + 1; i < N; i++)v[i]=A[i][k];}//计算P矩阵double P[N][N];for (int i = 0; i < N; i++)for (int j = 0; j < N; j++)if (i == j)P[i][j] = 1 - 2 * v[i] * v[j] / (norm * norm); elseP[i][j] = -2 * v[i] * v[j] / (norm * norm);}}//更新A和Q矩阵double temp[N][N] = {0};for (int i = 0; i < N; i++)for (int j = 0; j < N; j++)for (int m = 0; m < N; m++)temp[i][j] += P[i][m] * A[m][j];}}}for (int i = 0; i < N; i++)for (int j = 0; j < N; j++)A[i][j] = temp[i][j];Q[i][j]=Q[i][j]*P[i][j];}}}//计算矩阵Rfor (int i = 0; i < N; i++)for (int j = i; j < N; j++)R[i][j]=A[i][j];}}void eigenvalues(double A[N][N], double lambda[N])//通过矩阵R的对角线元素即可得到矩阵A的特征值for (int i = 0; i < N; i++)lambda[i] = A[i][i];}请注意,上述程序是一个简化版本的QR分解方法,并不适用于所有的矩阵。
用qr方法求矩阵的全部特征值例题矩阵的特征值问题是矩阵理论中的重要问题之一,QR方法是一种常用的求解矩阵特征值的方法。
本文将通过一个具体的例题,介绍如何使用QR方法求矩阵的全部特征值。
一、问题描述给定一个$n\timesn$矩阵$A$,我们需要求出其全部特征值。
矩阵的特征值通常可以通过求解矩阵的特征多项式来得到。
对于实对称矩阵,我们可以通过对角化矩阵的方法来求解特征值。
但对于一般矩阵,我们需要使用其他方法,如QR方法。
二、QR方法原理QR方法是基于矩阵的QR分解原理,将原矩阵分解为一个正交矩阵和一个上三角矩阵的乘积。
通过这个分解,我们可以将原矩阵的特征多项式转化为一个简单的多项式,从而方便地求解特征值。
三、例题及解答【例题】给定一个$3\times 3$矩阵:$A=\begin{bmatrix}1&2&3\\0&-2&4\\0&-1&2\end{bmatrix}$要求求出该矩阵的全部特征值。
【解法】1. 将矩阵A进行QR分解,得到正交矩阵$Q$和上三角矩阵$R$:$A=QR$2. 计算$Q^TAQ$的特征多项式,并求出全部特征值。
3. 将上三角矩阵$R$代入特征多项式中,得到原矩阵A的特征值。
【代码实现】(使用MATLAB)```matlab% 定义矩阵AA = [1 2 3; 0 -2 4; 0 -1 2];% 进行QR分解[Q, R] = qr(A);% 计算Q^TAQ的特征多项式,并求出全部特征值[eigvals,~,~] = eig(Q^TAQ);eigvals = real(eigvals); % 取实部作为特征值% 将上三角矩阵R代入特征多项式中,得到原矩阵A的特征值eigenvalues = diag(R) ./ (diag(R)+eigvals);```【结果】经过以上步骤,我们可以得到原矩阵A的全部特征值为:$\lambda_1=2.75+0. 866i,\lambda_2=2.75-0.866i,\lambda_3=4$。
Householder QR法求实对称矩阵特征值和特Householder+QR法求实对称矩阵特征值和特征向量(C程序代码)2010-12-01 10:29嗯,上一个雅可比法求n xn实对称矩阵的特征值和特征向量,是很有限制的。
虽然这个方法可靠,精度好,但是对于高于十几阶的矩阵,它就显得力不从心了,因为收敛速度比较慢。
所以对于高阶的实对称阵的特征值和特征向量的求解,目前用得比较多的还是先Householder法把对称阵变成"三对角阵",再用QR法求取特征值。
因为赶工的原因,这两个程序我没有亲自动手写(众:怎么每次都这么多理由捏)。
这里就顺便推荐一本书:《常用算法程序集(C语言描述)》,清华大学出版社出版。
这本书的内容比较广,虽然原理上讲得不详细,但是对于赶工的程序员来说那是一条不错的捷径。
有点可惜的是此书的代码写得比较凌乱,可读性极其糟糕,有空的话(又要等"有空")某U很想把里面的代码再改写封装一下…因为实在是太实用了…为了表示对作者的尊敬,这里贴出来的代码除了函数名之外,对内容就不作修改了。
#include stdlib.h#include math.h typedef double mydouble;int qr(int n,mydouble*b,mydouble*c,mydouble*q,mydouble eps,int l){int i,j,k,m,it,u,v;mydouble d,f,h,g,p,r,e,s;c[n-1]=0.0;d=0.0;f=0.0;for(j=0;j=n-1;j++){it=0;h=eps*(fabs(b[j])+fabs(c[j]));if(h d)d=h;m=j;while((m=n-1)&&(fabs(c[m])d))m++;if(m!=j){do{if(it==l)return 1;it=it+1;g=b[j];p=(b[j+1]-g)/(2.0*c[j]);r=sqrt(p*p+1.0);if(p=0.0)b[j]=c[j]/(p+r);else b[j]=c[j]/(p-r);h=g-b[j];for(i=j+1;i=n-1;i++)b[i]=b[i]-h;f=f+h;p=b[m];e=1.0;s=0.0;for(i=m-1;i=j;i--){g=e*c[i];h=e*p;if(fabs(p)=fabs(c[i])){e=c[i]/p;r=sqrt(e*e+1.0);c[i+1]=s*p*r;s=e/r;e=1.0/r;}else{e=p/c[i];r=sqrt(e*e+1.0);c[i+1]=s*c[i]*r;s=1.0/r;e=e/r;}p=e*b[i]-s*g;b[i+1]=h+s*(e*g+s*b[i]);for(k=0;k=n-1;k++){u=k*n+i+1;v=u-1;h=q[u];q[u]=s*q[v]+e*h;q[v]=e*q[v]-s*h;}}c[j]=s*p;b[j]=e*p;}while(fabs(c[j])d);}b[j]=b[j]+f;}for(i=0;i=n-1;i++){k=i;p=b[i];if(i+1=n-1){j=i+1;while((j=n-1)&&(b[j]=p)){ k=j;p=b[j];j=j+1;}}if(k!=i){b[k]=b[i];b[i]=p;for(j=0;j=n-1;j++){u=j*n+i;v=j*n+k;p=q[u];q[u]=q[v];q[v]=p;}}}return 0;}void householder(mydouble*a,int n,mydouble*q,mydouble*b,mydouble*c){int i,j,k,u;mydouble h,f,g,h2;for(i=0;i=n-1;i++){for(j=0;j=n-1;j++){u=i*n+j;q[u]=a[u];}}for(i=n-1;i=1;i--){h=0.0;if(i 1){for(k=0;k=i-1;k++){u=i*n+k;h=h+q[u]*q[u];}}if(h+1.0==1.0){c[i]=0.0;if(i==1)c[i]=q[i*n+i-1];b[i]=0.0;}else{c[i]=sqrt(h);u=i*n+i-1;if(q[u]0.0)c[i]=-c[i];h=h-q[u]*c[i];q[u]=q[u]-c[i];f=0.0;for(j=0;j=i-1;j++){q[j*n+i]=q[i*n+j]/h;g=0.0;for(k=0;k=j;k++)g=g+q[j*n+k]*q[i*n+k];if(j+1=i-1){for(k=j+1;k=i-1;k++)g=g+q[k*n+j]*q[i*n+k];}c[j]=g/h;f=f+g*q[j*n+i];}h2=f/(h+h);for(j=0;j=i-1;j++){f=q[i*n+j];g=c[j]-h2*f;c[j]=g;for(k=0;k=j;k++){u=j*n+k;q[u]=q[u]-f*c[k]-g*q[i*n+k];}}b[i]=h;}}for(i=0;i=n-2;i++)c[i]=c[i+1];c[n-1]=0.0;b[0]=0.0;for(i=0;i=n-1;i++){if((b[i]!=0.0)&&(i-1=0)){for(j=0;j=i-1;j++){g=0.0;for(k=0;k=i-1;k++)g=g+q[i*n+k]*q[k*n+j];for(k=0;k=i-1;k++){u=k*n+j;q[u]=q[u]-g*q[k*n+i];}}}u=i*n+i;b[i]=q[u];q[u]=1.0;if(i-1=0){for(j=0;j=i-1;j++){q[i*n+j]=0.0;q[j*n+i]=0.0;}}}}参数说明:void householder(mydouble*a,intn,mydouble*q,mydouble*b,mydouble*c);a-n xn实对称矩阵,用线性的方法存储(嗯,其实静态数组就是这样的)。
求矩阵特征值方法
求矩阵特征值的常见方法有以下几种:
1. 特征值分解:对于方阵A,特征值分解就是将其表示为特征向量的线性组合的形式,即A = PDP^-1,其中D是一个对角矩阵,对角线上的元素就是A的特征值,P是一个由A的特征向量组成的矩阵。
特征向量分解可通过求解矩阵的特征方程来实现。
2. 幂迭代法:幂迭代法是一种迭代计算特征值的方法,通过不断迭代向量与矩阵的乘积来逼近特征向量,从而得到特征值。
3. Jacobi方法:Jacobi方法通过不断迭代将矩阵A化为对角矩阵,每次迭代都通过旋转操作使得矩阵A中非对角线元素逐渐趋近于0,最终得到对角矩阵。
4. QR算法:QR算法通过迭代将矩阵A转化为上三角矩阵,从而得到对角线上的特征值。
5. 奇异值分解:奇异值分解将矩阵A分解为A=UDV^T的形式,其中U和V 是正交矩阵,D是一个对角矩阵,对角线上的元素即为A的奇异值。
这些方法具有不同的适用范围和计算复杂度,选择合适的方法取决于矩阵的特点和计算需求。
一种排序qr分解及其算法一、引言排序qr分解是一种常用的数值计算方法,用于将一个矩阵分解为正交矩阵和上三角矩阵的乘积。
它广泛应用于线性代数、数值分析和计算科学等领域。
本文将介绍排序qr分解的算法及其应用。
二、排序qr分解算法1. QR分解的定义QR分解是将一个$n×n$的矩阵$A$分解为两个矩阵的乘积:$A=QR$,其中$Q$是一个正交矩阵,$R$是一个上三角矩阵。
2. Gram-Schmidt方法Gram-Schmidt方法是QR分解的一种常用算法。
它通过对矩阵的列向量进行正交化处理,得到正交矩阵$Q$,然后通过矩阵乘法得到上三角矩阵$R$。
具体步骤如下:(1)对于矩阵$A$的第一列向量$a_1$,将其单位化,得到$Q$的第一列向量$q_1$。
(2)对于矩阵$A$的第二列向量$a_2$,将其在$q_1$上的投影减去,得到一个新的向量$b_2$,然后将$b_2$单位化,得到$Q$的第二列向量$q_2$。
(3)对于矩阵$A$的第三列向量$a_3$,将其在$q_1$和$q_2$上的投影减去,得到一个新的向量$b_3$,然后将$b_3$单位化,得到$Q$的第三列向量$q_3$。
以此类推,直到得到$Q$的所有列向量。
(4)对于上三角矩阵$R$,根据正交矩阵$Q$和矩阵$A$的关系,可以通过矩阵乘法得到。
3. Householder变换Householder变换是另一种常用的QR分解算法。
它通过对矩阵的行向量进行正交化处理,得到正交矩阵$Q$,然后通过矩阵乘法得到上三角矩阵$R$。
具体步骤如下:(1)选择一个向量$v$,使得$v$与矩阵$A$的第一列向量$a_1$平行,且长度为$||a_1||$。
(2)定义一个Householder矩阵$H=I-2vv^T/(v^Tv)$,其中$I$是单位矩阵。
(3)将矩阵$A$乘以$H$,得到一个新的矩阵$A_1$,使得$A_1$的第一列向量的其他元素都变为0。
QR分解迭代求特征值——原⽣python实现(不使⽤numpy)QR分解:有很多⽅法可以进⾏QR迭代,本⽂使⽤的是Schmidt正交化⽅法具体证明请参考链接迭代格式实际在进⾏QR分解之前⼀般将矩阵化为上hessnberg矩阵(奈何这个过程⽐较难以理解,本⼈智商不够,就不做这⼀步了哈哈哈)迭代终⽌条件看了很多⽂章都是设置⼀个迭代次数,感觉有些不是很合理,本来想采⽤A(k+1)-A(k)的对⾓线元素的⼆范数来作为误差的,但是我有没有⼀些严格的证明,所以本⽂也采⽤⽐较⼤众化的思路,设置迭代次数。
Python实现1 M = [[2, 4, 2], [-1, 0, -4], [2, 2, 1]]23import copy4import math567class QR(object):89def__init__(self, data):10 self.M = data11 self.degree = len(data)1213def get_row(self, index):14 res = []15for i in range(self.degree):16 res.append(self.M[i][index])17return res1819def get_col(self, index):20 res = []21for i in range(self.degree):22 res.append(self.M[i][index])23return res2425 @staticmethod26def dot(m1, m2):27 res = 028for i in range(len(m1)):29 res += m1[i] * m2[i]30return res3132 @staticmethod33def list_multi(k, lt):34 res = []35for i in range(len(lt)):36 res.append(k * lt[i])37return res3839 @staticmethod40def one_item(x, yArr):41 res = [0 for i in range(len(x))]42 temp_y_arr = []4344 n = len(yArr)45if n == 0:46 res = x47else:48for item in yArr:49 k = QR.dot(x, item) / QR.dot(item, item)50 temp_y_arr.append(QR.list_multi(-k, item))51 temp_y_arr.append(x)5253for item in temp_y_arr:54for i in range(len(item)):55 res[i] += item[i]56return res5758 @staticmethod59def normal(matrix):60 yArr = []61 yArr.append(matrix[0])6263for i in range(1, len(matrix)):64 yArr.append(QR.one_item(matrix[i], yArr))65return yArr6667 @staticmethod68def normalized(lt):69 res = []70 sm = 071for item in lt:72 sm += math.pow(item, 2)73 sm = math.sqrt(sm)74for item in lt:75 res.append(item / sm)76return res7778 @staticmethod79def matrix_T(matrix):80 mat = copy.deepcopy(matrix)81 m = len(mat[0])82 n = len(mat)83for i in range(m):84for j in range(n):85if i < j:86 temp = mat[i][j]87 mat[i][j] = mat[j][i]88 mat[j][i] = temp89return mat9091 @staticmethod92def matrix_multi(mat1, mat2):93 res = []94 rows = len(mat1[0])95 cols = len(mat1)96for i in range(rows):97 temp = [0 for i in range(cols)]98 res.append(temp)99100for i in range(rows):101for j in range(cols):102 sm = 0103for k in range(cols):104 sm += (mat1[k][i] * mat2[j][k])105 res[j][i] = sm106return res107108def execute(self):109 xArr = []110for i in range(self.degree):111 xArr.append(self.get_col(i))112 yArr = QR.normal(xArr)113 self.Q = []114for item in yArr:115 self.Q.append(QR.normalized(item))116117 self.R = QR.matrix_multi(QR.matrix_T(self.Q), xArr) 118return (self.Q, self.R)119120121# A = [122# [1, 0, -1, 2, 1],123# [3, 2, -3, 5, -3],124# [2, 2, 1, 4, -2],125# [0, 4, 3, 3, 1],126# [1, 0, 8, -11, 4]127# ]128# A = [129# [1, 2, 2],130# [2, 1, 2],131# [2, 2, 1]132# ]133 A = [134 [3, 2, 4],135 [2, 0, 2],136 [4, 2, 3]137 ]138139 temp = copy.deepcopy(A)140 val = [] # 特征值141 times = 20 # 迭代次数142for i in range(times):143 qr = QR(temp)144 (q, r) = qr.execute()145 temp = QR.matrix_multi(r, q)146 temp = QR.matrix_T(temp)147148for i in range(len(temp)):149for j in range(len(temp[0])):150if i == j:151 val.append(temp[i][j])152# 特征值153print(val)结果展⽰总结使⽤QR分解迭代求特征值,收敛的⽐较快,也可以求出所有的特征值,但是如果要求特征向量的话,还是需要求解线性⽅程组(感觉很⿇烦)。
以下为4个matlab子程序和一个主程序:实现了用隐式QR方法求方阵的特征值及特征向量实质上是将一个矩阵Hessenberg化,其中要用到Householder变换,然后使用QR迭代得到特征值最后用反幂法求解特征向量。
子程序1:function [v,b]=house(x)%householder变换H:满足H=I-bvv'且Hx只有第一个分量非零n=length(x);e=norm(x,inf);x=x/e;s=x(2:n)'*x(2:n);v=zeros(n,1);v(2:n)=x(2:n);if s==0b=0;elsea=sqrt(x(1)^2+s);if x(1)<=0v(1)=x(1)-a;elsev(1)=-s/(x(1)+a);endb=2*v(1)^2/(s+v(1)^2);v=v/v(1);end子程序2:function [A,d]=Hessen(A)%将A上Hessenberg化,次对角线以下用来存储Householder变换%次对角线及以上用来存A的上Hessenberg化n=length(A);d=zeros(n-1,1);for k=1:n-2[v,beta]=house(A(k+1:n,k));A(k+1:n,k:n)=A(k+1:n,k:n)-(beta*v)*(v'*A(k+1:n,k:n));A(1:n,k+1:n)=A(1:n,k+1:n)-(A(1:n,k+1:n)*v)*(beta*v');A(k+2:n,k)=v(2:n-k);d(k)=beta;end子程序3:function v=ipow(A,t)%反幂法计算对应于A的特征值为t的特征向量vn=length(A);[L,U,P]=lu(A-t*eye(n));v=ones(n,1);for i=2:nv(i)=v(i)-L(i,1:i-1)*v(1:i-1);endv(n)=v(n)/U(n,n);for i=n-1:-1:1v(i)=v(i)-U(i,i+1:n)*v(i+1:n);v(i)=v(i)/U(i,i);endv=v/norm(v);子程序4:function H=QR2shift(H)%双重步位移的QR方法,输入H并返回其一次双重步位移迭代后的矩阵%调用格式:循环H=QR2shift(H)n=length(H);m=n-1;s=H(m,m)+H(n,n);t=H(m,m)*H(n,n)-H(m,n)*H(n,m);x=H(1,1)*H(1,1)+H(1,2)*H(2,1)-s*H(1,1)+t;y=H(2,1)*(H(1,1)+H(2,2)-s);z=H(2,1)*H(3,2);for k=0:n-3[v,b]=house([x,y,z]');q=max(1,k);H(k+1:k+3,q:n)=H(k+1:k+3,q:n)-(b*v)*(v'*H(k+1:k+3,q:n));r=min(k+4,n);H(1:r,k+1:k+3)=H(1:r,k+1:k+3)-(H(1:r,k+1:k+3)*v)*(b*v');x=H(k+2,k+1);y=H(k+3,k+1);if k<n-3z=H(k+4,k+1);endend[v,b]=house([x,y]');H(n-1:n,n-2:n)=H(n-1:n,n-2:n)-(b*v)*(v'*H(n-1:n,n-2:n));H(1:n,n-1:n)=H(1:n,n-1:n)-(H(1:n,n-1:n)*v)*(b*v');主程序:function [eigval,V]=QReig(A)%计算A的特征值以及对应的特征向量%调用格式[eigval,V]=QReig(A),其中eigval为特征值,V的列为特征向量u=1e-15; %机器精度n=length(A);H=A;H=Hessen(H);%将A上Hessenberg化for i=1:n-2H(i+2:n,i)=zeros(n-i-1,1);endwhile(1)i1=0;i2=0;for i=2:nif(abs(H(i,i-1))<=(abs(H(i,i))+abs(H(i-1,i-1)))*u);H(i,i-1)=0;endend%将次对角线上足够小的元素置零for i=n:-1:2if((i~=2 && H(i,i-1)~=0) && H(i-1,i-2)~=0)i2=i;for j=i2:-1:2if(H(j,j-1)==0)i1=j;H(i1:i2,i1:i2)=QR2shift(H(i1:i2,i1:i2));break;endendif(i1==0)H(1:i2,1:i2)=QR2shift(H(1:i2,1:i2));endbreak;endendif(i==2)break;endend%以下求复特征值for i=2:nif(H(i,i-1)~=0)p=H(i,i)+H(i-1,i-1);q=H(i,i)*H(i-1,i-1)-H(i,i-1)*H(i-1,i);delta=p^2-4*q;H(i-1,i-1)=(p+sqrt(delta))/2;H(i,i)=(p-sqrt(delta))/2;endendeigval=diag(H);eigval=sort(eigval);%解得特征值后调用ipow解特征向量V=zeros(n);for i=1:nV(1:n,i)=ipow(A,eigval(i));end。
一种排序qr分解及其算法排序QR分解及其算法一、引言排序QR分解(Ordered QR Decomposition)是一种常用的矩阵分解方法,用于将一个矩阵分解为一个正交矩阵和一个上三角矩阵的乘积。
本文将介绍排序QR分解的算法原理和应用。
二、排序QR分解的算法原理1. QR分解的基本原理QR分解是将一个矩阵A分解为两个矩阵Q和R的乘积,其中Q 是正交矩阵,R是上三角矩阵。
分解的公式为A = QR。
2. 排序QR分解的特点排序QR分解是对传统的QR分解算法进行了改进,通过对矩阵A 的列向量进行排序,使得分解后的R矩阵的非零元素更加集中在主对角线上,从而提高了计算效率。
3. 排序QR分解的步骤(1)对矩阵A的列向量进行排序,得到排序后的矩阵A'。
(2)对矩阵A'进行QR分解,得到正交矩阵Q'和上三角矩阵R'。
(3)重新排列矩阵Q'的列向量和矩阵R'的行向量,得到最终的正交矩阵Q和上三角矩阵R。
三、排序QR分解的应用排序QR分解在很多领域中都有广泛的应用,下面我们将介绍几个典型的应用场景。
1. 线性方程组求解排序QR分解可以用于求解线性方程组。
假设有一个线性方程组Ax = b,其中A是已知的系数矩阵,b是已知的常数向量,x是未知的变量向量。
首先对A进行排序QR分解,得到Q和R两个矩阵,然后通过矩阵乘法和回代法求解线性方程组。
2. 特征值分解排序QR分解可以用于求解矩阵的特征值和特征向量。
首先对矩阵A进行排序QR分解,得到Q和R两个矩阵,然后通过QR算法迭代计算矩阵A的特征值和特征向量。
3. 最小二乘拟合排序QR分解可以用于最小二乘拟合问题。
假设有一组数据点(x, y),我们希望找到一个函数y = f(x)来拟合这些数据点。
首先将数据点组成矩阵A和向量b,然后对矩阵A进行排序QR分解,得到Q和R 两个矩阵,然后通过最小二乘法求解拟合函数的系数。
四、总结排序QR分解是一种常用的矩阵分解方法,可以用于线性方程组求解、特征值分解和最小二乘拟合等问题。
qr迭代法求特征值特征向量QR迭代法是一种求解矩阵特征值和特征向量的迭代算法。
它基于矩阵的QR分解,通过迭代得到一个上三角矩阵,该矩阵的对角线上的元素就是原矩阵的特征值,而对应的特征向量可以通过反向迭代求得。
QR分解是将一个矩阵分解为一个正交矩阵和一个上三角矩阵的乘积的过程。
对于一个实对称矩阵A,可以通过QR分解得到A=QR,其中Q是一个正交矩阵,R是一个上三角矩阵。
在QR迭代法中,我们通过反复进行QR分解,得到一系列上三角矩阵,最终收敛到一个上三角矩阵,其对角线上的元素就是原矩阵的特征值。
QR迭代法的具体步骤如下:1.将原矩阵A进行QR分解,得到A=QR。
2.计算RQ,得到新的矩阵A'。
3.将矩阵A'再进行QR分解,得到A'=Q'R'。
4.重复步骤2和3,直到矩阵A'收敛到一个上三角矩阵。
5.上三角矩阵的对角线上的元素就是原矩阵的特征值。
在QR迭代法中,每一次迭代的计算量较大,因为需要进行两次矩阵乘法。
为了降低计算量,可以使用对角线P的元素来近似表示特征向量。
具体方法是对矩阵A应用反向迭代,即计算A*x=b,其中x是特征向量的近似表示,b是一个长度为n的列向量,其元素全部为0,只有第i个元素为1。
通过迭代计算,可以得到逼近特征向量的值。
QR迭代法的优点是可以求解任意矩阵的特征值和特征向量,且收敛速度较快。
它不需要知道矩阵的特征值的个数和大小,适用于实对称矩阵和复矩阵。
但它对于存在特征值重复的矩阵求解时收敛速度会较慢,需要进行更多的迭代运算。
最后,我想提醒一下,QR迭代法是一种近似算法,无法保证得到的特征值和特征向量是精确解。
在实际应用中,如果需要得到精确解,可以使用其他更精确的算法,如雅可比迭代法或Hessenberg矩阵迭代法。
qr迭代法求特征值特征向量
摘要:
一、引言
二、QR 迭代法简介
1.QR 分解的概念
2.QR 迭代法的应用
三、QR 迭代法求特征值特征向量
1.方法步骤
2.举例说明
四、结论
正文:
一、引言
在线性代数中,特征值和特征向量是矩阵理论的重要概念。
对于给定的矩阵,如果存在非零向量和一个标量,使得矩阵乘以该向量等于向量乘以标量,那么该向量就称为特征向量,对应的标量称为特征值。
特征值和特征向量在很多实际问题中都有重要的应用,比如在信号处理、图像处理等领域。
求解特征值和特征向量的方法有很多,其中QR 迭代法是一种常用的方法。
二、QR 迭代法简介
1.QR 分解的概念
QR 分解是一种矩阵分解方法,可以将一个矩阵分解为一个正交矩阵和一个上三角矩阵的乘积。
对于给定的非奇异矩阵a,可以找到一个正交矩阵q
和一个上三角矩阵r,使得a = qr。
其中,正交矩阵q 的列向量是标准正交基,即满足qiqj = 0(i≠j)的向量;上三角矩阵r 的对角线元素为正,非对角线元素为0 或负数。
2.QR 迭代法的应用
QR 迭代法可以用于求解线性方程组、最小二乘问题等。
在实际应用中,QR 分解常用于求解最小二乘问题,即求解形如ax=b 的超定方程组,其中a 为mn 的矩阵,m>n,b 为m 维向量。
通过QR 分解,可以将最小二乘问题转化为求解形如rx=q 的方程组,从而降低问题的复杂度。
一、实验名称:用QR 算法求矩阵的特征值
二、实验目的:1、通过实验进一步熟悉掌握求矩阵特征值的QR 方法及原理。
2、理解QR 方法的计算流程。
3、能够编程实现QR 方法。
三、实验内容:给定矩阵 ⎪⎪⎪
⎭
⎫ ⎝⎛=111132126A , ⎪⎪
⎪⎪⎪
⎪⎭
⎫
⎝
⎛=0100098
20
087630
7654465432H ,采用QR 方法计算A 和H 矩阵的全部特征值。
四、实验要求:
(1) 根据QR 算法原理编写程序求矩阵A 及矩阵H 的全部特征值(要求误差<10
5
-)。
(2) 直接用MATLAB 的内部函数eig 求矩阵A 及矩阵H 的全部特征值,并与(1)的结果比较。
五、QR 方法计算矩阵特征值的程序: function
[namda,time,data_na]=qr_tz(A,tol) if nargin==1; tol=1e-5; end wucha=1; time=0;
while (wucha>tol)&(time<500) [q,r]=qr(A); A1=r*q; tz0=diag(A1); tz1=diag(A); wucha=norm(tz0-tz1);
A=A1; time=time+1; data_na(time,:)=tz1; end namda=tz1; disp(‘特征值为’) namda
disp(‘第一个特征在值’) time
n1=length(data_na); n2=(1:n1)’; temp1=[n2,data_na]; subplot(2,2,1:2)
plot(date_na(:,1))
title(‘迭代次数为’) grid
subplot(2,2,3)
plot(data-na(:,2))
title(‘第二个特征值’)grid
subplot(2,2,4)
plot(data-na(:,3))
title(‘第三个特征值’) grid
六、实验结果:
>> A=[6,2,1;2,3,1;1,1,1];[namda,time,data_na]=qr_tz(A,1e-5);特征值为
namda =
迭代次数为
time =
6
图 1
>> A=[6,2,1;2,3,1;1,1,1];[V,D]=eig(A,'nobalance'), V =
D =
0 0
0 0
0 0
>>
A=[2,3,4,5,6;4,4,5,6,7;0,3,6,7,8;0,0,2,8,9;0,0,0,1,0];[namda,time,data_na]=qr_t z(A,1e-5);
特征值为
namda =
迭代次数为
time =
22
图 2
>>
A=[2,3,4,5,6;4,4,5,6,7;0,3,6,7,8;0,0,2,8,9;0,0,0,1,0];[V,D]=eig(A,'nobalance'), V =
D =
0 0 0 0
0 0
0 0 0
0 0 0
0 0
0 0
表1 用两种方法求得矩阵A的全部特征值
表2 用两种方法求得矩阵H的全部特征值
七、实验结果分析:
从图1和图2中可以看出在迭代前几次可能会有一些波动,但逐渐趋于平稳,并且收敛速度快,算法稳定。
从表1和表2可以看出直接用MATLAB的内部函数eig求矩阵A及矩阵H的全部特征值与QR方法的结果进行比较时其特征值相差不大。
虽然用MATLAB内部函数eig直接求出矩阵的特征值的方法比较快,但是其计算结果和用QR方法相比存在一定的误差,即没有用QR 方法求得精确值高和特征值全。