Diagonalization and matrix powers

Repeating the same rule over and over is more common than you might think. You predict tomorrow's weather from today's, then the day after from that prediction. You compute next month's active and dormant users from this month's, then the month after that. Game economies, population changes and web page rankings are all "multiplying a vector by the same matrix again and again", that is, Aⁿx.

The trouble starts when n is large. Multiplying by a matrix 100 times is a lot of work, and it is hard to tell in advance where the result is heading. The key to this lesson is to use the eigenvectors from Lesson 11 as a basis (Lesson 10). Seen from the eigenvector directions, the matrix is just "stretching along each axis", and stretching 100 times only takes raising a number to the 100th power.

This one idea also explains why the golden ratio comes out of the Fibonacci sequence, and why the weather settles into fixed proportions in the long run. At the end you will write power iteration yourself, a way to find an eigenvalue by nothing more than repeated multiplication.

12-1. In the eigenvector basis it is just stretching (8 min)

Let's look again at $A = \begin{bmatrix} 2 & 1 \\ 1 & 2 \end{bmatrix}$ from Lesson 11. Its eigenvectors were $\mathbf{e}_1 = (1, 1)$ (eigenvalue 3) and $\mathbf{e}_2 = (1, -1)$ (eigenvalue 1). The two directions are different, so together they can build every vector in the plane (Lesson 4). In other words, the two eigenvectors form a basis.

Written in this basis, $\mathbf{v} = (3, 1)$ is $\mathbf{v} = 2\,\mathbf{e}_1 + 1\,\mathbf{e}_2$. Now multiply by $A$. All $A$ does to each eigenvector is "multiply by $\lambda$", so the calculation is very simple.

The eigenvector coordinates $(2, 1)$ became $(6, 1)$: each coordinate was simply multiplied by its own eigenvalue. The directions do not mix. Split into three steps, it looks like this.

1. Change coordinates. Turn the standard coordinates $(3, 1)$ into the eigenvector coordinates $(2, 1)$. As you saw in Lesson 10, this means multiplying by $P^{-1}$, the inverse of the basis matrix $P$ (columns = eigenvectors).
2. Stretch along each axis. Multiply the eigenvector coordinates by $\lambda_1$ and $\lambda_2$ respectively. This means multiplying by a diagonal matrix $D$ with only the eigenvalues on its diagonal.
3. Change coordinates back. Turn the eigenvector coordinates back into standard coordinates. This means multiplying by $P$.

As you learned in Lesson 6, matrix products are read from right to left, so chaining the three steps gives this.

Writing a matrix this way is called diagonalization. It means "seen in the right coordinate system, this matrix is a diagonal matrix".

12-2. Powers become easy (6 min)

Write $A^2$ using the diagonalization, and $P^{-1}P = I$ meets in the middle and disappears.

The same happens however many times you multiply, so $A^n = PD^nP^{-1}$. And to raise a diagonal matrix to a power, you just raise each diagonal entry to that power.

In words: "switch to eigenvector coordinates, multiply each coordinate by $\lambda^n$, and switch back". Multiply $\mathbf{v} = 2\,\mathbf{e}_1 + 1\,\mathbf{e}_2$ from 12-1 by $A$ ten times and you get this.

Instead of ten matrix products, one $3^{10}$ does it. More importantly, you can see where the result is heading. $3^n$ grows fast while $1^n$ stays put, so the larger $n$ gets, the more the $\mathbf{e}_1$ term dominates. That is why the direction of $A^n\mathbf{v}$ gets closer and closer to $(1, 1)$.

In general, it works like this.

Not every matrix can be diagonalized. When there is only one eigenvector line, like the shear matrix in Lesson 11, you cannot build a basis, so there is no $P^{-1}$. When the eigenvalues are complex, like a rotation matrix's, the matrix can be diagonalized in the world of complex numbers, but in this course we will only read it as a picture: "there is a rotation mixed in, so it spins while growing or shrinking".

The lab below starts from the vector x₀ (the teal handle) and plots x₀, Ax₀, A²x₀, …, Aⁿx₀ as points in order, joined by a line. Fainter points are earlier steps and darker ones are later, and the crimson arrow is the last xₙ. The yellow dotted line is the eigenvector direction of the eigenvalue with the largest absolute value. Move the slider n from 0 to 20 and watch whether the points get pulled toward the yellow dotted line.

In 'Weather Markov chain' and 'Fibonacci' the numbers are small, so a note tells you the picture is drawn 3.5× larger. The numbers on the right and in the NumPy panel below are the real values.

12-3. Markov chains: what will the weather be tomorrow? (7 min)

Suppose the weather is only ever "sunny" or "cloudy", and it follows these rules.

If you write today's state as a probability vector $\mathbf{x}$ = (chance of sun, chance of clouds), you get tomorrow's state by multiplying by one matrix.

Column 1 of the matrix says where you go tomorrow starting from "sunny today", and column 2 says where you go starting from "cloudy today" (Lesson 5: the columns are where the basis vectors land). That is why each column adds up to 1. If today is certainly sunny, $\mathbf{x}_0 = (1, 0)$; tomorrow is $A\mathbf{x}_0 = (0.9, 0.1)$, and the day after is $(0.86, 0.14)$. A process like this, where "the next state depends only on the current state", is called a Markov chain.

Let's find the eigenvalues. The trace is $1.4$ and the determinant is $0.45 - 0.05 = 0.4$.

Reading it with the rule from 12-2, the $0.4^n$ direction quickly fades away and only the direction of eigenvalue 1 remains. Since that eigenvalue is exactly 1, a vector in that direction does not change when multiplied by $A$. This is called the steady state. The first row of $(A - I)\mathbf{v} = \mathbf{0}$ is $-0.1s + 0.5c = 0$, that is $s = 5c$. Scaling it so the probabilities add up to 1 gives this.

The lab's "In the long run it settles at 83.3% sunny : 16.7% cloudy" is exactly this value. At n = 5 it has nearly arrived at $(0.835, 0.165)$, and at n = 20 it agrees to seven decimal places.

The surprising part is that it does not depend on where you start. Whether today is certainly cloudy or a fifty-fifty, in the long run you arrive at the same 83.3% : 16.7%. Every difference in the starting point is carried in the $0.4^n$ direction and fades away. PageRank, which ranks web pages, is a steady state with the same structure: "the share of time a user randomly following links spends on each page in the long run".

If you drag x₀ in the lab to a vector whose components do not add up to 1, the point it converges to changes in proportion to the sum. For example, with a sum of 2 it goes to $(1.667, 0.333)$. It is still on the yellow dotted line, with the same 5 : 1 ratio of sun to clouds.

12-4. Fibonacci and the golden ratio, plus a spiral (5 min)

The Fibonacci sequence $1, 1, 2, 3, 5, 8, 13, \dots$ follows the rule "add the previous two numbers to get the next". Bundle two consecutive numbers into a vector, and this rule also becomes a single matrix multiplication.

The first row is $F_{n+1} = F_n + F_{n-1}$, and the second just copies $F_n$ down one slot. The eigenvalues of this matrix are the solutions of $\lambda^2 - \lambda - 1 = 0$.

$\lambda_1$ is the golden ratio. Since $\lvert\lambda_2\rvert < 1$, that direction fades away, and the vector $(F_{n+1}, F_n)$ lines up with the eigenvector direction of $\lambda_1$. That eigenvector points along $(\lambda_1, 1)$, so the ratio of the two components, $F_{n+1} / F_n$, approaches the golden ratio $1.61803\dots$. In the lab's 'Fibonacci', raise n and the "Ratio x_n[0] / x_n[1]" wobbles through 1, 2, 1.5, 1.667, 1.6, 1.625, … as it closes in on 1.61803. It wobbles up and down because $\lambda_2$ is negative, so the sign of the remaining error flips every step.

In 'Fibonacci', 'Show direction only (scale to length 1)' is on from the start. The components grow by a factor of $1.618^n$ and would fly off the screen, so each step is scaled to length 1 to show only the direction. Turn it off and the points race off the screen along the yellow dotted line in no time.

The 'Spiral inward' matrix $\begin{bmatrix} 0.8 & -0.5 \\ 0.5 & 0.8 \end{bmatrix}$ has the complex eigenvalues $0.8 \pm 0.5i$. As you saw in Lesson 11, complex eigenvalues signal a rotation. The size of the complex number, $\sqrt{0.8^2 + 0.5^2} \approx 0.943$, is less than 1, so each multiplication shrinks the length to about 94%. That is why the points circle the origin and get sucked in like a whirlpool. There are no real eigenvectors, so no yellow dotted line appears either.

12-5. Diagonalization and powers with NumPy (11 min)

Reading the lab's NumPy panel

On the lab's starting screen ('Weather Markov chain', n = 5), the NumPy panel shows this code.

>>> import numpy as np
>>> A = np.array([[0.9, 0.5], [0.1, 0.5]])
>>> x0 = np.array([1, 0])
>>> np.linalg.matrix_power(A, 5)
array([[0.83504, 0.8248 ],
       [0.16496, 0.1752 ]])
>>> np.linalg.matrix_power(A, 5) @ x0
array([0.83504, 0.16496])
>>> np.round(np.linalg.eigvals(A), 6)
array([1. +0.j, 0.4+0.j])

np.linalg.matrix_power(A, n) is the matrix Aⁿ, A multiplied n times. It has the same value as A @ A @ A @ A @ A, but it computes it with fewer multiplications by repeated squaring. With n = 0 you get the identity matrix. Add @ x0 and you have x₅, the weather probabilities five days later. Stepping forward one step at a time gives the same answer.

>>> A @ x0
array([0.9, 0.1])
>>> A @ (A @ x0)
array([0.86, 0.14])
>>> np.linalg.matrix_power(A, 2) @ x0
array([0.86, 0.14])
>>> np.linalg.matrix_power(A, 0)
array([[1., 0.],
       [0., 1.]])

As you saw in Lesson 11, np.linalg.eigvals(A) returns only the eigenvalues, and the panel wraps it in np.round(…, 6) to show six decimal places. So when you pick 'Fibonacci' you see only six digits, as in array([ 1.618034+0.j, -0.618034+0.j]). The space in 1. +0.j is only there to line it up with 0.4+0.j.

Raise n to 20 or 50 and it sticks to the steady state (5/6, 1/6). Starting from a cloudy day lands in the same place. Starting from a vector whose sum is 2 lands at the place with the same ratio, twice as far out.

>>> np.linalg.matrix_power(A, 20) @ x0
array([0.83333334, 0.16666666])
>>> np.linalg.matrix_power(A, 50) @ x0
array([0.83333333, 0.16666667])
>>> np.linalg.matrix_power(A, 50) @ np.array([0, 1])
array([0.83333333, 0.16666667])
>>> np.linalg.matrix_power(A, 30) @ np.array([2, 0])
array([1.66666667, 0.33333333])

np.diag: building a diagonal matrix and taking out the diagonal

np.diag does two different jobs depending on what you pass it. Pass a 1-D array (a list of numbers) and it builds a diagonal matrix with those numbers on the diagonal; pass a 2-D array (a matrix) and it takes out just the diagonal.

>>> np.diag([3, 1])
array([[3, 0],
       [0, 1]])
>>> np.diag(np.array([[1, 2], [3, 4]]))
array([1, 4])

Checking that P D P⁻¹ gives the matrix back

Let's write out the diagonalization from 12-1 exactly and check it. The columns of P are the eigenvectors (1, 1) and (1, −1).

>>> S = np.array([[2, 1], [1, 2]])
>>> P = np.array([[1, 1], [1, -1]])
>>> D = np.diag([3, 1])
>>> P @ D @ np.linalg.inv(P)
array([[2., 1.],
       [1., 2.]])
>>> S @ np.array([3, 1])
array([7, 5])

np.linalg.inv works in real numbers, so the result came out as floats like 2.. Now do the same with the eigenvectors NumPy finds. S is symmetric, so use eigh from Lesson 11.

>>> vals, vecs = np.linalg.eigh(S)
>>> P = vecs
>>> D = np.diag(vals)
>>> np.allclose(P @ D @ np.linalg.inv(P), S)
True

Now use Aⁿ = P Dⁿ P⁻¹ from 12-2. Only the diagonal needs raising to the power, so vals ** 10 does it. ** on an array raises each element to the power (the element-wise operations from Lesson 2).

>>> P @ np.diag(vals ** 10) @ np.linalg.inv(P)
array([[29525., 29524.],
       [29524., 29525.]])
>>> np.linalg.matrix_power(S, 10)
array([[29525, 29524],
       [29524, 29525]])
>>> np.linalg.matrix_power(S, 10) @ np.array([3, 1])
array([118099, 118097])

The (118099, 118097) you found by hand in 12-2 comes out exactly.

One more thing as a bonus. The eigenvectors of a symmetric matrix are perpendicular and have length 1, so the inverse of P equals its transpose P.T. A matrix like this is called an orthogonal matrix, and it is the star of the next lesson.

>>> np.allclose(np.linalg.inv(P), P.T)
True

Finding the steady state directly from the eigenvector

The steady state of a Markov chain is the eigenvector for eigenvalue 1. np.argmax tells you which slot holds the largest value. Take out the column with that index (Lesson 11: vecs[:, k]) and divide by the sum of its components so it becomes probabilities.

>>> vals, vecs = np.linalg.eig(A)
>>> k = np.argmax(vals.real)
>>> k
np.int64(0)
>>> s = vecs[:, k].real
>>> s
array([0.98058068, 0.19611614])
>>> s / s.sum()
array([0.83333333, 0.16666667])
>>> A @ (s / s.sum())
array([0.83333333, 0.16666667])

The eigenvector NumPy gives is (0.98, 0.196), scaled to length 1, so it does not look like probabilities, but the ratio is 5 : 1. Divide by the sum and you get (5/6, 1/6), which does not change when multiplied by A. That is the steady state.

You check that each column adds up to 1 with sum(axis=0). axis=0 means "add along axis 0 (down the rows), top to bottom", that is, add up each column. Lesson 15 sums up axes again.

>>> A.sum(axis=0)
array([1., 1.])

Fibonacci

>>> F = np.array([[1, 1], [1, 0]])
>>> np.linalg.matrix_power(F, 10)
array([[89, 55],
       [55, 34]])
>>> x = np.array([1, 0])
>>> for n in range(1, 8):
...     x = F @ x
...     print(n, x, x[0] / x[1])
1 [1 1] 1.0
2 [2 1] 2.0
3 [3 2] 1.5
4 [5 3] 1.6666666666666667
5 [8 5] 1.6
6 [13  8] 1.625
7 [21 13] 1.6153846153846154
>>> np.linalg.eigvals(F)
array([ 1.61803399+0.j, -0.61803399+0.j])
>>> (1 + np.sqrt(5)) / 2
np.float64(1.618033988749895)

The entries of F¹⁰ are the Fibonacci numbers 89, 55 and 34. The ratio in the for loop wobbling up and down as it approaches the golden ratio is exactly what 12-4 described.

12-6. Power iteration: finding an eigenvalue by repetition alone (5 min)

Turn "the direction of the largest eigenvalue wins in the end" from 12-2 around, and you can find the largest eigenvalue without solving the characteristic equation. Start from any vector and keep multiplying by A: that is power iteration.

But if you only multiply, the numbers either grow without end, as with Fibonacci, or shrink to 0. So you scale to length 1 each time (the same thing as the lab's 'Show direction only (scale to length 1)'). How many times longer the vector gets in one multiplication is an approximation of the eigenvalue.

>>> x = np.array([1.0, 0.0])
>>> for k in range(1, 11):
...     y = F @ x
...     lam = np.linalg.norm(y) / np.linalg.norm(x)
...     x = y / np.linalg.norm(y)
...     print(k, round(float(lam), 6))
1 1.414214
2 1.581139
3 1.612452
4 1.617215
5 1.617914
6 1.618017
7 1.618031
8 1.618034
9 1.618034
10 1.618034

round(float(lam), 6) turns the value into a plain Python number, so each line does not get np.float64(…) attached, and shows six decimal places. After eight rounds the golden ratio is correct to six places. At that point x is an approximation of the eigenvector.

>>> x
array([0.85067229, 0.52569636])
>>> np.round(x[0] / x[1], 4)
np.float64(1.6182)

How fast it converges is set by the ratio of the sizes of the two eigenvalues, |λ₂| / |λ₁|. For Fibonacci it is 0.618 / 1.618 ≈ 0.38, so each step shrinks the error to about 0.38 times, and it converges quickly. When the two eigenvalues are close in size, it gets very slow.

In practice np.linalg.eig finds all the eigenvalues at once with a much more sophisticated method (the QR algorithm). Even so, power iteration is still used today when a matrix is so large, millions by millions, that eig is out of the question, and PageRank was originally computed this way. Calculations like this, which "repeat to close in on the answer", are covered properly in numerical analysis, the course after this one.

12-7. Common mistakes (3 min)

A ** n is not a matrix power. ** on an array raises each element to the power. It is the same story as A * B not being matrix multiplication in Lesson 6.

>>> S ** 3
array([[8, 1],
       [1, 8]])
>>> np.linalg.matrix_power(S, 3)
array([[14, 13],
       [13, 14]])

Integer overflow happens silently. NumPy's int64 only holds numbers up to about 9.2 × 10¹⁸. The Fibonacci number F₁₀₁ is far bigger than that, and you get a nonsense number with no error and no warning.

>>> np.linalg.matrix_power(F, 90) @ np.array([1, 0])
array([4660046610375530309, 2880067194370816120])
>>> np.linalg.matrix_power(F, 100) @ np.array([1, 0])
array([1298777728820984005, 3736710778780434371])

Up to the 90th power it is correct, but in the 100th power result the second number is bigger than the first, so these are not Fibonacci numbers. If you need very large integers, dtype=object makes the array hold Python integers (slow, but exact).

>>> Fo = np.array([[1, 1], [1, 0]], dtype=object)
>>> np.linalg.matrix_power(Fo, 100) @ np.array([1, 0])
array([573147844013817084101, 354224848179261915075], dtype=object)

Matrices that cannot be diagonalized. For the shear matrix, the two columns of the P that eig returns lie on the same line (Lesson 11), so its determinant is close to 0. Compute its inverse and instead of an error you get huge numbers around 10¹⁵. Do not trust results like that.

>>> vals, vecs = np.linalg.eig(np.array([[1, 1], [0, 1]]))
>>> np.linalg.det(vecs)
np.complex128(2.2204460492503185e-16+0j)

Exercises (5 min)

1. Use the diagonalization from 12-1 to find $\begin{bmatrix} 2 & 1 \\ 1 & 2 \end{bmatrix}^5$ by hand. (Hint: $3^5 = 243$. The entries of $A^n$ are $\tfrac{3^n + 1}{2}$ and $\tfrac{3^n - 1}{2}$.)
2. Suppose app users are in one of two states, "active" and "dormant". Next month 80% of active users stay active, and 30% of dormant users come back. Set up the matrix $T$ and find what share of users is active in the long run.
3. Predict which eigenvalue and eigenvector power iteration from 12-6 approaches when run 30 times on $\begin{bmatrix} 2 & 1 \\ 1 & 2 \end{bmatrix}$, then check.

Solutions

Question 1: $A^5 = P\,\mathrm{diag}(243, 1)\,P^{-1}$, and the entries are $(243 + 1)/2 = 122$ and $(243 - 1)/2 = 121$.

Question 2: the columns say "where you go next month starting from this month's state", so $T = \begin{bmatrix} 0.8 & 0.3 \\ 0.2 & 0.7 \end{bmatrix}$. With trace 1.5 and determinant $0.56 - 0.06 = 0.5$, the eigenvalues are 1 and 0.5. The eigenvector for eigenvalue 1 comes from $-0.2a + 0.3b = 0$, so $a = 1.5b$, and scaling it to add up to 1 gives $(0.6, 0.4)$. In the long run, 60% are active and 40% dormant.

Question 3: it approaches the largest eigenvalue 3 and its eigenvector direction $(1, 1)$, that is $(0.707, 0.707)$ scaled to length 1.

>>> np.linalg.matrix_power(S, 5)
array([[122, 121],
       [121, 122]])
>>> T = np.array([[0.8, 0.3], [0.2, 0.7]])
>>> np.linalg.eigvals(T)
array([1. +0.j, 0.5+0.j])
>>> vals, vecs = np.linalg.eig(T)
>>> s = vecs[:, np.argmax(vals.real)].real
>>> s / s.sum()
array([0.6, 0.4])
>>> np.linalg.matrix_power(T, 30) @ np.array([0, 1])
array([0.6, 0.4])
>>> x = np.array([1.0, 0.0])
>>> for _ in range(30):
...     x = S @ x
...     x = x / np.linalg.norm(x)
>>> x
array([0.70710678, 0.70710678])
>>> x @ S @ x
np.float64(3.0000000000000004)

The last line for question 2 confirms that even starting with every user dormant, you arrive at (0.6, 0.4) after 30 months. In question 3, x @ S @ x measures "how many times S stretches x" for an x of length 1, and it came out as the eigenvalue 3 with a little floating-point residue attached. For an exact comparison, use np.isclose(x @ S @ x, 3).

In the next lesson you will see why orthogonal vectors are so convenient, and learn least squares, which finds the line that best fits a set of data points. It is a way to find "the closest answer" to a system of equations that has no exact solution.