One matrix multiply is doing all the work.
Late-interaction search scores a query against every token of every document. Before class 04 replaces floats with bits, this class shows the float path properly: how MaxSim collapses into a single GEMM (general matrix multiply), why GEMM is the most optimized routine in all of computing, and the measured bar it sets — the one every clever kernel afterwards has to beat. One tiny worked example — a 4-token query and three toy documents — carries through every chapter.
MaxSim: every query token picks its favorite
A query is a bag of token vectors; so is a document. To score the pair: each query token finds the doc token it matches best (a max over dot products), and the matches are summed. That's the whole rule. Score a document and watch:
Try all three docs. Doc A wins (28), then B (20), then C (16) — because more of A's tokens are good matches for some query token. That per-token max is what makes late interaction sharper than squashing each document to a single vector: one great matching token can't be averaged away by fifty irrelevant ones.
The repackaging: S = Q · Dᵀ
Chapter 01 computed a grid of dot products — every query row against every doc column. There is a name for "grid of all dot products between two sets of vectors": a matrix multiply. Hover any cell of S:
Nothing new was computed — S = Q @ D.T contains exactly the 12 dot products
chapter 01 made one at a time. But the repackaging matters enormously, because the moment
your problem is a matrix multiply, you get to hand it to the most relentlessly
optimized routine in computing: GEMM, the heart of every BLAS library
(Accelerate, OpenBLAS, MKL…). The next four chapters are about why that hand-off is worth
10–100× — and what it still leaves undone.
The interpreter tax
Same math, three ways to ask numpy for it — on the repo's real corpus (full SciFact: 5,183 docs, 1.19M tokens, dim 128, 32-token query). Every strategy computes the identical 38 million dot products. Only the number of Python→C round trips changes:
Each Python-level call costs ~1 µs of ceremony — argument parsing, type dispatch, array bookkeeping — before a single multiplication happens. The actual math for this query is 18 ms (next chapters explain how). Strategy one spends two thousand times the entire answer on ceremony alone. And the per-document strategy isn't just taxed by calls: 5,183 tiny matmuls never let BLAS stretch its legs — no packing amortization, no threads, cold caches every time.
This is the first lesson of the float path: batch everything into one call, then make that one call unbeatable. Which raises the question — what makes GEMM itself so fast?
Wall-to-wall FMA
You already know GEMM's inner loop — it's class 01, chapter 06. Every cell of S is a dot product, and a dim-128 float dot product is 32 fused-multiply-add instructions, 4 lanes each, accumulators kept vertical until one fold at the very end:
// one register tile of S: 4 rows × 2 vector columns, // accumulated across the whole shared dimension k loop k: ld1 {a0.4s}, [Q_row0, k] // reuse ↓ ld1 {b0.4s, b1.4s}, [D, k] // reuse → fmla acc00.4s, a0.4s, b0.4s fmla acc01.4s, a0.4s, b1.4s // … 6 more fmla: every loaded value // used in MULTIPLE multiply-adds
one dot product: load 2 vectors → 1 fmla
1 load per multiply-add
4×2 register tile: load 6 vectors → 8 fmla
0.75 loads per multiply-add
bigger tile → more reuse per load →
the FMA units stay fed instead of
waiting on memoryThe instruction mix of a good GEMM is startling: it is almost nothing but
fmla (or AVX's vfmadd) back to back — the loads hide behind them.
A modern core can start multiple FMA instructions per cycle; keeping that pipeline
full is the entire game. Notice what makes the tile work: every value loaded into a register
is used in several multiply-adds, not one. That word — reuse — is
the real secret, and it goes far deeper than registers.
The cache game
RAM is ~100× slower than a register. In between sit small fast caches — and a matrix doesn't fit in them. Below, the same 8×8 multiply computed two ways; each cell of matrix B shows how many times it had to be fetched from slow memory. Watch the colors pile up:
Arithmetic intensity: FLOPs per byte
Why can GEMM be tiled into speed while the humble dot product can't? Count two things for an n×n multiply: arithmetic is 2n³ FLOPs, but the data is only 3 matrices ≈ 12n² bytes. Work grows a whole power of n faster than data — drag n:
Now the contrast that explains class 01's horizontal tax chapter: a single dot product does 2n FLOPs on 8n bytes — intensity 0.25, forever, no tiling can help because no value is ever reused. A dot product is doomed to run at memory speed; a matrix multiply is the rare computation that can saturate the arithmetic units of the machine. MaxSim scoring happens to be a bag of dot products that share both operand sets — which is exactly the property that lets chapter 02's repackaging turn doomed loops into a compute-bound GEMM.
What GEMM can't do: max
GEMM gives us every dot product. But MaxSim needs a max over each document's columns — and max isn't matrix algebra. nano-plaid's answer: concatenate all docs' tokens into one wide matrix, do ONE giant GEMM, then take segmented maxes at the document boundaries. Run it on our three docs:
The scores match chapter 01 exactly — 28, 20, 16 — but notice what changed: chapter 01
scored documents one at a time; this scored the whole corpus in one GEMM +
one reduceat + one sum, three numpy calls total. This is verbatim from
nanoplaid.py:
sim = q @ corpus.T # [nq, total_tokens] — THE GEMM per_doc = np.maximum.reduceat( # segmented max at doc starts sim, doc_offsets, axis=1) # [nq, n_docs] scores = per_doc.sum(axis=0) # one score per document
row of sim: [ 8 −5 −3 │ 4 5 │ 7 0 ]
offsets: 0 3 5
└─max─┘ └max┘ └max┘
per_doc row: [ 8, 5, 7 ]The measured bar
Full SciFact on an Apple M4: 5,183 docs, 1.19M tokens, dim 128 — that's
9.7 GFLOP per query for the exhaustive GEMM. Measured, end to end, p50 per
query (from the repo's eval.py --profile):
Read the bars in order. First: brute-force GEMM does 9.7 GFLOP in 18 ms — about half a teraflop per second, sustained, from three lines of Python. That is the bar. Second: the "smart" two-stage residual index — probe centroids, gather candidates, decompress, rescore — loses to brute force by 4–6× in numpy. Clever indexing has real constant costs (gathers, decode, Python glue), and it must pay them back before its asymptotic win shows up; at 5K docs it never does. Third: the only thing that beat the GEMM here didn't out-multiply it — it changed the problem: 1-bit documents scored by a fused SIMD kernel, 5.7 ms at 1/25th the memory.
That's the strategic lesson of the whole float path: you don't beat BLAS at its own game. You beat it by needing less — fewer candidates (pruning), fewer bits (quantization), or both.
The limit: GEMM only speaks float
So why not compress documents to 1 bit for storage, and decompress back to floats at query time — keeping our beloved GEMM for the scoring? Because of what decompression does to every doc token, on every query:
For 1-bit documents there is a better move: score the bits directly, without ever decompressing — an int8 query row against packed sign bits, through one line of algebra and the dot-product instructions from class 01. Different arithmetic, same fusion instincts, three CPUs. That kernel is class 04. But it raises the question the kernel classes quietly assume is settled: which bits should we store in the first place? Sign bits are one answer among several — class 03 lays out the whole menu (scalar, binary, centroid + residual, product quantization) and measures what each choice costs in ranking quality before we spend three classes making the survivors fast.
Six questions
Instant feedback, explanations included. Miss one? The chapter number is next to it.