inverse power method for smallest eigenvector calculation

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.

Answer 1

[ttopal]

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

Answer 2

[David Goodmanson]

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

[ttopal]

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?