{
 "cells": [
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "---\n",
    "execute:\n",
    "  enabled: false\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-svd-low-rank-md-1 fw=all -->\n",
    "\n",
    "# Singular Value Decomposition and Low-Rank Approximation\n",
    "\n",
    "The eigendecomposition of that section was powerful but\n",
    "*picky*: it needs a square matrix, and to be fully well-behaved (an orthonormal\n",
    "eigenbasis, real eigenvalues) it really wants a *symmetric* one. The defective\n",
    "shear there had a repeated eigenvalue but only a one-dimensional eigenspace, so it\n",
    "admitted no eigenbasis at all. The *singular value decomposition* (SVD) is the\n",
    "same idea made universal. It applies to **every** matrix (rectangular,\n",
    "rank-deficient, non-symmetric, defective) and factors it into a\n",
    "*rotate--scale--rotate* form whose scale factors, the *singular values*, are the\n",
    "right generalization of eigenvalues. The engine is one we already built: the\n",
    "spectral theorem of that section, applied to the\n",
    "symmetric positive semidefinite matrix $\\mathbf{A}^\\top\\mathbf{A}$. From this one\n",
    "factorization we obtain rank and the four fundamental subspaces (resolving the\n",
    "\"efficient rank\" promise of that section), the\n",
    "*optimal* low-rank approximation, principal component analysis, the pseudoinverse\n",
    "and least squares, and the condition number that predicts numerical trouble and\n",
    "gradient-descent zig-zag. It is also the tool behind LoRA and other low-rank\n",
    "adapters, the Muon optimizer, and the spectral diagnostics of the closing\n",
    "section.\n",
    "\n",
    "We use the following imports throughout the section."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "svd-imports"
   },
   "source": [
    "%matplotlib inline\n",
    "from d2l import torch as d2l\n",
    "import numpy as np\n",
    "import torch"
   ],
   "execution_count": 1,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-svd-low-rank-md-2 fw=all -->\n",
    "\n",
    "## The Singular Value Decomposition\n",
    "\n",
    "### Rotate--Scale--Rotate\n",
    "\n",
    "Recall the central picture of that section: a matrix\n",
    "sends the unit circle to an ellipse. For a *symmetric* matrix, the ellipse's axes\n",
    "lie along the eigenvectors, and one orthonormal frame does the whole job. A\n",
    "general matrix bends the circle to an ellipse too, but now the input directions\n",
    "that undergo pure stretching and the output directions they land on are *two\n",
    "different* orthonormal frames. The SVD names both.\n",
    "\n",
    "**Theorem (Singular Value Decomposition).** *Every real $m\\times n$ matrix\n",
    "$\\mathbf{A}$ can be written as*\n",
    "\n",
    "$$\n",
    "\\mathbf{A} = \\mathbf{U}\\,\\boldsymbol{\\Sigma}\\,\\mathbf{V}^\\top,\n",
    "$$\n",
    "\n",
    "*where $\\mathbf{U}\\in\\mathbb{R}^{m\\times m}$ and $\\mathbf{V}\\in\\mathbb{R}^{n\\times\n",
    "n}$ are orthogonal ($\\mathbf{U}^\\top\\mathbf{U}=\\mathbf{I}$,\n",
    "$\\mathbf{V}^\\top\\mathbf{V}=\\mathbf{I}$) and $\\boldsymbol{\\Sigma}\\in\\mathbb{R}^{m\\times\n",
    "n}$ is \"diagonal\" with non-negative entries\n",
    "$\\sigma_1\\ge\\sigma_2\\ge\\cdots\\ge0$ on its main diagonal and zeros elsewhere.* The\n",
    "columns $\\mathbf{v}_i$ of $\\mathbf{V}$ are the *right singular vectors*, the\n",
    "columns $\\mathbf{u}_i$ of $\\mathbf{U}$ are the *left singular vectors*, and the\n",
    "$\\sigma_i$ are the *singular values*.\n",
    "\n",
    "We prove existence in that section; first we read off what it\n",
    "*means*. Because $\\mathbf{V}$ and $\\mathbf{U}$ are orthogonal, they are pure\n",
    "rotations (possibly with a reflection), and the equation decomposes the\n",
    "action of $\\mathbf{A}$ on a vector $\\mathbf{x}$ into three stages applied right to\n",
    "left:\n",
    "\n",
    "1. $\\mathbf{V}^\\top$ **rotates** the input so that the right singular vectors\n",
    "   $\\mathbf{v}_i$ land on the coordinate axes;\n",
    "2. $\\boldsymbol{\\Sigma}$ **scales** axis $i$ by $\\sigma_i$ (and, if\n",
    "   $m\\neq n$, embeds into or projects onto the output dimension);\n",
    "3. $\\mathbf{U}$ **rotates** the scaled axes into the output frame, placing the\n",
    "   stretched axis $i$ along the left singular vector $\\mathbf{u}_i$.\n",
    "\n",
    "Equivalently, reading off column $i$ of $\\mathbf{A}\\mathbf{V}=\\mathbf{U}\\boldsymbol{\\Sigma}$,\n",
    "\n",
    "$$\n",
    "\\mathbf{A}\\mathbf{v}_i = \\sigma_i\\,\\mathbf{u}_i .\n",
    "$$\n",
    "\n",
    "The orthonormal input direction $\\mathbf{v}_i$ is sent to the orthonormal output\n",
    "direction $\\mathbf{u}_i$, stretched by $\\sigma_i$. The unit circle (sphere) of\n",
    "right singular vectors becomes an ellipse (ellipsoid) whose semi-axes are the\n",
    "$\\sigma_i\\mathbf{u}_i$. This is exactly the eigen-picture of\n",
    "that section, generalized: there one frame served as both\n",
    "input and output axes because the map was symmetric; here input frame $\\mathbf{V}$\n",
    "and output frame $\\mathbf{U}$ are allowed to differ, which is precisely what lets\n",
    "the SVD handle every matrix. the figure traces this for a\n",
    "non-symmetric $2\\times2$ matrix, so that $\\mathbf{U}\\neq\\mathbf{V}$ is visible: the\n",
    "unit circle and its right singular vectors $\\mathbf{v}_1,\\mathbf{v}_2$ are rotated\n",
    "by $\\mathbf{V}^\\top$ onto the axes, scaled by $\\boldsymbol{\\Sigma}$, then rotated\n",
    "by $\\mathbf{U}$ so the semi-axes land along the output frame\n",
    "$\\sigma_1\\mathbf{u}_1,\\sigma_2\\mathbf{u}_2$, a different orientation than the\n",
    "input frame, because the matrix rotates as well as stretches.\n",
    "\n",
    "![The action of a non-symmetric $\\mathbf{A}=\\mathbf{U}\\boldsymbol{\\Sigma}\\mathbf{V}^\\top$ on the unit circle, read right to left: rotate by $\\mathbf{V}^\\top$, scale by $\\boldsymbol{\\Sigma}$, rotate by $\\mathbf{U}$. The input frame $\\{\\mathbf{v}_i\\}$ and output frame $\\{\\mathbf{u}_i\\}$ point in different directions.](../img/mdl-la-svd-action.svg)\n",
    "\n",
    "**The polar decomposition: rotate--scale--rotate, rigorously.** The phrase\n",
    "\"rotate--scale--rotate\" can be made precise. For square\n",
    "$\\mathbf{A}$, insert $\\mathbf{V}\\mathbf{V}^\\top=\\mathbf{I}$:\n",
    "\n",
    "$$\n",
    "\\mathbf{A}\n",
    "  = \\mathbf{U}\\boldsymbol{\\Sigma}\\mathbf{V}^\\top\n",
    "  = (\\mathbf{U}\\mathbf{V}^\\top)\\,(\\mathbf{V}\\boldsymbol{\\Sigma}\\mathbf{V}^\\top)\n",
    "  = \\mathbf{Q}\\,\\mathbf{P}.\n",
    "$$\n",
    "\n",
    "Here $\\mathbf{Q}=\\mathbf{U}\\mathbf{V}^\\top$ is orthogonal (a product of\n",
    "orthogonal matrices) and $\\mathbf{P}=\\mathbf{V}\\boldsymbol{\\Sigma}\\mathbf{V}^\\top$\n",
    "is symmetric positive semidefinite (its eigendecomposition has the non-negative\n",
    "$\\sigma_i$ as eigenvalues). This *polar decomposition* says every linear map is a\n",
    "positive-semidefinite stretch $\\mathbf{P}$ followed by a rigid rotation\n",
    "$\\mathbf{Q}$: the exact matrix analog of writing a complex number as\n",
    "$z=re^{i\\theta}$, modulus times phase. It justifies the rotate--scale--rotate\n",
    "slogan: the \"scale\" is $\\mathbf{P}$ and the net \"rotate\" is $\\mathbf{Q}$.\n",
    "\n",
    "**Thin SVD and the dyadic sum.** If $r=\\operatorname{rank}\\mathbf{A}$, only the\n",
    "first $r$ singular values are nonzero. Discarding the zero columns of\n",
    "$\\boldsymbol{\\Sigma}$ and the matching singular vectors gives the *thin* (or\n",
    "*economy*) SVD $\\mathbf{A}=\\mathbf{U}_r\\boldsymbol{\\Sigma}_r\\mathbf{V}_r^\\top$ with\n",
    "$\\mathbf{U}_r\\in\\mathbb{R}^{m\\times r}$, $\\mathbf{V}_r\\in\\mathbb{R}^{n\\times r}$.\n",
    "Multiplying out, the SVD is a sum of rank-one *dyads*,\n",
    "\n",
    "$$\n",
    "\\mathbf{A} = \\sum_{i=1}^{r} \\sigma_i\\,\\mathbf{u}_i\\mathbf{v}_i^\\top,\n",
    "$$\n",
    "\n",
    "each term a single outer product weighted by its singular value, ordered from\n",
    "largest to smallest. This form, a *ranked* list of rank-one ingredients, is the\n",
    "key to the low-rank approximation of that section: keep the\n",
    "heavy terms, drop the light ones.\n",
    "\n",
    "The dyadic sum is also a one-line `einsum` in the Einstein notation of\n",
    "that section.\n",
    "Summing $\\sigma_i\\,\\mathbf{u}_i\\mathbf{v}_i^\\top$ over $i$ is the index pattern\n",
    "`'i,mi,in->mn'` (row $i$ of `Vt` *is* $\\mathbf{v}_i^\\top$), and it rebuilds the\n",
    "matrix exactly:"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "mdl-svd-low-rank-rotate-scale-rotate"
   },
   "source": [
    "rng = np.random.default_rng(5)\n",
    "A = rng.standard_normal((4, 3))\n",
    "U, s, Vt = np.linalg.svd(A, full_matrices=False)   # thin SVD\n",
    "A_dyads = np.einsum('i,mi,in->mn', s, U, Vt)       # sum_i sigma_i u_i v_i^T\n",
    "print('dyadic rebuild error:', round(float(np.linalg.norm(A_dyads - A)), 12))"
   ],
   "execution_count": 2,
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "dyadic rebuild error: 0.0\n"
     ]
    }
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-svd-low-rank-md-3 fw=all -->\n",
    "\n",
    "### Existence via $\\mathbf{A}^\\top\\mathbf{A}$\n",
    "\n",
    "The SVD is the spectral theorem of that section applied\n",
    "to the symmetric PSD matrix $\\mathbf{A}^\\top\\mathbf{A}$, which we built for\n",
    "exactly this purpose in the bridge at the end of that section. The crux of the\n",
    "construction is that orthonormality of the $\\mathbf{u}_i$ comes for free.\n",
    "\n",
    "**Proof of the equation (existence).** The matrix\n",
    "$\\mathbf{A}^\\top\\mathbf{A}$ is symmetric and positive semidefinite (recall\n",
    "that section:\n",
    "$\\mathbf{x}^\\top(\\mathbf{A}^\\top\\mathbf{A})\\mathbf{x} = \\|\\mathbf{A}\\mathbf{x}\\|^2 \\ge 0$).\n",
    "By the spectral theorem it has an orthonormal eigenbasis\n",
    "$\\mathbf{v}_1,\\ldots,\\mathbf{v}_n$ with real eigenvalues, non-negative by\n",
    "positive semidefiniteness; order them\n",
    "$\\lambda_1\\ge\\cdots\\ge\\lambda_n\\ge0$ and set\n",
    "\n",
    "$$\n",
    "\\sigma_i = \\sqrt{\\lambda_i}, \\qquad r = \\#\\{i : \\sigma_i > 0\\} .\n",
    "$$\n",
    "\n",
    "For each $i\\le r$ define $\\mathbf{u}_i = \\mathbf{A}\\mathbf{v}_i/\\sigma_i$. These\n",
    "left singular vectors are automatically orthonormal, and this single line is why:\n",
    "\n",
    "$$\n",
    "\\mathbf{u}_i^\\top\\mathbf{u}_j\n",
    "  = \\frac{1}{\\sigma_i\\sigma_j}\\,\\mathbf{v}_i^\\top\\mathbf{A}^\\top\\mathbf{A}\\,\\mathbf{v}_j\n",
    "  = \\frac{1}{\\sigma_i\\sigma_j}\\,\\mathbf{v}_i^\\top(\\lambda_j\\mathbf{v}_j)\n",
    "  = \\frac{\\lambda_j}{\\sigma_i\\sigma_j}\\,\\delta_{ij}\n",
    "  = \\delta_{ij} .\n",
    "$$\n",
    "\n",
    "Orthonormality of the *output* frame is inherited from orthonormality of the\n",
    "*input* frame, routed through $\\mathbf{A}^\\top\\mathbf{A}$. Extend\n",
    "$\\mathbf{u}_1,\\ldots,\\mathbf{u}_r$ to an orthonormal basis\n",
    "$\\mathbf{u}_1,\\ldots,\\mathbf{u}_m$ of $\\mathbb{R}^m$ (Gram--Schmidt on any\n",
    "completion). For the remaining right singular vectors, $i>r$, we have\n",
    "$\\|\\mathbf{A}\\mathbf{v}_i\\|^2=\\mathbf{v}_i^\\top\\mathbf{A}^\\top\\mathbf{A}\\mathbf{v}_i=\\lambda_i=0$,\n",
    "so $\\mathbf{A}\\mathbf{v}_i=\\mathbf 0$. In all cases, then,\n",
    "$\\mathbf{A}\\mathbf{v}_i=\\sigma_i\\mathbf{u}_i$ (with $\\sigma_i=0$ for $i>r$).\n",
    "Collecting these columns gives $\\mathbf{A}\\mathbf{V}=\\mathbf{U}\\boldsymbol{\\Sigma}$,\n",
    "and right-multiplying by $\\mathbf{V}^\\top=\\mathbf{V}^{-1}$ yields\n",
    "$\\mathbf{A}=\\mathbf{U}\\boldsymbol{\\Sigma}\\mathbf{V}^\\top$. $\\blacksquare$\n",
    "\n",
    "This proof also *computes* the SVD by hand: diagonalize the (small, symmetric)\n",
    "matrix $\\mathbf{A}^\\top\\mathbf{A}$, take square roots for the singular values, and\n",
    "push the eigenvectors through $\\mathbf{A}$ for the left singular vectors.\n",
    "\n",
    "**Uniqueness.** The singular *values* are unique: they are the square roots of\n",
    "the eigenvalues of $\\mathbf{A}^\\top\\mathbf{A}$, which are determined by\n",
    "$\\mathbf{A}$. The singular *vectors* are unique only up to the same ambiguities we\n",
    "met for eigenvectors in that section: a sign flip on a\n",
    "pair $(\\mathbf{u}_i,\\mathbf{v}_i)$ for a simple $\\sigma_i$, and an arbitrary\n",
    "orthonormal rotation within the subspace belonging to a repeated singular value.\n",
    "\n",
    "**Variational characterization of $\\sigma_1$.** There is a second route to the\n",
    "top singular value that needs no eigendecomposition and\n",
    "makes the *meaning* of $\\sigma_1$ plain. It is the Rayleigh quotient of\n",
    "that section in disguise:\n",
    "\n",
    "$$\n",
    "\\sigma_1 = \\max_{\\|\\mathbf{x}\\|=1}\\|\\mathbf{A}\\mathbf{x}\\| = \\|\\mathbf{A}\\|_2 ,\n",
    "$$\n",
    "\n",
    "because $\\|\\mathbf{A}\\mathbf{x}\\|^2=\\mathbf{x}^\\top\\mathbf{A}^\\top\\mathbf{A}\\mathbf{x}$\n",
    "is maximized over unit $\\mathbf{x}$ by the Rayleigh proposition at\n",
    "$\\lambda_1=\\sigma_1^2$, attained at $\\mathbf{x}=\\mathbf{v}_1$. So **$\\sigma_1$ is\n",
    "the most the matrix can stretch any unit vector**: its operator (spectral) norm\n",
    "$\\|\\mathbf{A}\\|_2$. The remaining singular values come from the same problem with\n",
    "deflation: $\\sigma_k=\\max\\{\\|\\mathbf{A}\\mathbf{x}\\| : \\|\\mathbf{x}\\|=1,\\ \\mathbf{x}\\perp\\mathbf{v}_1,\\ldots,\\mathbf{v}_{k-1}\\}$,\n",
    "the Courant--Fischer min-max principle of that section translated\n",
    "from Rayleigh quotients to singular values.\n",
    "This \"maximum stretch\" reading is the one that recurs in Eckart--Young, PCA,\n",
    "conditioning, and Lipschitz/spectral-norm arguments. The constructive proof is the\n",
    "better *first* proof (concrete, reusing machinery we have); the variational one is\n",
    "the better *meaning*.\n",
    "\n",
    "**Relationship to the eigendecomposition.** Substituting the equation into\n",
    "the two Gram matrices gives, since $\\mathbf{U}$ and $\\mathbf{V}$ are orthogonal,\n",
    "\n",
    "$$\n",
    "\\mathbf{A}^\\top\\mathbf{A} = \\mathbf{V}\\,\\boldsymbol{\\Sigma}^\\top\\boldsymbol{\\Sigma}\\,\\mathbf{V}^\\top,\n",
    "\\qquad\n",
    "\\mathbf{A}\\mathbf{A}^\\top = \\mathbf{U}\\,\\boldsymbol{\\Sigma}\\boldsymbol{\\Sigma}^\\top\\,\\mathbf{U}^\\top .\n",
    "$$\n",
    "\n",
    "These are *eigendecompositions*: the right singular vectors are the eigenvectors\n",
    "of $\\mathbf{A}^\\top\\mathbf{A}$, the left singular vectors are the eigenvectors of\n",
    "$\\mathbf{A}\\mathbf{A}^\\top$, and the squared singular values $\\sigma_i^2$ are the\n",
    "shared eigenvalues. For a symmetric PSD matrix the SVD and the eigendecomposition\n",
    "coincide. In general they differ, and the warning example against conflating\n",
    "singular values with eigenvalues is the scaled rotation\n",
    "\n",
    "$$\n",
    "\\mathbf{A} = \\begin{bmatrix} 0 & -2\\\\ 1 & \\phantom{-}0\\end{bmatrix},\n",
    "\\qquad\n",
    "\\mathbf{A}^\\top\\mathbf{A} = \\begin{bmatrix} 1 & 0\\\\ 0 & 4\\end{bmatrix},\n",
    "$$\n",
    "\n",
    "whose singular values are $\\{2,1\\}$ (the square roots of $\\{4,1\\}$), while its\n",
    "eigenvalues are $\\pm i\\sqrt2$ with modulus $|\\lambda|=\\sqrt2$. The singular values\n",
    "are *not* the eigenvalue magnitudes; here even their *geometric mean*\n",
    "$\\sqrt{\\sigma_1\\sigma_2}=\\sqrt2$ equals $|\\lambda|$, which is no accident: it is\n",
    "the matrix analog of the rotation-with-scaling reading from\n",
    "that section, with $|\\det\\mathbf{A}|=\\sigma_1\\sigma_2=|\\lambda_1\\lambda_2|=2$.\n",
    "The determinant identity holds in general: for any square $\\mathbf{A}$,\n",
    "$|\\det\\mathbf{A}|=\\prod_i\\sigma_i$, because orthogonal factors have determinant\n",
    "$\\pm1$, so the volume scaling of\n",
    "that section is carried entirely by the\n",
    "diagonal stretch $\\boldsymbol{\\Sigma}$.\n",
    "\n",
    "**The defective shear, finally decomposed.** The shear of\n",
    "that section\n",
    "\n",
    "$$\n",
    "\\mathbf{A} = \\begin{bmatrix} 1 & 1\\\\ 0 & 1\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "is *defective*: $\\lambda=1$ has algebraic multiplicity $2$ but only a\n",
    "one-dimensional eigenspace, so it has **no eigenbasis** and the eigendecomposition\n",
    "simply does not exist. The SVD has no such trouble. Form\n",
    "\n",
    "$$\n",
    "\\mathbf{A}^\\top\\mathbf{A} = \\begin{bmatrix} 1 & 1\\\\ 1 & 2\\end{bmatrix},\n",
    "$$\n",
    "\n",
    "a symmetric PSD matrix. Its characteristic polynomial is\n",
    "$\\lambda^2-3\\lambda+1$, with roots $\\lambda_{1,2}=(3\\pm\\sqrt5)/2$. Hence the\n",
    "singular values are\n",
    "\n",
    "$$\n",
    "\\sigma_1 = \\sqrt{\\tfrac{3+\\sqrt5}{2}} = \\frac{1+\\sqrt5}{2} = \\varphi \\approx 1.618,\n",
    "\\qquad\n",
    "\\sigma_2 = \\sqrt{\\tfrac{3-\\sqrt5}{2}} = \\frac{1}{\\varphi} = \\varphi-1 \\approx 0.618,\n",
    "$$\n",
    "\n",
    "the golden ratio and its reciprocal (consistent with\n",
    "$\\sigma_1\\sigma_2=|\\det\\mathbf{A}|=1$). Both are strictly positive, so the shear\n",
    "has full rank $2$ and a clean orthonormal factorization\n",
    "$\\mathbf{A}=\\mathbf{U}\\boldsymbol{\\Sigma}\\mathbf{V}^\\top$: there is nothing\n",
    "defective about it. The eigendecomposition stumbled because it insists on a\n",
    "*single* basis that both diagonalizes and is reused for input and output; the SVD\n",
    "succeeds because it is allowed two different orthonormal frames, which is exactly\n",
    "what the shear needs. The cell below verifies the singular values against\n",
    "the golden ratio."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "svd-defective-shear"
   },
   "source": [
    "A = np.array([[1., 1.], [0., 1.]])\n",
    "s = np.linalg.svd(A, compute_uv=False)\n",
    "phi = (1 + np.sqrt(5)) / 2\n",
    "print('singular values of the shear:', s.round(6))\n",
    "print('golden ratio phi, 1/phi:    ', np.array([phi, 1 / phi]).round(6))\n",
    "print('product sigma_1 * sigma_2 = |det A| =', round(float(s[0] * s[1]), 6))"
   ],
   "execution_count": 3,
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "singular values of the shear: [1.618034 0.618034]\ngolden ratio phi, 1/phi:     [1.618034 0.618034]\nproduct sigma_1 * sigma_2 = |det A| = 1.0\n"
     ]
    }
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-svd-low-rank-md-4 fw=all -->\n",
    "\n",
    "### The Four Fundamental Subspaces\n",
    "\n",
    "Because $\\mathbf{U}$ and $\\mathbf{V}$ are invertible, multiplying by them changes\n",
    "no dimensions, so $\\operatorname{rank}\\mathbf{A}=\\operatorname{rank}\\boldsymbol{\\Sigma}$,\n",
    "which is simply the number of nonzero diagonal entries:\n",
    "\n",
    "$$\n",
    "\\operatorname{rank}\\mathbf{A} = \\#\\{i : \\sigma_i > 0\\} = r .\n",
    "$$\n",
    "\n",
    "This is the *efficient* characterization of rank promised in\n",
    "that section, where counting independent\n",
    "columns by elimination was the only tool on offer. The dyadic form\n",
    "the equation makes the four fundamental subspaces equally explicit.\n",
    "Splitting the singular vectors at the index $r$,\n",
    "\n",
    "* the **range** (column space) of $\\mathbf{A}$ is\n",
    "  $\\operatorname{span}\\{\\mathbf{u}_1,\\ldots,\\mathbf{u}_r\\}$, since every output\n",
    "  $\\mathbf{A}\\mathbf{x}=\\sum_{i\\le r}\\sigma_i(\\mathbf{v}_i^\\top\\mathbf{x})\\mathbf{u}_i$\n",
    "  is a combination of these;\n",
    "* the **null space** of $\\mathbf{A}$ is\n",
    "  $\\operatorname{span}\\{\\mathbf{v}_{r+1},\\ldots,\\mathbf{v}_n\\}$, the input\n",
    "  directions sent to zero;\n",
    "* the **row space** (range of $\\mathbf{A}^\\top$) is\n",
    "  $\\operatorname{span}\\{\\mathbf{v}_1,\\ldots,\\mathbf{v}_r\\}$;\n",
    "* the **left null space** (null space of $\\mathbf{A}^\\top$) is\n",
    "  $\\operatorname{span}\\{\\mathbf{u}_{r+1},\\ldots,\\mathbf{u}_m\\}$.\n",
    "\n",
    "The picture is a clean bijection: $\\mathbf{A}$ maps the row space to the column\n",
    "space, sending each $\\mathbf{v}_i\\mapsto\\sigma_i\\mathbf{u}_i$ ($i\\le r$) one-to-one\n",
    "and onto, while crushing the null space to zero. Input space splits orthogonally\n",
    "as row space $\\oplus$ null space; output space splits orthogonally as column space\n",
    "$\\oplus$ left-null space. the figure draws this two-plane\n",
    "map; the pseudoinverse of that section will run the bijection\n",
    "backwards.\n",
    "\n",
    "![The four fundamental subspaces. $\\mathbf{A}$ maps the row space $\\operatorname{span}\\{\\mathbf{v}_1,\\dots,\\mathbf{v}_r\\}$ bijectively onto the column space $\\operatorname{span}\\{\\mathbf{u}_1,\\dots,\\mathbf{u}_r\\}$ via $\\mathbf{v}_i\\mapsto\\sigma_i\\mathbf{u}_i$, and crushes the null space $\\operatorname{span}\\{\\mathbf{v}_{r+1},\\dots,\\mathbf{v}_n\\}$ to zero.](../img/mdl-la-svd-subspaces.svg)\n",
    "\n",
    "**Numerical rank.** In floating-point arithmetic, a matrix that is mathematically\n",
    "rank-deficient rarely has exact zero singular values; rounding leaves tiny\n",
    "$\\sigma_i$ of size around $\\epsilon_{\\text{mach}}\\,\\sigma_1$ instead. The\n",
    "practical notion of rank therefore *thresholds*: count the singular values above a\n",
    "tolerance, which is exactly what `np.linalg.matrix_rank` does; its default cutoff is\n",
    "$\\sigma_1\\,\\max(m,n)\\,\\epsilon_{\\text{mach}}$, scaled by both the largest singular\n",
    "value and the matrix size. Building a deliberately rank-2\n",
    "matrix in $\\mathbb{R}^{4\\times4}$ shows two singular values collapse to near\n",
    "machine zero, so a tolerance recovers the true rank where a test for exact zeros\n",
    "would fail."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "svd-numerical-rank"
   },
   "source": [
    "rng = np.random.default_rng(0)\n",
    "A = rng.standard_normal((4, 2)) @ rng.standard_normal((2, 4))  # rank <= 2\n",
    "s = np.linalg.svd(A, compute_uv=False)\n",
    "print('singular values:', s.round(6))\n",
    "print('numerical rank :', int(np.linalg.matrix_rank(A)))"
   ],
   "execution_count": 4,
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "singular values: [4.351705 0.914865 0.       0.      ]\nnumerical rank : 2\n"
     ]
    }
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-svd-low-rank-md-5 fw=all -->\n",
    "\n",
    "## Low-Rank Approximation\n",
    "\n",
    "### Eckart--Young\n",
    "\n",
    "Here is the theorem that makes the SVD indispensable. The dyadic sum\n",
    "the equation lists the rank-one pieces of $\\mathbf{A}$ in order of\n",
    "importance. Keeping only the top $k$,\n",
    "\n",
    "$$\n",
    "\\mathbf{A}_k = \\sum_{i=1}^{k}\\sigma_i\\,\\mathbf{u}_i\\mathbf{v}_i^\\top,\n",
    "$$\n",
    "\n",
    "is the best rank-$k$ approximation of $\\mathbf{A}$, provably, in both the\n",
    "spectral and Frobenius norms. Eckart and Young proved the\n",
    "Frobenius case [@Eckart.Young.1936] and Mirsky extended it to every\n",
    "*unitarily invariant* norm (one unchanged by orthogonal multiplication on either\n",
    "side) [@Mirsky.1960]; we give the spectral-norm statement and proof\n",
    "in full.\n",
    "\n",
    "**Theorem (Eckart--Young--Mirsky).** *For every $k<r$,*\n",
    "\n",
    "$$\n",
    "\\min_{\\operatorname{rank}\\mathbf{B}\\le k}\\|\\mathbf{A}-\\mathbf{B}\\|_2 = \\sigma_{k+1},\n",
    "$$\n",
    "\n",
    "*and the minimum is attained at $\\mathbf{B}=\\mathbf{A}_k$. In the Frobenius norm,\n",
    "$\\min_{\\operatorname{rank}\\mathbf{B}\\le k}\\|\\mathbf{A}-\\mathbf{B}\\|_F^2=\\sum_{i>k}\\sigma_i^2$,\n",
    "again attained at $\\mathbf{A}_k$.*\n",
    "\n",
    "**Proof.** *The truncation achieves $\\sigma_{k+1}$.* The error\n",
    "$\\mathbf{A}-\\mathbf{A}_k=\\sum_{i>k}\\sigma_i\\mathbf{u}_i\\mathbf{v}_i^\\top$ is itself\n",
    "a (sub-)SVD whose largest singular value is $\\sigma_{k+1}$, so by the variational\n",
    "identity the equation,\n",
    "$\\|\\mathbf{A}-\\mathbf{A}_k\\|_2=\\sigma_{k+1}$. (And\n",
    "$\\|\\mathbf{A}-\\mathbf{A}_k\\|_F^2=\\sum_{i>k}\\sigma_i^2$, by the identity\n",
    "$\\|\\mathbf{A}\\|_F^2=\\operatorname{tr}(\\mathbf{A}^\\top\\mathbf{A})=\\sum_i\\sigma_i^2$:\n",
    "the Frobenius norm is the root-sum-of-squares of the singular values. Exercise 2\n",
    "asks you to prove this in two lines.)\n",
    "\n",
    "*No rank-$k$ matrix does better (the dimension-counting argument).* Let\n",
    "$\\mathbf{B}$ be any matrix with $\\operatorname{rank}\\mathbf{B}\\le k$. Its null\n",
    "space has dimension at least $n-k$. Consider also the\n",
    "$(k{+}1)$-dimensional subspace\n",
    "$\\mathcal{V}=\\operatorname{span}\\{\\mathbf{v}_1,\\ldots,\\mathbf{v}_{k+1}\\}$ spanned by\n",
    "the top $k{+}1$ right singular vectors. Two subspaces of $\\mathbb{R}^n$ whose\n",
    "dimensions add to more than $n$ must intersect nontrivially:\n",
    "\n",
    "$$\n",
    "\\dim(\\ker\\mathbf{B}\\cap\\mathcal{V}) \\ge (n-k) + (k+1) - n = 1 .\n",
    "$$\n",
    "\n",
    "So there is a *unit* vector $\\mathbf{x}\\in\\ker\\mathbf{B}\\cap\\mathcal{V}$. Write it\n",
    "in the top singular basis, $\\mathbf{x}=\\sum_{i\\le k+1}c_i\\mathbf{v}_i$ with\n",
    "$\\sum c_i^2=1$. Because $\\mathbf{x}\\in\\ker\\mathbf{B}$ we have $\\mathbf{B}\\mathbf{x}=\\mathbf 0$,\n",
    "and because $\\mathbf{A}\\mathbf{v}_i=\\sigma_i\\mathbf{u}_i$ with the $\\mathbf{u}_i$\n",
    "orthonormal,\n",
    "\n",
    "$$\n",
    "\\|\\mathbf{A}\\mathbf{x}\\|^2 = \\Bigl\\|\\sum_{i\\le k+1}c_i\\sigma_i\\mathbf{u}_i\\Bigr\\|^2\n",
    "   = \\sum_{i\\le k+1} c_i^2\\sigma_i^2\n",
    "   \\ge \\sigma_{k+1}^2 \\sum_{i\\le k+1} c_i^2\n",
    "   = \\sigma_{k+1}^2 ,\n",
    "$$\n",
    "\n",
    "using $\\sigma_i\\ge\\sigma_{k+1}$ for $i\\le k+1$. Therefore\n",
    "\n",
    "$$\n",
    "\\|\\mathbf{A}-\\mathbf{B}\\|_2^2\n",
    "  \\ge \\|(\\mathbf{A}-\\mathbf{B})\\mathbf{x}\\|^2\n",
    "  = \\|\\mathbf{A}\\mathbf{x}\\|^2 \\ge \\sigma_{k+1}^2 ,\n",
    "$$\n",
    "\n",
    "so $\\|\\mathbf{A}-\\mathbf{B}\\|_2\\ge\\sigma_{k+1}$ for every rank-$k$ $\\mathbf{B}$,\n",
    "with equality at $\\mathbf{A}_k$. $\\blacksquare$\n",
    "\n",
    "The kernel of $\\mathbf{B}$ and the top-$(k{+}1)$\n",
    "singular subspace *overfill* $\\mathbb{R}^n$, so they must share a direction. On\n",
    "that shared vector $\\mathbf{B}$ is blind ($\\mathbf{B}\\mathbf{x}=\\mathbf 0$) while\n",
    "$\\mathbf{A}$ still stretches by at least $\\sigma_{k+1}$, so no rank-$k$\n",
    "$\\mathbf{B}$ can track $\\mathbf{A}$ everywhere.\n",
    "\n",
    "For the Frobenius norm, the case PCA relies on, the same conclusion follows from\n",
    "a one-line singular-value inequality, which we state as a lemma and cite rather\n",
    "than prove.\n",
    "\n",
    "**Lemma (Weyl's inequality for singular values** [@Horn.Johnson.2012]**).**\n",
    "*For matrices $\\mathbf{X},\\mathbf{Y}$ of the same shape and all valid indices\n",
    "$i,j$,*\n",
    "\n",
    "$$\n",
    "\\sigma_{i+j-1}(\\mathbf{X}+\\mathbf{Y})\\le\\sigma_i(\\mathbf{X})+\\sigma_j(\\mathbf{Y}) .\n",
    "$$\n",
    "\n",
    "**Proof (Frobenius case).** The lemma\n",
    "applied with $\\mathbf{X}=\\mathbf{A}-\\mathbf{B}$, $\\mathbf{Y}=\\mathbf{B}$ and\n",
    "$j=k{+}1$ gives\n",
    "$\\sigma_{k+i}(\\mathbf{A})\\le\\sigma_i(\\mathbf{A}-\\mathbf{B})+\\sigma_{k+1}(\\mathbf{B})=\\sigma_i(\\mathbf{A}-\\mathbf{B})$,\n",
    "since $\\mathbf{B}$ has rank $\\le k$ and so $\\sigma_{k+1}(\\mathbf{B})=0$. Squaring\n",
    "and summing over $i\\ge1$,\n",
    "\n",
    "$$\n",
    "\\|\\mathbf{A}-\\mathbf{B}\\|_F^2=\\sum_{i\\ge1}\\sigma_i(\\mathbf{A}-\\mathbf{B})^2\n",
    "  \\ge\\sum_{i\\ge1}\\sigma_{k+i}(\\mathbf{A})^2=\\sum_{j>k}\\sigma_j(\\mathbf{A})^2 ,\n",
    "$$\n",
    "\n",
    "again with equality at the truncation $\\mathbf{A}_k$. $\\blacksquare$\n",
    "\n",
    "The error formulas give a quantitative *dial*. The fraction of \"energy\" retained\n",
    "by the rank-$k$ truncation is the *energy ratio*\n",
    "\n",
    "$$\n",
    "\\frac{\\sum_{i\\le k}\\sigma_i^2}{\\sum_{i}\\sigma_i^2}\n",
    "  = 1 - \\frac{\\|\\mathbf{A}-\\mathbf{A}_k\\|_F^2}{\\|\\mathbf{A}\\|_F^2} ,\n",
    "$$\n",
    "\n",
    "so choosing $k$ to capture, say, 95% of the energy is a principled way to set the\n",
    "rank. When the singular values decay quickly (as they do for images with\n",
    "large-scale structure and, empirically, for many trained weight matrices), a small\n",
    "$k$ captures almost everything, which is exactly the regime in which low-rank\n",
    "compression pays off.\n",
    "\n",
    "**Truncation as denoising.** There is a second, statistical reason to truncate.\n",
    "Suppose the matrix you observe is a low-rank signal plus noise,\n",
    "$\\mathbf{A}=\\mathbf{B}+\\mathbf{N}$, with $\\operatorname{rank}\\mathbf{B}=r$ small\n",
    "and $\\mathbf{N}$ having i.i.d. entries of standard deviation\n",
    "$\\sigma_{\\text{noise}}$. Then the spectrum of $\\mathbf{A}$ *splits*: the $r$\n",
    "signal values stand essentially where they were, while the noise contributes a\n",
    "floor of singular values clustered below roughly\n",
    "$(\\sqrt{m}+\\sqrt{n})\\,\\sigma_{\\text{noise}}$. Truncating just above that floor\n",
    "discards almost pure noise, so $\\mathbf{A}_k$ can be *closer to the truth*\n",
    "$\\mathbf{B}$ than the observed $\\mathbf{A}$ is. Gavish and Donoho made the cutoff\n",
    "precise: for an $n\\times n$ matrix with known noise level, the asymptotically\n",
    "optimal hard threshold is\n",
    "$\\tfrac{4}{\\sqrt3}\\sqrt{n}\\,\\sigma_{\\text{noise}}\\approx2.309\\,\\sqrt{n}\\,\\sigma_{\\text{noise}}$\n",
    "[@Gavish.Donoho.2014], a principled alternative to the 95%-energy dial when\n",
    "the data is noisy. Exercise 3 lets you watch the spectrum split.\n",
    "\n",
    "The same low-rank premise, with *missing* rather than noisy entries, is *matrix\n",
    "completion*: recommender systems fill in a sparsely observed ratings matrix by\n",
    "seeking the lowest-rank matrix consistent with the observed entries, in practice\n",
    "by minimizing the *nuclear norm* $\\|\\mathbf{A}\\|_*=\\sum_i\\sigma_i$, the convex\n",
    "stand-in for rank [@Candes.Recht.2009]. The nuclear norm plays the role for\n",
    "rank that the $\\ell_1$ norm plays for sparsity.\n",
    "\n",
    "the figure makes this concrete on a grayscale image. The\n",
    "left panel plots the singular-value spectrum on a log scale (note the rapid decay);\n",
    "the remaining panels reconstruct the image at ranks $k=1,5,20$ and full, each\n",
    "labeled with its relative Frobenius error\n",
    "$\\|\\mathbf{A}-\\mathbf{A}_k\\|_F/\\|\\mathbf{A}\\|_F=\\sqrt{\\sum_{i>k}\\sigma_i^2/\\sum_i\\sigma_i^2}$.\n",
    "A rank-20 truncation of this image already looks essentially correct while storing\n",
    "only a fraction of the numbers, because the discarded singular values carry little\n",
    "energy: a visual proof of Eckart--Young.\n",
    "\n",
    "![Eckart--Young on a grayscale image: the singular-value spectrum on a log scale (left, note the rapid decay) and the rank-$k$ reconstructions $\\mathbf{A}_k$ for $k=1,5,20$ and full rank, each labeled with its relative Frobenius error.](../img/mdl-la-eckart-young.svg)\n",
    "\n",
    "### Principal Component Analysis\n",
    "\n",
    "Principal component analysis is the most important single application of the SVD,\n",
    "and Eckart--Young turns it from a recipe into a theorem. Given data points as the\n",
    "rows of $\\mathbf{X}\\in\\mathbb{R}^{n\\times d}$, first *center* them by subtracting\n",
    "the mean row, $\\tilde{\\mathbf{X}}=\\mathbf{X}-\\mathbf 1\\bar{\\mathbf{x}}^\\top$. The\n",
    "empirical covariance is the symmetric PSD matrix\n",
    "$\\mathbf{C}=\\tfrac1n\\tilde{\\mathbf{X}}^\\top\\tilde{\\mathbf{X}}$ from\n",
    "that section (we use the maximum-likelihood $\\tfrac1n$ covariance; the\n",
    "unbiased $\\tfrac1{n-1}$ version rescales every eigenvalue by the same factor and so\n",
    "changes neither the principal directions nor the explained-variance ratio). PCA\n",
    "asks: which unit direction $\\mathbf{w}$ captures the most variance of the projected\n",
    "data?\n",
    "\n",
    "**Proposition (PCA via Rayleigh).** *The projected variance\n",
    "$\\tfrac1n\\|\\tilde{\\mathbf{X}}\\mathbf{w}\\|^2=\\mathbf{w}^\\top\\mathbf{C}\\mathbf{w}$ is\n",
    "maximized over unit $\\mathbf{w}$ by the top eigenvector of $\\mathbf{C}$, which is\n",
    "the top right singular vector $\\mathbf{v}_1$ of $\\tilde{\\mathbf{X}}$, with maximal\n",
    "variance $\\sigma_1^2/n$.*\n",
    "\n",
    "**Proof.** The projection of a centered point onto $\\mathbf{w}$ has coordinate\n",
    "$\\tilde{\\mathbf{x}}^\\top\\mathbf{w}$; the variance of these coordinates is\n",
    "$\\tfrac1n\\sum_i(\\tilde{\\mathbf{x}}_i^\\top\\mathbf{w})^2=\\tfrac1n\\|\\tilde{\\mathbf{X}}\\mathbf{w}\\|^2=\\mathbf{w}^\\top\\mathbf{C}\\mathbf{w}$.\n",
    "Maximizing this Rayleigh quotient over unit $\\mathbf{w}$ gives, by\n",
    "that section, the top eigenvector of $\\mathbf{C}$ and the value\n",
    "$\\lambda_1(\\mathbf{C})$. By the equation,\n",
    "$\\mathbf{C}=\\tfrac1n\\mathbf{V}\\boldsymbol{\\Sigma}^\\top\\boldsymbol{\\Sigma}\\mathbf{V}^\\top$,\n",
    "so its eigenvectors are the right singular vectors $\\mathbf{v}_i$ of\n",
    "$\\tilde{\\mathbf{X}}$ and its eigenvalues are $\\sigma_i^2/n$. $\\blacksquare$\n",
    "\n",
    "Iterating with deflation, the top-$k$ principal directions are\n",
    "$\\mathbf{v}_1,\\ldots,\\mathbf{v}_k$, and the variance *explained* by component $i$\n",
    "is $\\sigma_i^2/n$. The ranked list of these variances is the *scree plot*, and the\n",
    "explained-variance ratio is exactly the energy ratio\n",
    "the equation of $\\tilde{\\mathbf{X}}$. Moreover the rank-$k$\n",
    "projection $\\mathbf{z}=\\mathbf{V}_k^\\top(\\mathbf{x}-\\bar{\\mathbf{x}})$ is the\n",
    "*optimal* linear dimensionality reduction: minimizing reconstruction error over\n",
    "all rank-$k$ linear maps is Eckart--Young--Frobenius applied to $\\tilde{\\mathbf{X}}$,\n",
    "whose solution is the top-$k$ truncation. **Centering matters**: skip it and the\n",
    "first singular vector is dominated by the offset of the cloud from the origin,\n",
    "capturing where the data *is* rather than how it *varies*.\n",
    "\n",
    "the figure draws a 2-D correlated cloud with its two principal axes\n",
    "scaled by the per-axis standard deviation $\\sigma_i/\\sqrt{n}$: the first axis\n",
    "aligns with the long direction of the cloud, exactly the direction of maximal\n",
    "variance the proposition predicts. The cell below cross-checks the SVD principal\n",
    "axes against the eigenvectors of the covariance and reports the\n",
    "explained-variance ratio; the two routes agree, as the equation\n",
    "guarantees.\n",
    "\n",
    "![A 2-D correlated cloud with its principal axes drawn from the mean of the cloud, each scaled by the standard deviation $\\sigma_i/\\sqrt{n}$ along that direction. The first axis aligns with the direction of maximal variance.](../img/mdl-la-pca.svg)"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "svd-pca"
   },
   "source": [
    "rng = np.random.default_rng(1)\n",
    "theta = np.pi / 5\n",
    "rot = np.array([[np.cos(theta), -np.sin(theta)],\n",
    "                [np.sin(theta),  np.cos(theta)]])\n",
    "X = rng.standard_normal((300, 2)) * np.array([3.0, 0.8]) @ rot.T  # correlated\n",
    "Xc = X - X.mean(0)                                                # center\n",
    "s = np.linalg.svd(Xc, compute_uv=False)\n",
    "eigval = np.linalg.eigvalsh((Xc.T @ Xc) / len(Xc))[::-1]          # descending\n",
    "print('explained variance (sigma^2/n):', (s ** 2 / len(Xc)).round(4))\n",
    "print('eigenvalues of covariance     :', eigval.round(4))\n",
    "print('explained-variance ratio      :', (s ** 2 / (s ** 2).sum()).round(4))"
   ],
   "execution_count": 5,
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "explained variance (sigma^2/n): [8.0873 0.5705]\neigenvalues of covariance     : [8.0873 0.5705]\nexplained-variance ratio      : [0.9341 0.0659]\n"
     ]
    }
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-svd-low-rank-md-6 fw=all -->\n",
    "\n",
    "## Solving Linear Systems with the SVD\n",
    "\n",
    "### The Pseudoinverse and Least Squares\n",
    "\n",
    "When $\\mathbf{A}\\mathbf{x}=\\mathbf{b}$ has no solution (too many equations) or\n",
    "infinitely many (too few), the right formulation is the least-squares problem\n",
    "the equation of that section,\n",
    "and the SVD solves it in full generality. Define the\n",
    "*Moore--Penrose pseudoinverse* by inverting the nonzero singular values and\n",
    "transposing the rotations,\n",
    "\n",
    "$$\n",
    "\\mathbf{A}^{+} = \\mathbf{V}\\boldsymbol{\\Sigma}^{+}\\mathbf{U}^\\top,\n",
    "\\qquad\n",
    "\\boldsymbol{\\Sigma}^{+}_{ii} = \\begin{cases} 1/\\sigma_i & \\sigma_i>0,\\\\ 0 & \\sigma_i=0.\\end{cases}\n",
    "$$\n",
    "\n",
    "**Proposition (min-norm least squares).** *Among all minimizers of\n",
    "$\\|\\mathbf{A}\\mathbf{x}-\\mathbf{b}\\|_2$, the one of smallest norm is\n",
    "$\\hat{\\mathbf{x}}=\\mathbf{A}^{+}\\mathbf{b}$.*\n",
    "\n",
    "**Proof.** Rotate both the input and the target into the singular bases:\n",
    "$\\mathbf{y}=\\mathbf{V}^\\top\\mathbf{x}$ and $\\mathbf{c}=\\mathbf{U}^\\top\\mathbf{b}$.\n",
    "Since $\\mathbf{U}$ is orthogonal it preserves norms, so\n",
    "\n",
    "$$\n",
    "\\|\\mathbf{A}\\mathbf{x}-\\mathbf{b}\\|^2\n",
    "   = \\|\\mathbf{U}^\\top(\\mathbf{A}\\mathbf{x}-\\mathbf{b})\\|^2\n",
    "   = \\|\\boldsymbol{\\Sigma}\\mathbf{y}-\\mathbf{c}\\|^2\n",
    "   = \\sum_{i\\le r}(\\sigma_i y_i - c_i)^2 + \\sum_{i> r} c_i^2 .\n",
    "$$\n",
    "\n",
    "The rotation has *decoupled* the problem into independent scalar problems. The\n",
    "first sum is minimized term by term at $y_i=c_i/\\sigma_i$ for $i\\le r$; the second\n",
    "sum is a constant we cannot touch (it is the irreducible residual, the part of\n",
    "$\\mathbf{b}$ outside the column space). The coordinates $y_i$ for $i>r$ are *free*:\n",
    "they do not appear in the residual at all. Setting them to zero gives the unique\n",
    "solution of *smallest* norm (since $\\|\\mathbf{x}\\|=\\|\\mathbf{y}\\|$). That choice is\n",
    "exactly $\\mathbf{y}=\\boldsymbol{\\Sigma}^{+}\\mathbf{c}$, i.e.\n",
    "$\\hat{\\mathbf{x}}=\\mathbf{V}\\boldsymbol{\\Sigma}^{+}\\mathbf{U}^\\top\\mathbf{b}=\\mathbf{A}^{+}\\mathbf{b}$.\n",
    "$\\blacksquare$\n",
    "\n",
    "Rotating into the SVD basis turns a coupled\n",
    "least-squares problem into a list of one-dimensional problems, and \"small residual\"\n",
    "and \"small norm\" land on *disjoint* coordinate blocks ($i\\le r$ versus $i>r$), so\n",
    "both can be satisfied at once. Two corollaries follow. When $\\mathbf{A}$ is\n",
    "square and invertible, all $\\sigma_i>0$ and\n",
    "$\\mathbf{A}^{+}=\\mathbf{V}\\boldsymbol{\\Sigma}^{-1}\\mathbf{U}^\\top=\\mathbf{A}^{-1}$.\n",
    "And *truncating* the pseudoinverse (dropping terms with tiny $\\sigma_i$ instead of\n",
    "dividing by them) caps the dangerous $1/\\sigma_i$ blow-up; this is a form of\n",
    "regularization, closely related to ridge regression, which we revisit when we\n",
    "discuss weight decay in that section.\n",
    "\n",
    "The classical alternative, the *normal equations*\n",
    "$\\mathbf{A}^\\top\\mathbf{A}\\mathbf{x}=\\mathbf{A}^\\top\\mathbf{b}$, is mathematically\n",
    "equivalent for full-rank $\\mathbf{A}$ but numerically worse: forming\n",
    "$\\mathbf{A}^\\top\\mathbf{A}$ *squares* the condition number\n",
    "(that section), so an SVD- or QR-based solve is preferred.\n",
    "The cell below solves an overdetermined system via `pinv` and via the library\n",
    "`lstsq` (the two agree on the same minimum-residual fit) and prints `cond(A)`\n",
    "alongside `cond(A^T A)` to show the normal-equations matrix is far worse\n",
    "conditioned."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "svd-least-squares"
   },
   "source": [
    "A = np.array([[1., 1.], [1., 2.], [1., 3.], [1., 4.]])  # 4 eqns, 2 unknowns\n",
    "b = np.array([6., 5., 7., 10.])\n",
    "x_pinv = np.linalg.pinv(A) @ b\n",
    "x_lstsq, *_ = np.linalg.lstsq(A, b, rcond=None)\n",
    "print('pinv :', x_pinv.round(4), ' residual',\n",
    "      round(float(np.linalg.norm(A @ x_pinv - b)), 4))\n",
    "print('lstsq:', x_lstsq.round(4), ' residual',\n",
    "      round(float(np.linalg.norm(A @ x_lstsq - b)), 4))\n",
    "print('cond(A)     =', round(float(np.linalg.cond(A)), 3))\n",
    "print('cond(A^T A) =', round(float(np.linalg.cond(A.T @ A)), 3), '= cond(A)^2')"
   ],
   "execution_count": 6,
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "pinv : [3.5 1.4]  residual 2.0494\nlstsq: [3.5 1.4]  residual 2.0494\ncond(A)     = 7.469\ncond(A^T A) = 55.782 = cond(A)^2\n"
     ]
    }
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-svd-low-rank-md-7 fw=all -->\n",
    "\n",
    "### The Condition Number\n",
    "\n",
    "The *condition number* is the single SVD-derived scalar that predicts numerical\n",
    "pain. For an invertible (or full-rank) matrix it is the ratio of the largest to\n",
    "the smallest nonzero singular value,\n",
    "\n",
    "$$\n",
    "\\kappa(\\mathbf{A}) = \\frac{\\sigma_1}{\\sigma_r} .\n",
    "$$\n",
    "\n",
    "It measures how much a linear solve can amplify input error.\n",
    "\n",
    "**Proposition (error amplification, with tightness).** *Let $\\mathbf{A}$ be square\n",
    "invertible, $\\mathbf{A}\\mathbf{x}=\\mathbf{b}$, and let a perturbation\n",
    "$\\delta\\mathbf{b}$ change the solution to $\\mathbf{x}+\\delta\\mathbf{x}$. Then*\n",
    "\n",
    "$$\n",
    "\\frac{\\|\\delta\\mathbf{x}\\|}{\\|\\mathbf{x}\\|}\n",
    "   \\le \\kappa(\\mathbf{A})\\,\\frac{\\|\\delta\\mathbf{b}\\|}{\\|\\mathbf{b}\\|},\n",
    "$$\n",
    "\n",
    "*and the bound is attained for suitable $\\mathbf{b},\\delta\\mathbf{b}$.*\n",
    "\n",
    "**Proof.** From $\\mathbf{A}\\,\\delta\\mathbf{x}=\\delta\\mathbf{b}$ we get\n",
    "$\\delta\\mathbf{x}=\\mathbf{A}^{-1}\\delta\\mathbf{b}$, and the largest singular value\n",
    "of $\\mathbf{A}^{-1}$ is $1/\\sigma_n$, so by the equation\n",
    "$\\|\\delta\\mathbf{x}\\|\\le\\sigma_n^{-1}\\|\\delta\\mathbf{b}\\|$. From\n",
    "$\\mathbf{b}=\\mathbf{A}\\mathbf{x}$ we get $\\|\\mathbf{b}\\|\\le\\sigma_1\\|\\mathbf{x}\\|$.\n",
    "Multiplying the two inequalities,\n",
    "\n",
    "$$\n",
    "\\frac{\\|\\delta\\mathbf{x}\\|}{\\|\\mathbf{x}\\|}\n",
    "   \\le \\frac{\\sigma_1}{\\sigma_n}\\,\\frac{\\|\\delta\\mathbf{b}\\|}{\\|\\mathbf{b}\\|}\n",
    "   = \\kappa(\\mathbf{A})\\,\\frac{\\|\\delta\\mathbf{b}\\|}{\\|\\mathbf{b}\\|} .\n",
    "$$\n",
    "\n",
    "It is *tight*: align the signal with the direction the inverse amplifies least by\n",
    "taking $\\mathbf{b}=\\mathbf{u}_1$, so that $\\mathbf{x}=\\sigma_1^{-1}\\mathbf{v}_1$ is\n",
    "the *smallest*-norm solution any unit right-hand side can produce (a small\n",
    "$\\|\\mathbf{x}\\|$ in the denominator is exactly what maximizes the relative-error\n",
    "ratio), and align the error with the most-amplified direction\n",
    "$\\delta\\mathbf{b}=\\mathbf{u}_n$ (so $\\delta\\mathbf{x}=\\sigma_n^{-1}\\mathbf{v}_n$ is\n",
    "as large as a perturbation of its size can get);\n",
    "then both inequalities are equalities. $\\blacksquare$\n",
    "\n",
    "Three consequences follow. An orthogonal matrix has all singular\n",
    "values equal to $1$, so $\\kappa=1$: rotations and reflections are perfectly\n",
    "conditioned, which is *why* the SVD's $\\mathbf{U},\\mathbf{V}$ never amplify error\n",
    "and why orthogonal initialization is attractive. At the other extreme, forming the\n",
    "normal-equations matrix squares the conditioning, $\\kappa(\\mathbf{A}^\\top\\mathbf{A})=\\kappa(\\mathbf{A})^2$\n",
    "(its singular values are the $\\sigma_i^2$), which is the quantitative reason to\n",
    "prefer an SVD/QR least-squares solve over the normal equations\n",
    "(that section).\n",
    "\n",
    "Geometrically, a large $\\kappa$ means very *elongated* level sets. For a quadratic\n",
    "bowl $f(\\mathbf{x})=\\tfrac12\\mathbf{x}^\\top\\mathbf{M}\\mathbf{x}$ with $\\mathbf{M}$\n",
    "symmetric positive definite, the contours are ellipses whose axis ratio is exactly\n",
    "$\\sqrt{\\kappa(\\mathbf{M})}=\\sqrt{\\lambda_{\\max}/\\lambda_{\\min}}$; when\n",
    "$\\kappa\\gg1$ the bowl is a narrow valley, and gradient descent zig-zags across\n",
    "the steep walls while crawling along the flat floor. (For a least-squares\n",
    "objective $\\tfrac12\\|\\mathbf{A}\\mathbf{x}-\\mathbf{b}\\|^2$ the bowl's matrix\n",
    "is $\\mathbf{M}=\\mathbf{A}^\\top\\mathbf{A}$, so its geometric axis ratio is\n",
    "$\\kappa(\\mathbf A)$ while the Hessian condition number that controls gradient\n",
    "descent is $\\kappa(\\mathbf A)^2$.) This is the same picture that ended\n",
    "the Rayleigh discussion in that section, and it is no coincidence:\n",
    "with the best fixed step size, gradient descent's error contracts like\n",
    "$(\\kappa-1)/(\\kappa+1)$ per step on such a bowl: *one number, two consequences*,\n",
    "error amplification in a solve and slow convergence in optimization. We make this\n",
    "precise when we analyze gradient descent in that section\n",
    "and study numerical conditioning in that section.\n",
    "the figure contrasts a well-conditioned bowl ($\\kappa\\approx1$,\n",
    "near-circular contours, gradient descent heads almost straight to the minimum) with\n",
    "an ill-conditioned one ($\\kappa\\gg1$, elongated contours, a zig-zag trajectory). In\n",
    "the narrow valley, the step size is throttled by the steep direction (to stay\n",
    "stable) while the flat direction needs many such small steps to make progress;\n",
    "the zig-zag is the visible cost of a large condition number.\n",
    "\n",
    "![Gradient descent on a well-conditioned quadratic bowl ($\\kappa\\approx1$, near-circular contours, an almost straight path) versus an ill-conditioned one ($\\kappa\\gg1$, elongated contours, a zig-zag trajectory across the narrow valley).](../img/mdl-la-condition.svg)\n",
    "\n",
    "## The SVD in Modern Deep Learning\n",
    "\n",
    "Low-rank structure is everywhere in\n",
    "contemporary models, and Eckart--Young is the reason it can be exploited cheaply.\n",
    "\n",
    "**Low-rank adapters (LoRA).** Fine-tuning a large pretrained model by updating\n",
    "every weight is expensive. LoRA [@Hu.Shen.Wallis.ea.2021] freezes a\n",
    "pretrained weight $\\mathbf{W}\\in\\mathbb{R}^{m\\times n}$ and learns only a low-rank\n",
    "correction $\\Delta\\mathbf{W}=\\mathbf{B}\\mathbf{A}$ with\n",
    "$\\mathbf{B}\\in\\mathbb{R}^{m\\times r}$, $\\mathbf{A}\\in\\mathbb{R}^{r\\times n}$, and\n",
    "$r\\ll\\min(m,n)$. This trades $mn$ trainable parameters for $r(m+n)$, a ratio of\n",
    "\n",
    "$$\n",
    "\\frac{r(m+n)}{mn} ,\n",
    "$$\n",
    "\n",
    "which for a $4096\\times4096$ layer at rank $r=8$ is $8\\cdot8192/4096^2\\approx0.39\\%$\n",
    "of the parameters. the figure shows the architecture: the frozen\n",
    "weight and the trainable low-rank bypass act on the same input, and their outputs\n",
    "add. The hypothesis that fine-tuning updates are *intrinsically low\n",
    "rank* is what makes this work.\n",
    "\n",
    "![The LoRA architecture for one layer. The pretrained $\\mathbf{W}$ (all $mn = 16.8$M parameters for $m=n=4096$) stays frozen; a parallel bypass first compresses the input to $r=8$ dimensions with $\\mathbf{A}$ and then expands back with $\\mathbf{B}$, so only $r(m+n) \\approx 65.5$K parameters (about $0.39\\%$) are trained. The two branches add: $\\mathbf{h} = \\mathbf{W}\\mathbf{x} + \\mathbf{B}\\mathbf{A}\\mathbf{x}$. After training, $\\mathbf{B}\\mathbf{A}$ can be merged into $\\mathbf{W}$, so deployment costs nothing extra.](../img/mdl-la-lora.svg)\n",
    "\n",
    "Note what Eckart--Young does and does not say here:\n",
    "LoRA *learns* $\\mathbf{B}$ and $\\mathbf{A}$ by gradient descent rather than\n",
    "truncating a known matrix, so the theorem promises nothing about the learned\n",
    "adapter itself. What it quantifies is the *ceiling*: if the update that full\n",
    "fine-tuning would have made has singular values $\\sigma_1\\ge\\sigma_2\\ge\\cdots$,\n",
    "then no rank-$r$ adapter can come closer to it than $\\sigma_{r+1}$ in spectral\n",
    "norm ($\\sum_{i>r}\\sigma_i^2$ in squared Frobenius norm), so the approach can only succeed when\n",
    "the true update's spectrum decays fast, which is exactly the empirical finding\n",
    "that motivated LoRA. PiSSA makes the link to the SVD literal: it *initializes*\n",
    "$\\mathbf{B}$ and $\\mathbf{A}$ from the truncated SVD of the pretrained weight\n",
    "$\\mathbf{W}$ itself, so that fine-tuning starts by adapting the principal\n",
    "components [@Meng.Wang.Zhang.2024].\n",
    "\n",
    "**Orthogonalized updates (Muon).** The polar decomposition\n",
    "the equation is the engine of Muon [@Jordan.Jin.Boza.ea.2024], an\n",
    "optimizer adopted in recent large-scale language-model training. Write the\n",
    "*momentum* matrix of a weight (the running average of gradients that optimizers\n",
    "keep per weight; see the optimization chapters)\n",
    "as $\\mathbf{M}=\\mathbf{U}\\boldsymbol{\\Sigma}\\mathbf{V}^\\top$. A standard momentum\n",
    "step moves along $\\mathbf{M}$, which is dominated by its few largest dyads; Muon\n",
    "instead steps along the *polar factor*\n",
    "$\\mathbf{Q}=\\mathbf{U}\\mathbf{V}^\\top$, the rotation part of the update with\n",
    "every stretch factor equalized to $1$, so rare-but-consistent gradient\n",
    "directions are not drowned out by the dominant ones. (Equivalently, $\\mathbf{U}\\mathbf{V}^\\top$\n",
    "is the steepest-descent direction when distances between weight matrices are\n",
    "measured in the spectral norm rather than the Euclidean one\n",
    "[@Bernstein.Newhouse.2024].) Computing an SVD of every weight at every step\n",
    "would be hopeless; instead Muon runs a handful of *Newton--Schulz* iterations\n",
    "$\\mathbf{X}\\leftarrow\\tfrac12(3\\mathbf{X}-\\mathbf{X}\\mathbf{X}^\\top\\mathbf{X})$\n",
    "starting from $\\mathbf{X}=\\mathbf{M}/\\|\\mathbf{M}\\|_F$. Each iteration leaves\n",
    "$\\mathbf{U}$ and $\\mathbf{V}$ untouched and applies the odd polynomial\n",
    "$p(\\sigma)=\\tfrac12(3\\sigma-\\sigma^3)$ to every singular value, driving each\n",
    "nonzero $\\sigma_i$ toward the fixed point $1$. The Frobenius normalization is\n",
    "what makes that convergence guaranteed rather than lucky: it puts every singular\n",
    "value into $(0, 1]$, comfortably inside $(0, \\sqrt3\\,)$, the cubic's basin of\n",
    "attraction of $1$. Start a singular value beyond $\\sqrt3$ and the cubic maps it\n",
    "negative, after which the iterates either stick near the spurious fixed point\n",
    "$-1$ or diverge; either way the polar factor is lost. So the iterates converge to\n",
    "$\\mathbf{U}\\mathbf{V}^\\top$ using nothing but matrix multiplications (in practice\n",
    "Muon tunes the polynomial's coefficients so that about five iterations suffice).\n",
    "It is the same GPU-friendly bargain that spectral normalization strikes with\n",
    "power iteration below.\n",
    "\n",
    "The whole mechanism fits in a few lines.\n",
    "Frobenius-normalize a random \"momentum\" matrix, iterate the cubic, and every\n",
    "singular value marches to $1$ while the singular *vectors* stay put; the\n",
    "iterate lands on the polar factor $\\mathbf{U}\\mathbf{V}^\\top$ computed from a\n",
    "direct SVD. With the plain $\\tfrac12(3\\sigma-\\sigma^3)$ coefficients the\n",
    "smallest singular value (which starts closest to $0$, where $p$ only multiplies\n",
    "it by $\\tfrac32$) is the straggler, needing about ten iterations here; Muon's\n",
    "tuned polynomials accelerate exactly that small-$\\sigma$ regime."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "mdl-svd-low-rank-the-svd-in-modern-deep-learning-1"
   },
   "source": [
    "rng = np.random.default_rng(3)\n",
    "M = rng.standard_normal((6, 4))          # a stand-in momentum matrix\n",
    "X = M / np.linalg.norm(M)                # Frobenius normalize: all sigma <= 1\n",
    "for k in range(1, 11):\n",
    "    X = 1.5 * X - 0.5 * X @ X.T @ X      # applies p(sigma) = (3 sigma - sigma^3)/2\n",
    "    if k % 2 == 0:\n",
    "        print(f'iter {k:2d}: sigma =', np.linalg.svd(X, compute_uv=False).round(4))\n",
    "U, s, Vt = np.linalg.svd(M, full_matrices=False)\n",
    "print('distance to polar factor U V^T:',\n",
    "      round(float(np.linalg.norm(X - U @ Vt)), 6))"
   ],
   "execution_count": 7,
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "iter  2: sigma = [0.9982 0.8554 0.4648 0.1499]\niter  4: sigma = [1.     0.9987 0.8351 0.3291]\niter  6: sigma = [1.     1.     0.9978 0.6599]\niter  8: sigma = [1.     1.     1.     0.9663]\niter 10: sigma = [1. 1. 1. 1.]\ndistance to polar factor U V^T: 4e-06\n"
     ]
    }
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-svd-low-rank-md-8 fw=all -->\n",
    "\n",
    "**Spectral normalization.** Constraining the largest singular value of each weight\n",
    "matrix to $1$ caps the per-layer *Lipschitz constant*, the largest factor by\n",
    "which the layer can stretch any input difference\n",
    "(the equation: $\\sigma_1=\\|\\mathbf{W}\\|_2$).\n",
    "This stabilizes training, originally for GAN\n",
    "discriminators [@Miyato.Kataoka.Koyama.ea.2018]. $\\sigma_1$ is\n",
    "estimated by **power iteration on $\\mathbf{W}^\\top\\mathbf{W}$**, the very algorithm\n",
    "we analyzed in that section, since\n",
    "$\\sigma_1=\\sqrt{\\lambda_1(\\mathbf{W}^\\top\\mathbf{W})}$. A couple of iterations per\n",
    "training step suffice.\n",
    "\n",
    "**Weight and attention spectra.** Plotting the singular values of a trained layer\n",
    "(the spectrum panel of the figure, applied to $\\mathbf{W}$\n",
    "instead of an image) reveals its *effective rank* via the energy ratio\n",
    "the equation and exposes heavy-tailed spectra that correlate with\n",
    "generalization [@Martin.Mahoney.2021], heavy-tailed precisely against the\n",
    "Marchenko--Pastur noise baseline of that section.\n",
    "\n",
    "Attention matrices are often empirically near-low-rank, which\n",
    "motivates linear-attention approximations that replace the full softmax attention\n",
    "with a low-rank surrogate. The most prominent recent instance is low-rank\n",
    "key--value compression: DeepSeek-V2's multi-head latent attention stores a\n",
    "single low-rank latent in place of the full per-head key--value cache (the cached\n",
    "per-token keys and values that dominate inference memory) and\n",
    "expands it on the fly [@DeepSeek-AI.2024].\n",
    "\n",
    "The cell below makes the weight-side\n",
    "diagnostic concrete: it builds a synthetic weight matrix with a fast-decaying\n",
    "spectrum and reports the rank needed for 95% spectral energy and the parameter\n",
    "saving a LoRA of that rank would give; here rank 18 of 256, about 10.5% of the\n",
    "full parameter count. How small that rank is depends entirely on how fast the\n",
    "spectrum decays (compare the 0.39% headline above, which assumed $r=8$ suffices),\n",
    "which is why one measures the spectrum instead of assuming a rank."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "svd-weight-spectrum"
   },
   "source": [
    "rng = np.random.default_rng(2)\n",
    "m, n = 256, 512\n",
    "sigma = np.exp(-np.arange(min(m, n)) / 12.0)                 # fast-decaying spectrum\n",
    "Uw, _ = np.linalg.qr(rng.standard_normal((m, min(m, n))))\n",
    "Vw, _ = np.linalg.qr(rng.standard_normal((n, min(m, n))))\n",
    "s = np.linalg.svd((Uw * sigma) @ Vw.T, compute_uv=False)\n",
    "r95 = int(np.searchsorted(np.cumsum(s ** 2) / (s ** 2).sum(), 0.95) + 1)\n",
    "print(f'{m} x {n} matrix; rank for 95% energy: {r95}; '\n",
    "      f'LoRA params {r95 * (m + n) / (m * n):.1%} of full')"
   ],
   "execution_count": 8,
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "256 x 512 matrix; rank for 95% energy: 18; LoRA params 10.5% of full\n"
     ]
    }
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-svd-low-rank-md-9 fw=all -->\n",
    "\n",
    "**Scaling up.** For matrices too large to factor fully, *randomized SVD*\n",
    "[@Halko.Martinsson.Tropp.2011] computes an accurate rank-$k$ truncation at a\n",
    "fraction of the cost. The mechanism is a *range finder*: draw a random\n",
    "$n\\times(k{+}p)$ Gaussian matrix $\\boldsymbol\\Omega$ (a small oversampling $p$\n",
    "guards against unlucky draws), form $\\mathbf{Y}=\\mathbf{A}\\boldsymbol\\Omega$ (a\n",
    "handful of matrix--vector products that sketch the column space), and orthonormalize it,\n",
    "$\\mathbf{Q}=\\operatorname{qr}(\\mathbf{Y})$. Because $\\mathbf{A}$'s energy is\n",
    "concentrated in its top singular directions, the $k{+}p$ columns of $\\mathbf{Q}$\n",
    "capture them with high probability, so one then factors only the *tiny* matrix\n",
    "$\\mathbf{Q}^\\top\\mathbf{A}$ and maps its left singular vectors back through\n",
    "$\\mathbf{Q}$. The cost drops from the $O(mn\\min(m,n))$ of a full SVD to essentially\n",
    "a few passes over $\\mathbf{A}$; this is the standard tool when only the leading singular\n",
    "triples are needed. A related family, *CUR* decompositions,\n",
    "approximates $\\mathbf{A}$ with actual rows and columns of $\\mathbf{A}$, trading\n",
    "accuracy for interpretability.\n",
    "\n",
    "Let us verify the two facts the whole section rests on: that the SVD reconstructs\n",
    "$\\mathbf{A}$, and that $\\sigma_i^2$ are the eigenvalues of $\\mathbf{A}^\\top\\mathbf{A}$."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "svd-verify"
   },
   "source": [
    "A = np.array([[3., 1.], [1., 3.], [0., 2.]])   # a 3x2 matrix\n",
    "U, s, Vt = np.linalg.svd(A, full_matrices=False)\n",
    "recon = (U * s) @ Vt\n",
    "eig = np.sort(np.linalg.eigvalsh(A.T @ A))[::-1]\n",
    "print('reconstruction error:', round(float(np.linalg.norm(recon - A)), 12))\n",
    "print('sigma^2          :', (s ** 2).round(6))\n",
    "print('eig(A^T A) sorted:', eig.round(6))"
   ],
   "execution_count": 9,
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "reconstruction error: 0.0\nsigma^2          : [18.324555  5.675445]\neig(A^T A) sorted: [18.324555  5.675445]\n"
     ]
    }
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-svd-low-rank-md-10 fw=all -->\n",
    "\n",
    "The reconstruction error is at the level of floating-point round-off, and the\n",
    "squared singular values match the eigenvalues of $\\mathbf{A}^\\top\\mathbf{A}$:\n",
    "exactly the construction of that section.\n",
    "\n",
    "## Summary\n",
    "\n",
    "* Every matrix factors as $\\mathbf{A}=\\mathbf{U}\\boldsymbol{\\Sigma}\\mathbf{V}^\\top$\n",
    "  (rotate--scale--rotate; equivalently the polar form $\\mathbf{A}=\\mathbf{Q}\\mathbf{P}$),\n",
    "  with $\\mathbf{A}\\mathbf{v}_i=\\sigma_i\\mathbf{u}_i$ and the dyadic sum\n",
    "  $\\mathbf{A}=\\sum_i\\sigma_i\\mathbf{u}_i\\mathbf{v}_i^\\top$.\n",
    "* The SVD is the spectral theorem applied to the symmetric PSD matrix\n",
    "  $\\mathbf{A}^\\top\\mathbf{A}$, with $\\sigma_i=\\sqrt{\\lambda_i(\\mathbf{A}^\\top\\mathbf{A})}$;\n",
    "  it therefore *never fails*, even for the defective shear that had no\n",
    "  eigenbasis. Equivalently $\\sigma_1=\\max_{\\|\\mathbf{x}\\|=1}\\|\\mathbf{A}\\mathbf{x}\\|=\\|\\mathbf{A}\\|_2$.\n",
    "* Rank, range, and the four fundamental subspaces read off the spectrum;\n",
    "  numerical rank thresholds the $\\sigma_i$ at\n",
    "  $\\sim\\sigma_1\\max(m,n)\\,\\epsilon_{\\text{mach}}$.\n",
    "* **Eckart--Young--Mirsky:** the top-$k$ truncation $\\mathbf{A}_k$ is the optimal\n",
    "  rank-$k$ approximation, with $\\|\\mathbf{A}-\\mathbf{A}_k\\|_2=\\sigma_{k+1}$ and\n",
    "  $\\|\\mathbf{A}-\\mathbf{A}_k\\|_F^2=\\sum_{i>k}\\sigma_i^2$. The energy ratio is the\n",
    "  \"how much did we keep\" dial.\n",
    "* **PCA is Eckart--Young on centered data:** principal directions are the right\n",
    "  singular vectors $\\mathbf{v}_i$, explained variance is $\\sigma_i^2/n$.\n",
    "* The pseudoinverse $\\mathbf{A}^{+}=\\mathbf{V}\\boldsymbol{\\Sigma}^{+}\\mathbf{U}^\\top$\n",
    "  gives the minimum-norm least-squares solution; the condition number\n",
    "  $\\kappa=\\sigma_1/\\sigma_r$ predicts both error amplification and\n",
    "  gradient-descent speed, and $\\kappa(\\mathbf{A}^\\top\\mathbf{A})=\\kappa(\\mathbf{A})^2$\n",
    "  is why the normal equations are worse.\n",
    "* SVD powers PCA, LoRA (rank-$r$ updates at $r(m+n)$ parameters), Muon (steps\n",
    "  along the polar factor $\\mathbf{U}\\mathbf{V}^\\top$, computed by Newton--Schulz\n",
    "  iterations), spectral normalization (power iteration on\n",
    "  $\\mathbf{W}^\\top\\mathbf{W}$), and weight / attention spectral analysis.\n",
    "\n",
    "## Exercises\n",
    "\n",
    "1. Compute the SVD of $\\operatorname{diag}(3,1)$ by inspection. Then show the\n",
    "   scaled rotation $\\begin{bmatrix}0&-2\\\\1&0\\end{bmatrix}$ has singular values\n",
    "   $\\{2,1\\}$ even though its eigenvalue magnitudes are both $\\sqrt2$, i.e.\n",
    "   $\\sigma\\neq|\\lambda|$ for non-symmetric matrices.\n",
    "2. Prove that $\\|\\mathbf{A}\\|_F^2=\\sum_i\\sigma_i^2$. (*Hint:* first show\n",
    "   $\\|\\mathbf{A}\\|_F^2=\\operatorname{tr}(\\mathbf{A}^\\top\\mathbf{A})$ directly from\n",
    "   the definitions, then evaluate the trace as the sum of the eigenvalues of\n",
    "   $\\mathbf{A}^\\top\\mathbf{A}$ using the equation.) Conclude that\n",
    "   multiplying by orthogonal matrices on either side never changes the Frobenius\n",
    "   norm.\n",
    "3. Generate a rank-5 signal $\\mathbf{B}=\\mathbf{X}\\mathbf{Y}^\\top$ with\n",
    "   $\\mathbf{X},\\mathbf{Y}\\in\\mathbb{R}^{200\\times5}$ standard Gaussian, and add\n",
    "   noise $\\mathbf{N}$ with i.i.d. $\\mathcal{N}(0,\\sigma_{\\text{noise}}^2)$ entries,\n",
    "   $\\sigma_{\\text{noise}}=0.05$. Plot the singular values of\n",
    "   $\\mathbf{B}+\\mathbf{N}$ on a log scale and watch the spectrum split into five\n",
    "   signal values and a noise floor near\n",
    "   $(\\sqrt{m}+\\sqrt{n})\\,\\sigma_{\\text{noise}}$. Verify that hard-thresholding at\n",
    "   the Gavish--Donoho cutoff\n",
    "   $\\tfrac{4}{\\sqrt3}\\sqrt{n}\\,\\sigma_{\\text{noise}}$ recovers the true rank, and\n",
    "   check numerically whether the rank-5 truncation of $\\mathbf{B}+\\mathbf{N}$ is\n",
    "   closer to $\\mathbf{B}$ (in Frobenius norm) than $\\mathbf{B}+\\mathbf{N}$ itself\n",
    "   is. What happens as you keep more components than 5?\n",
    "4. Show that the singular values of any orthogonal matrix are all $1$, hence\n",
    "   $\\kappa=1$. Then prove $\\kappa(\\mathbf{A}^\\top\\mathbf{A})=\\kappa(\\mathbf{A})^2$\n",
    "   and explain why this makes the normal equations numerically inferior.\n",
    "5. Show $\\mathbf{A}$ and $\\mathbf{A}^\\top$ have the same nonzero singular values.\n",
    "   (*Hint:* relate $\\mathbf{A}^\\top\\mathbf{A}$ and $\\mathbf{A}\\mathbf{A}^\\top$ via\n",
    "   the equation.)\n",
    "6. Prove the polar decomposition $\\mathbf{A}=\\mathbf{Q}\\mathbf{P}$ has $\\mathbf{Q}$\n",
    "   orthogonal and $\\mathbf{P}\\succeq0$, and that $\\mathbf{P}$ is unique while\n",
    "   $\\mathbf{Q}$ is unique when $\\mathbf{A}$ is invertible.\n",
    "7. For a given weight matrix, find the LoRA rank achieving 95% spectral energy and\n",
    "   compute the resulting parameter saving relative to a full update (the\n",
    "   `#svd-weight-spectrum` cell). How does the answer change if the spectrum decays\n",
    "   more slowly?\n",
    "8. How much can an update $\\boldsymbol{\\Delta}$ (for instance a rank-$r$ LoRA\n",
    "   update $\\mathbf{B}\\mathbf{A}$) move the singular values of $\\mathbf{W}$?\n",
    "    1. Derive the *Weyl perturbation bound*\n",
    "       $|\\sigma_i(\\mathbf{W}+\\boldsymbol{\\Delta})-\\sigma_i(\\mathbf{W})|\\le\\sigma_1(\\boldsymbol{\\Delta})$\n",
    "       for every $i$ and *any* update $\\boldsymbol{\\Delta}$, by setting $j=1$ in\n",
    "       the Weyl inequality lemma of that section,\n",
    "       $\\sigma_{i+j-1}(\\mathbf{X}+\\mathbf{Y})\\le\\sigma_i(\\mathbf{X})+\\sigma_j(\\mathbf{Y})$,\n",
    "       applied in both directions.\n",
    "    1. Can a rank-*one* update move *all* the singular values? Take\n",
    "       $\\mathbf{W}=\\operatorname{diag}(1,2)$ and\n",
    "       $\\boldsymbol{\\Delta}=0.1\\,\\mathbf{1}\\mathbf{1}^\\top$, compute the singular\n",
    "       values of $\\mathbf{W}+\\boldsymbol{\\Delta}$ numerically, and check that\n",
    "       *both* moved, each by less than $\\sigma_1(\\boldsymbol{\\Delta})=0.2$, as the\n",
    "       bound demands.\n",
    "    1. Reconcile this with the fact that adding a rank-$r$ update changes\n",
    "       $\\operatorname{rank}\\mathbf{W}$ by at most $r$: a low-rank update can nudge\n",
    "       *every* singular value a little, but can create or destroy at most $r$\n",
    "       nonzero ones."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-svd-low-rank-md-11 fw=pytorch -->\n",
    "\n",
    "[Discussions](https://d2l.discourse.group/t/svd)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "name": "d2l-pytorch",
   "display_name": "d2l (PyTorch)",
   "language": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
