I'm trying to understand the power iteration to calculate the eigenvalues of a matrix.
I followed the algorithm from en.wikipedia.org/wiki/Power_iteration#The_method:
from math import sqrt
def powerIteration(A):
b = [random() for i in range(len(A))]
tmp = [0] * len(A)
for iteration in range(10000):
for i in range(0, len(A)):
tmp[i] = 0
for j in range(0, len(A)):
tmp[i] += A[i][j] * b[j]
normSq = 0
for k in range(0, len(A)):
normSq += tmp[k] * tmp[k]
norm = sqrt(normSq)
for i in range(len(A)):
b[i] = tmp[i] / norm
return b
When I run powerMethod([[0.0, 1.0], [1.0, 0.0]])
it returns random pair of numbers, such as: [0.348454142915605, 0.9373258293064111]
or [0.741752215683863, 0.6706740270266026]
Question #1 - why are those numbers random? Obviously I started with random vector b
but I hoped it would converge.
Question #2 - there is this Online Matrix Calculator to which when I feed:
0 1
1 0
it returns:
Eigenvalues:
( 1.000, 0.000i)
(-1.000, 0.000i)
Eigenvectors:
( 0.707, 0.000i) (-0.707, 0.000i)
( 0.707, 0.000i) ( 0.707, 0.000i)
If I understood correctly, returning b
should get one of those eigenvectors, but it does not. Why is the output so different?
Question #3 - what should I add to the above algorithm so that it returns one eigenvalue (In this example it is either 1 or -1)? (If understood correctly, the power iteration returns just one eigenvalue.) How do I actually calculate one eigenvalue?
The power method does not converge for your matrix.
From the wikipedia page:
Lambda_1 and lambda_2 are the two highest absolute value eigenvalues. In your case they are 1 and -1 so the convergence ratio is |1/-1| = 1. In other words the error stays the same at each iteration so the power method does not work.
Another way of understanding this is that your matrix takes a pair (a,b) and reverses it to become (b,a). The answer you get will simply depend on whether you do an even or odd number of iterations.