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 ;