Appendix A: Linear Algebra Meets Programming Languages

May 2026 (10512 Words, 59 Minutes)

1. Introduction

Lecture 5 ended with a working picture of the transformer: a stack of matrix multiplications wrapped around a softmax, applied many times, with positional information added at the start. The architecture rested on one operation — matrix multiplication — repeated billions of times per forward pass.

That lecture answered what the transformer computes. This appendix answers a different question. How does a matrix multiplication actually run? When you write A @ B in Python and the answer comes back in five milliseconds, what just happened? When the same A @ B written as a triple loop in pure Python takes twenty seconds for the same numbers, what’s the difference? Why does NumPy exist? Why does Eigen exist? Why does anyone build TPUs?

The honest answer threads through four levels of abstraction. Math notation writes \(C = A B\) and stops. Flat memory turns the rectangle into a one-dimensional buffer and gives us the index formula data[i * cols + j]. Cache-aware loops rearrange the same triple loop so memory loads come in the right order — buying 5–10× for free. Parallel hardware — SIMD lanes, GPU warps, TPU systolic arrays — does many of those loops at once and turns the engineering arc into the substrate that makes deep learning at all feasible.

That is the four-level story you’ll be able to tell after this appendix. None of it is mandatory for the rest of the course. None of it appears on any exam. But every single time you see model.forward(x), somebody did this work, and the gap between knowing roughly what they did and being lost is the difference between treating deep learning as a black box and treating it as a pile of well-understood engineering decisions.

This appendix has a hands-on companion. The advanced exercise Build Your Own Matrix Class (D5) walks you through building a Matrix class in Python and C++, implementing three matmul variants, and measuring the gap to NumPy and Eigen. The lecture is the conceptual map; the lab is the coordinate system. They are independent — you can read this without doing the lab, and you can do the lab without reading this — but they are designed to fit together.

We will make exactly one assumption from Lecture 5: matrix multiplication is the central operation of deep learning. Everything else we build from the ground.

2. From Math to Memory

A matrix in a textbook is a rectangle of numbers. A matrix in a computer is a contiguous block of memory plus two integers telling you how to interpret it. The first decision in any numerical library is how to put the rectangle into the block.

2.1 The shape contract

A Matrix has a shape \((\text{rows}, \text{cols})\) and a flat data buffer of length \(\text{rows} \times \text{cols}\). The element at logical position \((i, j)\) lives at flat position

That formula is the heart of every numerical library you will ever use. NumPy uses it. C arrays use it. BLAS uses it. The C++ skeleton in the lab uses it. Once you internalise it, indexing a matrix becomes mechanical — and once you can index a matrix, the rest of the appendix is about making the indexing happen fast.

2.2 The picture

A \(2 \times 3\) matrix lives in memory as six contiguous cells:

This layout is called row-major. We lay row 0 down first, then row 1, then row 2. Walking forward through data walks across a row, hits the right edge, and continues at the start of the next row. The forward formula reads naturally: i * cols skips past the full rows we’ve already seen, and j walks across the current row.

Worked examples on the diagram (\(\text{cols} = 3\)):

The reverse direction asks the inverse question: given a flat index \(k\), which row and column does it sit at? We need this whenever we walk through data linearly — printing the buffer, debugging, transposing in one pass — and want to know in human terms where we are.

The intuition comes out of one observation about row-major storage. Every cols cells you walk through data, you cross a row boundary. So for any flat index \(k\):

• The number of full rows already passed tells us the row index \(i\). One full row holds cols cells, so this count is \(k\) divided by cols, dropping the remainder.
• The cells past the start of the current row tell us the column index \(j\). That is the remainder.

Integer division and modulo, exactly:

In Python this is k // cols and k % cols. Both are O(1) integer operations on the CPU; you can do millions per second without thinking about cost.

Reading the formulas the same way we read the forward one — left to right, with the meaning of each piece spelled out:

i = ⌊ k / cols ⌋          ← which row are we on?
        │
        └── how many *full* rows we have already passed.
            One row is `cols` cells wide, so dividing
            k by cols (and dropping the remainder) tells
            us "how many complete rows lie before us in
            the buffer".

j = k mod cols            ← how far across that row are we?
        │
        └── what is *left over* after the rows we counted.
            The modulo is exactly the part of k that
            wasn't absorbed by a full row — the offset
            into the row we currently sit on.

Worked examples on the same \(2 \times 3\) matrix (cols = 3):

Watch where the row index changes. For \(k = 0, 1, 2\) the quotient is \(0\) — we have not filled a row yet. At \(k = 3\) the quotient flips to \(1\): we have just crossed the first row boundary. For \(k = 3, 4, 5\) the quotient stays at \(1\), and the remainder cycles \(0, 1, 2\) again — we are now walking across the second row.

That cycling is the whole story. The column index \(j = k \bmod \text{cols}\) is periodic with period cols — it counts \(0, 1, \ldots, \text{cols}-1, 0, 1, \ldots\) for as long as data has cells. The row index \(i = \lfloor k / \text{cols} \rfloor\) is the count of completed periods. Together they decompose the linear walk through data into a clean sequence of full rows.

2.2.1 Forward and reverse are inverses

The two formulas are not just a matched pair — they are the two halves of a bijection between the rectangular view and the flat view. Compose them and you land where you started:

That equality is just the division algorithm from elementary number theory: every non-negative integer \(k\) splits uniquely into a quotient \(q\) and a remainder \(r\) with \(0 \le r < \text{cols}\), and \(k = q \cdot \text{cols} + r\). The forward formula reassembles the two pieces \((q, r)\) into \(k\); the reverse formula splits \(k\) into the two pieces. Same arithmetic, applied in either direction.

The bijection is not a curiosity. It is what lets us hand the flat buffer to BLAS, OpenMP, a SIMD intrinsic, or a GPU kernel and have all of them agree on what they are operating on. None of those see rows and columns — they see a double* of length rows * cols. The forward and reverse formulas are the only things that tie the rectangular world we think in to the flat world the hardware reads. Nothing else in the buffer remembers the shape; the shape lives entirely in rows, cols, and these two formulas.

You will use the forward direction far more often than the reverse — most numerical code goes \((i, j) \to k\). The reverse shows up in three concrete places in the lab: in __repr__ (walk the buffer linearly, insert a row break exactly when \(j\) wraps back to \(0\)); in a transpose written as a one-pass index remap rather than a double loop; and in any debug tool that hands you a flat index and asks what cell that is.

2.3 Why row-major (and not column-major)

The same six cells could equally have been stored as [a, d, b, e, c, f] — column 0 first, then column 1, then column 2. That layout is called column-major and the formula becomes \(\text{flat_index}(i, j) = j \cdot \text{rows} + i\). Eigen and Fortran default to column-major. BLAS accepts either via a flag. Same data, different traversal.

We pick row-major in this appendix and in the lab for three reasons:

1. It’s NumPy’s default and C’s natural array layout. When you read or write a NumPy array without specifying order='F', you get row-major. The same default holds in C: double M[2][3] is laid out row by row.
2. Walking forward through data walks across a row. It matches the way we naturally write nested loops — for i in range(rows): for j in range(cols): — and the same pattern dominates code that consumes matrix data.
3. It sets up the cache-locality story in chapter 7. For fixed \(i\) and increasing \(j\), the addresses data[i * cols + j] are consecutive in memory, so the CPU’s prefetcher streams them one cache line at a time. Strided access — walking down a column — costs more. We will measure this.

The choice between row-major and column-major is conventional, not deep. But once you’ve picked, you have to be consistent everywhere. A bug where one function thinks the data is row-major and another thinks it’s column-major is silently wrong: the indexing produces some value, just the wrong one. The number doesn’t crash; it lies.

2.4 Why this is the most important design decision

Before we talk about algorithms, before we talk about cache, before we talk about hardware — we make this decision. The shape contract is the boundary between the mathematical object (a rectangle of numbers) and the computational object (a buffer plus two integers). Every line of code in the lab — and every line of code in NumPy and Eigen — sits on top of this contract. If the contract is murky, nothing above it is robust.

That is why the lab spends most of its first chapter on the contract before letting students write a single matrix multiplication. The lab is about building; this lecture is about understanding what you built.

3. Vectors as Matrices

Most libraries — NumPy included — treat vectors as a separate type from matrices. A 1-D NumPy array np.array([1, 2, 3]) has shape (3,), not (3, 1) or (1, 3). The library handles the dimensional ambiguity for you, mostly correctly.

Pedagogically and conceptually, this special-casing is unhelpful. The lab and this lecture take the opposite path: vectors are matrices, and there is no separate vector class.

A column vector is a matrix of shape \((n, 1)\). A row vector is a matrix of shape \((1, n)\). The dot product of two vectors \(\mathbf{a}\) and \(\mathbf{b}\) of length \(n\) is the matrix product \(\mathbf{a}^\top \mathbf{b}\) of shape \((1, n) \times (n, 1) = (1, 1)\) — a scalar dressed up as a one-by-one matrix. The outer product is \(\mathbf{a} \mathbf{b}^\top\) of shape \((n, 1) \times (1, n) = (n, n)\) — a rank-1 matrix.

Once Matrix @ Matrix works correctly, every vector operation falls out for free:

There is no special code for any of these. The same matmul you write for the general \((m, n) \times (n, p)\) case handles all four — the shapes happen to collapse to the trivial values, and the result is whatever the general formula produces.

3.1 Connection to Lecture 5

This is not just a pedagogical convenience — it is the same trick the transformer is using. In Lecture 5 §5.2 — Stacking vectors and right-multiplication, every token’s query is a row of a query matrix \(Q \in \mathbb{R}^{T \times d_k}\). The query for one token is a row-vector \(\mathbf{q}_i \in \mathbb{R}^{1 \times d_k}\). The score against one key is a dot product \(\mathbf{q}_i \mathbf{k}_j^\top\). The score against all keys at once is the matrix product \(Q K^\top \in \mathbb{R}^{T \times T}\).

The library does not write a separate “compute scores for one query” function. It writes one matmul. The single-query case is just \(T = 1\). That is the same elegance — same code path for vectors and matrices — that makes attention concise to express and efficient to compute.

4. The Three Views of Matrix Multiplication

Take \(A \in \mathbb{R}^{m \times k}\) and \(B \in \mathbb{R}^{k \times n}\). The product \(C = A B \in \mathbb{R}^{m \times n}\) is one mathematical object. It can be computed three different ways. The three are equivalent — same answer, every time — but they give very different stories about the structure of the operation. The story matters because each view leads to a different implementation, with different memory and cache behaviour.

4.1 View 1 — Entry-wise (the textbook formula)

The form you saw in linear algebra class:

Each output entry is the dot product of one row of \(A\) with one column of \(B\). The implementation is a triple loop:

for i in range(m):
    for j in range(n):
        s = 0.0
        for p in range(k):
            s += A[i, p] * B[p, j]
        C[i, j] = s

This is correct. It is also the slowest of the three when written in pure Python. We will return to why in chapter 7; for now, accept it as the baseline.

4.2 View 2 — Column-wise (the linear-map view)

Each column of \(C\) is the matrix \(A\) applied to one column of \(B\):

The view: \(A\) is a linear map; \(B\) is a stack of \(n\) input vectors; \(C\) is a stack of \(n\) output vectors. To compute \(C\), apply \(A\) to each input independently. The implementation is a double loop with an inner matrix–vector product:

for j in range(n):                  # for each column of B
    col = B[:, j]                   # one input vector
    out = matvec(A, col)            # one output vector
    C[:, j] = out

The flop count is identical to view 1 — \(m k n\) multiplications and adds. The arithmetic is the same; only the story is different. This is the right view when you want to think of the operation as “apply this linear transformation to many vectors”. Multi-head attention — applying the same head to every token — is exactly this picture.

4.3 View 3 — Outer-product (the rank decomposition)

The product \(C\) is the sum of \(k\) rank-1 matrices, one per “slot” of the inner dimension:

Each term \(A_{:,p} \cdot B_{p,:}^\top\) is an \(m \times n\) outer product — a matrix of rank at most one. The full \(C\) accumulates \(k\) of these. The implementation:

C = zeros(m, n)
for p in range(k):
    col = A[:, p]                   # m-vector
    row = B[p, :]                   # n-vector
    for i in range(m):
        for j in range(n):
            C[i, j] += col[i] * row[j]

This view is less common in introductory courses but it is closer to how some BLAS routines actually work internally — the inner loop accumulates into \(C\) using one column of \(A\) and one row of \(B\) at a time, both of which are contiguous in memory under the row-major contract for the right transposes.

4.4 Equivalence is theorem; equivalence in code is proof

Linear algebra textbooks prove the three views are equivalent by index manipulation: rearrange the sums and the same scalar comes out. The lab restates that proof experimentally — students implement all three, generate random \(A\) and \(B\), and check that \(C_1 = C_2 = C_3\) within floating-point tolerance. That equality is not a sanity check on the maths; it is a sanity check on the implementation. The maths is settled. The code might have bugs.

The reason we care about three views is that they have different memory access patterns. View 1 strides through \(B\) by column, which fights the row-major layout. View 3 walks contiguously through both \(A\) and \(B\) on every step. View 2 is in between. We will see exactly what this costs in chapter 7.

5. Building Matmul Bottom-Up

The textbook triple loop is correct, but it is a pedagogical dead end. If your first encounter with matrix multiplication in code is the three nested for loops, you have learned an algorithm but you have not learned a structure. The lab takes a different path: it builds matrix multiplication from primitives, one logical step at a time, and the three views fall out at the top.

5.1 The build, in dependency order

1. Scalar operations. Matrix.__add__(self, other) adds two matrices entry by entry. Matrix.__mul__(self, scalar) scales by a number. Two double loops; no surprises. These are the linear-algebra primitives — closure of \(\mathbb{R}^{m \times n}\) under the vector-space operations.

2. Matrix–vector product. matvec(A, v) returns \(A \mathbf{v}\). Crucially, in the lab implementation, matvec is built from the primitives above — it is a sum of scaled columns of \(A\):

   \[A \mathbf{v} = \sum_{p=1}^{k} v_p \cdot A_{:,p}.\]

   No nested loops. The body is out = sum(v[p] * A[:, p] for p in range(k)) — a linear combination of columns. The implementation looks like the definition of matrix–vector multiplication as a linear map.

3. Column-wise matmul. matmul_columnwise(A, B) computes \(C\) one column at a time, using matvec from step 2 for each. Loop over columns of \(B\); call matvec; assign to the corresponding column of \(C\). Again no triple-nested loops; the inner work delegates.

4. Entry-wise matmul. matmul_entrywise(A, B) is the textbook triple loop. We write it second, not first, because by now we know what it computes. It serves as a baseline — the slowest pure-Python variant, the one whose timing we will compare to the others.

5. Outer-product matmul. matmul_outerproduct(A, B) accumulates rank-1 outer products. Loop over the inner dimension; pull a column of \(A\) and a row of \(B\); add their outer product into \(C\) in place.

5.2 Why this order

The build is mathematically motivated. Each step is one logical idea: scaled add, linear combination, “apply the map per column”, “explicit triple sum”, “rank decomposition”. The Python (and C++) code mirrors the maths line for line. There are no nested for-loops in steps 2 and 3 because the linear-algebra definition does not have nested for-loops — it has a sum of scaled vectors. The implementation in the lab uses the operator overloads (__add__, __mul__) so the code of matvec reads literally as sum(v[p] * A.get_column(p) for p in range(k)).

This is also the order in which someone implementing a numerical library would think about the operation. You don’t reach for the triple loop and call it a day; you first understand the operation as a linear map (one matrix–vector product), then as a stack of linear maps (column-wise), then as a rank decomposition (outer-product). The triple loop is a flat unrolling of the column-wise view — fine to implement once you have the structure, dangerous to start from.

5.3 A note on naming

The lab is careful about naming. Class attributes — long-lived, public — get descriptive names: self.rows, self.cols, self.data. Local variables in numerical methods — short-lived, scoped to one function — match the maths: nrows, ncols, flat, col_p, row_p. Loop indices over mathematical objects are single letters: i, j, k, p. Loop indices over collections of “things” stay descriptive: for row in rows:.

This is not “lazy variable naming” — it is a style choice that reflects the audience of each identifier. When the code implements the formula \(C_{ij} = \sum_p A_{ip} B_{pj}\), the Python translation is more readable, not less, when it reuses the letters \(i, j, p\). Replacing them with row_index, column_index, inner_index breaks the visual mapping to the maths and makes the line harder to parse at a glance. Fortran (1957), MATLAB, R, NumPy, SciPy, JAX, PyTorch, BLAS, LAPACK, Eigen — every numerical library you will ever use — converges on this convention because it is the right one for code that implements maths.

6. Why NumPy Is So Fast

You have just spent four chapters building matrix multiplication from primitives. You can run the lab, compute a 256×256 matmul with your own three-views Python implementation, and watch it take about two seconds. Then you can run the same multiplication via numpy.matmul and watch it return in about two milliseconds. Same big-O. Same flop count. Same numbers in. Same numbers out. NumPy is a thousand times faster.

That gap is the core engineering puzzle of the appendix. It is not one trick; it is the cumulative effect of four engineering decisions, each of which costs Python.

6.1 Reason 1 — The interpreter is in the loop

Every data[i * cols + j] in your Python code goes through the CPython interpreter: bytecode dispatch, attribute lookup, type check, allocation of an intermediate float object. For a 256 × 256 matmul with the textbook triple loop, that is about \(256^3 \approx 16{,}777{,}216\) trips through the interpreter. Each trip costs hundreds of nanoseconds even when nothing real is happening.

NumPy’s inner loop is a tight C for loop. The interpreter is called once per matmul to enter numpy.matmul, and from there the work happens in compiled C — no bytecode, no attribute lookup, no boxing. The interpreter overhead does not scale with the matrix size; for a 256 × 256 matmul it amounts to roughly nothing. That alone is most of the gap.

6.2 Reason 2 — Packed doubles, not boxed Python floats

A Python float is a heap-allocated object. Its layout includes a refcount, a type pointer, the IEEE-754 double, and some padding — about 28 bytes per number on a typical 64-bit build. A NumPy float64 array is a contiguous block of doubles, eight bytes each. Same numerical precision, 3.5× less memory.

Memory is not just about footprint. Modern CPUs run far faster than they read from main memory; they hide the latency by streaming contiguous data through the L1 / L2 / L3 cache hierarchy. The cache works in 64-byte cache lines. With NumPy’s packed doubles, one cache line is exactly 8 numbers — one fetch buys you 8 doubles’ worth of work. With Python floats, the floats themselves live somewhere on the heap, indirected through pointers, and the contiguous “array” is really an array of pointers to heap objects. The CPU prefetcher cannot see the floats; it sees the pointer array. Cache locality is destroyed before the loop starts.

6.3 Reason 3 — BLAS dispatch

For matrices large enough to make the per-call overhead worth amortising, NumPy hands the multiplication off to a tuned BLAS (Basic Linear Algebra Subprograms) implementation — typically OpenBLAS, MKL, or, on Apple Silicon, the Accelerate framework. BLAS is a 50-year-old library specification; modern implementations are decades-tuned, hand-written assembly in their inner kernels, with cache blocking, register allocation, and architecture-aware code paths.

You will not out-write BLAS in an afternoon. That is fine. The point of BLAS is that you don’t have to.

6.4 Reason 4 — SIMD

A modern CPU has vector instructions — single instructions that operate on multiple doubles at once. AVX2, on most laptops since 2013, multiplies 4 doubles per cycle. AVX-512, on server-class chips, multiplies 8. ARM’s NEON does 2; Apple’s M-series processors do 2 with NEON and many more with the Accelerate framework’s matmul kernels. BLAS uses these. Your Python loop does not.

The factor here is bounded — you cannot beat the SIMD width — but a 4× to 8× factor is real, and it stacks on top of the other reasons. Combined with reasons 1, 2, and 3, you get the 100× to 1000× gap.

6.5 How much each reason matters

The factors are not equal in size and not independent of matrix size. For tiny matrices (8 × 8) the per-call overhead of NumPy itself dominates and the gap shrinks. For 256 × 256, the interpreter overhead and packed-doubles effects are the largest contributors. For 1024 × 1024 and up, BLAS dispatch + cache blocking + SIMD dominate, and the gap can exceed 1000×.

The lab measures all of this. The deliverable is a single log–log plot showing pure-Python triple-loop time versus NumPy time, across sizes. The shape of the plot tells the story: parallel lines, with NumPy two to three decades below. Same growth rate (cubic), vastly different constants.

7. Loop Order and Cache Locality

The four reasons in the previous chapter are why NumPy is fast. But you can move a smaller piece of the gap without leaving Python’s level of abstraction — let alone Python itself. The simplest such move is changing the order of the three nested loops in the textbook formula. We are about to get a 5–10× speedup for free, just by writing the loops in the order memory wants.

7.1 The setup: ijk and ikj

The textbook formula is

and the natural translation is the ijk loop:

for (i = 0; i < M; ++i) {
  for (j = 0; j < N; ++j) {
    for (p = 0; p < K; ++p) {
      C[i][j] += A[i][p] * B[p][j];
    }
  }
}

Equivalently, we can swap the inner two loops to ikj:

for (i = 0; i < M; ++i) {
  for (p = 0; p < K; ++p) {
    double a_ip = A[i][p];
    for (j = 0; j < N; ++j) {
      C[i][j] += a_ip * B[p][j];
    }
  }
}

These compute the same numbers. Same flop count. Same algorithmic complexity. Compile both at -O2 and you will see ikj run 5 to 10 times faster on most modern machines. There is no cleverness; the only difference is the order in which the multiplication consumes memory.

7.2 Why ikj wins: cache lines

Recall that B[p][j] lives in memory at flat index p * N + j. Walk through the inner loops of each version and watch which j’s and p’s vary together:

• ijk inner loop: p increments while i and j are fixed. The address B[p][j] jumps forward by N (one full row) every iteration. For N = 1024 doubles, that is 8 KB of stride per step. Each step crosses many cache lines; almost every access misses the L1 cache.
• ikj inner loop: j increments while i and p are fixed. The address B[p][j] increments by exactly 1 double each iteration — contiguous in memory. The CPU’s prefetcher sees the linear pattern, prefetches the next cache line before you need it, and the inner loop runs at memory bandwidth.

A modern x86 cache line is 64 bytes — 8 doubles. The ikj inner loop therefore amortises one cache miss across 8 multiplications. The ijk inner loop pays one cache miss per multiplication.

7.3 A tiny diagram

That is the whole story: row-major data wants to be walked along rows, not down columns. ikj walks rows of B in its inner loop; ijk walks down a column.

7.4 Free performance

Notice what we did not do. We did not change the algorithm. We did not parallelise. We did not use SIMD. We did not even change the language. We rewrote three nested for loops in a different order, and the program ran 5–10× faster.

This is the most memorable lesson in the lab. The compiler cannot fix it for you (-O3 will sometimes reorder loops, but only when it can prove the reorder is safe; for naive triple loops with arrays of pointers it usually cannot). The hardware cannot fix it for you. The libraries do it for you when you call them — but if you write the loops yourself, you have to write them in the right order.

The trade-off: ikj requires you to allocate C zeroed and accumulate into it (each output entry is touched many times), where ijk computes each output entry to completion before moving on. For numerical stability the difference is invisible. For performance the difference is everything.

8. Compiler Leverage: -O0, -O2, -O3

The previous chapter was about a change you made: rearranging loops. This chapter is about a change the compiler makes when you flip a flag. Both buy you orders of magnitude. Both are essentially free if you understand them.

C++ compilers ship a sequence of optimisation levels, named by their flag: -O0, -O1, -O2, -O3. Lower numbers mean less compiler effort and slower output; higher numbers mean more compiler effort and faster output. Same source code; the compiled binaries can differ by an order of magnitude.

8.1 What each level does (roughly)

• -O0 is “compile literally”. Every variable lives on the stack. Every load is from memory. No loop unrolling. No vectorisation. No inlining. The compiled binary is a faithful translation of your source code, line by line. This is what you want when stepping through code in a debugger; this is not what you want for performance.

• -O2 is “compile carefully”. The compiler unrolls small loops, allocates frequently-used variables to registers, inlines small functions across translation units, eliminates dead code, propagates constants, and starts to vectorise simple loops with SSE / AVX. This is the standard “release” build of most C++ projects.

• -O3 is “compile aggressively”. On top of -O2, the compiler does loop interchange (sometimes), more aggressive auto-vectorisation, function specialisation, and occasionally recognises dense matrix multiplication patterns and emits much tighter code. It also enables some optimisations that change floating-point semantics in tiny ways (-ffast-math is a separate flag, but -O3 is the level at which serious vectorisation kicks in for many compilers).

There is also -Os (“optimise for size”) and -Ofast (more aggressive than -O3, also enables -ffast-math). For this appendix -O0, -O2, -O3 are the three you care about.

8.2 What you measure

In the lab, you build the C++ matmul at three optimisation levels and benchmark a 1024 × 1024 matmul at each. Typical numbers (on a 2024-era laptop):

Numbers will differ on your hardware, but the shape of the table is robust. -O0 is the slow path; -O2 is the standard path; -O3 is the fastest path the compiler can give you for this code.

8.3 Why this is “free”

The same source code compiled at -O0 and -O3 produces different binaries. You did not change a line. You did not call a different library. You changed one flag in CMake (-DCMAKE_BUILD_TYPE=Release), and the compiler did the rest.

The catch — and there is one — is that the compiler can only optimise what it can prove safe. A naive triple loop over double arrays is straightforward to vectorise; a triple loop with pointer chasing through a std::vector<std::vector<double>> is much harder, because the compiler cannot prove the inner pointer is contiguous. The lab uses a flat std::vector<double> for exactly this reason — it gives the compiler the simplest possible memory model to reason about.

This is also why the gap between Python and naive C++ at -O3 is roughly 50× while the gap between Python and NumPy is 100–1000×: a lot of NumPy’s advantage is the C++ compiler’s advantage, applied to NumPy’s own inner kernels. NumPy is not magic; NumPy is a well-engineered C library that benefits from forty years of C compiler optimisation.

9. Eigen and BLAS

-O3 C++ takes you most of the way. It does not take you all the way. The remaining gap, between hand-written -O3 ikj and a serious linear-algebra library, is a small constant factor — usually 2–4× for medium sizes, sometimes more for large sizes when multi-threading kicks in. That gap is not free. It is the result of techniques the compiler does not know how to apply automatically. Eigen (the C++ library used in the lab) and BLAS (the C / Fortran library NumPy and Eigen both call for big enough operands) close it with four engineering tricks.

9.1 Trick 1 — Expression templates

When you write A * B in Eigen, the library does not immediately compute the product. It builds a small abstract syntax tree node — “I have an operation: A times B. The result type is a matrix of shape (rows of A) × (cols of B), but I have not allocated it yet.” When you then write C = A * B, Eigen recognises the assignment and generates code that fuses the multiplication and the assignment — writing directly into C’s buffer, with no intermediate temporary matrix.

This sounds like a small thing. It is not. A naive C++ Matrix operator*(const Matrix&, const Matrix&) returns a fresh matrix by value, which means an allocation, the multiplication into that fresh buffer, and (if the calling code has not been written carefully) a copy out. Eigen’s expression templates skip all of that. For chained operations — D = A * B + C — the savings compound.

9.2 Trick 2 — Cache blocking

The -O3 ikj loop streams along memory contiguously. That is good. But for matrices large enough that the whole matrix does not fit in L1 cache (32 KB on most CPUs, room for a ~64 × 64 double matrix), the inner loop still issues loads that miss the cache.

Cache blocking partitions the matrix into tiles small enough to live entirely in L1, multiplies tile-by-tile, and accumulates. The textbook three nested loops become six nested loops — three outer ones tiling the matrix into blocks, three inner ones doing the regular triple loop within each block. The inner triple loop now operates on data that fits in L1, so the reuse is high and the misses are rare.

Eigen does this. BLAS does this with multiple levels of blocking — L1, L2, sometimes L3 — and the inner kernels are written in assembly. This is the technique that gives big-matrix performance its 5–10× advantage over naive -O3 code.

9.3 Trick 3 — Hand-tuned SIMD kernels

The -O3 compiler can auto-vectorise, but its vectoriser is conservative — it only transforms loops it can prove safe. The innermost kernel of a serious BLAS implementation (the part that handles a 4 × 4 or 8 × 8 block at a time) is not written in C; it is written in intrinsics — function-shaped wrappers around the CPU’s vector instructions — or in raw assembly. The kernel author chose specific SIMD lanes, specific load patterns, specific accumulation registers, and laid out the inner loop to keep the CPU’s pipeline full.

Eigen ships such kernels for x86 (SSE2, AVX, AVX2, AVX-512) and ARM (NEON). NumPy gets the same kernels for free by linking against OpenBLAS or MKL. The C++ compiler at -O3 might auto-vectorise your hand-written ikj loop — but not as well, and not as reliably across sizes.

9.4 Trick 4 — BLAS dispatch

For matrices below a size threshold, Eigen does the work itself with its own kernels. Above the threshold (typically a few hundred rows), Eigen forwards the call to BLAS. BLAS is the seriously-tuned library; calling it is what NumPy effectively does for any matmul where the per-call overhead is amortised by the work.

The threshold exists because BLAS calls have non-trivial overhead — argument marshalling, function dispatch, sometimes thread pool startup. For a 16 × 16 matmul, going through BLAS is slower than just doing the work directly. Eigen knows this and routes accordingly.

9.5 How much each trick matters

For a 1024 × 1024 matmul, the breakdown is roughly:

• Hand-written -O3 ikj: 800 ms
• Add cache blocking: 300 ms
• Add hand-tuned SIMD kernels: 100 ms
• Add multi-threading + BLAS dispatch: 25 ms

The factors stack. The lab measures the endpoint — your -O3 ikj versus Eigen — and you will see the gap. Reproducing the cause of the gap (writing your own cache-blocked SIMD kernels) is a multi-week project. Understanding the cause is the work of this chapter.

10. Parallel Hardware: SIMD, GPUs, TPUs

We have stayed on a single CPU core for nine chapters. That is an enormous amount of engineering — packed memory layout, cache-aware loop order, compiler optimisation, expression templates, hand-tuned kernels — and it gets us within a factor of a few from theoretical CPU peak. To go further, we leave the CPU. Modern deep learning runs on parallel hardware: SIMD lanes inside the CPU, GPUs alongside it, and TPUs replacing it. This chapter is the bridge from “fast on one core” to “fast at scale”.

10.1 SIMD lanes — the inside view

Every modern CPU core has vector execution units alongside the scalar ones. AVX2 (Intel / AMD, 2013+) processes 4 doubles per cycle in a single instruction. AVX-512 (Intel server chips, some 2019+ desktop chips) processes 8 doubles. ARM’s NEON (Apple Silicon, modern phones) processes 2 doubles; SVE (newer ARM server chips) is variable-width and goes up to 32. Each vector instruction is one CPU-level operation; the cost is roughly the same as the scalar equivalent.

This is the smallest unit of parallelism. The compiler tries to use it automatically; serious libraries use it deliberately. The width matters: AVX-512 is roughly 2× faster than AVX2 for SIMD-friendly code. That factor explains why the same numerical library can be twice as fast on a server CPU as on a laptop, even at the same clock speed.

For matrix multiplication, SIMD is a near-perfect fit. The inner-most accumulation c += a * b is a fused multiply-add (FMA) instruction; the vector form does 4 or 8 of those per cycle. A serious BLAS kernel issues two FMA instructions per cycle (using two execution units), reaching 8–16 doubles’ worth of work per cycle on a single core. Multiply that by clock speed and you have CPU peak.

10.2 GPUs — many cores, lock-step execution

A GPU has thousands of cores. Not “thousands of CPU cores” — they are simpler than that. A GPU core is more like a single SIMD lane: it executes one operation per cycle. What makes a GPU enormous is how many of them there are (5,000+ on a high-end card) and how they are organised.

The organisational unit is the warp (NVIDIA) or wavefront (AMD) — 32 (or 64) threads that execute the same instruction in lockstep, on different data. This is sometimes called SIMT (single instruction, multiple threads) — it is morally SIMD with a thousand lanes. GPUs are extraordinarily good at workloads that look the same across many data points: matrix multiplication, convolution, element-wise operations. They are bad at workloads with branchy control flow, because divergence between threads in a warp serialises execution.

The other thing GPUs have is memory hierarchy designed for streaming. Global GPU memory (HBM, currently up to ~3 TB/s per accelerator) is much faster than a CPU’s main memory but still slow compared to compute. Each warp has access to a small, very fast on-chip shared memory (similar in role to a CPU’s L1 cache). Cache blocking on GPUs is even more important than on CPUs, because the bandwidth-to-compute ratio is more skewed.

For matmul specifically, modern GPUs ship dedicated matrix multiply hardware: NVIDIA’s Tensor Cores (Volta, 2017+) operate on 4×4 or 8×8 blocks in a single cycle, in mixed precision. AMD’s Matrix Cores are equivalent. These are not general-purpose units; they only do matmul. They are why a current data-centre GPU has tens of TFLOPS in standard precision but hundreds of TFLOPS in matmul-shaped operations.

10.3 TPUs — only matrix multiplication, but very fast

Google’s TPU is a different bet. Where GPUs are general parallel processors that happen to have matmul accelerators, a TPU is a matmul accelerator with some scalar bits attached. The core unit is a systolic array — a 128×128 grid of multiply-add cells, each connected to its neighbours. Data flows through the array in a wave: the top row consumes the operand, multiplies, and passes the result downward; the next row picks it up, multiplies, and passes downward again. In one cycle, a 128×128 systolic array does 16,384 multiply-adds. The whole TPU chip is one (or a few) of these.

This makes TPUs extraordinarily fast for matrix multiplication and almost useless for anything else. Convolutional networks and transformers — both dominated by dense matmul — run extremely well. Branchy code, irregular access, dynamic shapes — these the TPU does badly. The constraint is the design: the systolic array’s geometry only matches matmul; everything else goes through a much smaller scalar pipeline.

The other TPU design choice is bfloat16 — a 16-bit float format that has the same exponent range as float32 but fewer mantissa bits. Mantissa precision turns out to matter less than range for neural network training, so bfloat16 gives roughly 2× the throughput at very modest cost in accuracy. Modern GPUs adopted bfloat16 too; it is the dominant training format in 2026.

10.4 Why deep learning is matmul-bound

The architectural arc — CPU → SIMD → GPU → TPU — is the same story told four times: make matrix multiplication faster. Every accelerator is, at its core, a matmul accelerator. Convolutions reduce to matmul (im2col + matmul). Attention reduces to matmul (Q K^T and weights @ V). Feed-forward reduces to matmul. The dominant FLOPs of a forward pass — the dominant energy, the dominant time, the dominant cost — are matrix multiplications.

That is what we mean when we say deep learning is matmul-bound. The architecture of every accelerator on the market is a bet that the workload will be at least 80% matmul; the other 20% goes through a smaller, more general scalar / vector pipeline, and you choose the workload so that matmul stays in the majority.

This explains why “make this transformer faster” almost always reduces to “make matrix multiplication faster”, which reduces to either “use better hardware” or “use a less expensive variant of attention that is still matmul” (linear attention, sparse attention, low-rank approximations). The space of architectures is constrained by what your accelerator does well.

11. The Full Hierarchy

We have walked the engineering arc from raw Python to a TPU. It is worth seeing it in one place. Below is a schematic table — orders of magnitude only, current as of 2026, on representative hardware. Treat it as a map, not a benchmark. Your numbers will differ by hardware, problem size, and a thousand other things, but the shape — each step roughly an order of magnitude — is robust.

The first column is FLOPS — multiplications and adds per second — for dense matmul on the workload’s home turf (large matrices for the fast hardware, small ones for the slow). The “vs. pure Python” column is the cumulative speedup over the slowest layer.

11.1 What this table says

Each layer is roughly an order of magnitude faster than the previous, and the engineering arc tells you why:

• Python → NumPy (\(10^4\)): kill interpreter overhead, pack the doubles, dispatch to BLAS, use SIMD.
• NumPy → Eigen / -O3 C++ (\(\sim 1\)): not much, in fact. The two are comparable for many sizes. NumPy is a Python wrapper around C kernels; Eigen is a C++ library with similar kernels. The differences come from threading, dispatch overhead, and which BLAS implementation links in.
• CPU library → GPU (\(10^3\)): thousands of compute lanes instead of dozens, bigger memory bandwidth, dedicated matmul hardware (Tensor Cores).
• GPU → TPU (\(\sim 10\)): more matmul-specialised silicon, less general-purpose silicon. Faster for matmul, slower for everything else.

11.2 What this table does not say

The table omits energy, cost, and ergonomics. A TPU pod costs tens of millions of dollars and is impossible to access except through Google Cloud. A laptop GPU costs hundreds. NumPy costs your time installing it. The price-performance curve and the ergonomics curve do not move in lockstep with the throughput curve, and the right tool depends on what you are actually optimising. We come back to that in the next chapter.

The other thing the table does not capture is the engineering arc itself. None of these layers existed at full power on day one. NumPy started in 2006 as a successor to Numeric (1995). Eigen 1.0 was 2008. CUDA was 2007. The first TPU was 2015. Each layer is the work of hundreds of engineer-years. The fact that you can write import numpy as np; A @ B and get \\(10^{10}\\) FLOPS is the result of two decades of compounding engineering investment in the tools that turn linear algebra into running code.

That is the engineering arc that made deep learning possible. The transformer architecture is from 2017; the GPU computing stack that made training a transformer feasible at scale is from 2007; the matrix multiplication routines those GPUs run are descendants of BLAS routines from 1979. Every stack you use today is sitting on layers of work you did not see.

12. When You’d Actually Write This Yourself

You now have a conceptual map of the engineering arc. The natural follow-up question is: should you write your own matmul?

The honest answer is almost never. For 99% of working programmers, the right call is numpy.matmul (in Python), Eigen::Matrix (in C++), or torch.matmul / jax.numpy.matmul (in any language with PyTorch / JAX bindings). These libraries have absorbed every trick in this appendix and a thousand others. Beating them is a multi-month project, not a Tuesday afternoon.

When would you write your own? A short list, in roughly increasing order of plausibility:

• Embedded targets without BLAS. Microcontrollers. Specialised DSPs. Anywhere pip install numpy is not an option. You write a small, focused matmul tuned to your hardware and constraints.
• A new hardware architecture. TPU kernels are written by Google. NVIDIA writes the cuBLAS kernels for new generations of GPU. AMD writes ROCm. Each new accelerator generation requires a fresh round of kernel development. This is a job, and the people doing it are extraordinary engineers.
• Beating the library on a specific shape. Sometimes the library is general; your problem is a known specific shape; you can beat the library by 20% by writing a kernel that exploits the shape. Modern compiler infrastructures (Halide, TVM, Triton, JAX with jax.jit) automate much of this, but the human in the loop is still the deciding factor.
• Education. The lab.

For the remaining 1% of cases, you now have a map of where to look. You know what cache blocking is, what SIMD is, what expression templates buy you. You will not be lost when you read the OpenBLAS source code; you will know what each layer is doing.

For the 99%: use the library. The library is faster than what you would write. The library is also less buggy than what you would write — these libraries have had decades of fuzz testing, edge-case handling, and numerical-stability auditing that no individual project replicates. The right answer to “should I implement my own matrix class?” in production is, with extremely few exceptions, no.

12.1 The lab as the controlled exception

The hands-on companion to this lecture is the Build Your Own Matrix Class exercise (D5). It is the controlled exception — the one place in this course where you do implement matrix multiplication from scratch, in two languages, three ways each. The point is not that you will need this skill in production. The point is that having done it once, the words in this appendix mean something specific to you. “Cache locality” is not abstract; it is the 5× speedup you measured between ijk and ikj. “BLAS dispatch” is not folklore; it is the reason your -O3 C++ at 1024 × 1024 is twice as slow as Eigen.

The lab does in code what this lecture does in prose. Working through it, even partially, makes the conceptual map concrete. Skipping it and reading only the lecture is fine — the conceptual map stands alone. But if you have the time, the lab pays for itself.

12.2 Where this appendix sits in the course

This appendix is the only chapter of the course that goes below the level of numpy.matmul. The lectures before it treat matrix multiplication as a primitive — Lecture 5’s transformer is built out of matmuls and you are not asked to implement them. The lectures after it do the same. The appendix is here as a foundation — once, in one place, the course says: this is what is happening underneath. Everything else can refer back to it without re-deriving.

The next regular lecture continues the transformer story. You can return to this appendix any time you wonder why a model trains in 30 minutes on a GPU and 30 days on a CPU, why bfloat16 is the format everyone uses, why TPUs were built. The answers are all here. The map will be useful.

That is the appendix. Linear algebra meets programming languages, programming languages meet hardware, hardware enables deep learning. Every step is engineering. None of it is magic.