6-month learning plan · · session 18 / 13018 / 130

mathbeginner 55m read

S018 · Linear Algebra II — Matrices, Transforms, Eigenvalues

Matrices as linear transformations, matrix multiplication as function composition, and eigenvalues as the axes a transform respects — the shape of every neural net.

📐MathM02 · Math Foundations· Session 018 of 130 90 min

🎯 Read a matrix as a linear transformation, compose transformations by multiplying matrices, and know what eigenvalues/eigenvectors tell you about the transform's behaviour.

Why this session exists

A neural network is nothing but stacked matrix multiplications and non-linearities. A rendered 3-D scene is a sequence of matrix transformations. PCA, SVD, PageRank, JPEG, quantum states — all matrices. The reason matrices show up everywhere is that they are the only way to write down a linear transformation of a vector space. Get past "matrix = grid of numbers" to "matrix = function that stretches, rotates, and squashes space" and every downstream ML paper stops feeling like a wall of Greek letters.

You will be able to

Read a matrix as a function that maps input vectors to output vectors.
Multiply two matrices by hand for 2×2 and 3×3 cases (and know why row × column).
Recognise identity, transpose, inverse, and diagonal — and say what each does geometrically.
Explain eigenvalues and eigenvectors as the ‘axes the transform respects’ with a picture.
Use NumPy to multiply matrices, solve linear systems, and find eigenvalues.

Prerequisites

S017 · Vectors — dot product IS the atom of matrix multiplication.
Comfortable with element-wise math and Σ notation.

Pre-read · skim before the session (35 min total)

💡 Intuition 11 min

Linear transformations and matrices · Chapter 3, Essence of Linear Algebra

3Blue1Brown

The single most important 11 minutes of ML prep on the internet. Watch, then watch again.

🔬 Deep dive 17 min

Eigenvectors and eigenvalues · Chapter 14, Essence of Linear Algebra

3Blue1Brown

Eigenvalues finally make sense as ‘axes the transform doesn't rotate’. Beats every textbook.

🛠️ Hands-on 40 min

MIT 18.06 — Introduction to Linear Algebra · Lecture 1

MIT OpenCourseWare

Gilbert Strang's opening lecture. Watch parts you find confusing again in his other lectures.

📖 Reading 25 min

Mathematics for Machine Learning (Ch. 2.3–2.7 and 4.1–4.3)

Deisenroth, Faisal, Ong

Rigorous but readable. Chapter 4 is the eigen chapter — bookmark it.

📖 Reading 10 min

NumPy linalg reference

numpy.org

Every function you'll use: matmul, solve, inv, eig, svd. Read the top summaries.

📖 Reading 8 min

Matrix multiplication

Wikipedia

The definition and geometric intuition on one page. The animation in the article is gold.

(a) Intuition · 5 min

A matrix is a function that transforms space

🌍 Real world

Imagine a rubber sheet with a grid drawn on it. You pinch two corners and stretch it — the grid deforms. Straight lines stay straight, parallel lines stay parallel, the origin stays put. That constrained stretch is a linear transformation, and it is entirely described by where the two basis vectors [1,0] and [0,1] end up. Those two ending points are the two columns of a 2×2 matrix.

Multiply the matrix by any point on the original sheet and it tells you where that point moved to. Rotate the sheet 90°? Matrix. Stretch it 2× vertically? Matrix. Shear it? Matrix.

💻 Code world

In code, a matrix is a 2-D array. But its meaning is a function. A @ x in NumPy takes vector x and returns Ax — the transformed vector. Every operation you'll use — matrix-matrix multiply, inverse, eigendecomposition — is a fact about that function.

The columns of A are literally "where the basis vectors go." That's the one sentence that unlocks the whole subject.

The three views of a matrix

A matrix is always all three

A grid of numbers — m rows × n columns. This is the storage view.
A linear function ℝⁿ → ℝᵐ — takes a length-n vector, returns a length-m vector. This is the transformation view.
A composition of column vectors — the k-th column of A is where the k-th basis vector of the input space lands. This is the ‘how do I build one?’ view.

A quick history

200 BCE
The Nine Chapters (China)
Chinese mathematicians solve simultaneous linear equations with a proto-matrix method 2000 years before Europe.
1858
Cayley · Memoir on the Theory of Matrices
Arthur Cayley writes down matrix multiplication and inverse. The name ‘matrix’ was Sylvester's.
1929
von Neumann · quantum mechanics as matrices
John von Neumann formulates QM as linear operators on Hilbert space. Matrices become central to physics.
1947
SVD · Beltrami/Jordan much earlier, but computable now
The Singular Value Decomposition — the single most-used matrix factorisation in ML.
1998
PageRank · Brin & Page
The internet ranked by the principal eigenvector of a hyperlink matrix. Google is born.
2017+
Transformers · attention as matrix multiplies
GPT/BERT do trillions of matrix multiplications per training step. GPUs exist because matrices do.

(b) Visual walkthrough · 15 min

A matrix transforms space

The columns of A are exactly where the basis vectors land. Every other vector's destination is a weighted combination of those two column-destinations.

The matrix multiplication formula, unpacked

Check shapes

(m × k) · (k × n) → (m × n). Inner k must match; outer m and n are the output shape.

Row × column

Entry (i, j) of the result = dot product of row i of A with column j of B.

Result is a matrix of dot products

Every cell of the product is a scalar; it's the ‘answer’ to ‘what does row i of A × column j of B evaluate to?’

Alternative view · composition

A · B is the function ‘first apply B, then apply A.’ Order matters and reads right-to-left.

The special matrices worth memorising

Six matrices you'll see every day

Identity Iₙ

1s on the diagonal, 0s elsewhere. A · I = I · A = A. The ‘do nothing’ transform.

identity

Diagonal D

Non-zero only on the diagonal. Multiplying by D scales each coordinate independently. O(n) instead of O(n²).

diagonal

Transpose Aᵀ

Swap rows and columns. Preserves dot products in a sense: (Ax)·y = x·(Aᵀy). Central to gradient formulas.

transpose

Inverse A⁻¹

The transform that undoes A. A · A⁻¹ = I. Only exists if A is square and det(A) ≠ 0.

inverse

Orthogonal Q

Rows and columns are unit vectors, mutually perpendicular. Q · Qᵀ = I. Represents rotations/reflections — preserves lengths and angles.

orthogonal

Symmetric S

S = Sᵀ. All eigenvalues are real; eigenvectors are orthogonal. Covariance matrices are symmetric — that's why PCA works.

symmetric

Solving Ax = b: three interpretations, one answer

Row picture

Intersection of hyperplanes

Each row of A is a linear equation
Each equation is a hyperplane in ℝⁿ
Solution x is where they all meet
Sensitive to picture-blowing-up in high dims

Column picture

Linear combination of columns

b is a linear combination of columns of A
Solve for the weights
Cleaner geometric intuition
Preferred view for ML

Function picture

Undo the transform

A maps x → b
We want x such that A applies to it gives b
x = A⁻¹ b (when A is invertible)
Numerical: `np.linalg.solve(A, b)` (never invert!)

Eigenvalues and eigenvectors — geometric meaning

Eigenvector: a vector whose direction A does NOT change — A only stretches or shrinks it. Eigenvalue λ: the stretch factor.

Common misconception

✗ What most people think

"A matrix is a 2D array of numbers — a grid, like a spreadsheet or a table. Matrix multiplication is a slightly odd rule you memorise: rows times columns."

✓ What is actually true

A matrix is a linear transformation written in coordinates. Its columns are where the basis vectors land. Matrix multiplication is function composition — AB means "do B, then A" — and every property that looks arbitrary (why the inner dimensions must match, why it isn't commutative, why the identity looks like that) falls out of that single fact.

Why the myth is so sticky

Because the grid model is sufficient for storage and for the mechanics of computing a product, and because in data work you genuinely do use matrices as tables of records. The model fails the moment you need to reason rather than compute. Under the grid view, AB ≠ BA is a strange quirk to remember; under the transformation view it is obvious, because rotating-then-stretching is visibly not the same as stretching-then-rotating. Similarly "singular" sounds like a numerical accident, when it actually means the transformation collapsed a dimension and destroyed information — which is why it cannot be inverted, and why a near-singular matrix makes a regression's coefficients wildly unstable. Everything mysterious about linear algebra becomes mechanical once you stop seeing grids.

Prove it to yourself

The columns literally are the images of the basis vectors — check it:

import numpy as np
A = np.array([[0., -1.], [1., 0.]])     # 90 degree rotation
print(A @ np.array([1., 0.]))            # [0, 1]  <- column 1 of A
print(A @ np.array([0., 1.]))            # [-1, 0] <- column 2 of A

B = np.array([[2., 0.], [0., 1.]])       # stretch x by 2
print(A @ B)
print(B @ A)                             # different: order is composition order

From first principles

Start with the question

Why does a matrix fail to be invertible exactly when its determinant is zero? These look like unrelated facts — one about solving equations, one about a strange alternating sum.

1
A matrix maps input vectors to output vectors. Inverting it means recovering the unique input that produced a given output.
forced by · an inverse function must be well defined, so each output needs exactly one preimage
2
The determinant measures how the transformation scales volume: the unit cube maps to a parallelepiped whose signed volume is det(A).
forced by · linearity means the volume scaling factor is the same everywhere, so one number suffices
3
If det(A) = 0, the output volume is zero — the entire input space has been flattened onto a lower-dimensional subspace (a plane, a line, a point).
forced by · only a degenerate shape has zero volume, which means the columns became linearly dependent
4
Flattening is many-to-one: infinitely many distinct inputs map to the same output, because an entire direction was crushed to zero.
forced by · reducing dimension cannot be injective — that is the pigeonhole argument in continuous form
5
Therefore no inverse function can exist. There is no rule that recovers which of infinitely many preimages you started from; the information was genuinely destroyed, not merely hidden.
forced by · an inverse must map each output back to a single input, and there isn't one

⇒ Therefore

Therefore "determinant zero" and "not invertible" are the same statement viewed geometrically and algebraically: the transformation lost a dimension. And det(AB) = det(A)det(B) stops being a formula to memorise — composing two transformations multiplies their volume scalings, which is the only thing it could possibly do.

And note what this predicts, and it is the practically important part: a determinant near zero means nearly flattened, so a tiny change in output corresponds to a huge change in input. That is ill-conditioning, and it is exactly why perfectly collinear features make X'X singular and near-collinear features make regression coefficients enormous, unstable, and sign-flipping between refits. It also predicts the fix: ridge regularisation adds λI, pushing the eigenvalues away from zero and restoring invertibility by construction.

Mental modelMatrices are verbs, not tables

Read every matrix as an action on space. Its columns tell you where the basis vectors go, and because the transformation is linear, that's enough to determine where everything goes. Multiplying Ax is applying the verb to a vector; multiplying AB is chaining two verbs, applied right to left.

Once matrices are verbs, the vocabulary decodes itself: the identity is "do nothing", the inverse is "undo", the determinant is "how much does this squash or expand", singular is "collapsed a dimension", transpose swaps the roles of input and output space, and eigenvectors are the directions the verb only stretches without turning.

Inner dimensions must match because the output space of the first transformation must be the input space of the second. Shape errors are type errors about which space you're in.
AB ≠ BA because composition order matters. This is not a quirk; it is why the order of operations in a neural network layer stack is meaningful.
Never compute an explicit inverse to solve Ax = b. Use a solver (np.linalg.solve, or a decomposition) — it's faster and far more numerically stable, since forming the inverse amplifies error.
Eigen/SVD decompositions answer "what does this transformation do, in its own natural coordinates?" — which is why SVD underlies PCA, low-rank approximation, and recommender factorisation alike.

🔔 Fires when you see

Fire this model the moment you see: a shape mismatch error · LinAlgError: Singular matrix · regression coefficients that are huge or flip sign between runs · a condition number in a warning · PCA, SVD, or embedding factorisation · a transformer's attention as QK'.

The tradeoff

You need to solve a large least-squares problem Ax ≈ b. Normal equations, QR/SVD decomposition, or iterative methods?

Normal equations (A'A x = A'b)

+ you gain simplest to derive and implement; A'A is small when features are far fewer than rows, and it can be accumulated in one streaming pass over the data

− you pay forming A'A squares the condition number, so it roughly halves the precision you have to work with — near-collinear features that QR handles fine will produce garbage here

pick when features are few, well-conditioned, and you need a single pass over data that doesn't fit in memory

QR or SVD decomposition

+ you gain numerically stable — it works directly on A, so the condition number is never squared; SVD additionally reveals rank and provides the pseudo-inverse for genuinely rank-deficient problems

− you pay needs the full matrix (or a blocked variant), higher cost than the normal equations, and SVD in particular is expensive on very large or very wide matrices

pick when accuracy matters, the matrix might be ill-conditioned, or you need to diagnose rank deficiency rather than merely survive it — the default for any serious numerical work

Iterative methods (conjugate gradient, LSQR, SGD)

+ you gain never forms or factors the matrix — only needs matrix–vector products, so it exploits sparsity and scales to problems that cannot be stored densely; you can stop early at the accuracy you need

− you pay convergence depends on conditioning and often requires a preconditioner to be practical; the answer is approximate with an error you must monitor; more knobs to get wrong

pick when the matrix is huge and sparse, or exists only as an operator, or an approximate solution reached quickly beats an exact one reached slowly — which is the entire premise of training by gradient descent

What a senior engineer actually does

Use a library's least-squares routine (which will pick QR or SVD) unless scale forces you elsewhere, and check the condition number before trusting any coefficient. The recurring production mistake is not choosing the wrong solver — it's reporting coefficients from an ill-conditioned fit as though they were meaningful, when collinearity has made their individual values essentially arbitrary while the predictions remain fine.

That distinction is worth holding onto: ill-conditioning corrupts interpretation long before it corrupts prediction. A model can score well and still have coefficients you must not explain to a stakeholder.

(c) Hands-on · 25 min

You'll implement matrix multiplication from scratch, compare to NumPy, solve a real linear system, and compute eigenvalues to find principal directions of a 2-D dataset. Save as matrices_demo.py.

Run it:

uv run --with numpy python matrices_demo.py

Expected output:

What each block does

Anatomy of the code

matvec + matmul from scratch

Nested comprehensions make the formula obvious: each output entry is a dot product. Great teaching version; ~1000× slower than NumPy for real sizes.

primitives

compare_to_numpy

Differential test: for 200 random shape/value pairs, our result and NumPy's must agree. Catches off-by-one, wrong-axis bugs immediately.

verify

demo_composition

R first (rotate 90°), then S (scale x). Composition R∘S is the same as multiplying the matrices in reverse order. Order matters — this is where non-commutativity bites.

compose

np.linalg.solve(A, b)

NEVER compute A⁻¹ and multiply. `solve` uses LU decomposition — numerically stable and O(n³). Explicit inverse is worse in every way.

solve

np.linalg.eigh(cov)

`eigh` is for symmetric matrices — faster and always returns real eigenvalues. Covariance is always symmetric, so this is the right tool.

eig

PCA in 5 lines

Centre the data → covariance matrix → eigendecomposition → sort by eigenvalue → largest eigenvector is the principal direction. That IS PCA.

pca

Try itFeel that rotation is orthogonal

import numpy as np
R = np.array([[0.0, -1.0], [1.0, 0.0]])   # 90° rotation
print("R @ R.T =")
print(R @ R.T)                             # should be identity
print("det(R) =", np.linalg.det(R))        # should be 1
 
theta = np.pi / 4                          # 45°
R45 = np.array([[np.cos(theta), -np.sin(theta)],
                [np.sin(theta),  np.cos(theta)]])
v = np.array([1.0, 0.0])
print("R45 v =", R45 @ v)                  # should have length 1
print("length =", np.linalg.norm(R45 @ v)) # exactly 1 — rotations preserve length

💡 Hint · Verify by computing R @ R.T — you should get the identity to machine precision. Orthogonal matrices satisfy Q · Qᵀ = I. This is the defining property. Now try it with S (the scale matrix) — S @ S.T is NOT identity; S is not orthogonal.

(d) Production reality · 15 min

War story Google · original PageRank paper· 1998~150M web pages at launch

🔥 What broke

Larry Page and Sergey Brin needed to rank the whole web. Their idea: a page's importance is the sum of the importances of pages linking to it, weighted by how many links each source has. This is a system of equations: importance vector r = M · r, where M is a huge stochastic matrix of link probabilities.

The solution is the dominant eigenvector of M (eigenvalue 1). For 150M pages, computing that directly is impossible.

🧯 The fix

Use the power iteration method: start with an arbitrary vector, repeatedly multiply by M, and it converges to the dominant eigenvector. Each multiplication is a sparse matrix-vector product (most pages don't link to most other pages), so it's tractable. This is how Google's early ranking worked, and how every modern graph-embedding technique still works.

🎓 Lesson to steal

The dominant eigenvector of a suitable matrix IS the answer to many ranking / centrality / stability questions. Power iteration lets you compute it without ever storing an n × n matrix.

Post-mortem

War story Common failure · using inverse instead of solvesilent numerical disaster

🔥 What broke

A researcher writes x = np.linalg.inv(A) @ b to solve Ax = b. On tiny well-conditioned problems, it works. On a 5000 × 5000 problem with realistic data, they get gibberish answers or NaN. Blames the algorithm.

🧯 The fix

Never form the explicit inverse. Use np.linalg.solve(A, b) or (for repeated solves) lu, piv = scipy.linalg.lu_factor(A); scipy.linalg.lu_solve((lu, piv), b). The inverse amplifies floating-point noise; solve doesn't. For overdetermined systems, use np.linalg.lstsq. For symmetric positive definite (covariance matrices), use np.linalg.cholesky.

🎓 Lesson to steal

Numerical linear algebra has strict rules. Explicit inverses are almost always wrong — every good textbook says so, and every ML paper still occasionally violates it.

War story Common failure · matrix multiplication in the wrong ordersilent, hard-to-spot

🔥 What broke

A grad student rewrites a PyTorch loss function. They swap W @ x for x @ W. Shapes still match. Training runs, loss goes down more slowly than before, model accuracy drops 5 points. Nobody notices until much later.

🧯 The fix

Always annotate expected shapes as comments (or use einops.rearrange). Add shape assertions at every transformation. If you rely on broadcasting rules, write them out: y = W @ x # W: (d_out, d_in), x: (batch, d_in) → transpose? — use W.T.

🎓 Lesson to steal

Matrix multiplication isn't commutative. Shape errors that DO throw are gifts; shape errors that don't throw are bugs waiting to burn a paper deadline.

Where this shows up in the rest of the plan

Matrices are the atom of ML infrastructure

S019/S020 · Calculus & gradient descent

Gradients are vectors; Jacobians and Hessians are matrices; chain rule is matrix multiplication.

S078 · PCA & dimensionality reduction

SVD (matrix factorisation) is the workhorse — reduce a 10 000-D user vector to 50-D without losing much.

S085 · Convolutional nets

Convolutions unroll into large sparse matrix multiplies on hardware.

S095 · Word embeddings

Embedding matrix E: E[word_id] is the row for that word. All-of-vocab similarity = E · Eᵀ.

S115 · Transformers

Attention: softmax(Q · Kᵀ / sqrt(d)) · V. Three matrix multiplies per attention head.

S125 · Distributed training

Model parallelism = splitting big matrices across GPUs; all-reduce is a giant matrix-sum.

(e) Recall + stretch · 10 min

Recall — click each to reveal · click to reveal

★ = stretch question

Explain-out-loud test

What does a matrix do to a vector, in plain English?
Why is A @ B not the same as B @ A?
What is an eigenvector, and why does that word appear in every PCA paper?

What comes next

📐MathNext →

S019 · Calculus I — Derivatives & Chain Rule

Hub: The 6-Month Learning Plan

Part of a 130-session evergreen learning series. Session structure: intuition → visual → hands-on → production war stories → recall. Duration: 90 minutes.