Review

A symmetric matrix has real eigenvalues and perpendicular eigenvectors, which makes its eigendecomposition especially clean. More generally, whenever a matrix has a full set of eigenvectors we can write , and raising to a power only touches the diagonal core: . This chapter turns that fact into a surprisingly practical algorithm.

Power iteration (also called the power method) is a simple algorithm for finding the largest eigenvalue of a matrix and its corresponding eigenvector. It works by repeatedly applying the matrix to a vector, and it falls right out of what we already know about matrix powers.

---

17.1 Let's keep transforming a vector

Given a matrix and a starting vector , repeat these steps:

1. Normalize to unit length: .
2. Replace with the matrix-vector product: .
3. Go back to step 1.

Here is a small program that does exactly this:

Vector = [float]
Matrix = [[float]]

def scale_vector(scale: float, vector: Vector) -> Vector:
    return [i * scale for i in vector]

def add_vectors(array_of_vectors: Matrix) -> Vector:
    dimensionality = len(array_of_vectors[0])
    output = [0] * dimensionality
    for vector in array_of_vectors:
        output = [i + j for (i, j) in zip(output, vector)]
    return output

def get_matrix_vector_product(matrix: Matrix, vector: Vector) -> Vector:
    scaled_vectors = [scale_vector(scale, vector)
        for (scale, vector) in zip(vector, matrix)]
    return add_vectors(scaled_vectors)

def get_norm(vector: Vector) -> float:
    return (sum(i ** 2 for i in vector)) ** 0.5

def get_unit_vector(vector: Vector) -> Vector:
    norm = get_norm(vector)
    return [i / norm for i in vector]

vector = [1, 1]
matrix = [[1.5, 0.2], [0.5, 0.5]]
num_steps = 100

for _ in range(num_steps):
    vector = get_unit_vector(vector)
    vector = get_matrix_vector_product(matrix, vector)

---

17.1.1 When does it converge?

Let's watch the iteration run for two different matrices. The starting vector does not matter, so we begin with in both cases.

Start with

whose largest eigenvalue is .

First see how stretches the unit circle into an ellipse, then reveal the dominant eigenvector direction, the direction stretched the most.

Now run the iteration one step at a time. Start with , normalized. Then apply and renormalize, again, again, and once more. Each new arrow leans closer to the violet line, and after only a few steps has locked onto the dominant eigenvector: it converges to the direction .

---

Now try the scaled rotation

which has no real eigenvectors, since a rotation moves every direction.

See how maps the unit circle. Then iterate: start with , apply and renormalize, again, again, and again. The arrow turns by at every step and marches around the circle forever: it never converges.

So when does the iteration converge, and to what direction?

---

17.1.2 The answer

For an matrix , the iteration converges to the eigenvector belonging to the largest eigenvalue if all three of the following hold:

• is diagonalizable, meaning where has rank (its columns are the eigenvectors of ) and is diagonal (its entries are the eigenvalues of ).
• The largest absolute eigenvalue of is unique.
• The starting vector is not perpendicular to the eigenvector of that largest absolute eigenvalue.

The rotation fails the first two conditions (it has no real eigendecomposition), which is why its iteration never settles.

---

17.1.3 But why?

At each step we compute . This is similar to just repeatedly multiplying by , that is, computing where is the starting vector. Why similar? Because the loop

for _ in range(num_steps):
    vector = get_unit_vector(vector)
    vector = get_matrix_vector_product(matrix, vector)

differs from the plain power loop

for _ in range(num_steps):
    vector = get_matrix_vector_product(matrix, vector)

only by the renormalization at the start of each step. And scaling the input of a matrix only scales the output by the same factor: . It never changes the direction of the output.

So the normalized iteration and the raw power iteration always point in the same direction, differing only in length. If one converges in direction, so does the other. Now look at what a matrix power does through the eigendecomposition:

The key observation is that raising numbers to a power does not change which one is largest:

argmax(array) == argmax([i ** k for i in array])

Assume is the largest absolute eigenvalue, so is the largest entry of . As grows, pulls away from all the others, because each ratio (for ) blows up. So becomes dominated by its top-left entry:

which gives . Reading this from right to left:

• As grows, applying it to scales the first coordinate of by and crushes every other coordinate toward .
• So the middle vector keeps only its first entry: it looks like for some number .
• Multiplying that by picks out times the first column of , which is exactly the eigenvector of the largest absolute eigenvalue:

So : the raw power iteration lines up with the dominant eigenvector, and therefore so does the normalized power iteration. That is the direction converged to for matrix .

---

17.2 Finding the largest eigenvalue

This hands us an easy way to find the largest absolute eigenvalue and its eigenvector for any diagonalizable matrix : iterate until stops moving, and the amount by which then stretches is that eigenvalue.

It is especially valuable when is sparse, because a matrix-vector product with a sparse matrix is cheap, and power iteration only ever needs matrix-vector products, never the full eigendecomposition. According to this Wikipedia article, Google used exactly this idea to compute the dominant eigenvector of the enormous matrix representing links between web pages.

One of the most elegant consequences of these ideas is Binet's formula for Fibonacci numbers, which we turn to next.

In the next chapter we use eigendecomposition and matrix powers to derive a closed-form expression for the -th Fibonacci number, turning a recursive sequence into a single formula.

---

← 16. Eigenvectors of symmetric matrices · 18. Fibonacci numbers →