The SVD, Eckart–Young–Mirsky, polar decomposition and the pseudoinverse

Before this lesson: lesson 06 (the spectral theorem), lesson 07 ($A^\top A$ is PSD) and lesson 02 (the four subspaces). Time: about 95 minutes, plus the notebook.

---

Why this matters #

The spectral theorem only works for symmetric matrices. The SVD works for every matrix, of every shape. It gives the rank in floating point, orthonormal bases for all four subspaces, the best low-rank approximation, the pseudoinverse and the nearest orthogonal matrix. In ML it is behind low-rank compression, PCA (lesson 09), the spectral norm of a layer and the condition number $\sigma_{\max}/\sigma_{\min}$ (lesson 13). In science, the SVD of a space–time data matrix is called POD (proper orthogonal decomposition) in fluid dynamics and EOF analysis (empirical orthogonal functions) in climate science.

1. Predict first #

2. The picture #

Start with a diagonal matrix, $A = \operatorname{diag}(3, -2)$. It stretches the $x_1$-axis by 3, and it stretches the $x_2$-axis by 2 and flips it. The SVD separates the stretch from the flip: the singular values are $3$ and $2$ (always $\ge 0$), and the sign goes into an orthogonal factor.

Now take any $2 \times 2$ matrix $A$. It maps the unit circle to an ellipse. Call the semi-axes of the ellipse $\sigma_1 u_1$ and $\sigma_2 u_2$, with orthogonal unit vectors $u_1, u_2$ and $\sigma_1 \ge \sigma_2 \ge 0$. The unit vectors $v_1, v_2$ that land on these axes are also orthogonal. So an orthonormal input basis goes to an orthogonal output set:

Put the $v_i$ as columns of $V$ and the $u_i$ as columns of $U$. Then $AV = U\Sigma$, so $A = U\Sigma V^\top$. Read it from right to left:

1. $V^\top$ writes the input in coordinates along the $v_i$.
2. $\Sigma$ stretches coordinate $i$ by $\sigma_i$. If $A$ is not square, $\Sigma$ also changes the dimension from $n$ to $m$.
3. $U$ places the result along the $u_i$.

$U$ and $V$ are orthogonal, so each is a rotation or a reflection. For square $A$ with $\det A < 0$, at least one of them must contain a reflection, because $\det A = \det U \cdot \det \Sigma \cdot \det V$ and $\det \Sigma \ge 0$.

Interactive figure: loading…

Why such a pair always exists #

The claim "some orthonormal pair $v_1, v_2$ has orthogonal images" is the heart of the SVD. Here is why it is true for every $2 \times 2$ matrix, step by step.

1. Turn a pair. Every orthonormal pair can be written as $v_1 = (\cos t, \sin t)$ and $v_2 = (-\sin t, \cos t)$ for an angle $t$, up to the sign of $v_2$.
2. Measure the angle of the images. Let $f(t) = (Av_1) \cdot (Av_2) = v_1^\top A^\top A\, v_2$. The images are orthogonal exactly when $f(t) = 0$.
3. Turn by $90°$. Then $v_1$ becomes $v_2$, and $v_2$ becomes $-v_1$. So $f(t + 90°) = (Av_2) \cdot (-Av_1) = -f(t)$: the sign of $f$ flips.
4. So $f$ has a zero. $f$ is a continuous function of $t$. It changes sign between $t$ and $t + 90°$, so it is zero somewhere in between (intermediate value theorem).
5. What the zero means. $f(t) = 0$ says that $A^\top A v_1$ has no component along $v_2$. In the plane, that leaves only the direction $v_1$: $A^\top A v_1 = \lambda_1 v_1$. So $v_1$, and in the same way $v_2$, is an eigenvector of the symmetric matrix $A^\top A$.
6. Largest and smallest stretch. Since $\frac{d}{dt} v_1 = v_2$, we get $\frac{d}{dt}\|Av_1\|^2 = 2\, v_2^\top A^\top A\, v_1 = 2f(t)$. So $f = 0$ exactly where $\|Av_1\|$ is largest or smallest. These extreme lengths are $\sigma_1$ and $\sigma_2$.

In $n$ dimensions the same search is the spectral theorem for $A^\top A$ (lesson 06). The proof in section 3 uses it directly.

Interactive figure: loading…

The same factorization is a sum of rank-1 layers, sorted by importance:

3. Definitions and main results #

The SVD #

Let $A \in \mathbb{R}^{m \times n}$ have rank $r$, and let $p = \min(m, n)$.

• Full SVD: $A = U\Sigma V^\top$ with $U \in \mathbb{R}^{m\times m}$ and $V \in \mathbb{R}^{n \times n}$ orthogonal, and $\Sigma \in \mathbb{R}^{m \times n}$ zero except for the diagonal $\sigma_1 \ge \dots \ge \sigma_p \ge 0$. Exactly $r$ of the $\sigma_i$ are positive.
• Thin (reduced) SVD (NumPy's full_matrices=False): keep $p$ columns of $U$ and $V$ and the $p \times p$ block of $\Sigma$.
• Compact (rank-$r$) SVD: keep only $r$ columns: $A = U_r \Sigma_r V_r^\top$ with $U_r \in \mathbb{R}^{m \times r}$, $\Sigma_r \in \mathbb{R}^{r\times r}$, $V_r \in \mathbb{R}^{n \times r}$.
• The singular values $\sigma_i$ are unique. The left singular vectors $u_i$ and the right singular vectors $v_i$ are unique only up to sign (for distinct nonzero $\sigma_i$), and a sign flip of $u_i$ needs the same flip of $v_i$.
• Column by column: $Av_i = \sigma_i u_i$ and $A^\top u_i = \sigma_i v_i$.
• The four subspaces (lesson 02) get orthonormal bases: $u_1, \dots, u_r$ span $C(A)$; $u_{r+1}, \dots, u_m$ span $N(A^\top)$; $v_1, \dots, v_r$ span $C(A^\top)$; $v_{r+1}, \dots, v_n$ span $N(A)$.
• Norms: $\|A\|_2 = \sigma_1$ and $\|A\|_F = \sqrt{\sigma_1^2 + \dots + \sigma_r^2}$, because orthogonal factors do not change these norms.

The three shapes for one example. Take $m = 4$, $n = 3$ and rank $r = 2$, so $p = 3$ and $\sigma_3 = 0$.

Form	$U$	$\Sigma$	$V$	What is dropped
full	$4 \times 4$	$4 \times 3$	$3 \times 3$	nothing
thin	$4 \times 3$	$3 \times 3$, diagonal $(\sigma_1, \sigma_2, 0)$	$3 \times 3$	$u_4$, which only meets a zero row of $\Sigma$
compact	$4 \times 2$	$2 \times 2$	$3 \times 2$	also $u_3$, $v_3$, which only meet $\sigma_3 = 0$

All three give the same product $A$. The dropped columns are multiplied by zeros of $\Sigma$, so they do not change the product. They are still useful: $u_3, u_4$ span $N(A^\top)$ and $v_3$ spans $N(A)$.

Theorem (SVD existence) #

Every real $m\times n$ matrix $A$ has an SVD.

Proof.

1. $A^\top A$ is symmetric and PSD: $x^\top A^\top A x = \|Ax\|^2 \ge 0$ (lesson 07).
2. Spectral theorem: $A^\top A = V \Lambda V^\top$ with $V$ orthogonal and $\lambda_1 \ge \dots \ge \lambda_n \ge 0$. Since $N(A^\top A) = N(A)$, exactly $r$ eigenvalues are positive.
3. For $i \le r$, set $\sigma_i = \sqrt{\lambda_i}$ and $u_i = A v_i / \sigma_i$.
4. The $u_i$ are orthonormal: $u_i^\top u_j = v_i^\top A^\top A v_j/(\sigma_i\sigma_j) = \lambda_j v_i^\top v_j/(\sigma_i \sigma_j) = \delta_{ij}$.
5. For $i > r$: $\|Av_i\|^2 = \lambda_i = 0$, so $Av_i = 0$.
6. Extend $u_1, \dots, u_r$ to an orthonormal basis $u_1, \dots, u_m$ of $\mathbb{R}^m$ (the new vectors span $N(A^\top)$).
7. Now $AV = U\Sigma$ column by column ($Av_i = \sigma_i u_i$ for $i \le r$, $Av_i = 0$ for $i > r$). Multiply by $V^\top$ on the right: $A = U\Sigma V^\top$. $\blacksquare$

Steps 5 and 6 on a small example. Take the rank-1 matrix $A = \begin{pmatrix} 1 & 1 \\ 1 & 1 \\ 0 & 0 \end{pmatrix}$ ($m = 3$, $n = 2$, $r = 1$).

• $A^\top A = \begin{pmatrix} 2 & 2 \\ 2 & 2 \end{pmatrix}$ has eigenvalues $4$ and $0$, with $v_1 = (1, 1)/\sqrt2$ and $v_2 = (1, -1)/\sqrt2$.
• Step 3: $\sigma_1 = 2$ and $u_1 = Av_1/2 = (\sqrt2, \sqrt2, 0)/2 = (1, 1, 0)/\sqrt2$.
• Step 5: $Av_2 = (0, 0, 0)$. There is no $\sigma_2 > 0$ to divide by, so step 3 cannot give $u_2$.
• Step 6: $N(A^\top)$ is the set of $y$ with $y_1 + y_2 = 0$, a plane. Pick any orthonormal basis of it, for example $u_2 = (1, -1, 0)/\sqrt2$ and $u_3 = (0, 0, 1)$. Both are orthogonal to $u_1$.
• Result: $U = [u_1\ u_2\ u_3]$, $\Sigma = \begin{pmatrix} 2 & 0 \\ 0 & 0 \\ 0 & 0 \end{pmatrix}$, $V = [v_1\ v_2]$, and $U\Sigma V^\top = 2\,u_1 v_1^\top = A$.

The choice in step 6 is not unique: any orthonormal basis of the plane $N(A^\top)$ works. It does not matter, because $u_2$ and $u_3$ meet only zeros of $\Sigma$.

Two consequences. First, $A^\top A = V\Sigma^\top\Sigma V^\top$ and $AA^\top = U\Sigma\Sigma^\top U^\top$ have the same nonzero eigenvalues $\sigma_i^2$. Second, step 3 is the safe way to pair $U$ with $V$. If you eigendecompose $AA^\top$ separately, each $u_i$ comes with an arbitrary sign (or an arbitrary basis of a repeated eigenspace), and then $U\Sigma V^\top \ne A$.

In floating point, do not compute the SVD from $A^\top A$. Forming it squares the condition number, and singular values below about $\sqrt{\varepsilon}\,\sigma_1 \approx 10^{-8}\sigma_1$ drown in rounding error (section 5). Library SVDs work on $A$ directly (lesson 14).

Theorem (Eckart–Young–Mirsky) #

Let $A_k = \sum_{i \le k} \sigma_i u_i v_i^\top$ (keep the $k$ largest terms, $k < r$). This is the truncated SVD. For every $B$ with $\operatorname{rank} B \le k$:

Why this is true, in words ($n = 3$, $k = 1$). Let $B$ have rank 1. We show that $B$ misses $A$ by at least $\sigma_2$ somewhere.

1. $B$ has rank 1, so by rank–nullity its null space $N(B)$ has dimension $3 - 1 = 2$: a plane.
2. $W = \operatorname{span}(v_1, v_2)$ is also a plane. On $W$, $A$ stretches every unit vector by at least $\sigma_2$.
3. Two planes through $0$ in $\mathbb{R}^3$ always share a line, because $2 + 2 > 3$. Take a unit vector $z$ on that line.
4. $Bz = 0$, so $(A - B)z = Az$, and $\|Az\| \ge \sigma_2$ because $z \in W$.

This is the same dimension-counting idea as Courant–Fischer in lesson 06, which has a figure of two planes that share a line. The general proof replaces "plane" by the right dimensions.

Proof (spectral norm).

1. $A - A_k = \sum_{i > k} \sigma_i u_i v_i^\top$ is itself an SVD, with largest singular value $\sigma_{k+1}$. So $\|A - A_k\|_2 = \sigma_{k+1}$.
2. Let $\operatorname{rank} B \le k$. By rank–nullity, $\dim N(B) \ge n - k$.
3. $W = \operatorname{span}(v_1, \dots, v_{k+1})$ has dimension $k+1$. Since $(n-k) + (k+1) > n$, the two subspaces share a unit vector $z \in N(B) \cap W$.
4. $(A - B)z = Az = \sum_{i \le k+1} \sigma_i (v_i^\top z) u_i$, so $\|(A-B)z\|^2 = \sum_{i \le k+1} \sigma_i^2 (v_i^\top z)^2 \ge \sigma_{k+1}^2 \sum_{i\le k+1}(v_i^\top z)^2 = \sigma_{k+1}^2$.
5. Hence $\|A - B\|_2 \ge \|(A-B)z\| \ge \sigma_{k+1}$. $\blacksquare$

Proof of the Frobenius case (uses Weyl's inequality for singular values)

This proof uses one result that this track does not prove in full: Weyl's inequality for singular values,

$$\sigma_{i+j-1}(X + Y) \le \sigma_i(X) + \sigma_j(Y) \qquad (i, j \ge 1,\ i + j - 1 \le \min(m, n)).$$

See Horn and Johnson, Topics in Matrix Analysis (Cambridge University Press, 1991), Theorem 3.3.16. The special case $j = 1$ gives $|\sigma_k(A + E) - \sigma_k(A)| \le \|E\|_2$, which lesson 17 proves.

1. Let $\operatorname{rank} B \le k$. Then $\sigma_{k+1}(B) = 0$.
2. Use Weyl with $X = A - B$, $Y = B$ and $j = k + 1$. Since $X + Y = A$: $\sigma_{i+k}(A) \le \sigma_i(A - B) + \sigma_{k+1}(B) = \sigma_i(A - B)$.
3. Square and sum over $i = 1, \dots, p - k$: $\sum_{i > k}\sigma_i(A)^2 \le \sum_{i \le p-k} \sigma_i(A - B)^2 \le \|A - B\|_F^2$.
4. The left side is $\|A - A_k\|_F^2$. $\blacksquare$

Reading the formulas. The relative Frobenius error is $\sqrt{\sum_{i>k}\sigma_i^2 / \sum_i \sigma_i^2}$: the root of the fraction of "energy" that you drop. For an error budget, pick the smallest $k$ that meets it.

Other norms, and when the best approximation is unique

• Mirsky extended the theorem to every unitarily invariant norm (a norm with $\|Q_1 A Q_2\| = \|A\|$ for orthogonal $Q_1, Q_2$). For example, in the nuclear norm $\|A\|_* = \sum_i \sigma_i$ the error of $A_k$ is $\sigma_{k+1} + \dots + \sigma_r$, and it is again optimal.
• The Frobenius minimizer is unique when $\sigma_k > \sigma_{k+1}$. The spectral-norm minimizer is usually not unique. A simple case: for $A = \operatorname{diag}(3, 1)$ and $k = 1$, both $\operatorname{diag}(3, 0)$ and $\operatorname{diag}(3.5, 0)$ have rank 1 and spectral error $\|A - B\|_2 = 1 = \sigma_2$.

Interactive figure: loading…

Polar decomposition #

For square $A$, insert $V^\top V = I$ into the SVD:

For $m > n$, use the thin SVD; then $Q$ is $m \times n$ with orthonormal columns.

• Meaning. $P$ is symmetric PSD: it stretches by $\sigma_i$ along the orthogonal axes $v_i$. $Q$ is orthogonal: it rotates or reflects. This is the matrix version of $z = e^{i\theta}|z|$ for a complex number.
• Uniqueness. $P$ is always unique (the PSD square root of $A^\top A$ is unique). If $A$ is invertible, $Q = AP^{-1}$ is unique too.
• Rotation or not. For invertible square $A$, $\det Q = \operatorname{sign}(\det A)$. So $Q$ is a rotation only when $\det A > 0$. Example: the swap $A = \begin{pmatrix}0&1\\1&0\end{pmatrix}$ has $P = I$ and $Q = A$, a reflection. In the same way, $U$ and $V$ in the SVD are orthogonal, not always rotations (lesson 04 explains the difference).
• Nearest orthogonal matrix. $Q = UV^\top$ minimizes $\|A - W\|_F$ over all orthogonal $W$.

Why $UV^\top$ is the nearest orthogonal matrix

1. For orthogonal $W$, $\|W\|_F^2 = n$. Expand: $\|A - W\|_F^2 = \|A\|_F^2 + n - 2\operatorname{tr}(W^\top A)$. So minimizing the distance means maximizing $\operatorname{tr}(W^\top A)$.
2. Insert the SVD: $\operatorname{tr}(W^\top A) = \operatorname{tr}(W^\top U\Sigma V^\top) = \operatorname{tr}(V^\top W^\top U\,\Sigma) = \operatorname{tr}(Z\Sigma)$, with $Z = V^\top W^\top U$ (the trace does not change under cyclic shifts).
3. $Z$ is a product of orthogonal matrices, so it is orthogonal and every entry has $|z_{ii}| \le 1$.
4. $\operatorname{tr}(Z\Sigma) = \sum_i z_{ii}\sigma_i \le \sum_i \sigma_i$, with equality at $Z = I$, that is, $W = UV^\top$. $\blacksquare$

Interactive figure: loading…

Pseudoinverse #

where $\Sigma^+$ ($n \times m$) has $1/\sigma_i$ for each nonzero $\sigma_i$ and zeros everywhere else. This is the (Moore–Penrose) pseudoinverse. $A^+$ undoes $A$ from the column space back to the row space, and it sends $N(A^\top)$ to zero:

• $AA^+ = U_rU_r^\top$ is the orthogonal projector onto $C(A)$;
• $A^+A = V_rV_r^\top$ is the orthogonal projector onto $C(A^\top)$.

Special cases: $A^+ = (A^\top A)^{-1}A^\top$ for full column rank, $A^+ = A^\top(AA^\top)^{-1}$ for full row rank, and $A^+ = A^{-1}$ for invertible $A$.

Theorem. $x^+ = A^+ b$ minimizes $\|Ax - b\|_2$, and among all minimizers it has the smallest $\|x\|_2$. It is the minimum-norm least-squares solution.

Proof. Let $c = U^\top b$ and $y = V^\top x$. Orthogonal maps keep 2-norms, so $\|x\| = \|y\|$.

1. $\|Ax - b\|^2 = \|U\Sigma V^\top x - b\|^2 = \|\Sigma y - c\|^2 = \sum_{i \le r} (\sigma_i y_i - c_i)^2 + \sum_{i > r} c_i^2$.
2. The second sum does not depend on $x$. The first sum is zero exactly when $y_i = c_i/\sigma_i$ for $i \le r$; the $y_i$ with $i > r$ are free. These are all the least-squares solutions.
3. $\|x\|^2 = \|y\|^2 = \sum_{i\le r} (c_i/\sigma_i)^2 + \sum_{i>r} y_i^2$ is smallest when the free $y_i$ are $0$.
4. That $y$ is $\Sigma^+ c$, so $x = V\Sigma^+U^\top b = A^+ b$. $\blacksquare$

Geometrically, $Ax^+$ is the projection of $b$ onto $C(A)$, and $x^+$ lies in the row space. Any other solution is $x^+ + w$ with $w \in N(A)$, and $w \perp x^+$, so $\|x^+ + w\|^2 = \|x^+\|^2 + \|w\|^2$.

Numerics. $A^+$ is not a continuous function of $A$: $\operatorname{diag}(1, \varepsilon)^+ = \operatorname{diag}(1, 1/\varepsilon)$ blows up as $\varepsilon \to 0$, but $\operatorname{diag}(1, 0)^+ = \operatorname{diag}(1, 0)$. So libraries treat $\sigma_i \le \text{rcond}\cdot\sigma_1$ as zero. The number of singular values above this cut-off is the numerical rank. This truncation is a form of regularization, like ridge.

Interactive figure: loading…

Answers to section 1. (1) No. $A^\top A = \begin{pmatrix}1&1\\1&2\end{pmatrix}$ has eigenvalues $(3 \pm \sqrt5)/2$, so $\sigma = (\sqrt5 \pm 1)/2 \approx 1.618,\ 0.618$. Eigenvalues describe directions that stay on their own line; singular values describe stretching. The products agree: $\sigma_1\sigma_2 = |\det A| = 1$. (2) By Eckart–Young, $\sqrt{1^2 + 0.1^2} = \sqrt{1.01} \approx 1.005$.

4. Worked example #

SVD. $A^\top A = \begin{pmatrix} 3 & 1 \\ 1 & 3\end{pmatrix}$ has eigenvalues $4$ and $2$, with eigenvectors $v_1 = (1,1)/\sqrt2$ and $v_2 = (1,-1)/\sqrt2$. So $\sigma_1 = 2$, $\sigma_2 = \sqrt2$.

• $u_1 = Av_1/\sigma_1 = (2, 2, 0)/(2\sqrt2) = (1, 1, 0)/\sqrt2$.
• $u_2 = Av_2/\sigma_2 = (0, 0, 2)/(\sqrt2\cdot\sqrt2) = (0, 0, 1)$.
• For the full SVD, $u_3 = (1, -1, 0)/\sqrt2$ spans $N(A^\top)$ (check: $A^\top u_3 = 0$).

Best rank-1 approximation. $A_1 = \sigma_1 u_1 v_1^\top = \begin{pmatrix} 1 & 1 \\ 1 & 1 \\ 0 & 0\end{pmatrix}$ and $A - A_1 = \begin{pmatrix} 0&0\\0&0\\1&-1\end{pmatrix}$. Both errors are $\sqrt2 = \sigma_2$, as the theorem says (only one singular value is discarded, so the two norms agree).

Polar. $Q = u_1v_1^\top + u_2v_2^\top = \begin{pmatrix} 1/2 & 1/2 \\ 1/2 & 1/2 \\ 1/\sqrt2 & -1/\sqrt2\end{pmatrix}$ (orthonormal columns) and $P = 2\,v_1v_1^\top + \sqrt2\, v_2v_2^\top = \begin{pmatrix} 1 + \tfrac{1}{\sqrt2} & 1 - \tfrac{1}{\sqrt2} \\ 1 - \tfrac{1}{\sqrt2} & 1 + \tfrac{1}{\sqrt2}\end{pmatrix}$. Check: $P^2 = \begin{pmatrix}3&1\\1&3\end{pmatrix} = A^\top A$ and $QP = A$.

Pseudoinverse and least squares. $A^+ = \tfrac12 v_1u_1^\top + \tfrac{1}{\sqrt2} v_2u_2^\top = \begin{pmatrix} 1/4 & 1/4 & 1/2 \\ 1/4 & 1/4 & -1/2\end{pmatrix}$, the same as $(A^\top A)^{-1}A^\top$ (full column rank). For $b = (2, 4, 1)$: $x^+ = A^+b = (2, 1)$, $Ax^+ = (3, 3, 1)$, and the residual $b - Ax^+ = (-1, 1, 0)$ is a multiple of $u_3$, so it is orthogonal to $C(A)$. ✓

Minimum norm. Now solve with the rank-1 matrix $A_1$. Its pseudoinverse is $A_1^+ = \tfrac12 v_1u_1^\top = \tfrac14\begin{pmatrix} 1&1&0\\1&1&0\end{pmatrix}$, and $A_1^+ b = (3/2, 3/2)$. Since $A_1x = (x_1 + x_2)(1, 1, 0)$ and the best fit of $(2, 4)$ by one number is $3$, every $x$ with $x_1 + x_2 = 3$ is a least-squares solution. The pseudoinverse picks the one in the row space $\operatorname{span}\{(1,1)\}$, which has the smallest norm.

5. In code #

The libraries check the hand computations. The notebook builds each of these from the SVD (and the SVD itself from eigh).

import numpy as np
from scipy.linalg import polar

A = np.array([[1., 1.], [1., 1.], [1., -1.]])
U, s, Vt = np.linalg.svd(A)                # full SVD: U (3, 3), s (2,), Vt (2, 2)
print(s)                                   # [2.         1.41421356]
print(np.linalg.svd(A, full_matrices=False).U.shape)    # thin U: (3, 2)

A1 = np.array([[1., 1.], [1., 1.], [0., 0.]])           # rank-1 answer of section 4
print(np.linalg.norm(A - A1, 2), np.linalg.norm(A - A1, "fro"))   # both 1.414...

Q, P = polar(A)                            # A = Q P
print(Q.round(4))                          # [[0.5 0.5] [0.5 0.5] [0.7071 -0.7071]]
print(P.round(4))                          # [[1.7071 0.2929] [0.2929 1.7071]]

b = np.array([2., 4., 1.])
print(np.linalg.pinv(A) @ b)               # [2. 1.]
print(np.linalg.pinv(A1) @ b)              # [1.5 1.5]: the minimum-norm solution
print(np.linalg.lstsq(A1, b, rcond=None)[0])   # the same

print(np.linalg.svd(np.array([[1., 1.], [0., 1.]]), compute_uv=False))  # [1.618 0.618]

# Why not the SVD from eig(AᵀA): a 50 x 2 matrix with σ = 1 and 1e-12
rng = np.random.default_rng(0)
Qa, _ = np.linalg.qr(rng.standard_normal((50, 2)))
Qb, _ = np.linalg.qr(rng.standard_normal((2, 2)))
B = Qa @ np.diag([1.0, 1e-12]) @ Qb.T
print(np.linalg.svd(B, compute_uv=False)[1])           # 1.0000084e-12: correct to about ε·σ₁
print(np.sqrt(abs(np.linalg.eigvalsh(B.T @ B)[0])))    # 1.077e-09: rounding noise

np.linalg.svd never forms $A^\top A$. It calls LAPACK's gesdd: Householder bidiagonalization of $A$, then a divide-and-conquer SVD of the bidiagonal matrix (lesson 14). The QR-iteration driver is gesvd; SciPy offers it as scipy.linalg.svd(..., lapack_driver="gesvd").

6. Where it shows up #

• Low-rank compression with an error budget. A weight matrix $W \in \mathbb{R}^{m\times n}$ stored as $U_k\Sigma_k V_k^\top$ needs $k(m + n + 1)$ numbers instead of $mn$, with Frobenius error exactly $\sqrt{\sum_{i>k}\sigma_i^2}$, and no matrix of rank $\le k$ does better. LoRA (lesson 15) also restricts an update to rank $\le r$.
• POD and EOFs in science. Stack snapshots of a flow or a climate field as columns. The left singular vectors are the spatial modes, $\sigma_i^2$ is the energy of mode $i$, and truncation gives the basis of a reduced-order model or a surrogate. PCA (lesson 09) is the same on centered data.
• Least squares and inverse problems. np.linalg.lstsq uses the SVD and returns the minimum-norm solution for collinear features. Truncated SVD and ridge damp the small $\sigma_i$ that amplify noise in ill-posed problems (deblurring, tomography).
• Orthogonalization. Orthogonal Procrustes (aligning molecules, point clouds or embeddings): the orthogonal $W$ that minimizes $\|WX - Y\|_F$ is $UV^\top$ from the SVD of $YX^\top$. The Muon optimizer (2024) replaces each momentum matrix by an approximation of its polar factor $UV^\top$, computed with a Newton–Schulz iteration.
• Spectral norm. $\sigma_1(W)$ is the Lipschitz constant of $x \mapsto Wx$. Spectral normalization divides $W$ by an estimate of $\sigma_1$ from power iteration.

7. Common mistakes #

• Getting $U$ and $V$ from two separate eigendecompositions. The signs (or the bases of repeated eigenspaces) do not match. Use $u_i = Av_i/\sigma_i$.
• Treating singular values as $|$eigenvalues$|$. They agree only for normal matrices ($A^\top A = AA^\top$, for example symmetric or orthogonal). The shear has eigenvalues $1, 1$ and singular values $1.618, 0.618$.
• "$U$ and $V$ are rotations." They are orthogonal. A reflection is allowed, and it is forced when $\det A < 0$.
• Shapes. The full $U$ is $m\times m$, the thin one $m \times p$, the compact one $m \times r$. NumPy returns $V^\top$ (Vt), not $V$.
• The wrong Eckart–Young error. The spectral error is the single number $\sigma_{k+1}$; the Frobenius error is the root of the sum of all discarded $\sigma_i^2$.
• Inverting tiny singular values. $1/\sigma$ for $\sigma \approx 10^{-14}$ multiplies noise by $10^{14}$. Truncate with rcond.
• Writing $A^+ = (A^\top A)^{-1}A^\top$ for every $A$. It needs full column rank, and even then it is numerically worse than the SVD.

8. Check yourself #

1. A symmetric matrix has eigenvalues $3$ and $-2$. What are its singular values, and how are $u_i$ and $v_i$ related?
   Answer

   $3$ and $2$. $u_1 = v_1$ for $\lambda = 3$ and $u_2 = -v_2$ for $\lambda = -2$: $\sigma_i = |\lambda_i|$, and the sign goes into $U$.
2. $A$ is $1000 \times 500$ with rank $50$. Give the shapes of the factors of the full SVD and of the rank-$r$ SVD, and the storage of the rank-50 form.
   Answer

   Full: $U$ $1000\times1000$, $\Sigma$ $1000\times500$, $V$ $500\times500$. Rank-$r$: $1000\times50$, $50\times50$, $500\times50$. Storage: $50(1000 + 500 + 1) = 75{,}050$ numbers instead of $500{,}000$.
3. The singular values are $6, 3, 2, 1$. What is the smallest $k$ with a relative Frobenius error of at most 20%? What is the spectral error then?
   Answer

   $\|A\|_F^2 = 50$. $k = 2$: $\sqrt{5/50} \approx 0.316$, too big. $k = 3$: $\sqrt{1/50} \approx 0.141$. So $k = 3$, with spectral error $\sigma_4 = 1$ (relative $1/6$).
4. Compute the pseudoinverse of a nonzero column vector $a \in \mathbb{R}^m$ and interpret $a^+ b$.
   Answer

   The SVD is $a = (a/\|a\|)\,\|a\|\,(1)$, so $a^+ = a^\top/\|a\|^2$. $a^+b = a^\top b/\|a\|^2$ is the coefficient of the projection of $b$ onto the line through $a$.
5. Find the polar decomposition of $A = \begin{pmatrix}0 & -2\\ 1 & 0\end{pmatrix}$. Is $Q$ a rotation?
   Answer

   $A^\top A = \operatorname{diag}(1, 4)$, so $P = \operatorname{diag}(1, 2)$ and $Q = AP^{-1} = \begin{pmatrix}0&-1\\1&0\end{pmatrix}$, a rotation by $90°$ (consistent with $\det A = 2 > 0$).
6. Why does the spectral-norm proof of Eckart–Young use $k+1$ right singular vectors, not $k$?
   Answer

   To force an intersection with $N(B)$: $\dim N(B) \ge n - k$ and $(n - k) + (k + 1) > n$. On $\operatorname{span}(v_1, \dots, v_{k+1})$, $A$ stretches every unit vector by at least $\sigma_{k+1}$.

Recall cardsAnswer in your head first, then reveal. Add the cards to your review deck to see them again just before you would forget.

1. State the SVD of an $m \times n$ matrix $A$.

   Show answer

2. How are the singular values of $A$ related to $A^\top A$?

   Show answer

3. You know $v_i$ and $\sigma_i > 0$. How do you get $u_i$, and why not from a separate eigendecomposition of $AA^\top$?

   Show answer

4. Why does every $2 \times 2$ matrix have an orthonormal pair $v_1, v_2$ with $Av_1 \perp Av_2$?

   Show answer

5. Eckart–Young in the spectral norm: what is the smallest $\|A - B\|_2$ over $\operatorname{rank} B \le k$?

   Show answer

6. Eckart–Young in the Frobenius norm: what is the smallest $\|A - B\|_F$ over $\operatorname{rank} B \le k$?

   Show answer

7. $A$ has singular values $4, 2, 1$. What are the best rank-1 errors in the spectral and the Frobenius norm?

   Show answer

8. How do you get the polar decomposition $A = QP$ from the SVD?

   Show answer

9. Write the pseudoinverse with the SVD. What is special about $A^+b$?

   Show answer

10. Why should you not compute the SVD from eig(AᵀA) in floating point?

    Show answer

Add these cards to my reviewOpen the review →

9. Interview angle #

• "Explain the SVD in one sentence." — Every matrix maps some orthonormal basis $v_i$ to orthogonal vectors $\sigma_i u_i$: rotate or reflect, stretch, rotate or reflect.
• "How do singular values relate to eigenvalues?" — $\sigma_i^2$ are the eigenvalues of $A^\top A$. For symmetric $A$, $\sigma_i = |\lambda_i|$. For square $A$ in general, $\max_i|\lambda_i| \le \sigma_1$ and $\prod_i|\lambda_i| = \prod_i \sigma_i = |\det A|$.
• "How would you compress a layer to a given error?" — Compute the SVD, choose the smallest $k$ with $\sqrt{\sum_{i>k}\sigma_i^2}$ within the budget, and store $k(m+n+1)$ numbers. Eckart–Young says this is optimal for that rank.
• "Why not compute the SVD with eig(AᵀA)?" — It squares the condition number, so small singular values drown in rounding error.

10. Practice #

The notebook for this lesson runs in Colab. It has four exercises, each with a check cell and a folded solution:

1. svd_via_eig: build the thin SVD from np.linalg.eigh(A.T @ A) and $u_i = Av_i/\sigma_i$, the existence proof as code; see where it fails for a tiny singular value.
2. truncated_svd and eckart_young_error on a generated image: compare the measured errors with $\sigma_{k+1}$ and $\sqrt{\sum_{i>k}\sigma_i^2}$, and find the smallest $k$ for a 5% budget.
3. polar: $Q = UV^\top$, $P = V\Sigma V^\top$; check that $Q$ is the nearest orthogonal matrix and use it to undo a rotation of a point cloud (orthogonal Procrustes).
4. pinv: the pseudoinverse with an rcond cut-off; check the minimum-norm property on a rank-deficient least-squares problem.

11. Watch and read #

6. Singular Value Decomposition (SVD)Gilbert Strang · MIT 18.065, lecture 6 · 54 min

7. Eckart-Young: The Closest Rank k Matrix to AGilbert Strang · MIT 18.065, lecture 7 · 47 min

Singular Value Decomposition (SVD): Mathematical OverviewSteve Brunton · SVD series · 13 min

SVD and Optimal TruncationSteve Brunton · SVD series · 12 min

Least Squares Regression and the SVDSteve Brunton · SVD series · 6 min

• Gilbert Strang, MIT 18.065 (2018), lecture 6: Singular value decomposition (textbook section I.8) and lecture 7: Eckart–Young, the closest rank-$k$ matrix to $A$ (textbook section I.9).
• Nick Higham, What is the singular value decomposition? (Eckart–Young, pseudoinverse and polar decomposition on one page) and What is the polar decomposition?
• Brunton and Kutz, Data-Driven Science and Engineering, chapter 1, §1.1–1.4 (SVD, matrix approximation, pseudo-inverse and regression).