Linear algebra: the shape of data

You learned linear algebra as procedures on paper — row reduction, cofactor expansion, a list of identities to memorise. None of that survives contact with a dataset. What survives is the geometry: a vector is a point, a matrix is a transformation, a dot product is a measure of agreement, and a rank is how much of the space you can actually reach. This module rebuilds the subject that way, and attaches every idea to the exact line of ML where it does work.

Cold retrieval — answer first, then open the cards

Model answers

1 · Why linear algebra is the language of data

A dataset is a rectangle with a convention attached: one row is one sample, one column is one feature. $n$ customers, $d$ measurements each, and entry $(i,j)$ is feature $j$ of sample $i$:

Once you accept that convention, you can read the whole field as three operations. Almost every formula in machine learning is a sum over rows, a combination of columns, or a change of basis:

Say them in English and the ML appears: the loss averages an error over samples; a prediction is a weighted mix of the features; PCA rotates the data into coordinates where the axes mean something. A "sum over rows" is written $n$ times in a loop or once as a matrix product — and writing it once is what makes it run on a GPU. That is the entire practical difference between a script and a framework.

2 · Vectors: not lists, points with direction

A vector is a list of numbers that belong together as one object — but the useful picture is geometric: a point in space, or an arrow from the origin to that point. Two operations build everything else. Addition places the arrows head-to-tail (the diagonal of the parallelogram they span); scalar multiplication stretches ($c \gt 1$), shrinks ($0 \lt c \lt 1$) or, if $c$ is negative, flips the arrow through the origin. Feed a vector set to both operations and you get linear combinations, and the set of everything you can build is the span:

Two non-parallel vectors in the plane span the whole plane. Three vectors in 3-D that all lie in one plane span only that plane — no combination of them escapes it, no matter how you weight them. Hold that image; it is the whole of the next three sections.

1. A feature vector is a point. A customer with (age, income, clicks) = (34, 52000, 7) is a single point in $\mathbb{R}^3$. "Similar customers" is literally "nearby points" — every clustering, nearest-neighbour and recommendation method is a distance function applied to these points.
2. A linear combination is a reachable output. $f(x) = w^\top x + b$ is a weighted mix of the features; a model that only ever returns sums of its inputs can only ever reach the span of what you gave it.
3. Embeddings are vectors you can do arithmetic on. The classic (approximate) result in word embeddings: $\text{king} - \text{man} + \text{woman} \approx \text{queen}$. That only works because addition and scalar subtraction in the vector space correspond to real relationships in the data.

3 · The dot product: alignment, not just multiplication

The dot product is the one number that turns two vectors into one scalar, and it has a definition that looks mechanical and a meaning that is geometric — they are the same fact:

Take three readings of it, because all three show up in ML:

1. Projection. $\|u\|$ times the component of $v$ along $u$. How much of $v$ points the same way as $u$? Ask a question in one language, measure an answer in another — this is how a linear classifier scores a point against a weight direction.
2. Pairwise feature products. $\sum_i u_i v_i$: feature 1 with feature 1, feature 2 with feature 2… It is one number summarising two aligned feature vectors — the score $w^\top x + b$ of every linear model, the similarity $q \cdot k$ of every attention head.
3. Agreement. Positive when they lean the same way, exactly zero when they are perpendicular, negative when opposed. Orthogonality $\Rightarrow$ dot product zero $\Rightarrow$ "these two features carry no linear relationship" — which is precisely what PCA tries to manufacture.

Worked by hand: $(1,2,3) \cdot (4,5,6) = 4 + 10 + 18 = 32$, with $\|u\| = \sqrt{14} \approx 3.742$ and $\|v\| = \sqrt{77} \approx 8.775$. So $\cos\theta = 32 / \sqrt{14 \cdot 77} = 32/\sqrt{1078} \approx 0.975$ — the two vectors are about $13°$ apart, strongly aligned. Strip the lengths off and you have the similarity score that retrieval actually wants:

4 · Matrices are transformations

Stop reading a matrix as a table of numbers and start reading it as an instruction. Multiplying $A$ by a vector $x$ produces a weighted combination of $A$'s columns, with the entries of $x$ as the weights:

Three shapes to keep straight. The identity $I$ (ones on the diagonal) changes nothing, exactly like the number 1. The transpose $A^\top$ flips rows and columns, $(A^\top)_{ij} = A_{ji}$, and it satisfies $\left(AB\right)^\top = B^\top A^\top$ — note the order reversal. The inverse $A^{-1}$ satisfies $A^{-1}A = AA^{-1} = I$; for a two-sided inverse $A$ must be square ($n \times n$) and its columns must be independent.

Non-square matrices are not broken — they change the dimension of space. A $k \times d$ matrix maps $\mathbb{R}^d \to \mathbb{R}^k$: feed it $d$ numbers, get $k$ numbers back. When $k \lt d$ that is a compression, and it is the shape of every layer that reduces: 784 pixels into 10 hidden units, 512 features into 64, a whole vocabulary into a dense embedding. Reading shapes tells you what an operation does before you read a line of code.

The classical use is the linear system $Ax = b$: find the coefficients that build $b$ out of the columns of $A$. If $b$ lies in the span of $A$'s columns there are solutions (unique if the columns are independent, infinitely many if they are not); if $b$ is outside that span, there are none — only a least-squares compromise, which is exactly what regression computes. Singular matrices are precisely those whose columns do not span the space they live in:

5 · Matrix multiplication, properly

The definition is one line: take row $i$ of the left matrix and column $j$ of the right, dot them.

So $\begin{bmatrix} 1 & 2 \\ 3 & 4 \end{bmatrix}\begin{bmatrix} 5 & 6 \\ 7 & 8 \end{bmatrix}$ has top-left entry $1\cdot5 + 2\cdot7 = 19$, top-right $1\cdot6 + 2\cdot8 = 22$, and the product is $\begin{bmatrix} 19 & 22 \\ 43 & 50 \end{bmatrix}$. The definition is correct but useless for thinking. Use one of these three readings instead — they are all the same operation:

1. Columns. Column $j$ of $AB$ equals $A$ times column $j$ of $B$: each output column is a combination of $A$'s columns, mixed with the weights sitting in that column of $B$.
2. Rows. Row $i$ of $AB$ equals row $i$ of $A$ times $B$: each output row is a combination of $B$'s rows. If the rows of $B$ are data samples, row $i$ of $A$ decides how to blend them.
3. Composition. $(AB)x = A(Bx)$: apply $B$ first, then $A$. Two layers stacked in one matrix. This is why order matters and why $AB \ne BA$ in general.

Why the obsession with this one operation? Cost and hardware. Multiplying two square $1000 \times 1000$ matrices costs about $n^3 = 10^9$ multiply–add operations; training a deep network runs an enormous number of such products, and the overwhelming majority of wall-clock training time is spent inside general matrix multiply kernels. The other reason is batching. With $X$ an $n \times d$ batch of samples and $W$ a $d \times k$ weight matrix:

One giant multiply instead of $n$ tiny ones. GPUs have thousands of cheap cores and very high memory bandwidth, which is exactly the shape of $XW$ — and nothing else. Before BLAS-style batched matrix multiplication, a neural net trained on a handful of examples per second; after it, the same code runs on billions of tokens, because the arithmetic did not change, only its packaging.

6 · Rank, span and collinearity

The rank of a matrix is the number of linearly independent columns — equivalently, the dimension of the span they generate, equivalently the number of independent rows (rows and columns always agree). It is the honest size of the information in your data:

A rank deficiency is not an abstract defect — it is a duplicate or a derived feature hiding in your dataset. Two columns that are proportional (metres and centimetres), or one column that is the sum of others, and your matrix has collapsed. Three places you will meet it immediately:

• The dummy-variable trap. One-hot encode a 4-category feature into 4 indicator columns and add an intercept column of ones — the ones column is exactly the sum of the four indicators. $X^\top X$ is singular, the regression has no unique solution. Fix: drop one level, or drop the intercept.
• Multicollinearity. When two features are nearly (not exactly) collinear the model can still fit, but the coefficients become unstable and huge — tiny changes in the data swing them wildly. Ridge regression adds $\lambda I$ precisely to put a floor under those directions (Section 9 shows why that works).
• Why PCA exists. Real datasets are strongly correlated, so their effective rank is far below $d$: a 1 000-feature sensor rig driven by 8 physical quantities has rank 8 in the ways that matter. Low rank is not a problem to fix — it is structure to exploit, and it is what Section 11 and Section 12 harvest.

7 · Determinant: volume scaling and invertibility

The determinant is a single number that answers: if I feed the unit square (or cube, or hypercube) through this matrix, how does the volume change? Not the direction, just the factor. The two-dimensional case is worth memorising because it carries the whole intuition, $\det\begin{bmatrix} a & b \\ c & d \end{bmatrix} = ad - bc$: it is the signed area of the parallelogram spanned by the two columns.

Read the three cases as a rule. Absolute value = scaling factor: $|\det A| = 6$ means every volume in the space is multiplied by 6. Sign = orientation: $\det A \lt 0$ means space got flipped inside out, like a reflection — areas still grow, but handedness is reversed. Zero = collapse: the matrix sends a full-dimensional chunk of space onto a plane or a line, volume becomes 0, and no inverse can undo that (you cannot un-flatten a pancake). For square matrices, all three collapse into one equivalence:

Two caveats that prevent the usual damage. First, the determinant is only defined for square matrices — a $4 \times 7$ matrix has no determinant, so "the determinant is non-zero, therefore it is invertible" is meaningless until you know the shape. Second, do not compute it by cofactor expansion: computers use elimination (LU), and in practice you almost never compute a determinant directly — you ask a library whether the matrix is singular. What you do use is its meaning: the determinant of a covariance matrix controls the volume of the uncertainty ellipsoid in a Gaussian density (it appears as $\log|\Sigma|$ in Module 04), and it equals the product of the eigenvalues, $\det A = \prod_i \lambda_i$ — the volume scaling is just all the per-axis stretch factors multiplied together. That is the bridge to Section 10.

8 · Norms and distances

A norm is a length measurement — a function that turns a vector into a non-negative number, zero only for the zero vector, and that grows predictably when you scale the vector. There is a family of them, indexed by $p$:

Worked: $x = (3,-4)$ gives $\|x\|_2 = 5$, $\|x\|_1 = 7$, $\|x\|_\infty = 4$. To see why the three differ, draw the set of vectors with norm exactly 1 — the "unit circle" of each norm:

• $p = 2$: a circle. Smooth, no corners, every direction is equivalent. This is the length we mean when we say "length".
• $p = 1$: a diamond. A square stood on its corner, with its four vertices on the axes — points like $(1,0)$ and $(0,1)$. Every point on the diamond's edges is on the way to a vertex that lies on an axis.
• $p = \infty$: a square. Flat sides at $\pm 1$; the norm just watches the worst coordinate. The three nest: $\|x\|_\infty \le \|x\|_2 \le \|x\|_1$.

That diamond is why $L_1$ does something $L_2$ cannot. When you minimise a smooth objective over the unit ball, the answer tends to land where the constraint's corners are — and the diamond's corners sit exactly on the axes, where some coordinate is exactly zero. The circle has no corners, so its solutions land generically in the middle of nowhere: every coordinate survives, just a bit smaller. This is the entire mechanism of Lasso sparsity versus the gentle shrinkage of $L_2$ weight decay — you will derive it again in Modules 07 and 08, and it reappears as feature selection and as pruning.

For distances between data points, two metrics do nearly all the work. Euclidean distance $\|u - v\|_2$ measures how far apart two points are — k-means minimises exactly this (the "k" is not negotiable: change the norm and you have a different algorithm), and it is right when magnitude is meaningful: pixels, prices, sensor readings. Cosine distance, $1 - \cos(u,v)$, measures the angle and ignores the length — right for text and embeddings, where document length or user activity inflates the norm without adding meaning. The two are not rivals either — for unit vectors the whole distance collapses to the angle:

9 · Special matrices you will meet in ML

Four matrix types do almost all the work. Each has a one-line definition, a picture, and a specific place it shows up:

The last row deserves its working shown, because it answers a question you will hit the moment you add regularisation. Take any matrix $X$ (possibly rank-deficient) and any $\lambda \gt 0$:

Positive definiteness is "the function has a bowl shape": $f(x) = \tfrac{1}{2}x^\top A x$ with $A$ positive definite decreases in every direction you move away from zero, so zero is the unique minimum. That is why convexity arguments in Modules 03 and 07 keep asking whether a Hessian is positive definite, and why a covariance matrix must be at least positive semi-definite — a direction with zero variance is a direction the data never varies in.

10 · Eigenvalues and eigenvectors: the directions that do not turn

Feed a matrix a generic vector and it comes back rotated off its own span — no scalar multiple of the input can equal the output. Eigenvectors are the exceptions, the rare directions where the matrix behaves like plain arithmetic:

The length of $v$ is irrelevant — if $Av = \lambda v$ then $A(cv) = \lambda(cv)$ for any scalar $c$, so every line through an eigenvector is full of eigenvectors. Only the direction and the factor $\lambda$ carry meaning. The eigenvector is not "the longest" or "the most important" vector in any sense; importance, if you need it, comes from the associated $\lambda$ and from what the matrix is.

A quick intuition check: a rotation matrix in 2-D. If $A$ rotates every vector by $30°$, is there any direction it leaves on its own span? No — and the algebra agrees. The characteristic equation $\det(A - \lambda I) = 0$ for $R(\theta)$ expands to

For $\theta \ne 0, \pi$ there are no real eigenvectors: the roots are complex. That is a genuine diagnostic — complex eigenvalues mean the matrix contains rotation. The general recipe is the same for any $n \times n$ matrix: build $\det(A - \lambda I)$, a degree-$n$ polynomial, and its roots are the eigenvalues (counted with multiplicity, possibly complex). Work the trivial case to keep it grounded: for $A = \begin{bmatrix} 2 & 0 \\ 0 & 3 \end{bmatrix}$ the equation is $(2-\lambda)(3-\lambda) = 0$, giving $\lambda = 2, 3$ with eigenvectors $e_1, e_2$ — the axes, because a diagonal matrix never mixes them.

The payoff of finding them all is that the matrix becomes trivial in that basis. Stack the eigenvectors as columns of $P$ and put the eigenvalues on the diagonal of $D$, and you get a diagonalisation:

For $A = \operatorname{diag}(2,3)$, $A^{10} = \operatorname{diag}(2^{10}, 3^{10}) = \operatorname{diag}(1024, 59049)$ — one line instead of ten matrix multiplications. This is exactly how long-run behaviour is analysed. In a Markov chain with row-stochastic transition matrix $P$ (rows sum to 1, entries are probabilities), the distribution after many steps is $\pi_0 P^k$, and the steady state is the distribution that applying $P$ does not change:

Concretely, $P = \begin{bmatrix} 0.9 & 0.1 \\ 0.4 & 0.6 \end{bmatrix}$ has steady state $\pi = (0.8, 0.2)$: check $0.8 \times 0.9 + 0.2 \times 0.4 = 0.8$ and $0.8 \times 0.1 + 0.2 \times 0.6 = 0.2$. PageRank is the same idea with 10 million states, plus a damping trick to keep the matrix well behaved.

11 · Covariance matrix and PCA, derived

The covariance matrix collects how every pair of features co-varies across samples. Centre the data first ($X_c = X - \bar{x}$, each column minus its mean), then:

Now the question PCA answers. Take a unit direction $u$ (say, "the direction from the origin to (0.83, 0.55)") and project every sample onto it — that gives one number per sample, $X_c u$. How variable is that projection? Its variance is a dot-product identity away from $\Sigma$:

So PCA is an optimisation: pick the unit direction that maximises $u^\top \Sigma u$, subject to $u^\top u = 1$. Solve it with a Lagrange multiplier — watch how the constraint folds into the objective (this exact move returns in Module 07):

That is the whole derivation: maximise variance under a unit-length constraint and you get an eigenvector. Order the eigenvalues $\lambda_1 \ge \lambda_2 \ge \dots$, and the corresponding eigenvectors are the principal components — direction 1 holds the most variance, direction 2 the most of what is left, and they are orthogonal (because $\Sigma$ is symmetric). Worked example: for $\Sigma = \begin{bmatrix} 2 & 1.2 \\ 1.2 & 1 \end{bmatrix}$ the eigenvalues are exactly $2.8$ and $0.2$ (trace $3$, determinant $0.56$), and the first component is $u_1 = (3,2)/\sqrt{13} \approx (0.832,\, 0.555)$, explaining $2.8/3 = 93.3\%$ of the variance.

12 · SVD: the fundamental theorem of linear algebra

The eigendecomposition needs a square matrix. The singular value decomposition does not: every matrix, any shape, factors into an orthogonal basis on each side and a diagonal stretch in the middle:

The geometric story is three moves: $V^\top$ rotates the input so the coordinate axes line up with the data's natural directions; $\Sigma$ scales along those axes, stretching by $\sigma_1, \sigma_2, \dots$ and sending the extra dimensions to zero; $U$ rotates the result into the output's frame. Every matrix is a rotation–scaling–rotation (when the scaling sends dimensions to zero, it is also a collapse — which is how rank deficiency shows up). Written as a sum of rank-1 pieces, the structure is explicit:

Chopping the tail is not a hack; it is provably optimal. The Eckart–Young theorem states that the matrix built from the top $k$ singular values is the closest rank-$k$ matrix to $A$ in the Frobenius norm — no other matrix of rank $k$ does better:

Three applications, in decreasing order of how often you will meet them:

1. Image compression. A greyscale image is a matrix; store only the top $k$ singular triplets — $k(m + n + 1)$ numbers instead of $mn$ — and reconstruct. Rank 20 on a 128×128 image is a 3× saving with a few percent error (you will run this below).
2. Latent factors in recommender systems. A users × items matrix with a few hidden tastes behind it is low-rank: "this user is a 0.9 action fan", "this film is 0.8 action" — the product predicts the rating. Missing entries are filled from the low-rank reconstruction, which is exactly the exercise below.
3. PCA itself. Centre the data and take its SVD, $X_c = USV^\top$; then $\Sigma = \frac{1}{n-1}X_c^\top X_c = \frac{1}{n-1} V S^2 V^\top$, so the right singular vectors are the principal components and the explained variances are $s_i^2/(n-1)$. One algorithm, two entrances.

13 · The ML map

Close the loop: every concept in this module, and the machine it drives. Keep this table as the index — when an algorithm seems opaque, find its row here and name the load-bearing math before reading the code.

Four blocks: alignment, transformation, PCA, compression

Every output below is verified — run the code and check the printed numbers against the comments. Then break something: scale a vector, change the angle, drop a component, and predict what moves before you run it.

Block 1 — length bias. The same query against four candidates, scored two ways: raw dot product and cosine. Watch the rankings disagree.

import numpy as np

u = np.array([1., 2., 3.])          # the query embedding
v = np.array([4., 5., 6.])          # a candidate: long, and roughly the same direction
s = 0.1 * u                         # identical direction to u, 10x shorter
t = np.array([10., 0., 0.])         # long, but almost perpendicular to u

def dot(a, b):    return float(a @ b)
def cosine(a, b): return dot(a, b) / (np.linalg.norm(a) * np.linalg.norm(b))

print("u . v     =", dot(u, v))                 # 32.0
print("cos(u, v) =", round(cosine(u, v), 4))    # 0.9746  -> about 12.9 degrees apart
print()
print(f"{'candidate':10s} {'raw dot':>9s} {'cosine':>9s}")
for name, w in [("v", v), ("10 * v", 10 * v), ("s = 0.1 u", s), ("t", t)]:
    print(f"{name:10s} {dot(u, w):9.2f} {cosine(u, w):9.4f}")

# raw dot ranks: 10*v (320) > v (32) > t (10) > s (1.4)  -- size wins
# cosine ranks:    s (1.000) > v = 10*v (0.9746) > t (0.2673)  -- direction wins
# s and t swap places: the raw dot product sends the WRONG candidate to the top

The raw score and the cosine score disagree on which candidate is best. Neither is "wrong": the raw dot product measures alignment times magnitude, cosine measures alignment alone. If your ranking task cares about direction, you must normalise before you sort.

Block 2 — matrices as transformations. An orthogonal matrix (rotation) must preserve every length; a rank-deficient matrix must collapse space onto a plane and refuse to invert.

import numpy as np

# --- an orthogonal matrix: rotation by 30 degrees ---
th = np.pi / 6
R = np.array([[np.cos(th), -np.sin(th)],
              [np.sin(th),  np.cos(th)]])
print("R @ R.T = I :", np.allclose(R @ R.T, np.eye(2)))   # True -> inverse is the transpose
print("det R   =", round(float(np.linalg.det(R)), 12))     # 1.0   -> area preserved

pts = np.array([[1., 0.], [3., 4.], [-2., 5.]])
before = np.linalg.norm(pts, axis=1)
after  = np.linalg.norm(pts @ R.T, axis=1)   # rows are points, so rotate with R.T
print("lengths before:", np.round(before, 6)) # [1.  5.  5.385165]
print("lengths after :", np.round(after, 6))  # identical
print("max change    =", np.abs(before - after).max())   # 0.0

# --- a rank-deficient matrix: 3-D points collapse onto a plane ---
A = np.array([[1., 2., 3.],
              [2., 4., 6.],
              [1., 0., 1.]])                  # row 2 = 2 * row 1 -> rank 2
print("rank(A) =", np.linalg.matrix_rank(A))  # 2
print("det(A)  =", np.linalg.det(A))          # 0.0

y = A @ np.array([1., 2., 3.])
print("A x =", y, "-> 2nd coord is 2 x 1st:", np.isclose(y[1], 2 * y[0]))   # [14. 28.  4.] True
rand3 = np.random.default_rng(1).normal(size=(4, 3))
images = rand3 @ A.T
print("all images on the plane y2 = 2*y1:",
      np.allclose(images[:, 1], 2 * images[:, 0]))            # True
try:
    np.linalg.inv(A)
except np.linalg.LinAlgError as err:
    print("np.linalg.inv ->", err)                           # Singular matrix

Three facts in one run: $R^\top R = I$ with $\det R = 1$ (a rotation: no length is gained or lost), rank 2 with determinant 0 (a collapse — every image lands on the plane $y_2 = 2y_1$, a two-dimensional sliver inside 3-D), and np.linalg.inv refusing to work. Singular matrices are not an inconvenience libraries hide; they are a geometric fact you can watch.

Block 3 — PCA from scratch, then checked against scikit-learn. The covariance matrix is built with known eigenvalues ($2.8$ and $0.2$) so you can compare the answer to the truth.

import numpy as np
rng = np.random.default_rng(0)

# population covariance chosen so the true eigenvalues are exactly 2.8 and 0.2
Sigma_true = np.array([[2.0, 1.2],
                       [1.2, 1.0]])
X = rng.multivariate_normal([1.0, -2.0], Sigma_true, size=500)

# 1) centre: PCA is about spread, not location
Xc = X - X.mean(axis=0)
# 2) sample covariance, shape (d, d)
Sigma = Xc.T @ Xc / (len(Xc) - 1)
# 3) eigh is for symmetric matrices: ascending eigenvalues, orthonormal eigenvectors
lam, V = np.linalg.eigh(Sigma)
lam, V = lam[::-1], V[:, ::-1]                     # largest eigenvalue first

print("eigenvalues  :", np.round(lam, 3))          # [2.719 0.188]  (truth: 2.8, 0.2)
print("PC1          :", np.round(V[:, 0], 3))       # [0.833 0.555] up to sign (truth: (3,2)/sqrt(13))
print("explained    :", np.round(lam / lam.sum(), 3))  # [0.935 0.065]
print("eigen eq holds:", np.allclose(Sigma @ V[:, 0], lam[0] * V[:, 0]))   # True

# 4) project onto PC1, then rebuild a rank-1 version of the data
z = Xc @ V[:, 0]                                   # one score per sample
X_hat = np.outer(z, V[:, 0])                       # n x 2 reconstruction
print("reconstruction MSE:", round(float(np.mean((Xc - X_hat) ** 2)), 4))  # 0.0936 ~ lam2/2

# 5) sanity check: same direction as scikit-learn?
from sklearn.decomposition import PCA
p = PCA(n_components=1).fit(X)
print("sklearn PC1  :", np.round(p.components_[0], 3))       # [0.833 0.553]
print("same direction:", abs(float(p.components_[0] @ V[:, 0])) > 0.999)   # True

Read the output against the theory: the sample eigenvalues sit close to the true $2.8$ and $0.2$; PC1 is $\pm(0.833, 0.555)$, which is $(3,2)/\sqrt{13}$ — the direction of maximal spread; the explained ratio $0.935$ matches $2.8/3 = 0.933$; and the reconstruction error $0.094$ is almost exactly half of the discarded eigenvalue $0.188/2$ (half, because the MSE is averaged over both coordinates). The sign of an eigenvector is arbitrary — only the line matters.

Block 4 — truncated SVD on an image. Build a smooth synthetic texture (natural images have this frequency spectrum), then reconstruct it from the top 5, 20 and 100 singular values.

import numpy as np

# a synthetic 128x128 grayscale texture: smooth fields compress, noise does not
rng = np.random.default_rng(4)
n = 128
F = np.fft.fft2(rng.normal(size=(n, n)))
fy, fx = np.fft.fftfreq(n)[:, None], np.fft.fftfreq(n)[None, :]
r = np.hypot(fx, fy); r[0, 0] = 1e-6               # divide by frequency -> cloudy, low-rank-ish field
img = np.real(np.fft.ifft2(F / r ** 1.6))
img = (img - img.min()) / (img.max() - img.min())

U, sv, Vt = np.linalg.svd(img, full_matrices=False)
for k in (5, 20, 100):
    low = (U[:, :k] * sv[:k]) @ Vt[:k]             # best rank-k approximation (Eckart-Young)
    err = np.linalg.norm(img - low) / np.linalg.norm(img)
    cost = n * n / (k * (n + n + 1))               # numbers stored now vs. the full image
    print(f"rank {k:3d}   rel. error {err:.4f}   compression {cost:5.2f}x")
# rank   5   rel. error 0.0989   compression 12.75x
# rank  20   rel. error 0.0359   compression  3.19x
# rank 100   rel. error 0.0025   compression  0.64x   (worse than useless: keep the pixels)

The trade curve is the lesson: rank 5 costs 10% of the energy for a 12.75× saving, rank 20 costs 3.6% for 3.19×, and rank 100 is nearly lossless but costs more storage than the original (compression below 1× means you should not have bothered). Compression only pays when the singular values decay — which is a statement about your data, not about your algorithm.

import numpy as np

R = np.array([
    [5, 5, 4, 1, 1, np.nan],
    [5, 4, 5, 1, np.nan, 1],
    [4, 5, np.nan, 3, 3, np.nan],
    [np.nan, 1, 2, 5, 5, 4],
    [1, np.nan, 1, 5, 4, 5],
], float)                              # 5 users x 6 films, NaN = not yet seen
mask = np.isnan(R)
missing = list(zip(*(a.tolist() for a in np.where(mask))))

M = R.copy()
M[mask] = np.take(np.nanmean(M, axis=0), np.where(mask)[1])   # start from column means
for _ in range(200):
    U, s, Vt = np.linalg.svd(M, full_matrices=False)
    low = (U[:, :2] * s[:2]) @ Vt[:2]                         # rank-2 model of the ratings
    new = low.copy()
    new[~mask] = R[~mask]                                     # observed ratings are facts: never move them
    if np.allclose(new, M, atol=1e-10):
        M = new
        break
    M = new

print("singular values:", np.round(s, 2))   # [16.81 8.96 1.22 1.04 0.63] -> big gap after 2
for u, j in missing:
    print(f"user {u + 1} -> film {j + 1}:  {low[u, j]:.2f}")

# user 1 -> film 6:  0.66     user 3 -> film 6:  2.61
# user 2 -> film 5:  1.35     user 4 -> film 1:  1.43
# user 3 -> film 3:  4.53     user 5 -> film 2:  0.54

The singular values fall off a cliff — $16.81,\ 8.96$, then $1.22$ — so two hidden tastes really do explain this matrix, and rank 2 is a property of the data rather than a guess. The recommendation is unambiguous: user 3 → film 3 at 4.53 is the only missing pair the model likes; every other unseen combination is a rejection, and two of them (0.66, 0.54) are so negative the model has left the 1–5 scale entirely — clip predictions to the valid range in anything you ship. This is exactly how matrix completion works on real recommender systems, only with $10^6 \times 10^5$ instead of $5 \times 6$.

Multiple choice — instant feedback

Open questions — answer aloud, without notes

Re-answer these in 1 day, then 3 days, then 7 days — from memory, in writing. Ten minutes of recall on day three beats an hour of re-reading today.