Power iteration

Question

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?

Solution

The power method does not converge for your matrix.

From the wikipedia page:

The convergence is geometric, with ratio |lambda_2 / lambda_1|

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.

Comments

[Ecir Hana]

Um, that's...unfortunate. Is there anything I can do about it, besides using other algorithm? How do I calculate this lambda_{1, 2} anyway? When I uses different matrix, say [[0.0, 2.0], [1.0, 0.0]], it still does not converge..?

[David Eisenstat]

@EcirHana That one has two eigenvalues of norm sqrt(2), so the ratio is again 1.

[Lutz Lehmann]

There is one thing you can do and that is to apply a matrix shift that will break the symmetry. This does obviously not work for multiple eigenvalues. So try iterating on [[1.0, 2.0], [1.0, 1.0]] and remember to subtract 1 from the eigenvalues of this matrix to compensate for the added identity matrix.