Power Method in MATLAB

Question

I would like to implement the Power Method for determining the dominant eigenvalue and eigenvector of a matrix in MATLAB.

Here's what I wrote so far:

%function to implement power method to compute dominant
%eigenvalue/eigenevctor
function [m,y_final]=power_method(A,x);
m=0;
n=length(x);
y_final=zeros(n,1);
y_final=x;
tol=1e-3;
while(1)
    mold=m;
 y_final=A*y_final;
 m=max(y_final);
 y_final=y_final/m;
 if (m-mold)<tol
     break;
 end
end
end

With the above code, here is a numerical example:

 A=[1 1 -2;-1 2 1; 0 1 -1]

A =

     1     1    -2
    -1     2     1
     0     1    -1

>> x=[1 1 1];
>> x=x';
>> [m,y_final]=power_method(A,x);
>> A*x

ans =

     0
     2
     0

When comparing with the eigenvalues and eigenvectors of the above matrix in MATLAB, I did:

[V,D]=eig(A)

V =

    0.3015   -0.8018    0.7071
    0.9045   -0.5345    0.0000
    0.3015   -0.2673    0.7071


D =

    2.0000         0         0
         0    1.0000         0
         0         0   -1.0000

The eigenvalue coincides, but the eigenvector should be approaching [1/3 1 1/3]. Here, I get:

 y_final

y_final =

    0.5000
    1.0000
    0.5000

Is this acceptable to see this inaccuracy, or am I making some mistake?

Solution

You have the correct implementation, but you're not checking both the eigenvector and eigenvalue for convergence. You're only checking the eigenvalue for convergence. The power method estimates both the prominent eigenvector and eigenvalue, so it's probably a good idea to check to see if both converged. When I did that, I managed to get [1/3 1 1/3]. Here is how I modified your code to facilitate this:

function [m,y_final]=power_method(A,x)
m=0;
n=length(x);
y_final=x;
tol=1e-10; %// Change - make tolerance more small to ensure convergence
while(1)
     mold = m;
     y_old=y_final; %// Change - Save old eigenvector
     y_final=A*y_final;
     m=max(y_final);
     y_final=y_final/m;
     if abs(m-mold) < tol && norm(y_final-y_old,2) < tol %// Change - Check for both
         break;
     end
end
end

When I run the above code with your example input, I get:

>> [m,y_final]=power_method(A,x)

m =

     2


y_final =

    0.3333
    1.0000
    0.3333

---

On a side note with regards to eig, MATLAB most likely scaled that eigenvector using another norm. Remember that eigenvectors are not unique and are accurate up to scale. If you want to be sure, simply take the first column of V, which coincides with the dominant eigenvector, and divide by the largest value so that we can get one component to be normalized with the value of 1, just like the Power Method:

>> [V,D] = eig(A);
>> V(:,1) / max(abs(V(:,1)))

ans =

    0.3333
    1.0000
    0.3333

This agrees with what you have observed.