8. Eigenvalues and Eigenvectors

Most vectors change direction when a matrix acts on them, but a few special directions are only stretched, shrunk, or flipped; along those directions the matrix acts like a single number. This chapter finds those directions (eigenvectors) and numbers (eigenvalues), uses them to rewrite a matrix in its simplest possible form (diagonal, triangular, or Jordan), and applies the result to repeated processes, Markov chains, differential equations, population models, and PageRank. Eigenvalues are the main tool for answering "what happens in the long run?" about any linear process.

<a id="81-eigenvalue-fundamentals"></a>

8.1 Eigenvalue Fundamentals

This subchapter defines eigenvalues and eigenvectors, shows how to find them through the characteristic polynomial, and explains what happens when eigenvalues repeat or are complex.

Definition of Eigenvalues and Eigenvectors

Why it matters. A $2 \times 2$ matrix already mixes its inputs in a way that is hard to read: $A = \begin{bmatrix} 2 & 1 \ 1 & 2 \end{bmatrix}$ sends $(1, 0)$ to $(2, 1)$, turning it as well as stretching it. But it sends $(1, 1)$ to $(3, 3)$ and $(1, -1)$ to $(1, -1)$. Along these two directions $A$ is just "multiply by 3" and "multiply by 1". If every input can be split into pieces along such directions, then applying $A$, or applying it a hundred times, reduces to multiplying numbers. This is why eigenvalues govern vibration modes of bridges, energy levels in quantum mechanics, principal components in statistics, the stability of dynamical systems, and the ranking of web pages.

Analogy. Pull a rubber sheet so that it stretches along one diagonal and compresses along the other. A dot drawn on the sheet usually moves sideways as well as outward. Dots on the two diagonals only slide straight out from or in toward the center. Those two diagonals are the eigen-directions of the deformation, and the stretch factors along them are the eigenvalues. For a spinning globe, the axis is the only direction that does not move: it is an eigen-direction with eigenvalue 1.

Definition. Let $A$ be an $n \times n$ matrix over a field $F$ (or, more generally, let $T : V \to V$ be a linear operator). A scalar $\lambda \in F$ is an eigenvalue of $A$ if there is a nonzero vector $\mathbf{v}$ with

$$A\mathbf{v} = \lambda\mathbf{v}.$$

Such a $\mathbf{v}$ is an eigenvector of $A$ for the eigenvalue $\lambda$, and $(\lambda, \mathbf{v})$ is an eigenpair. The word eigen is German for "own" or "characteristic": these are the matrix's own directions.

Three details of the definition matter:

• $\mathbf{v} \ne \mathbf{0}$ is required. $A\mathbf{0} = \lambda\mathbf{0}$ holds for every $\lambda$, so allowing $\mathbf{0}$ would make every scalar an eigenvalue.
• $\lambda = 0$ is allowed. $A\mathbf{v} = 0\mathbf{v} = \mathbf{0}$ with $\mathbf{v} \ne \mathbf{0}$ means $\mathbf{v}$ is a nonzero vector in $\operatorname{Nul} A$. So $0$ is an eigenvalue of $A$ exactly when $A$ is singular, which is item (n) of the invertible matrix theorem (Section 4.3).
• Eigenvectors are directions, not single vectors. If $A\mathbf{v} = \lambda\mathbf{v}$, then $A(c\mathbf{v}) = c,A\mathbf{v} = \lambda(c\mathbf{v})$ for every scalar $c$. Any nonzero multiple of an eigenvector is again an eigenvector for the same eigenvalue. Software usually returns eigenvectors scaled to length 1, with an arbitrary sign.

What the eigenvalue says geometrically. For a real eigenvalue, the line through $\mathbf{v}$ is mapped onto itself, and $\lambda$ says how.

Eigenvalue	Effect on the eigen-direction
$\lambda > 1$	stretched away from the origin
$0 < \lambda < 1$	shrunk toward the origin
$\lambda = 1$	left fixed, point by point
$\lambda = 0$	collapsed to the origin (the direction is in the kernel)
$-1 < \lambda < 0$	flipped through the origin and shrunk
$\lambda < -1$	flipped and stretched
non-real $\lambda$	no real fixed direction; the map rotates (see Complex Eigenvalues below)

   A = [ 2  1 ]       most vectors turn,        eigenvectors keep their line
       [ 1  2 ]

        y                                     y
        |      A e1 = (2, 1)                  |        A(1, 1) = (3, 3)
        |     ^                               |       ^
        |    /                                |      /
        |   /                                 |     /
        |  /                                  |    /  (1, 1): stretched by 3
        | /                                   |   ^
        |/                                    |  /
        O---------> e1 = (1, 0)   x           O-------------------- x
                                              |\
                                              | \
                                              |  v (1, -1) = A(1, -1): stretched by 1

A catalogue of examples. Many eigenpairs can be found by thinking geometrically, before any computation.

Map	Eigenvalues and eigenvectors
projection of $\mathbb{R}^2$ onto a line $L$	$1$ on $L$; $0$ on the direction that is collapsed
reflection of $\mathbb{R}^2$ across a line $L$	$1$ on $L$; $-1$ on the perpendicular line
scaling $\mathbf{x} \mapsto c\mathbf{x}$	$c$, and every nonzero vector is an eigenvector
rotation of $\mathbb{R}^2$ by $90°$	no real eigenvalues: every nonzero vector changes direction
rotation of $\mathbb{R}^3$ about an axis	$1$ on the axis
shear $\begin{bmatrix} 1 & 1 \ 0 & 1 \end{bmatrix}$	$1$, with eigenvectors only along the $x$-axis
$D = \tfrac{d}{dx}$ on smooth functions	every $\lambda$: $D(e^{\lambda x}) = \lambda e^{\lambda x}$

The last row is the reason exponentials solve linear differential equations (Section 8.4). Eigenvectors of operators on function spaces are called eigenfunctions.

Checking a claimed eigenpair is cheap. Finding eigenvectors takes work (later topics), but verifying one takes a single multiplication: compute $A\mathbf{v}$ and see whether it is a multiple of $\mathbf{v}$.

import numpy as np

A = np.array([[2.0, 1.0],
              [1.0, 2.0]])

for v in [np.array([1.0, 0.0]), np.array([1.0, 1.0]), np.array([1.0, -1.0])]:
    Av = A @ v
    # Av is parallel to v exactly when the 2x2 determinant [v Av] is zero
    parallel = np.isclose(v[0] * Av[1] - v[1] * Av[0], 0.0)
    print(v, "->", Av, "eigenvector" if parallel else "turned")
# [1. 0.] -> [2. 1.] turned
# [1. 1.] -> [3. 3.] eigenvector
# [ 1. -1.] -> [ 1. -1.] eigenvector

# numpy returns the eigenvalues and unit-length eigenvectors (as columns).
# The arrays may have a complex dtype; .real drops the zero imaginary parts here.
vals, vecs = np.linalg.eig(A)
print(vals.real)                          # [3. 1.]
print(np.round(vecs.real, 4))             # [[ 0.7071 -0.7071]
                                          #  [ 0.7071  0.7071]]
for lam, v in zip(vals, vecs.T):          # iterate over columns
    print(np.allclose(A @ v, lam * v))    # True, True

# Eigenvalue 0 means singular
S = np.array([[1.0, 2.0], [2.0, 4.0]])
print(np.round(np.linalg.eigvals(S).real, 10), np.linalg.matrix_rank(S))   # [0. 5.] 1

Key takeaways

• $\lambda$ is an eigenvalue of $A$ when $A\mathbf{v} = \lambda\mathbf{v}$ for some nonzero $\mathbf{v}$; along $\mathbf{v}$ the matrix acts as multiplication by $\lambda$.
• Eigenvectors are directions: every nonzero multiple of an eigenvector is an eigenvector for the same eigenvalue.
• $\lambda = 0$ is an eigenvalue exactly when $A$ is singular.
• Projections, reflections, and scalings have obvious eigenvectors; rotations of the plane have none that are real.
• Verifying an eigenpair needs one matrix-vector product.

🧪 Practice

1. Without computing anything, give all eigenvalues and eigenvectors of the reflection of $\mathbb{R}^2$ across the line $y = x$, and of the projection of $\mathbb{R}^3$ onto the $xy$-plane.
2. Check whether $(1, 2)$, $(2, 1)$, and $(1, -1)$ are eigenvectors of $A = \begin{bmatrix} 3 & 1 \ 2 & 2 \end{bmatrix}$, and give the eigenvalue for each one that is.
3. Show that if $\mathbf{v}$ is an eigenvector of $A$ with eigenvalue $\lambda$, then it is also an eigenvector of $A^2$, of $A + 5I$, and (if $A$ is invertible) of $A^{-1}$. Find the corresponding eigenvalues.
4. Interview: Explain to a non-specialist what an eigenvector is, and why an engineer designing a footbridge cares about the eigenvalues of a stiffness matrix. Hint: an eigenvector is a pattern of motion that the system reproduces with only a change of scale; the eigenvalues give the natural frequencies, and footsteps that match one of them cause resonance (the Millennium Bridge in London, 2000).

Eigenspaces

Why it matters. An eigenvalue usually has a whole line of eigenvectors, and sometimes a plane or more. Collecting all of them, together with $\mathbf{0}$, gives a subspace on which $A$ is simply "multiply by $\lambda$". Splitting a space into such subspaces is the strategy behind diagonalization (Section 8.2): once you know these pieces and the scalars that act on them, you know $A$.

Definition. For an eigenvalue $\lambda$ of $A$, the eigenspace of $\lambda$ is

$$E_\lambda = {\mathbf{v} : A\mathbf{v} = \lambda\mathbf{v}} = \operatorname{Nul}(A - \lambda I).$$

The rewriting is just algebra: $A\mathbf{v} = \lambda\mathbf{v} \iff A\mathbf{v} - \lambda I\mathbf{v} = \mathbf{0} \iff (A - \lambda I)\mathbf{v} = \mathbf{0}$. Note that $\lambda\mathbf{v}$ must be written as $\lambda I\mathbf{v}$ before factoring, since "$A - \lambda$" (a matrix minus a number) is not defined.

Because $E_\lambda$ is a null space, it is a subspace (Section 6.6). It consists of all the eigenvectors for $\lambda$ together with the zero vector, which is not an eigenvector but is needed for a subspace. Its dimension,

$$\dim E_\lambda = n - \operatorname{rank}(A - \lambda I),$$

is called the geometric multiplicity of $\lambda$ (more in Algebraic vs. Geometric Multiplicity below). The scalar $\lambda$ is an eigenvalue exactly when $E_\lambda \ne {\mathbf{0}}$.

Eigenspaces are invariant. If $\mathbf{v} \in E_\lambda$, then $A\mathbf{v} = \lambda\mathbf{v} \in E_\lambda$. So $A$ maps each eigenspace into itself and acts on it as the scalar $\lambda$. This is the simplest case of an invariant subspace (Section 8.3).

Eigenvectors for different eigenvalues are independent.

Theorem. If $\mathbf{v}_1, \dots, \mathbf{v}_k$ are eigenvectors of $A$ for distinct eigenvalues $\lambda_1, \dots, \lambda_k$, then ${\mathbf{v}_1, \dots, \mathbf{v}_k}$ is linearly independent.

Proof. Suppose not, and let $m$ be the smallest number such that $\mathbf{v}_1, \dots, \mathbf{v}m$ are dependent. Then $\mathbf{v}m = c_1\mathbf{v}1 + \cdots + c{m-1}\mathbf{v}{m-1}$ (Section 6.3), and $\mathbf{v}1, \dots, \mathbf{v}{m-1}$ are independent. Apply $A$ to get $\lambda_m\mathbf{v}m = c_1\lambda_1\mathbf{v}1 + \cdots + c{m-1}\lambda{m-1}\mathbf{v}{m-1}$. Subtract $\lambda_m$ times the first equation:

$$\mathbf{0} = c_1(\lambda_1 - \lambda_m)\mathbf{v}1 + \cdots + c{m-1}(\lambda_{m-1} - \lambda_m)\mathbf{v}_{m-1}.$$

By independence every $c_i(\lambda_i - \lambda_m) = 0$, and since the eigenvalues are distinct, every $c_i = 0$. Then $\mathbf{v}_m = \mathbf{0}$, which is impossible for an eigenvector. $\blacksquare$

Corollary. The sum of the eigenspaces of distinct eigenvalues is a direct sum (Section 6.2): $E_{\lambda_1} + \cdots + E_{\lambda_k} = E_{\lambda_1} \oplus \cdots \oplus E_{\lambda_k}$. Combining bases of the separate eigenspaces always gives an independent set. In particular, an $n \times n$ matrix has at most $n$ distinct eigenvalues.

A worked example. Let

$$A = \begin{bmatrix} 4 & -1 & 6 \ 2 & 1 & 6 \ 2 & -1 & 8 \end{bmatrix}.$$

Its eigenvalues are $2$ and $9$ (the next two topics show how to find them). For $\lambda = 2$,

$$A - 2I = \begin{bmatrix} 2 & -1 & 6 \ 2 & -1 & 6 \ 2 & -1 & 6 \end{bmatrix} \longrightarrow \begin{bmatrix} 2 & -1 & 6 \ 0 & 0 & 0 \ 0 & 0 & 0 \end{bmatrix}.$$

One equation, $2x_1 - x_2 + 6x_3 = 0$, in three unknowns, with $x_2$ and $x_3$ free. So $E_2 = \operatorname{span}{(1, 2, 0), (-3, 0, 1)}$ is a plane. For $\lambda = 9$, $E_9 = \operatorname{span}{(1, 1, 1)}$ is a line. The three basis vectors together are independent, by the corollary.

           E_9 = span{(1, 1, 1)}: every vector here is multiplied by 9
              |
              |       E_2 = a plane through O: every vector here is doubled
          ____|_______________________
         /    |                      /
        /     O ------> v            /      A v = 2 v   (v in E_2)
       /                            /
      /____________________________/

import sympy as sp

A = sp.Matrix([[4, -1, 6],
               [2,  1, 6],
               [2, -1, 8]])
I3 = sp.eye(3)

for lam in [2, 9]:
    basis = (A - lam * I3).nullspace()           # eigenspace E_lam = Nul(A - lam I)
    print(lam, [list(b) for b in basis], "dim", len(basis))
# 2 [[1/2, 1, 0], [-3, 0, 1]] dim 2
# 9 [[1, 1, 1]] dim 1

print((A - 5 * I3).nullspace())                   # []: 5 is not an eigenvalue

# Any combination inside E_2 is again doubled; mixing eigenspaces is not
u = sp.Matrix([1, 2, 0]); w = sp.Matrix([-3, 0, 1]); z = sp.Matrix([1, 1, 1])
print((A * (3 * u - 4 * w)).T, (2 * (3 * u - 4 * w)).T)   # equal: both [30, 12, -8]
print((A * (u + z)).T)                                     # [11, 13, 9]: not a multiple of u + z = [2, 3, 1]

# Bases of different eigenspaces combine into an independent set
print(sp.Matrix.hstack(u, w, z).rank())           # 3

Key takeaways

• $E_\lambda = \operatorname{Nul}(A - \lambda I)$ is a subspace: all eigenvectors for $\lambda$ plus $\mathbf{0}$.
• $A$ maps $E_\lambda$ into itself and acts on it as multiplication by $\lambda$.
• $\dim E_\lambda = n - \operatorname{rank}(A - \lambda I)$ is the geometric multiplicity of $\lambda$.
• Eigenvectors for distinct eigenvalues are linearly independent, so eigenspaces form a direct sum.
• A sum of eigenvectors from different eigenspaces is not an eigenvector.

🧪 Practice

1. Find a basis of each eigenspace of $A = \begin{bmatrix} 3 & 0 & 0 \ 0 & 3 & 1 \ 0 & 0 & 5 \end{bmatrix}$ and state its dimension.
2. Show that if $\mathbf{u}$ and $\mathbf{w}$ are eigenvectors for distinct eigenvalues $\lambda \ne \mu$, then $\mathbf{u} + \mathbf{w}$ is not an eigenvector of $A$.
3. Let $P$ be the matrix of the projection of $\mathbb{R}^3$ onto the plane $x + y + z = 0$. Describe $E_1$ and $E_0$ geometrically and give their dimensions without computing $P$.
4. Interview: You are given 4 eigenvectors of a $4 \times 4$ matrix, each with a different eigenvalue. Can you conclude they form a basis of $\mathbb{R}^4$? What if two of them share an eigenvalue? Hint: distinct eigenvalues force independence; for a shared eigenvalue the two vectors might be parallel, so you must check that pair inside its eigenspace.

Characteristic Polynomial

Why it matters. The definition $A\mathbf{v} = \lambda\mathbf{v}$ has two unknowns, $\lambda$ and $\mathbf{v}$, multiplied together, so it is not a linear system. The determinant (Chapter 5) removes $\mathbf{v}$ from the problem and leaves a single polynomial equation in $\lambda$. This polynomial also explains why an $n \times n$ matrix has at most $n$ eigenvalues and links eigenvalues to the trace and determinant.

Derivation. For a scalar $\lambda$,

$$\lambda \text{ is an eigenvalue} \iff (A - \lambda I)\mathbf{v} = \mathbf{0} \text{ has a nonzero solution} \iff A - \lambda I \text{ is singular} \iff \det(A - \lambda I) = 0.$$

Definition. The characteristic polynomial of an $n \times n$ matrix $A$ is

$$p_A(\lambda) = \det(\lambda I - A).$$

The equation $p_A(\lambda) = 0$ is the characteristic equation, and its roots are the eigenvalues. Some books use $\det(A - \lambda I)$ instead; it equals $(-1)^n p_A(\lambda)$ and has the same roots. This guide uses $\det(\lambda I - A)$ because its leading coefficient is $1$ (it is monic).

Why it is a polynomial of degree $n$. Expanding the determinant as a sum of signed products of entries (Section 5.1), every product is a polynomial in $\lambda$, and only the product of the diagonal entries, $(\lambda - a_{11})(\lambda - a_{22})\cdots(\lambda - a_{nn})$, reaches degree $n$. So $p_A$ has degree exactly $n$ and, over any field, at most $n$ roots.

The $2 \times 2$ case. For $A = \begin{bmatrix} a & b \ c & d \end{bmatrix}$,

$$p_A(\lambda) = \det\begin{bmatrix} \lambda - a & -b \ -c & \lambda - d \end{bmatrix} = \lambda^2 - (a + d)\lambda + (ad - bc) = \lambda^2 - (\operatorname{tr}A)\lambda + \det A.$$

The quadratic formula then gives $\lambda = \tfrac{1}{2}\big(\operatorname{tr}A \pm \sqrt{(\operatorname{tr}A)^2 - 4\det A}\big)$. The sign of the discriminant $(\operatorname{tr}A)^2 - 4\det A$ decides between two real eigenvalues, one repeated eigenvalue, or a complex-conjugate pair.

Trace and determinant in general. Over $\mathbb{C}$, $p_A$ factors completely (Section 1.2): $p_A(\lambda) = (\lambda - \lambda_1)(\lambda - \lambda_2)\cdots(\lambda - \lambda_n)$, with eigenvalues repeated according to multiplicity. Comparing coefficients with the determinant expansion gives

$$p_A(\lambda) = \lambda^n - (\operatorname{tr}A),\lambda^{n-1} + \cdots + (-1)^n\det A,$$

$$\lambda_1 + \lambda_2 + \cdots + \lambda_n = \operatorname{tr}A, \qquad \lambda_1\lambda_2\cdots\lambda_n = \det A.$$

These are fast sanity checks on any eigenvalue computation. (Set $\lambda = 0$ for the constant term: $p_A(0) = \det(-A) = (-1)^n\det A$.)

Facts that follow immediately.

Fact	Reason
A triangular matrix has its diagonal entries as eigenvalues	$\lambda I - A$ is triangular, so $p_A(\lambda) = \prod(\lambda - a_{ii})$ (Section 5.1)
$A$ and $A^T$ have the same eigenvalues	$\det(\lambda I - A^T) = \det\big((\lambda I - A)^T\big) = \det(\lambda I - A)$
Similar matrices have the same characteristic polynomial	$\det(\lambda I - P^{-1}AP) = \det\big(P^{-1}(\lambda I - A)P\big) = \det(\lambda I - A)$
$AB$ and $BA$ have the same eigenvalues (for square $A$, $B$)	true even when neither is invertible; see Practice
$A$ is invertible exactly when $p_A(0) \ne 0$	$p_A(0) = (-1)^n\det A$

$A$ and $A^T$ share eigenvalues but generally not eigenvectors. The eigenvectors of $A^T$ are sometimes called the left eigenvectors of $A$, since $A^T\mathbf{w} = \lambda\mathbf{w}$ means $\mathbf{w}^TA = \lambda\mathbf{w}^T$.

Every monic polynomial is a characteristic polynomial. The companion matrix

$$C = \begin{bmatrix} 0 & 0 & \cdots & 0 & -c_0 \ 1 & 0 & \cdots & 0 & -c_1 \ 0 & 1 & \cdots & 0 & -c_2 \ \vdots & & \ddots & & \vdots \ 0 & 0 & \cdots & 1 & -c_{n-1} \end{bmatrix}$$

has $p_C(\lambda) = \lambda^n + c_{n-1}\lambda^{n-1} + \cdots + c_1\lambda + c_0$. So finding eigenvalues is exactly as hard as finding roots of polynomials. By the Abel-Ruffini theorem there is no formula in radicals for the roots of a general polynomial of degree 5 or more, so eigenvalues of large matrices cannot be computed by a finite formula; every practical algorithm is iterative (Section 12.3). The relationship also runs the other way: numpy.roots finds polynomial roots by building the companion matrix and computing its eigenvalues.

import numpy as np
import sympy as sp

lam = sp.symbols('lambda')
A = sp.Matrix([[4, -1, 6],
               [2,  1, 6],
               [2, -1, 8]])

p = (lam * sp.eye(3) - A).det()                     # det(lambda I - A)
print(sp.expand(p))                                 # lambda**3 - 13*lambda**2 + 40*lambda - 36
print(sp.factor(p))                                 # (lambda - 9)*(lambda - 2)**2
print(A.trace(), A.det())                           # 13 36: matches -(-13) and -(-36)
print(2 + 2 + 9, 2 * 2 * 9)                         # 13 36: sum and product of eigenvalues

# A triangular matrix shows its eigenvalues on the diagonal
T = sp.Matrix([[3, 1, 4], [0, -2, 5], [0, 0, 7]])
print(T.eigenvals())                                # {7: 1, -2: 1, 3: 1}: eigenvalue -> multiplicity

# numpy: np.poly gives characteristic polynomial coefficients from a matrix
An = np.array(A.tolist(), dtype=float)
print(np.round(np.poly(An), 6))                     # [  1. -13.  40. -36.]

# Companion matrix of p(x) = x^3 - 6x^2 + 11x - 6 = (x-1)(x-2)(x-3)
C = np.array([[0.0, 0.0,   6.0],
              [1.0, 0.0, -11.0],
              [0.0, 1.0,   6.0]])
print(np.sort(np.linalg.eigvals(C).real).round(10)) # [1. 2. 3.]
print(np.sort(np.roots([1, -6, 11, -6])).round(10)) # [1. 2. 3.]: same computation inside

Key takeaways

• $\lambda$ is an eigenvalue exactly when $\det(\lambda I - A) = 0$; the determinant eliminates the unknown eigenvector.
• $p_A(\lambda) = \det(\lambda I - A)$ is a monic polynomial of degree $n$, so there are at most $n$ eigenvalues.
• For $2 \times 2$ matrices, $p_A(\lambda) = \lambda^2 - (\operatorname{tr}A)\lambda + \det A$.
• Over $\mathbb{C}$, the eigenvalues (with multiplicity) sum to $\operatorname{tr}A$ and multiply to $\det A$.
• Triangular matrices show their eigenvalues on the diagonal; $A$, $A^T$, and every matrix similar to $A$ share eigenvalues.
• Every monic polynomial is the characteristic polynomial of its companion matrix, so eigenvalues require iterative algorithms for $n \ge 5$.

🧪 Practice

1. Compute the characteristic polynomial and eigenvalues of $\begin{bmatrix} 1 & 4 \ 2 & 3 \end{bmatrix}$ and $\begin{bmatrix} 1 & 2 & 3 \ 0 & 4 & 5 \ 0 & 0 & 6 \end{bmatrix}$. Check the trace and determinant rules.
2. A $3 \times 3$ matrix has trace $6$, determinant $6$, and eigenvalue $1$. Find its other two eigenvalues.
3. Show that $AB$ and $BA$ have the same nonzero eigenvalues: if $AB\mathbf{v} = \lambda\mathbf{v}$ with $\lambda \ne 0$, show that $B\mathbf{v}$ is a nonzero eigenvector of $BA$.
4. Interview: Why do numerical libraries not compute eigenvalues by forming the characteristic polynomial and finding its roots? Hint: tiny relative errors in the coefficients can move the roots a lot (Wilkinson's polynomial, next topic), and for $n \ge 5$ root-finding is itself iterative; libraries go the other way, from polynomial to companion matrix to eigenvalues.

Computing Eigenvalues and Eigenvectors

Why it matters. The previous topics give the two halves of a procedure: the characteristic equation finds the eigenvalues, and a null space computation finds the eigenvectors for each one. Doing small cases by hand builds the intuition needed to read the output of numerical software and to recognize when that output is suspect.

The hand procedure.

1. Compute $p_A(\lambda) = \det(\lambda I - A)$.
2. Find its roots $\lambda_1, \dots, \lambda_k$. These are the eigenvalues.
3. For each $\lambda_i$, row reduce $A - \lambda_iI$ and solve $(A - \lambda_iI)\mathbf{x} = \mathbf{0}$. A basis of the solution set is a basis of $E_{\lambda_i}$.
4. Check: $A\mathbf{v} = \lambda_i\mathbf{v}$ for each basis vector.

In step 3, $A - \lambda_iI$ must be singular, so the row reduction must produce at least one zero row. If it does not, the eigenvalue from step 2 is wrong.

A worked $3 \times 3$ example. Let

$$A = \begin{bmatrix} 1 & 2 & 1 \ 6 & -1 & 0 \ -1 & -2 & -1 \end{bmatrix}.$$

Step 1. Expanding along the first row (and simplifying) gives $p_A(\lambda) = \lambda^3 + \lambda^2 - 12\lambda = \lambda(\lambda - 3)(\lambda + 4)$. Check: $\operatorname{tr}A = -1$ and the coefficient of $\lambda^2$ is $-\operatorname{tr}A = 1$; $\det A = 0$, matching the constant term $0$.

Step 2. The eigenvalues are $0$, $3$, and $-4$.

Step 3. For $\lambda = 3$:

$$A - 3I = \begin{bmatrix} -2 & 2 & 1 \ 6 & -4 & 0 \ -1 & -2 & -4 \end{bmatrix} \longrightarrow \begin{bmatrix} 1 & 0 & 1 \ 0 & 1 & \tfrac{3}{2} \ 0 & 0 & 0 \end{bmatrix}, \qquad E_3 = \operatorname{span}{(2, 3, -2)}.$$

Similarly $E_0 = \operatorname{span}{(1, 6, -13)}$ and $E_{-4} = \operatorname{span}{(1, -2, -1)}$.

Step 4. $A(2, 3, -2) = (2 + 6 - 2,\ 12 - 3 + 0,\ -2 - 6 + 2) = (6, 9, -6) = 3(2, 3, -2)$.

Shortcuts worth knowing.

Situation	Eigen-information
$A$ triangular or diagonal	eigenvalues are the diagonal entries
$2 \times 2$ with $b \ne 0$	for eigenvalue $\lambda$, the vector $(b,\ \lambda - a)$ is an eigenvector
every row of $A$ sums to $s$	$s$ is an eigenvalue with eigenvector $(1, 1, \dots, 1)$
every column of $A$ sums to $s$	$s$ is an eigenvalue (apply the previous row to $A^T$)
$A = \mathbf{u}\mathbf{v}^T$ (rank one)	$\mathbf{v}^T\mathbf{u}$ with eigenvector $\mathbf{u}$; $0$ with multiplicity $n - 1$
block triangular $\begin{bmatrix} B & C \ 0 & D \end{bmatrix}$	eigenvalues of $B$ together with those of $D$ (Section 5.2)

Eigenvalues of related matrices. If $A\mathbf{v} = \lambda\mathbf{v}$, then with the same $\mathbf{v}$:

Matrix	Eigenvalue	Reason
$A^k$	$\lambda^k$	apply $A$ repeatedly
$A^{-1}$ (if it exists)	$1/\lambda$	multiply $A\mathbf{v} = \lambda\mathbf{v}$ by $\lambda^{-1}A^{-1}$
$A + cI$	$\lambda + c$	a shift; used heavily in algorithms (Section 12.3)
$cA$	$c\lambda$	scale both sides
$q(A) = a_0I + a_1A + \cdots + a_mA^m$	$q(\lambda)$	combine the rows above

There is no such rule for $A + B$ or $AB$ unless $A$ and $B$ share eigenvectors.

Numerically, never go through the polynomial. On a computer, eigenvalues are computed by the QR algorithm (Section 12.3), which works directly with the matrix: np.linalg.eig for general matrices and np.linalg.eigh for symmetric ones (faster, and guaranteed real and orthonormal output; Chapter 10). Forming $p_A$ first is unstable. The classic example is $W = \operatorname{diag}(1, 2, \dots, 20)$: its eigenvalues are obvious, but the coefficients of $p_W$ are as large as $20! \approx 2.4 \times 10^{18}$. Merely storing them in double precision moves some roots by $0.07$, and Wilkinson showed that changing the single coefficient $-210$ by $2^{-23} \approx 10^{-7}$ pushes several roots far off the real line.

import numpy as np
import sympy as sp

A = sp.Matrix([[1, 2, 1], [6, -1, 0], [-1, -2, -1]])
lam = sp.symbols('lambda')
print(sp.factor(A.charpoly(lam).as_expr()))        # lambda*(lambda - 3)*(lambda + 4)

# sympy: exact eigenpairs as (eigenvalue, algebraic multiplicity, eigenspace basis)
for val, mult, basis in A.eigenvects():
    v = basis[0] / basis[0][0]                       # scale so the first entry is 1
    print(val, mult, list(v), A * v == val * v)
# -4 1 [1, -2, -1] True
# 0 1 [1, 6, -13] True
# 3 1 [1, 3/2, -1] True

# numpy: floating point, unit-length eigenvectors in the columns
An = np.array(A.tolist(), dtype=float)
vals, vecs = np.linalg.eig(An)
print(np.round(vals.real, 10))                     # [-4.  3.  0.]: no guaranteed order
residual = np.linalg.norm(An @ vecs - vecs * vals) # A V - V diag(vals), column by column
print(residual < 1e-12)                            # True

# The rank-one shortcut: A = u v^T has eigenvalue v.u with eigenvector u
u, v = np.array([1.0, 2.0, 3.0]), np.array([4.0, 0.0, 1.0])
R = np.outer(u, v)
print(np.abs(np.round(np.sort(np.linalg.eigvals(R).real), 10)), v @ u)   # [0. 0. 7.] 7.0

# Wilkinson: the characteristic polynomial route loses the eigenvalues of diag(1..20)
W = np.diag(np.arange(1.0, 21.0))
exact = np.arange(1, 21)
print(np.max(np.abs(np.sort(np.linalg.eigvals(W).real) - exact)))   # 0.0
c = np.poly(W)                                      # coefficients, rounded to doubles
print(round(np.max(np.abs(np.sort(np.roots(c).real) - exact)), 3))  # 0.07: digits already lost
c[1] -= 2.0**-23                                    # change -210 by about 1e-7
print(round(np.max(np.abs(np.roots(c).imag)), 2))   # 2.81: complex "eigenvalues" of a diagonal matrix

Key takeaways

• By hand: compute $p_A$, find its roots, then solve $(A - \lambda I)\mathbf{x} = \mathbf{0}$ for each root, and check the answers.
• If $A - \lambda I$ row reduces without a zero row, $\lambda$ is not an eigenvalue: recheck the arithmetic.
• $A^k$, $A^{-1}$, $A + cI$, and $q(A)$ share the eigenvectors of $A$, with eigenvalues $\lambda^k$, $1/\lambda$, $\lambda + c$, $q(\lambda)$.
• Row sums, rank-one structure, and triangular blocks give eigenvalues without computation.
• In software, use eig or eigh; never compute eigenvalues as roots of a numerically formed characteristic polynomial.

🧪 Practice

1. Find all eigenvalues and a basis of each eigenspace of $\begin{bmatrix} 2 & 0 & 0 \ 1 & 3 & 0 \ 4 & -1 & 1 \end{bmatrix}$ and of $\begin{bmatrix} 5 & 4 \ 1 & 2 \end{bmatrix}$.
2. The matrix $J$ is the $n \times n$ matrix of all ones. Find its eigenvalues and eigenspaces without computing a determinant, then deduce the eigenvalues of $J - I$ and of $(J - I)^{-1}$.
3. Let $A$ have eigenvalues $1$, $2$, $-1$. Find the eigenvalues, trace, and determinant of $A^3 - 2A + I$.
4. Interview: A numerical eigensolver returns eigenvalues for a $1000 \times 1000$ matrix. Describe three inexpensive checks you would run on the output. Hint: compare the sum with the trace; compute residuals $\lVert A\mathbf{v} - \lambda\mathbf{v}\rVert$; for real input, verify that complex eigenvalues come in conjugate pairs; for symmetric input, verify they are real.

Algebraic vs. Geometric Multiplicity

Why it matters. When an eigenvalue is a repeated root of the characteristic polynomial, it may or may not have a correspondingly large eigenspace. The identity matrix $I_2$ and the shear $\begin{bmatrix} 1 & 1 \ 0 & 1 \end{bmatrix}$ both have characteristic polynomial $(\lambda - 1)^2$, yet every vector is an eigenvector of $I_2$, while the shear has only one eigen-direction. Telling the two situations apart decides whether a matrix can be diagonalized (Section 8.2), and when it cannot, the missing eigenvectors are what the Jordan form (Section 8.3) is built to handle.

Definitions. Let $\lambda$ be an eigenvalue of $A$.

• The algebraic multiplicity $\operatorname{am}(\lambda)$ is the multiplicity of $\lambda$ as a root of $p_A$: the exponent $m$ in the factor $(t - \lambda)^m$ of $p_A(t)$.
• The geometric multiplicity $\operatorname{gm}(\lambda) = \dim E_\lambda = n - \operatorname{rank}(A - \lambda I)$ is the number of independent eigenvectors for $\lambda$.

Analogy. The algebraic multiplicity is the number of seats reserved for $\lambda$ in the characteristic polynomial. The geometric multiplicity is how many independent eigenvectors actually show up to fill them. Some seats can stay empty, but there can never be more guests than seats.

Theorem. For every eigenvalue, $1 \le \operatorname{gm}(\lambda) \le \operatorname{am}(\lambda)$.

Proof. $\operatorname{gm} \ge 1$ because an eigenvalue has at least one eigenvector. For the upper bound, let $g = \operatorname{gm}(\lambda)$, pick a basis $\mathbf{v}_1, \dots, \mathbf{v}g$ of $E\lambda$, and extend it to a basis of $F^n$ (Section 6.4). Let $P$ have these basis vectors as columns. The first $g$ columns of $P^{-1}AP$ are $\lambda\mathbf{e}_1, \dots, \lambda\mathbf{e}_g$, so

$$P^{-1}AP = \begin{bmatrix} \lambda I_g & B \ 0 & C \end{bmatrix}, \qquad p_A(t) = p_{P^{-1}AP}(t) = (t - \lambda)^g,p_C(t),$$

using the block triangular determinant (Section 5.2) and the fact that similar matrices have the same characteristic polynomial. So $(t - \lambda)^g$ divides $p_A$, and $\operatorname{am}(\lambda) \ge g$. $\blacksquare$

An eigenvalue with $\operatorname{gm} < \operatorname{am}$ is called defective, and so is a matrix that has one. An eigenvalue with $\operatorname{am} = 1$ is simple, and then automatically $\operatorname{gm} = 1$.

Examples.

Matrix	Eigenvalue	am	gm	Defective?
$\begin{bmatrix} 1 & 0 \ 0 & 1 \end{bmatrix}$	$1$	2	2	no
$\begin{bmatrix} 1 & 1 \ 0 & 1 \end{bmatrix}$	$1$	2	1	yes
$\begin{bmatrix} 4 & -1 & 6 \ 2 & 1 & 6 \ 2 & -1 & 8 \end{bmatrix}$	$2$; $9$	2; 1	2; 1	no
$\begin{bmatrix} 2 & 1 & 0 \ 0 & 2 & 1 \ 0 & 0 & 2 \end{bmatrix}$	$2$	3	1	yes
$\begin{bmatrix} 2 & 1 & 0 \ 0 & 2 & 0 \ 0 & 0 & 2 \end{bmatrix}$	$2$	3	2	yes
$\begin{bmatrix} 0 & 1 \ 0 & 0 \end{bmatrix}$	$0$	2	1	yes

Over $\mathbb{C}$ the algebraic multiplicities add up to $n$, because $p_A$ splits into $n$ linear factors. The geometric multiplicities add up to at most $n$, with equality exactly when there is a basis of eigenvectors (Section 8.2).

Defective matrices are fragile. Perturbing the shear by a tiny $\varepsilon$ in its lower-left corner gives $\begin{bmatrix} 1 & 1 \ \varepsilon & 1 \end{bmatrix}$, with eigenvalues $1 \pm \sqrt{\varepsilon}$. A change of $10^{-10}$ in an entry moves the eigenvalues by $10^{-5}$, a hundred thousand times more. This is why floating-point software essentially never reports exactly repeated eigenvalues for a defective matrix, and why multiplicities are a question for exact arithmetic or for careful numerical rank decisions (Chapter 11).

   I_2: am = 2, gm = 2                  shear: am = 2, gm = 1

        y                                    y
        ^   every direction                  |   only the x-axis
        |   is an eigen-direction            |   is an eigen-direction
     \  |  /                                 |
      \ | /                                  |     points slide right
   ----\|/----> x                    --------O=========> x   (eigenvalue 1)
       /|\                                   |    ----> points above the axis
      / | \                                  |

import numpy as np
import sympy as sp

def multiplicities(M):
    """Return {eigenvalue: (algebraic, geometric)} using exact arithmetic."""
    n = M.shape[0]
    return {lam: (am, n - (M - lam * sp.eye(n)).rank()) for lam, am in M.eigenvals().items()}

print(multiplicities(sp.eye(2)))                                    # {1: (2, 2)}
print(multiplicities(sp.Matrix([[1, 1], [0, 1]])))                  # {1: (2, 1)}
print(multiplicities(sp.Matrix([[2, 1, 0], [0, 2, 1], [0, 0, 2]]))) # {2: (3, 1)}
print(multiplicities(sp.Matrix([[2, 1, 0], [0, 2, 0], [0, 0, 2]]))) # {2: (3, 2)}
print(multiplicities(sp.Matrix([[4, -1, 6], [2, 1, 6], [2, -1, 8]])))  # {9: (1, 1), 2: (2, 2)}

# A defective eigenvalue splits violently under perturbation: about sqrt(eps)
for eps in [1e-4, 1e-10]:
    S = np.array([[1.0, 1.0], [eps, 1.0]])
    print(eps, np.sort(np.linalg.eigvals(S).real) - 1)
# 0.0001 [-0.01  0.01]
# 1e-10 [-1.e-05  1.e-05]

Key takeaways

• Algebraic multiplicity counts repeated roots of $p_A$; geometric multiplicity counts independent eigenvectors, $\dim E_\lambda = n - \operatorname{rank}(A - \lambda I)$.
• Always $1 \le \operatorname{gm}(\lambda) \le \operatorname{am}(\lambda)$; simple eigenvalues have both equal to 1.
• A matrix is defective when some eigenvalue has fewer eigenvectors than its algebraic multiplicity; such a matrix has no basis of eigenvectors.
• Matrices with the same characteristic polynomial can differ in geometric multiplicity, so they need not be similar.
• Repeated eigenvalues of defective matrices are extremely sensitive to perturbation (roughly $\sqrt{\varepsilon}$ for a $2 \times 2$ block).

🧪 Practice

1. Find the algebraic and geometric multiplicity of each eigenvalue of $\begin{bmatrix} 3 & 1 & 0 \ 0 & 3 & 0 \ 0 & 0 & 5 \end{bmatrix}$ and of $\begin{bmatrix} 3 & 0 & 1 \ 0 & 3 & 0 \ 0 & 0 & 3 \end{bmatrix}$.
2. For which values of $a$ is $\begin{bmatrix} 2 & a \ 0 & 2 \end{bmatrix}$ defective? For which values of $b$ does $\begin{bmatrix} 1 & b \ 1 & 1 \end{bmatrix}$ have a repeated eigenvalue, and is it defective then?
3. Show that a nonzero nilpotent matrix ($N^k = 0$ for some $k$) is always defective. Start by showing that $0$ is its only eigenvalue.
4. Interview: An eigensolver returns eigenvalues $0.99999$ and $1.00001$ for a $2 \times 2$ matrix with nearly parallel eigenvectors. What might the exact matrix look like, and why are the computed eigenvectors unreliable? Hint: think of a shear perturbed by $\varepsilon \approx 10^{-10}$; the true matrix may be defective, and near a defective matrix the eigenvectors almost coincide, so tiny errors change them a lot.

Complex Eigenvalues

Why it matters. A rotation of the plane by $90°$ moves every nonzero vector off its line, so it has no real eigenvectors. Yet rotations, oscillations, and spirals are among the most common behaviors of linear systems: a mass on a spring, an AC circuit, a predator-prey cycle. Allowing complex eigenvalues restores the theory (every matrix has eigenvalues over $\mathbb{C}$), and the complex eigenvalue of a real matrix encodes a rotation angle and a scaling factor.

Existence over $\mathbb{C}$. By the fundamental theorem of algebra (Section 1.2), $p_A$ has exactly $n$ complex roots counted with multiplicity. So every square matrix, real or complex, has at least one complex eigenvalue and eigenvector. Over $\mathbb{R}$ this can fail; over $\mathbb{C}$ it never does. This existence is the first step of Schur's theorem and the Jordan form (Section 8.3).

The rotation matrix. For $R_\theta = \begin{bmatrix} \cos\theta & -\sin\theta \ \sin\theta & \cos\theta \end{bmatrix}$,

$$p(\lambda) = \lambda^2 - 2\cos\theta,\lambda + 1, \qquad \lambda = \cos\theta \pm i\sin\theta = e^{\pm i\theta},$$

with eigenvectors $(1, \mp i)$. Both eigenvalues have modulus 1, as they should for a map that preserves length, and their argument is the rotation angle.

Real matrices: conjugate pairs. If $A$ is real and $A\mathbf{v} = \lambda\mathbf{v}$, conjugating every entry gives $A\bar{\mathbf{v}} = \bar\lambda\bar{\mathbf{v}}$, since $\bar A = A$. So non-real eigenvalues of a real matrix come in conjugate pairs $a \pm bi$, with conjugate eigenvectors. Consequences:

• A real matrix of odd size has at least one real eigenvalue (a real cubic always crosses zero).
• Real symmetric matrices have only real eigenvalues (Chapter 10), so complex eigenvalues signal a non-symmetric, rotational component.

The rotation-scaling form. The real matrix $C = \begin{bmatrix} a & -b \ b & a \end{bmatrix}$ has eigenvalues $a \pm bi$. Writing $a + bi = r(\cos\varphi + i\sin\varphi)$ in polar form (Section 1.2),

$$C = r\begin{bmatrix} \cos\varphi & -\sin\varphi \ \sin\varphi & \cos\varphi \end{bmatrix}, \qquad r = \lvert\lambda\rvert = \sqrt{a^2 + b^2}, \quad \varphi = \arg\lambda.$$

So $C$ rotates by $\varphi$ and scales by $r$. This is exactly the matrix representation of multiplication by $a + bi$ on $\mathbb{C} \cong \mathbb{R}^2$ (Section 1.2).

Theorem. Let $A$ be a real $2 \times 2$ matrix with eigenvalue $\lambda = a - bi$ ($b \ne 0$) and eigenvector $\mathbf{v} \in \mathbb{C}^2$. Then

$$A = PCP^{-1}, \qquad P = \big[,\operatorname{Re}\mathbf{v}\ \ \operatorname{Im}\mathbf{v},\big], \qquad C = \begin{bmatrix} a & -b \ b & a \end{bmatrix}.$$

Proof sketch. Write $\mathbf{v} = \mathbf{x} + i\mathbf{y}$ and expand $A(\mathbf{x} + i\mathbf{y}) = (a - bi)(\mathbf{x} + i\mathbf{y})$. Real parts give $A\mathbf{x} = a\mathbf{x} + b\mathbf{y}$ and imaginary parts give $A\mathbf{y} = -b\mathbf{x} + a\mathbf{y}$. These say $AP = PC$ column by column. Since $b \ne 0$, $\mathbf{x}$ and $\mathbf{y}$ are independent (otherwise $\mathbf{v}$ would be a complex multiple of a real vector, which would make $\lambda$ real), so $P$ is invertible. $\blacksquare$

In words: in the real coordinates given by $\operatorname{Re}\mathbf{v}$ and $\operatorname{Im}\mathbf{v}$, every real $2 \times 2$ matrix with non-real eigenvalues is a rotation combined with a scaling. Iterating it traces an ellipse (if $r = 1$), an inward spiral ($r < 1$), or an outward spiral ($r > 1$).

$\lvert\lambda\rvert$	Orbit of $\mathbf{x}_0, A\mathbf{x}_0, A^2\mathbf{x}_0, \dots$ (real $2 \times 2$, non-real $\lambda$)
$< 1$	spirals in toward $\mathbf{0}$, turning by about $\arg\lambda$ per step
$= 1$	stays on a fixed ellipse forever
$> 1$	spirals out to infinity

import numpy as np

theta = np.pi / 2                                     # rotation by 90 degrees
R = np.array([[np.cos(theta), -np.sin(theta)],
              [np.sin(theta),  np.cos(theta)]])
vals, vecs = np.linalg.eig(R)
print(np.round(vals, 10))                             # [0.+1.j 0.-1.j]: e^{+i pi/2}, e^{-i pi/2}
print(np.abs(vals))                                   # [1. 1.]: rotations preserve length

# A real matrix with a conjugate pair of eigenvalues
A = np.array([[0.5, -0.6],
              [0.75, 1.1]])
vals, vecs = np.linalg.eig(A)
print(np.round(vals, 6))                              # [0.8+0.6j 0.8-0.6j]
lam, v = vals[1], vecs[:, 1]                          # take lambda = a - bi with b > 0
a, b = lam.real, -lam.imag
P = np.column_stack([v.real, v.imag])
C = np.array([[a, -b], [b, a]])
print(np.allclose(A, P @ C @ np.linalg.inv(P)))       # True: A = P C P^{-1}
print(abs(lam), np.degrees(np.angle(vals[0])))        # 1.0 36.86989764584402: rotate ~36.9 deg, no scaling

# Iterating: |lambda| = 1 keeps orbits bounded; scaling A by 0.9 gives an inward spiral
x = np.array([1.0, 0.0])
for k in range(4):
    x = (0.9 * A) @ x
    print(k + 1, np.round(x, 3), round(np.linalg.norm(x), 3))
# 1 [0.45  0.675] 0.811
# 2 [-0.162  0.972] 0.985
# 3 [-0.598  0.853] 1.042
# 4 [-0.73   0.441] 0.852

The norm in the last loop is not monotone because $P$ is not orthogonal: the orbit is an ellipse being shrunk, not a circle. Over every 10 steps it shrinks by a factor $0.9^{10} \approx 0.35$ in the $P$-coordinates.

Key takeaways

• Over $\mathbb{C}$ every $n \times n$ matrix has $n$ eigenvalues counted with multiplicity; over $\mathbb{R}$ it may have none.
• Rotation by $\theta$ has eigenvalues $e^{\pm i\theta}$, both of modulus 1.
• Non-real eigenvalues of a real matrix come in conjugate pairs with conjugate eigenvectors; odd-size real matrices have a real eigenvalue.
• A real $2 \times 2$ matrix with eigenvalue $a - bi$ equals $PCP^{-1}$ with $P = [\operatorname{Re}\mathbf{v}\ \operatorname{Im}\mathbf{v}]$ and $C$ a rotation-scaling by angle $\arg\lambda$ and factor $\lvert\lambda\rvert$.
• $\lvert\lambda\rvert < 1$, $= 1$, $> 1$ give inward spirals, closed ellipses, and outward spirals.

🧪 Practice

1. Find the eigenvalues and eigenvectors of $\begin{bmatrix} 1 & -2 \ 1 & 3 \end{bmatrix}$. Give the rotation angle and scaling factor, and say whether iterates spiral in or out.
2. Show that $\begin{bmatrix} a & -b \ b & a \end{bmatrix}$ has eigenvector $(1, -i)$ with eigenvalue $a + bi$, and find the eigenvector for $a - bi$ by conjugation.
3. Explain why a real $3 \times 3$ matrix cannot have eigenvalues $i$, $2i$, and $1$, and why a real $3 \times 3$ rotation always has a real axis.
4. Interview: In a recurrent neural network, the hidden state evolves roughly as $\mathbf{h}_{k+1} = W\mathbf{h}_k$. What do complex eigenvalues of $W$ with modulus slightly above 1 imply for long sequences, and how would you detect the problem? Hint: oscillating components that grow like $\lvert\lambda\rvert^k$ (exploding gradients); inspect the spectral radius $\max\lvert\lambda_i\rvert$ and keep it at or below 1.

<a id="82-diagonalization"></a>

8.2 Diagonalization

This subchapter answers the central question raised in Section 7.2: when can a basis be chosen so that the matrix of an operator is diagonal, how is that basis found, and what does a diagonal form make easy?

Similar Matrices and Invariants

Why it matters. Section 7.2 showed that similar matrices, $B = P^{-1}AP$, are the same operator written in two different bases. Diagonalization asks whether $A$ is similar to a diagonal matrix, so it is essential to know which properties similarity preserves and which it changes. Invariants also give fast proofs that two matrices are not similar, and they explain why eigenvalues are properties of the operator rather than of a particular matrix.

Analogy. A temperature can be recorded in Celsius or Fahrenheit. The numbers differ, but facts such as "it is above freezing" or "it rose overnight" do not depend on the scale. Similar matrices are the same operator recorded in different coordinate systems. A similarity invariant is a fact about the operator that every coordinate system reports in the same way.

Eigen-information under similarity. Let $B = P^{-1}AP$.

• Same characteristic polynomial (Section 8.1), hence the same eigenvalues with the same algebraic multiplicities, the same trace, and the same determinant.
• Eigenvectors are translated, not preserved. If $A\mathbf{v} = \lambda\mathbf{v}$, then $B(P^{-1}\mathbf{v}) = P^{-1}AP,P^{-1}\mathbf{v} = P^{-1}A\mathbf{v} = \lambda,P^{-1}\mathbf{v}$. So $P^{-1}\mathbf{v}$ is an eigenvector of $B$: it is the coordinate vector of the same eigenvector in the new basis.
• Same geometric multiplicities. $B - \lambda I = P^{-1}(A - \lambda I)P$, and multiplying by invertible matrices does not change rank, so $\dim\operatorname{Nul}(B - \lambda I) = \dim\operatorname{Nul}(A - \lambda I)$.
• More generally, $(B - \lambda I)^k = P^{-1}(A - \lambda I)^kP$, so $\operatorname{rank}(A - \lambda I)^k = \operatorname{rank}(B - \lambda I)^k$ for every $\lambda$ and every $k$.

Invariants, from coarse to complete.

Invariant	Detects	Fails to separate
trace, determinant	quick mismatches	many non-similar pairs
characteristic polynomial	the eigenvalues with algebraic multiplicities	$I_2$ and the shear $\begin{bmatrix} 1 & 1 \ 0 & 1 \end{bmatrix}$
geometric multiplicities $\dim E_\lambda$	$I_2$ versus the shear	some $4 \times 4$ nilpotent pairs (Practice)
minimal polynomial (Section 8.3)	the size of the largest Jordan block per eigenvalue	some larger nilpotent pairs
ranks of $(A - \lambda I)^k$ for all $\lambda$, $k$	everything	nothing: this is a complete invariant

The last row is the Jordan canonical form in disguise (Section 8.3): two complex matrices are similar if and only if all these ranks agree. For diagonalizable matrices the test is much simpler.

Theorem. If $A$ and $B$ are both diagonalizable, then $A \sim B$ if and only if they have the same eigenvalues with the same multiplicities.

Proof. If $A = PDP^{-1}$ and $B = QD'Q^{-1}$ with $D'$ a reordering of $D$, then $D' = R^{-1}DR$ for a permutation matrix $R$, and similarity is transitive. The converse holds because eigenvalues are invariants. $\blacksquare$

Similarity does not preserve everything. Entries, symmetry, orthogonality of eigenvectors, and norms such as $\lVert A\rVert$ can all change. For example, $\begin{bmatrix} 1 & 0 \ 0 & 2 \end{bmatrix}$ and $\begin{bmatrix} 1 & 100 \ 0 & 2 \end{bmatrix}$ are similar (both diagonalizable with eigenvalues 1 and 2), yet the second has much larger entries and nearly parallel eigenvectors. Transient behavior of $A^k\mathbf{x}$, before the eigenvalues take over, depends on these non-invariant features (Section 8.4).

   one operator T, many matrices                  what all of them share

         basis B1          basis B2               eigenvalues, am and gm,
     [T]_B1 = A  <---->  [T]_B2 = P^-1 A P        trace, det, rank,
            \                  /                  ranks of (A - lambda I)^k
             \     basis B3   /
              [T]_B3 = Q^-1 A Q                   what differs: entries,
                                                  symmetry, eigenvector coordinates

import sympy as sp

def similarity_profile(M):
    """Ranks of (M - lam I)^k for each eigenvalue and k = 1..n: a complete similarity invariant."""
    n = M.shape[0]
    return {lam: [((M - lam * sp.eye(n)) ** k).rank() for k in range(1, n + 1)]
            for lam in sorted(M.eigenvals())}

I2 = sp.eye(2)
S = sp.Matrix([[1, 1], [0, 1]])
print(I2.charpoly().as_expr() == S.charpoly().as_expr())   # True: same characteristic polynomial
print(similarity_profile(I2), similarity_profile(S))       # {1: [0, 0]} {1: [1, 0]}: not similar

# Same eigenvalues and both diagonalizable -> similar, even though the entries look different
A = sp.Matrix([[1, 0], [0, 2]])
B = sp.Matrix([[1, 100], [0, 2]])
print(similarity_profile(A) == similarity_profile(B))      # True
Q, _ = B.diagonalize()                                     # columns of Q: eigenvectors of B
print(Q.inv() * B * Q)                                     # Matrix([[1, 0], [0, 2]])

# Eigenvectors transform by P^{-1}
P = sp.Matrix([[2, 1], [1, 1]])
C = P.inv() * B * P
v = sp.Matrix([100, 1])                                    # eigenvector of B for eigenvalue 2
print(B * v == 2 * v, C * (P.inv() * v) == 2 * (P.inv() * v))   # True True

Key takeaways

• Similar matrices share the characteristic polynomial, eigenvalues, algebraic and geometric multiplicities, trace, determinant, and rank.
• If $A\mathbf{v} = \lambda\mathbf{v}$ and $B = P^{-1}AP$, then $P^{-1}\mathbf{v}$ is the corresponding eigenvector of $B$.
• The ranks of $(A - \lambda I)^k$ for all eigenvalues and powers form a complete invariant for similarity over $\mathbb{C}$.
• Two diagonalizable matrices are similar exactly when they have the same eigenvalues with the same multiplicities.
• Entries, symmetry, and eigenvector angles are not invariant, and they control short-term behavior.

🧪 Practice

1. Show that $\begin{bmatrix} 2 & 1 \ 0 & 3 \end{bmatrix}$ and $\begin{bmatrix} 3 & 0 \ 4 & 2 \end{bmatrix}$ are similar, and find an invertible $P$ with $P^{-1}AP = B$.
2. Prove that $\begin{bmatrix} 2 & 1 \ 0 & 2 \end{bmatrix}$ and $\begin{bmatrix} 2 & 0 \ 0 & 2 \end{bmatrix}$ are not similar using a single invariant.
3. Let $N_1$ have a Jordan structure of blocks of sizes $(2, 2)$ and $N_2$ of sizes $(3, 1)$, both with eigenvalue 0. Show that both have $\dim\operatorname{Nul} = 2$, and find the invariant that separates them.
4. Interview: Two engineers measure the same linear system in different sensor coordinates and get different matrices. How can they check that they modeled the same dynamics without finding the change of basis? Hint: compare eigenvalues (with multiplicities) and, if repeated eigenvalues occur, the ranks of $(A - \lambda I)^k$; numerically, use tolerances, since exact equality of computed eigenvalues is not reliable.

Diagonalizability Criteria

Why it matters. A diagonal matrix is the simplest possible operator: it scales each coordinate independently, and its powers, inverse, exponential, and every other function are computed entry by entry (Section 4.2). Knowing when a matrix can be brought to this form tells you when every problem about it reduces to $n$ separate scalar problems.

Definition. A square matrix $A$ is diagonalizable (over $F$) if $A = PDP^{-1}$ for some invertible $P$ and diagonal $D$ with entries in $F$. Equivalently, $A$ is similar to a diagonal matrix, or the operator $\mathbf{x} \mapsto A\mathbf{x}$ has a diagonal matrix in some basis.

The diagonalization theorem. An $n \times n$ matrix $A$ is diagonalizable if and only if $A$ has $n$ linearly independent eigenvectors. In that case $A = PDP^{-1}$, where the columns of $P$ are $n$ independent eigenvectors and the diagonal entries of $D$ are the corresponding eigenvalues, in the same order.

Proof. Write $P = [,\mathbf{p}_1 \cdots \mathbf{p}_n,]$ and $D = \operatorname{diag}(d_1, \dots, d_n)$. Column by column,

$$AP = [,A\mathbf{p}_1 \cdots A\mathbf{p}_n,], \qquad PD = [,d_1\mathbf{p}_1 \cdots d_n\mathbf{p}_n,].$$

So $AP = PD$ says exactly that $A\mathbf{p}_i = d_i\mathbf{p}_i$ for every $i$. $P$ is invertible exactly when its columns are independent, and then $AP = PD$ is equivalent to $A = PDP^{-1}$. Independent columns are nonzero, so each $\mathbf{p}_i$ is an eigenvector. $\blacksquare$

Equivalent and sufficient conditions.

Condition	Status
$A$ has $n$ linearly independent eigenvectors (an eigenbasis)	equivalent
$\sum_\lambda \dim E_\lambda = n$	equivalent
$p_A$ splits into linear factors over $F$ and $\operatorname{gm}(\lambda) = \operatorname{am}(\lambda)$ for every $\lambda$	equivalent
the minimal polynomial has no repeated roots and splits over $F$ (Section 8.3)	equivalent
$A$ has $n$ distinct eigenvalues in $F$	sufficient, not necessary ($I$ is diagonal)
$A$ is real symmetric	sufficient, even with an orthonormal eigenbasis (Chapter 10)
$A$ satisfies $A^2 = A$ (projection) or $A^2 = I$ (involution)	sufficient (Section 8.3)

The $n$-distinct-eigenvalues test follows from Section 8.1: one eigenvector per eigenvalue gives $n$ independent vectors.

The field matters. The rotation $\begin{bmatrix} 0 & -1 \ 1 & 0 \end{bmatrix}$ has no real eigenvalues, so it is not diagonalizable over $\mathbb{R}$. Over $\mathbb{C}$ it has the distinct eigenvalues $\pm i$ and is diagonalizable. Defectiveness, in contrast, cannot be fixed by enlarging the field: the shear has the single eigenvalue $1$ with one eigen-direction over every field.

Diagonalizability is unrelated to invertibility.

	Invertible	Singular
Diagonalizable	$\begin{bmatrix} 2 & 0 \ 0 & 3 \end{bmatrix}$	$\begin{bmatrix} 1 & 0 \ 0 & 0 \end{bmatrix}$
Not diagonalizable	$\begin{bmatrix} 1 & 1 \ 0 & 1 \end{bmatrix}$	$\begin{bmatrix} 0 & 1 \ 0 & 0 \end{bmatrix}$

Invertibility is about whether $0$ is an eigenvalue. Diagonalizability is about whether there are enough eigenvectors.

Most matrices are diagonalizable. Repeated eigenvalues require the discriminant of $p_A$ to vanish, a single polynomial equation in the entries. A random matrix, or one perturbed by rounding errors, almost surely has $n$ distinct complex eigenvalues and is therefore diagonalizable over $\mathbb{C}$. Defective matrices still matter: they arise from structure (shears, nilpotent operators, repeated roots in differential equations), and matrices near them have badly conditioned eigenvectors.

   Is the n x n matrix A diagonalizable over F?

   Does p_A split into linear factors over F? --- no ---> NOT diagonalizable over F
            |                                             (try a larger field, e.g. C)
           yes
            |
   n distinct eigenvalues? --- yes ---> diagonalizable
            |
            no
            |
   for each repeated lambda: n - rank(A - lambda I) == am(lambda)?
            |                                  |
           all yes                         some no
            |                                  |
     diagonalizable                    NOT diagonalizable (defective):
                                       use the Jordan form, Section 8.3

import sympy as sp

def diagonalizable(M):
    """Exact test: gm(lambda) == am(lambda) for every eigenvalue (over C)."""
    n = M.shape[0]
    return all(n - (M - lam * sp.eye(n)).rank() == am for lam, am in M.eigenvals().items())

tests = {
    "distinct eigenvalues": sp.Matrix([[1, 2], [3, 4]]),
    "identity":             sp.eye(3),
    "shear":                sp.Matrix([[1, 1], [0, 1]]),
    "rotation by 90 deg":   sp.Matrix([[0, -1], [1, 0]]),
    "repeated, enough eigenvectors": sp.Matrix([[4, -1, 6], [2, 1, 6], [2, -1, 8]]),
    "nilpotent":            sp.Matrix([[0, 1, 0], [0, 0, 1], [0, 0, 0]]),
}
for name, M in tests.items():
    print(f"{name:30s} {diagonalizable(M)!s:6s} sympy: {M.is_diagonalizable()}, "
          f"over R: {M.is_diagonalizable(reals_only=True)}")
# distinct eigenvalues           True   sympy: True, over R: True
# identity                       True   sympy: True, over R: True
# shear                          False  sympy: False, over R: False
# rotation by 90 deg             True   sympy: True, over R: False
# repeated, enough eigenvectors  True   sympy: True, over R: True
# nilpotent                      False  sympy: False, over R: False

Key takeaways

• $A$ is diagonalizable exactly when it has $n$ linearly independent eigenvectors; then $A = PDP^{-1}$ with eigenvectors in $P$ and eigenvalues in $D$.
• Equivalently, $p_A$ splits and every eigenvalue has $\operatorname{gm} = \operatorname{am}$.
• $n$ distinct eigenvalues guarantee diagonalizability; repeated eigenvalues require checking the eigenspace dimensions.
• Diagonalizability depends on the field (rotations need $\mathbb{C}$) but defectiveness does not go away over a larger field.
• Diagonalizable and invertible are independent properties; generic matrices are diagonalizable over $\mathbb{C}$.

🧪 Practice

1. Decide which matrices are diagonalizable over $\mathbb{R}$, and justify: $\begin{bmatrix} 3 & 1 \ 0 & 2 \end{bmatrix}$, $\begin{bmatrix} 3 & 1 \ 0 & 3 \end{bmatrix}$, $\begin{bmatrix} 1 & -1 \ 1 & 1 \end{bmatrix}$, $\begin{bmatrix} 2 & 0 & 0 \ 0 & 2 & 0 \ 1 & 0 & 2 \end{bmatrix}$.
2. For which values of $a$ and $b$ is $\begin{bmatrix} 1 & a & b \ 0 & 1 & 0 \ 0 & 0 & 2 \end{bmatrix}$ diagonalizable?
3. Show that if $A$ is diagonalizable and invertible, then $A^{-1}$ is diagonalizable, and that if $A$ is diagonalizable then so is $A^T$.
4. Interview: Your code calls np.linalg.eig, gets $P$ and $D$, and reconstructs $A$ as P @ D @ inv(P) with a large error. The matrix came from a model with a repeated root. What is going on, and what would you use instead? Hint: $A$ is probably defective or nearly so, which makes $P$ nearly singular (check np.linalg.cond(P)); use the Schur decomposition (Section 8.3) or the specific function you need, such as scipy.linalg.expm.

Eigenbasis Construction

Why it matters. A basis made of eigenvectors, an eigenbasis, is a coordinate system in which the operator does not mix coordinates at all. In eigen-coordinates, a coupled system of equations becomes a list of independent scalar equations. Building one is the practical content of diagonalization, and it works for abstract operators (on polynomials, matrices, or functions) as well as for matrices.

Analogy. A mixing desk where every slider affects several speakers at once is hard to tune. If you rewire it so that each slider controls exactly one speaker, every adjustment is independent. An eigenbasis is that rewiring: each eigen-coordinate is affected only by its own eigenvalue.

Construction.

1. Find the distinct eigenvalues $\lambda_1, \dots, \lambda_k$.
2. For each $\lambda_i$, find a basis $\mathcal{B}i$ of $E{\lambda_i} = \operatorname{Nul}(A - \lambda_iI)$.
3. Concatenate: $\mathcal{B} = \mathcal{B}_1 \cup \cdots \cup \mathcal{B}_k$.
4. If $\mathcal{B}$ has $n$ vectors, it is an eigenbasis. If it has fewer, no eigenbasis exists.

Step 3 needs no independence check: the eigenspaces form a direct sum (Section 8.1), so combining bases of different eigenspaces always gives an independent set. Within an eigenspace of dimension 2 or more, any basis can be used, which is one reason eigenbases are not unique. Eigenvectors can also be rescaled and reordered freely.

Decoupling. If $\mathcal{B} = (\mathbf{v}_1, \dots, \mathbf{v}_n)$ is an eigenbasis with $A\mathbf{v}_i = \lambda_i\mathbf{v}_i$ and $\mathbf{x} = c_1\mathbf{v}_1 + \cdots + c_n\mathbf{v}_n$, then

$$A\mathbf{x} = \lambda_1c_1\mathbf{v}1 + \cdots + \lambda_nc_n\mathbf{v}n, \qquad [A\mathbf{x}]{\mathcal{B}} = \operatorname{diag}(\lambda_1, \dots, \lambda_n),[\mathbf{x}]{\mathcal{B}}.$$

Each coordinate $c_i$ is multiplied by its own $\lambda_i$, independently of the others.

     standard coordinates                       eigen-coordinates

   x1' = 4 x1 -   x2 + 6 x3                 c1' = 2 c1
   x2' = 2 x1 +   x2 + 6 x3       P^-1      c2' = 2 c2
   x3' = 2 x1 -   x2 + 8 x3     -------->   c3' = 9 c3

   every output depends on every input      each coordinate evolves alone

Eigenbases for abstract operators. The construction does not need matrices with numbers in them, only an operator $T : V \to V$.

Operator	Eigenbasis	Eigenvalues
$T(p) = x,p'(x)$ on $\mathbb{P}_n$	$1, x, x^2, \dots, x^n$	$0, 1, 2, \dots, n$ (since $x \cdot kx^{k-1} = kx^k$)
$T(M) = M^T$ on $\mathbb{R}^{n \times n}$	symmetric matrices, then skew-symmetric matrices	$1$ with multiplicity $\tfrac{n(n+1)}{2}$; $-1$ with multiplicity $\tfrac{n(n-1)}{2}$
a projection $P$ onto $U$ along $W$	a basis of $U$, then a basis of $W$	$1$ on $U$, $0$ on $W$ (Section 7.1)
$D = \tfrac{d}{dx}$ on $\mathbb{P}_n$ ($n \ge 1$)	none	only $0$, with eigenspace the constants: defective

The transpose example explains the decomposition of every square matrix into symmetric and skew-symmetric parts (Section 6.2): it is the eigenspace decomposition of $M \mapsto M^T$.

import numpy as np
import sympy as sp

A = sp.Matrix([[4, -1, 6], [2, 1, 6], [2, -1, 8]])

# Concatenate eigenspace bases; independence across eigenspaces is automatic
basis = [v for lam, _, vs in A.eigenvects() for v in vs]
P = sp.Matrix.hstack(*basis)
print(P.shape[1] == 3, P.rank())                   # True 3: an eigenbasis

# Coordinates in the eigenbasis are scaled independently
x = sp.Matrix([1, 0, 0])
c = P.solve(x)                                     # [x]_B
print(c.T, (P.solve(A * x)).T)                     # Matrix([[-2/7, -2/7, 2/7]]) Matrix([[-4/7, -4/7, 18/7]])
                                                   # coordinates scaled by 2, 2, 9 independently

# The transpose operator on 2x2 matrices, written in the basis E11, E12, E21, E22
T = sp.Matrix([[1, 0, 0, 0],
               [0, 0, 1, 0],
               [0, 1, 0, 0],
               [0, 0, 0, 1]])                      # vec(M^T) = T vec(M)
for lam, mult, vs in T.eigenvects():
    print(lam, mult, [list(v.reshape(2, 2)) for v in vs])
# -1 1 [[0, -1, 1, 0]]
# 1 3 [[1, 0, 0, 0], [0, 1, 1, 0], [0, 0, 0, 1]]

In the last output each eigenvector is printed as a flattened $2 \times 2$ matrix: $\begin{bmatrix} 0 & -1 \ 1 & 0 \end{bmatrix}$ is skew-symmetric, and the three matrices for $\lambda = 1$ span the symmetric matrices.

Key takeaways

• An eigenbasis is built by concatenating bases of all the eigenspaces; independence across eigenspaces is automatic.
• An eigenbasis exists exactly when the eigenspace dimensions add up to $n$.
• In eigen-coordinates the operator multiplies each coordinate by its own eigenvalue, so the problem decouples.
• Eigenbases are not unique: vectors can be rescaled and reordered, and any basis of a multi-dimensional eigenspace works.
• The construction applies to operators on polynomials, matrices, and functions, not only to numerical matrices.

🧪 Practice

1. Find an eigenbasis of $\mathbb{R}^3$ for $\begin{bmatrix} 2 & 0 & 0 \ 0 & 1 & 1 \ 0 & 1 & 1 \end{bmatrix}$ and give the coordinates of $(1, 2, 3)$ in it.
2. Find an eigenbasis for $T(p) = p(2x)$ on $\mathbb{P}_2$ and give the eigenvalues. Then do the same for $S(p) = p(-x)$.
3. Show that $D = \tfrac{d}{dx}$ on $\mathbb{P}_2$ has no eigenbasis, but that $D$ on $\operatorname{span}{e^x, e^{2x}, e^{3x}}$ does.
4. Interview: You need to apply the same $1000 \times 1000$ diagonalizable matrix to a vector a million times with different numbers of repetitions $k$. How would an eigenbasis help, and what would you precompute? Hint: precompute $P$, $P^{-1}$ (or a factorization), and the eigenvalues once; then $A^k\mathbf{x} = P,(\lambda^k \odot P^{-1}\mathbf{x})$, two matrix-vector products and $n$ scalar powers per query, provided $P$ is well conditioned.

Diagonalization Procedure

Why it matters. This topic collects the theory into a reliable recipe that ends with a verifiable factorization $A = PDP^{-1}$. Most errors in practice come from mismatched orderings, eigenvectors placed as rows instead of columns, or skipped checks, so the recipe also lists the checks.

The procedure.

1. Eigenvalues. Solve $p_A(\lambda) = 0$, listing each eigenvalue with its algebraic multiplicity.
2. Eigenspaces. For each eigenvalue, find a basis of $\operatorname{Nul}(A - \lambda I)$.
3. Count. If the eigenspace dimensions add up to less than $n$, stop: $A$ is not diagonalizable.
4. Assemble. Put the eigenvectors as the columns of $P$, and put the matching eigenvalues on the diagonal of $D$ in the same order.
5. Check $AP = PD$. This avoids computing $P^{-1}$ and catches ordering mistakes.

Worked example. For $A = \begin{bmatrix} 4 & -1 & 6 \ 2 & 1 & 6 \ 2 & -1 & 8 \end{bmatrix}$, Section 8.1 found $p_A(\lambda) = (\lambda - 2)^2(\lambda - 9)$, $E_2 = \operatorname{span}{(1, 2, 0), (-3, 0, 1)}$, and $E_9 = \operatorname{span}{(1, 1, 1)}$. The dimensions add to $2 + 1 = 3$, so

$$P = \begin{bmatrix} 1 & -3 & 1 \ 2 & 0 & 1 \ 0 & 1 & 1 \end{bmatrix}, \qquad D = \begin{bmatrix} 2 & 0 & 0 \ 0 & 2 & 0 \ 0 & 0 & 9 \end{bmatrix}.$$

Check, using the first column: $A(1, 2, 0) = (4 - 2,\ 2 + 2,\ 2 - 2) = (2, 4, 0) = 2(1, 2, 0)$, which is the first column of $PD$. The other columns are checked the same way.

Reading $A = PDP^{-1}$ from right to left.

                      A
      x  ---------------------------->  A x
      |                                  ^
      | P^{-1}: standard -> eigen        | P: eigen -> standard
      v                                  |
   [x]_B  ---------------------------->  [Ax]_B
                      D
              (scale coordinate i by lambda_i)

To apply $A$: express $\mathbf{x}$ in eigen-coordinates, scale each coordinate by its eigenvalue, and convert back. This is the same picture as Section 7.2, with the special choice of basis that makes the middle matrix diagonal.

Common mistakes.

Mistake	Symptom	Fix
eigenvalues in $D$ in a different order from the columns of $P$	$AP \ne PD$	reorder one to match the other
eigenvectors written as rows of $P$	$AP \ne PD$	eigenvectors are columns
writing $A = P^{-1}DP$	correct only if $P$ is redefined	with eigenvectors in the columns of $P$, the formula is $A = PDP^{-1}$
a repeated eigenvalue listed once in $D$	$P$ is not square	list it $\operatorname{gm}$ times, one per basis vector
not checking $\operatorname{gm} = \operatorname{am}$	$P$ is singular	count eigenvectors before assembling

Non-uniqueness. Any reordering of the eigenvalues (with the same reordering of the columns), any nonzero rescaling of columns, and any change of basis inside an eigenspace give another valid pair $(P, D)$. $D$ is unique up to order; $P$ is far from unique.

Numerical diagonalization and conditioning. np.linalg.eig(A) returns the eigenvalues and a matrix of unit-length eigenvectors, which is a valid $P$ when $A$ is diagonalizable. How useful the factorization is depends on the condition number of $P$ (Chapter 11). If $P$ is nearly singular, which happens near defective matrices, then $P^{-1}$ amplifies rounding errors and $PDP^{-1}$ can be inaccurate. Symmetric matrices are the ideal case: $P$ can be chosen orthogonal, with condition number 1 (Chapter 10).

import numpy as np
import sympy as sp

A = sp.Matrix([[4, -1, 6], [2, 1, 6], [2, -1, 8]])
P = sp.Matrix([[1, -3, 1], [2, 0, 1], [0, 1, 1]])
D = sp.diag(2, 2, 9)
print(A * P == P * D)                      # True: the check, without inverting P
print(A == P * D * P.inv())                # True

# The wrong order is caught immediately
print(A * P == P * sp.diag(9, 2, 2))       # False

# sympy finds P and D itself (here the same ones; in general any valid choice)
P2, D2 = A.diagonalize()
print(P2.tolist(), D2.tolist())            # [[1, -3, 1], [2, 0, 1], [0, 1, 1]] [[2, 0, 0], [0, 2, 0], [0, 0, 9]]

# numpy: floating point, with a conditioning check on P
An = np.array(A.tolist(), dtype=float)
vals, V = np.linalg.eig(An)
print(np.allclose(An, V @ np.diag(vals) @ np.linalg.inv(V)))   # True
print(round(np.linalg.cond(V), 2))                              # 3.49: modest, so the factorization is trustworthy

# Near a defective matrix the eigenvector matrix is nearly singular
N = np.array([[1.0, 1.0], [1e-12, 1.0]])
_, V = np.linalg.eig(N)
print(f"{np.linalg.cond(V):.1e}")                               # 1.0e+06: columns almost parallel

Key takeaways

• Recipe: eigenvalues, eigenspace bases, count to $n$, eigenvectors as columns of $P$, eigenvalues in matching order in $D$.
• Verify with $AP = PD$, which needs no inverse and catches ordering errors.
• $A = PDP^{-1}$ means: convert to eigen-coordinates, scale each coordinate, convert back.
• $D$ is unique up to order; $P$ can be rescaled, reordered, and rebased within eigenspaces.
• Numerically, the condition number of $P$ measures how trustworthy the diagonalization is.

🧪 Practice

1. Diagonalize $\begin{bmatrix} 1 & 2 \ 2 & 1 \end{bmatrix}$ and $\begin{bmatrix} 2 & 2 & 1 \ 1 & 3 & 1 \ 1 & 2 & 2 \end{bmatrix}$ (the second has eigenvalues 1, 1, 5). Verify $AP = PD$.
2. Construct a $3 \times 3$ matrix with eigenvalues $1, 2, 3$ and eigenvectors $(1, 0, 0)$, $(1, 1, 0)$, $(1, 1, 1)$.
3. Suppose $A = PDP^{-1}$. Find a diagonalization of $A^T$ and of $A + 3I$ in terms of $P$ and $D$.
4. Interview: Why are symmetric matrices the "nice case" for diagonalization in numerical work? Hint: their eigenvalues are real, eigenvectors for distinct eigenvalues are orthogonal, and $P$ can be chosen orthogonal, so $P^{-1} = P^T$ costs nothing and has condition number 1 (Chapter 10).

Powers of Diagonalizable Matrices

Why it matters. Many processes apply the same linear step repeatedly: $\mathbf{x}_{k+1} = A\mathbf{x}_k$ (Section 4.1). After $k$ steps the state is $A^k\mathbf{x}_0$, and computing $A^k$ by repeated multiplication says little about what happens as $k$ grows. Diagonalization turns $A^k$ into $k$-th powers of numbers, which gives closed-form formulas and long-run predictions at once.

The formula. If $A = PDP^{-1}$, then the inner factors cancel in pairs:

$$A^k = (PDP^{-1})(PDP^{-1})\cdots(PDP^{-1}) = PD^kP^{-1}, \qquad D^k = \operatorname{diag}(\lambda_1^k, \dots, \lambda_n^k).$$

(Formally, by induction on $k$: $A^{k+1} = A^kA = PD^kP^{-1}PDP^{-1} = PD^{k+1}P^{-1}$.)

The eigenvector expansion. Writing $\mathbf{x}_0 = c_1\mathbf{v}_1 + \cdots + c_n\mathbf{v}_n$ in the eigenbasis,

$$A^k\mathbf{x}_0 = c_1\lambda_1^k\mathbf{v}_1 + c_2\lambda_2^k\mathbf{v}_2 + \cdots + c_n\lambda_n^k\mathbf{v}_n.$$

This one line explains long-run behavior. If $\lvert\lambda_1\rvert > \lvert\lambda_i\rvert$ for all $i \ge 2$ (a dominant eigenvalue) and $c_1 \ne 0$, the first term eventually outweighs all others:

$$A^k\mathbf{x}_0 \approx c_1\lambda_1^k\mathbf{v}_1 \quad \text{for large } k.$$

The state lines up with the dominant eigenvector, and its size changes by a factor of about $\lambda_1$ per step. The ratio $\lvert\lambda_2\rvert/\lvert\lambda_1\rvert$ sets the speed of convergence. This is the basis of the power method (Section 12.3), Markov chain steady states, and PageRank (Section 8.4).

Largest $\lvert\lambda_i\rvert$ (the spectral radius $\rho(A)$)	Behavior of $A^k\mathbf{x}_0$
$\rho(A) < 1$	decays to $\mathbf{0}$ for every $\mathbf{x}_0$
$\rho(A) = 1$	bounded if diagonalizable; may settle, cycle, or oscillate
$\rho(A) > 1$	grows without bound for almost every $\mathbf{x}_0$

Example: a closed formula for Fibonacci numbers. With $F_0 = 0$, $F_1 = 1$, $F_{k+1} = F_k + F_{k-1}$, the step is (Section 4.1)

$$\begin{bmatrix} F_{k+1} \ F_k \end{bmatrix} = \begin{bmatrix} 1 & 1 \ 1 & 0 \end{bmatrix}\begin{bmatrix} F_k \ F_{k-1} \end{bmatrix}.$$

The matrix has $p(\lambda) = \lambda^2 - \lambda - 1$, eigenvalues $\varphi = \tfrac{1 + \sqrt5}{2} \approx 1.618$ and $\psi = \tfrac{1 - \sqrt5}{2} \approx -0.618$, and eigenvectors $(\varphi, 1)$ and $(\psi, 1)$. Expanding $(F_1, F_0) = (1, 0) = \tfrac{1}{\sqrt5}\big[(\varphi, 1) - (\psi, 1)\big]$ and applying the formula gives Binet's formula:

$$F_k = \frac{\varphi^k - \psi^k}{\sqrt5}.$$

Since $\lvert\psi\rvert < 1$, the second term vanishes quickly, so $F_k$ is the integer nearest to $\varphi^k/\sqrt5$, and $F_{k+1}/F_k \to \varphi$.

Other functions of $A$. The same cancellation works for any function defined through powers:

Expression	Diagonal form
polynomial $q(A)$	$P,q(D),P^{-1}$, with $q(D) = \operatorname{diag}(q(\lambda_i))$
inverse $A^{-1}$ (all $\lambda_i \ne 0$)	$PD^{-1}P^{-1}$
square root (all $\lambda_i \ge 0$)	$PD^{1/2}P^{-1}$ squares to $A$
limit $\lim_k A^k$	$P,\big(\lim_k D^k\big),P^{-1}$, when the limit exists
exponential $e^{A}$	$Pe^{D}P^{-1}$ (next topic)

Two practical warnings. First, for exact integer results, such as Fibonacci numbers modulo a prime, exponentiation by squaring (Section 4.1) computes $A^k$ with about $2\log_2k$ matrix products and no irrational numbers. Second, eigenvalue-based formulas lose accuracy when $P$ is ill conditioned or eigenvalues of very different sizes are mixed.

import numpy as np

F = np.array([[1.0, 1.0], [1.0, 0.0]])
vals, P = np.linalg.eig(F)
vals, P = vals.real, P.real
print(np.round(vals, 6))                               # [ 1.618034 -0.618034]

def fib_eigen(k):
    Dk = np.diag(vals ** k)                            # powers of the eigenvalues only
    return (P @ Dk @ np.linalg.inv(P))[1, 0]           # A^k = [[F(k+1), F(k)], [F(k), F(k-1)]]

def fib_squaring(k):
    """Exact: exponentiation by squaring with Python integers (about log2 k products)."""
    def mul(X, Y):
        return [[X[0][0]*Y[0][0] + X[0][1]*Y[1][0], X[0][0]*Y[0][1] + X[0][1]*Y[1][1]],
                [X[1][0]*Y[0][0] + X[1][1]*Y[1][0], X[1][0]*Y[0][1] + X[1][1]*Y[1][1]]]
    R, B = [[1, 0], [0, 1]], [[1, 1], [1, 0]]
    while k:
        if k & 1:
            R = mul(R, B)
        B, k = mul(B, B), k >> 1
    return R[0][1]

phi = (1 + 5 ** 0.5) / 2
for k in [10, 30, 80]:
    print(k, round(fib_eigen(k)), round(phi ** k / 5 ** 0.5), fib_squaring(k))
# 10 55 55 55
# 30 832040 832040 832040
# 80 23416728348467740 23416728348467744 23416728348467685   (both floating-point results are wrong)

# Long run: the ratio of consecutive terms approaches the dominant eigenvalue
x = np.array([1.0, 0.0])
for _ in range(25):
    x = F @ x
print(round(x[0] / x[1], 10))                         # 1.6180339887

# Matrix square root through the eigen-decomposition
A = np.array([[5.0, 4.0], [1.0, 2.0]])                 # eigenvalues 6 and 1
w, V = np.linalg.eig(A)
root = V @ np.diag(np.sqrt(w.real)) @ np.linalg.inv(V)
print(np.allclose(root @ root, A))                     # True

Key takeaways

• $A^k = PD^kP^{-1}$: powers of a diagonalizable matrix reduce to powers of its eigenvalues.
• $A^k\mathbf{x}_0 = \sum c_i\lambda_i^k\mathbf{v}_i$; a dominant eigenvalue makes the state align with its eigenvector.
• The spectral radius decides decay ($\rho < 1$), boundedness ($\rho = 1$, diagonalizable), or growth ($\rho > 1$); $\lvert\lambda_2\rvert/\lvert\lambda_1\rvert$ sets the convergence rate.
• Binet's formula for Fibonacci numbers is the eigenvector expansion of a $2 \times 2$ recurrence.
• Polynomials, inverses, roots, and limits of $A$ are computed through $D$; for exact integer powers, use exponentiation by squaring.

🧪 Practice

1. Compute $A^{10}$ for $A = \begin{bmatrix} 4 & -2 \ 1 & 1 \end{bmatrix}$ by diagonalization, and find a formula for $A^k$.
2. Let $A = \begin{bmatrix} 0.9 & 0.2 \ 0.1 & 0.8 \end{bmatrix}$. Find $\lim_{k \to \infty} A^k$ and explain what its columns have in common.
3. Derive a closed form for the sequence $a_{k+1} = 3a_k - 2a_{k-1}$ with $a_0 = 0$, $a_1 = 1$ by diagonalizing a $2 \times 2$ matrix.
4. Interview: Explain why the ratio $F_{k+1}/F_k$ converges to the golden ratio, and how fast, using only eigenvalues. Hint: the error is governed by $(\psi/\varphi)^k \approx (-0.382)^k$, so each step gains about $\log_{10}(1/0.382) \approx 0.42$ correct decimal digits, and the sign alternates.

Matrix Exponential

Why it matters. The scalar equation $x'(t) = ax(t)$ has solution $x(t) = e^{at}x(0)$. Systems of linear differential equations, $\mathbf{x}'(t) = A\mathbf{x}(t)$, describe coupled springs, circuits, chemical reactions, epidemic models, and continuous-time Markov chains. The matrix exponential $e^{tA}$ gives the same one-line solution, $\mathbf{x}(t) = e^{tA}\mathbf{x}(0)$, and diagonalization is the easiest way to compute it.

Definition. For a square matrix $A$,

$$e^{A} = \sum_{k=0}^{\infty}\frac{A^k}{k!} = I + A + \frac{A^2}{2!} + \frac{A^3}{3!} + \cdots.$$

The series converges for every matrix, because $\lVert A^k\rVert/k! \le \lVert A\rVert^k/k!$ and $\sum\lVert A\rVert^k/k! = e^{\lVert A\rVert}$ is finite. It is not the entrywise exponential: $e^{A}$ applied entry by entry is a different, rarely useful matrix.

Diagonalizable case. If $A = PDP^{-1}$, every term of the series is $P\frac{D^k}{k!}P^{-1}$, so

$$e^{tA} = P,e^{tD},P^{-1}, \qquad e^{tD} = \operatorname{diag}(e^{t\lambda_1}, \dots, e^{t\lambda_n}).$$

Expanding $\mathbf{x}(0) = \sum c_i\mathbf{v}_i$ gives the solution of $\mathbf{x}' = A\mathbf{x}$ as a sum of independent exponential modes (Section 8.4):

$$\mathbf{x}(t) = c_1e^{\lambda_1t}\mathbf{v}_1 + \cdots + c_ne^{\lambda_nt}\mathbf{v}_n.$$

The real parts of the eigenvalues decide stability in continuous time: if every $\operatorname{Re}\lambda_i < 0$, all solutions decay to $\mathbf{0}$. Compare with the discrete case, where the test was $\lvert\lambda_i\rvert < 1$; the two are linked by $\lvert e^{\lambda}\rvert = e^{\operatorname{Re}\lambda}$.

Properties.

Property	Comment
$e^{0} = I$	the series starts with $I$
$\frac{d}{dt}e^{tA} = Ae^{tA} = e^{tA}A$	differentiate the series term by term; this is why $e^{tA}\mathbf{x}_0$ solves $\mathbf{x}' = A\mathbf{x}$
$e^{(s + t)A} = e^{sA}e^{tA}$	$sA$ and $tA$ commute
$(e^{A})^{-1} = e^{-A}$	so $e^{A}$ is always invertible
$\det e^{A} = e^{\operatorname{tr}A}$	over $\mathbb{C}$: $\prod e^{\lambda_i} = e^{\sum\lambda_i}$
$e^{A + B} = e^{A}e^{B}$ if $AB = BA$	false in general; see the code below
$e^{P^{-1}AP} = P^{-1}e^{A}P$	the exponential respects similarity
eigenvalues of $e^{A}$ are $e^{\lambda_i}$	with the same eigenvectors

Two examples where the series can be summed directly.

• Nilpotent: for $N = \begin{bmatrix} 0 & 1 \ 0 & 0 \end{bmatrix}$, $N^2 = 0$, so the series stops: $e^{tN} = I + tN = \begin{bmatrix} 1 & t \ 0 & 1 \end{bmatrix}$. Polynomials in $t$ (not just exponentials) appear for defective matrices, which is why $te^{\lambda t}$ shows up for repeated roots of differential equations. Since $\lambda I$ commutes with $N$, $e^{t(\lambda I + N)} = e^{\lambda t}(I + tN)$.
• Rotation generator: for $J = \begin{bmatrix} 0 & -1 \ 1 & 0 \end{bmatrix}$, $J^2 = -I$, so the even and odd terms of the series regroup into the cosine and sine series: $e^{tJ} = \cos t,I + \sin t,J = R_t$, the rotation by angle $t$. This is the matrix version of Euler's formula $e^{it} = \cos t + i\sin t$, with $J$ playing the role of $i$ (Section 1.2).

Computing it in practice. Use scipy.linalg.expm, which uses scaling and squaring with Padé approximants: it computes $e^{A/2^s}$ for a small matrix $A/2^s$ and squares the result $s$ times. The eigen-decomposition formula is fine for symmetric matrices and for well-conditioned $P$, but fails for defective or nearly defective matrices, and truncating the Taylor series directly is inaccurate when $\lVert A\rVert$ is large.

import numpy as np
from math import factorial
from scipy.linalg import expm

A = np.array([[-2.0, 1.0],
              [ 1.0, -2.0]])                       # eigenvalues -1 and -3: a decaying system

w, V = np.linalg.eig(A)
via_eig = V @ np.diag(np.exp(w)) @ np.linalg.inv(V)
via_series = sum(np.linalg.matrix_power(A, k) / factorial(k) for k in range(30))
print(np.allclose(expm(A), via_eig), np.allclose(expm(A), via_series))  # True True
print(np.allclose(expm(A), np.exp(A)))             # False: expm is not the entrywise exp

# Solve x' = A x with x(0) = (1, 0): x(t) = expm(t A) x0
x0 = np.array([1.0, 0.0])
for t in [0.0, 1.0, 3.0]:
    print(t, np.round(expm(t * A) @ x0, 4))
# 0.0 [1. 0.]
# 1.0 [0.2088 0.159 ]
# 3.0 [0.025  0.0248]   the slow mode e^{-t} (1, 1)/2 dominates

# det(e^A) = e^{tr A}
print(np.isclose(np.linalg.det(expm(A)), np.exp(np.trace(A))))   # True

# Rotation generator and nilpotent shear
J = np.array([[0.0, -1.0], [1.0, 0.0]])
print(np.round(expm(np.pi / 2 * J), 10))           # [[ 0. -1.]  [ 1.  0.]]: rotation by 90 degrees
N = np.array([[0.0, 1.0], [0.0, 0.0]])
print(expm(3 * N))                                 # [[1. 3.]  [0. 1.]]: I + 3N, the series stops

# e^{A+B} != e^A e^B when AB != BA
B = np.array([[0.0, 0.0], [1.0, 0.0]])
print(np.allclose(expm(N + B), expm(N) @ expm(B)))  # False

Key takeaways

• $e^{A} = \sum A^k/k!$ converges for every square matrix and is not the entrywise exponential.
• $\mathbf{x}(t) = e^{tA}\mathbf{x}(0)$ solves $\mathbf{x}' = A\mathbf{x}$; for diagonalizable $A$, $e^{tA} = Pe^{tD}P^{-1}$ with entries $e^{\lambda_it}$.
• Continuous-time stability is decided by the signs of $\operatorname{Re}\lambda_i$, discrete-time stability by $\lvert\lambda_i\rvert$ versus 1.
• $e^{A}$ is always invertible, $\det e^{A} = e^{\operatorname{tr}A}$, and $e^{A + B} = e^{A}e^{B}$ only when $A$ and $B$ commute.
• Nilpotent parts produce polynomial factors such as $te^{\lambda t}$; use scipy.linalg.expm in practice.

🧪 Practice

1. Compute $e^{tA}$ for $A = \begin{bmatrix} 1 & 0 \ 0 & -2 \end{bmatrix}$, $A = \begin{bmatrix} 0 & 1 \ 1 & 0 \end{bmatrix}$ (diagonalize first), and $A = \begin{bmatrix} 3 & 1 \ 0 & 3 \end{bmatrix}$ (split as $3I + N$).
2. Show from the series that $e^{tJ} = \cos t,I + \sin t,J$ when $J^2 = -I$.
3. Find $2 \times 2$ matrices $A$, $B$ with $e^{A + B} \ne e^{A}e^{B}$, and verify the inequality by hand using nilpotent matrices.
4. Interview: A continuous-time Markov chain has rate matrix $Q$ (rows sum to 0). Why is $e^{tQ}$ a transition matrix for every $t \ge 0$, and why would you call expm rather than diagonalize $Q$? Hint: $Q\mathbf{1} = \mathbf{0}$ gives $e^{tQ}\mathbf{1} = \mathbf{1}$ (rows sum to 1) and nonnegativity follows from $Q + cI \ge 0$; rate matrices are often non-symmetric and nearly defective, so the eigenvector matrix can be badly conditioned.

<a id="83-advanced-eigenstructure"></a>

8.3 Advanced Eigenstructure

Not every matrix is diagonalizable. This subchapter develops what can always be done: triangularize with the Schur decomposition, describe which polynomials annihilate a matrix, and reach the Jordan and rational canonical forms, which decide exactly when two matrices are similar.

Invariant Subspaces

Why it matters. Diagonalization splits a space into one-dimensional pieces that the operator maps into themselves. When that is impossible, the next best thing is to split the space into larger pieces that are still mapped into themselves. Such pieces let a big problem be broken into smaller independent ones, and they are the building blocks of every canonical form in this subchapter. They also power the Krylov subspace methods (Arnoldi, Lanczos, GMRES) that solve the largest linear systems and eigenproblems in practice (Chapter 12).

Analogy. Think of a building where $T$ is "take one step". A floor is an invariant subspace if no single step can take you off it: you can walk around on it forever, but never leave. An eigen-line is a floor so narrow that you can only move along one corridor, toward or away from the origin.

Definition. A subspace $U \subseteq V$ is invariant under $T : V \to V$ (or $T$-invariant) if $T(U) \subseteq U$, that is, $T\mathbf{u} \in U$ for every $\mathbf{u} \in U$. The restriction $T|_U : U \to U$ is then a linear operator on $U$.

Examples. For every operator $T$:

Subspace	Why it is invariant
${\mathbf{0}}$ and $V$	trivially
$\ker T$	$T\mathbf{u} = \mathbf{0} \in \ker T$
$\operatorname{im} T$	$T(T\mathbf{v}) \in \operatorname{im} T$
an eigenspace $E_\lambda$, or any span of eigenvectors	$T\mathbf{v} = \lambda\mathbf{v}$ (Section 8.1)
$\ker q(T)$ and $\operatorname{im} q(T)$ for any polynomial $q$	$T$ commutes with $q(T)$
the Krylov subspace $\mathcal{K}(\mathbf{x}) = \operatorname{span}{\mathbf{x}, T\mathbf{x}, T^2\mathbf{x}, \dots}$	applying $T$ shifts the list by one; it is the smallest invariant subspace containing $\mathbf{x}$

A one-dimensional invariant subspace is exactly a line spanned by an eigenvector: $T\mathbf{v} \in \operatorname{span}{\mathbf{v}}$ means $T\mathbf{v} = \lambda\mathbf{v}$. Invariant subspaces are therefore the higher-dimensional generalization of eigenvectors.

Invariant subspaces create zero blocks. Choose a basis $\mathbf{u}_1, \dots, \mathbf{u}_k$ of an invariant $U$ and extend it to a basis of $V$ (Section 6.4). Since $T\mathbf{u}_j \in U$, the first $k$ columns of $[T]$ have zeros below row $k$:

$$[T] = \begin{bmatrix} A_{11} & A_{12} \ 0 & A_{22} \end{bmatrix}, \qquad A_{11} = [T|U], \qquad p_T = p{A_{11}},p_{A_{22}}.$$

$A_{22}$ is the matrix of the operator that $T$ induces on the quotient space $V/U$ (Section 7.5). If $V = U \oplus W$ with both $U$ and $W$ invariant, choosing a basis of $W$ for the remaining vectors also makes $A_{12} = 0$, so $[T]$ is block diagonal and $T$ splits into two independent operators. Diagonalization is the extreme case in which all the blocks are $1 \times 1$.

   U invariant:                    U and W both invariant (V = U + W, direct):

   [ T|_U   *   ]                  [ T|_U    0   ]
   [   0  T_V/U ]                  [   0   T|_W  ]

   block upper triangular          block diagonal: two independent problems

Real operators always have small invariant subspaces. A real matrix may have no real eigenvector, but every real operator on $\mathbb{R}^n$ ($n \ge 1$) has an invariant subspace of dimension 1 or 2. If $\lambda = a - bi$ ($b \ne 0$) is a complex eigenvalue with eigenvector $\mathbf{v} = \mathbf{x} + i\mathbf{y}$, then $A\mathbf{x} = a\mathbf{x} + b\mathbf{y}$ and $A\mathbf{y} = -b\mathbf{x} + a\mathbf{y}$ (Section 8.1), so $\operatorname{span}{\mathbf{x}, \mathbf{y}}$ is a real invariant plane on which $A$ acts as a rotation-scaling.

Testing invariance. If the columns of $U$ form a basis of the subspace, then the subspace is invariant exactly when every column of $AU$ lies in $\operatorname{Col} U$, that is, when $\operatorname{rank}[,U\ \ AU,] = \operatorname{rank}U$. Equivalently, $AU = UB$ for some small matrix $B$, which is then the matrix of the restriction.

import numpy as np

A = np.array([[2.0, 0.0, 1.0, -1.0],
              [0.0, 2.0, 1.0, -1.0],
              [3.0, -3.0, 5.0, -6.0],
              [3.0, -3.0, 3.0, -4.0]])

def is_invariant(A, U):
    """Col U is A-invariant iff appending A U does not increase the rank."""
    return np.linalg.matrix_rank(np.hstack([U, A @ U])) == np.linalg.matrix_rank(U)

U = np.array([[1.0, 1.0, 0.0, 0.0],
              [0.0, 1.0, 1.0, 0.0]]).T              # columns (1,1,0,0) and (0,1,1,0)
print(is_invariant(A, U))                          # True
B, *_ = np.linalg.lstsq(U, A @ U, rcond=None)      # solve U B = A U: the matrix of A restricted to Col U
print(np.round(B, 10))                             # [[2. 1.]  [0. 2.]]: a Jordan block, see below
print(is_invariant(A, np.array([[1.0, 0.0, 0.0, 0.0]]).T))   # False: e1 is not an eigenvector

# Krylov subspaces: the smallest invariant subspace containing x
def krylov_dim(A, x):
    K = np.column_stack([np.linalg.matrix_power(A, j) @ x for j in range(A.shape[0])])
    return np.linalg.matrix_rank(K)

print(krylov_dim(A, np.array([1.0, 1.0, 0.0, 0.0])))   # 1: an eigenvector
print(krylov_dim(A, np.array([0.0, 1.0, 1.0, 0.0])))   # 2
print(krylov_dim(A, np.array([0.0, 0.0, 1.0, 0.0])))   # 3: the maximum; never 4 for this A (see Minimal Polynomial)

Key takeaways

• $U$ is invariant under $T$ when $T(U) \subseteq U$; one-dimensional invariant subspaces are eigen-lines.
• Kernels, images, eigenspaces, $\ker q(T)$, $\operatorname{im} q(T)$, and Krylov subspaces are always invariant.
• An invariant subspace gives a block upper triangular matrix; two complementary invariant subspaces give a block diagonal one.
• Every real operator has an invariant line or plane; complex eigenvalues give invariant planes $\operatorname{span}{\operatorname{Re}\mathbf{v}, \operatorname{Im}\mathbf{v}}$.
• Test invariance with $\operatorname{rank}[U\ AU] = \operatorname{rank}U$, and read off the restriction from $AU = UB$.

🧪 Practice

1. List all invariant subspaces of the shear $\begin{bmatrix} 1 & 1 \ 0 & 1 \end{bmatrix}$ on $\mathbb{R}^2$, and of the rotation by $90°$.
2. Show that $\mathbb{P}_k$ is an invariant subspace of $D = \tfrac{d}{dx}$ on $\mathbb{P}_n$ for each $k \le n$, and write the block form of $D$ on $\mathbb{P}_3$ adapted to $\mathbb{P}_1 \subseteq \mathbb{P}_3$.
3. Prove that if $AB = BA$, then every eigenspace of $A$ is invariant under $B$. Deduce that commuting complex matrices share at least one eigenvector.
4. Interview: Explain why Krylov methods such as GMRES or Lanczos can solve or approximate problems with a $10^6 \times 10^6$ sparse matrix using only matrix-vector products. Hint: they work inside $\operatorname{span}{\mathbf{b}, A\mathbf{b}, \dots, A^{m-1}\mathbf{b}}$ for $m \ll n$, which needs only $m$ products; if this subspace becomes invariant, the restricted problem is solved exactly.

Triangularization and Schur Decomposition

Why it matters. Diagonalization fails for defective matrices and is numerically fragile near them. Triangularization always succeeds over $\mathbb{C}$, and the Schur decomposition does it with a unitary (orthonormal) change of basis, which is perfectly conditioned. That makes Schur the form that numerical software actually computes: eig, expm, matrix square roots, and the solution of Sylvester equations all go through it. It is also the standard tool for proving facts about eigenvalues.

Theorem (Schur). For every complex $n \times n$ matrix $A$ there is a unitary $U$ ($U^*U = I$, Section 4.2) and an upper triangular $T$ with

$$A = UTU^*.$$

The diagonal entries of $T$ are the eigenvalues of $A$, in any order you choose.

Proof. By induction on $n$; $n = 1$ is trivial. Over $\mathbb{C}$, $A$ has an eigenvalue $\lambda_1$ with a unit eigenvector $\mathbf{u}_1$ (Section 8.1). Extend $\mathbf{u}_1$ to an orthonormal basis $\mathbf{u}_1, \dots, \mathbf{u}_n$ of $\mathbb{C}^n$ (Gram-Schmidt, Chapter 9) and let $U_1$ have these columns. The first column of $U_1^AU_1$ is $U_1^(\lambda_1\mathbf{u}_1) = \lambda_1\mathbf{e}_1$, so

$$U_1^*AU_1 = \begin{bmatrix} \lambda_1 & \ast \ \mathbf{0} & A_1 \end{bmatrix}$$

with $A_1$ of size $n - 1$. By induction $A_1 = U_2T_1U_2^*$; then $U = U_1\begin{bmatrix} 1 & 0 \ 0 & U_2 \end{bmatrix}$ is unitary and $U^*AU$ is upper triangular. The diagonal of a triangular matrix holds its eigenvalues, and similar matrices share eigenvalues (Section 8.2). $\blacksquare$

In terms of invariant subspaces: the first $k$ columns of $U$ span an invariant subspace for every $k$, giving a nested chain ${0} \subset U_1 \subset U_2 \subset \cdots \subset \mathbb{C}^n$ (a flag), each one dimension larger and each mapped into itself.

Variants.

Form	Statement	When it exists
triangularization	$A = PTP^{-1}$, $P$ invertible, $T$ upper triangular	whenever $p_A$ splits over $F$ (always over $\mathbb{C}$)
complex Schur	$A = UTU^*$, $U$ unitary	every complex square matrix
real Schur	$A = QTQ^T$, $Q$ real orthogonal, $T$ block upper triangular with $1 \times 1$ and $2 \times 2$ diagonal blocks	every real square matrix; each $2 \times 2$ block carries a complex-conjugate pair
spectral theorem	$A = UDU^*$, $D$ diagonal	exactly when $A$ is normal, $AA^* = A^*A$ (Chapter 10)

Consequences proved in one line each.

• $\operatorname{tr}A = \operatorname{tr}T = \sum\lambda_i$ and $\det A = \det T = \prod\lambda_i$ for every complex matrix, diagonalizable or not.
• The eigenvalues of $A^k$ are $\lambda_i^k$ and those of $q(A)$ are $q(\lambda_i)$, with multiplicities, since $T^k$ and $q(T)$ are triangular with those diagonals.
• A normal matrix has a diagonal Schur form: a normal triangular matrix must be diagonal (compare the squared lengths of the first row and first column of $T^T = TT^$, then recurse).
• The Cayley-Hamilton theorem (next topic).

Why software prefers Schur to diagonalization. In $A = PDP^{-1}$, errors are amplified by the condition number of $P$, which can be huge. In $A = UTU^*$ the change of basis is unitary with condition number exactly 1, so rounding errors are not amplified by the change of basis. The price is a triangular rather than diagonal $T$. The QR algorithm (Section 12.3) computes the Schur form iteratively, and eigenvectors, if needed, are recovered afterwards from $T$ by back substitution.

import numpy as np
from scipy.linalg import schur

A = np.array([[2.0, 0.0, 1.0, -1.0],
              [0.0, 2.0, 1.0, -1.0],
              [3.0, -3.0, 5.0, -6.0],
              [3.0, -3.0, 3.0, -4.0]])            # defective: not diagonalizable

T, U = schur(A, output='complex')                  # A = U T U^*
print(np.allclose(A, U @ T @ U.conj().T))          # True
print(np.allclose(np.tril(T, -1), 0))              # True: T is upper triangular
print(np.round(np.diag(T).real, 6))                # [ 2. -1.  2.  2.]: eigenvalues, in an order chosen by LAPACK
print(round(np.linalg.cond(U), 12))                # 1.0: unitary, perfectly conditioned

# trace and determinant from the diagonal, even for a defective matrix
d = np.diag(T)
print(np.isclose(d.sum(), np.trace(A)), np.isclose(np.prod(d), np.linalg.det(A)))   # True True

# Real Schur form: 2x2 diagonal blocks for complex-conjugate pairs
R = np.array([[1.0, -2.0, 0.0],
              [1.0,  3.0, 0.0],
              [0.0,  1.0, 5.0]])                    # eigenvalues 2 + i, 2 - i, 5
Tr, Q = schur(R)                                   # real output by default
print(np.round(Tr, 4))                             # a 2x2 block (eigenvalues 2 +- i) and a 1x1 block (5)
print(np.allclose(Q @ Q.T, np.eye(3)), np.allclose(R, Q @ Tr @ Q.T))   # True True

The order of the eigenvalues along the diagonal and the entries above it depend on the library and its version. What is guaranteed is the structure checked by the True lines.

Key takeaways

• Every complex square matrix is unitarily similar to an upper triangular matrix: $A = UTU^*$, with the eigenvalues on the diagonal of $T$.
• Over $\mathbb{R}$, the real Schur form uses an orthogonal $Q$ and $2 \times 2$ blocks for complex pairs.
• The columns of $U$ give a nested chain of invariant subspaces.
• Schur proves $\operatorname{tr}A = \sum\lambda_i$, $\det A = \prod\lambda_i$, and the eigenvalue rules for $A^k$ and $q(A)$ for all matrices.
• Software computes Schur forms rather than diagonalizations because unitary changes of basis do not amplify errors.

🧪 Practice

1. Find a Schur decomposition of $\begin{bmatrix} 1 & 1 \ -1 & 3 \end{bmatrix}$ by hand: take a unit eigenvector as the first column of $U$ and a perpendicular unit vector as the second.
2. Use the Schur form to prove that if $A$ is nilpotent then $\operatorname{tr}A^k = 0$ for every $k \ge 1$.
3. Show that if $T$ is upper triangular and normal ($T^T = TT^$), then $T$ is diagonal. Start with the $(1, 1)$ entries of both sides.
4. Interview: Why does scipy.linalg.expm or a matrix square root routine work through the Schur form instead of an eigendecomposition? Hint: Schur exists for every matrix and uses a unitary basis change (condition number 1); functions of a triangular matrix can be computed stably by recurrences (Parlett's method), while $PDP^{-1}$ fails or loses accuracy when $P$ is ill conditioned.

Cayley-Hamilton Theorem

Why it matters. The characteristic polynomial was built to find eigenvalues, but it says much more: every matrix is a root of its own characteristic polynomial. As a result, every power $A^k$ and every inverse $A^{-1}$ can be written as a polynomial in $A$ of degree less than $n$. This underlies efficient formulas for powers, the theory of the minimal polynomial and of Krylov methods, and control theory's controllability tests.

Theorem (Cayley-Hamilton). If $p_A(\lambda) = \lambda^n + c_{n-1}\lambda^{n-1} + \cdots + c_1\lambda + c_0$ is the characteristic polynomial of the $n \times n$ matrix $A$, then

$$p_A(A) = A^n + c_{n-1}A^{n-1} + \cdots + c_1A + c_0I = 0.$$

The tempting false proof. "Substitute $\lambda = A$ into $p_A(\lambda) = \det(\lambda I - A)$ to get $\det(A - A) = 0$." This is wrong: $p_A(A)$ is a matrix, while $\det(\cdot)$ is a number; the substitution confuses the two. A correct proof needs more work.

Proof for diagonalizable $A$. $p_A(A) = P,p_A(D),P^{-1}$ and $p_A(D) = \operatorname{diag}(p_A(\lambda_1), \dots, p_A(\lambda_n)) = 0$, since every $\lambda_i$ is a root.

Proof in general. By Schur, $A = UTU^$ with $T$ upper triangular and diagonal $t_{11}, \dots, t_{nn}$, so $p_A(\lambda) = \prod_k(\lambda - t_{kk})$ and $p_A(A) = U,p_A(T),U^$. Let $W_k = \operatorname{span}{\mathbf{e}1, \dots, \mathbf{e}k}$. Column $k$ of $T - t{kk}I$ has zeros in position $k$ and below, so $(T - t{kk}I)$ maps $W_k$ into $W_{k-1}$. In

$$p_A(T) = (T - t_{11}I)(T - t_{22}I)\cdots(T - t_{nn}I),$$

the rightmost factor maps $W_n = \mathbb{C}^n$ into $W_{n-1}$, the next into $W_{n-2}$, and so on, until $(T - t_{11}I)$ maps $W_1$ to $W_0 = {\mathbf{0}}$. So $p_A(T) = 0$, and $p_A(A) = 0$. $\blacksquare$

The $2 \times 2$ case. Cayley-Hamilton reads $A^2 - (\operatorname{tr}A)A + (\det A)I = 0$. For $A = \begin{bmatrix} 1 & 2 \ 3 & 4 \end{bmatrix}$: $A^2 = \begin{bmatrix} 7 & 10 \ 15 & 22 \end{bmatrix}$, and $A^2 - 5A - 2I = \begin{bmatrix} 7 - 5 - 2 & 10 - 10 \ 15 - 15 & 22 - 20 - 2 \end{bmatrix} = 0$.

Applications.

Task	How Cayley-Hamilton helps
inverse as a polynomial	if $c_0 = (-1)^n\det A \ne 0$: $A^{-1} = -\tfrac{1}{c_0}\big(A^{n-1} + c_{n-1}A^{n-2} + \cdots + c_1I\big)$
reducing high powers	divide $\lambda^k = q(\lambda)p_A(\lambda) + r(\lambda)$ with $\deg r < n$; then $A^k = r(A)$
every function of $A$	$e^{A}$, $A^{1/2}$, ... all equal polynomials in $A$ of degree $< n$ (with matrix-dependent coefficients)
Krylov spaces	$A^n\mathbf{x}$ is a combination of $\mathbf{x}, \dots, A^{n-1}\mathbf{x}$, so $\mathcal{K}(\mathbf{x})$ has dimension $\le n$
control theory	reachability of $\mathbf{x}_{k+1} = A\mathbf{x}_k + B\mathbf{u}_k$ is decided by $[,B\ AB\ \cdots\ A^{n-1}B,]$, because higher powers add nothing

The inverse formula shows that $A^{-1}$, when it exists, commutes with $A$ and with everything that commutes with $A$. It is a theoretical tool, not an efficient algorithm: computing $A^{-1}$ by elimination (Section 4.3) is cheaper and more stable.

import sympy as sp

lam = sp.symbols('lambda')
A = sp.Matrix([[2, 0, 1, -1],
               [0, 2, 1, -1],
               [3, -3, 5, -6],
               [3, -3, 3, -4]])
n = A.shape[0]

p = A.charpoly(lam)                                  # PurePoly in lambda
coeffs = p.all_coeffs()                              # [1, c_{n-1}, ..., c_0]
print(sp.factor(p.as_expr()), coeffs)                # (lambda - 2)**3*(lambda + 1) [1, -5, 6, 4, -8]

def poly_at(cs, M):
    """Evaluate a polynomial (coefficients high -> low) at a matrix with Horner's rule."""
    R = sp.zeros(*M.shape)
    for c in cs:
        R = R * M + c * sp.eye(M.shape[0])
    return R

print(poly_at(coeffs, A) == sp.zeros(n, n))          # True: Cayley-Hamilton

# Inverse as a polynomial in A: A^{-1} = -(A^3 - 5A^2 + 6A + 4I) / c0 with c0 = -8
Ainv = -poly_at(coeffs[:-1], A) / coeffs[-1]
print(Ainv == A.inv())                               # True

# A^100 = r(A), where r is the remainder of lambda^100 divided by p (degree < 4)
r = sp.rem(sp.Poly(lam**100, lam), sp.Poly(p.as_expr(), lam))
print(r.degree(), poly_at(r.all_coeffs(), A) == A**100)   # 3 True

Key takeaways

• Every square matrix satisfies its characteristic polynomial: $p_A(A) = 0$.
• "Substitute $\lambda = A$ into $\det(\lambda I - A)$" is not a proof; the triangular (Schur) argument is.
• For $2 \times 2$ matrices: $A^2 = (\operatorname{tr}A)A - (\det A)I$.
• Every power, the inverse, and every analytic function of $A$ is a polynomial in $A$ of degree less than $n$.
• High powers reduce to $A^k = r(A)$ where $r$ is the remainder of $\lambda^k$ modulo $p_A$.

🧪 Practice

1. Verify Cayley-Hamilton for $\begin{bmatrix} 2 & 1 \ -1 & 0 \end{bmatrix}$, and use it to write $A^{-1}$ and $A^5$ as linear combinations of $I$ and $A$.
2. A $3 \times 3$ matrix has characteristic polynomial $\lambda^3 - 2\lambda^2 + \lambda - 3$. Express $A^{-1}$ and $A^4$ as polynomials in $A$ of degree at most 2.
3. Explain precisely where the "substitute $\lambda = A$" argument breaks down, and give an analogous false argument that $\operatorname{tr}(\lambda I - A)$ evaluated at $\lambda = A$ is zero, which is easily seen to fail.
4. Interview: You need $A^{10^{9}}\mathbf{x}$ for a $3 \times 3$ integer matrix modulo a prime. Compare exponentiation by squaring with the Cayley-Hamilton reduction. Hint: both cost $O(\log k)$ steps; squaring multiplies $3 \times 3$ matrices, the reduction computes $\lambda^k \bmod p_A(\lambda)$ by squaring polynomials of degree $< 3$, which is cheaper for large $n$ and is the idea behind fast linear-recurrence algorithms.

Minimal Polynomial

Why it matters. The characteristic polynomial annihilates $A$, but often a polynomial of smaller degree does too: the identity satisfies $A - I = 0$, not merely $(A - I)^n = 0$. The smallest such polynomial, the minimal polynomial, records exactly how far from diagonalizable a matrix is, gives a clean diagonalizability test, and determines the dimensions of Krylov subspaces and the number of steps some iterative solvers need.

Definition. The minimal polynomial $m_A$ of a square matrix $A$ is the monic polynomial of least degree with $m_A(A) = 0$.

It exists because Cayley-Hamilton gives at least one annihilating polynomial, of degree $n$. It is unique: if two monic polynomials of the same least degree annihilate $A$, their difference annihilates $A$ and has lower degree, so it must be zero.

Theorem. If $q(A) = 0$, then $m_A$ divides $q$. In particular $m_A$ divides $p_A$.

Proof. Divide with remainder: $q = s,m_A + r$ with $\deg r < \deg m_A$. Then $r(A) = q(A) - s(A)m_A(A) = 0$. By minimality $r$ cannot be a nonzero polynomial, so $r = 0$. $\blacksquare$

Theorem. $m_A$ and $p_A$ have the same roots: the eigenvalues of $A$.

Proof. Roots of $m_A$ are roots of $p_A$ because $m_A$ divides $p_A$. Conversely, if $A\mathbf{v} = \lambda\mathbf{v}$ with $\mathbf{v} \ne \mathbf{0}$, then $\mathbf{0} = m_A(A)\mathbf{v} = m_A(\lambda)\mathbf{v}$, so $m_A(\lambda) = 0$. $\blacksquare$

So, over $\mathbb{C}$, if $p_A(t) = \prod(t - \lambda_i)^{a_i}$, then $m_A(t) = \prod(t - \lambda_i)^{s_i}$ with $1 \le s_i \le a_i$. The exponent $s_i$ turns out to be the size of the largest Jordan block for $\lambda_i$ (Jordan Canonical Form below).

Matrix	$p_A(t)$	$m_A(t)$
$I_n$	$(t - 1)^n$	$t - 1$
$\begin{bmatrix} 1 & 1 \ 0 & 1 \end{bmatrix}$	$(t - 1)^2$	$(t - 1)^2$
$\operatorname{diag}(2, 2, 3)$	$(t - 2)^2(t - 3)$	$(t - 2)(t - 3)$
$3 \times 3$ Jordan block $J_3(\lambda)$	$(t - \lambda)^3$	$(t - \lambda)^3$
a projection $P \ne 0, I$	$t^r(t - 1)^{n - r}$	$t(t - 1)$
the $4 \times 4$ running example $A$ of this subchapter	$(t - 2)^3(t + 1)$	$(t - 2)^2(t + 1)$

Diagonalizability criterion. $A$ is diagonalizable over $F$ if and only if $m_A$ splits over $F$ into distinct linear factors, $m_A(t) = (t - \lambda_1)\cdots(t - \lambda_k)$.

Proof. ($\Rightarrow$) If $A$ has an eigenbasis, $(A - \lambda_1I)\cdots(A - \lambda_kI)$ kills every eigenvector (the factors commute, so the one for its own eigenvalue can be applied first), hence kills a basis, hence is $0$. ($\Leftarrow$) Let $\ell_i(t) = \prod_{j \ne i}\frac{t - \lambda_j}{\lambda_i - \lambda_j}$ be the Lagrange polynomials for the distinct $\lambda_i$. Then $\sum_i\ell_i(t) = 1$ (a polynomial of degree $< k$ equal to 1 at $k$ points), so every $\mathbf{x} = \sum_i\ell_i(A)\mathbf{x}$. Each $(A - \lambda_iI)\ell_i(A)$ is a multiple of $m_A(A) = 0$, so $\ell_i(A)\mathbf{x} \in E_{\lambda_i}$. Thus eigenvectors span $F^n$. $\blacksquare$

Immediate consequences.

• Projections ($P^2 = P$, annihilated by $t(t - 1)$) are diagonalizable.
• Involutions ($A^2 = I$, annihilated by $(t - 1)(t + 1)$) are diagonalizable when $1 \ne -1$ in $F$, that is, outside characteristic 2.
• If $A^k = I$ over $\mathbb{C}$, then $A$ is diagonalizable, since $t^k - 1$ has $k$ distinct roots (the $k$-th roots of unity). Permutation matrices are examples.
• A nonzero nilpotent matrix is never diagonalizable: $m_A = t^s$ with $s \ge 2$.

Computing $m_A$. Look for the first linear dependence among $I, A, A^2, \dots$: flatten each power into a vector of length $n^2$ and stop at the first power that is a combination of the earlier ones. The coefficients of that combination give $m_A$. (Over $\mathbb{C}$, if the eigenvalues are known exactly, one can instead test exponents $s_i$ directly.)

import sympy as sp

def minimal_polynomial(M, t):
    """Find the first power M^d that is a combination of I, M, ..., M^(d-1)."""
    n = M.shape[0]
    powers = [sp.eye(n)]
    while True:
        powers.append(powers[-1] * M)
        cols = sp.Matrix.hstack(*[P.reshape(n * n, 1) for P in powers])   # vec(I), vec(M), ...
        null = cols.nullspace()
        if null:                                       # first dependence found
            c = null[0] / null[0][-1]                  # make it monic
            return sp.factor(sum(ci * t**i for i, ci in enumerate(c)))

t = sp.symbols('t')
A = sp.Matrix([[2, 0, 1, -1],
               [0, 2, 1, -1],
               [3, -3, 5, -6],
               [3, -3, 3, -4]])
print(sp.factor(A.charpoly(t).as_expr()))             # (t - 2)**3*(t + 1)
print(minimal_polynomial(A, t))                       # (t - 2)**2*(t + 1)
print(minimal_polynomial(sp.diag(2, 2, 3), t))        # (t - 3)*(t - 2): distinct roots -> diagonalizable
print(minimal_polynomial(sp.Matrix([[0, 0, 1], [1, 0, 0], [0, 1, 0]]), t))   # (t - 1)*(t**2 + t + 1)
print(A.is_diagonalizable())                          # False: the factor (t - 2) is repeated

The permutation matrix in the fourth line satisfies $A^3 = I$; its minimal polynomial $t^3 - 1 = (t - 1)(t^2 + t + 1)$ has three distinct complex roots, so it is diagonalizable over $\mathbb{C}$ but not over $\mathbb{R}$.

Key takeaways

• $m_A$ is the monic annihilating polynomial of least degree; it divides every annihilating polynomial, including $p_A$.
• $m_A$ and $p_A$ have the same roots, but $m_A$ may have smaller exponents.
• $A$ is diagonalizable exactly when $m_A$ splits into distinct linear factors.
• Projections, involutions, and matrices with $A^k = I$ (over $\mathbb{C}$) are diagonalizable; nonzero nilpotent matrices are not.
• Compute $m_A$ by finding the first linear dependence among $I, A, A^2, \dots$.

🧪 Practice

1. Find the minimal polynomials of $\begin{bmatrix} 3 & 1 & 0 \ 0 & 3 & 0 \ 0 & 0 & 3 \end{bmatrix}$, $\begin{bmatrix} 3 & 1 & 0 \ 0 & 3 & 1 \ 0 & 0 & 3 \end{bmatrix}$, and $\begin{bmatrix} 3 & 0 & 0 \ 0 & 3 & 0 \ 0 & 0 & 4 \end{bmatrix}$, and say which are diagonalizable.
2. Show that if $A^3 = A$ then $A$ is diagonalizable over $\mathbb{R}$, and list its possible eigenvalues.
3. Prove that the dimension of every Krylov subspace $\operatorname{span}{\mathbf{x}, A\mathbf{x}, A^2\mathbf{x}, \dots}$ is at most $\deg m_A$.
4. Interview: A matrix satisfies $A^2 = 5A - 6I$. Without knowing $A$, what can you say about its eigenvalues, diagonalizability, invertibility, and $A^{-1}$? Hint: $t^2 - 5t + 6 = (t - 2)(t - 3)$ annihilates $A$, so $m_A$ divides it: eigenvalues in ${2, 3}$, diagonalizable, invertible, and $A^{-1} = (5I - A)/6$.

Generalized Eigenvectors

Why it matters. A defective matrix does not have enough eigenvectors to form a basis. Generalized eigenvectors fill the gap: they are not mapped to multiples of themselves, but they reach an eigenvector after a few applications of $A - \lambda I$. They always supply a full basis, they explain the $te^{\lambda t}$ terms in differential equations with repeated roots, and they are the columns of the matrix $P$ in the Jordan form.

Analogy. An eigenvector is a ground-floor room: one application of $A - \lambda I$ sends it to $\mathbf{0}$. A generalized eigenvector of rank 2 is one floor up: one application takes it down to the ground floor (an eigenvector), and a second takes it to $\mathbf{0}$. A Jordan chain is a staircase that ends at the ground.

Definitions. Let $\lambda$ be an eigenvalue of $A$ and $N = A - \lambda I$.

• A nonzero $\mathbf{v}$ is a generalized eigenvector of rank $k$ for $\lambda$ if $N^k\mathbf{v} = \mathbf{0}$ but $N^{k-1}\mathbf{v} \ne \mathbf{0}$. Rank-1 generalized eigenvectors are ordinary eigenvectors.
• The generalized eigenspace is $G_\lambda = {\mathbf{v} : N^k\mathbf{v} = \mathbf{0} \text{ for some } k} = \operatorname{Nul}(A - \lambda I)^{n}$.
• A Jordan chain of length $k$ starting from a rank-$k$ vector $\mathbf{v}_k$ is $\mathbf{v}k$, $\mathbf{v}{k-1} = N\mathbf{v}_k$, ..., $\mathbf{v}_1 = N^{k-1}\mathbf{v}_k$. The bottom vector $\mathbf{v}_1$ is an eigenvector.

   a Jordan chain of length 3 for eigenvalue lambda   (N = A - lambda I)

      v3  --N-->  v2  --N-->  v1  --N-->  0
     rank 3      rank 2      rank 1 (an ordinary eigenvector)

   in terms of A:   A v1 = lambda v1
                    A v2 = lambda v2 + v1
                    A v3 = lambda v3 + v2

Nested kernels. $\operatorname{Nul}N \subseteq \operatorname{Nul}N^2 \subseteq \operatorname{Nul}N^3 \subseteq \cdots$, since $N^k\mathbf{v} = \mathbf{0}$ implies $N^{k+1}\mathbf{v} = \mathbf{0}$. The dimensions grow until they stop, and once two consecutive kernels are equal, all later ones are equal too. The final kernel is $G_\lambda$, and it is reached by the power $s$ = the exponent of $(t - \lambda)$ in the minimal polynomial.

Theorem (generalized eigenspace decomposition). Over $\mathbb{C}$ (or any field where $p_A$ splits), with distinct eigenvalues $\lambda_1, \dots, \lambda_k$:

1. $\dim G_{\lambda_i} = \operatorname{am}(\lambda_i)$;
2. each $G_{\lambda_i}$ is invariant, and the only eigenvalue of $A$ restricted to it is $\lambda_i$;
3. $F^n = G_{\lambda_1} \oplus G_{\lambda_2} \oplus \cdots \oplus G_{\lambda_k}$.

Idea of proof. Write $p_A = \prod(t - \lambda_i)^{a_i}$. The factors $(t - \lambda_i)^{a_i}$ are pairwise coprime, so by the extended Euclidean algorithm there are polynomials $h_i$ with $\sum_i h_i(t)\prod_{j \ne i}(t - \lambda_j)^{a_j} = 1$. Substituting $A$ and using $p_A(A) = 0$ (Cayley-Hamilton) gives projections onto the $\operatorname{Nul}(A - \lambda_iI)^{a_i}$ that sum to $I$, which is the same argument as the Lagrange proof in the previous topic. Restricting to $G_{\lambda_i}$, the operator $A - \lambda_iI$ is nilpotent, so $\lambda_i$ is the only eigenvalue there, and comparing characteristic polynomials of the block-diagonal form gives the dimensions. $\blacksquare$

This is the key fact for defective matrices: eigenvectors may fail to span, but generalized eigenvectors always do.

Worked example. $A = \begin{bmatrix} 1 & 1 \ -1 & 3 \end{bmatrix}$ has $p_A(\lambda) = \lambda^2 - 4\lambda + 4 = (\lambda - 2)^2$, and

$$N = A - 2I = \begin{bmatrix} -1 & 1 \ -1 & 1 \end{bmatrix}, \qquad N^2 = 0.$$

$\operatorname{Nul}N = \operatorname{span}{(1, 1)}$ has dimension 1, so $A$ is defective, while $\operatorname{Nul}N^2 = \mathbb{R}^2 = G_2$. Choose any $\mathbf{v}_2 \notin \operatorname{Nul}N$, say $\mathbf{v}_2 = (1, 0)$. Then $\mathbf{v}_1 = N\mathbf{v}_2 = (-1, -1)$ is an eigenvector, and $(\mathbf{v}_1, \mathbf{v}_2)$ is a basis of $\mathbb{R}^2$ in which $A$ has the matrix $\begin{bmatrix} 2 & 1 \ 0 & 2 \end{bmatrix}$.

Why $te^{\lambda t}$ appears. For $\mathbf{x}' = A\mathbf{x}$, an eigenvector gives the solution $e^{\lambda t}\mathbf{v}_1$. For the next vector in a chain, $e^{tA}\mathbf{v}_2 = e^{\lambda t}e^{tN}\mathbf{v}_2 = e^{\lambda t}(\mathbf{v}_2 + tN\mathbf{v}_2) = e^{\lambda t}(\mathbf{v}_2 + t\mathbf{v}_1)$, since $N^2\mathbf{v}_2 = \mathbf{0}$ cuts the exponential series short. This is the matrix origin of the $te^{\lambda t}$ solutions for repeated roots in differential equations courses.

import sympy as sp

A = sp.Matrix([[2, 0, 1, -1],
               [0, 2, 1, -1],
               [3, -3, 5, -6],
               [3, -3, 3, -4]])
N = A - 2 * sp.eye(4)

# Nested kernels grow and then stop: dims of Nul N^k
print([len((N ** k).nullspace()) for k in range(1, 5)])   # [2, 3, 3, 3]: G_2 has dim 3 = am(2), reached at k = 2

# Build a chain of length 2: pick v2 in Nul N^2 but not in Nul N
v2 = next(v for v in (N ** 2).nullspace() if N * v != sp.zeros(4, 1))
v1 = N * v2
print(v2.T, v1.T)                                 # Matrix([[-1, 0, 1, 0]]) Matrix([[1, 1, 0, 0]])
print(A * v1 == 2 * v1, A * v2 == 2 * v2 + v1)    # True True

# The 2x2 example from the text, and the te^{2t} solution
t = sp.symbols('t')
B = sp.Matrix([[1, 1], [-1, 3]])
w2 = sp.Matrix([1, 0]); w1 = (B - 2 * sp.eye(2)) * w2
x = sp.exp(2 * t) * (w2 + t * w1)                 # claimed solution of x' = B x
print(w1.T, sp.simplify(x.diff(t) - B * x).T)     # Matrix([[-1, -1]]) Matrix([[0, 0]])

Key takeaways

• $\mathbf{v}$ is a generalized eigenvector of rank $k$ when $(A - \lambda I)^k\mathbf{v} = \mathbf{0}$ but $(A - \lambda I)^{k-1}\mathbf{v} \ne \mathbf{0}$.
• A Jordan chain $\mathbf{v}k \to \mathbf{v}{k-1} \to \cdots \to \mathbf{v}_1 \to \mathbf{0}$ under $N = A - \lambda I$ ends in an ordinary eigenvector, and $A\mathbf{v}_j = \lambda\mathbf{v}j + \mathbf{v}{j-1}$.
• The kernels of $N^k$ increase and stabilize at $G_\lambda$, whose dimension is $\operatorname{am}(\lambda)$.
• Generalized eigenspaces of distinct eigenvalues form a direct sum equal to the whole space, even for defective matrices.
• Chains of length 2 or more produce solutions such as $te^{\lambda t}$ in differential equations.

🧪 Practice

1. For $A = \begin{bmatrix} 3 & 1 \ -1 & 1 \end{bmatrix}$, find the eigenvalue, an eigenvector, and a generalized eigenvector of rank 2, and write $A$ in the resulting basis.
2. For $J_3(0) = \begin{bmatrix} 0 & 1 & 0 \ 0 & 0 & 1 \ 0 & 0 & 0 \end{bmatrix}$, compute $\operatorname{Nul}J^k$ for $k = 1, 2, 3$ and identify a Jordan chain.
3. Solve $\mathbf{x}' = A\mathbf{x}$ with $A = \begin{bmatrix} 1 & 1 \ -1 & 3 \end{bmatrix}$ and $\mathbf{x}(0) = (0, 1)$, using the chain from the worked example.
4. Interview: Explain why generalized eigenvectors always span $\mathbb{C}^n$ even when eigenvectors do not, and why numerical software still avoids computing them. Hint: the generalized eigenspaces have dimensions equal to the algebraic multiplicities, which sum to $n$; but the chain structure is discontinuous under perturbation (a tiny change makes the matrix diagonalizable with almost parallel eigenvectors), so Schur vectors are used instead.

Jordan Canonical Form

Why it matters. The Jordan form is the answer to the question "what is the simplest matrix similar to $A$?" for every complex matrix. It is diagonal when $A$ is diagonalizable, and otherwise it is as close to diagonal as possible, with 1s just above the diagonal where eigenvectors are missing. It classifies matrices up to similarity completely, and it gives exact formulas for $A^k$ and $e^{tA}$ for defective matrices.

Jordan blocks. The $k \times k$ Jordan block for $\lambda$ is

$$J_k(\lambda) = \begin{bmatrix} \lambda & 1 & & \ & \lambda & \ddots & \ & & \ddots & 1 \ & & & \lambda \end{bmatrix} = \lambda I_k + N_k,$$

where $N_k$ has 1s on the superdiagonal. $N_k$ shifts coordinates ($\mathbf{e}j \mapsto \mathbf{e}{j-1}$), so $N_k^k = 0$. A Jordan block is exactly the matrix of $A$ on one Jordan chain $(\mathbf{v}_1, \dots, \mathbf{v}_k)$, by the relations $A\mathbf{v}_j = \lambda\mathbf{v}j + \mathbf{v}{j-1}$ from the previous topic. It has a single eigenvalue $\lambda$ with $\operatorname{am} = k$, $\operatorname{gm} = 1$, and minimal polynomial $(t - \lambda)^k$.

Theorem (Jordan canonical form). Every complex square matrix $A$ is similar to a block diagonal matrix

$$J = P^{-1}AP = \begin{bmatrix} J_{k_1}(\lambda_1) & & \ & \ddots & \ & & J_{k_r}(\lambda_r) \end{bmatrix},$$

where the $\lambda_i$ need not be distinct. $J$ is unique up to the order of the blocks. The columns of $P$ are Jordan chains, each listed from its eigenvector upward.

The existence proof combines the previous topic (split into generalized eigenspaces) with a combinatorial argument that each nilpotent restriction $N = (A - \lambda I)|{G\lambda}$ has a basis of chains. Uniqueness follows from the counting formula below, which expresses the block sizes through ranks, which are similarity invariants.

Reading the Jordan form. For each eigenvalue $\lambda$:

Quantity	Read from $J$
algebraic multiplicity	total size of all blocks for $\lambda$
geometric multiplicity	number of blocks for $\lambda$
exponent of $(t - \lambda)$ in $m_A$	size of the largest block for $\lambda$
diagonalizable	every block is $1 \times 1$

Counting the blocks from ranks. Let $r_k = \operatorname{rank}(A - \lambda I)^k$, with $r_0 = n$. Each block for $\lambda$ of size $\ge k$ loses exactly one unit of rank between powers $k - 1$ and $k$, and blocks for other eigenvalues keep full rank. So

$$#{\text{blocks for } \lambda \text{ of size} \ge k} = r_{k-1} - r_k, \qquad #{\text{blocks of size exactly } k} = r_{k-1} - 2r_k + r_{k+1}.$$

Running example. For the $4 \times 4$ matrix $A$ of this subchapter and $\lambda = 2$: $r_0 = 4$, $r_1 = 2$, $r_2 = 1$, $r_3 = 1$. Blocks of size $\ge 1$: $4 - 2 = 2$. Size $\ge 2$: $2 - 1 = 1$. Size $\ge 3$: $1 - 1 = 0$. So $\lambda = 2$ has one block of size 2 and one of size 1, and $\lambda = -1$ has a single $1 \times 1$ block:

          [ 2  1 |   |    ]
          [ 0  2 |   |    ]        dot diagram for lambda = 2:
   J  =   [------+---+----]           * *      <- block of size 2 (a chain of length 2)
          [      | 2 |    ]           *        <- block of size 1 (a lone eigenvector)
          [------+---+----]        rows = blocks, columns = rank of the chain vectors
          [      |   | -1 ]

All Jordan types of a $3 \times 3$ matrix with a single eigenvalue $\lambda$.

Block sizes	$J$	gm	$m_A$	Diagonalizable
$1, 1, 1$	$\lambda I$	3	$t - \lambda$	yes
$2, 1$	$\begin{bmatrix} \lambda & 1 & 0 \ 0 & \lambda & 0 \ 0 & 0 & \lambda \end{bmatrix}$	2	$(t - \lambda)^2$	no
$3$	$\begin{bmatrix} \lambda & 1 & 0 \ 0 & \lambda & 1 \ 0 & 0 & \lambda \end{bmatrix}$	1	$(t - \lambda)^3$	no

For $3 \times 3$ matrices the geometric multiplicity alone distinguishes the types; from $4 \times 4$ on, the partitions $(2, 2)$ and $(3, 1)$ have the same $\operatorname{gm} = 2$ and need the ranks of $(A - \lambda I)^2$.

Powers and exponentials of Jordan blocks. Since $\lambda I$ and $N$ commute and $N^k = 0$, the binomial theorem and the exponential series stop after $k$ terms:

$$J_k(\lambda)^m = \sum_{j=0}^{k-1}\binom{m}{j}\lambda^{m-j}N^j, \qquad e^{tJ_k(\lambda)} = e^{\lambda t}\Big(I + tN + \frac{t^2}{2!}N^2 + \cdots + \frac{t^{k-1}}{(k-1)!}N^{k-1}\Big).$$

For example, $J_2(\lambda)^m = \begin{bmatrix} \lambda^m & m\lambda^{m-1} \ 0 & \lambda^m \end{bmatrix}$. The off-diagonal term $m\lambda^{m-1}$ means that a defective matrix with $\lvert\lambda\rvert = 1$ produces polynomial growth, and one with $\lvert\lambda\rvert < 1$ can grow for a while before it decays.

A theoretical tool, not a numerical one. The Jordan form depends discontinuously on the entries: $\begin{bmatrix} 0 & 1 \ \varepsilon & 0 \end{bmatrix}$ is diagonalizable for every $\varepsilon \ne 0$ but is a Jordan block at $\varepsilon = 0$. In floating point, rounding errors destroy Jordan structure, so it should be computed only in exact arithmetic (as sympy does). Numerical work uses the Schur form instead.

import sympy as sp

A = sp.Matrix([[2, 0, 1, -1],
               [0, 2, 1, -1],
               [3, -3, 5, -6],
               [3, -3, 3, -4]])
P, J = A.jordan_form()
print(J)                              # Matrix([[-1, 0, 0, 0], [0, 2, 1, 0], [0, 0, 2, 0], [0, 0, 0, 2]])
print(P * J * P.inv() == A)           # True

def jordan_block_sizes(M, lam):
    """Block sizes for eigenvalue lam from the ranks r_k = rank (M - lam I)^k."""
    n = M.shape[0]
    r = [n] + [((M - lam * sp.eye(n)) ** k).rank() for k in range(1, n + 2)]
    sizes = []
    for k in range(1, n + 1):
        sizes += [k] * (r[k - 1] - 2 * r[k] + r[k + 1])   # number of blocks of size exactly k
    return sorted(sizes, reverse=True)

print(jordan_block_sizes(A, 2), jordan_block_sizes(A, -1))   # [2, 1] [1]

# Powers of a Jordan block: binomial coefficients appear above the diagonal
lam, m = sp.symbols('lambda m')
J3 = sp.Matrix([[lam, 1, 0], [0, lam, 1], [0, 0, lam]])
print(sp.simplify(J3 ** 5))
# Matrix([[lambda**5, 5*lambda**4, 10*lambda**3], [0, lambda**5, 5*lambda**4], [0, 0, lambda**5]])

# |lambda| < 1 but transient growth: J_2(0.9)^m peaks before decaying
J2 = sp.Matrix([[sp.Rational(9, 10), 1], [0, sp.Rational(9, 10)]])
print([round(float((J2 ** k)[0, 1]), 2) for k in [1, 5, 10, 20, 60]])   # [1.0, 3.28, 3.87, 2.7, 0.12]

Key takeaways

• Every complex matrix is similar to a Jordan matrix, a block diagonal of Jordan blocks $J_k(\lambda) = \lambda I + N$, unique up to block order.
• The columns of $P$ are Jordan chains; each block corresponds to one chain.
• Number of blocks for $\lambda$ = geometric multiplicity; total size = algebraic multiplicity; largest size = exponent in $m_A$.
• Block sizes are determined by the ranks of $(A - \lambda I)^k$, so two matrices are similar exactly when their Jordan forms agree.
• Powers and exponentials of Jordan blocks contain $\binom{m}{j}\lambda^{m-j}$ and $t^je^{\lambda t}/j!$; the form is unstable numerically and is computed only exactly.

🧪 Practice

1. Find the Jordan form of $\begin{bmatrix} 1 & 1 \ -1 & 3 \end{bmatrix}$ and of $\begin{bmatrix} 5 & 1 & 0 \ 0 & 5 & 0 \ 0 & 0 & 5 \end{bmatrix}$, including the matrix $P$.
2. A $6 \times 6$ matrix has single eigenvalue $3$ with $\operatorname{rank}(A - 3I) = 3$, $\operatorname{rank}(A - 3I)^2 = 1$, $\operatorname{rank}(A - 3I)^3 = 0$. Find its Jordan form and minimal polynomial.
3. List all possible Jordan forms (up to block order) of a $4 \times 4$ matrix with characteristic polynomial $(t - 1)^2(t + 2)^2$, and give the minimal polynomial of each.
4. Interview: A discrete-time control system has all eigenvalues inside the unit disk, yet simulations show the state growing by a factor of 50 before decaying. How is that possible? Hint: non-normality: Jordan blocks (or nearly defective, non-orthogonal eigenvectors) give terms like $m\lambda^{m-1}$ that grow before $\lambda^m$ wins; eigenvalues only govern the long run.

Rational Canonical Form

Why it matters. The Jordan form requires all eigenvalues to lie in the field. Over $\mathbb{Q}$ or $\mathbb{R}$, or over a finite field, they often do not: $\begin{bmatrix} 0 & 2 \ 1 & 0 \end{bmatrix}$ has eigenvalues $\pm\sqrt2$, which are not rational. The rational canonical form is a canonical form that uses only arithmetic in the original field, so it exists for every matrix over every field, and two matrices over $F$ are similar exactly when their rational canonical forms are equal. It is computed with exact polynomial arithmetic, without finding any roots.

Building block: the companion matrix. For a monic $f(t) = t^d + c_{d-1}t^{d-1} + \cdots + c_0$, the companion matrix $C(f)$ (Section 8.1) satisfies $p_{C(f)} = m_{C(f)} = f$. It has a cyclic vector: $\mathbf{e}_1 \mapsto \mathbf{e}_2 \mapsto \cdots \mapsto \mathbf{e}_d$, so $\mathbf{e}_1, C\mathbf{e}_1, \dots, C^{d-1}\mathbf{e}_1$ is a basis. Conversely, if $A$ has a cyclic vector $\mathbf{x}$ (a vector whose Krylov subspace is all of $F^n$), then in the basis $\mathbf{x}, A\mathbf{x}, \dots, A^{n-1}\mathbf{x}$ the matrix of $A$ is $C(p_A)$. A matrix has a cyclic vector exactly when $m_A = p_A$.

Theorem (rational canonical form). Every square matrix $A$ over a field $F$ is similar over $F$ to a unique block diagonal matrix

$$R = \begin{bmatrix} C(f_1) & & \ & \ddots & \ & & C(f_r) \end{bmatrix}, \qquad f_1 \mid f_2 \mid \cdots \mid f_r,$$

where the $f_i$ are monic, nonconstant polynomials over $F$, called the invariant factors of $A$. Moreover $f_r = m_A$ and $f_1f_2\cdots f_r = p_A$.

The idea is to split $F^n$ into cyclic invariant subspaces, each a Krylov subspace on which $A$ acts as a companion matrix. The divisibility condition makes the decomposition unique. The invariant factors can be computed from the Smith normal form of the polynomial matrix $tI - A$, by row and column operations over $F[t]$, which is Gaussian elimination with polynomial entries.

Running example. For the $4 \times 4$ matrix $A$ of this subchapter, $p_A = (t - 2)^3(t + 1)$ and $m_A = (t - 2)^2(t + 1) = t^3 - 3t^2 + 4$. The invariant factors must divide each other, end with $m_A$, and multiply to $p_A$, which forces $f_1 = t - 2$, $f_2 = t^3 - 3t^2 + 4$. Therefore

$$R = \begin{bmatrix} 2 & 0 & 0 & 0 \ 0 & 0 & 0 & -4 \ 0 & 1 & 0 & 0 \ 0 & 0 & 1 & 3 \end{bmatrix}.$$

The last column of the companion block holds $-c_0, -c_1, -c_2 = -4, 0, 3$.

Comparison of canonical and standard forms.

Form	Field requirement	Unique?	Classifies similarity?	Numerically stable?
diagonal $PDP^{-1}$	eigenvalues in $F$ and enough eigenvectors	up to order	only among diagonalizable matrices	depends on $\operatorname{cond}(P)$
Jordan $PJP^{-1}$	eigenvalues in $F$	up to block order	yes	no
rational canonical	any field	yes	yes, over $F$ itself	no (exact arithmetic)
Schur $UTU^*$	$\mathbb{C}$ (real Schur for $\mathbb{R}$)	no	no	yes

A useful consequence: similarity does not depend on the field. If two real (or rational) matrices are similar over $\mathbb{C}$, they are similar over $\mathbb{R}$ (or $\mathbb{Q}$), because the invariant factors are computed inside the smaller field.

import sympy as sp
from sympy.matrices.normalforms import invariant_factors

t = sp.symbols('t')

def companion(f):
    """Companion matrix of a monic polynomial f (ones below the diagonal, -coefficients in the last column)."""
    c = sp.Poly(f, t).all_coeffs()[::-1][:-1]           # c_0, c_1, ..., c_{d-1}
    d = len(c)
    C = sp.zeros(d, d)
    for i in range(1, d):
        C[i, i - 1] = 1
    for i in range(d):
        C[i, d - 1] = -c[i]
    return C

def rational_canonical_form(M):
    n = M.shape[0]
    fs = [f for f in invariant_factors(t * sp.eye(n) - M, domain=sp.QQ[t]) if sp.degree(f, t) > 0]
    return fs, sp.diag(*[companion(f) for f in fs])

A = sp.Matrix([[2, 0, 1, -1],
               [0, 2, 1, -1],
               [3, -3, 5, -6],
               [3, -3, 3, -4]])
fs, R = rational_canonical_form(A)
print(fs)                                       # [t - 2, t**3 - 3*t**2 + 4]
print(R.tolist())                               # [[2, 0, 0, 0], [0, 0, 0, -4], [0, 1, 0, 0], [0, 0, 1, 3]]
print(R.jordan_form()[1] == A.jordan_form()[1]) # True: same Jordan form, so A and R are similar

# Irrational eigenvalues: the rational form stays in Q, the Jordan form does not
B = sp.Matrix([[1, 2], [1, -1]])                 # eigenvalues +- sqrt(3)
print(rational_canonical_form(B))               # ([t**2 - 3], Matrix([[0, 3], [1, 0]]))
print(B.jordan_form()[1])                       # Matrix([[-sqrt(3), 0], [0, sqrt(3)]])

Key takeaways

• The rational canonical form uses only the field's own arithmetic, so it exists over every field, including $\mathbb{Q}$ and finite fields.
• Its blocks are companion matrices of the invariant factors $f_1 \mid f_2 \mid \cdots \mid f_r$, with $f_r = m_A$ and $\prod f_i = p_A$.
• A matrix has a cyclic vector exactly when $m_A = p_A$, and then it is similar to the single companion matrix $C(p_A)$.
• Invariant factors come from the Smith normal form of $tI - A$; no roots are needed.
• Similarity over a larger field implies similarity over the original field.

🧪 Practice

1. Find the rational canonical form of $\operatorname{diag}(1, 1, 2)$ and of the $3 \times 3$ Jordan block $J_3(1)$. Which one has a cyclic vector?
2. A $5 \times 5$ rational matrix has $p_A = (t^2 + 1)^2(t - 3)$ and $m_A = (t^2 + 1)(t - 3)$. Find its invariant factors and write its rational canonical form.
3. Show that $\mathbf{x} = \mathbf{e}_1$ is a cyclic vector for $\begin{bmatrix} 0 & 2 \ 1 & 0 \end{bmatrix}$, and that $\mathbf{x} = \mathbf{e}_1$ is not cyclic for $2I$. Explain the second fact using the minimal polynomial.
4. Interview: In coding theory or cryptography, a linear map over $\mathbb{F}_2$ is used repeatedly (for example, a linear feedback shift register). Why is the rational canonical form, rather than the Jordan form, the natural way to understand its period? Hint: eigenvalues may lie only in extension fields of $\mathbb{F}_2$; the period of an LFSR is the order of $t$ modulo its (companion) characteristic polynomial, a statement purely in $\mathbb{F}_2[t]$.

<a id="84-applications-of-eigenvalues"></a>

8.4 Applications of Eigenvalues

This subchapter applies the eigen-machinery to processes that evolve in steps or in continuous time: dynamical systems, random walks, differential equations, structured populations, and rankings of nodes in a network.

Discrete Dynamical Systems

Why it matters. A discrete linear dynamical system $\mathbf{x}_{k+1} = A\mathbf{x}_k$ models anything updated once per period: interacting populations, a loan portfolio, a digital filter, a linearized economy, a recurrent network's hidden state. The questions are always the same. Will the state die out, settle, oscillate, or explode? Which pattern dominates in the long run, and how fast does it take over? The eigenvalues answer all of them.

The solution. If $A$ has an eigenbasis $\mathbf{v}_1, \dots, \mathbf{v}_n$ (Section 8.2) and $\mathbf{x}_0 = c_1\mathbf{v}_1 + \cdots + c_n\mathbf{v}_n$, then

$$\mathbf{x}_k = A^k\mathbf{x}_0 = c_1\lambda_1^k\mathbf{v}_1 + c_2\lambda_2^k\mathbf{v}_2 + \cdots + c_n\lambda_n^k\mathbf{v}_n.$$

Each eigen-direction evolves on its own as a geometric sequence. Order the eigenvalues so that $\lvert\lambda_1\rvert \ge \lvert\lambda_2\rvert \ge \cdots$.

• Stability. $\mathbf{x}_k \to \mathbf{0}$ for every $\mathbf{x}_0$ exactly when the spectral radius $\rho(A) = \max\lvert\lambda_i\rvert$ is less than 1. This holds for defective matrices too, by the Jordan form (Section 8.3): $k^j\lambda^k \to 0$ when $\lvert\lambda\rvert < 1$.
• Dominance. If $\lvert\lambda_1\rvert > \lvert\lambda_2\rvert$ and $c_1 \ne 0$, then $\mathbf{x}_k \approx c_1\lambda_1^k\mathbf{v}_1$: the direction of $\mathbf{x}_k$ approaches $\mathbf{v}_1$, and the state is multiplied by about $\lambda_1$ per step.
• Rate. The error in that approximation shrinks like $\lvert\lambda_2/\lambda_1\rvert^k$.

Classifying the origin for $2 \times 2$ systems.

Eigenvalues	Name	Trajectories
real, both $\lvert\lambda\rvert < 1$	attractor (sink)	all approach $\mathbf{0}$, eventually along the eigenvector of the larger $\lvert\lambda\rvert$
real, both $\lvert\lambda\rvert > 1$	repeller (source)	all move away from $\mathbf{0}$
real, $\lvert\lambda_1\rvert > 1 > \lvert\lambda_2\rvert$	saddle point	approach along $\mathbf{v}_2$, then leave along $\mathbf{v}_1$
complex, $\lvert\lambda\rvert < 1$	spiral attractor	spiral inward (Section 8.1)
complex, $\lvert\lambda\rvert > 1$	spiral repeller	spiral outward
complex, $\lvert\lambda\rvert = 1$	rotation	stay on an ellipse

A negative real eigenvalue makes the corresponding component flip sign every step, so trajectories jump from one side of the origin to the other.

   saddle point: lambda_1 = 1.5 along v1, lambda_2 = 0.5 along v2

                  v2
           \       |       /
            \      v      /          points start near the v2-axis,
       <---- \     |     / ---->     are pulled in along v2,
              `--. | .--'            then pushed out along v1
   v1 <----------- O -----------> v1
              .--' | `--.
       <---- /     ^     \ ---->
            /      |      \
           /       |       \

Example: owls and rats. Let $O_k$ be the number of owls and $R_k$ the number of rats (in thousands) in month $k$, with

$$\begin{bmatrix} O_{k+1} \ R_{k+1} \end{bmatrix} = \begin{bmatrix} 0.5 & 0.4 \ -0.104 & 1.1 \end{bmatrix}\begin{bmatrix} O_k \ R_k \end{bmatrix}.$$

Without rats, half the owls die each month; rats feed the owls ($0.4$); owls eat rats ($-0.104$); without owls, rats grow 10% a month. The characteristic polynomial $\lambda^2 - 1.6\lambda + 0.5916$ has roots $1.02$ and $0.58$, with eigenvectors $(10, 13)$ and $(5, 1)$. In the long run, both populations grow by 2% per month, in the ratio of 10 owls to every 13 thousand rats. The $0.58^k$ component dies out within a few months; after 20 months it is about $0.58^{20} \approx 2 \times 10^{-5}$ of its starting size.

Affine systems. For $\mathbf{x}_{k+1} = A\mathbf{x}_k + \mathbf{b}$, the fixed point $\mathbf{x}^* = (I - A)^{-1}\mathbf{b}$ exists when $1$ is not an eigenvalue. The deviation $\mathbf{e}_k = \mathbf{x}k - \mathbf{x}^*$ satisfies $\mathbf{e}{k+1} = A\mathbf{e}_k$, so the same eigenvalue analysis decides whether the system settles at $\mathbf{x}^*$.

import numpy as np

A = np.array([[0.5, 0.4],
              [-0.104, 1.1]])
vals, vecs = np.linalg.eig(A)
vals, vecs = vals.real, vecs.real
print(np.round(vals, 4))                               # [0.58 1.02]
v1 = vecs[:, np.argmax(np.abs(vals))]
print(np.round(v1 / v1[0] * 10, 4))                    # [10. 13.]: owls : rats = 10 : 13

x = np.array([20.0, 100.0])                            # start: 20 owls, 100 thousand rats
for k in range(1, 61):
    x_prev, x = x, A @ x
    if k in (1, 10, 60):
        print(k, np.round(x, 1), np.round(x / x_prev, 4), round(x[1] / x[0], 4))
# 1 [ 50.  107.9] [2.5    1.0792] 2.1584
# 10 [106.1 138.2] [1.0221 1.0203] 1.303
# 60 [286.3 372.2] [1.02 1.02] 1.3

def classify(M):
    lam = np.linalg.eigvals(M)
    r = np.abs(lam)
    if np.iscomplexobj(lam) and np.any(np.abs(lam.imag) > 1e-12):
        return "spiral attractor" if r[0] < 1 else "spiral repeller" if r[0] > 1 else "rotation"
    if np.all(r < 1):
        return "attractor"
    if np.all(r > 1):
        return "repeller"
    return "saddle" if r.min() < 1 < r.max() else "borderline"

print(classify(A))                                     # saddle
print(classify(np.array([[0.8, -0.3], [0.3, 0.8]])))   # spiral attractor
print(classify(np.array([[0.9, 0.0], [0.0, 0.5]])))    # attractor

Key takeaways

• $\mathbf{x}_k = \sum c_i\lambda_i^k\mathbf{v}_i$: each eigen-direction evolves independently as a geometric sequence.
• The system is stable ($\mathbf{x}_k \to \mathbf{0}$ for all starts) exactly when $\rho(A) < 1$.
• A dominant eigenvalue determines the long-run growth factor and the long-run proportions (its eigenvector).
• The ratio $\lvert\lambda_2/\lambda_1\rvert$ sets how fast the dominant pattern takes over.
• Real eigenvalues give attractors, repellers, and saddles; complex ones give spirals and rotations.

🧪 Practice

1. Classify the origin for $\begin{bmatrix} 1.2 & 0 \ 0 & 0.6 \end{bmatrix}$, $\begin{bmatrix} 0.6 & 0.3 \ 0.3 & 0.6 \end{bmatrix}$, and $\begin{bmatrix} 0.6 & -0.8 \ 0.8 & 0.6 \end{bmatrix}$, and sketch typical trajectories.
2. In the owl-rat model, change the predation rate $0.104$ to $p$. For which $p$ do both populations eventually die out, and for which $p$ do they stabilize at constant levels?
3. Find the fixed point of $\mathbf{x}_{k+1} = \begin{bmatrix} 0.5 & 0.2 \ 0.1 & 0.7 \end{bmatrix}\mathbf{x}_k + \begin{bmatrix} 1 \ 2 \end{bmatrix}$ and decide whether iterates converge to it.
4. Interview: A linearized model of an economy has a $50 \times 50$ update matrix with spectral radius $0.97$. Your manager asks how long until a shock decays to 1% of its size. What do you answer, and what caveat do you add? Hint: $0.97^k = 0.01$ gives $k \approx 151$ periods; but if $A$ is far from normal (badly conditioned eigenvectors), transient growth can delay this, so simulate or bound $\lVert A^k\rVert$ directly.

Markov Chains and Steady States

Why it matters. A Markov chain models a system that moves randomly between states, where the next move depends only on the current state: customers switching brands, a web surfer clicking links, a molecule changing conformation, a word generator choosing the next token. The long-run fraction of time spent in each state is an eigenvector for the eigenvalue 1, and the second eigenvalue says how fast the chain forgets its starting point. Markov chain Monte Carlo, used throughout Bayesian statistics, depends on both facts.

Setup. With $n$ states, let $p_{ij}$ be the probability of moving to state $i$ from state $j$ in one step. The transition matrix $P = [p_{ij}]$ has nonnegative entries and columns that sum to 1 (it is column-stochastic). A probability vector $\mathbf{x}_k$ (nonnegative, entries summing to 1) gives the distribution at step $k$, and

$$\mathbf{x}_{k+1} = P\mathbf{x}_k, \qquad \mathbf{x}_k = P^k\mathbf{x}_0.$$

Many statistics and probability texts use the transposed convention: rows sum to 1 and distributions are row vectors, $\boldsymbol{\pi}_{k+1}^T = \boldsymbol{\pi}_k^TP$. Everything below transfers by transposing.

A steady state always exists. A steady-state (or stationary) vector is a probability vector $\mathbf{q}$ with $P\mathbf{q} = \mathbf{q}$: an eigenvector for eigenvalue 1.

• 1 is an eigenvalue. Column sums equal 1 means $\mathbf{1}^TP = \mathbf{1}^T$, so $\mathbf{1}$ is an eigenvector of $P^T$ with eigenvalue 1, and $P$ and $P^T$ have the same eigenvalues (Section 8.1).
• No eigenvalue is larger. If $P\mathbf{v} = \lambda\mathbf{v}$, summing absolute values gives $\lvert\lambda\rvert\sum_i\lvert v_i\rvert \le \sum_i\sum_jp_{ij}\lvert v_j\rvert = \sum_j\lvert v_j\rvert$, so $\lvert\lambda\rvert \le 1$.
• A nonnegative eigenvector for 1 exists (Perron-Frobenius theory, Section 13.2), so it can be scaled into a probability vector.

Convergence for regular chains. A stochastic matrix is regular if some power $P^m$ has all entries strictly positive (from every state you can reach every state in exactly $m$ steps). For a regular chain:

1. the steady state $\mathbf{q}$ is unique and strictly positive;
2. every other eigenvalue satisfies $\lvert\lambda\rvert < 1$;
3. $\mathbf{x}_k \to \mathbf{q}$ for every starting distribution, and $P^k \to \mathbf{q}\mathbf{1}^T$, the matrix whose columns all equal $\mathbf{q}$;
4. the distance to $\mathbf{q}$ shrinks roughly like $\lvert\lambda_2\rvert^k$, so about $\ln(1/\varepsilon)/\ln(1/\lvert\lambda_2\rvert)$ steps bring it below $\varepsilon$. The quantity $1 - \lvert\lambda_2\rvert$ is called the spectral gap.

The proof of 3 is the eigenvector expansion: $\mathbf{x}_0 = \mathbf{q} + c_2\mathbf{v}_2 + \cdots$ (the coefficient of $\mathbf{q}$ is 1 because all the other eigenvectors have entries summing to 0, since $\mathbf{1}^T\mathbf{v} = \mathbf{1}^TP\mathbf{v}/\lambda$ forces $\mathbf{1}^T\mathbf{v} = 0$ when $\lambda \ne 1$), and every $\lambda_i^k \to 0$ for $i \ge 2$.

When things go wrong.

Chain	Eigenvalues on the unit circle	Behavior
regular	only $1$, simple	converges to a unique $\mathbf{q}$
periodic, e.g. $\begin{bmatrix} 0 & 1 \ 1 & 0 \end{bmatrix}$	$1$ and $-1$	oscillates forever; $\mathbf{q} = (\tfrac12, \tfrac12)$ is only a time average
reducible, e.g. two separate groups of states	$1$ with multiplicity 2	several steady states; the limit depends on the start
absorbing states (e.g. "game over")	$1$ for each absorbing state	ends in an absorbing state; eigen-analysis gives absorption probabilities

Example: brand switching. Three brands A, B, C. Each month, customers switch according to

$$P = \begin{bmatrix} 0.8 & 0.1 & 0.2 \ 0.1 & 0.7 & 0.3 \ 0.1 & 0.2 & 0.5 \end{bmatrix} \qquad (\text{column } j: \text{where brand } j\text{'s customers go}).$$

Solving $(P - I)\mathbf{q} = \mathbf{0}$ and normalizing gives $\mathbf{q} = \big(\tfrac{9}{22}, \tfrac{8}{22}, \tfrac{5}{22}\big) \approx (0.409, 0.364, 0.227)$ (computed below). Whatever the current market shares, they approach these values. The other two eigenvalues are the roots of $\lambda^2 - \lambda + 0.22$, namely $\tfrac{1}{2}(1 \pm \sqrt{0.12}) \approx 0.673$ and $0.327$, so the distance to the steady state shrinks by a factor of about $0.673$ per month and falls below $10^{-6}$ of its initial size after about 35 months.

   transition probabilities (column = from, row = to)        long-run shares

              from A   from B   from C
     to A      0.8      0.1      0.2                        A  ####################  40.9%
     to B      0.1      0.7      0.3                        B  ##################    36.4%
     to C      0.1      0.2      0.5                        C  ###########           22.7%

import numpy as np

P = np.array([[0.8, 0.1, 0.2],
              [0.1, 0.7, 0.3],
              [0.1, 0.2, 0.5]])
print(P.sum(axis=0))                                  # [1. 1. 1.]: column-stochastic

vals, vecs = np.linalg.eig(P)
i = np.argmin(np.abs(vals - 1))                       # locate the eigenvalue 1
q = (vecs[:, i] / vecs[:, i].sum()).real              # scale the eigenvector to sum to 1
print(np.round(q, 4), np.round(q * 22, 4))            # [0.4091 0.3636 0.2273] [9. 8. 5.]
print(np.round(np.sort(np.abs(vals)), 4))             # [0.3268 0.6732 1.    ]: |lambda_2| ~ 0.673

# Convergence from two very different starts
for x in [np.array([1.0, 0.0, 0.0]), np.array([0.0, 0.0, 1.0])]:
    for _ in range(40):
        x = P @ x
    print(np.round(x, 4))                             # both [0.4091 0.3636 0.2273]

# P^k tends to the rank-one matrix with every column equal to q
print(np.allclose(np.linalg.matrix_power(P, 50), np.outer(q, np.ones(3))))   # True

# A periodic chain never settles
F = np.array([[0.0, 1.0], [1.0, 0.0]])
print(np.linalg.eigvals(F).real, np.linalg.matrix_power(F, 7) @ np.array([1.0, 0.0]))  # [ 1. -1.] [0. 1.]

Key takeaways

• A column-stochastic $P$ always has eigenvalue 1, and every eigenvalue satisfies $\lvert\lambda\rvert \le 1$.
• A steady state is an eigenvector for eigenvalue 1, scaled so its entries sum to 1.
• For a regular chain the steady state is unique and positive, and every starting distribution converges to it.
• The convergence rate is governed by $\lvert\lambda_2\rvert$; the spectral gap $1 - \lvert\lambda_2\rvert$ measures how fast the chain mixes.
• Periodic chains (eigenvalue $-1$ or other roots of unity) oscillate; reducible chains have several steady states.

🧪 Practice

1. Find the steady state of $\begin{bmatrix} 0.9 & 0.5 \ 0.1 & 0.5 \end{bmatrix}$ and the second eigenvalue. Roughly how many steps does it take for any start to be within $10^{-6}$ of the steady state?
2. Show that for a $2 \times 2$ stochastic matrix $\begin{bmatrix} 1 - a & b \ a & 1 - b \end{bmatrix}$, the eigenvalues are $1$ and $1 - a - b$, and the steady state is $\tfrac{1}{a + b}(b, a)$.
3. Construct a $3 \times 3$ reducible stochastic matrix with two different steady states, and describe which one the chain approaches from each starting state.
4. Interview: An MCMC sampler for a Bayesian model is "slow to mix". What does that mean in terms of the eigenvalues of its transition operator, and what are typical remedies? Hint: $\lvert\lambda_2\rvert$ is close to 1 (a small spectral gap), often because of nearly disconnected regions of high probability; remedies include better proposals, reparameterization, tempering, or Hamiltonian Monte Carlo.

Systems of Linear Differential Equations

Why it matters. Continuous-time processes such as heat flow between rooms, vibrating structures, electrical circuits, chemical kinetics, and the linearization of any smooth system near an equilibrium are modeled by $\mathbf{x}'(t) = A\mathbf{x}(t)$. Section 8.2 gave the formal solution $e^{tA}\mathbf{x}_0$. Eigenvalues turn that formula into an explicit sum of exponential modes and decide stability from their real parts.

Eigen-solutions. Try $\mathbf{x}(t) = e^{\lambda t}\mathbf{v}$ with a constant vector $\mathbf{v}$. Then $\mathbf{x}' = \lambda e^{\lambda t}\mathbf{v}$ and $A\mathbf{x} = e^{\lambda t}A\mathbf{v}$, so this is a solution exactly when $A\mathbf{v} = \lambda\mathbf{v}$. By linearity, if $A$ has an eigenbasis, the general solution is

$$\mathbf{x}(t) = c_1e^{\lambda_1t}\mathbf{v}_1 + c_2e^{\lambda_2t}\mathbf{v}_2 + \cdots + c_ne^{\lambda_nt}\mathbf{v}_n,$$

where the $c_i$ are the coordinates of $\mathbf{x}(0)$ in the eigenbasis. In eigen-coordinates, the coupled system becomes $n$ separate equations $c_i' = \lambda_ic_i$.

Example: heat exchange between two rooms. Room temperatures $x_1, x_2$ (relative to the outside) evolve as

$$\mathbf{x}' = \begin{bmatrix} -2 & 1 \ 1 & -2 \end{bmatrix}\mathbf{x}.$$

Each room loses heat to the outside (rate 1) and exchanges heat with the other (rate 1). The eigenpairs are $\lambda = -1$, $\mathbf{v} = (1, 1)$, and $\lambda = -3$, $\mathbf{v} = (1, -1)$. So

$$\mathbf{x}(t) = c_1e^{-t}\begin{bmatrix} 1 \ 1 \end{bmatrix} + c_2e^{-3t}\begin{bmatrix} 1 \ -1 \end{bmatrix}.$$

The difference between the rooms (the mode $(1, -1)$) disappears three times faster than the common temperature (the mode $(1, 1)$): the rooms first equalize, then cool together.

Complex eigenvalues give oscillations. For a real $A$ with eigenvalue $\lambda = a + bi$ and eigenvector $\mathbf{v}$, the complex solution $e^{\lambda t}\mathbf{v} = e^{at}(\cos bt + i\sin bt)\mathbf{v}$ has real and imaginary parts that are both real solutions:

$$e^{at}\big(\cos(bt)\operatorname{Re}\mathbf{v} - \sin(bt)\operatorname{Im}\mathbf{v}\big), \qquad e^{at}\big(\sin(bt)\operatorname{Re}\mathbf{v} + \cos(bt)\operatorname{Im}\mathbf{v}\big).$$

The imaginary part $b$ is the angular frequency, and the real part $a$ is the growth or decay rate of the envelope.

Higher-order equations become first-order systems. The damped oscillator $y'' + 2\zeta\omega y' + \omega^2y = 0$ becomes, with $\mathbf{x} = (y, y')$,

$$\mathbf{x}' = \begin{bmatrix} 0 & 1 \ -\omega^2 & -2\zeta\omega \end{bmatrix}\mathbf{x},$$

whose characteristic polynomial $\lambda^2 + 2\zeta\omega\lambda + \omega^2$ is the familiar characteristic equation of the ODE (the matrix is a companion matrix, Section 8.1). Underdamping ($\zeta < 1$) gives complex eigenvalues, critical damping ($\zeta = 1$) a repeated, defective eigenvalue with solutions $te^{-\omega t}$ (Section 8.3), and overdamping ($\zeta > 1$) two negative real eigenvalues.

Stability and the trace-determinant plane. The equilibrium $\mathbf{0}$ is asymptotically stable (all solutions tend to $\mathbf{0}$) exactly when every eigenvalue has $\operatorname{Re}\lambda < 0$. For $2 \times 2$ systems, the eigenvalues are determined by $\tau = \operatorname{tr}A$ and $\delta = \det A$, which gives a map of all possible behaviors:

                          det A = delta
                               ^
      spiral source            |            spiral sink
      (Re > 0, complex)        |       (Re < 0, complex)
                  .  `  .      |      .  `  .
   node source  .         `  . | .  `         .  node sink
   (real > 0)  .     tau^2 = 4 delta    .  (real < 0)
              .                |  (parabola)   .
   ------------------------------+------------------------------> tr A = tau
                               |
                 saddle (delta < 0): one positive and one negative real eigenvalue

   center: tau = 0, delta > 0 (purely imaginary eigenvalues, closed orbits)
   stable exactly when tau < 0 and delta > 0

Continuous versus discrete time.

	Discrete $\mathbf{x}_{k+1} = A\mathbf{x}_k$	Continuous $\mathbf{x}' = A\mathbf{x}$
mode	$\lambda^k\mathbf{v}$	$e^{\lambda t}\mathbf{v}$
stable mode	$\lvert\lambda\rvert < 1$	$\operatorname{Re}\lambda < 0$
neutral mode	$\lvert\lambda\rvert = 1$	$\operatorname{Re}\lambda = 0$
oscillation	$\lambda$ non-real (or negative)	$\lambda$ non-real
link	sampling every $h$ time units gives $A_{\text{discrete}} = e^{hA}$, eigenvalues $e^{h\lambda}$

import numpy as np
from scipy.linalg import expm
from scipy.integrate import solve_ivp

A = np.array([[-2.0, 1.0],
              [1.0, -2.0]])
lam, V = np.linalg.eigh(A)                       # symmetric: real eigenvalues, orthonormal eigenvectors
print(lam, np.round(V * np.sqrt(2), 4))          # [-3. -1.] [[ 1.  1.]  [-1.  1.]] (column signs may vary)

x0 = np.array([10.0, 2.0])                       # room 1 at +10, room 2 at +2 degrees
c = V.T @ x0                                     # eigen-coordinates (V is orthogonal)

def x_eigen(t):
    return V @ (np.exp(lam * t) * c)             # sum of c_i e^{lambda_i t} v_i

sol = solve_ivp(lambda t, x: A @ x, (0, 2), x0, rtol=1e-10, atol=1e-12, dense_output=True)
for t in [0.5, 2.0]:
    print(t, np.round(x_eigen(t), 4), np.allclose(x_eigen(t), expm(t * A) @ x0),
          np.allclose(x_eigen(t), sol.sol(t), atol=1e-7))
# 0.5 [4.5317 2.7467] True True
# 2.0 [0.8219 0.8021] True True     the rooms have nearly equalized

# Damped oscillator y'' + 2 zeta w y' + w^2 y = 0: eigenvalues for three damping ratios
w = 2.0
for zeta in [0.1, 1.0, 2.0]:
    M = np.array([[0.0, 1.0], [-w**2, -2 * zeta * w]])
    print(zeta, np.round(np.linalg.eigvals(M), 4))
# 0.1 [-0.2+1.99j -0.2-1.99j]       underdamped: decaying oscillation at frequency ~1.99
# 1.0 [-2.+0.j -2.+0.j]             critically damped: repeated, defective
# 2.0 [-0.5359+0.j -7.4641+0.j]     overdamped: two decaying exponentials

Key takeaways

• $e^{\lambda t}\mathbf{v}$ solves $\mathbf{x}' = A\mathbf{x}$ exactly when $(\lambda, \mathbf{v})$ is an eigenpair; an eigenbasis gives the general solution $\sum c_ie^{\lambda_it}\mathbf{v}_i$.
• Eigen-coordinates decouple the system into scalar equations $c_i' = \lambda_ic_i$.
• $\operatorname{Re}\lambda$ is the growth or decay rate and $\operatorname{Im}\lambda$ the angular frequency; stability means every $\operatorname{Re}\lambda < 0$.
• An $n$-th order linear ODE is a first-order system with a companion matrix; repeated roots correspond to Jordan blocks and $te^{\lambda t}$ terms.
• For $2 \times 2$ systems, the trace and determinant classify the equilibrium; sampling a continuous system gives $A_{\text{discrete}} = e^{hA}$.

🧪 Practice

1. Solve $\mathbf{x}' = \begin{bmatrix} 1 & 2 \ 2 & 1 \end{bmatrix}\mathbf{x}$ with $\mathbf{x}(0) = (3, 1)$, and classify the origin.
2. Find real solutions of $\mathbf{x}' = \begin{bmatrix} -1 & -4 \ 1 & -1 \end{bmatrix}\mathbf{x}$, giving the decay rate and the period of oscillation.
3. Write $y''' - 6y'' + 11y' - 6y = 0$ as a first-order system, find the eigenvalues of its matrix, and write the general solution for $y$.
4. Interview: A control engineer linearizes a drone's dynamics near hover and gets a $12 \times 12$ matrix with one eigenvalue $0.3 + 4i$. What does that mean physically, and what does a controller $\mathbf{u} = -K\mathbf{x}$ do in eigenvalue terms? Hint: an oscillation at about $4/(2\pi) \approx 0.64$ Hz whose amplitude grows like $e^{0.3t}$; feedback replaces $A$ by $A - BK$, and $K$ is chosen to move all eigenvalues into the left half-plane (pole placement).

Population and Growth Models

Why it matters. Many populations are structured: newborns cannot reproduce, adults reproduce at different rates by age, and survival depends on life stage. Ecologists, demographers, and fisheries managers use matrix models to predict long-run growth, to find the stable age structure, and, most importantly, to decide which life stage to protect for the greatest effect. All three answers come from the dominant eigenvalue and its eigenvectors.

The Leslie matrix. Divide females into age classes $0, 1, \dots, n - 1$ (each one time step wide). Let $f_i$ be the average number of female offspring per individual in class $i$, and $s_i$ the probability of surviving from class $i$ to class $i + 1$. Then $\mathbf{x}_{k+1} = L\mathbf{x}_k$ with

$$L = \begin{bmatrix} f_0 & f_1 & f_2 & \cdots & f_{n-1} \ s_0 & 0 & 0 & \cdots & 0 \ 0 & s_1 & 0 & \cdots & 0 \ \vdots & & \ddots & & \vdots \ 0 & 0 & \cdots & s_{n-2} & 0 \end{bmatrix}.$$

The first row creates newborns; the subdiagonal moves survivors up one age class.

The dominant eigenvalue is the growth rate. $L$ has nonnegative entries. If it is irreducible (every age class eventually contributes to every other) and primitive (some power is strictly positive; for example, two consecutive age classes reproduce), the Perron-Frobenius theorem (Section 13.2) guarantees:

1. a real eigenvalue $\lambda_1 > 0$ larger in modulus than every other eigenvalue;
2. a corresponding eigenvector $\mathbf{v}_1$ with all entries positive;
3. $\mathbf{x}_k \approx c\lambda_1^k\mathbf{v}_1$ for large $k$.

So $\lambda_1$ is the long-run growth factor per time step ($\lambda_1 > 1$ grows, $< 1$ declines, $= 1$ is stationary), and $\mathbf{v}_1$, scaled to sum to 1, is the stable age distribution: the long-run proportions in each class, regardless of the starting mix.

The characteristic equation. Expanding $\det(\lambda I - L)$ along the first row gives the Euler-Lotka equation

$$1 = \frac{f_0}{\lambda} + \frac{s_0f_1}{\lambda^2} + \frac{s_0s_1f_2}{\lambda^3} + \cdots + \frac{s_0\cdots s_{n-2}f_{n-1}}{\lambda^n}.$$

At $\lambda = 1$ the right side is $R_0 = f_0 + s_0f_1 + s_0s_1f_2 + \cdots$, the expected lifetime number of offspring per newborn. Since the right side decreases in $\lambda$, $\lambda_1 > 1$ exactly when $R_0 > 1$: a population grows when each individual more than replaces itself.

Sensitivity: which parameter matters most? Let $\mathbf{w}$ be the left eigenvector ($\mathbf{w}^TL = \lambda_1\mathbf{w}^T$, an eigenvector of $L^T$), which measures the reproductive value of each class. Then

$$\frac{\partial\lambda_1}{\partial\ell_{ij}} = \frac{w_iv_j}{\mathbf{w}^T\mathbf{v}}.$$

This formula holds for any simple eigenvalue of any matrix. Conservation biologists famously used it for loggerhead sea turtles: improving the survival of large juveniles (with devices that keep turtles out of shrimp nets) raised $\lambda_1$ far more than protecting eggs on beaches, which had been the focus of earlier efforts.

Example. A species with three age classes:

$$L = \begin{bmatrix} 0 & 1.5 & 1.2 \ 0.6 & 0 & 0 \ 0 & 0.5 & 0 \end{bmatrix}, \qquad R_0 = 0.6 \cdot 1.5 + 0.6 \cdot 0.5 \cdot 1.2 = 1.26 > 1.$$

The computation below gives $\lambda_1 \approx 1.107$: the population grows about 10.7% per step, with about 56% newborns, 30% one-year-olds, and 14% two-year-olds in the long run. The other two eigenvalues form a complex pair of modulus $0.57$, so the age structure settles with damped oscillations.

import numpy as np

L = np.array([[0.0, 1.5, 1.2],
              [0.6, 0.0, 0.0],
              [0.0, 0.5, 0.0]])

vals, V = np.linalg.eig(L)
i = np.argmax(np.abs(vals))
lam1 = vals[i].real
v = np.abs(V[:, i].real); v /= v.sum()                 # stable age distribution
print(round(lam1, 4), np.round(v, 4))                  # 1.1069 [0.5596 0.3033 0.137 ]
print(np.round(np.abs(vals), 4))                       # [1.1069 0.5703 0.5703]: the dominant one is unique

R0 = 1.5 * 0.6 + 1.2 * 0.6 * 0.5
print(round(R0, 4), lam1 > 1)                          # 1.26 True

# Simulation from a population of only newborns
x = np.array([100.0, 0.0, 0.0])
for k in range(1, 41):
    x_prev, x = x, L @ x
print(np.round(x.sum() / x_prev.sum(), 4), np.round(x / x.sum(), 4))   # 1.1069 [0.5596 0.3033 0.137 ]

# Sensitivity of lambda_1 to each entry: w_i v_j / (w . v)
wl, W = np.linalg.eig(L.T)
w = np.abs(W[:, np.argmax(np.abs(wl))].real)
S = np.outer(w, v) / (w @ v)
print(np.round(S, 3))
# [[0.441 0.239 0.108]
#  [0.814 0.441 0.199]
#  [0.479 0.259 0.117]]
print("survival s0:", round(S[1, 0], 3), " fecundity f1:", round(S[0, 1], 3))
# survival s0: 0.814  fecundity f1: 0.239

The sensitivity matrix shows that raising first-year survival $s_0$ (entry $(2, 1)$) by $0.01$ raises $\lambda_1$ by about $0.008$, more than three times the effect of raising the fecundity $f_1$ by the same amount. Entries of $S$ at positions where $L$ has a structural zero (for example, survival from class 0 straight to class 2) have no biological meaning and are ignored.

Key takeaways

• A Leslie matrix has fecundities in its first row and survival probabilities on its subdiagonal.
• For irreducible, primitive nonnegative matrices, the dominant eigenvalue $\lambda_1$ is real, positive, and strictly dominant (Perron-Frobenius).
• $\lambda_1$ is the long-run growth factor and its positive eigenvector is the stable age distribution.
• $\lambda_1 > 1$ exactly when the net reproductive number $R_0 > 1$.
• The sensitivity $\partial\lambda_1/\partial\ell_{ij} = w_iv_j/(\mathbf{w}^T\mathbf{v})$, from left and right eigenvectors, tells managers which life stage to protect.

🧪 Practice

1. For $L = \begin{bmatrix} 0 & 2 \ 0.5 & 0 \end{bmatrix}$, find the eigenvalues and explain why the population oscillates rather than settling into a stable age distribution. Which Perron-Frobenius hypothesis fails?
2. A two-class model has $L = \begin{bmatrix} 0.5 & 2 \ s & 0 \end{bmatrix}$. Find the value of $s$ for which the population is exactly stationary ($\lambda_1 = 1$).
3. Derive the Euler-Lotka equation for a $3 \times 3$ Leslie matrix by expanding $\det(\lambda I - L)$.
4. Interview: You are asked to choose between two interventions for an endangered species: increase egg survival by 50% or adult survival by 5%. How would you decide quantitatively? Hint: compare the changes in $\lambda_1$ using sensitivities, or better elasticities $\tfrac{\ell_{ij}}{\lambda_1}\tfrac{\partial\lambda_1}{\partial\ell_{ij}}$ (proportional changes), and check with a direct recomputation of $\lambda_1$ for each scenario.

PageRank and Network Centrality

Why it matters. Which web page is most important? Which person in a social network is most influential? Which protein is central to a cell's interaction network? Counting links is a poor answer, because a link from an important node should count more than a link from an obscure one. That definition is circular: importance depends on importance. Eigenvectors are exactly the solutions of such circular definitions, and Google's original PageRank algorithm made one of them the basis of web search.

Eigenvector centrality. For an undirected graph with adjacency matrix $A$, define a score $x_i$ for each node as proportional to the sum of its neighbors' scores:

$$x_i = \frac{1}{\lambda}\sum_ja_{ij}x_j \iff A\mathbf{x} = \lambda\mathbf{x}.$$

Requiring all scores to be positive selects the dominant eigenvector (Perron-Frobenius, for a connected graph). This is eigenvector centrality.

PageRank: the random surfer. On the directed web graph, imagine a surfer who, at each step, clicks a uniformly random outgoing link of the current page. The importance of a page is the long-run fraction of time the surfer spends there: the steady state of a Markov chain (previous topics). Let $n$ be the number of pages and define the column-stochastic link matrix

$$H_{ij} = \begin{cases} 1/\operatorname{outdeg}(j) & \text{if page } j \text{ links to page } i, \ 0 & \text{otherwise.} \end{cases}$$

Two problems prevent this chain from having a unique steady state:

1. Dangling nodes (pages with no outgoing links, such as PDFs) give zero columns. Fix: replace each such column by $\tfrac1n\mathbf{1}$ (the surfer jumps to a random page). Call the result $S$.
2. Disconnected or trapping regions make the chain reducible or periodic. Fix: with probability $1 - \alpha$, the surfer teleports to a uniformly random page.

The Google matrix is

$$G = \alpha S + \frac{1 - \alpha}{n}\mathbf{1}\mathbf{1}^T, \qquad \alpha \approx 0.85.$$

Every entry of $G$ is positive, so $G$ is regular, and the PageRank vector is its unique steady state: $G\mathbf{r} = \mathbf{r}$, $\mathbf{r} > 0$, $\sum r_i = 1$.

Computing it: the power method. For billions of pages, $G$ is never formed (it is dense). Instead iterate

$$\mathbf{r}_{k+1} = \alpha S\mathbf{r}_k + \frac{1 - \alpha}{n}\mathbf{1},$$

which uses only the sparse link structure ($\mathbf{1}^T\mathbf{r}_k = 1$ makes this equal to $G\mathbf{r}_k$). Every eigenvalue of $G$ other than 1 satisfies $\lvert\lambda\rvert \le \alpha$, so the error shrinks at least like $0.85^k$: about 50 iterations give 3 to 4 correct digits, independent of the size of the web. The damping factor is a trade-off: $\alpha$ closer to 1 follows the link structure more faithfully but converges more slowly.

Example. A five-page web:

        +-----+          +-----+
        |  1  | -------> |  2  |
        +-----+          +-----+
         ^   \            |   ^
         |    \           |   |
         |     v          v   |
        +-----+          +-----+          +-----+
        |  4  | <------- |  3  | -------> |  5  |   page 5 links nowhere (dangling)
        +-----+          +-----+          +-----+

   links: 1 -> 2, 1 -> 3, 2 -> 3, 3 -> 2, 3 -> 4, 3 -> 5, 4 -> 1

A family of centrality measures.

Measure	Definition	Eigen-connection
degree centrality	number of links	none: ignores who links to you
eigenvector centrality	$A\mathbf{x} = \lambda_1\mathbf{x}$	dominant eigenvector of the adjacency matrix
Katz centrality	$\mathbf{x} = \alpha A\mathbf{x} + \mathbf{1}$	$(I - \alpha A)^{-1}\mathbf{1}$; needs $\alpha < 1/\rho(A)$
PageRank	$G\mathbf{r} = \mathbf{r}$	steady state of the random-surfer chain
HITS hubs and authorities	$\mathbf{h} \propto A\mathbf{a}$, $\mathbf{a} \propto A^T\mathbf{h}$	dominant eigenvectors of $AA^T$ and $A^TA$ (singular vectors, Chapter 11)

import numpy as np

n, alpha = 5, 0.85
links = [(1, 2), (1, 3), (2, 3), (3, 2), (3, 4), (3, 5), (4, 1)]   # (from, to), pages 1..5

H = np.zeros((n, n))
for j, i in links:
    H[i - 1, j - 1] = 1.0                                # column j = outgoing links of page j
out = H.sum(axis=0)
S = H.copy()
S[:, out > 0] /= out[out > 0]                            # each column sums to 1
S[:, out == 0] = 1.0 / n                                 # dangling page 5: jump anywhere
G = alpha * S + (1 - alpha) / n * np.ones((n, n))        # dense Google matrix (small n only)

# Power method using only S (what a search engine would do)
r = np.full(n, 1.0 / n)
for k in range(100):
    r_new = alpha * S @ r + (1 - alpha) / n
    if np.abs(r_new - r).sum() < 1e-12:
        break
    r = r_new
print(k, np.round(r_new, 4))                             # 45 [0.1769 0.2192 0.316  0.144  0.144 ]

# Same answer from the eigenvector for eigenvalue 1
vals, vecs = np.linalg.eig(G)
q = vecs[:, np.argmin(np.abs(vals - 1))].real
print(np.allclose(q / q.sum(), r_new))                   # True
print(np.round(sorted(np.abs(vals))[-2], 4) <= alpha)    # True: |lambda_2| <= alpha

print("ranking:", [int(p) for p in np.argsort(-r_new, kind="stable") + 1])   # ranking: [3, 2, 1, 4, 5]

Page 3 ranks first: it receives links from both 1 and 2, and page 2 passes on its entire vote. Page 2 is second because one of its incoming links comes from page 3, the most important page. Page 1 has a single incoming link, but page 4 gives it its whole vote, so page 1 beats page 4. Pages 4 and 5 tie exactly: each receives one third of page 3's vote and nothing else.

Key takeaways

• Eigenvector centrality and PageRank resolve the circular definition "important nodes are linked from important nodes" with a dominant eigenvector.
• PageRank is the steady state of a random surfer who follows links with probability $\alpha$ and teleports otherwise.
• Dangling-node fixes and teleportation make the Google matrix positive, so Perron-Frobenius guarantees a unique positive ranking.
• The power method converges at rate $\alpha^k$ using only sparse matrix-vector products, regardless of the graph size.
• Katz, HITS, and PageRank are variations of the same eigenvector idea; HITS leads to singular vectors (Chapter 11).

🧪 Practice

1. For the five-page example, compute PageRank with $\alpha = 0.5$ and $\alpha = 0.99$ and compare the rankings. Which pages gain or lose?
2. Compute the eigenvector centrality of the four-node graph with edges $1\text{-}2$, $1\text{-}3$, $2\text{-}3$, $3\text{-}4$ (Section 4.1), and compare it with degree centrality.
3. Show that if $S$ is column-stochastic and $\mathbf{1}^T\mathbf{r} = 1$, then $G\mathbf{r} = \alpha S\mathbf{r} + \tfrac{1 - \alpha}{n}\mathbf{1}$, so the dense matrix $G$ never needs to be formed.
4. Interview: How would you compute PageRank for a graph with 10 billion nodes and 100 billion edges that does not fit on one machine, and how many iterations would you budget? Hint: store the sparse $S$ distributed by columns or edges, do one distributed sparse matrix-vector product per iteration (MapReduce or Pregel style), handle dangling mass as a single scalar, and budget $\log(\varepsilon)/\log(\alpha)$ iterations, about 50 to 100 for $\alpha = 0.85$.

<a id="9-inner-product-spaces-and-orthogonality"></a>