Eigenvectors

The last chapter ended with an observation: a linear transformation can have directions that it doesn’t turn but only stretches. Let’s check on the matrix $A = \left(\begin{smallmatrix} 2 & 1 \\ 1 & 2 \end{smallmatrix}\right)$. It takes the vector $(1, 0)$ to $(2, 1)$: the arrow has turned. The vector $(0, 1)$ goes to $(1, 2)$ and turns too. But the vector $(1, 1)$ goes to $(3, 3)$: the same direction, three times the length. And the vector $(1, -1)$ stays exactly where it is.

Such directions will be our quarry. They’re worth hunting, because they hide the answers to three problems that at first sight have nothing in common. How do you find the hundredth Fibonacci number without computing the ninety-nine before it? How does a search engine decide which of billions of pages is more important? And what is the curve $5x^2 + 4xy + 2y^2 = 6$, with the very product $xy$ that the chapter on conic sections left for us to sort out?

The answer in all three cases is the same. Find the directions the matrix doesn’t turn, and look at it along them. Along these directions even the most tangled transformation turns out to be a simple stretch.

Tracks on the grid

Take a $2 \times 2$ matrix and start turning a vector $\mathbf v$ of length one. Its image $A\mathbf v$ turns too, but at its own speed: in places it lags behind $\mathbf v$, in places it runs ahead. If in one position $A\mathbf v$ lags and in another it overtakes, then somewhere in between the arrows are bound to lie on one line. Those are the positions we’re trying to catch.

A square matrix: the transformation we’re studying. The eigenvector. It has to be non-zero: the zero vector satisfies the equation for any $\lambda$ and says nothing about the matrix. The eigenvalue: the factor by which the vector is stretched. It can be negative (the vector turns round but stays on its line) or zero (the vector is squashed to a point). Example: for $A = \left(\begin{smallmatrix} 2 & 1 \\ 1 & 2 \end{smallmatrix}\right)$ and $\mathbf v = (1, 1)$ we get $A\mathbf v = (3, 3) = 3\mathbf v$, that is, $\lambda = 3$. For $\mathbf v = (1, -1)$ it comes out as $\lambda = 1$.

As in the previous chapters, vectors are written in bold, and coordinates in the text are written in a row: $(1, -1)$ means the column $\left(\begin{smallmatrix} 1 \\ -1 \end{smallmatrix}\right)$.

If $\mathbf v$ is an eigenvector, so are $2\mathbf v$, $-\mathbf v$ and in general $c\mathbf v$ for any $c \ne 0$: $A(c\mathbf v) = c\,A\mathbf v = \lambda\,(c\mathbf v)$. So what we’re hunting are lines, not individual arrows. The sum of two eigenvectors with the same $\lambda$ again satisfies $A\mathbf v = \lambda\mathbf v$, so all such vectors together with the zero vector form a subspace, the eigenspace. Most often it’s a line, but it can be bigger: for the scaling $\mathbf v \mapsto 2\mathbf v$ the whole space is an eigenspace.

Here’s what the hunt gives for the familiar transformations from chapter 35:

Next to the plane in the widget a sonar is running. It shows the signed area of the parallelogram on $\mathbf v$ and $A\mathbf v$: plus if $A\mathbf v$ is turned anticlockwise from $\mathbf v$, minus if clockwise. Where the curve crosses zero, the vectors lie on one line. For a rotation the whole curve is above zero, and the quarry escapes. For a shear it only touches zero at one point.

Why the sonar draws a sine wave

For $A = \left(\begin{smallmatrix} a & b \\ c & d \end{smallmatrix}\right)$ and $\mathbf v = (\cos\theta, \sin\theta)$ the signed area equals

$$S(\theta) = \det(\mathbf v, A\mathbf v) = \frac{c - b}{2} + \frac{c + b}{2}\cos 2\theta + \frac{d - a}{2}\sin 2\theta.$$

This is a sine wave with period $180°$, raised by $\frac{c - b}{2}$. The rise is the contribution of the “rotational” part of the matrix, and the amplitude of the sine wave is $\frac12\sqrt{(c + b)^2 + (d - a)^2}$. It reaches zero when the rise is no bigger than the amplitude: $(c - b)^2 \le (c + b)^2 + (d - a)^2$, that is, $(a - d)^2 + 4bc \ge 0$. On the left is the discriminant of the characteristic polynomial from the next section. For a symmetric matrix $b = c$, there’s no rise, and there are always two eigendirections (or all of them at once). We’ll come back to this at the end of the chapter.

How many eigendirections does a rotation of the plane through $180°$ about the origin have?

None: it’s a rotation, after all Two perpendicular ones All directions One

A rotation through $180°$ takes every vector $\mathbf v$ to $-\mathbf v$, that is, it coincides with the matrix $-I$. Every non-zero vector is an eigenvector, with $\lambda = -1$: the arrow turns round but stays on its line.

The trap: the characteristic polynomial

Walking an arrow round a circle is possible in the plane. In three-dimensional space the directions form a whole sphere, and where there are billions of coordinates (a search engine works with vectors like that), trying them one by one is hopeless. We need a trap that catches all the quarry at once.

Let’s rewrite the definition. The equation $A\mathbf v = \lambda\mathbf v$ means the same as $A\mathbf v - \lambda I\mathbf v = \mathbf 0$, that is,

$$(A - \lambda I)\,\mathbf v = \mathbf 0.$$

Here $I$ is the identity matrix (some books write $E$); subtracting $\lambda I$ means subtracting $\lambda$ from each diagonal entry. We’ve got a homogeneous system of equations, and an eigenvector is a non-zero solution of it. Whether there is such a solution, the determinant will tell us.

An $n \times n$ matrix whose eigenvalues we’re looking for. An unknown number. Each $\lambda$ gives its own matrix $A - \lambda I$, and we’re interested in those that squash space. The identity matrix of the same size: $\lambda I$ has $\lambda$ on the diagonal and zeros elsewhere. Example: for $A = \left(\begin{smallmatrix} 4 & 1 \\ 2 & 3 \end{smallmatrix}\right)$ we get $\det\left(\begin{smallmatrix} 4 - \lambda & 1 \\ 2 & 3 - \lambda \end{smallmatrix}\right) = (4 - \lambda)(3 - \lambda) - 2 = \lambda^2 - 7\lambda + 10 = 0$, so $\lambda_1 = 2$, $\lambda_2 = 5$.

The idea: the equation $A\mathbf v = \lambda\mathbf v$ means that the matrix $A - \lambda I$ sends a non-zero vector to zero. Such a matrix squashes the plane onto a line, and the determinant shows how area changes, so it’s zero. The picture is for $2 \times 2$ matrices; what changes in $n$ dimensions is said at the end.

Let $\lambda$ be an eigenvalue: $A\mathbf v = \lambda\mathbf v$ for some $\mathbf v \ne \mathbf 0$. In the picture $\mathbf v$ and $A\mathbf v$ lie on one line. Since $I\mathbf v = \mathbf v$, the right side can be written as $\lambda I\mathbf v$, and by the distributive law for matrices $A\mathbf v - \lambda I\mathbf v = (A - \lambda I)\mathbf v$. So $(A - \lambda I)\mathbf v = A\mathbf v - \lambda\mathbf v = \mathbf 0$: the matrix $B = A - \lambda I$ takes the non-zero vector $\mathbf v$ to zero. Hence $\det B = 0$. Extend $\mathbf v$ by some vector $\mathbf w$ to a basis of the plane. Any vector equals $s\mathbf v + t\mathbf w$, and by linearity $B(s\mathbf v + t\mathbf w) = sB\mathbf v + tB\mathbf w = tB\mathbf w$. Everything $B$ produces lies on one line, the line of the vector $B\mathbf w$, so the unit square goes to a segment of zero area. The area of the image of the square is $|\det B|$ (chapter 35), so $\det B = 0$. In the picture you can see how the columns $B\mathbf e_1 = A\mathbf e_1 - \lambda\mathbf e_1$ and $B\mathbf e_2 = A\mathbf e_2 - \lambda\mathbf e_2$ are obtained by moving the tips of the columns of $A$ by $\lambda$ along the axes, and fall onto one line. Conversely, let $\det B = 0$. Then the parallelogram on the vectors $B\mathbf e_1$ and $B\mathbf e_2$ has zero area, that is, they lie on one line. If $B\mathbf e_1 = \mathbf 0$, take $\mathbf v = \mathbf e_1$. Otherwise $B\mathbf e_2 = kB\mathbf e_1$ for some number $k$, and the vector $\mathbf v = (k, -1)$ gives $B\mathbf v = kB\mathbf e_1 - B\mathbf e_2 = \mathbf 0$. In both cases we’ve found $\mathbf v = (v_1, v_2) \ne \mathbf 0$ with $v_1B\mathbf e_1 + v_2B\mathbf e_2 = B\mathbf v = \mathbf 0$: in the picture two arrows come back to the origin. Reading the rearrangement from the second step backwards, from $(A - \lambda I)\mathbf v = \mathbf 0$ we get $A\mathbf v = \lambda\mathbf v$ with $\mathbf v \ne \mathbf 0$, that is, $\lambda$ is an eigenvalue. In $n$ dimensions the argument is the same, only with volume instead of area, and “the columns are linearly dependent” instead of “they lie on one line”. That $\det B = 0$ exactly when the equation $B\mathbf v = \mathbf 0$ has a non-zero solution is a theorem from chapter 36; for complex numbers Gaussian elimination works just the same.

For an $n \times n$ matrix the characteristic polynomial has degree $n$: when the determinant is expanded, the $\lambda$’s from the diagonal entries are multiplied together at most $n$ times, and the leading term is $(-\lambda)^n$. So there are at most $n$ distinct eigenvalues. Some books write $\det(\lambda I - A)$ so that the leading coefficient is one; the roots don’t change.

Once the eigenvalues are found, the eigenvectors come from solving a homogeneous system by Gaussian elimination from chapter 36.

For a $2 \times 2$ matrix the determinant can be expanded in general: $\det(A - \lambda I) = (a - \lambda)(d - \lambda) - bc = \lambda^2 - (a + d)\lambda + (ad - bc)$. The constant term is the determinant of $A$, and the coefficient of $\lambda$ is minus the sum of the diagonal entries. This sum is called the trace of the matrix and written $\operatorname{tr} A$.

The trace $a + d$. By Vieta’s formulas from chapter 10 it equals the sum of the roots: $\lambda_1 + \lambda_2 = \operatorname{tr} A$. The determinant $ad - bc$ equals the product of the roots: $\lambda_1\lambda_2 = \det A$. Example: $A = \left(\begin{smallmatrix} 4 & 1 \\ 2 & 3 \end{smallmatrix}\right)$ has trace $7$ and determinant $10$; the roots are $2$ and $5$, and indeed $2 + 5 = 7$, $2 \cdot 5 = 10$. If the discriminant $(\operatorname{tr} A)^2 - 4\det A$ is negative, there are no real eigenvalues.

Both equations hold for matrices of any size too.

Hence a useful test: a matrix is singular exactly when zero is among its eigenvalues. The whole set of eigenvalues is called the spectrum of the matrix.

Find the eigenvalues of the matrix $\left(\begin{smallmatrix} 1 & 2 \\ 3 & 2 \end{smallmatrix}\right)$.

The trace is $1 + 2 = 3$, the determinant $1 \cdot 2 - 2 \cdot 3 = -4$. The equation $\lambda^2 - 3\lambda - 4 = 0$ factors as $(\lambda - 4)(\lambda + 1) = 0$, so $\lambda_1 = 4$, $\lambda_2 = -1$. Check: $4 + (-1) = 3$ and $4 \cdot (-1) = -4$.

The trophy: a diagonal matrix

Suppose the hunt was a success: an $n \times n$ matrix turned out to have $n$ linearly independent eigenvectors $\mathbf v_1, \dots, \mathbf v_n$ with eigenvalues $\lambda_1, \dots, \lambda_n$. They form a basis, and in this basis the matrix is as simple as it gets. The vector with coordinates $c_1, \dots, c_n$ in the new basis is $c_1\mathbf v_1 + \dots + c_n\mathbf v_n$, and the matrix turns it into $\lambda_1c_1\mathbf v_1 + \dots + \lambda_nc_n\mathbf v_n$. Each coordinate is multiplied by its own number, and nothing gets mixed.

Let’s write this with matrices. Put the eigenvectors as columns of a matrix $P$ (in chapter 37 such a matrix was called the change-of-basis matrix), and the eigenvalues on the diagonal of the matrix $D = \operatorname{diag}(\lambda_1, \dots, \lambda_n)$, which has zeros off the diagonal. Then

The return to the original coordinates: the columns of $P$ are the eigenvectors themselves. The diagonal matrix of eigenvalues: in eigencoordinates $A$ simply stretches the $i$th axis by a factor of $\lambda_i$. The move to eigencoordinates: $P^{-1}\mathbf x$ is the coordinates of the vector $\mathbf x$ in the basis of eigenvectors. The product is read from right to left, like a composition of transformations. Example: $\left(\begin{smallmatrix} 4 & 1 \\ 2 & 3 \end{smallmatrix}\right) = \left(\begin{smallmatrix} 1 & 1 \\ -2 & 1 \end{smallmatrix}\right)\left(\begin{smallmatrix} 2 & 0 \\ 0 & 5 \end{smallmatrix}\right)\cdot\frac13\left(\begin{smallmatrix} 1 & -1 \\ 2 & 1 \end{smallmatrix}\right)$. The columns of the first matrix are the vectors $(1, -2)$ and $(1, 1)$ found above.

The idea: in a basis of eigenvectors the matrix stretches each coordinate separately, and $P^{-1}$ and $P$ are merely the move to this basis and back. The picture is for $n = 2$, but not a single step uses the fact that there are two coordinates.

Let $A\mathbf p_i = d_i\mathbf p_i$ for all $i$. In the picture $\mathbf p_1 = \mathbf v_1$, $\mathbf p_2 = \mathbf v_2$, $d_1 = 2$, $d_2 = \frac12$: on the basis vectors the matrix acts most simply, stretching $\mathbf v_1$ twofold and shrinking $\mathbf v_2$ by half. The skewed grid is the coordinate grid of this basis. Take any vector $\mathbf x$ and expand it in the basis: $\mathbf x = c_1\mathbf v_1 + c_2\mathbf v_2$. This equation means $\mathbf x = P\mathbf c$ for the column $\mathbf c = (c_1, c_2)$ (multiplying a matrix by a column means adding its columns with these coefficients). The matrix $P$ is invertible, so $\mathbf c = P^{-1}\mathbf x$: $P^{-1}$ turns a vector into its coordinates in the eigenbasis. By linearity $A\mathbf x = c_1A\mathbf v_1 + c_2A\mathbf v_2 = d_1c_1\mathbf v_1 + d_2c_2\mathbf v_2$. Each coordinate has simply been multiplied by its own number, and the new column of coordinates is $(d_1c_1, d_2c_2) = D\mathbf c = DP^{-1}\mathbf x$. Assemble the vector from the new coordinates, which is again multiplication by $P$: $A\mathbf x = P(DP^{-1}\mathbf x) = (PDP^{-1})\mathbf x$ (matrix multiplication is associative). The equation holds for every $\mathbf x$, in particular for $\mathbf e_1, \dots, \mathbf e_n$, and the columns of a matrix are the images of the $\mathbf e_i$. So the matrices $A$ and $PDP^{-1}$ have the same columns: $A = PDP^{-1}$. Conversely, let $A = PDP^{-1}$. Multiply both sides on the right by $P$: $AP = PD$. On the left the $i$th column is $A\mathbf p_i$, on the right $d_i\mathbf p_i$, because multiplying by a diagonal matrix on the right multiplies the $i$th column by $d_i$. So $A\mathbf p_i = d_i\mathbf p_i$; the columns of an invertible matrix are non-zero and linearly independent, so $\mathbf p_1, \dots, \mathbf p_n$ is a basis of eigenvectors. The whole parallelogram of the grid on $\mathbf v_1$, $\mathbf v_2$ goes to the parallelogram on $2\mathbf v_1$, $\frac12\mathbf v_2$.

The trophy justifies the hunt when the matrix has to be applied many times. In the product $(PDP^{-1})(PDP^{-1})\cdots(PDP^{-1})$ all the inner pairs $P^{-1}P$ cancel.

For our example

$$A^n = \frac13\begin{pmatrix} 2^n + 2\cdot 5^n & 5^n - 2^n \\ 2\cdot 5^n - 2^{n+1} & 5^n + 2^{n+1} \end{pmatrix}.$$

At $n = 1$ we get $A$ itself, and at $n = 10$ the matrix $\left(\begin{smallmatrix} 6{,}510{,}758 & 3{,}254{,}867 \\ 6{,}509{,}734 & 3{,}255{,}891 \end{smallmatrix}\right)$, without having to multiply $A$ by itself nine times. Something else shows too: for large $n$ the terms with $5^n$ swamp everything else, so $A^n$ is almost proportional to the matrix $\left(\begin{smallmatrix} 2 & 1 \\ 2 & 1 \end{smallmatrix}\right)$, both of whose columns point along $(1, 1)$. We’ll come back to this observation.

Find the top left entry of the matrix $A^5$ for $A = \left(\begin{smallmatrix} 4 & 1 \\ 2 & 3 \end{smallmatrix}\right)$.

By the formula for $A^n$: $\frac{2^5 + 2 \cdot 5^5}{3} = \frac{32 + 6250}{3} = 2094$.

When is a basis of eigenvectors guaranteed? A sufficient condition is simple.

Who escapes the hunter

The hunt doesn’t always succeed, and failures come in two kinds. A rotation escapes the hunter entirely, while a shear leaves less quarry than a basis needs.

Rotation: quarry in the complex plane

Rotation through the angle $\varphi$ has the matrix $\left(\begin{smallmatrix} \cos\varphi & -\sin\varphi \\ \sin\varphi & \cos\varphi \end{smallmatrix}\right)$, trace $2\cos\varphi$ and determinant $1$. The equation $\lambda^2 - 2\lambda\cos\varphi + 1 = 0$ has discriminant $4\cos^2\varphi - 4 = -4\sin^2\varphi$, negative for every angle except $0°$ and $180°$. The roots are complex:

$$\lambda = \cos\varphi \pm i\sin\varphi = e^{\pm i\varphi}.$$

In the chapter on complex numbers we saw that multiplication by $a + bi$ is a rotation with a stretch. Since $(a + bi)(x + yi) = (ax - by) + (bx + ay)i$, in the coordinates $(x, y)$ it’s written with the matrix $\left(\begin{smallmatrix} a & -b \\ b & a \end{smallmatrix}\right)$, and its eigenvalues are $a \pm bi$. It turns out that every real $2 \times 2$ matrix with complex eigenvalues looks like this in a suitable basis.

Example: $A = \left(\begin{smallmatrix} 1 & -1 \\ 1 & 1 \end{smallmatrix}\right)$ has eigenvalues $1 \pm i = \sqrt2\,e^{\pm i\pi/4}$. So $A$ rotates through $45°$ and stretches by a factor of $\sqrt2$, and eight such steps make a full turn: $A^8 = (\sqrt2)^8I = 16I$. There are eigenvectors too, only complex ones: $A\left(\begin{smallmatrix} 1 \\ -i \end{smallmatrix}\right) = \left(\begin{smallmatrix} 1 + i \\ 1 - i \end{smallmatrix}\right) = (1 + i)\left(\begin{smallmatrix} 1 \\ -i \end{smallmatrix}\right)$.

Over the complex numbers the hunter never goes home empty-handed: by the fundamental theorem of algebra the characteristic polynomial always has a root, so every square matrix has at least one eigenvalue and an eigenvector in $\mathbb C^n$.

Why every rotation of a ball has an axis

Let $A$ be a linear transformation of three-dimensional space that preserves the lengths of vectors and orientation ($\det A = 1$). Then there’s a non-zero vector $\mathbf n$ with $A\mathbf n = \mathbf n$: the transformation is a rotation about the axis $\mathbf n$.

The idea: work out what the eigenvalues can be, and see that $1$ must be among them. The characteristic polynomial of a $3 \times 3$ matrix has degree $3$ and real coefficients, and such a polynomial always has a real root: for large $|\lambda|$ it takes values of opposite signs, and a continuous function that changes sign passes through zero (the intermediate value theorem, chapter 25; a rigorous proof is in chapter 53). The other two roots are either real or a conjugate pair $\mu$, $\bar\mu$.

All the eigenvalues have absolute value $1$. For a real eigenvector this follows directly from $\|A\mathbf v\| = \|\mathbf v\|$ and $A\mathbf v = \lambda\mathbf v$. For a complex $\mathbf w = \mathbf u + i\mathbf v$ put $\|\mathbf w\|^2 = \|\mathbf u\|^2 + \|\mathbf v\|^2$; since $A\mathbf w = A\mathbf u + iA\mathbf v$ and $A$ preserves the lengths of $\mathbf u$ and $\mathbf v$, we get $\|A\mathbf w\| = \|\mathbf w\|$, while $A\mathbf w = \mu\mathbf w$ gives $\|A\mathbf w\| = |\mu|\,\|\mathbf w\|$, that is, $|\mu| = 1$.

The product of the three eigenvalues equals $\det A = 1$ (the statement about trace and determinant above). If there’s a conjugate pair, then $\mu\bar\mu = |\mu|^2 = 1$, and the third, real value is also $1$. If all three are real, they’re the numbers $\pm1$ with product $1$: either $1, 1, 1$ or $1, -1, -1$, and in both cases there’s a $1$. The eigenvector $\mathbf n$ with $\lambda = 1$ stays put, and the line through it is the axis. That’s how the modern language proves the theorem Euler obtained in 1775: any displacement of a rigid body with a fixed point is a rotation about some axis.

Shear: less quarry than needed

The shear $\left(\begin{smallmatrix} 1 & 1 \\ 0 & 1 \end{smallmatrix}\right)$ behaves differently. Its characteristic polynomial $(1 - \lambda)^2$ has the double root $\lambda = 1$, but the system $(A - I)\mathbf v = \mathbf 0$, that is, $\left(\begin{smallmatrix} 0 & 1 \\ 0 & 0 \end{smallmatrix}\right)\mathbf v = \mathbf 0$, leaves only the vectors with $y = 0$: one line. The root is counted twice, but there’s only one direction. There’s no basis of eigenvectors, and the matrix isn’t diagonalisable.

Here two notions of multiplicity part ways. The algebraic multiplicity of an eigenvalue is how many times it occurs among the roots of the characteristic polynomial. The geometric multiplicity is the dimension of its eigenspace. For the shear the algebraic multiplicity is $2$ and the geometric $1$. The geometric multiplicity is never larger than the algebraic one.

If the entries of a matrix are chosen at random, multiple roots almost certainly won’t occur, so non-diagonalisable matrices are the exception. But they turn up precisely where a system’s behaviour changes, for example on the boundary between oscillation and smooth decay (see the note below).

The matrices $\left(\begin{smallmatrix} 3 & 0 \\ 0 & 3 \end{smallmatrix}\right)$ and $\left(\begin{smallmatrix} 3 & 1 \\ 0 & 3 \end{smallmatrix}\right)$ have the same characteristic polynomial, $(\lambda - 3)^2$. Which of them is diagonalisable?

Both: the polynomial is the same Only the first Only the second Neither: the root is double

The first is already diagonal, and its eigenspace is the whole space. For the second the system $(A - 3I)\mathbf v = \mathbf 0$ gives one line, as for the shear. The characteristic polynomial doesn’t determine the matrix, and a double root isn’t yet a death sentence.

Fibonacci’s rabbits

The Fibonacci numbers $0, 1, 1, 2, 3, 5, 8, 13, \dots$, going back to the rabbit problem in the Liber Abaci (1202), obey the rule $F_{n+1} = F_n + F_{n-1}$ (chapter 13). To get the next number you have to remember the last two. Put them in a vector, and a step of the sequence becomes multiplication by a matrix:

$$\begin{pmatrix} F_{n+1} \\ F_n \end{pmatrix} = \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix}\begin{pmatrix} F_n \\ F_{n-1} \end{pmatrix}.$$

Call this matrix $M$. Starting from $(F_1, F_0) = (1, 0)$ and multiplying $n$ times, we get $(F_{n+1}, F_n) = M^n(1, 0)$ (at $n = 0$ this holds if we take $M^0 = I$). Along the way we get a pretty formula for the powers themselves.

Let’s hunt. The trace of $M$ is $1$, the determinant $-1$, and the characteristic equation $\lambda^2 - \lambda - 1 = 0$. Its roots are the golden ratio and its “twin”:

$$\varphi = \frac{1 + \sqrt5}{2} \approx 1.618, \qquad \psi = \frac{1 - \sqrt5}{2} \approx -0.618.$$

The eigenvectors are $(\varphi, 1)$ and $(\psi, 1)$. Expanding the starting vector in them, we get the famous formula.

The idea: expand the starting vector along the two eigendirections. After that each multiplication by $M$ simply stretches both parts, each by its own factor.

By the lemma the vector $(F_{n+1}, F_n)$ is obtained from $(1, 0)$ by multiplying by $M$ $n$ times. The picture shows these vectors for $n = 0, 1, 2, 3, 4$: $(1, 0)$, $(1, 1)$, $(2, 1)$, $(3, 2)$, $(5, 3)$. They hug one line ever more closely. $M$ has two eigendirections: $(\varphi, 1)$ with $\lambda = \varphi$ and $(\psi, 1)$ with $\lambda = \psi$. Let’s check the first: $M(\varphi, 1) = (\varphi + 1, \varphi)$, and this is $\varphi\,(\varphi, 1)$, because $\varphi^2 = \varphi + 1$; that’s how $\varphi$ is defined, as a root of the equation $\lambda^2 - \lambda - 1 = 0$. For $\psi$, the second root, the check is the same. Expand the starting vector: $(1, 0) = \frac{1}{\sqrt5}(\varphi, 1) - \frac{1}{\sqrt5}(\psi, 1)$. Check: the first coordinate is $\frac{\varphi - \psi}{\sqrt5} = \frac{\sqrt5}{\sqrt5} = 1$, the second $\frac{1 - 1}{\sqrt5} = 0$. Multiplication by $M$ is linear, so the parts can be multiplied separately. An eigenvector is only stretched: $M\mathbf v = \lambda\mathbf v$ gives $M^2\mathbf v = \lambda M\mathbf v = \lambda^2\mathbf v$ and so on, $M^n\mathbf v = \lambda^n\mathbf v$. So $(F_{n+1}, F_n) = \frac{\varphi^n}{\sqrt5}(\varphi, 1) - \frac{\psi^n}{\sqrt5}(\psi, 1)$. In the picture the first part grows along its line by a factor of $\varphi \approx 1.618$ per step, while the second, multiplied by $\psi \approx -0.618$, melts away and jumps to the other side each time. Compare the second coordinates: $F_n = \frac{\varphi^n}{\sqrt5} \cdot 1 - \frac{\psi^n}{\sqrt5} \cdot 1 = \frac{\varphi^n - \psi^n}{\sqrt5}$. In the picture these are the heights of the points: $0, 1, 1, 2, 3$.

The growing component: the eigenvalue $\varphi \approx 1.618$ to the power $n$. The fading component: $|\psi| \approx 0.618 < 1$, so $\psi^n \to 0$, changing sign at every step. A normalisation that came from expanding the starting vector: $\sqrt5 = \varphi - \psi$. Example: $\varphi^{10} \approx 122.9919$ and $\psi^{10} \approx 0.0081$, so $F_{10} \approx \frac{122.9919 - 0.0081}{2.2361} \approx 55.000$. Since $|\psi^n| / \sqrt5 < \frac12$ for every $n$, the number $F_n$ is $\frac{\varphi^n}{\sqrt5}$ rounded to the nearest integer: $\frac{\varphi^{10}}{\sqrt5} \approx 55.004$.

The hundredth Fibonacci number no longer has to be dug out with ninety-nine additions: $F_{100}$ is the integer nearest to $\varphi^{100}/\sqrt5$, that is, $354{,}224{,}848{,}179{,}261{,}915{,}075$. True, for the rounding to be honest, $\varphi^{100}$ has to be computed to more than twenty digits. The formula is named after Jacques Binet, who published it in 1843, although Abraham de Moivre and Daniel Bernoulli knew it more than a hundred years before him.

The matrix view also gives unexpected identities.

This identity was found by Giovanni Domenico Cassini in 1680. An old trick rests on it: an $8 \times 8$ square is cut into four pieces and reassembled into a $5 \times 13$ rectangle, and the area seems to grow from $64$ to $65$. The extra cell is smeared along a narrow, almost invisible gap along the diagonal of the rectangle: $5 \cdot 13 - 8^2 = 1$.

The Lucas numbers $2, 1, 3, 4, 7, 11, \dots$ obey the same rule but start with $L_0 = 2$, $L_1 = 1$. Their starting vector expands in the eigenvectors even more simply, and $L_n = \varphi^n + \psi^n$. Find $L_{10}$.

By addition: $2, 1, 3, 4, 7, 11, 18, 29, 47, 76, 123$. By the formula: $\varphi^{10} + \psi^{10} \approx 122.9919 + 0.0081 = 123$. This also shows why $\varphi^{10}$ is so close to an integer: only the tiny $\psi^{10}$ separates it from $L_{10}$.

Hunting without a trap: the power method

Look at the ratios of neighbouring Fibonacci numbers: $\frac21 = 2$, $\frac32 = 1.5$, $\frac53 \approx 1.667$, $\frac85 = 1.6$, $\frac{13}{8} = 1.625$, $\frac{21}{13} \approx 1.615$, $\frac{34}{21} \approx 1.619$. They converge to $\varphi$, and Binet’s formula explains why. In the vector $M^n(1, 0)$ the component along $(\psi, 1)$ shrinks compared with the main one as $(\psi/\varphi)^n \approx (-0.382)^n$, and the vector lies ever more precisely on the line $(\varphi, 1)$.

This observation works for any diagonalisable matrix with one eigenvalue strictly larger in absolute value than the rest. Expand the starting vector in eigenvectors: $\mathbf x = c_1\mathbf v_1 + c_2\mathbf v_2 + \dots$, where $|\lambda_1| > |\lambda_2| \ge \dots$. After $k$ multiplications

A common factor: the vector grows or shrinks, but its direction doesn’t depend on that. In practice the vector is divided by its length after each step, so that the numbers don’t overflow. The error: the share of the second eigenvector is multiplied by $\lambda_2/\lambda_1$ at each step. We need $|\lambda_1| > |\lambda_2|$ and $c_1 \ne 0$: the starting vector has to have at least some component along $\mathbf v_1$. Example: for the Fibonacci matrix $|\psi/\varphi| \approx 0.382$, and in ten steps the error shrinks by a factor of about $1/0.382^{10} \approx 15{,}000$.

That’s how we get the power method: take almost any vector, multiply by $A$ again and again, normalising after each step, and its direction converges to the eigenvector with the eigenvalue of largest absolute value. $\lambda_1$ itself is obtained at the end as $\langle \mathbf x, A\mathbf x\rangle$ for a unit $\mathbf x$. The method was described by Richard von Mises and Hilda Pollaczek-Geiringer in 1929. It needs neither determinants nor polynomials, only multiplication of a matrix by a vector, so it works even for matrices with billions of rows.

The idea: in eigencoordinates multiplication by $A^k$ multiplies $c_i$ by $\lambda_i^k$, and compared with the first, all the other parts melt away as $|\lambda_i/\lambda_1|^k$. In the picture $n = 2$, $\lambda_1 = 2$, $\lambda_2 = 1$, so the ratio of the parts halves at each step.

Expand $\mathbf x$ in the eigenbasis: $\mathbf x = c_1\mathbf v_1 + \dots + c_n\mathbf v_n$; in the picture the parts $c_1\mathbf v_1$ and $c_2\mathbf v_2$ are the sides of a parallelogram. The number $\lambda_1$ is real: the non-real eigenvalues of a real matrix come in conjugate pairs $\mu$, $\bar\mu$ with the same absolute value, while $|\lambda_1|$ is strictly larger than the absolute values of all the others. Multiplication by $A$ multiplies each part by its eigenvalue: $A\mathbf x = \lambda_1c_1\mathbf v_1 + \dots + \lambda_nc_n\mathbf v_n$. After $k$ steps $A^k\mathbf x = \lambda_1^k\,(c_1\mathbf v_1 + \mathbf r_k)$, where $\mathbf r_k = \sum_{i \ge 2} c_i\,(\lambda_i/\lambda_1)^k\,\mathbf v_i$. The normalisation $\mathbf x_k = A^k\mathbf x/\|A^k\mathbf x\|$ changes the length but not the direction; in the picture the normalised vectors lie on a circle. The remainder is small: by the triangle inequality $\|\mathbf r_k\| \le \sum_{i \ge 2}|c_i|\,|\lambda_i/\lambda_1|^k\,\|\mathbf v_i\| \le K\,|\lambda_2/\lambda_1|^k$, where $K = \sum_{i \ge 2}|c_i|\,\|\mathbf v_i\|$, because $|\lambda_i| \le |\lambda_2|$ for $i \ge 2$. In the picture the ratio of the parts, $|c_2/c_1|\,(1/2)^k$, halves at each step, and the points on the circle creep towards the line of $\mathbf v_1$. Let’s turn this into a distance between directions. For any non-zero $\mathbf y$, $\mathbf z$ we have $\bigl\|\frac{\mathbf y}{\|\mathbf y\|} - \frac{\mathbf z}{\|\mathbf z\|}\bigr\| \le \frac{2\|\mathbf y - \mathbf z\|}{\|\mathbf z\|}$: the difference equals $\frac{\mathbf y - \mathbf z}{\|\mathbf z\|} + \mathbf y\bigl(\frac{1}{\|\mathbf y\|} - \frac{1}{\|\mathbf z\|}\bigr)$, and the second term has length $\frac{|\|\mathbf z\| - \|\mathbf y\||}{\|\mathbf z\|} \le \frac{\|\mathbf y - \mathbf z\|}{\|\mathbf z\|}$. Take $\mathbf y = c_1\mathbf v_1 + \mathbf r_k$ and $\mathbf z = c_1\mathbf v_1$: the vector $\mathbf x_k$ differs from $\mathbf y/\|\mathbf y\|$ at most by its sign (the factor $\lambda_1^k$ may be negative), so the distance from $\mathbf x_k$ to $\pm\mathbf v_1/\|\mathbf v_1\|$ is at most $\frac{2\|\mathbf r_k\|}{|c_1|\,\|\mathbf v_1\|} \le \frac{2K}{|c_1|\,\|\mathbf v_1\|}\,|\lambda_2/\lambda_1|^k$. This is where $c_1 \ne 0$ is needed: put $\mathbf x$ on the line of $\mathbf v_2$, and the first part disappears, leaving nothing to converge to. So $\mathbf x_k$ comes arbitrarily close to $\pm\mathbf u$, where $\mathbf u = \mathbf v_1/\|\mathbf v_1\|$. For a unit eigenvector $\langle \mathbf u, A\mathbf u\rangle = \lambda_1\langle \mathbf u, \mathbf u\rangle = \lambda_1$ (and for $-\mathbf u$ too), and the dot product $\langle \mathbf x, A\mathbf x\rangle$ depends continuously on $\mathbf x$, so $\langle \mathbf x_k, A\mathbf x_k\rangle \to \lambda_1$.

Where the method breaks down can be seen from the same two cases as before. If the two leading eigenvalues are equal in absolute value, like $1$ and $-1$ for a reflection, no direction wins. If they’re complex, the fan spins for ever. And if $|\lambda_2/\lambda_1|$ is close to one, the method works in principle, but very slowly.

Scooters spread across the city

A company rents out $600$ e-scooters. For the company the city splits into two zones: the centre and the outskirts. The statistics for the season are these: of the scooters that start the day in the centre, $90\,\%$ end it in the centre and $10\,\%$ go off to the outskirts. Of those parked in the outskirts, by evening half end up in the centre and half stay put. Where will the scooters be in a week if, after an all-night festival, all $600$ have gathered in the outskirts?

Write the shares of scooters in the centre and the outskirts as a vector $\mathbf p = (p_1, p_2)$. One day is multiplication by a matrix:

$$\mathbf p_{k+1} = P\,\mathbf p_k, \qquad P = \begin{pmatrix} 0.9 & 0.5 \\ 0.1 & 0.5 \end{pmatrix}.$$

Column $j$ records where the scooters from zone $j$ scatter in a day; each column adds up to $1$, because scooters don’t vanish. The expected number of scooters in the centre, day after day: $0$, $300$, $420$, $468$, $487.2$, $494.9$, and after that it closes in on $500$. The shortfall from $500$ is $500$, $200$, $80$, $32$, $12.8$, shrinking by a factor of exactly $2.5$ each day.

This is the power method in its pure form. The first eigenvalue of $P$ is $1$ (why, see below), the trace is $1.4$, so the second eigenvalue is $1.4 - 1 = 0.4$: exactly the factor the shortfall is multiplied by each day. And the limit is the eigenvector with $\lambda = 1$. The equation $P\boldsymbol\pi = \boldsymbol\pi$ comes down to a balance of flows, $0.1\,\pi_1 = 0.5\,\pi_2$: as many scooters leave the centre as arrive. Together with the condition $\pi_1 + \pi_2 = 1$ we get $\boldsymbol\pi = (5/6, 1/6)$: $500$ scooters in the centre and $100$ in the outskirts, whatever distribution we start from.

As with the scooters, convergence to equilibrium is guaranteed if all the transition probabilities are positive. To prove it, we’ll measure the distance between distributions by the sum of the absolute values of the differences of the coordinates: $\|\mathbf y\|_1 = |y_1| + \dots + |y_n|$. The key is a single lemma.

The idea: the matrix $P$ maps the set of all distributions inside itself, and each application shrinks everything by a factor of at least $1/q$. The images, nested inside one another, shrink to a single point, and that point is $\boldsymbol\pi$. The picture is for three states.

A distribution on three states is a triple $(p_1, p_2, p_3)$ with $p_i \ge 0$ adding up to $1$. Draw it as the point $p_1\mathbf e_1 + p_2\mathbf e_2 + p_3\mathbf e_3$ of the triangle with vertices $\mathbf e_1, \mathbf e_2, \mathbf e_3$: it’s the centre of mass of weights $p_1, p_2, p_3$ placed at the vertices. The vertex $\mathbf e_j$ means the chain is certainly in state $j$. $P\mathbf e_j$ is the $j$th column of $P$, a distribution too; its coordinates are at least $\varepsilon$, so the point lies strictly inside the triangle. The matrix is linear: $P(p_1\mathbf e_1 + p_2\mathbf e_2 + p_3\mathbf e_3) = p_1P\mathbf e_1 + p_2P\mathbf e_2 + p_3P\mathbf e_3$, so centres of mass go to centres of mass with the same weights. Hence the whole triangle goes to the triangle with vertices at the columns of $P$. Apply $P$ again: $P^2(\Delta) = P(P(\Delta)) \subset P(\Delta)$, and so on: the triangles are nested inside one another. Moreover they shrink: the difference of two distributions has coordinates adding up to $0$, and by the lemma the distance between the images is at most $q$ times the distance between the originals. At the start any two distributions are at most $2$ apart (the sum of the absolute values of the difference is at most $1 + 1$), so any two points of $P^k(\Delta)$ are at most $2q^k$ apart. Take any distribution $\mathbf p$ and the sequence $\mathbf p_k = P^k\mathbf p$. Neighbouring terms are close: $\|\mathbf p_{k+1} - \mathbf p_k\|_1 = \|P^k(P\mathbf p - \mathbf p)\|_1 \le 2q^k$, and the series $2 + 2q + 2q^2 + \dots$ converges. So each coordinate of $\mathbf p_k$ is a partial sum of an absolutely convergent series and has a limit (chapter 30). The limit $\boldsymbol\pi$ is again a distribution (non-negativity and the sum $1$ survive in the limit), and passing to the limit in the equation $\mathbf p_{k+1} = P\mathbf p_k$, we get $\boldsymbol\pi = P\boldsymbol\pi$. This is the point the triangles shrink to. Uniqueness and positivity. If $P\boldsymbol\pi' = \boldsymbol\pi'$ for another distribution, then $\|\boldsymbol\pi - \boldsymbol\pi'\|_1 = \|P(\boldsymbol\pi - \boldsymbol\pi')\|_1 \le q\,\|\boldsymbol\pi - \boldsymbol\pi'\|_1$, and with $q < 1$ this is possible only if $\boldsymbol\pi' = \boldsymbol\pi$. Each coordinate $\pi_i = \sum_j p_{ij}\pi_j \ge \varepsilon\sum_j \pi_j = \varepsilon$: the point $\boldsymbol\pi$ lies in the inner triangle where all the coordinates are at least $\varepsilon$. Finally, for any starting $\mathbf p$: $\|P^k\mathbf p - \boldsymbol\pi\|_1 = \|P^k(\mathbf p - \boldsymbol\pi)\|_1 \le q^k\,\|\mathbf p - \boldsymbol\pi\|_1 \le 2q^k$: step by step the point $\mathbf p$ moves into $\boldsymbol\pi$. For $n$ states the argument is the same, only with the polytope of distributions in $n$-dimensional space instead of a triangle.

For the scooters the smallest entry of $P$ is $0.1$, there are two states, and the theorem guarantees a factor of $q = 1 - 2 \cdot 0.1 = 0.8$ per day. In fact the shortfall shrank faster, by a factor of $0.4$: that’s the second eigenvalue. The lemma gives an estimate that works for every positive matrix at once, while the exact rate is set by $|\lambda_2|$.

In books on probability and in chapter 50 the transition matrix is usually written transposed, so that the rows add up to one, and the distribution is written as a row: $\boldsymbol\pi P = \boldsymbol\pi$. It’s the same eigenvector problem, only the multiplication is on the left. Andrei Markov began studying such chains in 1906, and how he tested them on the text of Eugene Onegin is told in that chapter.

Another company works in two districts, $A$ and $B$. In a day $30\,\%$ of the scooters in $A$ go to $B$, and $20\,\%$ of the scooters in $B$ go to $A$. What share of all the scooters will be in district $A$ after many days? Give the answer as a fraction or a decimal.

The transition matrix is $P = \left(\begin{smallmatrix} 0.7 & 0.2 \\ 0.3 & 0.8 \end{smallmatrix}\right)$. The balance of flows: $0.3\,\pi_A = 0.2\,\pi_B$, that is, $\pi_B = 1.5\,\pi_A$. Together with $\pi_A + \pi_B = 1$ we get $\pi_A = 0.4$. The second eigenvalue is $\operatorname{tr}P - 1 = 0.5$, so the deviation from equilibrium halves every day.

You can get your hand in with eigenvalues, eigenvectors and Markov chains in the trainer. And the solver will take a matrix of your own apart step by step.

How Google ranked the pages

In 1998 two Stanford graduate students, Sergey Brin and Larry Page, published the paper The Anatomy of a Large-Scale Hypertextual Web Search Engine, and in the autumn of the same year they founded Google. One of the paper’s main ideas was PageRank, a way of judging the importance of a page by the links to it. A page is important if important pages link to it. The definition looks like a vicious circle, but in fact it’s an equation for an eigenvector.

Let page $j$ have $L_j$ outgoing links. It shares its importance $r_j$ equally among the pages it links to, and each gets $r_j/L_j$. A page’s importance is the sum of everything it receives:

$$r_i = \sum_{j \to i} \frac{r_j}{L_j}, \qquad \text{that is,} \qquad \mathbf r = M\mathbf r,$$

where column $j$ of the matrix $M$ has the numbers $1/L_j$ opposite the pages that $j$ links to, and zeros elsewhere. The columns of $M$ add up to $1$, so this is the transition matrix of a Markov chain, and $\mathbf r$ is its stationary distribution. It has a vivid meaning: it’s the share of time spent on each page by a random surfer who clicks random links for ever.

With the real internet it’s not so simple. There are pages with no outgoing links, and the random surfer gets stuck on them. There are groups of pages that link only to one another: like a trap, they gradually take all the importance for themselves. And if the web falls into unconnected pieces, there are many stationary distributions. Brin and Page fixed all this with a single trick: with probability $d$ the surfer clicks a random link, and with probability $1 - d$ gets bored and opens a random page out of all $N$. The paper says that $d = 0.85$ is usually taken. From a page with no links we’ll assume the surfer always jumps to a random one, which is what’s usually done.

The probability of following a link, usually $0.85$. The link matrix: column $j$ shares one equally among the pages that $j$ links to. The share each of the $N$ pages gets from “teleportation”; $\mathbf 1$ is the vector of ones. Since $r_1 + \dots + r_N = 1$, the equation can be written as $\mathbf r = G\mathbf r$ with the matrix $G = dM + \frac{1 - d}{N}J$, where $J$ is the matrix of all ones. Example: for the three pages above with $d = 0.85$ we get $\mathbf r = \frac{1}{1769}(686, 380, 703) \approx (0.388, 0.215, 0.397)$. Now $C$ is slightly ahead of $A$: part of the importance that $C$ passed to $A$ spreads over all the pages.

Why does this equation have a solution, and a unique one? We already have the answer in hand. The matrix $G$ is a transition matrix: the columns of $M$ add up to $1$ (from a page with no links the surfer goes to any of the $N$ with probability $1/N$), so do the columns of $\frac1N J$, and $d + (1 - d) = 1$. Moreover, all the entries of $G$ are at least $\varepsilon = \frac{1 - d}{N}$. By the theorem on the stationary distribution of a positive chain with $n = N$ we get $q = 1 - N\varepsilon = d$: the ranking exists and is unique, each page gets a weight of at least $\frac{1 - d}{N}$, and the power method converges, with error after $k$ steps at most $2d^k$. At $d = 0.85$, after $50$ iterations that’s $2 \cdot 0.85^{50} \approx 6 \cdot 10^{-4}$. And $G$ can be multiplied by a vector in a single pass through the list of links, even when there are billions of pages.

Our theorem on positive chains is a special case of a general result about matrices with positive entries.

The perfect quarry: symmetric matrices

Back to the curve $5x^2 + 4xy + 2y^2 = 6$. Its left side is a quadratic form, a homogeneous polynomial of degree two in the coordinates. Any such form can be written with a symmetric matrix:

$$5x^2 + 4xy + 2y^2 = \begin{pmatrix} x & y \end{pmatrix}\begin{pmatrix} 5 & 2 \\ 2 & 2 \end{pmatrix}\begin{pmatrix} x \\ y \end{pmatrix} = \langle \mathbf x, S\mathbf x\rangle.$$

The coefficient of $xy$ is split equally between the two off-diagonal cells: $4xy = 2xy + 2yx$. A matrix that doesn’t change under transposition ($S^{\mathsf T} = S$) is called symmetric. With such matrices the hunt always succeeds.

Let’s start with the plane, where everything can be seen in a picture. For a symmetric matrix $S$ it’s convenient to look at the function $Q(\mathbf x) = \langle \mathbf x, S\mathbf x\rangle$ on the unit circle: as we’ll see, its largest and smallest values are the eigenvalues.

The idea: find the point on the unit circle where the form $Q$ is largest, and notice that there its gradient points along the radius. And that’s exactly what $S\mathbf u \parallel \mathbf u$ means.

Write the unit vectors as $\mathbf u(\theta) = (\cos\theta, \sin\theta)$ and consider the function $f(\theta) = Q(\mathbf u(\theta)) = a\cos^2\theta + 2b\sin\theta\cos\theta + c\sin^2\theta$. It’s continuous on the interval $[0, 2\pi]$ and so takes a largest value there (Weierstrass’s theorem, chapter 53). Call the point of maximum $\mathbf u_1 = \mathbf u(\theta_1)$. In the picture the ring is darker where $Q$ is larger. Geometrically, at this point the level curve $Q = f(\theta_1)$ touches the circle: the whole circle lies in the region $Q \le f(\theta_1)$, and at the point $\mathbf u_1$ the boundary of that region touches the circle. Let’s check this by calculation. The function $f$ is periodic, so the point of maximum can be taken to be an interior point of the interval (shift the interval of length $2\pi$ so that $\theta_1$ is in its middle), and at an interior point of maximum the derivative is zero (chapter 27). By the product rule $f' = \langle \mathbf u', S\mathbf u\rangle + \langle \mathbf u, S\mathbf u'\rangle$. For a symmetric $S$ the second term equals the first: $\langle \mathbf u, S\mathbf u'\rangle = \mathbf u^{\mathsf T}S\mathbf u' = (S\mathbf u)^{\mathsf T}\mathbf u' = \langle S\mathbf u, \mathbf u'\rangle$, because $S^{\mathsf T} = S$. So $f'(\theta) = 2\langle \mathbf u'(\theta), S\mathbf u(\theta)\rangle$, and at the point of maximum $S\mathbf u_1$ is perpendicular to the vector $\mathbf u'(\theta_1) = (-\sin\theta_1, \cos\theta_1)$. In the plane only multiples of $\mathbf u_1$ are perpendicular to this vector, so $S\mathbf u_1 = \lambda_1\mathbf u_1$ for some real $\lambda_1$. That’s what the picture shows: the gradient $2S\mathbf u_1$ points along the radius. Taking the dot product of the equation with $\mathbf u_1$, we find $\lambda_1 = \langle \mathbf u_1, S\mathbf u_1\rangle = \max Q$. Take the unit vector $\mathbf u_2 = (-\sin\theta_1, \cos\theta_1)$, perpendicular to $\mathbf u_1$. Again by symmetry $\langle S\mathbf u_2, \mathbf u_1\rangle = \langle \mathbf u_2, S\mathbf u_1\rangle = \lambda_1\langle \mathbf u_2, \mathbf u_1\rangle = 0$. So $S\mathbf u_2$ is perpendicular to $\mathbf u_1$, that is, a multiple of $\mathbf u_2$: $S\mathbf u_2 = \lambda_2\mathbf u_2$ with real $\lambda_2 = \langle \mathbf u_2, S\mathbf u_2\rangle$. Any unit vector can be written as $\mathbf x = \cos\alpha\,\mathbf u_1 + \sin\alpha\,\mathbf u_2$, and then $Q(\mathbf x) = \lambda_1\cos^2\alpha + \lambda_2\sin^2\alpha$ (the mixed terms have gone, because $\mathbf u_1 \perp \mathbf u_2$). This number is between $\lambda_2$ and $\lambda_1$, so $\lambda_2 = \min Q$: the inner level curve $Q = \lambda_2$ touches the circle at the points $\pm\mathbf u_2$. We’ve got two perpendicular unit eigenvectors with real eigenvalues.

In any dimension the argument is the same, only the circle becomes a sphere, and the move to the perpendicular direction becomes an induction.

The idea: an eigenvector gives the maximum of the form on the sphere, and the matrix $S$ doesn’t let the subspace perpendicular to it out, so the same can be repeated there.

Induction on $n$; for $n = 1$ there’s nothing to prove. Let $n \ge 2$. The function $Q(\mathbf x) = \langle \mathbf x, S\mathbf x\rangle$ is continuous on the unit sphere $\{\mathbf x : \|\mathbf x\| = 1\}$, a closed bounded set, so it attains its largest value on it at some point $\mathbf u_1$ (Weierstrass’s theorem for functions of several variables, chapter 53). $\mathbf u_1$ is an eigenvector. Take any unit $\mathbf w \perp \mathbf u_1$. The circle $\cos\theta\,\mathbf u_1 + \sin\theta\,\mathbf w$ lies on the sphere, and on it $Q$ is largest at $\theta = 0$. The same calculation as in the plane gives $\langle \mathbf w, S\mathbf u_1\rangle = 0$. The vector $S\mathbf u_1$ is perpendicular to all vectors perpendicular to $\mathbf u_1$, so it’s a multiple of $\mathbf u_1$: split $S\mathbf u_1 = \alpha\mathbf u_1 + \mathbf w_0$ with $\mathbf w_0 \perp \mathbf u_1$; then $0 = \langle \mathbf w_0, S\mathbf u_1\rangle = \|\mathbf w_0\|^2$, that is, $\mathbf w_0 = \mathbf 0$. So $S\mathbf u_1 = \lambda_1\mathbf u_1$, and $\lambda_1 = \langle \mathbf u_1, S\mathbf u_1\rangle$ is a real number. The subspace $W = \{\mathbf w : \langle \mathbf w, \mathbf u_1\rangle = 0\}$ has dimension $n - 1$, and $S$ takes it to itself: if $\mathbf w \perp \mathbf u_1$, then $\langle S\mathbf w, \mathbf u_1\rangle = \langle \mathbf w, S\mathbf u_1\rangle = \lambda_1\langle \mathbf w, \mathbf u_1\rangle = 0$. Choose an orthonormal basis $\mathbf w_1, \dots, \mathbf w_{n-1}$ in $W$ (it’s built by the Gram–Schmidt process from chapter 39). In this basis the transformation $S$ on $W$ is written with the matrix with entries $\langle \mathbf w_i, S\mathbf w_j\rangle$, and it’s symmetric again: $\langle \mathbf w_i, S\mathbf w_j\rangle = \langle S\mathbf w_i, \mathbf w_j\rangle$. By the inductive hypothesis $W$ has an orthonormal basis $\mathbf u_2, \dots, \mathbf u_n$ of eigenvectors with real eigenvalues. The vectors $\mathbf u_1, \dots, \mathbf u_n$ are unit vectors and pairwise perpendicular, that is, they form an orthonormal basis of $\mathbb R^n$ of eigenvectors. By the diagonalisation criterion $S = UDU^{-1}$, and for a matrix $U$ with orthonormal columns $U^{\mathsf T}U = I$ (the entries of the product are the dot products of the columns), so $U^{-1} = U^{\mathsf T}$. The characteristic polynomial of $S$ is $(\lambda_1 - \lambda)\cdots(\lambda_n - \lambda)$ with real $\lambda_i$, so $S$ has no other eigenvalues, complex or otherwise.

Matrices with orthonormal columns are called orthogonal; more about them in chapter 39. For $2 \times 2$ the reality of the eigenvalues can be seen without any maxima: the discriminant of the characteristic polynomial of $S = \left(\begin{smallmatrix} a & b \\ b & c \end{smallmatrix}\right)$ is $(a + c)^2 - 4(ac - b^2) = (a - c)^2 + 4b^2 \ge 0$.

In our example the trace is $7$, the determinant $6$, the eigenvalues $6$ and $1$, and the eigenvectors $(2, 1)$ and $(1, -2)$ are perpendicular, as promised: $2 \cdot 1 + 1 \cdot (-2) = 0$. Draw new axes $X$ and $Y$ along the unit vectors $\frac{1}{\sqrt5}(2, 1)$ and $\frac{1}{\sqrt5}(1, -2)$. In these coordinates the matrix is diagonal, the product $XY$ disappears, and the form becomes a sum of squares.

The first eigenvalue; $X$ is the coordinate along its unit eigenvector. The second eigenvalue; $Y$ is the coordinate along the perpendicular direction. There’s no mixed term, because in eigencoordinates $S$ is diagonal. Example: $5x^2 + 4xy + 2y^2 = 6$ turns into $6X^2 + Y^2 = 6$, that is, $X^2 + \frac{Y^2}{6} = 1$. This is an ellipse with semi-axes $1$ along $(2, 1)$ and $\sqrt6 \approx 2.449$ along $(1, -2)$. Its long axis is tilted to the $x$-axis at the angle $\arctan 2 \approx 63.4°$, and its area is $\pi\sqrt6 \approx 7.70$.

The lines along the eigenvectors of a symmetric matrix are called the principal axes of the quadratic form. The signs of the eigenvalues tell at once what kind of curve we have.

For our curve $\det S = 5 \cdot 2 - 2^2 = 6 > 0$ and the trace is $7 > 0$, so both eigenvalues are positive: an ellipse, just as we found.

The picture from the proof also works as a tool for finding the axes. The normal to the curve $\langle \mathbf x, S\mathbf x\rangle = 6$ at a point $\mathbf x$ points along the gradient, and the gradient of a quadratic form is $2S\mathbf x$ (chapter 32). The ends of the axes are the points where the radius vector is perpendicular to the curve, that is, points along the normal: $S\mathbf x \parallel \mathbf x$, the eigenvector equation again. For an ellipse the length of the radius is largest and smallest at these points.

Principal axes turn up wherever there’s a symmetric matrix. For a rigid body they’re the axes of inertia: any body, be it a book or a space station, has three mutually perpendicular axes about which it can spin without wobbling. They’re the eigenvectors of a symmetric matrix called the inertia tensor. Rotation about the axes with the largest and smallest moments of inertia is stable, but about the middle one it isn’t. Toss a book held shut with a rubber band, spinning it about each of the three axes in turn: in one of the three cases it starts to tumble. In 1985 the cosmonaut Vladimir Dzhanibekov noticed aboard the Salyut 7 station that a wing nut spinning in weightlessness flips over at regular intervals. Since then this has been called the Dzhanibekov effect, or the tennis racket theorem.

Where next

Take a cloud of points, for example the heights and weights of a thousand people, and try to draw through it the line that describes it best. The points don’t lie on a line: the system “all the points on one line” is inconsistent, and no solution fits exactly. Which approximation counts as the best? The answer will again lead to a symmetric matrix and its eigenvectors. And a rectangular matrix has no eigenvectors at all: it takes vectors from one space to another, and $A\mathbf v$ can’t even be compared with $\mathbf v$. What takes their place, and how it’s used to compress photographs, is the subject of chapter 39.