Chapter 8 The Eigenvalue Problem

Counting the ways to combine long and short syllables and tracking how rabbits breed seem to have nothing to do with each other, yet both lead to recursively defined sequences. The tradition of Indian prosody and Fibonacci’s Liber Abaci of 1202 are different threads in the history of such problems. For the Fibonacci sequence, the ratio of consecutive terms approaches the golden ratio $\varphi=(1+\sqrt5)/2$ ; Binet’s formula, in turn, writes a sequence of integers as a combination of powers of irrational numbers. What structure is hidden behind all this?

Modern linear algebra lays a remarkably clean card on the table: the golden ratio $\varphi$ is, in essence, the dominant eigenvalue of a $2 \times 2$ state transition matrix. Once you have mastered the geometric viewpoint of diagonalizing a matrix, the square root of five and the irrational powers in Binet’s formula reveal themselves as naturally as a finely made gift being unwrapped—a number-theoretic coincidence spanning two thousand years turns out, in a linear space, to be nothing more than a completely standard eigenvalue computation.

This is precisely the ultimate lens that this chapter sets out to forge for you.

If a nonzero vector $\mathbf{v}$ satisfies $\mathbf{A}\mathbf{v}=\lambda\mathbf{v}$ , then the one-dimensional subspace it spans is left invariant by the transformation; $\mathbf{v}$ is called an eigenvector and $\lambda$ an eigenvalue. In the real case this can be pictured as stretching, compressing, reversing, or collapsing to zero along a single line. But a real matrix need not have any real eigendirection—a rotation of the plane by ninety degrees is an example—while the identity matrix has every nonzero vector as an eigenvector. Over the complex numbers every nonempty square matrix has at least one eigenvalue, yet it need not have enough linearly independent eigenvectors to be diagonalized.

The eigen in the German word Eigenwert means “own” or “proper.” Eigenvalues are unchanged by a similarity change of basis, which expresses an intrinsic property of the operator; but a list of eigenvalues alone is usually not enough to reconstruct the whole operator.

This language shows a remarkable universality across modern science and engineering:

The quantum world. The bound-state eigenvalues of the Hamiltonian operator give the discrete energy levels; the frequency of a spectral line is determined by an energy difference $h\nu=E_i-E_j$ and is not itself an eigenvalue.
Information networks. PageRank models the web with a damped random transition matrix, and the weights of the web pages are given by the nonnegative stationary eigenvector for the eigenvalue 1 whose components sum to 1.
High-dimensional data. Principal component analysis (PCA) in machine learning takes the covariance matrix of high-dimensional data and uses the eigenvector belonging to the largest eigenvalue to lock onto the principal axis of greatest variance (the eigenvalue is a scalar and only gives the amount of variance in that direction); it is a cornerstone of modern dimensionality reduction.

Chapter 9 will reveal why the eigenvalues of symmetric systems must be purely real; Chapter 10 will explain why the observables of quantum mechanics are described by self-adjoint operators, with the measurement rules given jointly by the eigenvalues and the corresponding spectral projections; and the singular value decomposition (SVD) of Chapter 11 will carry this beautiful idea of the eigen-spectrum all the way to matrices of arbitrary shape and to finite-dimensional linear maps.

Different worlds, read through the same eigenspace.

Chapter Structure and Learning Objectives¶

The story of the Fibonacci sequence pays off in full in §8.3: how the recurrence is written in matrix form, how the eigenvalues emerge as the golden ratio, and how diagonalization yields Binet’s formula. But to follow that story, we first need to lay the foundation in the first two sections of this chapter.

§8.1 begins by stating the central problem clearly. How does the geometric intuition of an “unchanged direction” become an algebraic condition that can be computed? The answer lies in the characteristic polynomial, in the equation $\det(\mathbf{A} - \lambda\mathbf{I}) = 0$ —finding eigenvalues is equivalent to finding the roots of a polynomial, and the fundamental theorem of algebra guarantees that the roots exist over the complex numbers. Here geometric intuition and algebraic tools shake hands for the first time.

§8.2 goes deeper into the structure of the characteristic polynomial. An eigenvalue can be a repeated root, and the gap between algebraic multiplicity and geometric multiplicity is precisely the watershed that decides whether a matrix can be “fully understood.” This section also introduces the Cayley–Hamilton theorem—every matrix satisfies its own characteristic equation. At first glance it looks like magic; on reflection it is nothing more than an inevitable consequence of linear dependence.

§8.3 reaches the central goal: change the basis so that the matrix becomes as simple as possible. A matrix being diagonalizable means that, in some basis, the linear transformation is nothing but independent scaling along each direction. This section answers “when can this be done” and “what can we do once it is done”—the closed-form formula for the Fibonacci sequence is the most concise demonstration of the power of diagonalization.

§8.4 faces an unavoidable reality: not every matrix can be diagonalized. The Jordan canonical form gives the most general answer—even when a matrix cannot be fully diagonalized, it can be brought to an almost diagonal block structure. This is one of the deepest theorems of finite-dimensional linear algebra, and a stage that cannot be bypassed in understanding the inner structure of matrices.

§8.5 verifies the conclusions of this chapter in Python and implements numerical methods such as power iteration, so that you can feel the distance between “diagonalizing by hand” and “efficient solution by computer”—both near and far.

Once you have read this chapter, you will have a far deeper answer than before to the question “what, really, is a matrix?”; and the sequence that left traces both in the Indian tradition of prosody and in Liber Abaci, and whose structure linear algebra finally brought into view, will become your most concrete memory of the power of eigenvalues.

8.1 Fundamentals of Eigenvalues and Eigenvectors¶

At its core, the eigenvalue problem is a search for the invariant subspaces and invariant directions of a matrix. When a matrix is viewed as a linear transformation, eigenvalues and eigenvectors give us the simplest way to represent the mapping and reveal the intrinsic structure and properties of the matrix. This section starts from the definition, introduces the basic concepts of eigenvalues and eigenvectors, and explores how to compute them and what they mean geometrically.

8.1.1 The Characteristic Equation and the Characteristic Polynomial¶

The concepts of eigenvalue and eigenvector grew out of a close study of the basic behavior of linear transformations. When a linear transformation acts on certain special vectors, the direction of these vectors is preserved and only their magnitude is scaled.

This definition has three key ingredients:

Square matrices only. The eigenvalue problem applies only to square matrices, that is, matrices whose number of rows equals their number of columns.
Nonzero vectors. An eigenvector must be a nonzero vector; otherwise the equation would hold automatically for every $\lambda$ .
Invariance of direction. After the transformation, an eigenvector still lies on its original invariant line, and its length is multiplied by $|\lambda|$ ; a negative value reverses it, and the value zero sends it to zero (here we mean real eigenvalues).

To solve for eigenvalues and eigenvectors, we need to turn the defining equation of an eigenvalue into a form that is easier to work with.

The definition of an eigenvalue can be rewritten as: $\mathbf{A}\mathbf{v} = \lambda\mathbf{v}$ , $\mathbf{A}\mathbf{v} - \lambda\mathbf{v} = \mathbf{0}$ , $(\mathbf{A} - \lambda\mathbf{I})\mathbf{v} = \mathbf{0}$

where $\mathbf{I}$ is the $n \times n$ identity matrix. Since we require the eigenvector $\mathbf{v}$ to be nonzero, $\mathbf{A} - \lambda\mathbf{I}$ must be a singular matrix, that is, its determinant must be zero:

\det(\mathbf{A} - \lambda\mathbf{I}) = 0

(2)

This equation is called the characteristic equation of the matrix $\mathbf{A}$ . When the determinant on the left is expanded, it is a polynomial in $\lambda$ , called the characteristic polynomial of the matrix $\mathbf{A}$ .

Historical significance and applications.

This is one of the major achievements of the mathematician Gauss in the early nineteenth century, and a central theorem of algebra. It tells us that the field of complex numbers is algebraically closed: every polynomial equation in one variable with complex coefficients can be solved within the complex numbers.

It is a concise existence theorem; in practice, other constructive methods (such as Newton’s method) are still needed to actually compute the roots, and the relevant Python packages can usually find approximate solutions that meet the required precision.

A direct corollary is that a polynomial of degree $n$ has exactly $n$ roots in the complex numbers (counted with multiplicity), so an $n \times n$ matrix has exactly $n$ eigenvalues (possibly repeated). These eigenvalues may be real or complex, even when all the entries of the matrix are real.

Proof 1 (A Proof by Real Analysis)

We use proof by contradiction, combined with the extreme-value property of continuous functions.

Without loss of generality, let

f(z) = z^n + a_{n-1} z^{n-1} + \cdots + a_1 z + a_0

(4)

(every term has been divided by $a_n$ ). Assume that $f(z)$ has no complex root, that is, $f(z) \ne 0$ for all $z \in \mathbb{C}$ .

Consider the function $g(z) := |f(z)|^2$ . Since $f(z)$ is a polynomial, $g(z)$ is a continuous function on $\mathbb{R}^2$ (regarding $\mathbb{C}$ as $\mathbb{R}^2$ ).
As $|z| \to \infty$ ,

|f(z)| = |z|^n \left|1 + \frac{a_{n-1}}{z} + \cdots + \frac{a_0}{z^n}\right| \to \infty

(5)

so the minimum of $g$ is not attained at infinity. There exists a sufficiently large $R > 0$ such that $g(z) > g(0)$ for all $|z| \ge R$ .

On the closed disk $\{z \in \mathbb{C} : |z| \le R\}$ , the continuous function $g$ must attain a minimum by the extreme value theorem; say it is attained at $z_0$ , so that $g(z_0) \le g(0)$ . Combined with the behavior on the boundary, $z_0$ is a global minimum point.
Let $h(w)=f(z_0+w)/f(z_0)$ . Since $f(z_0)\neq0$ by assumption, we can write
$h(w)=1+b_kw^k+\cdots+b_nw^n,$
(6)
where $1\le k\le n$ is the degree of the first nonzero nonconstant term and $b_k\neq0$ ; we do not assume that the leading coefficient is still 1.
Choose $w=re^{i\theta}$ with $e^{ik\theta}=-|b_k|/b_k$ . Then
$h(re^{i\theta})=1-|b_k|r^k+O(r^{k+1}),\qquad |h(re^{i\theta})|^2=1-2|b_k|r^k+O(r^{k+1})<1$
(7)
holds for all sufficiently small $r>0$ ; when $k=n$ the higher-order remainder is zero and the inequality still holds. Hence $|f(z_0+w)|<|f(z_0)|$ , contradicting the global minimality.

This contradiction shows that the assumption is false, so

\underline{f(z) \text{ must have a complex root.}}

(8)

Remark. The key to this proof is the extreme-value property of continuous functions on compact sets (a central result of real analysis), together with the asymptotic behavior of polynomials at infinity. The proof deliberately avoids tools of complex analysis such as contour integrals, but it still needs the basic properties of complex numbers. The fundamental theorem of algebra has many proofs, yet every one of them ultimately relies on properties from analysis or topology.

Standard procedure for computing eigenvalues and eigenvectors.

Compute the characteristic polynomial $p_A(\lambda) = \det(\mathbf{A} - \lambda\mathbf{I})$ .
Solve the characteristic equation $p_A(\lambda) = 0$ to obtain the eigenvalues.
For each eigenvalue $\lambda_i$ , solve the homogeneous linear system $(\mathbf{A} - \lambda_i\mathbf{I})\mathbf{v} = \mathbf{0}$ to obtain the corresponding eigenvectors.

Example 1

Compute the eigenvalues and eigenvectors of the matrix $\mathbf{A} = \begin{pmatrix} 3 & 1 \\ 1 & 3 \end{pmatrix}$ .

Compute the characteristic polynomial

\begin{align} p_A(\lambda) &= \det(\mathbf{A} - \lambda\mathbf{I}) \\ &= \det\begin{pmatrix} 3-\lambda & 1 \\ 1 & 3-\lambda \end{pmatrix} \\ &= (3-\lambda)(3-\lambda) - 1 \cdot 1 \\ &= (3-\lambda)^2 - 1 \\ &= \lambda^2 - 6\lambda + 9 - 1 \\ &= \lambda^2 - 6\lambda + 8 \end{align}

(9)

Solve the characteristic equation
$\lambda^2 - 6\lambda + 8 = 0$
(10)

Using factorization or the quadratic formula:

(\lambda - 2)(\lambda - 4) = 0

(11)

we obtain the eigenvalues $\lambda_0 = 2$ and $\lambda_1 = 4$ .

Compute the eigenvectors

For $\lambda_0 = 2$ , solve $(\mathbf{A} - 2\mathbf{I})\mathbf{v} = \mathbf{0}$ :

\begin{pmatrix} 1 & 1 \\ 1 & 1 \end{pmatrix} \begin{pmatrix} v_0 \\ v_1 \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \end{pmatrix}

(12)

This gives the equation $v_0 + v_1 = 0$ , that is, $v_0 = -v_1$ . Choosing $v_1 = 1$ gives $v_0 = -1$ . The eigenvector is $\mathbf{v}_0 = \begin{pmatrix} -1 \\ 1 \end{pmatrix}$ (or any nonzero multiple of it).

For $\lambda_1 = 4$ , solve $(\mathbf{A} - 4\mathbf{I})\mathbf{v} = \mathbf{0}$ :

\begin{pmatrix} -1 & 1 \\ 1 & -1 \end{pmatrix} \begin{pmatrix} v_0 \\ v_1 \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \end{pmatrix}

(13)

This gives the equation $-v_0 + v_1 = 0$ , that is, $v_0 = v_1$ . Choosing $v_0 = 1$ gives $v_1 = 1$ . The eigenvector is $\mathbf{v}_1 = \begin{pmatrix} 1 \\ 1 \end{pmatrix}$ (or any nonzero multiple of it).

Geometric interpretation.

For the matrix $\mathbf{A} = \begin{pmatrix} 3 & 1 \\ 1 & 3 \end{pmatrix}$ :

The eigenvector $\mathbf{v}_0 = \begin{pmatrix} -1 \\ 1 \end{pmatrix}$ points in the “northwest–southeast” direction, and the transformation stretches along this direction by a factor of 2.
The eigenvector $\mathbf{v}_1 = \begin{pmatrix} 1 \\ 1 \end{pmatrix}$ points in the “northeast” direction, and the transformation stretches along this direction by a factor of 4.

Note that these two eigenvectors are perpendicular to each other ( $\mathbf{v}_0 \cdot \mathbf{v}_1 = 0$ ); this is because $\mathbf{A}$ is a symmetric matrix (we discuss this property in detail in Chapter 9).

In the coordinate system formed by these eigenvectors, the action of the linear transformation $\mathbf{A}$ becomes extremely simple: it merely scales along two perpendicular directions separately, with no rotation or shear at all.

import sympy as sp

# --- Visual hierarchy helpers ---
def print_header(title):
    print("=" * 60)
    print(f"  {title}")
    print("=" * 60)

def print_step(step, desc):
    print(f"\n▶ Step {step}: {desc}")
    print("-" * 40)

# ============================================================
print_header("Example 8.1 | Symbolic Computation of the Characteristic Polynomial")
# ============================================================

# --- Step 1: Define the matrix A ---
print_step("1", "Define the matrix A")
A = sp.Matrix([[3, 1],
               [1, 3]])
sp.pprint(A)

# --- Step 2: Compute the characteristic polynomial p_A(λ) = det(A - λI) (consistent with def-characteristic-polynomial) ---
# Note: numerical packages often use the monic convention det(λI - A); the two agree for even n and differ by a sign (-1)^n for odd n,
#       while the set of roots (the eigenvalues) stays the same.
print_step("2", "Compute the characteristic polynomial p_A(λ) = det(A - λI)")
lam = sp.Symbol('lambda')
char_poly = (A - lam * sp.eye(2)).det()
print(f"  Expanded: p_A(λ) = {sp.expand(char_poly)}")
print(f"  Factored: p_A(λ) = {sp.factor(char_poly)}")

# --- Step 3: Find the eigenvalues ---
print_step("3", "Solve p_A(λ) = 0")
eigenvalues = sp.solve(char_poly, lam)
print(f"  Eigenvalues: λ₀ = {eigenvalues[0]},  λ₁ = {eigenvalues[1]}")

============================================================
  Example 8.1 | Symbolic Computation of the Characteristic Polynomial
============================================================

▶ Step 1: Define the matrix A
----------------------------------------
⎡3  1⎤
⎢    ⎥
⎣1  3⎦

▶ Step 2: Compute the characteristic polynomial p_A(λ) = det(A - λI)
----------------------------------------
  Expanded: p_A(λ) = lambda**2 - 6*lambda + 8
  Factored: p_A(λ) = (lambda - 4)*(lambda - 2)

▶ Step 3: Solve p_A(λ) = 0
----------------------------------------
  Eigenvalues: λ₀ = 2,  λ₁ = 4

8.1.2 Geometric Interpretation of Eigenvalues and Eigenvectors¶

Eigenvalues and eigenvectors carry deep geometric meaning: they reveal the essential behavior of a linear transformation. To search for invariant vectors yourself, see Experiment 1.

Basic geometric interpretation.

Eigenvectors represent the invariant directions of a linear transformation.
Eigenvalues represent the stretching or compression ratios along these directions.

When a matrix $\mathbf{A}$ acts on a vector space as a linear transformation, most vectors are not only scaled but also change direction. Eigenvectors are the exception—they are only stretched or compressed along their original direction, and the scaling ratio is exactly the corresponding eigenvalue. An interesting and important question: does a rotation matrix still have invariant eigenvectors?

Depending on the eigenvalue, we can observe several typical situations:

When $\lambda > 1$ , the eigenvector is stretched.
When $0 < \lambda < 1$ , the eigenvector is compressed.
When $\lambda < 0$ , the eigenvector’s direction is reversed and it is scaled by $|\lambda|$ .
When $\lambda = 0$ , the eigenvector is mapped to the zero vector, indicating that these directions “vanish” under the transformation.
When $\lambda = 1$ , the eigenvector is left completely unchanged.

As shown in Chapter 3, several typical transformations have the following eigenvalues and eigenvectors:

The identity matrix $\mathbf{I}$ :
- Every nonzero vector is an eigenvector, and the corresponding eigenvalue is always 1.
- Geometric interpretation: the identity transformation, which leaves every vector unchanged.
A projection matrix $\mathbf{P}$ (projecting onto some subspace):
- The only eigenvalues are 0 and 1.
- The eigenvectors for the eigenvalue 1 lie in the subspace and are left unchanged.
- The eigenvectors for the eigenvalue 0 are perpendicular to the subspace and are mapped to zero.
A two-dimensional rotation matrix $\mathbf{R}$ :
- If the angle of rotation is not an integer multiple of $\pi$ , there are no real eigenvalues.
- The eigenvalues are complex, of the form $e^{±i\theta}$ ( $\theta$ is the angle of rotation).
- Geometric interpretation: under this condition on the angle there is no real invariant line.
- This also shows that complex matrices arise naturally.

Example 2 (Eigenvalues and Eigenvectors of Rotation Matrices)

Consider the two-dimensional rotation matrix

\mathbf{R}_{\theta} = \begin{pmatrix} \cos\theta & -\sin\theta \\ \sin\theta & \cos\theta \end{pmatrix}

(14)

Case 1: $\theta = 0$ (no rotation)

\mathbf{R}_0 = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} = \mathbf{I}

(15)

Characteristic polynomial:

p(\lambda) = \det(\mathbf{R}_0 - \lambda\mathbf{I}) = (1-\lambda)^2 = 0 \implies \lambda = 1 \text{ (multiplicity 2)}

(16)

Every nonzero vector $\mathbf{v} \in \mathbb{R}^2$ satisfies $\mathbf{I}\mathbf{v} = 1 \cdot \mathbf{v}$ , so every nonzero vector is an eigenvector, and the only eigenvalue is $\lambda = 1$ .

Case 2: $\theta = \pi$ (rotation by 180°)

\mathbf{R}_\pi = \begin{pmatrix} -1 & 0 \\ 0 & -1 \end{pmatrix} = -\mathbf{I}

(17)

Characteristic polynomial:

p(\lambda) = \det(-\mathbf{I} - \lambda\mathbf{I}) = (-1-\lambda)^2 = 0 \implies \lambda = -1 \text{ (multiplicity 2)}

(18)

Every nonzero vector satisfies $(-\mathbf{I})\mathbf{v} = -1 \cdot \mathbf{v}$ , so every nonzero vector is an eigenvector, and the only eigenvalue is $\lambda = -1$ .

Case 3: $0 < \theta < \pi$ (a general rotation)

Step 1. Compute the characteristic polynomial.

p(\lambda) = \det(\mathbf{R}_\theta - \lambda\mathbf{I}) = \det\begin{pmatrix} \cos\theta - \lambda & -\sin\theta \\ \sin\theta & \cos\theta - \lambda \end{pmatrix}

(19)

= (\cos\theta - \lambda)^2 + \sin^2\theta = \lambda^2 - 2\cos\theta\,\lambda + (\cos^2\theta + \sin^2\theta) = \lambda^2 - 2\cos\theta\,\lambda + 1

(20)

Step 2. Solve the characteristic equation $\lambda^2 - 2\cos\theta\,\lambda + 1 = 0$ .

\lambda = \frac{2\cos\theta \pm \sqrt{4\cos^2\theta - 4}}{2} = \cos\theta \pm i\sin\theta

(21)

Since $\sin\theta > 0$ for $0 < \theta < \pi$ , the discriminant $4\cos^2\theta - 4 < 0$ , so

\lambda_0 = e^{i\theta} = \cos\theta + i\sin\theta, \qquad \lambda_1 = e^{-i\theta} = \cos\theta - i\sin\theta

(22)

There is no solution in the real numbers; in the complex numbers, the two eigenvalues are complex conjugates of each other, and both have modulus 1 (reflecting the fact that rotations preserve length).

Step 3. Find the eigenvector for $\lambda_0 = e^{i\theta}$ by solving $(\mathbf{R}_\theta - \lambda_0 \mathbf{I})\mathbf{v} = \mathbf{0}$ :

\begin{pmatrix} \cos\theta - e^{i\theta} & -\sin\theta \\ \sin\theta & \cos\theta - e^{i\theta} \end{pmatrix} \begin{pmatrix} v_0 \\ v_1 \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \end{pmatrix}

(23)

Note that $\cos\theta - e^{i\theta} = -i\sin\theta$ , so the equation from row 0 is

-i\sin\theta\, v_0 - \sin\theta\, v_1 = 0 \implies v_1 = -i\, v_0

(24)

Taking $v_0 = 1$ gives

\mathbf{v}_0 = \begin{pmatrix} 1 \\ -i \end{pmatrix}

(25)

Step 4. Find the eigenvector for $\lambda_1 = e^{-i\theta}$ . Similarly, $\cos\theta - e^{-i\theta} = i\sin\theta$ , and the equation from row 0 is

i\sin\theta\, v_0 - \sin\theta\, v_1 = 0 \implies v_1 = i\, v_0

(26)

Taking $v_0 = 1$ gives

\mathbf{v}_1 = \begin{pmatrix} 1 \\ i \end{pmatrix}

(27)

Step 5. Geometric interpretation.

Note that $\mathbf{v}_0$ and $\mathbf{v}_1$ form a complex-conjugate pair: for a real matrix, complex eigenvalues and eigenvectors always occur in conjugate pairs. No real vector keeps its direction under the rotation precisely because a pure rotation ( $0 < \theta < \pi$ ) changes the direction of every real vector; the complex eigenvectors $\begin{pmatrix}1\\-i\end{pmatrix}$ and $\begin{pmatrix}1\\i\end{pmatrix}$ encode the periodicity and the phase information of the rotation.

Through geometric intuition we have seen that eigenvalues and eigenvectors reveal the “principal axes” and the “scaling factors” of a linear transformation. But this is only the tip of the iceberg. In Section 8.2 we explore the deeper theoretical properties of eigenvalues, including eigenspaces and the notions of algebraic and geometric multiplicity, which provide the theoretical foundation for understanding the complete structure of a matrix.

import numpy as np

def print_header(title):
    print("=" * 60)
    print(f"  {title}")
    print("=" * 60)

def print_step(step, desc):
    print(f"\n▶ Step {step}: {desc}")
    print("-" * 40)

np.set_printoptions(precision=6, suppress=True, linewidth=100)

# ============================================================
print_header("Verifying the Eigenvalues of Rotation Matrices")
# ============================================================

for theta_label, theta in [("θ = 0", 0.0),
                             ("θ = π", np.pi),
                             ("θ = π/3", np.pi / 3)]:
    print(f"\n{'━'*50}")
    print(f"  {theta_label}")
    print(f"{'━'*50}")

    # Step 0: Build the rotation matrix
    print_step(0, "Build the rotation matrix R_θ")
    R = np.array([[np.cos(theta), -np.sin(theta)],
                  [np.sin(theta),  np.cos(theta)]])
    print(R)

    # Step 1: Compute the eigenvalues and eigenvectors (over the complex numbers)
    print_step(1, "NumPy eigendecomposition (over the complex numbers)")
    eigvals, eigvecs = np.linalg.eig(R)
    for k in range(2):
        lam = eigvals[k]
        v   = eigvecs[:, k]
        # Verify R @ v = λ * v
        lhs = R @ v
        rhs = lam * v
        residual = np.max(np.abs(lhs - rhs))
        print(f"  λ_{k} = {lam:.6f}  |λ_{k}| = {abs(lam):.6f}")
        print(f"  v_{k} = {v}")
        print(f"  ‖R·v_{k} − λ_{k}·v_{k}‖∞ = {residual:.2e}")

    # Step 2: Verify |λ| = 1
    print_step(2, "Verify that every eigenvalue has modulus 1 (rotations preserve length)")
    print(f"  |λ_0| = {abs(eigvals[0]):.8f}")
    print(f"  |λ_1| = {abs(eigvals[1]):.8f}")

    # Step 3: Verify the trace and the determinant
    print_step(3, "Verify tr = 2cosθ, det = 1")
    print(f"  tr(R)   = {np.trace(R):.6f},  2cos(θ) = {2*np.cos(theta):.6f}")
    print(f"  det(R)  = {np.linalg.det(R):.6f},  expected = 1.000000")

§8.1 Summary

This section answered a most basic question: does a linear transformation have vectors whose “direction does not change”?

The answer lies among the zeros of the characteristic polynomial $p_{\mathbf{A}}(\lambda) = \det(\mathbf{A} - \lambda\mathbf{I})$ . The fundamental theorem of algebra guarantees that a polynomial of degree $n$ has exactly $n$ roots in the complex numbers, so an $n \times n$ matrix must have $n$ eigenvalues (counted with multiplicity)—even when all the entries of the matrix are real, its eigenvalues may be complex. Rotation matrices are the most intuitive example: a two-dimensional rotation through an angle that is not an integer multiple of $\pi$ has no real invariant line, and this is reflected precisely in the complex nature of its eigenvalues.

Geometrically, an eigenvalue describes the scaling ratio along an eigendirection: $|\lambda|>1$ means stretching, $0<|\lambda|<1$ means compression, a negative value comes with a reversal of direction, and zero means that the direction is “flattened.” In the theory of diagonalization in §8.3, these cases become a precise classification of matrix structure.

Computationally, the standard procedure is: find the characteristic polynomial → find its roots to obtain the eigenvalues → for each eigenvalue, solve for the null space to obtain the eigenvectors. The procedure looks simple, but in higher dimensions it hides a subtle gap between algebraic multiplicity and geometric multiplicity—and that will be the central topic of §8.2.

8.2 Theoretical Properties of Eigenvalues¶

This section explores the important theoretical properties of eigenvalues in depth, relates eigenvalues to other properties of a matrix, and introduces the concept of an eigenspace. These theoretical foundations provide a solid framework for understanding matrix structure and for applying eigenvalue analysis.

8.2.1 Relations Between Eigenvalues and Matrix Properties¶

Eigenvalues are closely connected with several basic properties of a matrix, and these connections reveal the essential features of the matrix’s intrinsic structure.

Exercise 1 (Basic Relations Between Eigenvalues and Matrix Properties)

For an $n \times n$ matrix $\mathbf{A}$ with eigenvalues $\lambda_0, \lambda_1, \ldots, \lambda_{n-1}$ (repeated values included), prove the following relations:

Trace and the sum of the eigenvalues. $\text{tr}(\mathbf{A}) = \sum_{i=0}^{n-1} \lambda_i = \lambda_0 + \lambda_1 + \cdots + \lambda_{n-1}$
Determinant and the product of the eigenvalues. $\det(\mathbf{A}) = \prod_{i=0}^{n-1} \lambda_i = \lambda_0 \times \lambda_1 \times \cdots \times \lambda_{n-1}$
Powers. If $\lambda$ is an eigenvalue of $\mathbf{A}$ , then $\lambda^k$ is an eigenvalue of $\mathbf{A}^k$ .
Inverse. If $\lambda \neq 0$ is an eigenvalue of $\mathbf{A}$ , then $\frac{1}{\lambda}$ is an eigenvalue of $\mathbf{A}^{-1}$ (if the inverse exists).
Shift. If $\lambda$ is an eigenvalue of $\mathbf{A}$ , then $\lambda + c$ is an eigenvalue of $\mathbf{A} + c\mathbf{I}$ .

Solution to Exercise 1

1. Trace and the sum of the eigenvalues. By the fundamental theorem of algebra (Theorem 1), the characteristic polynomial factors completely as

p_A(\lambda) = \det(\mathbf{A} - \lambda\mathbf{I}) = (\lambda_0 - \lambda)(\lambda_1 - \lambda)\cdots(\lambda_{n-1} - \lambda)

(28)

Compare the coefficients of $\lambda^{n-1}$ on both sides. On the right, after expansion, it is $(-1)^{n-1}\sum_{i=0}^{n-1}\lambda_i$ . On the left, in the expansion of the determinant, the terms containing $\lambda^{n-1}$ can come only from the product of the main diagonal $\prod_{i=0}^{n-1}(a_{ii} - \lambda)$ (every other permutation term contains at most $n-2$ diagonal entries and so has degree at most $n-2$ ), and its coefficient is $(-1)^{n-1}\sum_{i=0}^{n-1} a_{ii} = (-1)^{n-1}\,\text{tr}(\mathbf{A})$ . The two sides are equal, so $\text{tr}(\mathbf{A}) = \sum_{i=0}^{n-1} \lambda_i$ .

2. Determinant and the product of the eigenvalues. Substitute $\lambda = 0$ into the factorization above:

p_A(0) = \det(\mathbf{A}) = \lambda_0 \lambda_1 \cdots \lambda_{n-1}

(29)

3. Powers. Let $\mathbf{A}\mathbf{v} = \lambda\mathbf{v}$ with $\mathbf{v} \neq \mathbf{0}$ . We use induction on $k$ : the case $k=1$ is the hypothesis; if $\mathbf{A}^{k-1}\mathbf{v} = \lambda^{k-1}\mathbf{v}$ , then

\mathbf{A}^k\mathbf{v} = \mathbf{A}(\mathbf{A}^{k-1}\mathbf{v}) = \mathbf{A}(\lambda^{k-1}\mathbf{v}) = \lambda^{k-1}(\mathbf{A}\mathbf{v}) = \lambda^k\mathbf{v}

(30)

4. Inverse. If $\mathbf{A}$ is invertible and $\lambda \neq 0$ , multiplying both sides of $\mathbf{A}\mathbf{v} = \lambda\mathbf{v}$ on the left by $\mathbf{A}^{-1}$ gives $\mathbf{v} = \lambda\,\mathbf{A}^{-1}\mathbf{v}$ , that is, $\mathbf{A}^{-1}\mathbf{v} = \frac{1}{\lambda}\mathbf{v}$ , so $\frac{1}{\lambda}$ is an eigenvalue of $\mathbf{A}^{-1}$ (with the same eigenvector).

5. Shift. $(\mathbf{A} + c\mathbf{I})\mathbf{v} = \mathbf{A}\mathbf{v} + c\mathbf{v} = (\lambda + c)\mathbf{v}$ , and $\mathbf{v} \neq \mathbf{0}$ . $\square$

The relation between the trace and the sum of the eigenvalues can be proved from the expansion of the characteristic polynomial. After the characteristic polynomial $p_A(\lambda) = \det(\mathbf{A} - \lambda\mathbf{I})$ is expanded, the coefficient of $\lambda^{n-1}$ is exactly $(-1)^{n-1}\text{tr}(\mathbf{A})$ .

The relation between the determinant and the product of the eigenvalues explains the invertibility condition for a square matrix: a square matrix $\mathbf{A}$ is invertible if and only if all of its eigenvalues are nonzero, that is, $\det(\mathbf{A}) \neq 0$ .

Exercise 2 (Eigenvalues of Special Matrices)

Prove the following properties:

Triangular matrices. The eigenvalues of an upper or lower triangular matrix are exactly its main diagonal entries.
Idempotent matrices. If $\mathbf{A}^2 = \mathbf{A}$ , then its eigenvalues can only be 0 or 1 (if $\lambda$ is an eigenvalue, then $\lambda^2 = \lambda$ , that is, $\lambda(\lambda-1) = 0$ ).
Orthogonal matrices. If $\mathbf{A}^\top \mathbf{A} = \mathbf{I}$ , or $\mathbf{A}^{*}\mathbf{A} =\mathbf{I}$ , then all of its eigenvalues have absolute value 1.

Example 3

Consider the matrix $\mathbf{A} = \begin{pmatrix} 4 & 2 \\ 1 & 3 \end{pmatrix}$ .

Compute the characteristic polynomial:

\begin{align} p_A(\lambda) &= \det(\mathbf{A} - \lambda\mathbf{I}) \\ &= \det\begin{pmatrix} 4-\lambda & 2 \\ 1 & 3-\lambda \end{pmatrix} \\ &= (4-\lambda)(3-\lambda) - 2 \cdot 1 \\ &= \lambda^2 - 7\lambda + 10 \end{align}

(34)

Solving the characteristic equation $\lambda^2 - 7\lambda + 10 = 0$ gives the eigenvalues $\lambda_0 = 5$ and $\lambda_1 = 2$ .

Verify the matrix properties:

Trace: $\text{tr}(\mathbf{A}) = 4 + 3 = 7 = \lambda_0 + \lambda_1$
Determinant: $\det(\mathbf{A}) = 4 \cdot 3 - 2 \cdot 1 = 10 = \lambda_0 \cdot \lambda_1$

import numpy as np

# --- Visual hierarchy helpers ---

def compare_print(label, actual, expected_desc):
    print(f"\n[{label}]")
    print(f"  Computed value: {actual}")
    print(f"  Expected value: {expected_desc}")

np.set_printoptions(precision=4, suppress=True, linewidth=100)

# ============================================================
print_header("Example 8.2 | Verifying Basic Properties of Eigenvalues")
# ============================================================

# --- Step 1: Define the matrix A ---
print_step("1", "Define the matrix A")
A = np.array([[4, 2],
              [1, 3]])
print(A)

# --- Step 2: Compute the eigenvalues λ₀, λ₁ ---
print_step("2", "Compute the eigenvalues λ₀, λ₁")
eigenvalues = np.linalg.eigvals(A)
print(f"  λ₀ = {eigenvalues[0]:.4f},  λ₁ = {eigenvalues[1]:.4f}")

# --- Step 3: Verify the trace and the determinant ---
print_step("3", "Verify tr(A) = Σλᵢ, det(A) = Πλᵢ")
compare_print(
    "tr(A) vs λ₀ + λ₁",
    f"{np.trace(A):.4f}",
    f"λ₀ + λ₁ = {np.sum(eigenvalues):.4f}"
)
compare_print(
    "det(A) vs λ₀ × λ₁",
    f"{np.linalg.det(A):.4f}",
    f"λ₀ × λ₁ = {np.prod(eigenvalues):.4f}"
)

# --- Step 4: Verify the eigenvalues of A² ---
print_step("4", "Verify that the eigenvalues of A² are λᵢ²")
compare_print(
    "eig(A²)",
    np.linalg.eigvals(A @ A),
    f"λᵢ² = {eigenvalues**2}"
)

# --- Step 5: Verify the eigenvalues of A⁻¹ ---
print_step("5", "Verify that the eigenvalues of A⁻¹ are 1/λᵢ")
compare_print(
    "eig(A⁻¹)",
    np.linalg.eigvals(np.linalg.inv(A)),
    f"1/λᵢ = {1/eigenvalues}"
)

# --- Step 6: Verify the eigenvalues of (A + cI) ---
print_step("6", "Verify that the eigenvalues of (A + cI) are λᵢ + c")
c = 3
compare_print(
    f"eig(A + {c}I)",
    np.linalg.eigvals(A + c * np.eye(2)),
    f"λᵢ + {c} = {eigenvalues + c}"
)

8.2.2 Eigenspaces, Algebraic Multiplicity, and Geometric Multiplicity¶

For each eigenvalue, the corresponding eigenvectors (together with the zero vector) form a vector subspace, called the eigenspace.

Eigenspaces have the following important properties:

Subspace property. An eigenspace is a subspace of the vector space: it is closed under addition and under scalar multiplication.
Invariant subspace. If $\mathbf{v} \in E_\lambda$ , then $\mathbf{A}\mathbf{v} \in E_\lambda$ . That is, the matrix $\mathbf{A}$ maps the eigenspace into itself.
Eigenvectors for distinct eigenvalues are linearly independent. If $\lambda_0, \lambda_1, \ldots, \lambda_{k-1}$ are pairwise distinct, then the corresponding nonzero eigenvectors $\mathbf{v}_0, \ldots, \mathbf{v}_{k-1}$ must be linearly independent (see the theorem below).

This theorem is one of the cornerstones of eigenvalue theory: it guarantees that eigenvectors taken from distinct eigenvalues always form a linearly independent set, which lays the theoretical foundation for the diagonalization of matrices that follows.

Proof 2

We use mathematical induction on the number $m$ of eigenvectors.

Base case ( $m = 1$ ). Since an eigenvector is by definition a nonzero vector, $\{\mathbf{v}_0\}$ is clearly linearly independent.

Inductive step. Suppose $\textcolor{#57068C}{\mathbf{v}_0, \mathbf{v}_1, \ldots, \mathbf{v}_{m-1}}$ are linearly independent, and consider the set obtained by adding $\textcolor{#D62728}{\mathbf{v}_m}$ .

Suppose $\textcolor{#57068C}{\mathbf{v}_0, \ldots, \mathbf{v}_{m-1}},\, \textcolor{#D62728}{\mathbf{v}_m}$ are linearly dependent, that is, there exist coefficients, not all zero, such that

c_0\textcolor{#57068C}{\mathbf{v}_0} + c_1\textcolor{#57068C}{\mathbf{v}_1} + \cdots + c_{m-1}\textcolor{#57068C}{\mathbf{v}_{m-1}} + c_m\textcolor{#D62728}{\mathbf{v}_m} = \mathbf{0}.

(36)

By the induction hypothesis, $\textcolor{#57068C}{\mathbf{v}_0, \ldots, \mathbf{v}_{m-1}}$ are linearly independent, so $c_m \ne 0$ (otherwise (36) would reduce to a vanishing linear combination of the first $m$ vectors with coefficients not all zero, contradicting their independence).
Applying $\mathbf{A}$ to both sides of (36) and using $\mathbf{A}\mathbf{v}_i = \lambda_i \mathbf{v}_i$ , we obtain

c_0\lambda_0\textcolor{#57068C}{\mathbf{v}_0} + \cdots + c_{m-1}\lambda_{m-1}\textcolor{#57068C}{\mathbf{v}_{m-1}} + c_m\lambda_m\textcolor{#D62728}{\mathbf{v}_m} = \mathbf{0}.

(37)

Multiplying (36) by $\lambda_m$ , we obtain

c_0\lambda_m\textcolor{#57068C}{\mathbf{v}_0} + \cdots + c_{m-1}\lambda_m\textcolor{#57068C}{\mathbf{v}_{m-1}} + c_m\lambda_m\textcolor{#D62728}{\mathbf{v}_m} = \mathbf{0}.

(38)

Subtracting (38) from (37) eliminates the terms containing $\textcolor{#D62728}{\mathbf{v}_m}$ and gives
$c_0(\lambda_0 - \lambda_m)\textcolor{#57068C}{\mathbf{v}_0} + \cdots + c_{m-1}(\lambda_{m-1} - \lambda_m)\textcolor{#57068C}{\mathbf{v}_{m-1}} = \mathbf{0}.$
(39)
By the induction hypothesis, every coefficient is zero; and since the eigenvalues are pairwise distinct, $\lambda_i - \lambda_m \ne 0$ , so
$c_0 = c_1 = \cdots = c_{m-1} = 0.$
(40)
Substituting back into (36) gives $c_m\textcolor{#D62728}{\mathbf{v}_m} = \mathbf{0}$ ; but $c_m \ne 0$ and $\textcolor{#D62728}{\mathbf{v}_m} \ne \mathbf{0}$ , a contradiction.

This contradiction shows that the assumption is false, so

\underline{\textcolor{#57068C}{\mathbf{v}_0, \mathbf{v}_1, \ldots, \mathbf{v}_{m-1}},\, \textcolor{#D62728}{\mathbf{v}_m} \text{ are linearly independent.}}

(41)

By induction, the statement holds for all $k \ge 1$ .

A short dimension count, together with the finite-dimensional abstract linear algebra of Chapter 4, readily shows that if the eigenvalues of a square matrix $\mathbf{A}$ are pairwise distinct, then $\mathbf{A}$ can be diagonalized using the corresponding eigenvectors; this means that collecting all the eigenspaces recovers the entire domain. But things do not always go so smoothly—eigenvalues sometimes repeat.

When an eigenvalue $\lambda$ is a repeated root, “repeated” has two quite different meanings, and the gap between these two meanings decides whether the matrix can be diagonalized.

Algebraic multiplicity answers: how “heavy” is this eigenvalue in the polynomial?
Geometric multiplicity answers: how many independent invariant axes can this eigenvalue hold up?

Proof 3

Since $\lambda$ is an eigenvalue, there is at least one nonzero eigenvector, so $g_\lambda \geq 1$ .
Let $\{\mathbf{v}_0, \ldots, \mathbf{v}_{g_\lambda-1}\}$ be a basis of $E_\lambda$ . Fix the field $\mathbb{F}=\mathbb{R}$ or $\mathbb{C}$ with $\lambda\in\mathbb{F}$ , and extend this basis to a full basis of $\mathbb{F}^n$ :
$\mathcal{B} = \{\mathbf{v}_0, \ldots, \mathbf{v}_{g_\lambda-1},\, \mathbf{u}_0, \ldots, \mathbf{u}_{n-g_\lambda-1}\}.$
(44)
In the basis $\mathcal{B}$ , since $\mathbf{A}\mathbf{v}_i = \lambda\mathbf{v}_i$ ( $i < g_\lambda$ ), the matrix $\mathbf{A}$ takes block upper triangular form:
$\mathbf{A}_{\mathcal{B}} = \begin{pmatrix} \lambda\mathbf{I}_{g_\lambda} & \mathbf{B} \\ \mathbf{0} & \mathbf{C} \end{pmatrix}.$
(45)
Since similar matrices have the same characteristic polynomial,
$p_{\mathbf{A}}(t) = p_{\mathbf{A}_{\mathcal{B}}}(t) = (\lambda-t)^{g_\lambda}p_{\mathbf{C}}(t).$
(46)
Hence $\lambda$ occurs at least $g_\lambda$ times in the characteristic polynomial, that is, $\underline{g_\lambda \leq a_\lambda}$ .

Side-by-side comparison.

	Algebraic multiplicity $a_\lambda$	Geometric multiplicity $g_\lambda$
Defined via	the characteristic polynomial (algebra)	the kernel (geometry)
How to compute	the multiplicity of $\lambda$ as a root	$n - \operatorname{rank}(\mathbf{A}-\lambda\mathbf{I})$
Intuitive meaning	how many directions are “expected”	how many independent directions there “actually” are

Over $\mathbb{C}$ (or, more generally, when the characteristic polynomial splits completely over $\mathbb{F}$ ), if $g_\lambda = a_\lambda$ holds for every eigenvalue, the matrix is diagonalizable. If any eigenvalue has $g_\lambda < a_\lambda$ , the Jordan canonical form (§8.4) is needed to handle those “missing directions.”

The next two examples use the same eigenvalue $\lambda = 2$ to illustrate two extreme situations.

Example 4 (Algebraic and Geometric Multiplicities of an Upper Triangular Matrix)

Consider the matrix

\mathbf{A} = \begin{pmatrix} 2 & 1 & 0 \\ 0 & 2 & 0 \\ 0 & 0 & 3 \end{pmatrix}

(47)

Step 1. Read off the eigenvalues. The eigenvalues of an upper triangular matrix are its diagonal entries: $\lambda = 2$ (repeated twice) and $\lambda = 3$ . Hence $a_2 = 2$ and $a_3 = 1$ .

Step 2. Compute the geometric multiplicity of $\lambda = 2$ . Solve $(\mathbf{A} - 2\mathbf{I})\mathbf{v} = \mathbf{0}$ :

\begin{pmatrix} 0 & 1 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} v_0 \\ v_1 \\ v_2 \end{pmatrix} = \mathbf{0} \quad\Rightarrow\quad v_1 = 0,\; v_2 = 0.

(48)

The kernel is $E_2 = \operatorname{span}\!\left\{\begin{pmatrix}1\\0\\0\end{pmatrix}\right\}$ , so $g_2 = 1 < 2 = a_2$ , and the matrix is not diagonalizable.

Step 3. Compute the geometric multiplicity of $\lambda = 3$ . Similarly, $E_3 = \operatorname{span}\!\left\{\begin{pmatrix}0\\0\\1\end{pmatrix}\right\}$ and $g_3 = 1 = a_3$ , so this eigenvalue has no gap.

Example 5 (Algebraic and Geometric Multiplicities of a Nilpotent Matrix)

Consider the nilpotent matrix

\mathbf{A} = \begin{pmatrix} 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \end{pmatrix}

(49)

Step 1. Compute the eigenvalues. $p_A(\lambda) = \lambda^4$ , so the only eigenvalue is $\lambda = 0$ , with $a_0 = 4$ .

Step 2. Compute the geometric multiplicity. The solution space of $\mathbf{A}\mathbf{v} = \mathbf{0}$ is $E_0 = \operatorname{span}\!\left\{\begin{pmatrix}0\\0\\0\\1\end{pmatrix}\right\}$ , and $g_0 = 1 \ll a_0 = 4$ ; the gap is very large, and the matrix is severely non-diagonalizable.

Step 3. Observe the generalized eigenspaces (the pirate-ship analogy). The action of $\mathbf{A}$ on the standard basis is

\mathbf{A}\mathbf{e}_0 = \mathbf{e}_1,\quad \mathbf{A}\mathbf{e}_1 = \mathbf{e}_2,\quad \mathbf{A}\mathbf{e}_2 = \mathbf{e}_3,\quad \mathbf{A}\mathbf{e}_3 = \mathbf{0}.

(50)

Each application of $\mathbf{A}$ moves a vector “down” one position, and a vector that moves off the bottom becomes zero. Picture four people $\mathbf{e}_0, \mathbf{e}_1, \mathbf{e}_2, \mathbf{e}_3$ standing on the deck of a pirate ship, from bow to stern. In each round everyone “steps forward” (one application of $\mathbf{A}$ ), and whoever steps off the stern falls overboard:

Round	Falls overboard this round	$\ker(\mathbf{A}^k)$ (overboard so far)	Dimension
$k=1$	$\mathbf{e}_3$	$\operatorname{span}\{\mathbf{e}_3\}$	1
$k=2$	$\mathbf{e}_2$	$\operatorname{span}\{\mathbf{e}_2,\mathbf{e}_3\}$	2
$k=3$	$\mathbf{e}_1$	$\operatorname{span}\{\mathbf{e}_1,\mathbf{e}_2,\mathbf{e}_3\}$	3
$k=4$	$\mathbf{e}_0$	$\mathbb{R}^4$	4

“Overboard so far” is exactly the chain of generalized eigenspaces:

\ker(\mathbf{A}) \subset \ker(\mathbf{A}^2) \subset \ker(\mathbf{A}^3) \subset \ker(\mathbf{A}^4) = \mathbb{R}^4.

(51)

This round-by-round structure of “falling overboard” is precisely where Jordan blocks come from; §8.4 treats it in detail.

Besides illustrating the situation $g_\lambda < a_\lambda$ , these two examples are highly representative in form—they are exactly the prototypes of the Jordan canonical form, and they preview the conclusions of our study of non-diagonalizable matrices.

Remark 1 (How the Gap Between Multiplicities Shows Up in Difference Equations)

The difference between algebraic and geometric multiplicity has an extremely intuitive manifestation in the structure of the solutions of difference equations. (For further discussion, see Experiment 7.)

Consider the difference equation $\mathbf{x}_{k+1} = \mathbf{A}\mathbf{x}_k$ , where

\mathbf{A} = \begin{pmatrix} 2 & 1 \\ 0 & 2 \end{pmatrix}.

(52)

$\mathbf{A}$ is a $2\times 2$ Jordan block with $p_A(\lambda) = (\lambda-2)^2$ and $a_2 = 2$ , but $g_2 = 1$ —one independent eigendirection is missing.

Writing $\mathbf{A} = 2\mathbf{I} + \mathbf{N}$ (where $\mathbf{N}^2 = \mathbf{0}$ ) and using the binomial expansion:

\mathbf{A}^k = \begin{pmatrix} 2^k & k\cdot 2^{k-1} \\ 0 & 2^k \end{pmatrix}, \qquad \mathbf{x}_k = \begin{pmatrix} 2^k x_0^{(0)} + k\cdot 2^{k-1} x_1^{(0)} \\ 2^k x_1^{(0)} \end{pmatrix}.

(53)

The term $k \cdot 2^{k-1}$ appears in the solution, that is, a term of the form “polynomial $\times$ exponential.” Compare this with the case $g_2 = a_2 = 2$ (that is, $\mathbf{A} = 2\mathbf{I}$ ), where the solution is the pure exponential $2^k\mathbf{x}_0$ with no polynomial factor.

Conclusion.

$g_\lambda = a_\lambda$ (diagonalizable) $\Rightarrow$ the solutions are pure exponentials $c\lambda^k$ .
$g_\lambda < a_\lambda$ (a Jordan block is present) $\Rightarrow$ when $\lambda\neq0$ , the general solution admits “polynomial $\times$ exponential” terms $k^j\lambda^k$ , whose coefficients depend on the initial values; a Jordan block for $\lambda=0$ dies out after finitely many steps.

The size of the largest Jordan block, minus 1, determines the highest polynomial degree that the corresponding solutions can have; the multiplicity gap $a_\lambda-g_\lambda$ by itself is not enough to determine it. This is the most direct manifestation of the Jordan canonical form in applications.

8.2.3 ◆The Cayley–Hamilton Theorem¶

The Cayley–Hamilton theorem links a square matrix to its characteristic equation and gives the concise result $p_{\mathbf{A}}(\mathbf{A}) = \mathbf{0}$ . For this expression to make sense, we first need to define polynomial functions of a matrix.

Remark 2 (A Distinction in Notation)

Be careful to distinguish:

$p(\lambda)$ : a scalar polynomial in the variable $\lambda$
$p(\mathbf{A})$ : the matrix obtained by replacing the powers of $\lambda$ with the powers of $\mathbf{A}$

In particular, $\det(\mathbf{A} - \lambda\mathbf{I})$ expands to a polynomial in $\lambda$ ,

p_{\mathbf{A}}(\lambda) = c_n\lambda^n + \cdots + c_0

(55)

and $p_{\mathbf{A}}(\mathbf{A})$ then means $c_n\mathbf{A}^n + \cdots + c_0\mathbf{I}$ . The polynomial function sends the original constant term $c_0$ to $c_0\mathbf{I}$ , so one cannot simply substitute $\mathbf{A}$ into the original determinant expression to get $\det(\mathbf{A}-\mathbf{AI})= \det(\mathbf{0})=0$ ; these are two different things.

Example 6 (A Pencil-and-Paper Check of the Cayley–Hamilton Theorem)

Take

\mathbf{A} = \begin{pmatrix} 4 & 2 \\ 1 & 3 \end{pmatrix}

(56)

and verify $p_{\mathbf{A}}(\mathbf{A}) = \mathbf{0}$ step by step.

Step 1. Compute the characteristic polynomial.

p_{\mathbf{A}}(\lambda) = \det(\mathbf{A} - \lambda\mathbf{I}) = \det\begin{pmatrix} 4-\lambda & 2 \\ 1 & 3-\lambda \end{pmatrix} = (4-\lambda)(3-\lambda) - 2 \cdot 1

(57)

Expanding:

= \lambda^2 - 7\lambda + 12 - 2 = \lambda^2 - 7\lambda + 10

(58)

Hence $p_{\mathbf{A}}(\lambda) = \lambda^2 - 7\lambda + 10$ .

Step 2. Compute $\mathbf{A}^2$ .

\mathbf{A}^2 = \begin{pmatrix} 4 & 2 \\ 1 & 3 \end{pmatrix} \begin{pmatrix} 4 & 2 \\ 1 & 3 \end{pmatrix} = \begin{pmatrix} 4\cdot4+2\cdot1 & 4\cdot2+2\cdot3 \\ 1\cdot4+3\cdot1 & 1\cdot2+3\cdot3 \end{pmatrix} = \begin{pmatrix} 18 & 14 \\ 7 & 11 \end{pmatrix}

(59)

Step 3. Compute $-7\mathbf{A}$ .

-7\mathbf{A} = \begin{pmatrix} -28 & -14 \\ -7 & -21 \end{pmatrix}

(60)

Step 4. Compute $10\mathbf{I}$ .

10\mathbf{I} = \begin{pmatrix} 10 & 0 \\ 0 & 10 \end{pmatrix}

(61)

Step 5. Combine the pieces and verify $p_{\mathbf{A}}(\mathbf{A}) = \mathbf{A}^2 - 7\mathbf{A} + 10\mathbf{I}$ .

p_{\mathbf{A}}(\mathbf{A}) = \begin{pmatrix} 18 & 14 \\ 7 & 11 \end{pmatrix} + \begin{pmatrix} -28 & -14 \\ -7 & -21 \end{pmatrix} + \begin{pmatrix} 10 & 0 \\ 0 & 10 \end{pmatrix} = \begin{pmatrix} 18-28+10 & 14-14+0 \\ 7-7+0 & 11-21+10 \end{pmatrix} = \begin{pmatrix} 0 & 0 \\ 0 & 0 \end{pmatrix} = \mathbf{0}

(62)

$\mathbf{A}$ indeed satisfies its own characteristic equation.

import numpy as np

np.set_printoptions(precision=4, suppress=True, linewidth=100)

# ============================================================
print_header("Example 8.3 | Verifying the Cayley-Hamilton Theorem")
# ============================================================

# --- Step 1: Define the matrix A ---
print_step("1", "Define the matrix A")
A = np.array([[4, 2],
              [1, 3]])
print(A)

# --- Step 2: Compute the coefficients of the characteristic polynomial ---
print_step("2", "Compute the coefficients of the characteristic polynomial p(λ)")
coeffs = np.poly(A)   # returns [c_n, c_{n-1}, ..., c_0], e.g. [1, -7, 10]
print(f"  coefficients of p(λ) (in descending powers): {coeffs}")
print(f"  that is, p(λ) = λ² - 7λ + 10")

# --- Step 3: Compute the matrix polynomial p(A) ---
print_step("3", "Compute the matrix polynomial p(A) = A² - 7A + 10I")
n = len(A)
p_A = sum(c * np.linalg.matrix_power(A, n - i) for i, c in enumerate(coeffs))

# --- Step 4: Verify the Cayley-Hamilton theorem ---
print_step("4", "Verify p(A) ≈ O (Cayley-Hamilton theorem)")
print("  p(A) =")
print(np.round(p_A, 10))   # ← split over two lines to avoid a Typst PDF tag error
compare_print(
    "Is p(A) the zero matrix?",
    np.allclose(p_A, 0),
    "the Cayley-Hamilton theorem guarantees p(A) = O"
)

Proof 4

Set up the key identity. Let $\mathbf{B}(\lambda) = \mathbf{A} - \lambda\mathbf{I}$ . By the property of the adjugate matrix from Chapter 7 (see Theorem 15), every square matrix $\mathbf{M}$ satisfies $\text{adj}(\mathbf{M})\,\mathbf{M} = \det(\mathbf{M})\,\mathbf{I}$
This identity holds whether $\mathbf{M}$ is invertible or singular. Taking $\mathbf{M} = \mathbf{B}(\lambda)$ , we obtain an identity that holds for all $\lambda \in \mathbb{R}$ :
$\text{adj}(\mathbf{B}(\lambda))\,\mathbf{B}(\lambda) \;=\; p_{\mathbf{A}}(\lambda)\,\mathbf{I}$
(64)
Expand the adjugate as a polynomial in $\lambda$ . The entries of $\mathbf{B}(\lambda)$ are

[\mathbf{B}(\lambda)]_{ij} = \begin{cases} a_{ii} - \lambda, & i = j \\ a_{ij}, & i \neq j \end{cases}

(65)

That is, the diagonal entries are polynomials of degree one in $\lambda$ , and the off-diagonal entries are constants. Each entry of the adjugate matrix is the (signed) determinant of some $(n{-}1)\times(n{-}1)$ submatrix; expanding that determinant over permutations gives a sum of products of $n-1$ entries, and since each entry of $\mathbf{B}(\lambda)$ has degree at most 1 as a polynomial in $\lambda$ , each product has degree at most $n-1$ ; therefore each cofactor, as a polynomial in $\lambda$ , has degree at most $n-1$ .

Hence the adjugate matrix can be written as $\text{adj}(\mathbf{B}(\lambda)) \;=\; \mathbf{C}_{n-1}\lambda^{n-1} + \mathbf{C}_{n-2}\lambda^{n-2} + \cdots + \mathbf{C}_1\lambda + \mathbf{C}_0$

where $\mathbf{C}_0, \mathbf{C}_1, \ldots, \mathbf{C}_{n-1} \in M_{n \times n}(\mathbb{F})$ are constant matrices.

Expand the product and compare coefficients degree by degree. Substituting this into the identity in (2): $\left(\sum_{k=0}^{n-1}\mathbf{C}_k\,\lambda^k\right)\!\left(\mathbf{A} - \lambda\mathbf{I}\right) \;=\; \sum_{j=0}^{n} c_j\,\lambda^j\,\mathbf{I}$

Expand the product on the left—multiply each term in the parentheses by $\mathbf{A}$ and by $(-\lambda\mathbf{I})$ separately:

\underbrace{\sum_{k=0}^{n-1}\mathbf{C}_k\mathbf{A}\,\lambda^k}_{\text{times }\mathbf{A}} \;-\; \underbrace{\sum_{k=0}^{n-1}\mathbf{C}_k\,\lambda^{k+1}}_{\text{times }-\lambda\mathbf{I}}

(66)

In the second sum, let $j = k+1$ (aligning $\lambda^{k+1}$ as $\lambda^j$ ) to obtain

\sum_{k=0}^{n-1}\mathbf{C}_k\mathbf{A}\,\lambda^k \;-\; \sum_{j=1}^{n}\mathbf{C}_{j-1}\,\lambda^j \;=\; \sum_{j=0}^{n} c_j\,\lambda^j\,\mathbf{I}

(67)

Since the identity in (2) holds for all $\lambda \in \mathbb{R}$ , the two sides, as polynomials in $\lambda$ , must have equal coefficients in every degree (each entry of the difference of the two sides is a polynomial in $\lambda$ , and a nonzero polynomial has only finitely many roots). Matching the coefficients of $\lambda^0, \lambda^1, \ldots, \lambda^n$ one by one gives $n+1$ matrix equations:

\begin{cases} \lambda^0\colon & \mathbf{C}_0\mathbf{A} = c_0\,\mathbf{I} & \cdots\; (E_0) \\[4pt] \lambda^j\ (1 \leq j \leq n{-}1)\colon & \mathbf{C}_j\mathbf{A} - \mathbf{C}_{j-1} = c_j\,\mathbf{I} & \cdots\; (E_j) \\[4pt] \lambda^n\colon & -\mathbf{C}_{n-1} = c_n\,\mathbf{I} & \cdots\; (E_n) \end{cases}

(68)

Telescoping sum.

Our goal is to extract $p_{\mathbf{A}}(\mathbf{A}) = \sum_{j=0}^n c_j\,\mathbf{A}^j$ from the system above. Observe that the right side of $(E_j)$ is $c_j\,\mathbf{I}$ ; multiplying it on the right by $\mathbf{A}^j$ gives $c_j\,\mathbf{A}^j$ , which is exactly the $j$ -th term of the target expansion. With this motivation, multiply both sides of $(E_j)$ on the right by $\mathbf{A}^j$ (using the associativity of matrix multiplication to simplify the left side, e.g., $\mathbf{C}_j\mathbf{A}\cdot\mathbf{A}^j = \mathbf{C}_j\mathbf{A}^{j+1}$ ):

\begin{aligned} (E_0)\times\mathbf{I}\;\colon & \quad \mathbf{C}_0\mathbf{A} &&= c_0\,\mathbf{I} \\[2pt] (E_1)\times\mathbf{A}\;\colon & \quad \mathbf{C}_1\mathbf{A}^2 - \mathbf{C}_0\mathbf{A} &&= c_1\,\mathbf{A} \\[2pt] (E_2)\times\mathbf{A}^2\colon & \quad \mathbf{C}_2\mathbf{A}^3 - \mathbf{C}_1\mathbf{A}^2 &&= c_2\,\mathbf{A}^2 \\[2pt] & \quad \vdots && \vdots \\[2pt] (E_{n-1})\times\mathbf{A}^{n-1}\colon & \quad \mathbf{C}_{n-1}\mathbf{A}^n - \mathbf{C}_{n-2}\mathbf{A}^{n-1} &&= c_{n-1}\,\mathbf{A}^{n-1} \\[2pt] (E_n)\times\mathbf{A}^n\colon & \quad -\mathbf{C}_{n-1}\mathbf{A}^n &&= c_n\,\mathbf{A}^n \end{aligned}

(69)

Add these $n+1$ equations together. The left side forms a telescoping sum: each $\mathbf{C}_k\mathbf{A}^{k+1}$ ( $0 \leq k \leq n-1$ ) appears with a plus sign in line $k$ and with a minus sign in line $k+1$ , so the terms cancel in turn, every term disappears, and the left side $= \mathbf{0}$ . The right side collects directly into

c_0\,\mathbf{I} + c_1\,\mathbf{A} + c_2\,\mathbf{A}^2 + \cdots + c_n\,\mathbf{A}^n = p_{\mathbf{A}}(\mathbf{A})

(70)

Therefore $\underline{p_{\mathbf{A}}(\mathbf{A}) = \mathbf{0}}$ .

This proof never substitutes $\mathbf{A}$ for $\lambda$ directly: if we set $\lambda = \mathbf{A}$ directly in the identity in (2), then $\mathbf{B}(\mathbf{A}) = \mathbf{A} - \mathbf{A}\mathbf{I} = \mathbf{0}$ , and the equation collapses to $\mathbf{0}=\mathbf{0}$ , yielding no conclusion at all. The real strategy is to compare the coefficients of each degree first, while $\lambda$ is a scalar, breaking the identity into constant matrix equations, and then to piece these together with a telescoping sum.

Here $\mathbb{F}$ denotes the field containing the entries of the matrix; the Cayley–Hamilton theorem does not require the characteristic polynomial to split over that field.

Example 7 (A Concrete

3\times 3

Illustration of the Proof Strategy for the Cayley–Hamilton Theorem)

This example uses

\mathbf{A} = \begin{pmatrix} 2 & 1 & 0 \\ 0 & 2 & 1 \\ 0 & 0 & 3 \end{pmatrix}

(71)

to illustrate concretely the two core steps of the proof of the theorem: the adjugate identity and the telescoping sum.

Step 1. Write out $\mathbf{B}(\lambda) = \mathbf{A} - \lambda\mathbf{I}$ .

\mathbf{B}(\lambda) = \begin{pmatrix} 2-\lambda & 1 & 0 \\ 0 & 2-\lambda & 1 \\ 0 & 0 & 3-\lambda \end{pmatrix}

(72)

Step 2. Compute the characteristic polynomial $p_\mathbf{A}(\lambda) = \det(\mathbf{B}(\lambda))$ .

Since $\mathbf{B}(\lambda)$ is upper triangular, its determinant equals the product of its diagonal entries:

p_\mathbf{A}(\lambda) = (2-\lambda)^2(3-\lambda) = -\lambda^3 + 7\lambda^2 - 16\lambda + 12

(73)

Hence the coefficients are $c_3 = -1,\; c_2 = 7,\; c_1 = -16,\; c_0 = 12$ .

Step 3. Compute $\text{adj}(\mathbf{B}(\lambda))$ and verify the adjugate identity.

For a $3\times 3$ upper triangular matrix, the adjugate is still upper triangular. Computing the cofactors one by one gives

\text{adj}(\mathbf{B}(\lambda)) = \begin{pmatrix} \lambda^2 - 5\lambda + 6 & \lambda - 3 & 1 \\ 0 & \lambda^2 - 5\lambda + 6 & \lambda - 2 \\ 0 & 0 & \lambda^2 - 4\lambda + 4 \end{pmatrix}

(74)

Check the $(0,0)$ entry of the product (the other entries are similar):

\bigl[\text{adj}(\mathbf{B})\,\mathbf{B}\bigr]_{0,0} = (\lambda^2-5\lambda+6)(2-\lambda) + (\lambda-3)\cdot 0 + 1\cdot 0 = (2-\lambda)^2(3-\lambda) = p_\mathbf{A}(\lambda) \checkmark

(75)

The off-diagonal entries are all 0, which confirms

\text{adj}(\mathbf{B}(\lambda))\,\mathbf{B}(\lambda) = p_\mathbf{A}(\lambda)\,\mathbf{I}

(76)

Step 4. Expand $\text{adj}(\mathbf{B}(\lambda))$ in powers of $\lambda$ and read off the coefficient matrices $\mathbf{C}_2, \mathbf{C}_1, \mathbf{C}_0$ .

\text{adj}(\mathbf{B}(\lambda)) = \mathbf{C}_2\,\lambda^2 + \mathbf{C}_1\,\lambda + \mathbf{C}_0

(77)

where

\mathbf{C}_2 = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix}, \quad \mathbf{C}_1 = \begin{pmatrix} -5 & 1 & 0 \\ 0 & -5 & 1 \\ 0 & 0 & -4 \end{pmatrix}, \quad \mathbf{C}_0 = \begin{pmatrix} 6 & -3 & 1 \\ 0 & 6 & -2 \\ 0 & 0 & 4 \end{pmatrix}

(78)

Substituting into (76) and comparing the coefficients of $\lambda^0, \lambda^1, \lambda^2, \lambda^3$ gives four matrix equations:

\begin{cases} \lambda^0\colon & \mathbf{C}_0\mathbf{A} = c_0\,\mathbf{I} = 12\,\mathbf{I} & \cdots\;(E_0) \\[4pt] \lambda^1\colon & \mathbf{C}_1\mathbf{A} - \mathbf{C}_0 = c_1\,\mathbf{I} = -16\,\mathbf{I} & \cdots\;(E_1) \\[4pt] \lambda^2\colon & \mathbf{C}_2\mathbf{A} - \mathbf{C}_1 = c_2\,\mathbf{I} = 7\,\mathbf{I} & \cdots\;(E_2) \\[4pt] \lambda^3\colon & {-}\mathbf{C}_2 = c_3\,\mathbf{I} = -\mathbf{I} & \cdots\;(E_3) \end{cases}

(79)

Check them one by one:

\mathbf{C}_0\mathbf{A} = \begin{pmatrix} 6 & -3 & 1 \\ 0 & 6 & -2 \\ 0 & 0 & 4 \end{pmatrix} \begin{pmatrix} 2 & 1 & 0 \\ 0 & 2 & 1 \\ 0 & 0 & 3 \end{pmatrix} = \begin{pmatrix} 12 & 0 & 0 \\ 0 & 12 & 0 \\ 0 & 0 & 12 \end{pmatrix} = 12\,\mathbf{I} \checkmark

(80)

\mathbf{C}_1\mathbf{A} - \mathbf{C}_0 = \begin{pmatrix} -10+0 & -5+2 & \cdots \\ \cdots & \cdots & \cdots \\ \cdots & \cdots & \cdots \end{pmatrix} - \mathbf{C}_0 = -16\,\mathbf{I} \checkmark

(81)

\mathbf{C}_2\mathbf{A} - \mathbf{C}_1 = \mathbf{A} - \mathbf{C}_1 = \begin{pmatrix} 7 & 0 & 0 \\ 0 & 7 & 0 \\ 0 & 0 & 7 \end{pmatrix} = 7\,\mathbf{I} \checkmark \qquad {-}\mathbf{C}_2 = -\mathbf{I} \checkmark

(82)

Step 5. Telescoping sum.

Multiply the four equations on the right by $\mathbf{A}^0, \mathbf{A}^1, \mathbf{A}^2, \mathbf{A}^3$ respectively, and add them up:

\begin{aligned} (E_0)\times\mathbf{I}\colon\quad & \mathbf{C}_0\mathbf{A} &&= 12\,\mathbf{I} \\[2pt] (E_1)\times\mathbf{A}\colon\quad & \mathbf{C}_1\mathbf{A}^2 - \mathbf{C}_0\mathbf{A} &&= -16\,\mathbf{A} \\[2pt] (E_2)\times\mathbf{A}^2\colon\quad & \mathbf{C}_2\mathbf{A}^3 - \mathbf{C}_1\mathbf{A}^2 &&= 7\,\mathbf{A}^2 \\[2pt] (E_3)\times\mathbf{A}^3\colon\quad & {-}\mathbf{C}_2\mathbf{A}^3 &&= -\mathbf{A}^3 \end{aligned}

(83)

Adding the left sides, each term $\mathbf{C}_k\mathbf{A}^{k+1}$ cancels against its negative:

\underbrace{\mathbf{C}_0\mathbf{A}}_{} \underbrace{- \mathbf{C}_0\mathbf{A} + \mathbf{C}_1\mathbf{A}^2}_{} \underbrace{- \mathbf{C}_1\mathbf{A}^2 + \mathbf{C}_2\mathbf{A}^3}_{} \underbrace{- \mathbf{C}_2\mathbf{A}^3}_{} = \mathbf{0}

(84)

The right sides collect into $p_\mathbf{A}(\mathbf{A})$ , so

-\mathbf{A}^3 + 7\mathbf{A}^2 - 16\mathbf{A} + 12\,\mathbf{I} = \mathbf{0}

(85)

in complete agreement with the characteristic polynomial of Step 2.

import numpy as np
from sympy import symbols, Matrix, det, factor, expand, latex, eye, pprint

np.set_printoptions(precision=4, suppress=True, linewidth=100)

def print_header(title):
    print("=" * 60)
    print(f"  {title}")
    print("=" * 60)

def print_step(step, desc):
    print(f"\n▶ Step {step}: {desc}")
    print("-" * 40)

def compare_print(label, actual, expected_desc):
    print(f"\n[{label}]")
    print(f"  Computed value: {actual}")
    print(f"  Expected value: {expected_desc}")

# ============================================================
print_header("Example: The Proof Strategy of the Cayley-Hamilton Theorem (3×3)")
# ============================================================

lam = symbols('lambda')

# --- Step 1: Define the matrix A and B(λ) ---
print_step("1", "Define the matrix A and B(λ) = A - λI")
A = Matrix([[2, 1, 0],
            [0, 2, 1],
            [0, 0, 3]])
B = A - lam * eye(3)
print("A ="); pprint(A)
print("\nB(λ) = A - λI ="); pprint(B)

# --- Step 2: Compute the characteristic polynomial p_A(λ) ---
print_step("2", "Compute the characteristic polynomial p_A(λ) = det(B(λ))")
p = expand(det(B))
p_factored = factor(p)
print(f"  p_A(λ) = {p}")
compare_print(
    "Characteristic polynomial (expanded vs factored)",
    p,
    f"factored = {p_factored}"
)

# --- Step 3: Compute adj(B(λ)) ---
print_step("3", "Compute adj(B(λ)) and verify adj(B)·B = p_A(λ)·I")
adj_B = B.adjugate()
print("adj(B(λ)) ="); pprint(adj_B)

product = expand(adj_B * B)
expected = expand(p * eye(3))
compare_print(
    "adj(B(λ))·B(λ) vs p_A(λ)·I",
    "(see the matrices below)",
    "the two should be equal"
)
print("  adj(B)·B ="); pprint(product)
print("  p_A(λ)·I ="); pprint(expected)
print(f"\n  Equality check: {product == expected}")

# --- Step 4: Expand adj(B(λ)) into the coefficient matrices C2, C1, C0 ---
print_step("4", "Expand adj(B(λ)) = C₂λ² + C₁λ + C₀ and read off the coefficient matrices")
adj_expanded = adj_B.applyfunc(expand)

def coeff_matrix(M, var, deg):
    """Extract the coefficient of var^deg from every entry of the matrix M"""
    return M.applyfunc(lambda expr: expr.coeff(var, deg))

C2 = coeff_matrix(adj_expanded, lam, 2)
C1 = coeff_matrix(adj_expanded, lam, 1)
C0 = coeff_matrix(adj_expanded, lam, 0)

print("C₂ ="); pprint(C2)
print("C₁ ="); pprint(C1)
print("C₀ ="); pprint(C0)

# Check the reconstruction: C2*λ² + C1*λ + C0 should equal adj_expanded
reconstructed = expand(C2 * lam**2 + C1 * lam + C0)
compare_print(
    "Reconstruction check of C₂λ² + C₁λ + C₀",
    reconstructed == adj_expanded.applyfunc(expand),
    "True (identical to the original adj(B(λ)))"
)

# --- Step 5: Extract the coefficients c3, c2, c1, c0 and check the four matrix equations ---
print_step("5", "Compare the coefficients of each power of λ and check the system (E0)–(E3)")
from sympy import Poly
coeffs = Poly(p, lam).all_coeffs()  # [c3, c2, c1, c0]
c3, c2, c1, c0 = coeffs
print(f"  coefficients of the characteristic polynomial: c₃={c3}, c₂={c2}, c₁={c1}, c₀={c0}")

I3 = eye(3)
E0 = C0 * A - c0 * I3
E1 = C1 * A - C0 - c1 * I3
E2 = C2 * A - C1 - c2 * I3
E3 = -C2 - c3 * I3

compare_print("(E0): C₀A = c₀I", E0 == Matrix.zeros(3), "should be the zero matrix")
compare_print("(E1): C₁A - C₀ = c₁I", E1 == Matrix.zeros(3), "should be the zero matrix")
compare_print("(E2): C₂A - C₁ = c₂I", E2 == Matrix.zeros(3), "should be the zero matrix")
compare_print("(E3): -C₂ = c₃I", E3 == Matrix.zeros(3), "should be the zero matrix")

# --- Step 6: Telescoping sum; check p_A(A) = 0 ---
print_step("6", "Telescoping sum: check p_A(A) = c₃A³ + c₂A² + c₁A + c₀I = 0")
A2 = A ** 2
A3 = A ** 3
p_A = c3 * A3 + c2 * A2 + c1 * A + c0 * I3
print("p_A(A) ="); pprint(p_A)
compare_print(
    "Cayley-Hamilton theorem: p_A(A) = 0",
    p_A == Matrix.zeros(3),
    "True (every square matrix satisfies its own characteristic equation)"
)

Example 8 (Characteristic and Minimal Polynomials of Jordan Blocks)

Take the matrix (made up of one $2\times 2$ Jordan block and one $1\times 1$ Jordan block)

\mathbf{J} = \begin{pmatrix} 2 & 1 & 0 \\ 0 & 2 & 0 \\ 0 & 0 & 2 \end{pmatrix}

(86)

Step 1. Compute the characteristic polynomial.

$\mathbf{J}$ is upper triangular, so its determinant equals the product of its diagonal entries:

p_\mathbf{J}(\lambda) = (2-\lambda)^3 = -(\lambda-2)^3

(87)

By the Cayley–Hamilton theorem, $p_\mathbf{J}(\mathbf{J}) = -(\mathbf{J}-2\mathbf{I})^3 = \mathbf{0}$ , so the degree of the minimal polynomial is at most 3.

Step 2. Introduce the nilpotent matrix $\mathbf{N} = \mathbf{J} - 2\mathbf{I}$ and compute its powers one at a time.

\mathbf{N} = \begin{pmatrix} 0 & 1 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}

(88)

\mathbf{N}^2 = \begin{pmatrix} 0 & 1 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} 0 & 1 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix} = \begin{pmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix} = \mathbf{0}

(89)

Step 3. Determine the minimal polynomial.

From the above, $\mathbf{N}^1 \neq \mathbf{0}$ and $\mathbf{N}^2 = \mathbf{0}$ , so the monic polynomial of lowest degree with $q(\mathbf{J}) = \mathbf{0}$ is

m_\mathbf{J}(\lambda) = (\lambda - 2)^2

(90)

Step 4. Compare the two polynomials.

p_\mathbf{J}(\lambda) = -(\lambda-2)^3 \quad (\text{degree 3}), \qquad m_\mathbf{J}(\lambda) = (\lambda-2)^2 \quad (\text{degree 2})

(91)

The minimal polynomial properly divides the characteristic polynomial: $m_\mathbf{J} \mid p_\mathbf{J}$ , and

\frac{p_\mathbf{J}(\lambda)}{m_\mathbf{J}(\lambda)} = -(\lambda - 2)

(92)

The two degrees are not equal. The sets of roots are still the same (every root is $\lambda = 2$ ), in accordance with property 3 of Exercise 3.

Step 5. Understand where the difference in degree comes from.

Split $\mathbf{J}$ into blocks according to its Jordan block structure:

\mathbf{J} = \underbrace{\begin{pmatrix} 2 & 1 \\ 0 & 2 \end{pmatrix}}_{\mathbf{J}_2,\; 2\times 2 \text{ block}} \oplus \underbrace{\begin{pmatrix} 2 \end{pmatrix}}_{\mathbf{J}_1,\; 1\times 1 \text{ block}}

(93)

The minimal polynomial of $\mathbf{J}_2$ is $(\lambda-2)^2$ (it takes 2 steps to reach zero).
The minimal polynomial of $\mathbf{J}_1$ is $(\lambda-2)^1$ (1 step suffices).

The minimal polynomial of the whole matrix $\mathbf{J}$ is the least common multiple of the minimal polynomials of the blocks:

m_\mathbf{J}(\lambda) = \text{lcm}\bigl[(\lambda-2)^2,\,(\lambda-2)^1\bigr] = (\lambda-2)^2

(94)

whereas the characteristic polynomial is the product of the characteristic polynomials of the blocks:

p_\mathbf{J}(\lambda) = -(\lambda-2)^2 \cdot (\lambda-2)^1 = -(\lambda-2)^3

(95)

This is the fundamental reason why the degree of the minimal polynomial (2) is strictly smaller than that of the characteristic polynomial (3): the nilpotency depth of the $1\times 1$ block is “covered” by that of the $2\times 2$ block.

Exercise 3 (◆Properties of the Minimal Polynomial)

Use division with remainder and the factor theorem to prove the following:

$m_{\mathbf{A}}(\lambda) \mid p_{\mathbf{A}}(\lambda)$ (the minimal polynomial divides the characteristic polynomial)
For every polynomial $q(\lambda)$ : $q(\mathbf{A}) = \mathbf{0} \Leftrightarrow m_{\mathbf{A}}(\lambda) \mid q(\lambda)$
The roots of $m_{\mathbf{A}}(\lambda)$ are exactly the distinct eigenvalues of $\mathbf{A}$

§8.2 Summary

This section asked: besides being “the scaling ratio along an invariant direction,” what global information about a matrix do eigenvalues give us?

The answer is richer than one might expect. At the quantitative level, the trace equals the sum of all the eigenvalues, and the determinant equals their product—these two equations anchor global properties of the matrix directly to its eigenvalue structure, and they also give invertibility a concise equivalent condition: all eigenvalues are nonzero.

At the structural level, the inequality $1 \leq g_\lambda \leq a_\lambda$ between geometric and algebraic multiplicity is the central result of this section. The geometric multiplicity is the actual dimension of the eigenspace, and the algebraic multiplicity is the multiplicity of the eigenvalue as a root of the characteristic polynomial; the gap between them is precisely the measure of “how far a matrix is from being diagonalizable”—if the characteristic polynomial splits completely over the field in use (which always holds over $\mathbb{C}$ ) and $g_\lambda = a_\lambda$ for every eigenvalue, the matrix is diagonalizable; otherwise the Jordan theory of §8.4 is needed to deal with what remains.

Another high point of this section is the Cayley–Hamilton theorem: every matrix is a solution of its own characteristic equation, that is, $p_{\mathbf{A}}(\mathbf{A}) = \mathbf{0}$ . At first glance the result looks like a magic trick, but in fact it reflects a deep connection between the algebra of matrix polynomials and eigenstructure. The introduction of the minimal polynomial shows further that the set of roots of the lowest-degree polynomial annihilating a matrix is exactly the set of all distinct eigenvalues of the matrix. This thread is developed in full in the Jordan decomposition of §8.4.

8.3 Similarity Transformations and Diagonalization¶

In Sections 8.1 and 8.2 we studied eigenvalues and their eigenspaces one at a time. Now we ask a global question: taken together, can these eigenspaces rebuild the entire vector space?

If the answer is yes, we can find a basis made up entirely of eigenvectors, and in that basis the matrix representation of the linear transformation becomes extremely simple—this is diagonalization. To reach this goal we need a tool for “switching viewpoints”: the similarity transformation. It lets us move freely between different bases and thus look for the basis that best reveals the structure of the matrix.

8.3.1 Similarity Transformations: Looking at the Same Thing from Another Angle¶

Geometric meaning. $\mathbf{A}$ and $\mathbf{B}$ are the matrix representations of the same linear transformation $T: V \to V$ in different bases.

Suppose we have two bases:

The standard basis $\mathcal{B}_{\text{std}} = \{\mathbf{e}_0, \mathbf{e}_1, \ldots, \mathbf{e}_{n-1}\}$ : the matrix representation of the linear transformation $T$ is $\mathbf{A}$ .
A new basis $\mathcal{B}_{\text{new}} = \{\mathbf{v}_0, \mathbf{v}_1, \ldots, \mathbf{v}_{n-1}\}$ : the matrix representation of the same $T$ is $\mathbf{B}$ .

Here $\mathbf{P} = [\mathbf{v}_0, \mathbf{v}_1, \ldots, \mathbf{v}_{n-1}]$ is the change-of-basis matrix.

Why is the formula $\mathbf{B} = \mathbf{P}^{-1}\mathbf{A}\mathbf{P}$ ?

\begin{array}{ccccc} [\mathbf{x}]_{\text{new basis}} & \overset{\mathbf{P}}{\rightarrow} & [\mathbf{x}]_{\text{standard basis}} \\[10pt] {\mathbf{B}}↓ & & ↓{\mathbf{A}} \\[10pt] [\mathbf{y}]_{\text{new basis}} & \overset{\mathbf{P}^{-1}}{\leftarrow} & [\mathbf{y}]_{\text{standard basis}} \end{array}

(111)

Consider the coordinates $[\mathbf{x}]_{\mathcal{B}_{\text{new}}}$ of a vector $\mathbf{x}$ in the new basis. To compute the coordinates of $T(\mathbf{x})$ in the new basis, we go through the following steps:

[\mathbf{x}]_{\text{new}} \overset{\mathbf{P}}{\rightarrow} [\mathbf{x}]_{\text{std}} \overset{\mathbf{A}}{\rightarrow} [\mathbf{y}]_{\text{std}} \overset{\mathbf{P}^{-1}}{\rightarrow}[\mathbf{y}]_{\text{new}}

(112)

Therefore, the matrix representation of the linear transformation in the new basis $\mathcal{B}_{\text{new}}$ is

\underline{\mathbf{B} = \mathbf{P}^{-1}\mathbf{A}\mathbf{P}}

(113)

Example 9 (Similarity Transformation and Diagonalization)

Consider the matrix $\mathbf{A} = \begin{pmatrix} 4 & 2 \\ 1 & 3 \end{pmatrix}$ (the example from Section 8.2).

Recap. In Section 8.2.1 we already computed:

Eigenvalues: 5 and 2 (denoted $\lambda_0 = 5, \lambda_1 = 2$ in Section 8.2.1; in this example we relabel them as $\lambda_0 = 2, \lambda_1 = 5$ )
Trace: $\text{tr}(\mathbf{A}) = 7$ ; determinant: $\det(\mathbf{A}) = 10$

Diagonalization steps.

Compute the eigenvectors.

For $\lambda_0 = 2$ , solve $(\mathbf{A} - 2\mathbf{I})\mathbf{v} = \mathbf{0}$ :

\begin{pmatrix} 2 & 2 \\ 1 & 1 \end{pmatrix}\begin{pmatrix} v_0 \\ v_1 \end{pmatrix} = \mathbf{0} \quad \Rightarrow \quad v_0 + v_1 = 0

(114)

Take the eigenvector $\mathbf{v}_0 = \begin{pmatrix} 1 \\ -1 \end{pmatrix}$ .

For $\lambda_1 = 5$ , solve $(\mathbf{A} - 5\mathbf{I})\mathbf{v} = \mathbf{0}$ :

\begin{pmatrix} -1 & 2 \\ 1 & -2 \end{pmatrix}\begin{pmatrix} v_0 \\ v_1 \end{pmatrix} = \mathbf{0} \quad \Rightarrow \quad v_0 = 2v_1

(115)

Take the eigenvector $\mathbf{v}_1 = \begin{pmatrix} 2 \\ 1 \end{pmatrix}$ .

Construct the diagonalizing matrices.

\mathbf{P} = [\mathbf{v}_0, \mathbf{v}_1] = \begin{pmatrix} 1 & 2 \\ -1 & 1 \end{pmatrix}, \quad \mathbf{D} = \begin{pmatrix} 2 & 0 \\ 0 & 5 \end{pmatrix}

(116)

\mathbf{P}^{-1} = \frac{1}{3}\begin{pmatrix} 1 & -2 \\ 1 & 1 \end{pmatrix}

(117)

Verify $\mathbf{P}^{-1}\mathbf{A}\mathbf{P} = \mathbf{D}$

\mathbf{P}^{-1}\mathbf{A}\mathbf{P} = \frac{1}{3}\begin{pmatrix} 1 & -2 \\ 1 & 1 \end{pmatrix} \begin{pmatrix} 4 & 2 \\ 1 & 3 \end{pmatrix} \begin{pmatrix} 1 & 2 \\ -1 & 1 \end{pmatrix} = \begin{pmatrix} 2 & 0 \\ 0 & 5 \end{pmatrix}

(118)

Geometric interpretation.

In the standard basis, the effect of $\mathbf{A}$ is complicated (it both stretches and rotates).
In the eigenvector basis $\{\mathbf{v}_0, \mathbf{v}_1\}$ , the transformation becomes pure stretching:
- stretching by a factor of 2 along the direction of $\mathbf{v}_0 = \begin{pmatrix} 1 \\ -1 \end{pmatrix}$
- stretching by a factor of 5 along the direction of $\mathbf{v}_1 = \begin{pmatrix} 2 \\ 1 \end{pmatrix}$

Verify the invariants.

Eigenvalues: the diagonal of $\mathbf{D}$ is $2, 5$ ✓
Trace: $\text{tr}(\mathbf{D}) = 2 + 5 = 7 = \text{tr}(\mathbf{A})$ ✓
Determinant: $\det(\mathbf{D}) = 2 \times 5 = 10 = \det(\mathbf{A})$ ✓

import numpy as np

# --- Visual hierarchy helpers ---
def print_header(title):
    print("=" * 60)
    print(f"  {title}")
    print("=" * 60)

def print_step(step, desc):
    print(f"\n▶ Step {step}: {desc}")
    print("-" * 40)

def compare_print(label, actual, expected_desc):
    print(f"\n[{label}]")
    print(f"  Computed value: {actual}")
    print(f"  Expected value: {expected_desc}")

np.set_printoptions(precision=4, suppress=True, linewidth=100)

# ============================================================
print_header("Example 8.4 | Diagonalization and Verification of the Invariants")
# ============================================================

# --- Step 1: Define the matrix A ---
print_step("1", "Define the matrix A")
A = np.array([[4, 2],
              [1, 3]])
print(A)

# --- Step 2: Compute the eigenvalues λᵢ and the eigenvector matrix P ---
print_step("2", "Compute the eigenvalues λᵢ and the eigenvector matrix P")
eigenvalues, eigenvectors = np.linalg.eig(A)
P = eigenvectors
print(f"  Eigenvalues: λ₀ = {eigenvalues[0]:.4f},  λ₁ = {eigenvalues[1]:.4f}")
print("\n  Eigenvector matrix P =")
print(P)

# --- Step 3: Construct the diagonal matrix D ---
print_step("3", "Construct the diagonal matrix D = diag(λ₀, λ₁)")
D = np.diag(eigenvalues)
print(D)

# --- Step 4: Verify the diagonalization P⁻¹AP = D ---
print_step("4", "Verify the diagonalization: P⁻¹AP = D")
P_inv_A_P = np.linalg.inv(P) @ A @ P
print("  P⁻¹AP =")
print(P_inv_A_P)
compare_print(
    "‖P⁻¹AP − D‖ (Frobenius error)",
    f"{np.linalg.norm(P_inv_A_P - D):.2e}",
    "a similarity transformation does not change the eigenvalues; the error should be close to 0"
)

# --- Step 5: Verify the similarity invariants ---
print_step("5", "Verify the similarity invariants: tr and det")
compare_print(
    "tr(A) vs tr(D)",
    f"{np.trace(A):.4f}",
    f"tr(D) = {np.trace(D):.4f}"
)
compare_print(
    "det(A) vs det(D)",
    f"{np.linalg.det(A):.4f}",
    f"det(D) = {np.linalg.det(D):.4f}"
)

Proof 5

Let $\mathbf{B} = \mathbf{P}^{-1}\mathbf{A}\mathbf{P}$ .

1. Characteristic polynomial.

\begin{align} p_B(\lambda) &= \det(\mathbf{B} - \lambda\mathbf{I}) \\ &= \det(\mathbf{P}^{-1}\mathbf{A}\mathbf{P} - \lambda\mathbf{I}) \\ &= \det(\mathbf{P}^{-1}\mathbf{A}\mathbf{P} - \lambda\mathbf{P}^{-1}\mathbf{I}\mathbf{P}) \\ &= \det(\mathbf{P}^{-1}(\mathbf{A} - \lambda\mathbf{I})\mathbf{P}) \\ &= \det(\mathbf{P}^{-1})\det(\mathbf{A} - \lambda\mathbf{I})\det(\mathbf{P}) \\ &= \det(\mathbf{A} - \lambda\mathbf{I}) = p_A(\lambda) \end{align}

(119)

2. Eigenvalues. Since the characteristic polynomials are the same, their roots (the eigenvalues) are naturally the same.

3. Determinant. Set $\lambda = 0$ in the characteristic polynomial, or compute directly:

\det(\mathbf{B}) = \det(\mathbf{P}^{-1}\mathbf{A}\mathbf{P}) = \det(\mathbf{P}^{-1})\det(\mathbf{A})\det(\mathbf{P}) = \det(\mathbf{A})

(120)

4. Trace. Using the cyclic property of the trace, $\text{tr}(\mathbf{XYZ}) = \text{tr}(\mathbf{YZX})$ :

\text{tr}(\mathbf{B}) = \text{tr}(\mathbf{P}^{-1}\mathbf{A}\mathbf{P}) = \text{tr}(\mathbf{A}\mathbf{P}\mathbf{P}^{-1}) = \text{tr}(\mathbf{A})

(121)

5–6. Rank and nullity. A similarity transformation does not change the dimensions of the image and the kernel of the linear map, because the change of basis is invertible.

Intuition. These invariants are all intrinsic properties of the linear transformation itself and do not change with the viewpoint (the choice of basis):

Eigenvalues: the stretching factors along the eigendirections, independent of the coordinate system
Determinant: in real geometry, the signed factor by which volume is scaled (its absolute value is the ordinary volume factor), independent of the coordinate system
Trace: the sum of all the eigenvalues, reflecting the overall “strength” of the transformation
Rank: the dimension of the image, that is, the number of “effective directions” after the transformation
Nullity: the dimension of the kernel, that is, the number of directions “flattened to zero”

8.3.2 Definition of and Conditions for Diagonalization¶

With the concept of a similarity transformation in hand, we can now pose an important question: can we find a special basis in which the matrix takes diagonal form?

Diagonalization is exactly such a special similarity transformation: the new basis is made up of eigenvectors.

Geometric meaning of diagonalization. Diagonalization means that, in the basis formed by the eigenvectors, the linear transformation reduces to simple scaling along each coordinate axis, with no rotation or shear at all.

Proof 6

Forward direction (sufficiency). Let $\mathbf{P} = [\mathbf{v}_0, \mathbf{v}_1, \ldots, \mathbf{v}_{n-1}]$ consist of $n$ linearly independent eigenvectors, with corresponding eigenvalues $\lambda_0, \lambda_1, \ldots, \lambda_{n-1}$ .

Then

\mathbf{A}\mathbf{P} = \mathbf{A}[\mathbf{v}_0, \mathbf{v}_1, \ldots, \mathbf{v}_{n-1}] = [\lambda_0\mathbf{v}_0, \lambda_1\mathbf{v}_1, \ldots, \lambda_{n-1}\mathbf{v}_{n-1}]

(123)

On the other hand,

\mathbf{P}\mathbf{D} = [\mathbf{v}_0, \mathbf{v}_1, \ldots, \mathbf{v}_{n-1}] \begin{pmatrix} \lambda_0 & & \\ & \ddots & \\ & & \lambda_{n-1} \end{pmatrix} = [\lambda_0\mathbf{v}_0, \lambda_1\mathbf{v}_1, \ldots, \lambda_{n-1}\mathbf{v}_{n-1}]

(124)

Hence $\mathbf{A}\mathbf{P} = \mathbf{P}\mathbf{D}$ .

Since the eigenvectors are linearly independent, $\mathbf{P}$ is invertible; multiplying both sides on the left by $\mathbf{P}^{-1}$ gives

\underline{\mathbf{P}^{-1}\mathbf{A}\mathbf{P} = \mathbf{D}}

(125)

Converse (necessity). Suppose $\mathbf{P}^{-1}\mathbf{A}\mathbf{P} = \mathbf{D}$ , where $\mathbf{D} = \text{diag}(\lambda_0, \ldots, \lambda_{n-1})$ .

Multiplying both sides on the left by $\mathbf{P}$ : $\mathbf{A}\mathbf{P} = \mathbf{P}\mathbf{D}$

Write $\mathbf{P} = [\mathbf{v}_0, \ldots, \mathbf{v}_{n-1}]$ and expand the equation above column by column:

\mathbf{A}\mathbf{v}_i = \lambda_i\mathbf{v}_i, \quad i = 0, 1, \ldots, n-1

(126)

This shows that $\mathbf{v}_i$ is an eigenvector corresponding to the eigenvalue $\lambda_i$ . Since $\mathbf{P}$ is invertible, its columns are linearly independent, so $\mathbf{A}$ has $n$ linearly independent eigenvectors.

This corollary follows directly from Theorem Theorem 2: eigenvectors corresponding to pairwise distinct eigenvalues are necessarily linearly independent.

Important

Important. Having $n$ pairwise distinct eigenvalues is a sufficient condition for diagonalizability, but not a necessary one. A matrix with repeated eigenvalues may still be diagonalizable (for example, the identity matrix $\mathbf{I}$ has only one eigenvalue $\lambda = 1$ , yet it is diagonalizable).

A matrix is diagonalizable over a field $\mathbb{F}$ if and only if its characteristic polynomial splits completely over $\mathbb{F}$ and for every eigenvalue $\lambda$ , the geometric multiplicity equals the algebraic multiplicity, that is,

g_\lambda = a_\lambda, \quad \forall \lambda

(127)

Equivalently, the dimensions of all the eigenspaces add up to $n$ :

\sum_{\lambda} \dim(E_\lambda) = n

(128)

This guarantees that $n$ linearly independent eigenvectors can be taken from the eigenspaces to form a complete basis.

Example 10

Diagonalize the matrix $\mathbf{A} = \begin{pmatrix} 3 & 1 \\ 1 & 3 \end{pmatrix}$ .

From the computation in Section 8.1 we know:

Eigenvalues: $\lambda_0 = 2, \lambda_1 = 4$
Eigenvectors: $\mathbf{v}_0 = \begin{pmatrix} -1 \\ 1 \end{pmatrix}, \mathbf{v}_1 = \begin{pmatrix} 1 \\ 1 \end{pmatrix}$

Construct the matrices:

\mathbf{P} = [\mathbf{v}_0, \mathbf{v}_1] = \begin{pmatrix} -1 & 1 \\ 1 & 1 \end{pmatrix}, \quad \mathbf{D} = \begin{pmatrix} 2 & 0 \\ 0 & 4 \end{pmatrix}

(129)

Compute $\mathbf{P}^{-1}$ :

\mathbf{P}^{-1} = \frac{1}{-1-1}\begin{pmatrix} 1 & -1 \\ -1 & -1 \end{pmatrix} = \begin{pmatrix} -0.5 & 0.5 \\ 0.5 & 0.5 \end{pmatrix}

(130)

Verify the diagonalization:

\mathbf{P}^{-1}\mathbf{A}\mathbf{P} = \begin{pmatrix} -0.5 & 0.5 \\ 0.5 & 0.5 \end{pmatrix} \begin{pmatrix} 3 & 1 \\ 1 & 3 \end{pmatrix} \begin{pmatrix} -1 & 1 \\ 1 & 1 \end{pmatrix} = \begin{pmatrix} 2 & 0 \\ 0 & 4 \end{pmatrix} = \mathbf{D}

(131)

8.3.3 Applications of Diagonalization¶

Diagonalization turns complicated matrix operations into simple operations on diagonal matrices. Diagonal matrices have the following features:

Powers are easy. For $\mathbf{D}^k$ , simply raise each diagonal entry to the $k$ -th power.
Functions are easy. For $f(\mathbf{D})$ , simply apply the function to each diagonal entry.
Decoupling. The diagonal form decomposes the vector space into independent one-dimensional subspaces.

These properties make diagonalization extremely useful in the following settings:

Application 1: Computing powers of a matrix

If $\mathbf{A} = \mathbf{P}\mathbf{D}\mathbf{P}^{-1}$ , then

\mathbf{A}^k = (\mathbf{P}\mathbf{D}\mathbf{P}^{-1})^k = \mathbf{P}\mathbf{D}^k\mathbf{P}^{-1}

(132)

and powers of a diagonal matrix are very easy to compute:

\mathbf{D}^k = \begin{pmatrix} \lambda_0^k & 0 & \cdots & 0 \\ 0 & \lambda_1^k & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \lambda_{n-1}^k \end{pmatrix}

(133)

Application 2: Solving difference equations

Consider the system of linear difference equations

\mathbf{x}_{k+1} = \mathbf{A}\mathbf{x}_k

(136)

If $\mathbf{A}$ is diagonalizable, then

\mathbf{x}_k = \mathbf{A}^k\mathbf{x}_0 = \mathbf{P}\mathbf{D}^k\mathbf{P}^{-1}\mathbf{x}_0

(137)

This reduces the problem to simple operations on a diagonal matrix.

Example 12 (The Fibonacci Sequence)

The Fibonacci sequence is defined by $F_0 = 0, F_1 = 1$ and $F_{n+1} = F_n + F_{n-1}$ for all $n \geq 1$ . We use diagonalization to derive a closed-form formula for the Fibonacci sequence.

Step 1. Set up the matrix form

Write the recurrence in vector form:

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

(138)

Let $\mathbf{A} = \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix}$ and $\mathbf{x}_n = \begin{pmatrix} F_{n+1} \\ F_n \end{pmatrix}$ ; then $\mathbf{x}_n = \mathbf{A}\mathbf{x}_{n-1} = \mathbf{A}^n\mathbf{x}_0$ , where $\mathbf{x}_0 = \begin{pmatrix} 1 \\ 0 \end{pmatrix}$ .

Step 2. Compute the eigenvalues

\det(\mathbf{A} - \lambda\mathbf{I}) = \det\begin{pmatrix} 1-\lambda & 1 \\ 1 & -\lambda \end{pmatrix} = (1-\lambda)(-\lambda) - 1 = \lambda^2 - \lambda - 1 = 0

(139)

By the quadratic formula,

\lambda = \frac{1 \pm \sqrt{5}}{2}

(140)

Write $\phi = \frac{1+\sqrt{5}}{2} \approx 1.618$ (the golden ratio) and $\hat{\phi} = \frac{1-\sqrt{5}}{2} \approx -0.618$ .

Note that $\phi + \hat{\phi} = 1$ , $\phi \cdot \hat{\phi} = -1$ , and $\phi^2 = \phi + 1$ .

Step 3. Compute the eigenvectors

For $\lambda_0 = \phi$ , solve $(\mathbf{A} - \phi\mathbf{I})\mathbf{v} = \mathbf{0}$ :

\begin{pmatrix} 1-\phi & 1 \\ 1 & -\phi \end{pmatrix}\begin{pmatrix} v_0 \\ v_1 \end{pmatrix} = \mathbf{0}

(141)

Row 1 gives $v_0 = \phi v_1$ ; taking $v_1 = 1$ gives $\mathbf{v}_0 = \begin{pmatrix} \phi \\ 1 \end{pmatrix}$ .

Similarly, for $\lambda_1 = \hat{\phi}$ we obtain $\mathbf{v}_1 = \begin{pmatrix} \hat{\phi} \\ 1 \end{pmatrix}$ .

Step 4. Construct the diagonalization

\mathbf{P} = \begin{pmatrix} \phi & \hat{\phi} \\ 1 & 1 \end{pmatrix}, \quad \mathbf{D} = \begin{pmatrix} \phi & 0 \\ 0 & \hat{\phi} \end{pmatrix}

(142)

Compute $\mathbf{P}^{-1}$ :

\det(\mathbf{P}) = \phi - \hat{\phi} = \frac{1+\sqrt{5}}{2} - \frac{1-\sqrt{5}}{2} = \sqrt{5}

(143)

\mathbf{P}^{-1} = \frac{1}{\sqrt{5}}\begin{pmatrix} 1 & -\hat{\phi} \\ -1 & \phi \end{pmatrix}

(144)

Step 5. Derive the closed-form formula

\begin{pmatrix} F_{n+1} \\ F_n \end{pmatrix} = \mathbf{A}^n \begin{pmatrix} 1 \\ 0 \end{pmatrix} = \mathbf{P}\mathbf{D}^n\mathbf{P}^{-1} \begin{pmatrix} 1 \\ 0 \end{pmatrix}

(145)

The intermediate steps:

\mathbf{P}^{-1} \begin{pmatrix} 1 \\ 0 \end{pmatrix} = \frac{1}{\sqrt{5}}\begin{pmatrix} 1 \\ -1 \end{pmatrix}

(146)

\mathbf{D}^n\mathbf{P}^{-1} \begin{pmatrix} 1 \\ 0 \end{pmatrix} = \frac{1}{\sqrt{5}}\begin{pmatrix} \phi^n & 0 \\ 0 & \hat{\phi}^n \end{pmatrix}\begin{pmatrix} 1 \\ -1 \end{pmatrix} = \frac{1}{\sqrt{5}}\begin{pmatrix} \phi^n \\ -\hat{\phi}^n \end{pmatrix}

(147)

\mathbf{P}\mathbf{D}^n\mathbf{P}^{-1} \begin{pmatrix} 1 \\ 0 \end{pmatrix} = \frac{1}{\sqrt{5}}\begin{pmatrix} \phi & \hat{\phi} \\ 1 & 1 \end{pmatrix}\begin{pmatrix} \phi^n \\ -\hat{\phi}^n \end{pmatrix} = \frac{1}{\sqrt{5}}\begin{pmatrix} \phi^{n+1} - \hat{\phi}^{n+1} \\ \phi^n - \hat{\phi}^n \end{pmatrix}

(148)

Component 1 ( $F_n$ ) gives the famous Binet’s formula:

\underline{F_n = \frac{1}{\sqrt{5}}\left[\left(\frac{1+\sqrt{5}}{2}\right)^n - \left(\frac{1-\sqrt{5}}{2}\right)^n\right] = \frac{\phi^n - \hat{\phi}^n}{\sqrt{5}}}

(149)

Step 6. Checks and observations

For $n=0$ : $F_0 = \frac{\phi^0 - \hat{\phi}^0}{\sqrt{5}} = \frac{1-1}{\sqrt{5}} = 0$ ✓
For $n=1$ : $F_1 = \frac{\phi - \hat{\phi}}{\sqrt{5}} = \frac{\sqrt{5}}{\sqrt{5}} = 1$ ✓
Since $|\hat{\phi}| < 1$ , we have $\hat{\phi}^n \to 0$ for large $n$ , so $F_n \approx \frac{\phi^n}{\sqrt{5}}$

This shows that the Fibonacci sequence grows exponentially with the golden ratio as its base!

The ratio of consecutive terms tends to the golden ratio:

\lim_{n \to \infty} \frac{F_{n+1}}{F_n} = \lim_{n \to \infty} \frac{\phi^{n+1} - \hat{\phi}^{n+1}}{\phi^n - \hat{\phi}^n} = \lim_{n \to \infty} \frac{\phi^{n+1}}{\phi^n} = \phi = \frac{1+\sqrt{5}}{2}

(150)

Application 3: Simplifying the computation of matrix functions

For a function $f(x)$ , if $\mathbf{A}= \mathbf{P}\mathbf{D}\mathbf{P}^{-1}$ is diagonalizable, we can define the matrix function $f(\mathbf{A})$ :

f(\mathbf{A}) := \mathbf{P}f(\mathbf{D})\mathbf{P}^{-1} = \mathbf{P}\text{diag}(f(\lambda_0), f(\lambda_1), \ldots, f(\lambda_{n-1}))\mathbf{P}^{-1}

(151)

In particular, the matrix exponential is

e^{\mathbf{A}} = \mathbf{P}\text{diag}(e^{\lambda_0}, e^{\lambda_1}, \ldots, e^{\lambda_{n-1}})\mathbf{P}^{-1}

(152)

Exercise 4 (Computing a Matrix Exponential)

Consider the matrix $\mathbf{A} = \begin{pmatrix} 1 & 1 \\ 0 & 2 \end{pmatrix}$ .

Verify that $\mathbf{A}$ is diagonalizable, and find $\mathbf{P}$ and $\mathbf{D}$ such that $\mathbf{A} = \mathbf{P}\mathbf{D}\mathbf{P}^{-1}$ .
Use the diagonalization to compute $e^{\mathbf{A}}$ .
Verify that for a diagonal matrix $\mathbf{D} = \text{diag}(\lambda_0, \lambda_1)$ ,
$e^{\mathbf{D}} = \text{diag}(e^{\lambda_0}, e^{\lambda_1})$
(153)
◆Challenge. Compute $\sin(\mathbf{A})$ and $\cos(\mathbf{A})$ (Hint. Use the definitions of $\sin(\mathbf{D})$ and $\cos(\mathbf{D})$ .)

Application 4: Solving systems of differential equations

Diagonalization has an important application to solving systems of linear differential equations. Consider the first-order linear system

\frac{d\mathbf{x}}{dt} = \mathbf{A}\mathbf{x}

(158)

where $\mathbf{x}(t) \in \mathbb{R}^n$ is the unknown vector-valued function and $\mathbf{A}$ is an $n \times n$ constant coefficient matrix.

If $\mathbf{A}$ can be diagonalized as $\mathbf{A} = \mathbf{P}\mathbf{D}\mathbf{P}^{-1}$ , then the solution can be written as

\mathbf{x}(t) = e^{\mathbf{A}t}\mathbf{x}_0 = \mathbf{P}e^{\mathbf{D}t}\mathbf{P}^{-1}\mathbf{x}_0

(159)

where $\mathbf{x}_0 = \mathbf{x}(0)$ is the initial condition and

e^{\mathbf{D}t} = \text{diag}(e^{\lambda_0 t}, e^{\lambda_1 t}, \ldots, e^{\lambda_{n-1} t})

(160)

Example 13 (Decoupling a System of Differential Equations by a Similarity Transformation)

Consider a model of a two-tank chemical reactor. Let $x_0(t)$ and $x_1(t)$ be the concentrations of a substance in the two tanks; their rates of change are described by the coupled system of differential equations

\frac{d}{dt}\begin{pmatrix} x_0 \\ x_1 \end{pmatrix} = \underbrace{\begin{pmatrix} -3 & 1 \\ 2 & -2 \end{pmatrix}}_{\mathbf{A}} \begin{pmatrix} x_0 \\ x_1 \end{pmatrix}, \quad \begin{pmatrix} x_0(0) \\ x_1(0) \end{pmatrix} = \begin{pmatrix} 3 \\ 1 \end{pmatrix}

(161)

The difficulty of solving directly. The equation for $x_0$ contains $x_1$ , and the equation for $x_1$ contains $x_0$ ; the two are coupled and cannot be solved independently. What a similarity transformation does is find a “decoupled language.”

Step 1. Diagonalize $\mathbf{A}$

Characteristic polynomial: $p_A(\lambda) = \lambda^2 + 5\lambda + 4 = (\lambda+1)(\lambda+4)$

Eigenvalues: $\lambda_0 = -1$ , $\lambda_1 = -4$ (distinct, so the matrix is diagonalizable).

For $\lambda_0 = -1$ , solve $(\mathbf{A}+\mathbf{I})\mathbf{v} = \mathbf{0}$ :

\begin{pmatrix} -2 & 1 \\ 2 & -1 \end{pmatrix}\mathbf{v} = \mathbf{0} \quad\Rightarrow\quad \mathbf{v}_0 = \begin{pmatrix} 1 \\ 2 \end{pmatrix}

(162)

For $\lambda_1 = -4$ , solve $(\mathbf{A}+4\mathbf{I})\mathbf{v} = \mathbf{0}$ :

\begin{pmatrix} 1 & 1 \\ 2 & 2 \end{pmatrix}\mathbf{v} = \mathbf{0} \quad\Rightarrow\quad \mathbf{v}_1 = \begin{pmatrix} 1 \\ -1 \end{pmatrix}

(163)

\mathbf{P} = \begin{pmatrix} 1 & 1 \\ 2 & -1 \end{pmatrix}, \quad \mathbf{D} = \begin{pmatrix} -1 & 0 \\ 0 & -4 \end{pmatrix}, \quad \mathbf{P}^{-1} = \frac{1}{3}\begin{pmatrix} 1 & 1 \\ 2 & -1 \end{pmatrix}

(164)

(In this example $\mathbf{P}^{-1} = \frac{1}{3}\begin{pmatrix}1&1\\2&-1\end{pmatrix}$ ; verify this yourself.)

Step 2. Change variables: let $\mathbf{y}(t) = \mathbf{P}^{-1}\mathbf{x}(t)$ , that is, describe the system in the eigenvector basis

Multiply both sides of the original equation $\dot{\mathbf{x}} = \mathbf{A}\mathbf{x}$ on the left by $\mathbf{P}^{-1}$ :

\mathbf{P}^{-1}\dot{\mathbf{x}} = \mathbf{P}^{-1}\mathbf{A}\mathbf{x} = \underbrace{\mathbf{P}^{-1}\mathbf{A}\mathbf{P}}_{\mathbf{D}}\,\underbrace{\mathbf{P}^{-1}\mathbf{x}}_{\mathbf{y}}

(165)

That is, $\dot{\mathbf{y}} = \mathbf{D}\mathbf{y}$ , which expands into two independent scalar equations:

\dot{y}_0 = -y_0, \qquad \dot{y}_1 = -4\,y_1

(166)

The coupling has disappeared completely!

Step 3. Solve the decoupled equations

y_0(t) = y_0(0)\,e^{-t}, \qquad y_1(t) = y_1(0)\,e^{-4t}

(167)

Transform the initial condition:

\mathbf{y}(0) = \mathbf{P}^{-1}\mathbf{x}(0) = \frac{1}{3}\begin{pmatrix}1&1\\2&-1\end{pmatrix}\begin{pmatrix}3\\1\end{pmatrix} = \frac{1}{3}\begin{pmatrix}4\\5\end{pmatrix}

(168)

So $y_0(0) = \tfrac{4}{3}$ and $y_1(0) = \tfrac{5}{3}$ .

Step 4. Return to the original variables, $\mathbf{x}(t) = \mathbf{P}\mathbf{y}(t)$

\begin{pmatrix} x_0(t) \\ x_1(t) \end{pmatrix} = \begin{pmatrix} 1 & 1 \\ 2 & -1 \end{pmatrix} \begin{pmatrix} \tfrac{4}{3}e^{-t} \\ \tfrac{5}{3}e^{-4t} \end{pmatrix} = \begin{pmatrix} \tfrac{4}{3}e^{-t} + \tfrac{5}{3}e^{-4t} \\[4pt] \tfrac{8}{3}e^{-t} - \tfrac{5}{3}e^{-4t} \end{pmatrix}

(169)

Physical interpretation. The solution is a superposition of two “pure decay modes”: the $e^{-t}$ mode decays along the eigendirection $\mathbf{v}_0 = (1,2)^\top$ , and the $e^{-4t}$ mode decays faster along $\mathbf{v}_1 = (1,-1)^\top$ . What the similarity transformation does is precisely to decompose “mixed dynamics” into “independent decays.”

Exercise 5 (◆A Radioactive Decay Chain)

This exercise requires basic knowledge of ordinary differential equations. If you have not yet studied differential equations, you may skip it and come back after taking the relevant course.

Consider the decay of two radioactive isotopes A and B:

A decays into B with rate constant $k_A = 2$
B decays into the stable nuclide C with rate constant $k_B = 1$

Let $x_0(t)$ and $x_1(t)$ denote the relative amounts of A and B at time $t$ . The decay process is then described by the system

\frac{d}{dt}\begin{pmatrix} x_0 \\ x_1 \end{pmatrix} = \begin{pmatrix} -2 & 0 \\ 2 & -1 \end{pmatrix} \begin{pmatrix} x_0 \\ x_1 \end{pmatrix}

(170)

The initial condition is $\mathbf{x}(0) = \begin{pmatrix} 1 \\ 0 \end{pmatrix}$ (initially there is only A and no B).

Find the eigenvalues and eigenvectors of the matrix $\mathbf{A} = \begin{pmatrix} -2 & 0 \\ 2 & -1 \end{pmatrix}$ .
Verify that $\mathbf{A}$ is diagonalizable, and find $\mathbf{P}$ and $\mathbf{D}$ .
Compute $e^{\mathbf{A}t}$ .
Solve $\mathbf{x}(t) = e^{\mathbf{A}t}\mathbf{x}_0$ to obtain explicit expressions for $x_0(t)$ and $x_1(t)$ .
Explain the physical meaning of the result:
- How should $x_0(t)$ change over time?
- Does $x_1(t)$ first increase and then decrease? Why?
- As $t \to \infty$ , what values do $x_0(t)$ and $x_1(t)$ approach?
Food for thought. If $k_A = k_B$ (that is, the two decay rates are equal), is the matrix still diagonalizable? This question is discussed in Section 8.4.

Solution to Exercise 5

1. $\mathbf{A}$ is lower triangular, so its eigenvalues are its main diagonal entries: $\lambda_0 = -2$ and $\lambda_1 = -1$ .

$\lambda_0 = -2$ : solve $(\mathbf{A} + 2\mathbf{I})\mathbf{v} = \mathbf{0}$ , that is, $2v_0 + v_1 = 0$ ; take $\mathbf{v}_0 = \begin{pmatrix} 1 \\ -2 \end{pmatrix}$
$\lambda_1 = -1$ : solve $(\mathbf{A} + \mathbf{I})\mathbf{v} = \mathbf{0}$ , that is, $v_0 = 0$ ; take $\mathbf{v}_1 = \begin{pmatrix} 0 \\ 1 \end{pmatrix}$

2. The two eigenvalues are distinct, so $\mathbf{A}$ is diagonalizable:

\mathbf{P} = \begin{pmatrix} 1 & 0 \\ -2 & 1 \end{pmatrix}, \qquad \mathbf{D} = \begin{pmatrix} -2 & 0 \\ 0 & -1 \end{pmatrix}, \qquad \mathbf{P}^{-1} = \begin{pmatrix} 1 & 0 \\ 2 & 1 \end{pmatrix}

(171)

3.

e^{\mathbf{A}t} = \mathbf{P}\begin{pmatrix} e^{-2t} & 0 \\ 0 & e^{-t} \end{pmatrix}\mathbf{P}^{-1} = \begin{pmatrix} e^{-2t} & 0 \\ 2e^{-t} - 2e^{-2t} & e^{-t} \end{pmatrix}

(172)

4.

\mathbf{x}(t) = e^{\mathbf{A}t}\begin{pmatrix} 1 \\ 0 \end{pmatrix} = \begin{pmatrix} e^{-2t} \\ 2e^{-t} - 2e^{-2t} \end{pmatrix}, \qquad x_0(t) = e^{-2t}, \quad x_1(t) = 2e^{-t} - 2e^{-2t}

(173)

5. Physical meaning:

$x_0(t) = e^{-2t}$ decays exponentially and monotonically—A only decays and is never replenished.
$x_1(t)$ first increases and then decreases: at first A decays in large amounts and replenishes B (production outpaces consumption); later A is gradually exhausted, and B is left only with its own decay. Setting $x_1'(t) = -2e^{-t} + 4e^{-2t} = 0$ gives the peak time $t = \ln 2$ , with peak value $x_1(\ln 2) = \frac{1}{2}$ .
As $t \to \infty$ , $x_0(t), x_1(t) \to 0$ : every nuclide eventually decays into the stable nuclide C.

6. If $k_A = k_B = k$ , then $\mathbf{A} = \begin{pmatrix} -k & 0 \\ k & -k \end{pmatrix}$ has only one eigenvalue $\lambda = -k$ , with algebraic multiplicity $a_\lambda = 2$ ; but $\mathbf{A} + k\mathbf{I} = \begin{pmatrix} 0 & 0 \\ k & 0 \end{pmatrix}$ has rank 1, so the geometric multiplicity is $g_\lambda = 1 < a_\lambda$ , and the matrix is not diagonalizable—this is exactly the situation in which the Jordan canonical form of §8.4 enters, and the solution will contain “polynomial × exponential” terms of the form $t e^{-kt}$ .

§8.3 Summary

This section answered a global question: can we find a basis in which the matrix representation of a given linear transformation is as simple as possible?

The similarity transformation $\mathbf{B} = \mathbf{P}^{-1}\mathbf{A}\mathbf{P}$ is the mathematical language for “looking at the same thing from another angle.” Similar matrices share the characteristic polynomial, the eigenvalues, the trace, the determinant, and the rank; these similarity invariants are intrinsic attributes of the linear transformation itself and do not change with the choice of basis.

Diagonalization is the best case of a similarity transformation: when $\mathbf{A}$ has $n$ linearly independent eigenvectors, taking them as the basis reduces the transformation to independent scaling along each direction, and the effects of rotation and shear disappear completely. The necessary and sufficient condition is that the characteristic polynomial splits completely over the field in use (which always holds over $\mathbb{C}$ ) and that $g_\lambda = a_\lambda$ for every eigenvalue; pairwise distinct eigenvalues form a sufficient but not a necessary condition.

The power of diagonalization shows fully in its applications. The matrix power $\mathbf{A}^k = \mathbf{P}\mathbf{D}^k\mathbf{P}^{-1}$ reduces to entrywise scalar powers; the solutions of linear difference equations and systems of differential equations can be written explicitly in terms of exponential functions of the eigenvalues—Binet’s formula behind the closed form of the Fibonacci sequence is the most elegant demonstration of this idea. The non-diagonalizable case is left to §8.4.

8.4 The Jordan Canonical Form¶

Not every matrix is diagonalizable. When a matrix is not diagonalizable, we need a more general canonical form—the Jordan canonical form.

8.4.1 Overview of the Jordan Canonical Form¶

In Section 8.2.2 we saw examples of non-diagonalizable matrices (Examples Example 4 and Example 5). What these matrices have in common is that some eigenvalue has geometric multiplicity smaller than its algebraic multiplicity, that is, $g_\lambda < a_\lambda$ . This means that the eigenvalue does not have enough linearly independent eigenvectors to form a complete diagonalizing basis.

Jordan blocks: intuition and cause

To understand the alternative when diagonalization is impossible, we return to the “pirate ship” analogy of Example Example 5:

Repeated application of the nilpotent matrix $\mathbf{A}$ makes vectors “fall level by level”; mathematically, this corresponds to a chain of generalized eigenvectors.

When diagonalization fails, the operator $(\mathbf{A} - \lambda\mathbf{I})$ induces this kind of chained filtration structure in the space (because of the minimal polynomial, the chain cannot grow indefinitely):

\{\mathbf{0}\} \subset \ker(\mathbf{A} - \lambda\mathbf{I}) \subset \ker(\mathbf{A} - \lambda\mathbf{I})^2 \subset \cdots \subset \ker(\mathbf{A} - \lambda\mathbf{I})^k

(174)

Merely extending the basis layer by layer along the kernels is not enough; we must choose a chain basis satisfying $(\mathbf{A}-\lambda\mathbf{I})\mathbf{v}_0=\mathbf{0}$ and $(\mathbf{A}-\lambda\mathbf{I})\mathbf{v}_j=\mathbf{v}_{j-1}$ . The same eigenvalue may require several chains, and each chain of length $k$ corresponds to a Jordan block $\mathbf{J}_k(\lambda)$ :

\mathbf{J}_k(\lambda) = \begin{pmatrix} \lambda & 1 & 0 & \cdots & 0 \\ 0 & \lambda & 1 & \cdots & 0 \\ 0 & 0 & \lambda & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & \lambda \end{pmatrix}_{k \times k}

(175)

Main diagonal: every entry is the eigenvalue $\lambda$ .
Superdiagonal: every entry is 1; this is precisely the mathematical trace of vectors being “pushed level by level” along the kernel chain.
Structural essence: a Jordan block captures exactly the property of being “almost diagonalizable”: it differs from the diagonal matrix $\lambda\mathbf{I}_k$ only by a nilpotent matrix consisting entirely of the 1s on the superdiagonal.

The Jordan canonical form theorem

Since each eigenvalue breaks down into a number of Jordan blocks according to its kernel chains, assembling these blocks along the diagonal like pieces of a jigsaw puzzle produces the ultimate form in the similarity classification of arbitrary matrices.

Theoretical value and numerical limitations

Although the Jordan canonical form is perfect in theory, there is a huge gap between its algebraic structure and numerical computation:

A perfect answer about algebraic structure:
Geometric multiplicity $g_\lambda$ : corresponds to the number of Jordan blocks for the eigenvalue $\lambda$ (that is, $\dim\ker(\mathbf{A}-\lambda\mathbf{I})$ ).
Algebraic multiplicity $a_\lambda$ : corresponds to the total size of all the Jordan blocks for the eigenvalue $\lambda$ .
Exponent of the factor $(t-\lambda)$ in the minimal polynomial: corresponds to the size of the largest Jordan block for the eigenvalue $\lambda$ .
A numerical disaster:
Extremely unstable structure: how the Jordan blocks split is extremely sensitive to tiny perturbations of the matrix entries. For example, $\begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix}$ has only one block, but adding just an $\epsilon$ in the lower right corner immediately makes the matrix diagonalizable.
High computational cost: determining exactly the order at which the generalized eigenspaces stabilize (that is, the end of the kernel chain) involves high powers of matrices, and numerical errors grow rapidly. In practical numerical ODE computations, other, more stable algorithms are used.

Where this topic sits in the course For these reasons, the numerical procedure of computing generalized eigenvector chains by hand lies outside the practical scope of this course, but the theoretical framework of the existence and uniqueness of the Jordan canonical form is worth understanding. Our main purpose in studying it is to understand how it brings the theory of matrix similarity to a perfect close. In practical numerical applications we turn to more stable tools, such as the singular value decomposition (SVD) of Chapter 11.

The Concavity Theorem for Kernel Dimension Increments¶

Underlying the theory of the Jordan canonical form is a geometrically very elegant property. For a linear map $N: V \to V$ , the layer-by-layer dimension increments of its kernel chain

\{\mathbf{0}\} \subset \ker(N) \subset \ker(N^2) \subset \ker(N^3) \subset \cdots

(178)

namely

\underbrace{\dim\ker(N) - 0}_{d_1} \;\ge\; \underbrace{\dim\ker(N^2) - \dim\ker(N)}_{d_2} \;\ge\; \underbrace{\dim\ker(N^3) - \dim\ker(N^2)}_{d_3} \;\ge\; \cdots

(179)

are monotonically nonincreasing. This concavity is the behind-the-scenes guarantee that the distribution of Jordan block sizes is unique.

Proof 7

Without using the theory of quotient spaces, we use the Steinitz exchange lemma (see Theorem Theorem 2): if a linear space $W$ contains a linearly independent set $L$ , and $S$ is a spanning set of $W$ , then $|L| \le |S|$ .

Below we construct, in the space $W = \ker(N^k)$ , a linearly independent set and a spanning set, and then apply the lemma to obtain the conclusion.

Notation. Let
$\textcolor{#D62728}{m} = d_k = \dim\ker(N^k) - \dim\ker(N^{k-1})$
(182)
$\textcolor{#57068C}{n} = d_{k+1} = \dim\ker(N^{k+1}) - \dim\ker(N^k)$
(183)
The goal is to prove $\textcolor{#57068C}{n} \le \textcolor{#D62728}{m}$ .
Construct the linearly independent set $L$ (from the higher layer). Since $d_{k+1} = \textcolor{#57068C}{n}$ , $\ker(N^{k+1})$ has exactly $\textcolor{#57068C}{n}$ more dimensions than $\ker(N^k)$ . Let $\{\mathbf{u}_0, \dots, \mathbf{u}_{a+m-1}\}$ be a basis of $\ker(N^k)$ , and extend it to a basis of $\ker(N^{k+1})$ :
$\{\mathbf{u}_0, \dots, \mathbf{u}_{a+m-1},\; \textcolor{#57068C}{\mathbf{w}_0, \dots, \mathbf{w}_{n-1}}\}$
(184)
where $\textcolor{#57068C}{\mathbf{w}_0, \dots, \mathbf{w}_{n-1}} \in \ker(N^{k+1}) \setminus \ker(N^k)$ are the $\textcolor{#57068C}{n}$ new basis vectors added in the extension. By the construction of the basis extension, these $\textcolor{#57068C}{n}$ vectors are linearly independent modulo $\ker(N^k)$ .
Apply $N^k$ to these $\textcolor{#57068C}{n}$ vectors: since $\mathbf{w}_i ∉ \ker(N^k)$ , we have $N^k(\mathbf{w}_i) \ne \mathbf{0}$ . And since $\mathbf{w}_i \in \ker(N^{k+1})$ , we have $N(N^k(\mathbf{w}_i)) = N^{k+1}(\mathbf{w}_i) = \mathbf{0}$ , so $N^k(\mathbf{w}_i) \in \ker(N) \subset \ker(N^k)$ .
If $\sum_i c_i\, N^k(\mathbf{w}_i) = \mathbf{0}$ , then $N^k\!\left(\sum_i c_i \mathbf{w}_i\right) = \mathbf{0}$ , that is, $\sum_i c_i \mathbf{w}_i \in \ker(N^k)$ . But the $\mathbf{w}_i$ are the part of the basis of $\ker(N^{k+1})$ lying outside $\ker(N^k)$ , so $\sum_i c_i \mathbf{w}_i \in \ker(N^k)$ if and only if $c_i = 0$ for all $i$ .
Hence $\{N^k(\mathbf{w}_i)\}$ is linearly independent. Apply $N$ to each $\mathbf{w}_i$ and let
$L = \{N(\mathbf{w}_0),\, \dots,\, N(\mathbf{w}_{n-1})\}$
(185)
- Location. Since $N^k(N(\mathbf{w}_i)) = N^{k+1}(\mathbf{w}_i) = \mathbf{0}$ , we have $L \subset \ker(N^k)$ . And since $\mathbf{w}_i ∉ \ker(N^k)$ , we have $N(\mathbf{w}_i) ∉ \ker(N^{k-1})$ , that is, $L \cap \ker(N^{k-1}) = \varnothing$ .
$L$ together with a basis of $\ker(N^{k-1})$ is linearly independent. Let $\{\mathbf{u}_0, \dots, \mathbf{u}_{a-1}\}$ be a basis of $\ker(N^{k-1})$ , and consider the linear combination
$\sum_{i} c_i\, N(\mathbf{w}_i) + \sum_{j} b_j\, \mathbf{u}_j = \mathbf{0}$
(186)
Apply $N^{k-1}$ to both sides: the right side is $\mathbf{0}$ ; the second term on the left vanishes because $\mathbf{u}_j \in \ker(N^{k-1})$ ; the first term becomes $\sum_i c_i\, N^k(\mathbf{w}_i)$ . Since $\mathbf{w}_i ∉ \ker(N^k)$ , we have $N^k(\mathbf{w}_i) \ne \mathbf{0}$ , and since the $\mathbf{w}_i$ were chosen so that $\{N^k(\mathbf{w}_i)\}$ is linearly independent, it follows that $c_i = 0$ for all $i$ . Substituting back gives $\sum_j b_j \mathbf{u}_j = \mathbf{0}$ , and the linear independence of the basis gives $b_j = 0$ for all $j$ .
Hence $L \cup \{\mathbf{u}_0, \dots, \mathbf{u}_{a-1}\}$ is a linearly independent set in $\ker(N^k)$ with $\textcolor{#57068C}{n} + a$ vectors.
Compare dimensions. The basis of $\ker(N^k)$ is extended from the basis of $\ker(N^{k-1})$ ; it contains $\textcolor{#D62728}{m} + a$ vectors and spans $\ker(N^k)$ . A linearly independent set is no larger than a spanning set (the Steinitz lemma), so
$\textcolor{#57068C}{n} + a \le \textcolor{#D62728}{m} + a$
(187)
Cancelling $a$ gives $\underline{\textcolor{#57068C}{n} \le \textcolor{#D62728}{m}}$ , that is, $d_{k+1} \le d_k$ , which proves concavity. $\square$

Remark 3 (Geometric Intuition for Concavity)

Geometrically, the concavity theorem says: each step deeper into the kernel chain adds fewer and fewer new dimensions; it never suddenly expands.

If we plot the power $k$ on the horizontal axis and $\dim\ker(N^k)$ on the vertical axis, the growth curve of $\ker(N^k)$ is an upward-convex polygonal line—its slope (that is, $d_k$ ) decreases monotonically to zero.

\dim\ker(N^0) = 0 \;\le\; \dim\ker(N^1) \;\le\; \dim\ker(N^2) \;\le\; \cdots

(188)

For each eigenvalue $\lambda$ , the sequence of Jordan block sizes in the Jordan canonical form is completely determined by the shape of this convex polygonal line.

Example 14 (Reading Off the Jordan Structure from the Kernel Chain)

Suppose the characteristic and minimal polynomials of a linear operator $T: \mathbb{R}^4 \to \mathbb{R}^4$ are

p_T(\lambda) = (\lambda - 2)^4, \qquad m_T(\lambda) = (\lambda - 2)^2.

(189)

Rather than applying the Jordan canonical form theorem directly, we use the kernel chain of a nilpotent operator to derive all possible Jordan structures of $T$ from the point of view of dimensions.

Step 1. Reduce to a nilpotent operator

Let $N = T - 2I$ . From the given conditions:

$p_T(\lambda) = (\lambda-2)^4$ and the Cayley–Hamilton theorem $\Rightarrow$ $N^4 = \mathbf{0}$ , and hence $\ker(N^4) = \mathbb{R}^4$ .
$m_T(\lambda) = (\lambda-2)^2$ $\Rightarrow$ $N^2 = \mathbf{0}$ and $N \neq \mathbf{0}$ .

Hence $N$ is a nilpotent operator with $N^2 = \mathbf{0}$ and $N \neq \mathbf{0}$ .

Step 2. Set up the kernel chain

Since $N^2 = \mathbf{0}$ , the kernel chain must stabilize at the second step:

\{\mathbf{0}\} \subset \ker(N) \subset \ker(N^2) = \mathbb{R}^4.

(190)

Let $d_1 = \dim\ker(N)$ . By the rank-nullity theorem, $\dim\operatorname{Im}(N) = 4 - d_1$ .

Step 3. Use the concavity of dimensions to constrain $d_1$

The dimension increments of the kernel chain are concave (the increment never increases from one step to the next):

\underbrace{d_1 - 0}_{\text{first increment}} \;\geq\; \underbrace{4 - d_1}_{\text{second increment}} \implies 2d_1 \geq 4 \implies d_1 \geq 2.

(191)

Moreover, $N \neq \mathbf{0}$ gives $d_1 < 4$ , and $N$ is not invertible, so $d_1 > 0$ . Altogether:

d_1 \in \{2, 3\}.

(192)

Step 4. Case 1: $d_1 = 2$ (two $2\times 2$ Jordan blocks)

Here $\dim\ker(N) = \dim\operatorname{Im}(N) = 2$ , and $N^2 = \mathbf{0}$ guarantees $\operatorname{Im}(N) \subset \ker(N)$ . Since the dimensions are equal,

\operatorname{Im}(N) = \ker(N).

(193)

Choose two vectors $\mathbf{v}_0, \mathbf{v}_2$ in $\mathbb{R}^4$ , not in $\ker(N)$ , that are linearly independent modulo $\ker(N)$ , and let $\mathbf{v}_1 = N\mathbf{v}_0$ and $\mathbf{v}_3 = N\mathbf{v}_2$ . From $N^2 = \mathbf{0}$ we know $\mathbf{v}_1, \mathbf{v}_3 \in \ker(N)$ , and $\{\mathbf{v}_0, \mathbf{v}_1, \mathbf{v}_2, \mathbf{v}_3\}$ is a basis of $\mathbb{R}^4$ ; in this basis the action of $N$ consists of two independent chains of generalized eigenvectors:

\mathbf{v}_0 \overset{\mathbf{N}}{\rightarrow} \mathbf{v}_1 \overset{\mathbf{N}}{\rightarrow} \mathbf{0}, \qquad \mathbf{v}_2 \overset{\mathbf{N}}{\rightarrow} \mathbf{v}_3 \overset{\mathbf{N}}{\rightarrow} \mathbf{0}.

(194)

Restoring $T = N + 2I$ , the corresponding Jordan canonical form is the direct sum of two $2\times 2$ blocks:

\mathbf{J} = \mathbf{J}_2(2) \oplus \mathbf{J}_2(2) = \begin{pmatrix} 2 & 1 & 0 & 0 \\ 0 & 2 & 0 & 0 \\ 0 & 0 & 2 & 1 \\ 0 & 0 & 0 & 2 \end{pmatrix}.

(195)

Here the geometric multiplicity is $g_2 = \dim\ker(N) = 2$ (the number of Jordan blocks), and the largest block has size 2 (in agreement with $\deg m_T = 2$ ).

Step 5. Case 2: $d_1 = 3$ (one $2\times 2$ block plus two $1\times 1$ blocks)

Here $\dim\operatorname{Im}(N) = 1$ , so the range is spanned by a single vector. Choose $\mathbf{v}_0 ∉ \ker(N)$ so that $\mathbf{v}_1 = N\mathbf{v}_0 \neq \mathbf{0}$ ; then $\mathbf{v}_1$ spans all of $\operatorname{Im}(N)$ , and $N\mathbf{v}_1 = \mathbf{0}$ , that is, $\mathbf{v}_1 \in \ker(N)$ .

In $\ker(N)$ (which is three-dimensional), add two vectors $\mathbf{v}_2, \mathbf{v}_3$ so that $\{\mathbf{v}_1, \mathbf{v}_2, \mathbf{v}_3\}$ is a basis of $\ker(N)$ . The action of $N$ is then:

\mathbf{v}_0 \overset{\mathbf{N}}{\rightarrow} \mathbf{v}_1 \overset{\mathbf{N}}{\rightarrow} \mathbf{0}, \qquad \mathbf{v}_2 \overset{\mathbf{N}}{\rightarrow} \mathbf{0}, \qquad \mathbf{v}_3 \overset{\mathbf{N}}{\rightarrow} \mathbf{0}.

(196)

Restoring $T = N + 2I$ , the corresponding Jordan canonical form is the direct sum of one $2\times 2$ block and two $1\times 1$ blocks:

\mathbf{J} = \mathbf{J}_2(2) \oplus \mathbf{J}_1(2) \oplus \mathbf{J}_1(2) = \begin{pmatrix} 2 & 1 & 0 & 0 \\ 0 & 2 & 0 & 0 \\ 0 & 0 & 2 & 0 \\ 0 & 0 & 0 & 2 \end{pmatrix}.

(197)

Here the geometric multiplicity is $g_2 = 3$ , and the largest block still has size 2 (consistent with $\deg m_T = 2$ ).

Step 6. Conclusion

Given the characteristic polynomial $(\lambda-2)^4$ and the minimal polynomial $(\lambda-2)^2$ , the Jordan structure of $T$ has exactly two possibilities:

$\dim\ker(N)$	Partition into Jordan blocks	Jordan canonical form
2	$2+2$	$\mathbf{J}_2(2)\oplus\mathbf{J}_2(2)$
3	$2+1+1$	$\mathbf{J}_2(2)\oplus\mathbf{J}_1(2)\oplus\mathbf{J}_1(2)$

Both cases are compatible with $m_T(\lambda) = (\lambda-2)^2$ (the largest block has size 2 in each), but their geometric multiplicities differ. To tell them apart, one must actually compute $\dim\ker(T - 2I)$ .

Key insight. The Jordan structure of a nilpotent operator is completely determined by how the dimension increments are distributed across the layers $\ker(N^k)$ of the kernel chain. The concavity of the sequence of increments is the geometric root of the decreasing arrangement of Jordan block sizes, and the degree of the minimal polynomial corresponds to the length of the longest chain of generalized eigenvectors.

8.4.2 Applications of the Jordan Canonical Form¶

Although computing the Jordan canonical form is beyond the scope of this course, the application formulas for 2×2 Jordan blocks are very important and relatively simple. These formulas are indispensable for solving difference equations and differential equations with non-diagonalizable matrices. For the complete derivations, visual comparisons, and application exercises, see Experiment 7.

Recall the form of a 2×2 Jordan block:

\mathbf{J}_2(\lambda) = \begin{pmatrix} \lambda & 1 \\ 0 & \lambda \end{pmatrix}

(198)

Key observation: in the non-diagonalizable case, the general solution admits terms of the “polynomial × exponential” type (for particular initial values the coefficients of the polynomial terms may vanish); this is the essential difference from the diagonalizable case.

Application 1: Computing powers of a matrix¶

The formula for powers of a Jordan block is

\mathbf{J}_2(\lambda)^n = \begin{pmatrix} \lambda^n & n\lambda^{n-1} \\ 0 & \lambda^n \end{pmatrix}

(199)

This power formula is stated for $\lambda\neq0$ . In addition, $\mathbf{J}^0=\mathbf{I}$ ; when $\mathbf{J}_k(0)=\mathbf{N}$ , use $\mathbf{J}_k(0)^n=\mathbf{N}^n$ directly, which is zero for $n\ge k$ ; this avoids negative powers and the ambiguity of 0⁰.

Comparison with the diagonalizable case.

For the diagonal matrix, $\begin{pmatrix} \lambda & 0 \\ 0 & \lambda \end{pmatrix}^n = \begin{pmatrix} \lambda^n & 0 \\ 0 & \lambda^n \end{pmatrix}$ ; the structure is unchanged by taking powers.
For the Jordan block, the term $n\lambda^{n-1}$ appears in the upper right corner, that is, a “polynomial × exponential” term.

Where does this extra term in the upper right corner come from? It comes from the fact that the binomial expansion of $\mathbf{J}_2(\lambda) = \lambda\mathbf{I} + \mathbf{N}$ (where the nilpotent matrix satisfies $\mathbf{N}^2 = \mathbf{0}$ ) terminates after finitely many terms—for the derivation and a numerical check, see Experiment 7, §1.1.

Application 2: Solving difference equations¶

Consider the system of linear difference equations $\mathbf{x}_{n+1} = \mathbf{A}\mathbf{x}_n$ .

If $\mathbf{A}$ is similar to a Jordan canonical form, $\mathbf{A} = \mathbf{P}\mathbf{J}\mathbf{P}^{-1}$ , then

\mathbf{x}_n = \mathbf{A}^n\mathbf{x}_0 = \mathbf{P}\mathbf{J}^n\mathbf{P}^{-1}\mathbf{x}_0

(200)

By the power formula of Application 1, when $\lambda\neq0$ the general solution admits terms of the form $n\lambda^{n-1}$ , that is, a polynomial ( $n$ ) times an exponential ( $\lambda^n$ ); for particular initial values (for example, ones lying in the eigenspace) their coefficients may be zero. For a complete comparison of the trajectories in the diagonalizable case ( $g_\lambda = a_\lambda$ ) and the non-diagonalizable case ( $g_\lambda < a_\lambda$ ), see Part 1 of Experiment 7.

Exercise 6 (A Second-Order Difference Equation with Constant Coefficients)

Consider the difference equation

a_n = 2r a_{n-1} - r^2 a_{n-2}, \quad n \geq 2

(201)

where $r \neq 0$ is a constant and the initial conditions are $a_0 = 0, a_1 = 1$ .

Write this equation in the matrix form $\mathbf{x}_n = \mathbf{A}\mathbf{x}_{n-1}$ , where $\mathbf{x}_n = \begin{pmatrix} a_n \\ a_{n-1} \end{pmatrix}$ .
Compute the eigenvalues of the matrix $\mathbf{A}$ .
Verify that $\mathbf{A}$ is not diagonalizable.
Find the Jordan canonical form of $\mathbf{A}$ .
Find the closed-form formula for $a_n$ .
Verify $a_0 = 0, a_1 = 1$ , and compute $a_2, a_3, a_4$ .

Solution to Exercise 6

1. Take $\mathbf{x}_n = \begin{pmatrix} a_n \\ a_{n-1} \end{pmatrix}$ ; the recurrence then becomes

\mathbf{x}_n = \mathbf{A}\mathbf{x}_{n-1}, \qquad \mathbf{A} = \begin{pmatrix} 2r & -r^2 \\ 1 & 0 \end{pmatrix}, \qquad \mathbf{x}_1 = \begin{pmatrix} a_1 \\ a_0 \end{pmatrix} = \begin{pmatrix} 1 \\ 0 \end{pmatrix}

(202)

2. $p_A(\lambda) = \det(\mathbf{A} - \lambda\mathbf{I}) = \lambda^2 - 2r\lambda + r^2 = (\lambda - r)^2$ , so the only eigenvalue is $\lambda = r$ , with algebraic multiplicity $a_r = 2$ .

3. $\mathbf{A} - r\mathbf{I} = \begin{pmatrix} r & -r^2 \\ 1 & -r \end{pmatrix}$ ; since $r \neq 0$ , its two rows are proportional and its rank is 1, so the geometric multiplicity is $g_r = \dim\ker(\mathbf{A} - r\mathbf{I}) = 2 - 1 = 1 < 2 = a_r$ , and $\mathbf{A}$ is not diagonalizable.

4. Take the eigenvector $\mathbf{v}_0 = \begin{pmatrix} r \\ 1 \end{pmatrix}$ ; solving the generalized eigenvector equation $(\mathbf{A} - r\mathbf{I})\mathbf{v}_1 = \mathbf{v}_0$ gives $\mathbf{v}_1 = \begin{pmatrix} 1 \\ 0 \end{pmatrix}$ . Let $\mathbf{P} = \begin{pmatrix} \mathbf{v}_0 & \mathbf{v}_1 \end{pmatrix}$ ; a direct check gives

\mathbf{P} = \begin{pmatrix} r & 1 \\ 1 & 0 \end{pmatrix}, \qquad \mathbf{P}^{-1}\mathbf{A}\mathbf{P} = \mathbf{J} = \begin{pmatrix} r & 1 \\ 0 & r \end{pmatrix}

(203)

5. By the structure of powers of a $2 \times 2$ Jordan block (Application 1), the general solution has the “polynomial × exponential” form $a_n = (C_0 + C_1 n)\,r^n$ . Substituting the initial conditions: $a_0 = C_0 = 0$ , and $a_1 = C_1 r = 1$ gives $C_1 = \frac{1}{r}$ , so

a_n = n\,r^{n-1}

(204)

6. Check: $a_0 = 0$ and $a_1 = 1 \cdot r^0 = 1$ hold; moreover $a_2 = 2r$ , $a_3 = 3r^2$ , $a_4 = 4r^3$ . Substituting back into the recurrence as a check:

2r\,a_{n-1} - r^2 a_{n-2} = 2(n-1)r^{n-1} - (n-2)r^{n-1} = n\,r^{n-1} = a_n \;\checkmark

(205)

Application 3: The matrix exponential¶

The matrix exponential of a 2×2 Jordan block is

e^{\mathbf{J}_2(\lambda)t} = e^{\lambda t}\begin{pmatrix} 1 & t \\ 0 & 1 \end{pmatrix}

(206)

This formula for the matrix exponential holds for all $\lambda$ (including $\lambda=0$ ) and is not subject to the restriction $\lambda\neq0$ on the discrete powers in Application 1.

Comparison with the diagonalizable case.

For the diagonal matrix, $e^{\operatorname{diag}(\lambda,\lambda)t} = \begin{pmatrix} e^{\lambda t} & 0 \\ 0 & e^{\lambda t} \end{pmatrix}$ ; the upper right corner is zero.
For the Jordan block, the term $t\,e^{\lambda t}$ appears in the upper right corner, which is how “polynomial × exponential” shows up in continuous time.

This formula is the core tool for solving systems of differential equations with non-diagonalizable matrices. The $t$ in the upper right corner comes from the fact that the Taylor expansion $e^{\mathbf{N}t} = \mathbf{I} + \mathbf{N}t$ of the nilpotent part $\mathbf{N}$ ( $\mathbf{N}^2=\mathbf{0}$ ) terminates after finitely many terms—for the derivation, see Part 3 of Experiment 7.

Application 4: Solving systems of differential equations¶

Consider the first-order linear system

\frac{d\mathbf{x}}{dt} = \mathbf{A}\mathbf{x}

(207)

If $\mathbf{A} = \mathbf{P}\mathbf{J}\mathbf{P}^{-1}$ , then the solution is

\mathbf{x}(t) = e^{\mathbf{A}t}\mathbf{x}_0 = \mathbf{P}e^{\mathbf{J}t}\mathbf{P}^{-1}\mathbf{x}_0

(208)

By the matrix exponential formula of Application 3, the general solution admits terms of the form $te^{\lambda t}$ , that is, a polynomial ( $t$ ) times an exponential ( $e^{\lambda t}$ ). For the complete computations in the two cases—decoupling by a similarity transformation ( $g_\lambda = a_\lambda$ ) and a Jordan block ( $g_\lambda < a_\lambda$ )—together with phase-space trajectories and the application to radioactive decay chains, see Part 2 of Experiment 7 and its Exercises 1 and 2.

Generalization to Jordan Blocks of Any Size¶

The formulas above generalize to Jordan blocks of any size. For a $k \times k$ Jordan block

\mathbf{J}_k(\lambda) = \begin{pmatrix} \lambda & 1 & & \\ & \lambda & \ddots & \\ & & \ddots & 1 \\ & & & \lambda \end{pmatrix}

(209)

the general solution admits polynomial factors of degree up to $k-1$ ; for particular initial values the coefficients may be zero:

Difference equations: when $\lambda\neq0$ , the general solution admits $n^j \lambda^n$ , $j = 0, 1, \ldots, k-1$
Differential equations: the general solution admits $t^j e^{\lambda t}$ , $j = 0, 1, \ldots, k-1$

The 2×2 case ( $k=2$ ) consists of exactly the two terms $j=0, 1$ , that is, a constant term and a linear term. Larger Jordan blocks admit higher-degree polynomials, but the essential “polynomial × exponential” structure is unchanged.

In practice, Jordan blocks of high dimension are computed with a computer algebra system (such as Python’s sympy.jordan_form()) rather than by hand.

The four applications side by side (2×2 Jordan block; the polynomial-times-exponential form for difference equations is restricted to $\lambda\neq0$ )

	Diagonalizable ( $g_\lambda = a_\lambda$ )	Jordan block ( $g_\lambda < a_\lambda$ )
$\mathbf{A}^n$	$\mathbf{P}\operatorname{diag}(\lambda^n)\mathbf{P}^{-1}$	$n\lambda^{n-1}$ appears in the upper right corner (Application 1)
Difference equation $\mathbf{x}_n$	superposition of pure exponentials $\lambda^n$	terms of the form $n\lambda^n$ appear (Application 2)
$e^{\mathbf{A}t}$	$\mathbf{P}\operatorname{diag}(e^{\lambda t})\mathbf{P}^{-1}$	$te^{\lambda t}$ appears in the upper right corner (Application 3)
Differential equation $\mathbf{x}(t)$	superposition of pure exponentials $e^{\lambda t}$	terms of the form $te^{\lambda t}$ appear (Application 4)

The common root: the 2×2 Jordan block here equals $\lambda\mathbf{I} + \mathbf{N}$ with $\mathbf{N}^2 = \mathbf{0}$ , so the binomial expansion or the Taylor expansion terminates after the linear term, which admits a polynomial factor; particular initial values can make its coefficient zero.

§8.4 Summary

This section answered a practical question that cannot be avoided: when diagonalization fails, how far can a matrix still be simplified?

The answer is given by the Jordan canonical form theorem: every complex square matrix is similar to a Jordan canonical form, assembled from Jordan blocks along the diagonal. A Jordan block $\mathbf{J}_k(\lambda)$ is a nearly diagonal matrix of size $k$ whose diagonal entries are all the same eigenvalue $\lambda$ , with a line of 1s on the superdiagonal immediately above—these 1s are precisely the algebraic fingerprint of a chain of generalized eigenvectors. From the point of view of the pirate-ship analogy, the $\lambda$ on the diagonal describes the “scaling” at each step, and the 1 above it describes the “passing along” of vectors on the chain; the length of the chain equals the size of the Jordan block, which is exactly the depth of the corresponding generalized eigenspace.

The Jordan canonical form is one of the ultimate classification theorems of finite-dimensional complex linear algebra: it tells us that, up to similarity, the inner structure of a matrix is completely determined by two sets of data—the set of eigenvalues, and the distribution of Jordan block sizes for each eigenvalue. Diagonalizable matrices are just the special case in which every Jordan block is $1 \times 1$ .

In actual computation, the Jordan canonical form is numerically sensitive and is usually replaced by the Schur decomposition; but as a theoretical tool it provides the clearest theoretical framework for matrix functions, for the general solutions of differential equations, and for the convergence analysis of the numerical algorithms in §8.5.

8.5 Numerical Methods for Computing Eigenvalues¶

The earlier sections of this chapter have already shown NumPy’s basic tools for computing eigenvalues (such as np.linalg.eig()). In this section we explore a few basic iterative numerical methods: the power method and the inverse power method. These methods are especially useful for large matrices or when only particular eigenvalues are needed.

8.5.1 The Power Method: Finding the Largest Eigenvalue Iteratively¶

When we only need the largest eigenvalue of a matrix (the eigenvalue of largest modulus) and its eigenvector, the power method provides a simple and effective iterative algorithm.

Basic principle¶

Suppose the matrix $\mathbf{A}$ has $n$ linearly independent eigenvectors $\mathbf{v}_0, \mathbf{v}_1, \ldots, \mathbf{v}_{n-1}$ , with corresponding eigenvalues satisfying

|\lambda_0| > |\lambda_1| \geq |\lambda_2| \geq \cdots \geq |\lambda_{n-1}|

(210)

For any initial vector $\mathbf{x}_0 = \sum_{i=0}^{n-1} c_i \mathbf{v}_i$ (assuming $c_0 \neq 0$ ), multiply repeatedly:

\mathbf{x}_{k+1} = \frac{\mathbf{A}\mathbf{x}_k}{\|\mathbf{A}\mathbf{x}_k\|}

(211)

First normalize the eigenvectors to unit length. The unimodular coefficient in the equation below is a sign $\pm1$ in the real case and a phase in the complex case; the conclusion is convergence in direction, not necessarily convergence of the vectors themselves. One can show that

\mathbf{x}_k \approx \frac{\mathbf{A}^k \mathbf{x}_0}{\|\mathbf{A}^k \mathbf{x}_0\|} \approx \frac{c_0 \lambda_0^k \mathbf{v}_0}{\|c_0 \lambda_0^k \mathbf{v}_0\|} = \frac{c_0\lambda_0^k}{|c_0\lambda_0^k|}\mathbf{v}_0

(212)

that is, the sequence $\{\mathbf{x}_k\}$ converges to the direction of the dominant eigenvector $\mathbf{v}_0$ .

Estimating the eigenvalue with the Rayleigh quotient¶

We estimate the eigenvalue with the Rayleigh quotient; for real symmetric or Hermitian matrices, the error of the eigenvalue estimate can be one order higher than the error of the vector, but this does not change the vector iteration here or its number of iterations:

\lambda \approx R(\mathbf{A}, \mathbf{x}) = \frac{\mathbf{x}^* \mathbf{A} \mathbf{x}}{\mathbf{x}^* \mathbf{x}}

(213)

For a unit vector $\mathbf{x}$ , this simplifies to $R(\mathbf{A}, \mathbf{x}) = \mathbf{x}^* \mathbf{A} \mathbf{x}$ .

Convergence analysis¶

The rate of convergence of the power method depends on the eigenvalue ratio $|\lambda_1/\lambda_0|$ :

If $|\lambda_1| \ll |\lambda_0|$ (the dominant eigenvalue is much larger than the others), convergence is fast.
If $|\lambda_1| \approx |\lambda_0|$ (there are eigenvalues close to it), convergence is slow.
The number of iterations is about $O(\log(1/\epsilon) / \log(|\lambda_0/\lambda_1|))$ , where $\epsilon$ is the error tolerance.

Python implementation¶

import numpy as np

# --- Visual hierarchy helpers ---
def print_header(title):
    print("=" * 60)
    print(f"  {title}")
    print("=" * 60)

def print_step(step, desc):
    print(f"\n▶ Step {step}: {desc}")
    print("-" * 40)

def compare_print(label, actual, expected_desc):
    print(f"\n[{label}]")
    print(f"  Computed value: {actual}")
    print(f"  Expected value: {expected_desc}")

np.set_printoptions(precision=6, suppress=True, linewidth=100)

def power_method(A, tol=1e-10, max_iter=1000):
    """Stop on the scaled eigen-residual; return the estimate, a unit vector, the iteration count, and whether the tolerance was met.
    A small residual confirms an approximate eigenpair; singling out the dominant / smallest-modulus one still requires a spectral gap and suitable initial conditions.
    """
    A = np.asarray(A)
    if A.ndim != 2 or A.shape[0] != A.shape[1] or A.shape[0] == 0:
        raise ValueError("A must be a nonempty square matrix")
    if tol <= 0 or max_iter < 1:
        raise ValueError("tol and max_iter must be positive")
    rng = np.random.default_rng(42)
    x = rng.random(A.shape[0])
    if np.iscomplexobj(A):
        x = x + 1j * rng.random(A.shape[0])
    x /= np.linalg.norm(x)
    scale = np.linalg.norm(A, 2)
    for k in range(1, max_iter + 1):
        y = A @ x
        ynorm = np.linalg.norm(y)
        if ynorm == 0:
            return np.vdot(x, A @ x) / np.vdot(x, x), x, k, True
        x = y / ynorm
        value = np.vdot(x, A @ x) / np.vdot(x, x)
        residual = np.linalg.norm(A @ x - value * x)
        if residual <= tol * (scale + abs(value)):
            return value, x, k, True
    return value, x, max_iter, False

# ============================================================
print_header("Example 8.5 | Finding the Dominant Eigenvalue by the Power Method")
# ============================================================

# --- Step 1: Define the matrix A ---
print_step("1", "Define the matrix A")
A = np.array([[4, -3,  0],
              [-3,  4,  0],
              [ 0,  0,  5]])
print(A)

# --- Step 2: Run the power method ---
print_step("2", "Run the power method (tol = 1e-10)")
lambda_max, v_max, iterations, converged = power_method(A)
print(f"  Dominant eigenvalue: λ_max    = {lambda_max:.10f}")
print(f"  Eigenvector:         v_max    = {v_max}")
print(f"  Iterations: {iterations}; tolerance met: {converged}")

# --- Step 3: Verify A v = λ v ---
print_step("3", "Verify the eigenvalue equation Av = λv")
print("  Av =")
print(f"  {A @ v_max}")
print("  λv =")
print(f"  {lambda_max * v_max}")
compare_print(
    "Residual ‖Av − λv‖",
    f"{np.linalg.norm(A @ v_max - lambda_max * v_max):.2e}",
    "for an exact eigenvector the residual should approach 0"
)

# --- Step 4: Compare with all the eigenvalues from NumPy ---
print_step("4", "Compare with all the eigenvalues from NumPy")
eigs_numpy = sorted(np.linalg.eigvals(A), key=abs, reverse=True)
compare_print(
    "Power method λ_max vs NumPy λ₀",
    f"{lambda_max:.10f}",
    f"NumPy λ₀ = {eigs_numpy[0]:.10f} (sorted by modulus)"
)
print(f"\n  All NumPy eigenvalues (in decreasing modulus): {[round(v, 6) for v in eigs_numpy]}")

8.5.2 The Inverse Power Method: Finding the Smallest Eigenvalue¶

The inverse power method is used to compute the smallest eigenvalue of a matrix (smallest in modulus) and its eigenvector. Its core idea is very clever:

Key observation¶

If $\lambda$ is an eigenvalue of $\mathbf{A}$ , then $1/\lambda$ is an eigenvalue of $\mathbf{A}^{-1}$ , with the same eigenvector:

\mathbf{A}\mathbf{v} = \lambda\mathbf{v} \quad \Rightarrow \quad \mathbf{A}^{-1}\mathbf{v} = \frac{1}{\lambda}\mathbf{v}

(214)

Therefore:

the smallest eigenvalue $\lambda_{\min}$ of $\mathbf{A}$ corresponds to the largest eigenvalue $1/\lambda_{\min}$ of $\mathbf{A}^{-1}$ ;
applying the power method to $\mathbf{A}^{-1}$ yields $\lambda_{\min}$ .

The algorithm¶

The iteration formula of the inverse power method is

\mathbf{A}\mathbf{x}_{k+1} = \mathbf{x}_k \quad \Rightarrow \quad \mathbf{x}_{k+1} = \mathbf{A}^{-1}\mathbf{x}_k

(215)

In actual computation, we do not compute $\mathbf{A}^{-1}$ explicitly; instead, in each iteration we solve the linear system

\mathbf{A}\mathbf{y} = \mathbf{x}_k, \quad \mathbf{x}_{k+1} = \frac{\mathbf{y}}{\|\mathbf{y}\|}

(216)

Key advantages¶

The key advantages of the inverse power method are:

No explicit inverse is needed. Each iteration only requires solving the linear system $\mathbf{A}\mathbf{y} = \mathbf{x}$ .
Numerical stability. For well-conditioned matrices, solving a linear system is more stable than inverting the matrix.
Rate of convergence. It depends on $|\lambda_{n-1}/\lambda_{n-2}|$ (with the eigenvalues arranged in decreasing order of modulus, the ratio of the eigenvalue of smallest modulus to the one of second-smallest modulus); this conclusion assumes that the matrix is invertible and diagonalizable, that the eigenvalue of smallest modulus is unique, and that the initial vector has a nonzero component along its eigendirection.

Python implementation¶

import numpy as np

# --- Visual hierarchy helpers ---
def print_header(title):
    print("=" * 60)
    print(f"  {title}")
    print("=" * 60)

def print_step(step, desc):
    print(f"\n▶ Step {step}: {desc}")
    print("-" * 40)

def compare_print(label, actual, expected_desc):
    print(f"\n[{label}]")
    print(f"  Computed value: {actual}")
    print(f"  Expected value: {expected_desc}")

np.set_printoptions(precision=6, suppress=True, linewidth=100)

# --- Reuse the power method of Example 8.5 ---
def power_method(A, tol=1e-10, max_iter=1000):
    """Stop on the scaled eigen-residual; return the estimate, a unit vector, the iteration count, and whether the tolerance was met.
    A small residual confirms an approximate eigenpair; singling out the dominant / smallest-modulus one still requires a spectral gap and suitable initial conditions.
    """
    A = np.asarray(A)
    if A.ndim != 2 or A.shape[0] != A.shape[1] or A.shape[0] == 0:
        raise ValueError("A must be a nonempty square matrix")
    if tol <= 0 or max_iter < 1:
        raise ValueError("tol and max_iter must be positive")
    rng = np.random.default_rng(42)
    x = rng.random(A.shape[0])
    if np.iscomplexobj(A):
        x = x + 1j * rng.random(A.shape[0])
    x /= np.linalg.norm(x)
    scale = np.linalg.norm(A, 2)
    for k in range(1, max_iter + 1):
        y = A @ x
        ynorm = np.linalg.norm(y)
        if ynorm == 0:
            return np.vdot(x, A @ x) / np.vdot(x, x), x, k, True
        x = y / ynorm
        value = np.vdot(x, A @ x) / np.vdot(x, x)
        residual = np.linalg.norm(A @ x - value * x)
        if residual <= tol * (scale + abs(value)):
            return value, x, k, True
    return value, x, max_iter, False

def inverse_power_method(A, tol=1e-10, max_iter=1000):
    """Stop on the scaled eigen-residual; return the estimate, a unit vector, the iteration count, and whether the tolerance was met.
    A small residual confirms an approximate eigenpair; singling out the dominant / smallest-modulus one still requires a spectral gap and suitable initial conditions.
    """
    A = np.asarray(A)
    if A.ndim != 2 or A.shape[0] != A.shape[1] or A.shape[0] == 0:
        raise ValueError("A must be a nonempty square matrix")
    if tol <= 0 or max_iter < 1:
        raise ValueError("tol and max_iter must be positive")
    rng = np.random.default_rng(42)
    x = rng.random(A.shape[0])
    if np.iscomplexobj(A):
        x = x + 1j * rng.random(A.shape[0])
    x /= np.linalg.norm(x)
    scale = np.linalg.norm(A, 2)
    for k in range(1, max_iter + 1):
        y = np.linalg.solve(A, x)
        ynorm = np.linalg.norm(y)
        if ynorm == 0:
            return np.vdot(x, A @ x) / np.vdot(x, x), x, k, True
        x = y / ynorm
        value = np.vdot(x, A @ x) / np.vdot(x, x)
        residual = np.linalg.norm(A @ x - value * x)
        if residual <= tol * (scale + abs(value)):
            return value, x, k, True
    return value, x, max_iter, False

# ============================================================
print_header("Example 8.6 | Finding the Smallest Eigenvalue by the Inverse Power Method")
# ============================================================

# --- Step 1: Define the matrix A ---
print_step("1", "Define the matrix A")
A = np.array([[4, -3,  0],
              [-3,  4,  0],
              [ 0,  0,  5]], dtype=float)
print(A)

# --- Step 2: Run the inverse power method ---
print_step("2", "Run the inverse power method (tol = 1e-10)")
lambda_min, v_min, iter_min, converged_min = inverse_power_method(A)
print(f"  Smallest eigenvalue: λ_min    = {lambda_min:.10f}")
print(f"  Eigenvector:         v_min    = {v_min}")
print(f"  Iterations: {iter_min}; tolerance met: {converged_min}")

# --- Step 3: Verify A v = λ v ---
print_step("3", "Verify the eigenvalue equation Av = λv")
print("  Av =")
print(f"  {A @ v_min}")
print("  λv =")
print(f"  {lambda_min * v_min}")
compare_print(
    "Residual ‖Av − λv‖",
    f"{np.linalg.norm(A @ v_min - lambda_min * v_min):.2e}",
    "for an exact eigenvector the residual should approach 0"
)

# --- Step 4: Compare with all the eigenvalues from NumPy ---
print_step("4", "Compare with all the eigenvalues from NumPy")
eigs_numpy = sorted(np.linalg.eigvals(A), key=abs)
compare_print(
    "Inverse power method λ_min vs NumPy λ₀ (smallest modulus)",
    f"{lambda_min:.10f}",
    f"NumPy λ₀ = {eigs_numpy[0]:.10f}"
)
print(f"\n  All NumPy eigenvalues (in increasing modulus): {[round(v, 6) for v in eigs_numpy]}")

# --- Step 5: Power method vs inverse power method ---
print_step("5", "Power method vs inverse power method")
lambda_max, _, iter_max, converged_max = power_method(A)
compare_print(
    "Dominant eigenvalue λ_max (power method)",
    f"{lambda_max:.6f}  ({iter_max} iterations, tolerance met: {converged_max})",
    "the eigenvalue of largest modulus"
)
compare_print(
    "Smallest eigenvalue λ_min (inverse power method)",
    f"{lambda_min:.6f}  ({iter_min} iterations, tolerance met: {converged_min})",
    "the eigenvalue of smallest modulus"
)
print(f"\n  Convergence ratio |λ_min / λ_next| = {abs(eigs_numpy[0] / eigs_numpy[1]):.6f}")
print(f"  (the smaller this ratio, the faster the inverse power method converges)")

§8.5 Summary

This section asked: how can eigenvalues, defined in theory, be computed efficiently on a computer?

Finding the roots of the characteristic polynomial is both unstable and expensive in high dimensions; in its place come iterative algorithms that act directly on the matrix. The core idea of the power method is remarkably simple: apply the matrix repeatedly and normalize, and the direction of the vector gradually converges under the “gravitational pull” of the dominant eigenvector—the rate of convergence is governed by the eigenvalue ratio $|\lambda_1/\lambda_0|$ . The inverse power method turns the problem around: power iteration on the inverse matrix locates the eigenvalue of smallest modulus, and each iteration only requires solving one linear system, with no need at all to compute the inverse explicitly. In both methods the Rayleigh quotient is used to estimate the eigenvalue; it does not change the vector iteration, and it is a different algorithm from Rayleigh quotient iteration, which updates the shift at every step.

The beauty of these algorithms is that they use only matrix multiplication (or the equivalent solution of linear equations), which makes them especially well suited to large sparse matrices—and the matrices in real applications (such as the link matrix of PageRank) are often exactly of this kind. Each multiplication in power iteration is, in essence, also a display of the Jordan structure: convergence is governed mainly by the spectral ratio and the initial vector, and nontrivial Jordan blocks may also introduce polynomial factors. This is the most direct bridge between the theory and numerical computation.

8.6 Summary: From Invariant Directions to a Complete Picture of the Matrix¶

This chapter started from a question that can be stated in almost a single sentence: does a linear transformation have vectors whose “direction stays fixed”? Starting from this question, we followed a tightly linked logical path through five sections and arrived at one of the deepest structure theorems of finite-dimensional linear algebra. Before closing the chapter, it is worth looking back along this path and confirming its place in the book as a whole.

Review of the Theoretical Thread¶

§8.1 set up the language: the characteristic equation $\det(\mathbf{A} - \lambda\mathbf{I}) = 0$ translates the geometric question “which directions are invariant” into the problem of finding the roots of the characteristic polynomial. The fundamental theorem of algebra guarantees that the roots exist over the complex numbers, and geometric intuition tells us how eigenvalues encode stretching, compression, and even reversal along each direction—a two-dimensional rotation through an angle that is not an integer multiple of $\pi$ has no real eigenvector, which is precisely the algebraic statement that in this case there is no real invariant line.

§8.2 deepened our understanding of eigenvalues: the trace equals the sum of the eigenvalues, and the determinant equals their product; these two equations anchor global properties of the matrix directly to its eigenvalue structure. More important is the inequality $1 \leq g_\lambda \leq a_\lambda$ between geometric and algebraic multiplicity, which measures precisely how far a matrix is from being diagonalizable. The Cayley–Hamilton theorem builds a bridge between matrices and the algebra of polynomials: every matrix is a solution of its own characteristic equation, and this seemingly magical conclusion is the algebraic foundation of the Jordan theory that follows.

§8.3 answered the global question: if the characteristic polynomial splits completely over the field in use (for example, when working over $\mathbb{C}$ ) and $g_\lambda = a_\lambda$ for all eigenvalues, then the matrix can be brought to the diagonal form $\mathbf{D} = \mathbf{P}^{-1}\mathbf{A}\mathbf{P}$ using eigenvectors as the basis. The similarity invariants (trace, determinant, characteristic polynomial) ensure that a change of basis does not change the essence of the linear transformation. Once diagonalization is achieved, powers of the matrix, functions of the matrix, and even difference equations and systems of differential equations all revolve around scalar operations on the eigenvalues; Binet’s formula for the Fibonacci sequence is the most elegant demonstration of this idea.

§8.4 faced the case of failure honestly: not every matrix is diagonalizable. The Jordan canonical form theorem tells us that even when diagonalization fails, every complex square matrix is similar to a nearly diagonal matrix assembled from Jordan blocks—the sizes of the blocks record the depths of the chains of generalized eigenvectors, and diagonalizability is just the special case in which every Jordan block degenerates to $1 \times 1$ . This is one of the ultimate classification results of finite-dimensional linear algebra.

§8.5 brought the theory back down to earth: power iteration and the inverse power method use the most basic matrix multiplication (or the solution of linear equations) to approach the dominant eigenvalue and the smallest eigenvalue step by step; the Rayleigh quotient provides an estimate of the eigenvalue; convergence is governed mainly by the spectral ratio and the initial vector; and nontrivial Jordan blocks may also introduce polynomial factors. This path from theory to algorithm embodies the enduring spirit of linear algebra, in which abstraction and computation complement each other.

Connections to Other Chapters¶

This chapter does not stand alone. It starts from the determinants of Chapter 7: the characteristic equation $\det(\mathbf{A} - \lambda\mathbf{I}) = 0$ is itself an application of the determinant, and the identity “the determinant equals the product of the eigenvalues” is the most direct bridge between the two chapters. In Chapter 7 the determinant described the factor by which a real linear map scales signed volume, with the ordinary volume factor given by its absolute value; in this chapter that factor is decomposed into the product of the eigenvalues, and its geometric meaning becomes clearer as a result—the matrix scales space separately along each eigendirection, and the determinant is the total product of these scalings.

Chapter 9 introduces inner product structure, and there the theory of eigenvalues takes a qualitative leap. Symmetric (Hermitian) matrices not only have real eigenvalues; their eigenvectors can also be chosen to form a mutually orthogonal basis—diagonalization is then no longer merely an algebraic possibility but the most natural geometric decomposition. Once the spectral decomposition of §9.3 is available, one can see directly that the values of the Rayleigh quotient $R(\mathbf{A}, \mathbf{x}) = \mathbf{x}^* \mathbf{A}\mathbf{x} / \|\mathbf{x}\|^2$ lie exactly in $[\lambda_{\min}, \lambda_{\max}]$ —the simplest form of the min-max theorem. The similarity invariants discussed in this chapter are refined further, in the context of the symmetric matrices of Chapter 9, into the spectral theorem, which gives a structurally stronger guarantee.

Further on, the singular value decomposition (SVD) of Chapter 11 can be viewed as the natural generalization of eigenvalue theory to non-square matrices: for a general $m \times n$ matrix $\mathbf{A}$ , the matrices $\mathbf{A}^*\mathbf{A}$ and $\mathbf{A}\mathbf{A}^*$ have the same nonzero eigenvalues, equal to the squares of the nonzero singular values; in their full spectra, the numbers of zero eigenvalues are $n-\operatorname{rank}\mathbf{A}$ and $m-\operatorname{rank}\mathbf{A}$ respectively, and the left and right singular vectors are eigenvectors of these two positive semidefinite matrices. The whole language built in this chapter—eigenspaces, algebraic and geometric multiplicity, the minimal polynomial—will be borrowed there, but its object expands from the eigenvalue problem of square matrices to more general matrix factorizations.

The Role of This Chapter in the Book¶

If the abstract theory of linear spaces in Chapter 4 laid down the language of the whole book, and the determinants of Chapter 7 supplied the geometric tool of “volume,” then what this chapter does is use eigenvalues and eigenvectors to open up the internal structure of matrices.

Before this, a matrix was still a black box to us: vectors go in, vectors come out, but it is hard to say clearly what the geometric essence of the transformation itself is. Eigenvalue theory gives us a key—find the “principal axis directions” of the matrix and the corresponding scaling ratios, and the geometric behavior of the linear transformation decomposes from a complicated mixture of effects into independent scaling along each eigendirection. From diagonalization to the Jordan canonical form, this decomposition extends from the most ideal case all the way to the most general one, and finally gives a complete picture of the classification in finite-dimensional complex linear algebra.

The influence of this picture reaches far beyond linear algebra itself. The long-term behavior of dynamical systems is determined by the dominant eigenvalue; in quantum mechanics, the eigenvalues of the Hamiltonian operator are the energies of the system; Google’s PageRank algorithm essentially computes the dominant eigenvector of the transition matrix of a Markov chain; and in deep learning, the spectral structure of attention matrices determines a model’s ability to capture long-range dependencies. The eigenvalue problem has flourished for three hundred years precisely because “searching for invariant directions” is itself one of the most pervasive structures in nature and in the world of information.

In this sense, this chapter is not only the end point of one theory but also the starting point of the chapters that follow: with the perspective of eigenvalues, look again at the symmetric matrices of Chapter 9, the quantum operators of Chapter 10, and the singular value decomposition of Chapter 11, and you will find that those more complicated theories are simply different answers, in different settings, to the same core question.

Chapter Structure and Learning Objectives¶

8.1 Fundamentals of Eigenvalues and Eigenvectors¶

8.1.1 The Characteristic Equation and the Characteristic Polynomial¶

8.1.2 Geometric Interpretation of Eigenvalues and Eigenvectors¶

8.2 Theoretical Properties of Eigenvalues¶

8.2.1 Relations Between Eigenvalues and Matrix Properties¶

8.2.2 Eigenspaces, Algebraic Multiplicity, and Geometric Multiplicity¶

8.2.3 ◆The Cayley–Hamilton Theorem¶

8.3 Similarity Transformations and Diagonalization¶

8.3.1 Similarity Transformations: Looking at the Same Thing from Another Angle¶

8.3.2 Definition of and Conditions for Diagonalization¶

8.3.3 Applications of Diagonalization¶

8.4 The Jordan Canonical Form¶

8.4.1 Overview of the Jordan Canonical Form¶

The Concavity Theorem for Kernel Dimension Increments¶

8.4.2 Applications of the Jordan Canonical Form¶

Application 1: Computing powers of a matrix¶

Application 2: Solving difference equations¶

Application 3: The matrix exponential¶

Application 4: Solving systems of differential equations¶

Generalization to Jordan Blocks of Any Size¶

8.5 Numerical Methods for Computing Eigenvalues¶

8.5.1 The Power Method: Finding the Largest Eigenvalue Iteratively¶

Basic principle¶

Estimating the eigenvalue with the Rayleigh quotient¶

Convergence analysis¶

Python implementation¶

8.5.2 The Inverse Power Method: Finding the Smallest Eigenvalue¶

Key observation¶

The algorithm¶

Key advantages¶

Python implementation¶

8.6 Summary: From Invariant Directions to a Complete Picture of the Matrix¶

Review of the Theoretical Thread¶

Connections to Other Chapters¶

The Role of This Chapter in the Book¶

Concept Map¶