Dictionary of Applied Machine Learning · singular value decomposition

singular value decomposition — Python demo

Numerical companion to the entry singular value decomposition: 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 with python3 svd.py, from any directory — it writes its output files into the current directory. Requires NumPy and Matplotlib only, and uses fixed seeds, so the printed numbers reproduce exactly. Download svd.py

The script, block by block

One cell per block of the script: the code, and what that code printed when it last ran here

setup

"""
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
         returned by np.linalg.svd reconstruct A to machine precision,
         V and U are orthonormal (V^T V = I, U^T U = I), and — in
         contrast to an EVD — the factorization exists for every matrix,
         including the rectangular 5 x 3 example and the defective
         matrix [[0, 1], [0, 0]] that admits no EVD.
[P-sing] Lambda is nonzero only on its main diagonal and its diagonal
         entries, the singular values, are nonnegative.

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

rng = np.random.default_rng(42)
report = []


def check(name, ok):
    report.append((name, bool(ok)))
    print(f"  [{'ok' if ok else 'FAIL'}] {name}")

P-def

The SVD A = V Lambda U^T of a rectangular matrix: the factors returned by np.linalg.svd reconstruct A to machine precision, V and U are orthonormal (V^T V = I, U^T U = I), and — in contrast to an EVD — the factorization exists for every matrix, including the rectangular 5 x 3 example and the defective matrix [[0, 1], [0, 0]] that admits no EVD.

print("[P-def] A = V Lambda U^T with orthonormal V, U — any matrix")
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)))
# exists even for the defective matrix that has no EVD
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)

[P-def] A = V Lambda U^T with orthonormal V, U — any matrix
  [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] defective [[0,1],[0,0]] (no EVD) still has an exact SVD

P-sing

Lambda is nonzero only on its main diagonal and its diagonal entries, the singular values, are nonnegative.

print("[P-sing] Lambda diagonal, singular values nonnegative")
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", np.all(s >= 0))
check("np.linalg.svd returns them sorted (descending)",
      np.all(np.diff(s) <= 0))

# ------------------------------------------------------------ preview
fig, ax = plt.subplots(1, 3, figsize=(9, 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)
fig.suptitle("[P-def] SVD of a rectangular matrix")
fig.tight_layout()
fig.savefig("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-sing] Lambda diagonal, singular values nonnegative
  [ok] Lambda vanishes off the main diagonal
  [ok] singular values are nonnegative
  [ok] np.linalg.svd returns them sorted (descending)

7/7 checks passed