P-pinv writes when the script runsDictionary of Applied Machine Learning · singular value decomposition (SVD)
Numerical companion to the entry singular value decomposition (SVD): 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/svd.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 svd.py into a folder and run uv run svd.py there. With NumPy and Matplotlib already installed, python3 svd.py, from any directory — it writes its output files into the current directory. Fixed seeds, so the printed numbers reproduce exactly. Download svd.py · Notebook · Open in Colab
One cell per block of the script: the code, and what that code printed when it last ran here
"""
svd.py — numerical companion to the glossary entry
'singular value decomposition (SVD)'.
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 SVD A = V Lambda U^T of a rectangular matrix: the factors
delivered by np.linalg.svd reconstruct A with error below
1e-12, V and U are orthonormal (V^T V = I, U^T U = I),
Lambda is nonzero only on its main diagonal with
nonnegative entries sorted in descending order, and
A u^(j) = lambda_j v^(j) for every j.
[P-exist] An SVD exists for every matrix, including the defective
matrix [[0, 1], [0, 0]] that admits no EVD; for a symmetric
psd matrix the SVD coincides with the EVD; the largest
singular value equals the spectral norm and the ratio of
largest to smallest nonzero singular value equals the
condition number.
[P-lowrank] Eckart-Young: truncating after the k largest singular
values gives error ||A - A_k||_2 = lambda_{k+1}, and no
random rank-k matrix among 300 samples does better; keeping
k singular values of an m x d image matrix stores
k(m + d + 1) numbers.
[P-pinv] The pseudoinverse from the SVD (invert the nonzero singular
values) matches np.linalg.pinv, and X^+ y solves the
least-squares problem for linear regression.
Outputs
-------
svd.png : preview figure (checking only).
Data generated by pythondemos/svd.py.
"""
import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
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 SVD A = V Lambda U^T of a rectangular matrix: the factors delivered by np.linalg.svd reconstruct A with error below 1e-12, V and U are orthonormal (V^T V = I, U^T U = I), Lambda is nonzero only on its main diagonal with nonnegative entries sorted in descending order, and A u^(j) = lambda_j v^(j) for every j.
print("[P-def] A = V Lambda U^T with orthonormal V, U")
A = rng.normal(size=(5, 3)) # rectangular
Vf, s, Ut = np.linalg.svd(A) # full SVD
Lam = np.zeros((5, 3)); Lam[:3, :3] = np.diag(s)
check("reconstruction V Lambda U^T = A (err < 1e-12)",
np.max(np.abs(Vf @ Lam @ Ut - A)) < 1e-12)
check("V orthonormal: V^T V = I", np.allclose(Vf.T @ Vf, np.eye(5)))
check("U orthonormal: U^T U = I", np.allclose(Ut @ Ut.T, np.eye(3)))
off_diag = Lam.copy(); np.fill_diagonal(off_diag, 0.0)
check("Lambda vanishes off the main diagonal", np.all(off_diag == 0))
check("singular values are nonnegative and sorted (descending)",
np.all(s >= 0) and np.all(np.diff(s) <= 0))
check("A u^(j) = lambda_j v^(j) for every j",
all(np.linalg.norm(A @ Ut[j] - s[j] * Vf[:, j]) < 1e-12
for j in range(3)))
[P-def] A = V Lambda U^T with orthonormal V, U [ok] reconstruction V Lambda U^T = A (err < 1e-12) [ok] V orthonormal: V^T V = I [ok] U orthonormal: U^T U = I [ok] Lambda vanishes off the main diagonal [ok] singular values are nonnegative and sorted (descending) [ok] A u^(j) = lambda_j v^(j) for every j
An SVD exists for every matrix, including the defective matrix [[0, 1], [0, 0]] that admits no EVD; for a symmetric psd matrix the SVD coincides with the EVD; the largest singular value equals the spectral norm and the ratio of largest to smallest nonzero singular value equals the condition number.
print("[P-exist] an SVD exists for every matrix; special cases")
D = np.array([[0.0, 1.0], [0.0, 0.0]])
Vd, sd, Utd = np.linalg.svd(D)
check("defective [[0,1],[0,0]] (no EVD) still has an exact SVD",
np.max(np.abs(Vd @ np.diag(sd) @ Utd - D)) < 1e-14)
S = A.T @ A # symmetric psd 3 x 3
lam_evd = np.sort(np.linalg.eigvalsh(S))[::-1]
sing_S = np.linalg.svd(S, compute_uv=False)
check("symmetric psd: singular values equal the eigenvalues "
"(SVD = EVD)", np.allclose(sing_S, lam_evd))
check("largest singular value = spectral norm ||A||_2",
np.isclose(s[0], np.linalg.norm(A, 2)))
check("lambda_1 / lambda_min = condition number of A",
np.isclose(s[0] / s[-1], np.linalg.cond(A)))
[P-exist] an SVD exists for every matrix; special cases [ok] defective [[0,1],[0,0]] (no EVD) still has an exact SVD [ok] symmetric psd: singular values equal the eigenvalues (SVD = EVD) [ok] largest singular value = spectral norm ||A||_2 [ok] lambda_1 / lambda_min = condition number of A
Eckart-Young: truncating after the k largest singular values gives error ||A - A_k||_2 = lambda_{k+1}, and no random rank-k matrix among 300 samples does better; keeping k singular values of an m x d image matrix stores k(m + d + 1) numbers.
print("[P-lowrank] Eckart-Young: the truncated SVD is the best "
"low-rank approximation")
k = 1
A_k = s[0] * np.outer(Vf[:, 0], Ut[0])
err_trunc = np.linalg.norm(A - A_k, 2)
check("truncation error ||A - A_1||_2 equals lambda_2",
np.isclose(err_trunc, s[1]))
beaten = 0
for _ in range(300):
a, b = rng.normal(size=5), rng.normal(size=3)
B = np.outer(a, b)
B *= np.trace(B.T @ A) / np.trace(B.T @ B) # best scale for this B
if np.linalg.norm(A - B, 2) < err_trunc - 1e-9:
beaten += 1
check("no random rank-1 matrix among 300 samples beats the truncation",
beaten == 0)
m_px, d_px = 5, 3
check("storing the rank-k truncation takes k(m + d + 1) numbers",
k * (m_px + d_px + 1) == k * m_px + k * d_px + k)
[P-lowrank] Eckart-Young: the truncated SVD is the best low-rank approximation [ok] truncation error ||A - A_1||_2 equals lambda_2 [ok] no random rank-1 matrix among 300 samples beats the truncation [ok] storing the rank-k truncation takes k(m + d + 1) numbers
The pseudoinverse from the SVD (invert the nonzero singular values) matches np.linalg.pinv, and X^+ y solves the least-squares problem for linear regression.
print("[P-pinv] the pseudoinverse from the SVD solves least squares")
Lam_pinv = np.zeros((3, 5)); Lam_pinv[:3, :3] = np.diag(1.0 / s)
A_pinv = Ut.T @ Lam_pinv @ Vf.T
check("U Lambda^+ V^T matches np.linalg.pinv",
np.allclose(A_pinv, np.linalg.pinv(A)))
yv = rng.normal(size=5)
w_hat = A_pinv @ yv
w_lstsq = np.linalg.lstsq(A, yv, rcond=None)[0]
check("A^+ y equals the least-squares solution",
np.allclose(w_hat, w_lstsq))
# ------------------------------------------------------------ preview
fig, ax = plt.subplots(1, 4, figsize=(12, 2.8))
for a, M, t in ((ax[0], A, "A (5 x 3)"), (ax[1], Lam, "Lambda"),
(ax[2], Vf @ Lam @ Ut - A, "reconstruction error")):
im = a.imshow(M, cmap="gray"); a.set_title(t)
fig.colorbar(im, ax=a, shrink=0.75)
ks = np.arange(0, 3)
errs = [np.linalg.norm(A - sum(s[j] * np.outer(Vf[:, j], Ut[j])
for j in range(kk)), 2) if kk else
np.linalg.norm(A, 2) for kk in ks]
ax[3].plot(ks, errs, "o-", c="k", label="$\\|A - A_k\\|_2$")
ax[3].plot(ks[:-1] + 1, s[1:], "s", mfc="white", mec="k",
label="$\\lambda_{k+1}$")
ax[3].set_xlabel("rank $k$"); ax[3].set_ylabel("error")
ax[3].set_xticks(ks)
ax[3].set_title("[P-lowrank] truncation error")
ax[3].legend(frameon=False)
fig.tight_layout()
fig.savefig(OUT_DIR / "svd.png", dpi=110)
print(f"\n{sum(ok for _, ok in report)}/{len(report)} checks passed")
assert all(ok for _, ok in report)
[P-pinv] the pseudoinverse from the SVD solves least squares [ok] U Lambda^+ V^T matches np.linalg.pinv [ok] A^+ y equals the least-squares solution 15/15 checks passed
P-pinv writes when the script runs