Jacobi算法和根据定义求特征值和特征向

Jacobi算法和根据定义求特征值和特征向 Jacobi算法和根据定义求特征值和特 征向 【转】 [Matlab]Jacobi算法和根据定义求特征值和特征向 量 2011年08月23日 matlab内建的求特征值和特征向量的函数是eig。 下面是网上搜集和改写的部分求特征值和特征向量的算法 一、Jacobi算法求特征值和特征向量 1、这个比较简单，不过没加容错处理，遇到非对称矩阵或者含有复数的矩阵会死循环。 % Jacobi方法求特征值 % [D,V] = YKB(A) % 输入:A为方阵 % 输出:D为特征值，V为特征向量 function [d,v]=ykb(a) [m,n]=size(a); if ndims(a)~=2 | m~=n, error('Matrix must be square');end a1=a; k=1; b=a1-diag(diag(a1)); [m1,i]=max(abs(b)); [m2,j]=max(m1); i=i(j); v=eye(m,n); while sum(sum(b.*b))>eps k=k+1; x=(a1(i,i)-a1(j,j))/(2.0*a1(i,j)); u=acot(x)/2; v1=eye(m,n); v1(i,i)=cos(u); v1(j,j)=cos(u); v1(i,j)=-sin(u); v1(j,i)=sin(u); v=v*v1; a1=v1'*a1*v1; b=a1-diag(diag(a1)); [m1,i]=max(abs(b)); [m2,j]=max(m1); i=i(j); end v=v'; d=diag(a1); 2、这个也是jacobi算法，容错处理做的比较好 function [d,v,history,historyend,numrot]=jacobi(a_in,iterma x) % [d,v,history,historyend,numrot]=jacobi(a_in,iterma x) % computes the eigenvalues d and % eigenvectors v of the real symmetric matrix a_in, % using rutishausers modfications of the classical % jacobi rotation method with treshold pivoting. % history(1:historyend) is a column vector of the length of % total sweeps used containing the sum of squares of % strict upper diagonal elements of a. a is not % touched but copied locally % the upper triangle is used only % itermax is the maximum number of total sweeps allowed % numrot is the number of rotations applied in total % check arguments siz=size(a_in); if siz(1) ~= siz(2) error('jacobi : matrix must be square ' ); end if norm(a_in-a_in',inf) ~= 0 error('jacobi ; matrix must be symmetric '); end if ~isreal(a_in) error(' jacobi : valid for real matrices only'); end if nargin==1 itermax=20; end n=siz(1); v=eye(n); a=a_in; history=zeros(itermax,1); d=diag(a); bw=d; zw=zeros(n,1); iter=0; numrot=0; while iter < itermax iter=iter+1; history(iter)=sqrt(sum(sum(triu(a,1).^2))); historyend=iter; tresh=history(iter)/(4*n); if tresh ==0 return; end for p=1:n for q=p+1:n gapq=10*abs(a(p,q)); termp=gapq+abs(d(p)); termq=gapq+abs(d(q)); if iter>4 && termp==abs(d(p)) && termq==abs(d(q)) % annihilate tiny elements a(p,q)=0; else if abs(a(p,q)) >= tresh %apply rotation h=d(q)-d(p); term=abs(h)+gapq; if term == abs(h) t=a(p,q)/h; else theta=0.5*h/a(p,q); t=1/(abs(theta)+sqrt(1+theta^2)); if theta < 0 t=-t; end end c=1/sqrt(1+t^2); s=t*c; tau=s/(1+c); h=t*a(p,q); zw(p)=zw(p)-h; %accumulate corrections to diagonal elements zw(q)=zw(q)+h; d(p)=d(p)-h; d(q)=d(q)+h; a(p,q)=0; %rotate, use information from the upper triangle of a only %for a pipelined cpu it may be better to work %on full rows and columns instead for j=1:p-1 g=a(j,p); h=a(j,q); a(j,p)=g-s*(h+g*tau); a(j,q)=h+s*(g-h*tau); end for j=p+1:q-1 g=a(p,j); h=a(j,q); a(p,j)=g-s*(h+g*tau); a(j,q)=h+s*(g-h*tau); end for j=q+1:n g=a(p,j); h=a(q,j); a(p,j)=g-s*(h+g*tau); a(q,j)=h+s*(g-h*tau); end % accumulate information in eigenvectormatrix for j=1:n g=v(j,p); h=v(j,q); v(j,p)=g-s*(h+g*tau); v(j,q)=h+s*(g-h*tau); end numrot=numrot+1; end end %if end % for q end % for p bw=bw+zw; d=bw; zw=zeros(n,1); end %while 二、根据定义求特征值和特征向量。 设A是一个n阶方阵，I是单位矩阵，如果存在一个数x及一个n维 非零列v，使得Av=xv,或 (A-xI)v=0，称数x为方阵A的一个值，称向量特征非零列向量v为方阵A的对应于(或属于)特征值x的特征向量。 设A是n阶方阵， I是单位矩阵， 如果存在一个数x，使得 A-xI 是奇异矩阵(即不可逆矩阵， 亦即行列式的值为零)， 那么x称为A的特征值。 第二个定义实际上是第一个定义的推论，由 (A-xI)v=0 推出行列式 |A-xI| = 0。 所以，求值，用 |A-xI| = 0 来解出值x。 特征特征 求特征向量，用 (A-xI)v = 0 来解出特征向量v。 下面这个是根据定义来的算法，效率低了点，不过可以计算非对称矩阵。 function [d,v]=defeig(a) % 根据定义求特征向量和特征值 % [d,v]=defeig(a) % a为方阵,输出d为特征值，v为特征向量 siz=size(a); if ndims(a)~=2 || siz(1)~=siz(2) error('Matrix must be square!'); end % 求特征值 syms x; b = a-x*eye(siz(1)); % 由特征值的定义，|A-xI|=0，解出的x就是 特征值。 b1= det(b); % 求行列式的值 d = subs(solve(b1)); % 用solve解出来是sym类型的，用subs转 换成数值型。 % 求特征向量 for k=1:siz(1) b = a-d(k)*eye(siz(1)); % (A-xI)v=O中的系数矩阵A-xI v(:,k) = null(b); % 解齐次线性方程组用null，v的每一列为一个 特征向量 end 三、这个不知道叫什么算法，不过非对称矩阵时结果不正确。 function [x, v] = findeigen(A) % Usage: % compute the subsequent eigenvalue and eigenvector % Input: % A orginal matrix % x0 initial eigen value % v0 initial eigen vector % Output: % x final eigen value % v final eigen vector % Author: % % Date: % % maximum iteration itermax = 100 ; % minimum error errmax = 1e-8 ; N = size(A, 1) ; xnew = 0 ; vnew = ones(N, 1) ; x = zeros(1, N) ; v = zeros(N, N) ; % calculate eigenvalue use The Deflation Method B = A ; for num1 = 1 : N if num1 > 1 B = B - xnew * vnew * vnew' ; else end [xnew, vnew] = powermethod(B, itermax, errmax) ; % call power method to obtain the eigenvalue x(num1) = xnew ; end % calculate eigenvalue use The Inverse iteration method u = 0.1 ; % shift value for num1 = 1 : N C = inv(A - (x(num1)-u)*eye(N)) ; [xnew, vnew] = powermethod(C, itermax, errmax) ; % call power method to obtain the eigenvector v(:, num1) = vnew ; end function [x, v] = powermethod(A, itermax, errmax) N = size(A, 1) ; xold = 0 ; vold = ones(N, 1) ; for num2 = 1 : itermax vnew = A * vold ; [xnew, i] = max(abs(vnew)) ; % get eigenvalue xnew = vnew(i) ; vnew = vnew/xnew ; % normlize errtemp = abs((xnew-xold)/xnew) ; % calculate the error if(errtemp < errmax) break ; end xold = xnew ; vold = vnew ; end x = xnew ; v = vnew ;