Power iteration

2019-04-13 01:57发布

问题:

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?

回答1:

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.