inverse power method for smallest eigenvector calculation

52 次查看(过去 30 天)
Hi,
I need to calculate the smallest eigenvector of a matrix. I use eigs(A,1,'sm') and I would like to compare the result with inverse power method and see how many iteration it takes to calculate the same result. Here is the code,
function [x,iter] = invitr(A, ep, numitr)
[m,n] = size(A);
if m~=n
disp('matrix A is not square') ;
return;
end;
x=rand(n,1);
for k = 1 : numitr
iter = k;
xhat = A \ x;
x = xhat/norm(xhat,2);
if norm((A)* x , inf) <= ep
break;
end;
end;
end
First of all after some point the eigenvector stops converging yet the result comes with a sign change. That is ;
A =
1 3 5
2 6 8
3 8 10
x1 = (by eigs)
0.8241
-0.5356
0.1844
x2 = (by inverse power method)
-0.8241
0.5356
-0.1844
iter =
100
and
x1 =
-0.8241
0.5356
-0.1844
x2 =
0.8241
-0.5356
0.1844
iter =
1000
I might have done something wrong with my function, yet I don't understand why the sign changes with eigs. I appreciate all comments.

采纳的回答

ttopal
ttopal 2017-2-22
Here is another version of inverse iteration method, where if statement works fine.
function [x,iter] = invitr(A, ep, numitr)
%INVITR Inverse iteration
%[x,iter] = invitr(A, ep, numitr) computes an approximation x, smallest
%eigenvector using inverse iteration. initial approximation is vector of ones,
%ep is the tolerance and numitr is the maximum number of iterations.
%If the iteration converged, iter is the number of iterations
%needed to converge. If the iteration did not converge,
%iter contains numitr.
%This program implements Algorithm in
%http://www.netlib.org/utk/people/JackDongarra/etemplates/node96.html
%input : Matrix A, ep and integer numitr
%output : vector x and integer iter
[m,n] = size(A);
if m~=n
disp('matrix A is not square');
return;
end;
y=ones(n,1);
for k = 1 : numitr
iter = k;
v = y/norm(y,2);
y = A\v;
th =v'*y;
if norm(y-th.*v,2) < ep*abs(th)
break;
end;
end;
x = y/th;
end

更多回答(1 个)

David Goodmanson
David Goodmanson 2017-2-22
编辑:David Goodmanson 2017-2-22
Hello Turker, There is nothing wrong here. The eigenvalue equation is
A*v = lambda*v
and so for the eigenvector, both v and -v are good solutions. The eigenvalue can't do that but it comes out correctly, which you can verify (since all components of your eigenvector are well away from equaling zero):
>> (A*x2)./x2
ans =
0.1689
0.1689
0.1689
compared to
>> eig(A)
ans =
17.5075
-0.6764
0.1689
  1 个评论
ttopal
ttopal 2017-2-22
Oh! Thanks David. Do you think that might be the reason why my if statement doesn't help? I mean if 100 iteration were enough to calculate good eigenvector why it would continue for 1000?

请先登录,再进行评论。

类别

Help CenterFile Exchange 中查找有关 Linear Algebra 的更多信息

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by