Part III: Decompositions | Chapter 1

Jordan Normal Form

The canonical form for matrices that cannot be diagonalized

Historical Context

Camille Jordan introduced the canonical form bearing his name in his 1870 Traité des substitutions et des équations algébriques, though Karl Weierstrass had independently developed equivalent results for matrix pencils. The Jordan form provides the finest invariant structure theory for linear operators: it classifies all matrices up to similarity, answering the fundamental question of when two matrices represent the same linear transformation in different bases.

While the Jordan form is theoretically crucial, it is numerically unstable—small perturbations can change the Jordan structure entirely. This led to the development of the Schur form and other numerically stable alternatives. Nevertheless, the Jordan form remains indispensable for theoretical analysis of matrix functions, differential equations, and control theory.

3.1 Generalized Eigenvectors

Definition: Generalized Eigenvector

A vector $v \neq 0$ is a generalized eigenvector of $A$ corresponding to $\lambda$ of rank $m$ if $(A - \lambda I)^m v = 0$ but $(A - \lambda I)^{m-1} v \neq 0$.

A Jordan chain of length $m$ starting from $v_m$ is the sequence$v_1, v_2, \ldots, v_m$ where $v_k = (A - \lambda I)^{m-k}v_m$, satisfying$(A - \lambda I)v_1 = 0$ and $(A - \lambda I)v_{k+1} = v_k$ for each $k$.

The generalized eigenspace is $G_\lambda = \ker(A - \lambda I)^n$. The dimension of$G_\lambda$ equals the algebraic multiplicity of $\lambda$, and the entire space decomposes as $V = G_{\lambda_1} \oplus \cdots \oplus G_{\lambda_k}$.

3.2 Jordan Blocks and Jordan Matrix

Definition: Jordan Block

A Jordan block $J_m(\lambda)$ of size $m$ is the $m \times m$ matrix with$\lambda$ on the diagonal and $1$s on the superdiagonal:

$$J_m(\lambda) = \begin{pmatrix} \lambda & 1 & & \\ & \lambda & \ddots & \\ & & \ddots & 1 \\ & & & \lambda \end{pmatrix}$$

Theorem (Jordan Normal Form)

Every matrix $A \in \mathbb{C}^{n \times n}$ is similar to a block-diagonal matrix$J = P^{-1}AP = \text{diag}(J_{m_1}(\lambda_1), \ldots, J_{m_r}(\lambda_r))$, where the Jordan blocks are unique up to permutation.

3.3 Computing the Jordan Form

The Jordan structure is determined by the ranks of successive powers of $(A - \lambda I)$:

Number of Jordan blocks for $\lambda$ = $\dim\ker(A - \lambda I)$ (geometric multiplicity)
Number of blocks of size $\geq k$ = $\dim\ker(A - \lambda I)^k - \dim\ker(A - \lambda I)^{k-1}$
Largest block size = smallest $m$ such that $\ker(A - \lambda I)^m = G_\lambda$

The minimal polynomial $m_A(\lambda)$ is the monic polynomial of smallest degree annihilating $A$. Its factorization $m_A(\lambda) = \prod(\lambda - \lambda_i)^{s_i}$ gives the sizes of the largest Jordan blocks: $s_i$ is the size of the largest block for $\lambda_i$.

3.4 Matrix Functions via Jordan Form

For any analytic function $f$, the value on a Jordan block is:

$$f(J_m(\lambda)) = \begin{pmatrix} f(\lambda) & f'(\lambda) & \frac{f''(\lambda)}{2!} & \cdots \\ & f(\lambda) & f'(\lambda) & \cdots \\ & & \ddots & \ddots \\ & & & f(\lambda) \end{pmatrix}$$

For the matrix exponential $e^{tJ_m(\lambda)}$, the $(i,j)$ entry (with $j \geq i$) is:

$$[e^{tJ_m(\lambda)}]_{ij} = e^{\lambda t} \frac{t^{j-i}}{(j-i)!}$$

This shows that defective eigenvalues produce polynomial-times-exponential behavior in ODE solutions $x(t) = e^{tA}x_0$, explaining the secular growth phenomenon in resonant systems.

3.5 Primary Decomposition Theorem

Theorem (Primary Decomposition)

If the characteristic polynomial factors as $\chi_A(\lambda) = \prod_{i=1}^k (\lambda - \lambda_i)^{n_i}$, then $V = G_{\lambda_1} \oplus \cdots \oplus G_{\lambda_k}$ where $G_{\lambda_i} = \ker(A - \lambda_i I)^{n_i}$and each $G_{\lambda_i}$ is $A$-invariant with $\dim G_{\lambda_i} = n_i$.

The primary decomposition is the conceptual foundation of the Jordan form. Each generalized eigenspace $G_{\lambda_i}$ is further decomposed into cyclic subspaces corresponding to individual Jordan chains, yielding the block diagonal structure. This decomposition is related to the Chinese Remainder Theorem in algebra and extends to modules over principal ideal domains.

Computational Laboratory

This simulation constructs Jordan blocks, computes their powers, demonstrates generalized eigenvectors, evaluates matrix exponentials via the Jordan form, and solves ODE systems with defective coefficient matrices.

Jordan Normal Form: Structure and Applications

Python

jordan_form.py153 lines

import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from scipy.linalg import expm

# ============================================================
# Jordan Normal Form: Generalized Eigenvectors and Applications
# ============================================================

# --- 1. Jordan block structure ---
print("=" * 60)
print("JORDAN BLOCKS")
print("=" * 60)

def jordan_block(lam, size):
    J = lam * np.eye(size)
    for i in range(size - 1):
        J[i, i+1] = 1
    return J

J1 = jordan_block(2, 3)
J2 = jordan_block(5, 2)
print(f"J(2, 3):\n{J1}")
print(f"\nJ(5, 2):\n{J2}")

# Full Jordan matrix
J = np.block([[J1, np.zeros((3, 2))],
              [np.zeros((2, 3)), J2]])
print(f"\nFull Jordan matrix J:\n{J}")

# --- 2. Powers of Jordan blocks ---
print("\n" + "=" * 60)
print("POWERS OF JORDAN BLOCKS")
print("=" * 60)

J_small = jordan_block(2, 3)
for k in [1, 2, 3, 5]:
    J_k = np.linalg.matrix_power(J_small.astype(int), k)
    print(f"J(2,3)^{k}:")
    print(f"{J_k}")
    print()

# --- 3. Generalized eigenvectors ---
print("=" * 60)
print("GENERALIZED EIGENVECTORS")
print("=" * 60)

A = np.array([[5, 4, 2, 1],
              [0, 1, -1, -1],
              [-1, -1, 3, 0],
              [1, 1, -1, 2]])

eigenvalues = np.linalg.eigvals(A)
print(f"Matrix A:\n{A}")
print(f"Eigenvalues: {np.round(eigenvalues, 4)}")

for lam in set(np.round(eigenvalues, 6)):
    lam = lam.real
    rank1 = np.linalg.matrix_rank(A - lam * np.eye(4))
    rank2 = np.linalg.matrix_rank(np.linalg.matrix_power((A - lam * np.eye(4)).astype(float), 2))
    gm = 4 - rank1
    print(f"\nλ = {lam:.4f}: rank(A-λI) = {rank1}, nullity = {gm}")
    print(f"  rank((A-λI)²) = {rank2}, nullity = {4 - rank2}")

# --- 4. Matrix exponential via Jordan form ---
print("\n" + "=" * 60)
print("MATRIX EXPONENTIAL VIA JORDAN FORM")
print("=" * 60)

J3 = jordan_block(1.0, 3)
print(f"J(1, 3):\n{J3}")

# exp(tJ) for Jordan block
def exp_jordan_block(lam, size, t):
    result = np.zeros((size, size))
    for i in range(size):
        for j in range(i, size):
            result[i, j] = np.exp(lam * t) * t**(j-i) / np.math.factorial(j-i)
    return result

for t in [0, 0.5, 1.0, 2.0]:
    exp_J = exp_jordan_block(1.0, 3, t)
    exp_J_scipy = expm(t * J3)
    print(f"t={t}: exp(tJ) match scipy? {np.allclose(exp_J, exp_J_scipy)}")
    if t == 1.0:
        print(f"  exp(J) =\n{np.round(exp_J, 4)}")

# --- 5. ODE system solution ---
print("\n" + "=" * 60)
print("ODE SYSTEM: dx/dt = Ax")
print("=" * 60)

A_ode = np.array([[-1, 1], [0, -1]])  # Jordan block
x0 = np.array([1.0, 1.0])
print(f"A = {A_ode.tolist()}, x0 = {x0}")
print(f"Solution: x(t) = exp(tA)x0")

for t in [0, 0.5, 1.0, 2.0, 5.0]:
    x_t = expm(t * A_ode) @ x0
    print(f"  t={t:.1f}: x = [{x_t[0]:.6f}, {x_t[1]:.6f}]")

# --- 6. Visualization ---
fig, axes = plt.subplots(1, 3, figsize=(16, 5))

# Plot 1: Powers of Jordan block entries
ax = axes[0]
ks = np.arange(1, 15)
J_base = jordan_block(0.9, 3)
entries = {"(1,1)": [], "(1,2)": [], "(1,3)": []}
for k in ks:
    Jk = np.linalg.matrix_power(J_base, k)
    entries["(1,1)"].append(abs(Jk[0,0]))
    entries["(1,2)"].append(abs(Jk[0,1]))
    entries["(1,3)"].append(abs(Jk[0,2]))
for (label, vals), color in zip(entries.items(), ["#3b82f6", "#0ea5e9", "#06b6d4"]):
    ax.semilogy(ks, vals, "o-", color=color, markersize=4, label=f"|J^k{label}|")
ax.set_xlabel("k", color="gray")
ax.set_title("Jordan Block Powers (λ=0.9)", fontsize=12, fontweight="bold")
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3)
ax.set_facecolor("#0f172a")

# Plot 2: ODE trajectories
ax = axes[1]
ts = np.linspace(0, 5, 200)
for x0_trial in [[1,0], [0,1], [1,1], [-1,1]]:
    x0_arr = np.array(x0_trial, dtype=float)
    traj = np.array([expm(t * A_ode) @ x0_arr for t in ts])
    ax.plot(traj[:, 0], traj[:, 1], linewidth=1.5)
ax.set_xlabel("x₁", color="gray")
ax.set_ylabel("x₂", color="gray")
ax.set_title("ODE Phase Portrait", fontsize=12, fontweight="bold")
ax.grid(True, alpha=0.3)
ax.set_facecolor("#0f172a")

# Plot 3: exp(tJ) components
ax = axes[2]
ts_plot = np.linspace(0, 3, 100)
for i, j, color, label in [(0,0,"#3b82f6","(1,1)"), (0,1,"#0ea5e9","(1,2)"), (0,2,"#06b6d4","(1,3)")]:
    vals = [exp_jordan_block(1.0, 3, t)[i,j] for t in ts_plot]
    ax.plot(ts_plot, vals, color=color, linewidth=2, label=f"exp(tJ){label}")
ax.set_xlabel("t", color="gray")
ax.set_title("exp(tJ) Components (λ=1)", fontsize=12, fontweight="bold")
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3)
ax.set_facecolor("#0f172a")

plt.tight_layout()
plt.savefig("output.png", dpi=150, bbox_inches="tight", facecolor="#0f172a")
plt.show()
print("\nVisualization saved.")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

← Spectral Theorem Singular Value Decomposition →