{
 "nbformat": 4,
 "nbformat_minor": 5,
 "metadata": {
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3",
   "language": "python"
  },
  "language_info": {
   "name": "python"
  },
  "colab": {
   "name": "svd.ipynb"
  }
 },
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "# singular value decomposition (SVD) \u2014 Python demo\n\nNumerical companion to the entry [singular value decomposition (SVD)](https://dictionaryofml.org/terms/svd.html) of the [Dictionary of Applied Machine Learning](https://dictionaryofml.org/): it recomputes what the entry states and prints one line per check.\n\nOne block per paragraph of the entry (marked [P...]): each block verifies numerically what the corresponding statement asserts. Self-contained (numpy/matplotlib only), fixed seed.\n\nRequires NumPy and Matplotlib only, and uses fixed seeds, so the printed numbers reproduce exactly. Generated from [`pythondemos/svd.py`](https://dictionaryofml.org/terms/svd.py); CC BY 4.0."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "# Notebook shim: the script resolves output paths relative to __file__,\n# which a notebook kernel does not define; everything lands in the\n# working directory instead.\nimport os\n__file__ = os.path.join(os.getcwd(), \"svd.py\")\nos.makedirs(\"pythondemos\", exist_ok=True)"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "\"\"\"\nsvd.py \u2014 numerical companion to the glossary entry\n'singular value decomposition (SVD)'.\n\nOne block per paragraph of the entry (marked [P...]): each block verifies\nnumerically what the corresponding statement asserts. Self-contained\n(numpy/matplotlib only), fixed seed.\n\nBlocks\n------\n[P-def]     The SVD A = V Lambda U^T of a rectangular matrix: the factors\n            delivered by np.linalg.svd reconstruct A with error below\n            1e-12, V and U are orthonormal (V^T V = I, U^T U = I),\n            Lambda is nonzero only on its main diagonal with\n            nonnegative entries sorted in descending order, and\n            A u^(j) = lambda_j v^(j) for every j.\n[P-exist]   An SVD exists for every matrix, including the defective\n            matrix [[0, 1], [0, 0]] that admits no EVD; for a symmetric\n            psd matrix the SVD coincides with the EVD; the largest\n            singular value equals the spectral norm and the ratio of\n            largest to smallest nonzero singular value equals the\n            condition number.\n[P-lowrank] Eckart-Young: truncating after the k largest singular\n            values gives error ||A - A_k||_2 = lambda_{k+1}, and no\n            random rank-k matrix among 300 samples does better; keeping\n            k singular values of an m x d image matrix stores\n            k(m + d + 1) numbers.\n[P-pinv]    The pseudoinverse from the SVD (invert the nonzero singular\n            values) matches np.linalg.pinv, and X^+ y solves the\n            least-squares problem for linear regression.\n\nOutputs\n-------\nsvd.png : preview figure (checking only).\n\nData generated by pythondemos/svd.py.\n\"\"\"\n\nimport numpy as np\nimport matplotlib\n\nmatplotlib.use(\"Agg\")\nimport matplotlib.pyplot as plt\nfrom pathlib import Path\n\nOUT_DIR = Path(__file__).parent\n\nrng = np.random.default_rng(42)\nreport = []\n\n\ndef check(name, ok):\n    report.append((name, bool(ok)))\n    print(f\"  [{'ok' if ok else 'FAIL'}] {name}\")"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[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."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "print(\"[P-def] A = V Lambda U^T with orthonormal V, U\")\nA = rng.normal(size=(5, 3))                     # rectangular\nVf, s, Ut = np.linalg.svd(A)                    # full SVD\nLam = np.zeros((5, 3)); Lam[:3, :3] = np.diag(s)\ncheck(\"reconstruction V Lambda U^T = A (err < 1e-12)\",\n      np.max(np.abs(Vf @ Lam @ Ut - A)) < 1e-12)\ncheck(\"V orthonormal: V^T V = I\", np.allclose(Vf.T @ Vf, np.eye(5)))\ncheck(\"U orthonormal: U^T U = I\", np.allclose(Ut @ Ut.T, np.eye(3)))\noff_diag = Lam.copy(); np.fill_diagonal(off_diag, 0.0)\ncheck(\"Lambda vanishes off the main diagonal\", np.all(off_diag == 0))\ncheck(\"singular values are nonnegative and sorted (descending)\",\n      np.all(s >= 0) and np.all(np.diff(s) <= 0))\ncheck(\"A u^(j) = lambda_j v^(j) for every j\",\n      all(np.linalg.norm(A @ Ut[j] - s[j] * Vf[:, j]) < 1e-12\n          for j in range(3)))"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[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."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "print(\"[P-exist] an SVD exists for every matrix; special cases\")\nD = np.array([[0.0, 1.0], [0.0, 0.0]])\nVd, sd, Utd = np.linalg.svd(D)\ncheck(\"defective [[0,1],[0,0]] (no EVD) still has an exact SVD\",\n      np.max(np.abs(Vd @ np.diag(sd) @ Utd - D)) < 1e-14)\nS = A.T @ A                                     # symmetric psd 3 x 3\nlam_evd = np.sort(np.linalg.eigvalsh(S))[::-1]\nsing_S = np.linalg.svd(S, compute_uv=False)\ncheck(\"symmetric psd: singular values equal the eigenvalues \"\n      \"(SVD = EVD)\", np.allclose(sing_S, lam_evd))\ncheck(\"largest singular value = spectral norm ||A||_2\",\n      np.isclose(s[0], np.linalg.norm(A, 2)))\ncheck(\"lambda_1 / lambda_min = condition number of A\",\n      np.isclose(s[0] / s[-1], np.linalg.cond(A)))"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[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."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "print(\"[P-lowrank] Eckart-Young: the truncated SVD is the best \"\n      \"low-rank approximation\")\nk = 1\nA_k = s[0] * np.outer(Vf[:, 0], Ut[0])\nerr_trunc = np.linalg.norm(A - A_k, 2)\ncheck(\"truncation error ||A - A_1||_2 equals lambda_2\",\n      np.isclose(err_trunc, s[1]))\nbeaten = 0\nfor _ in range(300):\n    a, b = rng.normal(size=5), rng.normal(size=3)\n    B = np.outer(a, b)\n    B *= np.trace(B.T @ A) / np.trace(B.T @ B)   # best scale for this B\n    if np.linalg.norm(A - B, 2) < err_trunc - 1e-9:\n        beaten += 1\ncheck(\"no random rank-1 matrix among 300 samples beats the truncation\",\n      beaten == 0)\nm_px, d_px = 5, 3\ncheck(\"storing the rank-k truncation takes k(m + d + 1) numbers\",\n      k * (m_px + d_px + 1) == k * m_px + k * d_px + k)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[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."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "print(\"[P-pinv] the pseudoinverse from the SVD solves least squares\")\nLam_pinv = np.zeros((3, 5)); Lam_pinv[:3, :3] = np.diag(1.0 / s)\nA_pinv = Ut.T @ Lam_pinv @ Vf.T\ncheck(\"U Lambda^+ V^T matches np.linalg.pinv\",\n      np.allclose(A_pinv, np.linalg.pinv(A)))\nyv = rng.normal(size=5)\nw_hat = A_pinv @ yv\nw_lstsq = np.linalg.lstsq(A, yv, rcond=None)[0]\ncheck(\"A^+ y equals the least-squares solution\",\n      np.allclose(w_hat, w_lstsq))\n\n# ------------------------------------------------------------ preview\nfig, ax = plt.subplots(1, 4, figsize=(12, 2.8))\nfor a, M, t in ((ax[0], A, \"A (5 x 3)\"), (ax[1], Lam, \"Lambda\"),\n                (ax[2], Vf @ Lam @ Ut - A, \"reconstruction error\")):\n    im = a.imshow(M, cmap=\"gray\"); a.set_title(t)\n    fig.colorbar(im, ax=a, shrink=0.75)\nks = np.arange(0, 3)\nerrs = [np.linalg.norm(A - sum(s[j] * np.outer(Vf[:, j], Ut[j])\n                               for j in range(kk)), 2) if kk else\n        np.linalg.norm(A, 2) for kk in ks]\nax[3].plot(ks, errs, \"o-\", c=\"k\", label=\"$\\\\|A - A_k\\\\|_2$\")\nax[3].plot(ks[:-1] + 1, s[1:], \"s\", mfc=\"white\", mec=\"k\",\n           label=\"$\\\\lambda_{k+1}$\")\nax[3].set_xlabel(\"rank $k$\"); ax[3].set_ylabel(\"error\")\nax[3].set_xticks(ks)\nax[3].set_title(\"[P-lowrank] truncation error\")\nax[3].legend(frameon=False)\nfig.tight_layout()\nfig.savefig(OUT_DIR / \"svd.png\", dpi=110)\nprint(f\"\\n{sum(ok for _, ok in report)}/{len(report)} checks passed\")\nassert all(ok for _, ok in report)"
  }
 ]
}