Power method in Matlab giving inconsistent results

Question

I'm trying to make a simple implementation of the power method for determining dominant eigenvectors of a matrix in Matlab. If you repeatedly multiply any vector by a matrix, normalize that, multiply again and repeat, it should converge to an eigenvector of the matrix corresponding to eigenvalue of greatest absolute value.

This is what I tried:

A1=[13 6;-12 -5]
x=[1;1]
y=(1/norm(x))*x

powerMethod(A1,y,4)

function yn=powerMethod(A,y0,n)
    if n>0
        x1=A*y0
        y1=(1/norm(x1))*x1
        powerMethod(A,y1,n-1)
    else
        yn=y0
    end
end

I thought it was working fine at first. When x=[-1;1] it converges to [-0.7071;0.7071] which is correct. The actual eigenvectors of the matrix are (-1,1) and (-1,2). But when I changed x to [1;1], now it converges to [0.7072;-0.7070] which isn't an eigenvector of the matrix at all. x=[1;1], [1;10], and [1;1000] all converge to basically the same vector. What am I doing wrong to make different values of x converge to a vector that's not even an eigenvector? I don't have a lot of experience with Matlab so I don't know if it's because I'm misunderstanding how the power method works or how Matlab works.

Answer 1

Any non-zero multiple of an eigenvector is also an eigenvector. So your answer is only unique up to a non-zero scalar multiple. In particular, if e is an eigenvector then so is -e. Your two answers are correct, they just converge to opposites because of your starting point. You can do the multiply by hand to prove it to yourself:

>> A1 * [-0.7071;0.7071]
ans =
  -4.949700000000000
   4.949699999999999
>> ans/norm(ans)
ans =
  -0.707106781186548
   0.707106781186547
>> A1 * [0.7072;-0.7070]
ans =
   4.951600000000000
  -4.951400000000000
>> ans/norm(ans)
ans =
   0.707121061700346
  -0.707092500384339

Also, your function doesn't need to be recursive. Just put the calculations in a loop. E.g.,

A1=[13 6;-12 -5]
x=[-1;1]
y=(1/norm(x))*x

e = powerMethod(A1,y,20)

function yn = powerMethod(A,y0,n)
yn = y0 / norm(y0);
for k=1:n
    yn = A * yn;
    yn = yn / norm(yn);
end
return
end