P-spec writes when the script runsDictionary of Applied Machine Learning · eigenvalue decomposition (EVD)
Numerical companion to the entry eigenvalue decomposition (EVD): it recomputes what the entry states and prints one line per check
One block per paragraph of the entry (marked [P...]): each block verifies numerically what the corresponding statement asserts. Self-contained (numpy/matplotlib only), fixed seed.
Run it without installing anything:uv run https://dictionaryofml.org/terms/evd.py
uv downloads this script and the pinned NumPy and Matplotlib it needs, then runs it; the script fetches any input file it uses. To keep the output files, download evd.py into a folder and run uv run evd.py there. With NumPy and Matplotlib already installed, python3 evd.py, from any directory — it writes its output files into the current directory. Fixed seeds, so the printed numbers reproduce exactly. Download evd.py · Notebook · Open in Colab
One cell per block of the script: the code, and what that code printed when it last ran here
"""
evd.py — numerical companion to the glossary entry
'eigenvalue decomposition (EVD)'.
One block per paragraph of the entry (marked [P...]): each block verifies
numerically what the corresponding statement asserts. Self-contained
(numpy/matplotlib only), fixed seed.
Blocks
------
[P-def] The EVD A = V Lambda V^{-1}: reconstructing A from the factors
returned by np.linalg.eig gives A back with error below 1e-12;
the columns of V are eigenvectors and Lambda is diagonal with
the matching eigenvalues. The entry's 2x2 example matrix
[[1.3, 0.8], [0.4, 0.9]] has eigenvalues 1.7 and 0.5 with
eigenvectors along (2, 1) and (-1, 1), as drawn in Fig. 1.
[P-diag] Matrices that admit an EVD are the diagonalizable ones: for the
defective matrix [[0, 1], [0, 0]] the eigenvector matrix is
singular (rank 1), so V^{-1} does not exist and no EVD is
possible, while the diagonalizable matrix A above passes.
[P-fast] The EVD speeds up computations: for GD on a linear regression
training set, transforming with the orthogonal eigenvector
matrix of X^T X decouples the update into element-wise
operations; the decoupled recursion reproduces the plain GD
iterates exactly.
[P-spec] Spectral clustering builds on the EVD of a graph Laplacian:
the Laplacian is symmetric and psd, so its EVD has an
orthogonal eigenvector matrix (V^{-1} = V^T) and real
nonnegative eigenvalues with lambda_1 = 0; the second-smallest
eigenvalue lambda_2 is positive iff the graph is connected,
and the signs of the entries of the corresponding eigenvector
(the Fiedler vector) split the two three-node clusters of the
entry's six-node example graph; the rounded entries match the
column vector displayed in the entry's Fig. 2.
Outputs
-------
evd.png : preview figure (checking only).
Data generated by pythondemos/evd.py.
"""
import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib.patches import Patch
from pathlib import Path
OUT_DIR = Path(__file__).parent
rng = np.random.default_rng(42)
report = []
def check(name, ok):
report.append((name, bool(ok)))
print(f" [{'ok' if ok else 'FAIL'}] {name}")
The EVD A = V Lambda V^{-1}: reconstructing A from the factors returned by np.linalg.eig gives A back with error below 1e-12; the columns of V are eigenvectors and Lambda is diagonal with the matching eigenvalues. The entry's 2x2 example matrix [[1.3, 0.8], [0.4, 0.9]] has eigenvalues 1.7 and 0.5 with eigenvectors along (2, 1) and (-1, 1), as drawn in Fig. 1.
print("[P-def] A = V Lambda V^{-1} reconstructs A")
A = np.array([[2.0, 1.0, 0.0], [1.0, 3.0, 0.5], [0.0, 0.5, 1.5]])
lam, V = np.linalg.eig(A)
Lam = np.diag(lam)
A_rec = V @ Lam @ np.linalg.inv(V)
check("reconstruction error < 1e-12",
np.max(np.abs(A_rec - A)) < 1e-12)
check("columns of V are eigenvectors (A v = lambda v)",
all(np.linalg.norm(A @ V[:, i] - lam[i] * V[:, i]) < 1e-12
for i in range(3)))
check("Lambda is diagonal", np.allclose(Lam, np.diag(np.diag(Lam))))
A2 = np.array([[1.3, 0.8], [0.4, 0.9]]) # the entry's 2x2 example
lam2 = np.sort(np.linalg.eigvals(A2))[::-1]
check("the 2x2 example matrix has eigenvalues 1.7 and 0.5",
np.allclose(lam2, [1.7, 0.5]))
check("its eigenvectors lie along (2, 1) and (-1, 1)",
np.linalg.norm(A2 @ [2, 1] - 1.7 * np.array([2, 1])) < 1e-12
and np.linalg.norm(A2 @ [-1, 1] - 0.5 * np.array([-1, 1])) < 1e-12)
[P-def] A = V Lambda V^{-1} reconstructs A
[ok] reconstruction error < 1e-12
[ok] columns of V are eigenvectors (A v = lambda v)
[ok] Lambda is diagonal
[ok] the 2x2 example matrix has eigenvalues 1.7 and 0.5
[ok] its eigenvectors lie along (2, 1) and (-1, 1)
Matrices that admit an EVD are the diagonalizable ones: for the defective matrix [[0, 1], [0, 0]] the eigenvector matrix is singular (rank 1), so V^{-1} does not exist and no EVD is possible, while the diagonalizable matrix A above passes.
print("[P-diag] EVD exists iff the matrix is diagonalizable")
D = np.array([[0.0, 1.0], [0.0, 0.0]]) # defective (not diagonalizable)
lamD, VD = np.linalg.eig(D)
rank_VD = np.linalg.matrix_rank(VD)
check("defective matrix: eigenvector matrix is singular (rank 1)",
rank_VD == 1)
sym = A # symmetric example above
check("diagonalizable matrix: eigenvector matrix invertible (rank 3)",
np.linalg.matrix_rank(V) == 3)
[P-diag] EVD exists iff the matrix is diagonalizable [ok] defective matrix: eigenvector matrix is singular (rank 1) [ok] diagonalizable matrix: eigenvector matrix invertible (rank 3)
The EVD speeds up computations: for GD on a linear regression training set, transforming with the orthogonal eigenvector matrix of X^T X decouples the update into element-wise operations; the decoupled recursion reproduces the plain GD iterates exactly.
print("[P-fast] EVD of X^T X decouples the GD update for linear regression")
m_tr, d = 50, 4
X = rng.normal(size=(m_tr, d))
yv = rng.normal(size=m_tr)
lamX, VX = np.linalg.eigh(X.T @ X)
eta = 0.02
w = np.zeros(d) # plain GD
wt = VX.T @ w # decoupled coordinates
b = VX.T @ (X.T @ yv)
for _ in range(30):
w = w - (2 * eta / m_tr) * (X.T @ (X @ w - yv))
wt = (1 - (2 * eta / m_tr) * lamX) * wt + (2 * eta / m_tr) * b
check("eigenvector matrix of the symmetric X^T X is orthogonal",
np.max(np.abs(VX.T @ VX - np.eye(d))) < 1e-12)
check("the element-wise recursion reproduces the GD iterates exactly",
np.max(np.abs(VX @ wt - w)) < 1e-10)
[P-fast] EVD of X^T X decouples the GD update for linear regression [ok] eigenvector matrix of the symmetric X^T X is orthogonal [ok] the element-wise recursion reproduces the GD iterates exactly
Spectral clustering builds on the EVD of a graph Laplacian: the Laplacian is symmetric and psd, so its EVD has an orthogonal eigenvector matrix (V^{-1} = V^T) and real nonnegative eigenvalues with lambda_1 = 0; the second-smallest eigenvalue lambda_2 is positive iff the graph is connected, and the signs of the entries of the corresponding eigenvector (the Fiedler vector) split the two three-node clusters of the entry's six-node example graph; the rounded entries match the column vector displayed in the entry's Fig. 2.
print("[P-spec] EVD of a graph Laplacian: spectral clustering")
n_nodes = 6
# two clusters of three nodes each, joined by a single edge (2, 3)
edge_list = [(0, 1), (1, 2), (0, 2), (3, 4), (4, 5), (3, 5), (2, 3)]
L = np.zeros((n_nodes, n_nodes))
for i, j in edge_list:
L[i, j] -= 1.0
L[j, i] -= 1.0
L[i, i] += 1.0
L[j, j] += 1.0
lamL, VL = np.linalg.eigh(L)
check("Laplacian symmetric: eigenvector matrix orthogonal (V^T V = I)",
np.max(np.abs(VL.T @ VL - np.eye(n_nodes))) < 1e-12)
check("eigenvalues nonnegative with smallest lambda_1 = 0",
lamL[0] > -1e-12 and abs(lamL[0]) < 1e-12)
check("graph connected: second-smallest eigenvalue lambda_2 > 0",
lamL[1] > 1e-8)
fiedler = VL[:, 1]
if fiedler[0] > 0: # the sign of an eigenvector
fiedler = -fiedler # is arbitrary; fix node 1 < 0
side = fiedler > 0
check("signs of the Fiedler vector recover the two clusters",
len(set(side[:3])) == 1 and len(set(side[3:])) == 1
and side[0] != side[3])
check("the entries match the column vector displayed in the entry's "
"Fig. 2 (rounded to two decimals)",
np.allclose(np.round(fiedler, 2),
[-0.46, -0.46, -0.26, 0.26, 0.46, 0.46]))
# ------------------------------------------------------------ preview
fig, ax = plt.subplots(1, 4, figsize=(12, 2.8))
for a, M, t in ((ax[0], A, "A"), (ax[1], Lam.real, "Lambda"),
(ax[2], (A_rec - A).real, "V Lambda V^{-1} - A")):
im = a.imshow(M, cmap="gray"); a.set_title(t)
fig.colorbar(im, ax=a, shrink=0.75)
nodes = np.arange(1, n_nodes + 1)
bars = ax[3].bar(nodes, fiedler, color=np.where(side, "C0", "C1"))
for b, s in zip(bars, side):
if s:
b.set_hatch("//")
legend_handles = [Patch(facecolor="C0", hatch="//", label="cluster 1"),
Patch(facecolor="C1", label="cluster 2")]
ax[3].set_xlim(0.4, n_nodes + 0.6)
ax[3].axhline(0.0, color="k", lw=0.8)
ax[3].set_xlabel("node $i$")
ax[3].set_ylabel("$v^{(2)}_i$")
ax[3].set_title("[P-spec] Fiedler vector")
ax[3].legend(handles=legend_handles, frameon=False)
fig.suptitle("EVD factors, reconstruction error, and spectral clustering")
fig.tight_layout()
fig.savefig(OUT_DIR / "evd.png", dpi=110)
print(f"\n{sum(ok for _, ok in report)}/{len(report)} checks passed")
assert all(ok for _, ok in report)
[P-spec] EVD of a graph Laplacian: spectral clustering [ok] Laplacian symmetric: eigenvector matrix orthogonal (V^T V = I) [ok] eigenvalues nonnegative with smallest lambda_1 = 0 [ok] graph connected: second-smallest eigenvalue lambda_2 > 0 [ok] signs of the Fiedler vector recover the two clusters [ok] the entries match the column vector displayed in the entry's Fig. 2 (rounded to two decimals) 14/14 checks passed
P-spec writes when the script runs