Multiply a vector by a matrix and two things happen at once: the arrow swings round to a new angle, and its length changes. Do that to a fan of arrows and the picture is a mess. Every one moves differently, and nothing about the grid of numbers tells you what the matrix is for.
A few directions escape the mess. Feed one of those in and the arrow comes back pointing exactly where it started, only longer or shorter — turned by nothing at all. Those directions are the eigenvectors of the matrix, and the amount each one stretches by is its eigenvalue. Find them and a grid of numbers collapses into a short list of directions and scale factors.
Most matrices have a few directions they only stretch
Everything below runs on one small matrix, \(A = \begin{bmatrix} 4 & 1 \\ 2 & 3 \end{bmatrix}\). Those four numbers are made up — hand-picked for this post, not measured off anything. Read \(A\) as one step of a two-group population model: two classes of something, and next year’s size of each is a mixture of this year’s two, with the four entries as the mixing rates. Whole numbers were chosen because they make the eigenvalues come out at exactly 5 and 2, which keeps the arithmetic out of the way of the geometry. A matrix fitted to real counts would run to hundreds of dimensions with decimal eigenvalues, and you could not draw any of it. Drawing it is what the rest of this post does.
A repeated map is also the honest home for this technique rather than a pretext for it: apply the same matrix year after year and its eigenvalues are what decide where the population ends up, which is where this post finishes. One thing to notice before that: \(A\) is not symmetric, and that matters shortly.
numpy.linalg.eig hands you the directions and the stretches at once:
Code
import numpy as npnp.set_printoptions(precision=3, suppress=True)A = np.array([[4.0, 1.0], [2.0, 3.0]])vals, vecs = np.linalg.eig(A)# numpy 2.x returns complex dtype from eig even when every eigenvalue is real,# as they are here. It also returns unit-length eigenvectors; since these are# defined only up to scale, rescale each to start with a 1 so the whole numbers# below match the text. (Fine here -- a zero first entry would need another# pivot, and a rotation matrix has no real eigenvector to rescale at all.)vals, vecs = vals.real, vecs.realvecs = vecs / vecs[0]for lam, v inzip(vals, vecs.T):print(f"lambda = {lam:5.2f} eigenvector = {v} A @ v = {A @ v}")
lambda = 5.00 eigenvector = [1. 1.] A @ v = [5. 5.]
lambda = 2.00 eigenvector = [ 1. -2.] A @ v = [ 2. -4.]
Both do what they promise: A @ v is v with every number multiplied by the same amount. Two directions, two numbers — and that pair is enough to rebuild \(A\) itself.
\(A = PDP^{-1}\) is that idea written as a formula
Put the eigenvectors in the columns of \(P\), their eigenvalues down the diagonal of \(D\), and you can rebuild the matrix:
Read it right to left, as three moves. \(P^{-1}\) redescribes an incoming vector in eigenvector terms — how much of it points along \(v_1\), how much along \(v_2\). \(D\) stretches each of those parts by its own eigenvalue. \(P\) converts the result back to ordinary coordinates. Nothing here is new; it is the sentence that opened this post, written so a computer can use it.
The eigenvalues themselves come from the characteristic equation \(\det(A - \lambda I) = 0\) — the condition for \(A\) to squash some direction to zero once you subtract \(\lambda\) from its diagonal. For this \(A\) it is the quadratic \(\lambda^2 - 7\lambda + 10 = 0\), whose roots are the 5 and 2 printed above.
Multiplying the three matrices back together gives \(A\) again. Drawing the map on a circle of unit vectors shows why the eigenvectors are the special ones:
Code
import matplotlib.pyplot as pltACCENT ="#4A3AA7"GREY ="#999999"P = vecsD = np.diag(vals)print(f"max |A - P D P^-1| = {np.abs(A - P @ D @ np.linalg.inv(P)).max():.2e}")theta = np.linspace(0, 2* np.pi, 400)circle = np.vstack([np.cos(theta), np.sin(theta)])image = A @ circlefig, ax = plt.subplots(figsize=(6.2, 6.2))ax.axhline(0, color="0.85", lw=1, zorder=0)ax.axvline(0, color="0.85", lw=1, zorder=0)ax.plot(*circle, color=GREY, lw=1.3, label="unit circle")ax.plot(*image, color="#444444", lw=1.6, label="image under $A$")for k, (lam, v) inenumerate(zip(vals, vecs.T)): u = v / np.linalg.norm(v) ax.annotate("", xy=lam * u, xytext=(0, 0), arrowprops=dict(arrowstyle="-|>", color=ACCENT, lw=2.2, shrinkA=0, shrinkB=0), zorder=3, )# Mark where the eigenvector meets the unit circle: the arrow is that same# point scaled by lambda, and its tip lands back on the ellipse. ax.plot(*u, marker="o", ms=6, color="white", mec=ACCENT, mew=1.8, zorder=5) ax.plot([], [], color=ACCENT, lw=2.2, label=f"eigenvector, $\\lambda = {lam:.0f}$") ax.annotate(f"$\\lambda_{k +1} = {lam:.0f}$", xy=lam * u, xytext=(16, 12* np.sign(u[1])), textcoords="offset points", color=ACCENT, fontsize=11, va="center", )ax.set_aspect("equal")ax.set_xlim(-6, 6)ax.set_ylim(-6, 6)ax.grid(False)for side in ("top", "right"): ax.spines[side].set_visible(False)ax.legend(frameon=False, loc="upper left", fontsize=10)ax.set_title("Dots: unit eigenvectors on the circle. Arrows: their images under $A$.", fontsize=10, color="0.3")plt.show()
max |A - P D P^-1| = 8.88e-16
The reconstruction error prints at machine precision, so the formula holds. In the picture, the circle becomes an ellipse and the purple arrows land near its long and short axes without quite reaching them. That near miss is real rather than a bug. The axes of the ellipse mark the directions \(A\) stretches by the most and the least — its singular vectors — and those coincide with the eigenvectors only when the matrix is symmetric. This one is not, which is the detail flagged earlier, and it comes back at the end.
Diagonal matrices turn repeated work into simple arithmetic
Rebuilding \(A\) from its own parts is a party trick until you apply it twice. Square \(PDP^{-1}\) and the inner \(P^{-1}P\) cancels, leaving \(A^2 = PD^2P^{-1}\); the same cancellation runs all the way up, so \(A^k = PD^kP^{-1}\). Raising a matrix to the hundredth power becomes raising two numbers to that power, and the larger of the two soon dwarfs the smaller — which is how one eigenvalue ends up dictating the long run of a population model. That is why eigenvalues turn up wherever the same map runs again and again.
Which is most places, under different names:
Concept
What it is
Where it shows up
Eigenvector
a direction \(v\) with \(Av = \lambda v\)
PCA’s principal axes
Eigenvalue
the scale factor \(\lambda\) along \(v\)
variance explained
Characteristic equation
\(\det(A - \lambda I) = 0\)
finding the spectrum
Decomposition
\(A = PDP^{-1}\)
Markov chain steady states
When it exists
\(n\) independent eigenvectors (diagonalizable)
vibration normal modes
Two kinds of matrix have no such viewpoint to offer. A rotation turns every direction without exception, so it has no real eigenvector to find. A defective matrix has eigenvectors but too few independent ones to fill the columns of \(P\), so there is nothing to invert. Neither failure is exotic, and both are why the decomposition is stated with a condition attached.
The asymmetry of this \(A\) costs something milder. Its eigenvectors are not at right angles to each other, so \(P^{-1} \neq P^{\top}\) and the change of viewpoint has to be inverted rather than transposed — an operation that goes numerically wobbly when the directions sit close together. That is why Matrix Factorizations as Optimization Problems reaches instead for the singular value decomposition, the SVD — the factorization built on those singular vectors rather than the eigenvectors, which buys a right angle for any matrix at all.
Inside those limits, the opening claim stands. The directions a matrix refuses to turn are the viewpoint from which it stops being a grid of numbers at all, and becomes a handful of arrows with a scale factor each.
