MATLAB: Inverse power method for smallest eigenvector calculation

eigenvectoreigsinverse iteration

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.

Best Answer

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
Related Question