USAAIO Lesson 16, from Week 6 on linear algebra, fully worked. It derives the SVD A=UΣVᵀ from AᵀA with no steps skipped, identifies the singular values as the non-negative roots of the eigenvalues of AᵀA, and proves the rotate-scale-rotate geometry on a running 2×3 example. It then computes low-rank approximation and Eckart-Young entry by entry, works a real compression on a decaying spectrum, re-derives PCA as the SVD of centered data, and builds the pseudoinverse A⁺=VΣ⁺Uᵀ by hand, matching it to np.linalg.pinv. Every snippet runs standalone, and every number came from real execution. The lesson runs to 63 slides.
Subject: Machine Learning · 106 slides · code lesson
Open the interactive version of this deck · Homework for this lesson
Title
USAAIO · Lesson 16 · Week 6 (Linear Algebra)
The crown jewel: every matrix factors as rotate · scale · rotate. We derive A = UΣVᵀ from AᵀA with no skipped steps, compute it by hand on one running matrix, then build low-rank compression, PCA, and the pseudoinverse on top of it.
Objectives
A = UΣVᵀ for any m×n matrix and name U, Σ, V and their shapesAᵀA, step by stepA vᵢ = σᵢ uᵢk approximation Aₖ = Σ σᵢ uᵢ vᵢᵀ and quote Eckart-Young: the error is the first dropped singular valueA⁺ = VΣ⁺Uᵀ, matching np.linalg.pinv and the normal equationsWarm-up
Discussion prompt
Before we open Lesson 16: Singular Value Decomposition: without looking back, what was the main idea of Modern Optimizers, and what could you do by the end of it that you could not do before?
Hint: One sentence for the idea, one for the skill. If the second one is blank, that is the part to revisit.
Answer:
SGD with momentum, RMSprop's per-parameter rates, Adam = momentum + RMSprop + bias correction, AdamW's decoupled weight decay, and learning-rate warmup. Implement momentum, RMSprop, and Adam from scratch and race them on an ill-conditioned bowl.
Section
Part 1 of 8
Concept
In Lesson 10 we diagonalized a matrix as A = QΛQ⁻¹ using eigenvectors. But that needs A square, and even then it can fail — a defective matrix has too few independent eigenvectors, and a non-square matrix has no eigenvalues at all.
Most data matrices are rectangular: m samples by n features, m ≠ n. We need a factorization that works for every shape and every rank. That factorization is the SVD.
Counterexample
Discussion prompt
Most data matrices are rectangular: m samples by n features, m ≠ n. We need a factorization that works for every shape and every rank. That factorization is the SVD.
That is stated as though it always holds. Do one of two things: produce a case where it fails, or say precisely what rules such a case out. "It just does" is not on the menu.
Hint: Hunt at the extremes first — zero, one, negative, empty, equal. If every extreme survives, the reason they survive is the proof.
Concept
Every m×n matrix A — square or not, symmetric or not, full-rank or not — factors into two orthogonal matrices and one diagonal matrix of non-negative numbers:
\[ A = U \Sigma V^\top, \qquad U^\top U = I,\quad V^\top V = I,\quad \Sigma = \operatorname{diag}(\sigma_1 \ge \sigma_2 \ge \dots \ge 0) \]
U's columns uᵢ are the left singular vectors, V's columns vᵢ are the right singular vectors, and the σᵢ on Σ's diagonal are the singular values, sorted large to small.
Analogy
Discussion prompt
Explain A = UΣVᵀ by analogy to something with no Machine Learning in it at all — a queue, a recipe, a map, a bank balance, whatever fits. Then say where your analogy breaks.
Hint: An analogy that never breaks is not an analogy, it is the same idea wearing a hat. Find the seam — that is the part that is actually new.
Answer:
Every m×n matrix A — square or not, symmetric or not, full-rank or not — factors into two orthogonal matrices and one diagonal matrix of non-negative numbers:
Concept
For an m×n matrix A, the factors have fixed shapes. In the full SVD, U is m×m, Σ is m×n (rectangular, zeros off the diagonal), and Vᵀ is n×n.
\[ \underset{m\times n}{A} \;=\; \underset{m\times m}{U}\;\underset{m\times n}{\Sigma}\;\underset{n\times n}{V^\top} \]
In the thin (economy) SVD — full_matrices=False in NumPy — we keep only the first r = \min(m,n) columns, so U is m×r, Σ is r×r, Vᵀ is r×n. Same product A, less storage. We'll use thin everywhere.
Explain it
Discussion prompt
Explain The three shapes to a student a year behind you. No notation, no jargon they have not met — and it still has to be true.
Hint: If your explanation needs a symbol they have never seen, you are describing the notation rather than the idea.
Answer:
For an m×n matrix A, the factors have fixed shapes. In the full SVD, U is m×m, Σ is m×n (rectangular, zeros off the diagonal), and Vᵀ is n×n.
Intuition
Applying A to a vector z is three clean moves read right-to-left: Vᵀ rotates z into the input axes, Σ stretches each axis by its σᵢ, then U rotates the result into the output space.
Orthogonal matrices are rigid rotations (or reflections): they never stretch, they only reorient. All the stretching a linear map does is packed into the diagonal Σ.
So the entire geometric content of any linear map is: pick the right input axes, scale, pick the right output axes. That is what the SVD reveals.
Picture it
Figure (svg): A unit circle on the left maps through A to a tilted ellipse on the right; the ellipse's long semi-axis is labelled sigma-1 and its short semi-axis sigma-2.
Discussion prompt
Read the picture before the words. What is this showing, and what is the one thing it is built to make obvious? Commit to an answer, then read on.
Hint: Name the parts, then say what changes between them — and if nothing changes, say what is being held still.
Answer:
Feed every unit vector (the unit circle) through A. Because Vᵀ and U only rotate, the shape you get out is an ellipse, and its semi-axis lengths are exactly the singular values σ₁ ≥ σ₂ ≥ ….
Intuition
Feed every unit vector (the unit circle) through A. Because Vᵀ and U only rotate, the shape you get out is an ellipse, and its semi-axis lengths are exactly the singular values σ₁ ≥ σ₂ ≥ ….
Figure (svg): A unit circle on the left maps through A to a tilted ellipse on the right; the ellipse's long semi-axis is labelled sigma-1 and its short semi-axis sigma-2.
Concept
We carry one matrix through the whole lesson so every slide compounds. It is 2×3 — deliberately non-square, so eigendecomposition can't touch it but the SVD can:
\[ A = \begin{bmatrix} 3 & 2 & 2 \\ 2 & 3 & -2 \end{bmatrix} \]
By the end we'll know its singular values, both sets of singular vectors, its best rank-1 summary, and its pseudoinverse — all by hand, all matching NumPy. Spoiler: σ = [5, 3].
Section
Part 2 of 8 — every step
Concept
We can't eigendecompose the rectangular A. But AᵀA is always square (n×n), always symmetric, and always positive semidefinite (Lesson 13). Symmetric PSD matrices have a full set of orthonormal eigenvectors and non-negative eigenvalues.
That is exactly the well-behaved object we know how to diagonalize. The plan: eigendecompose AᵀA, and read the SVD of A straight off it.
Ranking
Put in order
Put the moves of Substitute A = UΣVᵀ into AᵀA into the order they have to happen.
Why: These are the moves of the worked example in the order it makes them, and each one is set up by the one before it. Just substitution. Transpose of a product reverses order: (UΣVᵀ)ᵀ = V Σᵀ Uᵀ.
Worked example
Assume the factorization exists and see what AᵀA must be. Plug A = UΣVᵀ in and simplify, one move per line:
Write AᵀA with A = UΣVᵀ
Why: Just substitution. Transpose of a product reverses order: (UΣVᵀ)ᵀ = V Σᵀ Uᵀ.
\[ A^\top A = (U\Sigma V^\top)^\top (U\Sigma V^\top) = V\,\Sigma^\top U^\top\, U\,\Sigma\, V^\top \]
Cancel UᵀU = I
Why: U has orthonormal columns, so UᵀU is the identity — it vanishes from the middle.
\[ A^\top A = V\,\Sigma^\top \Sigma\, V^\top \]
ΣᵀΣ is diagonal with entries σᵢ²
Why: Σ is diagonal, so ΣᵀΣ is diagonal too, holding the squared singular values σ₁², σ₂², … on its diagonal.
\[ A^\top A = V\,\operatorname{diag}(\sigma_1^2,\sigma_2^2,\dots)\,V^\top \]
Notation
Annotate
From Substitute A = UΣVᵀ into AᵀA — read this one piece at a time. What is each part doing?
On: \( A^\top A = (U\Sigma V^\top)^\top (U\Sigma V^\top) = V\,\Sigma^\top U^\top\, U\,\Sigma\, V^\top \)
Concept
The last line AᵀA = V \operatorname{diag}(σᵢ²) Vᵀ is exactly the eigendecomposition of a symmetric matrix: orthogonal eigenvectors on the outside, eigenvalues on the diagonal inside.
\[ \lambda_i(A^\top A) = \sigma_i^2 \quad\Longrightarrow\quad \boxed{\;\sigma_i = \sqrt{\lambda_i(A^\top A)}\;} \]
So the right singular vectors V are the eigenvectors of AᵀA, and the singular values are the non-negative square roots of its eigenvalues. Because AᵀA is PSD, every λᵢ ≥ 0, so every σᵢ is real — the SVD always exists.
Concept
The same argument on the other side gives U. Compute AAᵀ = UΣVᵀ V ΣᵀUᵀ = U (ΣΣᵀ) Uᵀ, so the left singular vectors U are the eigenvectors of AAᵀ.
\[ V = \text{eigenvectors of } A^\top A, \qquad U = \text{eigenvectors of } A A^\top \]
AᵀA (n×n) and AAᵀ (m×m) are different sizes but share the same nonzero eigenvalues σᵢ² — the extra eigenvalues on the bigger one are just zeros. One spectrum, two viewpoints.
Intuition
It feels surprising that a 2×3 and a 3×3 matrix agree on eigenvalues. But A sends the v-frame to the u-frame with stretches σᵢ, so whichever side you measure the stretch from, you get the same σᵢ.
AᵀA measures the squared stretch in the input space; AAᵀ measures it in the output space. Same stretches, so the same σᵢ². The bigger matrix just has extra directions that get squashed — those show up as the zero eigenvalues.
Practical upshot: always eigendecompose the smaller of AᵀA and AAᵀ. For a tall data matrix (m ≫ n), that's the tiny n×n AᵀA.
Estimation
Predict first
Take our running A and confirm the derivation numerically: eigenvalues of AᵀA, then their square roots. This block runs on its own:
Commit before you compute: what does Build AᵀA and its eigenvalues — in code come out to? A rough magnitude and the right form is enough — the point is to have something concrete to be wrong about.
Correct: eig(AᵀA) = [25, 9, 0] → σ = [√25, √9] = [5, 3]
Why: A prediction you can defend turns the computation into a check rather than a leap of faith — and an answer that contradicts it is caught on the spot. AᵀA is 3×3 but rank 2, so its third eigenvalue is exactly 0 — that zero is why A only has TWO nonzero singular values.
Worked example
Take our running A and confirm the derivation numerically: eigenvalues of AᵀA, then their square roots. This block runs on its own:
import numpy as np
A = np.array([[3.,2.,2.],[2.,3.,-2.]])
AtA = A.T @ A
eigvals = np.sort(np.linalg.eigvalsh(AtA))[::-1]
print(eigvals.round(4)) # [25. 9. 0.]
print(np.sqrt(eigvals[:2])) # [5. 3.] = singular valueseig(AᵀA) = [25, 9, 0] → σ = [√25, √9] = [5, 3]
Why: AᵀA is 3×3 but rank 2, so its third eigenvalue is exactly 0 — that zero is why A only has TWO nonzero singular values. Verified by execution.
| quantity | value (verified) |
|---|---|
| A^T A (3×3) | [[13,12,2],[12,13,−2],[2,−2,8]] |
| eig(A^T A), sorted | [25, 9, 0] |
| σ = √(top-2 eig) | [5, 3] |
| rank(A) | 2 (one zero eigenvalue) |
Comparison
Comparison matrix
From Build AᵀA and its eigenvalues — in code: refill the value (verified) column from what you know. The rest of the table is as it appeared.
| quantity | value (verified) |
|---|---|
| A^T A (3×3) | [[13,12,2],[12,13,−2],[2,−2,8]] |
| eig(A^T A), sorted | [25, 9, 0] |
| σ = √(top-2 eig) | [5, 3] |
| rank(A) | 2 (one zero eigenvalue) |
Missing information
Discussion prompt
The 2×2 matrix AAᵀ should give the same nonzero eigenvalues {25, 9} — just without the extra zero. Compute both spectra side by side:
What do you need to know — or decide — before the first line can be written? List everything the problem has to hand you.
Hint: Anything you would have to invent to get started is a thing the problem must supply.
Answer:
By hand for a 2×2, eigenvalues satisfy λ = (trace ± √(trace² − 4det))/2 = (34 ± √(1156−900))/2 = (34 ± 16)/2 = {25, 9}. Matches the code exactly.
Worked example
The 2×2 matrix AAᵀ should give the same nonzero eigenvalues {25, 9} — just without the extra zero. Compute both spectra side by side:
import numpy as np
A = np.array([[3.,2.,2.],[2.,3.,-2.]])
print(np.sort(np.linalg.eigvalsh(A @ A.T))[::-1].round(4)) # [25. 9.]
print(np.sort(np.linalg.eigvalsh(A.T @ A))[::-1].round(4)) # [25. 9. 0.]AAᵀ = [[17,8],[8,17]]: trace 34, det 225
Why: By hand for a 2×2, eigenvalues satisfy λ = (trace ± √(trace² − 4det))/2 = (34 ± √(1156−900))/2 = (34 ± 16)/2 = {25, 9}. Matches the code exactly.
| matrix | shape | eigenvalues (verified) |
|---|---|---|
| A Aᵀ | 2×2 | [25, 9] |
| Aᵀ A | 3×3 | [25, 9, 0] |
| shared nonzero part | — | {25, 9} → σ = {5, 3} |
Trade off
Comparison matrix
From Cross-check with AAᵀ: every row here is a choice with a cost. Fill the shape column, then say which row you would actually pick and what you give up for it.
| matrix | shape | eigenvalues (verified) |
|---|---|---|
| A Aᵀ | 2×2 | [25, 9] |
| Aᵀ A | 3×3 | [25, 9, 0] |
| shared nonzero part | — | {25, 9} → σ = {5, 3} |
Anomaly
Predict first
A student writes this, and it looks reasonable:
Like eigendecomposition, the SVD only applies to square (or at least symmetric) matrices — a 2×3 has no decomposition.
It is wrong. Say what breaks — and say it before you turn the page.
Correct: This confuses SVD with EIGENdecomposition.
The SVD exists for every matrix, any shape, any rank — that universality is its whole point.
Why: This confuses SVD with EIGENdecomposition. A 2×3 matrix has no eigenvalues (Az = λz needs the output to have the same shape as z), so people wrongly conclude it has no factorization at all.
Trap
Like eigendecomposition, the SVD only applies to square (or at least symmetric) matrices — a 2×3 has no decomposition.
Refuse to decompose the 2×3 A
Why: This confuses SVD with EIGENdecomposition. A 2×3 matrix has no eigenvalues (Az = λz needs the output to have the same shape as z), so people wrongly conclude it has no factorization at all.
The SVD exists for every matrix, any shape, any rank — that universality is its whole point.
A(2×3) = U(2×2) Σ(2×3) Vᵀ(3×3), σ = [5, 3]
Why: We built σ from eig(AᵀA), which is always square/symmetric/PSD no matter what shape A is. SVD generalizes eigendecomposition to rectangular AND defective matrices — verified: reconstruction returns A exactly.
Break the constraint
Discussion prompt
The rule this trap just fixed:
The SVD exists for every matrix, any shape, any rank — that universality is its whole point.
Now break it on purpose. Build a case that violates it and follow the consequences until something visibly fails. Where does the failure first show up — and would you have noticed it if you had not been looking?
Hint: The dangerous rules are the ones whose violation still produces an answer. If yours fails loudly, try to find one that fails quietly.
Answer:
This confuses SVD with EIGENdecomposition. A 2×3 matrix has no eigenvalues (Az = λz needs the output to have the same shape as z), so people wrongly conclude it has no factorization at all.
Concept
The SVD reads the rank of A straight off Σ: the rank is simply how many singular values are strictly positive. Every zero σᵢ marks a direction A collapses to nothing.
\[ \operatorname{rank}(A) = \#\{\,i : \sigma_i > 0\,\} \]
Our A has σ = [5, 3, 0], so its rank is 2, not 3 — even though it has three columns. This is also the numerically honest way to find rank: count singular values above a small tolerance, rather than trusting an exact det = 0.
Section
Part 3 of 8
Concept
Multiply A = UΣVᵀ on the right by V. Since VᵀV = I, the Vᵀ and V collapse and we get AV = UΣ. Reading that one column at a time gives the identity that IS the geometry:
\[ A\,v_i = \sigma_i\, u_i \]
In words: the right singular vector vᵢ is an input direction that A sends to the output direction uᵢ, stretched by exactly σᵢ. A maps the orthonormal v-frame to the orthonormal u-frame, with a clean scale on each axis.
Pattern
Predict first
The table runs: ‖A v₁‖ | 5.0 (= σ₁) · (A v₁)/σ₁ | [−0.7071, −0.7071]
In Verify A v₁ = σ₁ u₁ in code, given the rows so far: what is the next one — the row where quantity is u₁ (col 0 of U)?
Correct: u₁ (col 0 of U) | [−0.7071, −0.7071] ✓ same
| quantity | value (verified) |
|---|---|
| ‖A v₁‖ | 5.0 (= σ₁) |
| (A v₁)/σ₁ | [−0.7071, −0.7071] |
| u₁ (col 0 of U) | [−0.7071, −0.7071] ✓ same |
Why: The relationship between the columns, not the individual numbers, is what generates the next row. A stretches the v₁ direction by exactly 5 and lands on u₁ — the identity A v₁ = σ₁ u₁ holds numerically.
Worked example
Take the top right singular vector v₁ (first row of Vt), push it through A, and check the length is σ₁ and the direction is u₁. Standalone:
import numpy as np
A = np.array([[3.,2.,2.],[2.,3.,-2.]])
U, S, Vt = np.linalg.svd(A, full_matrices=False)
v1 = Vt[0]
print(np.linalg.norm(A @ v1).round(4), S[0].round(4)) # 5.0 5.0
print((A @ v1 / S[0]).round(4), U[:,0].round(4)) # equal -> u1‖A v₁‖ = 5.0 = σ₁, and (A v₁)/σ₁ = u₁
Why: A stretches the v₁ direction by exactly 5 and lands on u₁ — the identity A v₁ = σ₁ u₁ holds numerically. This is the ellipse's long axis being traced out.
| quantity | value (verified) |
|---|---|
| ‖A v₁‖ | 5.0 (= σ₁) |
| (A v₁)/σ₁ | [−0.7071, −0.7071] |
| u₁ (col 0 of U) | [−0.7071, −0.7071] ✓ same |
Reverse engineer
Discussion prompt
Work backwards. The example finished here:
‖A v₁‖ = 5.0 = σ₁, and (A v₁)/σ₁ = u₁
What was it asked to do, and what must it have been given? Reconstruct the problem from its answer.
Hint: Every quantity in the result had to enter somewhere. Account for each one.
Answer:
Take the top right singular vector v₁ (first row of Vt), push it through A, and check the length is σ₁ and the direction is u₁. Standalone:
Intuition
Among all unit input vectors, which one does A stretch the most? Answer: v₁, and the stretch is σ₁. The largest singular value is the operator's maximum gain — the spectral norm ‖A‖₂ = σ₁.
Likewise σ₂ is the most A can stretch anything perpendicular to v₁, and so on down. That ordering — biggest stretch first — is why we sort Σ large to small, and it's what makes truncation meaningful later.
Sorting
Sort into buckets
These are the pieces of Lesson 16: Singular Value Decomposition, out of order. Put each one back under the part of the lesson it belongs to.
Section
Part 4 of 8 — our 2×3, start to finish
Worked example
(AᵀA)ⱼₖ is column j of A dotted with column k. The three columns are [3,2], [2,3], [2,−2]. Fill the symmetric 3×3:
Diagonal: c₁·c₁ = 9+4 = 13, c₂·c₂ = 4+9 = 13, c₃·c₃ = 4+4 = 8
Why: Each diagonal entry is a column dotted with itself — the squared length of that column.
Off-diagonals: c₁·c₂ = 6+6 = 12, c₁·c₃ = 6−4 = 2, c₂·c₃ = 4−6 = −2
Why: Symmetric, so we only compute the upper triangle and mirror it.
\[ A^\top A = \begin{bmatrix} 13 & 12 & 2 \\ 12 & 13 & -2 \\ 2 & -2 & 8 \end{bmatrix} \]
| entry | columns dotted | value |
|---|---|---|
| (1,1) | [3,2]·[3,2] | 13 |
| (1,2) | [3,2]·[2,3] | 12 |
| (1,3) | [3,2]·[2,−2] | 2 |
| (3,3) | [2,−2]·[2,−2] | 8 |
Notation
Annotate
From Step 1 — form AᵀA, entry by entry — read this one piece at a time. What is each part doing?
On: \( A^\top A = \begin{bmatrix} 13 & 12 & 2 \\ 12 & 13 & -2 \\ 2 & -2 & 8 \end{bmatrix} \)
Worked example
The eigenvalues of that 3×3 are 25, 9, 0 (trace = 13+13+8 = 46 = 25+9+0 checks out; the matrix is rank 2 so one eigenvalue must be 0). Take non-negative roots of the nonzero ones:
\[ \sigma_1 = \sqrt{25} = 5, \qquad \sigma_2 = \sqrt{9} = 3, \qquad \sigma_3 = \sqrt{0} = 0 \]
Keep the two nonzero singular values: σ = [5, 3]
Why: A zero singular value means a direction A collapses entirely — it carries no information, so the thin SVD drops it. Two nonzero σ ⇒ rank 2.
| λᵢ(AᵀA) | σᵢ = √λᵢ | kept? |
|---|---|---|
| 25 | 5 | yes |
| 9 | 3 | yes |
| 0 | 0 | no (rank-deficient direction) |
Translation
\( \sigma_1 = \sqrt{25} = 5, \qquad \sigma_2 = \sqrt{9} = 3, \qquad \sigma_3 = \sqrt{0} = 0 \)
Draw it
Translate both ways. First write the expression above as a sentence with no symbols in it at all. Then cover it, and write your sentence back as notation. If the two versions disagree, the disagreement is the thing to fix.
Fill the middle
Fill in the blanks
From Step 3 — decompose and reconstruct — one line has had its right-hand side removed. Put it back.
import numpy as np
A = np.array([[3.,2.,2.],[2.,3.,-2.]])
U, S, Vt = np.linalg.svd(A, full_matrices=False)
print(S.round(4)) # [5. 3.]
print(U.round(4))
print(Vt.round(4))
print(np.allclose(U @ np.diag(S) @ Vt, A)) # True
Why: A is what everything below it consumes, so the wrong expression here fails later and somewhere else. S is the vector of singular values, not a matrix — rebuild Σ with np.diag(S).
Worked example
Let NumPy produce U, Σ, Vᵀ, print them, and verify UΣVᵀ rebuilds A exactly. Standalone block:
import numpy as np
A = np.array([[3.,2.,2.],[2.,3.,-2.]])
U, S, Vt = np.linalg.svd(A, full_matrices=False)
print(S.round(4)) # [5. 3.]
print(U.round(4))
print(Vt.round(4))
print(np.allclose(U @ np.diag(S) @ Vt, A)) # Truenp.linalg.svd returns Vt (already transposed) and a 1-D S
Why: S is the vector of singular values, not a matrix — rebuild Σ with np.diag(S). The rows of Vt are the vᵢ; the columns of U are the uᵢ. allclose confirms the factorization is exact.
| object | value (verified) |
|---|---|
| S | [5, 3] |
| U | [[−0.7071, −0.7071], [−0.7071, 0.7071]] |
| Vt (rows = vᵢ) | [[−0.7071, −0.7071, 0], [−0.2357, 0.2357, −0.9428]] |
| UΣVᵀ == A | True |
Concept
If A vᵢ = σᵢ uᵢ, then flipping both vᵢ → −vᵢ and uᵢ → −uᵢ still satisfies it — the product σᵢ uᵢ vᵢᵀ is unchanged. So NumPy's U and Vt may come back with different signs on another machine, but the reconstruction is identical.
The singular values Σ are always the same and always non-negative. Never test U/V for equality against a printout — test the reconstruction UΣVᵀ ≈ A, which is what actually matters.
Section
Part 5 of 8 — Eckart-Young
Concept
Expand A = UΣVᵀ as a sum of outer products, one per singular triple. Each term σᵢ uᵢ vᵢᵀ is a rank-1 matrix scaled by its singular value:
\[ A = \sum_{i=1}^{r} \sigma_i\, u_i v_i^\top = \sigma_1 u_1 v_1^\top + \sigma_2 u_2 v_2^\top + \dots \]
The terms are sorted by importance: σ₁ biggest first. The matrix is built mostly from its leading terms — the tail contributes little. That is the whole idea behind compression.
Faded example
Fill in the blanks
The pieces really do add back to A, with the scaffolding fading: two lines are gone now — fill both.
import numpy as np
A = np.array([[3.,2.,2.],[2.,3.,-2.]])
U, S, Vt = np.linalg.svd(A, full_matrices=False)
A1 = *S[0]np.outer(U[:,0], Vt[0]) # rank-1 piece
A2 = S[1]np.outer(U[:,1], Vt[1])* # rank-2 piece
print((A1 + A2).round(4)) # == A
print(np.allclose(A1 + A2, A)) # True
Why: Reproducing these unaided, rather than reading them, is what tells you the method has transferred. Each np.outer(U[:,i], Vt[i]) is the rank-1 outer product uᵢvᵢᵀ; scaling by σᵢ and summing rebuilds A.
Worked example
Confirm the sum-of-outer-products claim on our running matrix: build both rank-1 pieces σ₁u₁v₁ᵀ and σ₂u₂v₂ᵀ and add them. Standalone:
import numpy as np
A = np.array([[3.,2.,2.],[2.,3.,-2.]])
U, S, Vt = np.linalg.svd(A, full_matrices=False)
A1 = S[0]*np.outer(U[:,0], Vt[0]) # rank-1 piece
A2 = S[1]*np.outer(U[:,1], Vt[1]) # rank-2 piece
print((A1 + A2).round(4)) # == A
print(np.allclose(A1 + A2, A)) # TrueA₁ + A₂ = A exactly (rank 2 has only two pieces)
Why: Each np.outer(U[:,i], Vt[i]) is the rank-1 outer product uᵢvᵢᵀ; scaling by σᵢ and summing rebuilds A. Since rank(A)=2, two pieces suffice — the (zeroed) third triple adds nothing.
| piece | matrix (verified) |
|---|---|
| σ₁u₁v₁ᵀ | [[2.5, 2.5, 0], [2.5, 2.5, 0]] |
| σ₂u₂v₂ᵀ | [[0.5, −0.5, 2], [−0.5, 0.5, −2]] |
| sum | [[3, 2, 2], [2, 3, −2]] = A |
Concept
Keep only the top k triples and throw the rest away. The result Aₖ is a rank-k matrix that approximates A:
\[ A_k = \sum_{i=1}^{k} \sigma_i\, u_i v_i^\top \]
Eckart-Young theorem: Aₖ is the closest rank-k matrix to A in both the spectral and Frobenius norms. No other rank-k matrix does better — the SVD truncation is provably optimal, not just convenient.
Concept
Dropping triples k+1, k+2, … leaves exactly those terms as the error A − Aₖ. Its spectral norm is the largest thing you dropped — the first discarded singular value:
\[ \lVert A - A_k \rVert_2 = \sigma_{k+1}, \qquad \lVert A - A_k \rVert_F = \sqrt{\sigma_{k+1}^2 + \sigma_{k+2}^2 + \dots} \]
So the error is governed entirely by the tail of the singular-value spectrum. If the σᵢ decay fast, a tiny k captures almost everything.
Estimation
Predict first
For our A with σ = [5, 3], the rank-1 approximation keeps only σ₁ u₁ v₁ᵀ, and Eckart-Young predicts the error is exactly σ₂ = 3. Build it and measure:
Commit before you compute: what does Rank-1 approximation of our A come out to? A rough magnitude and the right form is enough — the point is to have something concrete to be wrong about.
Correct: A₁ = [[2.5, 2.5, 0], [2.5, 2.5, 0]], error = 3.0 = σ₂
Why: A prediction you can defend turns the computation into a check rather than a leap of faith — and an answer that contradicts it is caught on the spot. np.outer(U[:,0], Vt[0]) builds the rank-1 outer product u₁v₁ᵀ; scaling by σ₁ gives the best rank-1 piece.
Worked example
For our A with σ = [5, 3], the rank-1 approximation keeps only σ₁ u₁ v₁ᵀ, and Eckart-Young predicts the error is exactly σ₂ = 3. Build it and measure:
import numpy as np
A = np.array([[3.,2.,2.],[2.,3.,-2.]])
U, S, Vt = np.linalg.svd(A, full_matrices=False)
A1 = S[0] * np.outer(U[:,0], Vt[0])
print(A1.round(4))
print(round(np.linalg.norm(A - A1, 2), 4)) # 3.0
print(round(S[1], 4)) # 3.0A₁ = [[2.5, 2.5, 0], [2.5, 2.5, 0]], error = 3.0 = σ₂
Why: np.outer(U[:,0], Vt[0]) builds the rank-1 outer product u₁v₁ᵀ; scaling by σ₁ gives the best rank-1 piece. Dropping σ₂ costs exactly its size — Eckart-Young made concrete.
| quantity | value (verified) |
|---|---|
| A₁ (rank-1) | [[2.5, 2.5, 0], [2.5, 2.5, 0]] |
| ‖A − A₁‖₂ (spectral) | 3.0 |
| ‖A − A₁‖_F (Frobenius) | 3.0 |
| σ₂ | 3.0 ✓ matches |
Anomaly
Predict first
A student writes this, and it looks reasonable:
To compress, drop the largest singular values — surely the big numbers are the ones eating all the storage.
It is wrong. Say what breaks — and say it before you turn the page.
Correct: Backwards. σ₁ carries the MOST structure — for our A it holds A₁ with error only 3 out of the total.
Keep the largest singular values; discard the smallest tail.
Why: Backwards. σ₁ carries the MOST structure — for our A it holds A₁ with error only 3 out of the total. Removing σ₁ instead would leave error 5 (= σ₁), the worst possible rank-1 choice.
Trap
To compress, drop the largest singular values — surely the big numbers are the ones eating all the storage.
Truncate σ₁ first, keep the tail
Why: Backwards. σ₁ carries the MOST structure — for our A it holds A₁ with error only 3 out of the total. Removing σ₁ instead would leave error 5 (= σ₁), the worst possible rank-1 choice.
Keep the largest singular values; discard the smallest tail.
Keep σ₁…σₖ (largest), drop σₖ₊₁…
Why: The largest σ capture the most stretch/energy, so keeping them minimizes the error (= σₖ₊₁). Image compression, PCA, and LoRA all keep the top few and drop the rest. Big σ = keep; small σ = throw away.
Section
Part 6 of 8 — a decaying spectrum
Intuition
A full m×n matrix costs m·n numbers. A rank-k truncation stores only the top k triples: k columns of U (m·k), k singular values (k), and k rows of Vᵀ (k·n) — about k(m+n+1) numbers.
When k is much smaller than m and n, that is a huge saving — and if the spectrum decays fast, the reconstruction is nearly perfect. This is exactly how SVD image compression works.
Fill the middle
Fill in the blanks
From A matrix with a decaying spectrum — one line has had its right-hand side removed. Put it back.
import numpy as np
np.random.seed(0)
Q1, _ = np.linalg.qr(np.random.randn(5, 4))
Q2, _ = np.linalg.qr(np.random.randn(4, 4))
B = Q1 @ np.diag([9.,5.,2.,0.5]) @ Q2.T # known singular values
U, S, Vt = np.linalg.svd(B, full_matrices=False)
print(S.round(4)) # [9. 5. 2. 0.5]
energy = np.cumsum(S2) / np.sum(S2)
print((energy*100).round(2)) # cumulative % energy
Why: energy is what everything below it consumes, so the wrong expression here fails later and somewhere else. Energy is Frobenius energy = Σσᵢ².
Worked example
Build a 5×4 matrix whose singular values we chose to be [9, 5, 2, 0.5], then read the cumulative energy — the fraction of total Σσᵢ² captured by the top k. Seeded, standalone:
import numpy as np
np.random.seed(0)
Q1, _ = np.linalg.qr(np.random.randn(5, 4))
Q2, _ = np.linalg.qr(np.random.randn(4, 4))
B = Q1 @ np.diag([9.,5.,2.,0.5]) @ Q2.T # known singular values
U, S, Vt = np.linalg.svd(B, full_matrices=False)
print(S.round(4)) # [9. 5. 2. 0.5]
energy = np.cumsum(S**2) / np.sum(S**2)
print((energy*100).round(2)) # cumulative % energyTop-1 already captures 73.5% of the energy; top-2 captures 96.2%
Why: Energy is Frobenius energy = Σσᵢ². Because σ₁=9 dominates, one component holds most of the matrix. The last component (σ₄=0.5) adds a mere 0.23% — cheap to drop.
| k | singular values kept | cumulative energy |
|---|---|---|
| 1 | [9] | 73.47% |
| 2 | [9, 5] | 96.15% |
| 3 | [9, 5, 2] | 99.77% |
| 4 | [9, 5, 2, 0.5] | 100.00% |
Pattern
Step through it
Step through A matrix with a decaying spectrum one row at a time. What is driving the change, and what would the row after the last one be?
Worked example
Truncate B to rank 2 and confirm both Eckart-Young formulas: spectral error = σ₃ = 2, Frobenius error = √(σ₃² + σ₄²). Standalone:
import numpy as np
np.random.seed(0)
Q1, _ = np.linalg.qr(np.random.randn(5, 4))
Q2, _ = np.linalg.qr(np.random.randn(4, 4))
B = Q1 @ np.diag([9.,5.,2.,0.5]) @ Q2.T
U, S, Vt = np.linalg.svd(B, full_matrices=False)
B2 = U[:, :2] @ np.diag(S[:2]) @ Vt[:2] # best rank-2
print(round(np.linalg.norm(B - B2, 2), 4)) # 2.0 = sigma3
print(round(np.linalg.norm(B - B2, 'fro'), 4)) # 2.0616
print(round(np.sqrt(S[2]**2 + S[3]**2), 4)) # 2.0616Spectral error 2.0 = σ₃; Frobenius error 2.0616 = √(2² + 0.5²)
Why: Slicing [:, :2] and [:2] keeps the first two columns/rows — the top-2 triples. Both norms match their Eckart-Young predictions to the decimal. The theory is not approximate; it is exact.
| norm | measured | theory |
|---|---|---|
| ‖B − B₂‖₂ | 2.0 | σ₃ = 2.0 |
| ‖B − B₂‖_F | 2.0616 | √(σ₃²+σ₄²) = √4.25 = 2.0616 |
Comparison
Comparison matrix
From Rank-2 error matches the theory: refill the theory column from what you know. The rest of the table is as it appeared.
| norm | measured | theory |
|---|---|---|
| ‖B − B₂‖₂ | 2.0 | σ₃ = 2.0 |
| ‖B − B₂‖_F | 2.0616 | √(σ₃²+σ₄²) = √4.25 = 2.0616 |
Section
Part 7 of 8
Concept
In Lesson 10 PCA meant: center the data, form the covariance C = XcᵀXc/(n−1), and take its top eigenvectors. But Xc = UΣVᵀ already, so C = V Σ²Vᵀ/(n−1) — the principal components V fall straight out of the SVD of Xc. No covariance matrix needed.
\[ C = \frac{X_c^\top X_c}{n-1} = V\,\frac{\Sigma^2}{n-1}\,V^\top \;\Longrightarrow\; \text{PCs} = V,\quad \text{variances} = \frac{\sigma_i^2}{n-1} \]
The right singular vectors of centered X are the principal directions, and each σᵢ²/(n−1) is the variance captured along that direction. SVD is the numerically stable way to do PCA.
Fill the middle
Fill in the blanks
From PCA of a small dataset via SVD — one line has had its right-hand side removed. Put it back.
import numpy as np
Xdata = np.array([[2.,0.],[3.,1.],[4.,3.],[5.,4.],[6.,5.],[4.,2.]])
Xc = Xdata - Xdata.mean(axis=0) # center
U, S, Vt = np.linalg.svd(Xc, full_matrices=False)
print(Vt.round(4)) # rows = principal components
print((S**2/(len(Xdata)-1)).round(4)) # variance along each PC
Why: Xc is what everything below it consumes, so the wrong expression here fails later and somewhere else. Subtracting the column mean centers the cloud at the origin.
Worked example
Six 2-D points, correlated. Center them, take the SVD, and read the principal components off Vt. Standalone:
import numpy as np
Xdata = np.array([[2.,0.],[3.,1.],[4.,3.],[5.,4.],[6.,5.],[4.,2.]])
Xc = Xdata - Xdata.mean(axis=0) # center
U, S, Vt = np.linalg.svd(Xc, full_matrices=False)
print(Vt.round(4)) # rows = principal components
print((S**2/(len(Xdata)-1)).round(4)) # variance along each PCPC1 = [0.6012, 0.7991]; variances = [5.456, 0.044]
Why: Subtracting the column mean centers the cloud at the origin. The first row of Vt is the direction of maximum spread, and σ₁²/(n−1)=5.456 is the variance along it — 99.2% of the total.
| quantity | value (verified) |
|---|---|
| mean removed | [4.0, 2.5] |
| PC1 (row 0 of Vt) | [0.6012, 0.7991] |
| PC2 (row 1 of Vt) | [−0.7991, 0.6012] |
| variances σ²/(n−1) | [5.456, 0.044] |
Pattern
Predict first
The table runs: eig of covariance | [5.456, 0.044] · σ²/(n−1) from SVD | [5.456, 0.044]
In Cross-check against the covariance eigenvectors, given the rows so far: what is the next one — the row where method is agree??
Correct: agree? | yes — identical
| method | variances (verified) |
|---|---|
| eig of covariance | [5.456, 0.044] |
| σ²/(n−1) from SVD | [5.456, 0.044] |
| agree? | yes — identical |
Why: The relationship between the columns, not the individual numbers, is what generates the next row. The two paths land on identical variances.
Worked example
The slow way — eigendecompose the covariance matrix — must give the same variances. Confirm the two routes agree:
import numpy as np
Xdata = np.array([[2.,0.],[3.,1.],[4.,3.],[5.,4.],[6.,5.],[4.,2.]])
Xc = Xdata - Xdata.mean(axis=0)
U, S, Vt = np.linalg.svd(Xc, full_matrices=False)
cov = (Xc.T @ Xc) / (len(Xdata)-1)
w, V = np.linalg.eigh(cov)
print(np.sort(w)[::-1].round(4)) # [5.456 0.044]
print((S**2/(len(Xdata)-1)).round(4)) # [5.456 0.044] -- matchCovariance eigenvalues [5.456, 0.044] = SVD's σ²/(n−1)
Why: The two paths land on identical variances. SVD skips forming XcᵀXc — which squares the condition number — so it's the numerically preferred route to the exact same PCA.
| method | variances (verified) |
|---|---|
| eig of covariance | [5.456, 0.044] |
| σ²/(n−1) from SVD | [5.456, 0.044] |
| agree? | yes — identical |
Section
Part 8 of 8 — A⁺ = VΣ⁺Uᵀ
Intuition
If A is rotate (Vᵀ) → scale (Σ) → rotate (U), then undoing A should just run those steps backward: undo the last rotation, undo the scaling, undo the first rotation.
Un-rotating is transposing (U → Uᵀ, Vᵀ → V), and un-scaling is dividing by each σᵢ. Reversed order gives V · Σ⁺ · Uᵀ — that's the pseudoinverse, before we write a single formula.
The only step that can't be undone is a zero stretch: once a direction is squashed to nothing, no inverse can bring it back. So Σ⁺ leaves those zeros as zeros — the one honest compromise.
Concept
A rectangular or rank-deficient A has no true inverse. But the SVD gives the next best thing — the Moore-Penrose pseudoinverse — by inverting each nonzero singular value and transposing the rotations:
\[ A^+ = V\,\Sigma^+ U^\top, \qquad \Sigma^+ = \operatorname{diag}\!\left(\tfrac{1}{\sigma_1}, \tfrac{1}{\sigma_2}, \dots\right) \]
For each nonzero σᵢ we store 1/σᵢ; zero singular values stay zero (you can't invert a collapse). Undo the two rotations by transposing them. That's the whole recipe.
Explain it
Discussion prompt
Explain Inverting the un-invertible to a student a year behind you. No notation, no jargon they have not met — and it still has to be true.
Hint: If your explanation needs a symbol they have never seen, you are describing the notation rather than the idea.
Answer:
For each nonzero σᵢ we store 1/σᵢ; zero singular values stay zero (you can't invert a collapse). Undo the two rotations by transposing them. That's the whole recipe.
Concept
For an overdetermined system Xw = y (Lesson 7), w = X⁺y is the least-squares solution — the same w the normal equations XᵀXw = Xᵀy give when XᵀX is invertible, and the minimum-norm solution when it isn't.
So the pseudoinverse is a single, always-defined formula that covers OLS, ridge-free rank deficiency, and underdetermined systems at once. It never chokes on a singular XᵀX the way solve(XᵀX, Xᵀy) does.
Missing information
Discussion prompt
Reuse Lesson 7's five-student data. np.linalg.pinv(X) @ y must return the least-squares line [2.2, 0.6]. Standalone:
What do you need to know — or decide — before the first line can be written? List everything the problem has to hand you.
Hint: Anything you would have to invent to get started is a thing the problem must supply.
Answer:
np.c_[np.ones(5), x] prepends the bias column. pinv builds X⁺ from X's SVD and applies it to y — the SVD-based least-squares solver, stable even when XᵀX is near-singular. Same answer as Lesson 7.
Worked example
Reuse Lesson 7's five-student data. np.linalg.pinv(X) @ y must return the least-squares line [2.2, 0.6]. Standalone:
import numpy as np
x = np.array([1.,2.,3.,4.,5.]); y = np.array([2.,4.,5.,4.,5.])
X = np.c_[np.ones(5), x]
w = np.linalg.pinv(X) @ y # A+ = V Σ+ Uᵀ, applied to y
print(w.round(4)) # [2.2 0.6]
print(np.linalg.svd(X, compute_uv=False).round(4)) # singular values of Xpinv(X) @ y = [2.2, 0.6] — identical to the normal equations
Why: np.c_[np.ones(5), x] prepends the bias column. pinv builds X⁺ from X's SVD and applies it to y — the SVD-based least-squares solver, stable even when XᵀX is near-singular. Same answer as Lesson 7.
| method | w (verified) |
|---|---|
| normal equations (Lesson 7) | [2.2, 0.6] |
| pinv(X) @ y (SVD) | [2.2, 0.6] |
| singular values of X | [7.6912, 0.9194] |
Estimation
Predict first
Prove pinv is nothing magic: assemble X⁺ = VΣ⁺Uᵀ yourself from svd, and check it equals np.linalg.pinv(X) and solves the same system. Standalone:
Commit before you compute: what does Build A⁺ by hand from the SVD come out to? A rough magnitude and the right form is enough — the point is to have something concrete to be wrong about.
Correct: Vt.T @ diag(1/S) @ U.T equals np.linalg.pinv(X) exactly
Why: A prediction you can defend turns the computation into a check rather than a leap of faith — and an answer that contradicts it is caught on the spot. We transpose Vt back to V, invert each singular value with 1/S on the diagonal, and transpose U.
Worked example
Prove pinv is nothing magic: assemble X⁺ = VΣ⁺Uᵀ yourself from svd, and check it equals np.linalg.pinv(X) and solves the same system. Standalone:
import numpy as np
x = np.array([1.,2.,3.,4.,5.]); y = np.array([2.,4.,5.,4.,5.])
X = np.c_[np.ones(5), x]
U, S, Vt = np.linalg.svd(X, full_matrices=False)
Xplus = Vt.T @ np.diag(1.0/S) @ U.T # V Σ+ Uᵀ, built by hand
print(np.allclose(Xplus, np.linalg.pinv(X))) # True
print((Xplus @ y).round(4)) # [2.2 0.6]Vt.T @ diag(1/S) @ U.T equals np.linalg.pinv(X) exactly
Why: We transpose Vt back to V, invert each singular value with 1/S on the diagonal, and transpose U. The hand-built X⁺ matches the library and returns [2.2, 0.6] — the pseudoinverse IS this formula.
| check | result (verified) |
|---|---|
| Σ⁺ = diag(1/σᵢ) | diag(1/7.6912, 1/0.9194) |
| VΣ⁺Uᵀ == np.linalg.pinv(X) | True |
| (VΣ⁺Uᵀ) @ y | [2.2, 0.6] |
Reverse engineer
Discussion prompt
Work backwards. The example finished here:
Vt.T @ diag(1/S) @ U.T equals np.linalg.pinv(X) exactly
What was it asked to do, and what must it have been given? Reconstruct the problem from its answer.
Hint: Every quantity in the result had to enter somewhere. Account for each one.
Answer:
Prove pinv is nothing magic: assemble X⁺ = VΣ⁺Uᵀ yourself from svd, and check it equals np.linalg.pinv(X) and solves the same system. Standalone:
Anomaly
Predict first
A student writes this, and it looks reasonable:
The formula is Σ⁺ = diag(1/σᵢ), so just take the reciprocal of every diagonal entry — including any zeros.
It is wrong. Say what breaks — and say it before you turn the page.
Correct: A zero singular value marks a direction A collapses to nothing — there's no information to recover, so 1/0 is meaningless and would poison the whole matrix with infinities.
Invert only the nonzero singular values; leave zeros (and numerically tiny σ below a tolerance) as 0 in Σ⁺.
Why: A zero singular value marks a direction A collapses to nothing — there's no information to recover, so 1/0 is meaningless and would poison the whole matrix with infinities. Our A has σ₃ = 0 exactly.
Trap
The formula is Σ⁺ = diag(1/σᵢ), so just take the reciprocal of every diagonal entry — including any zeros.
1/σ₃ = 1/0 → ∞ (or a NaN blow-up)
Why: A zero singular value marks a direction A collapses to nothing — there's no information to recover, so 1/0 is meaningless and would poison the whole matrix with infinities. Our A has σ₃ = 0 exactly.
Invert only the nonzero singular values; leave zeros (and numerically tiny σ below a tolerance) as 0 in Σ⁺.
Σ⁺ = diag(1/σ₁, 1/σ₂, 0)
Why: This is exactly what np.linalg.pinv does: it zeros out σ below rcond·σ_max before inverting. The result is the minimum-norm least-squares solution — finite and well-defined even for rank-deficient A.
Two truths and a lie
Sort into buckets
Some of these hold up and some are the exact mistakes this lesson is built to prevent. Sort them.
m×n matrix A — square or not, symmetric or not, full-rank or not — factors into two orthogonal matrices and one diagonal matrix of non-negative numbers:; For an m×n matrix A, the factors have fixed shapes. In the full SVD, U is m×m, Σ is m×n (rectangular, zeros off the diagonal), and Vᵀ is n×n.; Orthogonal matrices are rigid rotations (or reflections): they never stretch, they only reorient. All the stretching a linear map does is packed into the diagonal Σ.2×3 has no decomposition.; To compress, drop the largest singular values — surely the big numbers are the ones eating all the storage.Constraint
Discussion prompt
Run The SVD toolkit with this step confiscated:
Compress: keep the top-k triples → best rank-k (Eckart-Young); spectral error = σₖ₊₁, Frobenius error = √(Σ_{i>k} σᵢ²)
Is it still possible? If it is, say what takes its place and what it costs you. If it is not, say exactly what that step was providing that nothing else does.
Hint: A step you can drop for free was never load-bearing. If you cannot drop it, name the thing that goes wrong the moment it is gone.
Answer:
U, S, Vt = np.linalg.svd(A, full_matrices=False) — works for any shape; S is a 1-D vector, Vt is already transposedσᵢ = √λᵢ(AᵀA), always real and ≥ 0; V = eigenvectors of AᵀA, U = eigenvectors of AAᵀA vᵢ = σᵢ uᵢ; σ₁ = ‖A‖₂ is the maximum stretchk triples → best rank-k (Eckart-Young); spectral error = σₖ₊₁, Frobenius error = √(Σ_{i>k} σᵢ²)X; variances = σᵢ²/(n−1)A⁺ = VΣ⁺Uᵀ (invert only nonzero σ); np.linalg.pinv does it stablyPattern
U, S, Vt = np.linalg.svd(A, full_matrices=False) — works for any shape; S is a 1-D vector, Vt is already transposedσᵢ = √λᵢ(AᵀA), always real and ≥ 0; V = eigenvectors of AᵀA, U = eigenvectors of AAᵀA vᵢ = σᵢ uᵢ; σ₁ = ‖A‖₂ is the maximum stretchk triples → best rank-k (Eckart-Young); spectral error = σₖ₊₁, Frobenius error = √(Σ_{i>k} σᵢ²)X; variances = σᵢ²/(n−1)A⁺ = VΣ⁺Uᵀ (invert only nonzero σ); np.linalg.pinv does it stablyEdge cases
Discussion prompt
The SVD toolkit works on the cases you have just seen. Push it to the edge: what is the most degenerate input it still handles — empty, zero, one item, everything equal — and what is the first case where it stops being true? Name the case, not just "it breaks".
Hint: Try the smallest legal input, then the largest, then the one where two things collide. Methods are specified at their edges; the middle takes care of itself.
Answer:
U, S, Vt = np.linalg.svd(A, full_matrices=False) — works for any shape; S is a 1-D vector, Vt is already transposedσᵢ = √λᵢ(AᵀA), always real and ≥ 0; V = eigenvectors of AᵀA, U = eigenvectors of AAᵀA vᵢ = σᵢ uᵢ; σ₁ = ‖A‖₂ is the maximum stretchk triples → best rank-k (Eckart-Young); spectral error = σₖ₊₁, Frobenius error = √(Σ_{i>k} σᵢ²)X; variances = σᵢ²/(n−1)A⁺ = VΣ⁺Uᵀ (invert only nonzero σ); np.linalg.pinv does it stablyCheck
Trace the derivation back to its source. What object are the σᵢ built from?
Check your understanding
The singular values of a matrix A are…
Answer: A
Why: Substituting A = UΣVᵀ into AᵀA gives V diag(σᵢ²) Vᵀ, so λᵢ(AᵀA) = σᵢ² and σᵢ = √λᵢ(AᵀA). Since AᵀA is PSD its eigenvalues are ≥ 0, so the σᵢ are real and non-negative — and exist for any shape of A.
Elimination
Eliminate the wrong options
A has singular values [10, 6, 2]. The spectral-norm error ‖A − A₂‖₂ of the best rank-2 approximation is:
3 of these 4 are wrong. Strike them one at a time, and say what rules each one out before you strike the next. The survivor is the answer.
Survives elimination: A
Why: Eckart-Young: the best rank-k approximation's spectral-norm error equals the first discarded singular value, σ_{k+1}. Keeping the top 2 (10 and 6) discards σ₃ = 2, so the error is exactly 2.
Check
Reason directly from the singular-value spectrum — no computation needed.
Check your understanding
A has singular values [10, 6, 2]. The spectral-norm error ‖A − A₂‖₂ of the best rank-2 approximation is:
Answer: A
Why: Eckart-Young: the best rank-k approximation's spectral-norm error equals the first discarded singular value, σ_{k+1}. Keeping the top 2 (10 and 6) discards σ₃ = 2, so the error is exactly 2.
Prediction
Predict first
In A = UΣVᵀ, what does the identity A vᵢ = σᵢ uᵢ say?
Answer it in your own words, now, with nothing to choose from. The options are on the next slide — and picking the right one off a list is an easier skill than producing it.
Correct: A sends the right singular vector vᵢ to the left singular vector uᵢ, stretched by σᵢ
Why: Right-multiplying A = UΣVᵀ by V gives AV = UΣ; column i reads A vᵢ = σᵢ uᵢ. A maps the input frame {vᵢ} to the output frame {uᵢ}, scaling axis i by σᵢ — the circle-to-ellipse picture.
Check
Picture the unit circle mapping to an ellipse.
Check your understanding
In A = UΣVᵀ, what does the identity A vᵢ = σᵢ uᵢ say?
Answer: A
Why: Right-multiplying A = UΣVᵀ by V gives AV = UΣ; column i reads A vᵢ = σᵢ uᵢ. A maps the input frame {vᵢ} to the output frame {uᵢ}, scaling axis i by σᵢ — the circle-to-ellipse picture.
Elimination
Eliminate the wrong options
For an overdetermined system Xw = y (more rows than columns), np.linalg.pinv(X) @ y returns:
3 of these 4 are wrong. Strike them one at a time, and say what rules each one out before you strike the next. The survivor is the answer.
Survives elimination: A
Why: X⁺ = VΣ⁺Uᵀ gives the minimum-norm least-squares solution. For full-rank X it equals (XᵀX)⁻¹Xᵀy — the normal-equations answer, [2.2, 0.6] on Lesson 7's data.
Check
What does A⁺ actually compute for a tall system?
Check your understanding
For an overdetermined system Xw = y (more rows than columns), np.linalg.pinv(X) @ y returns:
Answer: A
Why: X⁺ = VΣ⁺Uᵀ gives the minimum-norm least-squares solution. For full-rank X it equals (XᵀX)⁻¹Xᵀy — the normal-equations answer, [2.2, 0.6] on Lesson 7's data.
Section
The project
Concept
Three milestones on the running matrices: SVD A and verify reconstruction, build its rank-1 approximation and confirm the Eckart-Young error, then solve Lesson 7's OLS with the pseudoinverse. You've derived every piece — now assemble it.
| # | requirement | tool |
|---|---|---|
| 1 | SVD A and verify UΣVᵀ = A | np.linalg.svd |
| 2 | rank-1 approx; error == σ₂ | np.outer + np.linalg.norm |
| 3 | pseudoinverse OLS → [2.2, 0.6] | np.linalg.pinv |
Build rules: type every line yourself, use full_matrices=False so shapes line up, remember svd returns Vt (already transposed) and a 1-D S, and reconstruct Σ with np.diag(S).
Analogy
Discussion prompt
Explain Project: decompose, compress, invert by analogy to something with no Machine Learning in it at all — a queue, a recipe, a map, a bank balance, whatever fits. Then say where your analogy breaks.
Hint: An analogy that never breaks is not an analogy, it is the same idea wearing a hat. Find the seam — that is the part that is actually new.
Answer:
Build rules: type every line yourself, use full_matrices=False so shapes line up, remember svd returns Vt (already transposed) and a 1-D S, and reconstruct Σ with np.diag(S).
Worked example
Your turn: SVD A = [[3,2,2],[2,3,−2]] and verify UΣVᵀ = A. Predict the singular values out loud before you print — you built them from AᵀA.
Hint: U, S, Vt = np.linalg.svd(A, full_matrices=False); reconstruct with U @ np.diag(S) @ Vt and test with np.allclose.
import numpy as np
A = np.array([[3.,2.,2.],[2.,3.,-2.]])
U, S, Vt = np.linalg.svd(A, full_matrices=False)
print(S.round(4))
print(np.allclose(U @ np.diag(S) @ Vt, A))| output | expected value |
|---|---|
| S | [5, 3] |
| U shape | (2, 2) |
| Vt shape | (2, 3) |
| reconstruct == A | True |
Worked example
Your turn: build the rank-1 approximation σ₁ u₁ v₁ᵀ and check its spectral-norm error equals σ₂. Say what number you expect before running.
Hint: A1 = S[0] * np.outer(U[:,0], Vt[0]); then np.linalg.norm(A - A1, 2) should equal S[1].
import numpy as np
A = np.array([[3.,2.,2.],[2.,3.,-2.]])
U, S, Vt = np.linalg.svd(A, full_matrices=False)
A1 = S[0] * np.outer(U[:,0], Vt[0])
print(round(np.linalg.norm(A - A1, 2), 4))
print(round(S[1], 4))| quantity | expected value |
|---|---|
| A₁ | [[2.5, 2.5, 0], [2.5, 2.5, 0]] |
| ‖A − A₁‖₂ | 3.0 |
| σ₂ | 3.0 |
Worked example
Your turn: solve Lesson 7's regression with pinv and confirm it matches the normal-equations answer. Predict the coefficients first.
Hint: build X = np.c_[np.ones(5), x], then np.linalg.pinv(X) @ y.
import numpy as np
x = np.array([1.,2.,3.,4.,5.]); y = np.array([2.,4.,5.,4.,5.])
X = np.c_[np.ones(5), x]
print((np.linalg.pinv(X) @ y).round(4))| method | w (expected) |
|---|---|
| pinv(X) @ y | [2.2, 0.6] |
| normal equations (Lesson 7) | [2.2, 0.6] |
| agree? | yes |
Trade off
Comparison matrix
From Milestone 3 — pseudoinverse OLS: every row here is a choice with a cost. Fill the w (expected) column, then say which row you would actually pick and what you give up for it.
| method | w (expected) |
|---|---|
| pinv(X) @ y | [2.2, 0.6] |
| normal equations (Lesson 7) | [2.2, 0.6] |
| agree? | yes |
Concept
import numpy as np
A = np.array([[3.,2.,2.],[2.,3.,-2.]])
U, S, Vt = np.linalg.svd(A, full_matrices=False)
A1 = S[0]*np.outer(U[:,0], Vt[0]) # best rank-1
err = np.linalg.norm(A - A1, 2) # == S[1]
print('sigma =', S.round(4), ' rank-1 err =', round(err,4))
x = np.array([1.,2.,3.,4.,5.]); y = np.array([2.,4.,5.,4.,5.])
X = np.c_[np.ones(5), x]
print('pinv OLS w =', (np.linalg.pinv(X) @ y).round(4))| printed line | value (verified) |
|---|---|
| sigma = [5. 3.] rank-1 err = | 3.0 |
| pinv OLS w = | [2.2 0.6] |
If σ = [5, 3], the rank-1 error is exactly 3 (= σ₂), and pinv recovers [2.2, 0.6] — you've wielded the crown jewel of linear algebra from the math up.
Comparison
Comparison matrix
From The full program: refill the value (verified) column from what you know. The rest of the table is as it appeared.
| printed line | value (verified) |
|---|---|
| sigma = [5. 3.] rank-1 err = | 3.0 |
| pinv OLS w = | [2.2 0.6] |
Concept
Slides closed, out loud: explain (1) why every matrix has an SVD but not every matrix has an eigendecomposition, (2) why the rank-1 error equals σ₂, and (3) what Σ⁺ does to a zero singular value and why.
Stretch: prove U's columns are eigenvectors of AAᵀ and V's of AᵀA from A = UΣVᵀ. SVD returns as PCA (Lesson 19), LoRA low-rank adapters (Week 38), and attention-matrix analysis (Week 30).
Counterexample
Discussion prompt
Stretch: prove U's columns are eigenvectors of AAᵀ and V's of AᵀA from A = UΣVᵀ. SVD returns as PCA (Lesson 19), LoRA low-rank adapters (Week 38), and attention-matrix analysis (Week 30).
That is stated as though it always holds. Do one of two things: produce a case where it fails, or say precisely what rules such a case out. "It just does" is not on the menu.
Hint: Hunt at the extremes first — zero, one, negative, empty, equal. If every extreme survives, the reason they survive is the proof.
Connect it up
Draw it
One page, no notation unless you need it: draw how these connect — Why every matrix needs this · Where singular values come from · The geometry: A vᵢ = σᵢ uᵢ · Compute the SVD by hand · Low-rank approximation · Compression that actually saves. Put an arrow wherever one of them is what makes another possible, and label the arrow with why.
Recap
A = UΣVᵀ (rotate · scale · rotate), and name the shapes of U, Σ, Vσᵢ = √λᵢ(AᵀA) step by step, with V = eigenvectors of AᵀA and U = eigenvectors of AAᵀA vᵢ = σᵢ uᵢ and read σ₁ = ‖A‖₂ as the maximum stretchσₖ₊₁ (spectral) / √Σσᵢ² (Frobenius)A⁺ = VΣ⁺Uᵀ, matching pinv and the normal equations| idea | the one thing to remember |
|---|---|
| existence | every matrix has an SVD — any shape, any rank |
| singular values | √eig(AᵀA), always real and ≥ 0 |
| geometry | A vᵢ = σᵢ uᵢ; σ₁ is the max stretch = ‖A‖₂ |
| low rank | keep the largest σ; error = σₖ₊₁ |
| PCA | right singular vectors of centered X; var = σ²/(n−1) |
| pseudoinverse | A⁺ = VΣ⁺Uᵀ, invert only nonzero σ |
Want this taught 1-on-1? Alexander tutors Machine Learning — $55/session, free consultation.