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.
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