P14 hands-on — Hardware-aware ML, block by block

Two measured numbers predict a workload, and two modelling bugs get caught.

Source: handson/h14_hardware.py --- run it with python3 handson/h14_hardware.py
Full project spec: P14 — Hardware-Aware ML System

A roofline is two measured numbers and a claim: no kernel can run faster than its arithmetic intensity times the machine's bandwidth, or than the machine's peak compute, whichever is lower. The claim is falsifiable, which is the whole point --- a measurement that beats it is proof that the model is wrong.

This page catches two such bugs. The first version measured peak throughput in fp64 and priced an fp32 workload against it, a factor of four, and the assembly duly reported kernels running at twice the speed of light. The second is subtler and lives in the byte count: a 25 MB working set that stays in cache never pays the DRAM cost the model charges it, and the ratio wanders across 1.0 from run to run as a result.

Both were invisible in the absolute timings. Both were unmissable the moment the number was divided by a bound it could not legally cross.

Contents

• Block 1 — Measure the machine, not the spec sheet
• Block 2 — Arithmetic intensity
• Block 3 — Predict, then measure
• Block 4 — Tiling
• Block 5 — Quantisation
• Block 6 — Batching and Little's Law
• Block 7 — The decode wall
• The assembly
• The design space
• Two modelling bugs this page catches
• The hardware, level by level
• Optimisation levers, in the order they pay
• Advanced topics
• How this connects to the rest of the track
• Failure modes at scale
• Primary sources
• Running it
• What to do with this

How to read this page

Each block below is a self-contained lego piece: it builds one mechanism, proves it works on its own, and returns what the next block needs. The code is the real source, sliced out of the script. The output underneath it is the real output, captured by running that script --- not transcribed, not idealised. Where a measurement contradicted what I expected, the contradiction is in the output and the prose says so.

The assembly at the end wires every block into one working thing and measures it.

Block 1 — Measure the machine, not the spec sheet

Teaches: the roofline is per-dtype, and that bites

The problem. A roofline is two numbers and a falsifiable claim. Getting either number wrong makes the model report impossible results — which is the most useful failure a performance model can have.

@block(1, "Measure the machine, not the spec sheet", "the roofline is per-dtype, and that bites")
def b1(s, show):
    n = 1024
    peaks = {}
    for dt in (np.float64, np.float32):
        A = np.random.rand(n, n).astype(dt); B = np.random.rand(n, n).astype(dt)
        peaks[np.dtype(dt).name] = 2 * n**3 / best(lambda: A @ B, 5)
    # peak bandwidth: a streaming triad, far larger than any cache
    N = 24_000_000
    x = np.random.rand(N); y = np.random.rand(N); z = np.empty(N)
    def triad(): np.multiply(x, 2.0, out=z); np.add(z, y, out=z)
    tb = best(triad, 5)
    bw = 3 * N * 8 / tb
    if show:
        print(f"  peak bandwidth : {bw/1e9:>8.1f} GB/s   "
              f"(triad over {N*8/1e6:.0f} MB in {tb*1000:.1f} ms)")
        print(f"  {'dtype':<10}{'peak GFLOP/s':>15}{'ridge point':>15}")
        for k, v in peaks.items():
            print(f"  {k:<10}{v/1e9:>15.1f}{v/bw:>13.2f}")
        print(f"  fp32 is {peaks['float32']/peaks['float64']:.2f}x fp64 -- two SIMD lanes")
        print("  per register instead of one, and on this machine a little more than")
        print("  the theoretical 2x because the fp32 kernel also gets better cache")
        print("  reuse per byte.\n")
        print("  THERE IS NO SUCH THING AS 'the' ROOFLINE. There is one per dtype, and")
        print(f"  the ridge moves with it: {peaks['float64']/bw:.1f} FLOP/byte in fp64, "
              f"{peaks['float32']/bw:.1f} in fp32. Mixing")
        print("  them is not a rounding error -- the first version of this file")
        print("  measured peak in fp64 and predicted an fp32 workload, and the")
        print("  assembly reported kernels running at 2x the speed of light. A ratio")
        print("  below 1.0 against a roofline is never a fast kernel; it is always a")
        print("  broken model, and it is the most useful bug the roofline can produce")
        print("  because it is impossible to rationalise away. There are exactly two")
        print("  ways to get one: the wrong ceiling (this bug) or the wrong byte count")
        print("  (a working set that never left cache -- the assembly measures that")
        print("  one too).")
        print("  Both numbers are MEASURED. The vendor's peak assumes an FMA on every")
        print("  port every cycle, which no real kernel reaches; numbers.md records")
        print("  this machine's bandwidth at 57.5-99.1 GB/s depending on working set,")
        print("  and this triad sits at the streaming end of that range.")
    return {"flops": peaks["float64"], "flops32": peaks["float32"], "bw": bw,
            "ridge": peaks["float64"] / bw, "ridge32": peaks["float32"] / bw}

Reading the implementation

• Peak compute from a large square matmul, because that is the closest any real kernel gets to the hardware ceiling. A 1024³ matmul has arithmetic intensity ~170 FLOP/byte, well above any ridge point, so it is compute-bound and measures the FMA units rather than the memory system.
• Peak bandwidth from a streaming triad over 192 MB — deliberately far larger than any cache, so the measurement is DRAM and not L3. Getting this wrong is the single most common roofline error, and the out= parameters matter: without them numpy allocates a fresh array per operation and the measurement includes allocation.

The ridge point is where the two lines cross: \(I_{\text{ridge}} = \text{peak FLOP/s} / \text{peak bytes/s}\). Above it a kernel can be compute-bound; below it, no amount of arithmetic optimisation helps.

And there is one ridge per dtype. fp32 peak here is ~4× fp64 peak — two SIMD lanes per register instead of one, plus better cache reuse per byte — so the ridge moves from ~9 to ~39 FLOP/byte. Mixing them is not a rounding error: the first version of this file measured peak in fp64, priced an fp32 workload against it, and the assembly reported kernels running at twice the speed of light.

What the numbers say

Output:

  peak bandwidth :     54.9 GB/s   (triad over 192 MB in 10.5 ms)
  dtype        peak GFLOP/s    ridge point
  float64             453.5         8.26
  float32            2085.8        38.01
  fp32 is 4.60x fp64 -- two SIMD lanes
  per register instead of one, and on this machine a little more than
  the theoretical 2x because the fp32 kernel also gets better cache
  reuse per byte.

  THERE IS NO SUCH THING AS 'the' ROOFLINE. There is one per dtype, and
  the ridge moves with it: 8.3 FLOP/byte in fp64, 38.0 in fp32. Mixing
  them is not a rounding error -- the first version of this file
  measured peak in fp64 and predicted an fp32 workload, and the
  assembly reported kernels running at 2x the speed of light. A ratio
  below 1.0 against a roofline is never a fast kernel; it is always a
  broken model, and it is the most useful bug the roofline can produce
  because it is impossible to rationalise away. There are exactly two
  ways to get one: the wrong ceiling (this bug) or the wrong byte count
  (a working set that never left cache -- the assembly measures that
  one too).
  Both numbers are MEASURED. The vendor's peak assumes an FMA on every
  port every cycle, which no real kernel reaches; numbers.md records
  this machine's bandwidth at 57.5-99.1 GB/s depending on working set,
  and this triad sits at the streaming end of that range.

Beyond the toy

A ratio below 1.0 against a roofline is never a fast kernel. It is always a broken model, and there are exactly two ways to break it:

1. Wrong ceiling — a dtype mismatch, or using a vendor spec sheet instead of a measurement. Nobody reaches vendor peak; the gap is workload-specific.
2. Wrong byte count — a working set that never left cache, so the model charged DRAM traffic the kernel never paid. The assembly measures this one too.

Both are invisible in absolute timings and unmissable the moment the number is divided by a bound it cannot legally cross. That is what the model is for.

Block 2 — Arithmetic intensity

Teaches: the property that decides which wall you hit

The problem. Arithmetic intensity is the single property that decides which wall a kernel hits, and it is computable on paper before writing any code.

@block(2, "Arithmetic intensity", "the property that decides which wall you hit")
def b2(s, show):
    if show:
        print(f"  fp64 ridge = {s['ridge']:.2f} FLOP/byte on this machine "
              f"(fp32 ridge = {s['ridge32']:.2f})\n")
        print(f"  {'kernel':<28}{'FLOPs':>12}{'bytes':>12}{'intensity':>12}{'bound by':>11}")
        N = 4_000_000
        cases = [
            ("vector add (a+b)",        N,        3*N*8),
            ("scale (a*2)",             N,        2*N*8),
            ("dot product",            2*N,       2*N*8),
            ("matmul 64x64",       2*64**3,    3*64*64*8),
            ("matmul 1024x1024", 2*1024**3, 3*1024*1024*8),
            ("attention, seq 1024",  4*1024**2*64, 3*1024*64*8),
        ]
        for lbl, f, b in cases:
            ai = f / b
            print(f"  {lbl:<28}{f/1e6:>10.1f}M{b/1e6:>10.1f}M{ai:>12.2f}"
                  f"{('COMPUTE' if ai > s['ridge'] else 'MEMORY'):>11}")
        print("  Elementwise operations have intensity below 1 and can never be")
        print("  compute-bound on any machine built since about 1990 -- the ridge has")
        print("  been climbing for thirty years while DRAM latency barely moved.")
        print("  Matmul's intensity grows as O(n), which is the entire reason deep")
        print("  learning runs on hardware designed for it.")
    return {}

Reading the implementation

\[ I = \frac{\text{FLOPs}}{\text{bytes moved}} \]

The table is worth reading as a hierarchy rather than a list:

• Elementwise (a+b): 1 FLOP per 3 × 8 bytes = 0.042. Memory-bound on every machine built since about 1990, and no optimisation of the arithmetic can possibly help.
• Dot product: 2 FLOPs per 16 bytes = 0.125. Same regime — and this is P02's inner loop, which is why that project is latency-bound rather than compute-bound.
• Matmul \(n \times n\): \(2n^3\) FLOPs over \(3n^2 \times 8\) bytes, so \(I = n/12\). Intensity grows linearly with \(n\), which is the entire reason deep learning runs on hardware designed for matmul.

That last line deserves emphasis. Matmul is the only common primitive whose intensity improves with size, which is why the field bent itself toward architectures expressible as large matmuls, and why attention's \(T^2 d\) term is tolerable while an equivalent amount of elementwise work would not be.

What the numbers say

Output:

  fp64 ridge = 8.26 FLOP/byte on this machine (fp32 ridge = 38.01)

  kernel                             FLOPs       bytes   intensity   bound by
  vector add (a+b)                   4.0M      96.0M        0.04     MEMORY
  scale (a*2)                        4.0M      64.0M        0.06     MEMORY
  dot product                        8.0M      64.0M        0.12     MEMORY
  matmul 64x64                       0.5M       0.1M        5.33     MEMORY
  matmul 1024x1024                2147.5M      25.2M       85.33    COMPUTE
  attention, seq 1024              268.4M       1.6M      170.67    COMPUTE
  Elementwise operations have intensity below 1 and can never be
  compute-bound on any machine built since about 1990 -- the ridge has
  been climbing for thirty years while DRAM latency barely moved.
  Matmul's intensity grows as O(n), which is the entire reason deep
  learning runs on hardware designed for it.

Beyond the toy

The byte count is the part people get wrong, because it depends on what is in cache. A matmul that streams \(3n^2\) bytes from DRAM has intensity \(n/12\); the same matmul with all operands resident in L2 has effectively infinite DRAM intensity and is bounded by L2 bandwidth instead. Hierarchical roofline models this with a separate ceiling per cache level, and it is the right tool when the single-level model gives a ratio you cannot explain.

Block 3 — Predict, then measure

Teaches: a roofline is a falsifiable claim

The problem. A model that is never checked against measurement is a decoration. This block predicts five kernels and compares — and the direction of each error is diagnostic.

@block(3, "Predict, then measure", "a roofline is a falsifiable claim")
def b3(s, show):
    if show:
        print(f"  {'kernel':<26}{'intensity':>11}{'predicted':>13}{'measured':>12}"
              f"{'ratio':>8}")
        N = 8_000_000
        a = np.random.rand(N); b = np.random.rand(N); c = np.empty(N)
        tests = []
        t = best(lambda: np.add(a, b, out=c), 5)
        tests.append(("vector add", N, 3*N*8, t))
        t = best(lambda: np.multiply(a, 2.0, out=c), 5)
        tests.append(("scale by constant", N, 2*N*8, t))
        t = best(lambda: float(a @ b), 5)
        tests.append(("dot product", 2*N, 2*N*8, t))
        for n in (256, 1024):
            A = np.random.rand(n, n); B = np.random.rand(n, n)
            t = best(lambda: A @ B, 5)
            tests.append((f"matmul {n}x{n}", 2*n**3, 3*n*n*8, t))
        for lbl, f, byt, t in tests:
            ai = f / byt
            pred = min(s["flops"], ai * s["bw"])          # the roofline itself
            meas = f / t
            print(f"  {lbl:<26}{ai:>11.2f}{pred/1e9:>11.1f}G{meas/1e9:>10.1f}G"
                  f"{meas/pred:>8.2f}")
        print("  A ratio near 1.0 means the roofline explained the kernel. Below 1.0")
        print("  means something else is the limit -- latency, a missing")
        print("  vectorisation, an unaligned access. Above 1.0 means a modelling")
        print("  error: usually the kernel hit cache and never touched DRAM, so the")
        print("  'bytes' figure is fiction. The model earns trust by being wrong in")
        print("  ways you can explain.")
    return {}

Reading the implementation

For each kernel: count FLOPs and bytes on paper, compute intensity, take \(\min(\text{peak}, I \times \text{bandwidth})\) as the predicted rate, measure the actual, and report the ratio.

Reading the ratios:

• ≈1.0 — the roofline explained the kernel. Nothing more to find.
• <1.0 (measured faster than the bound) — impossible; the model is wrong. Usually the byte count (cache residency) or the ceiling (dtype).
• >1.0 (measured slower) — something the roofline does not model is the limit: latency rather than bandwidth, poor vectorisation, unaligned access, or insufficient memory-level parallelism.

That third case is the common one and it is where the model earns its keep — it tells you there is another effect, which is far more actionable than a raw timing.

What the numbers say

Output:

  kernel                      intensity    predicted    measured   ratio
  vector add                       0.04        2.3G       2.8G    1.23
  scale by constant                0.06        3.4G       5.4G    1.58
  dot product                      0.12        6.9G       8.0G    1.16
  matmul 256x256                  21.33      453.5G     331.8G    0.73
  matmul 1024x1024                85.33      453.5G     463.7G    1.02
  A ratio near 1.0 means the roofline explained the kernel. Below 1.0
  means something else is the limit -- latency, a missing
  vectorisation, an unaligned access. Above 1.0 means a modelling
  error: usually the kernel hit cache and never touched DRAM, so the
  'bytes' figure is fiction. The model earns trust by being wrong in
  ways you can explain.

Beyond the toy

Little's Law is the missing half of the model: \(L = \lambda W\). To sustain 50 GB/s at 64 bytes per cache line with 121 ns latency requires \(50\text{e}9/64 \times 121\text{e-}9 \approx 94\) concurrent line fetches. A core has ~10--16 line-fill buffers, so a single thread physically cannot reach peak bandwidth — which is why single-threaded STREAM numbers are always far below the machine's rating, and why memory-level parallelism (multiple threads, software prefetch, independent load chains) is the lever rather than a faster loop.

Block 4 — Tiling

Teaches: the same FLOPs, a different traffic pattern

The problem. Same FLOPs, different traffic pattern. Tiling is the canonical demonstration that how you touch memory dominates how much arithmetic you do.

@block(4, "Tiling", "the same FLOPs, a different traffic pattern")
def b4(s, show):
    def naive(A, B):
        n = A.shape[0]; C = np.zeros((n, n))
        for i in range(n):
            for k in range(n):
                C[i] += A[i, k] * B[k]
        return C
    def tiled(A, B, T=64):
        n = A.shape[0]; C = np.zeros((n, n))
        for i0 in range(0, n, T):
            for k0 in range(0, n, T):
                for j0 in range(0, n, T):
                    C[i0:i0+T, j0:j0+T] += (A[i0:i0+T, k0:k0+T]
                                            @ B[k0:k0+T, j0:j0+T])
        return C
    if show:
        n = 512
        A = np.random.rand(n, n); B = np.random.rand(n, n)
        ref = A @ B
        t_bl = best(lambda: A @ B, 3)
        t_na = best(lambda: naive(A, B), 1)
        rows = [("row-at-a-time (python loop)", t_na, naive(A, B))]
        for T in (32, 64, 128, 256):
            tt = best(lambda: tiled(A, B, T), 3)
            rows.append((f"tiled, T={T}", tt, tiled(A, B, T)))
        rows.append(("numpy (BLAS)", t_bl, ref))
        print(f"  C = A @ B, {n}x{n}, {2*n**3/1e9:.2f} GFLOP")
        print(f"  {'implementation':<30}{'time':>10}{'GFLOP/s':>10}"
              f"{'vs BLAS':>10}{'correct':>9}")
        for lbl, tt, out in rows:
            print(f"  {lbl:<30}{tt*1000:>8.1f}ms{2*n**3/tt/1e9:>10.1f}"
                  f"{t_bl/tt:>9.2f}x"
                  f"{str(bool(np.allclose(out, ref))):>9}")
        wss = lambda T: 3 * T * T * 8 / 1024
        print(f"  working set per tile: T=32 -> {wss(32):.0f} KB, "
              f"T=64 -> {wss(64):.0f} KB, T=128 -> {wss(128):.0f} KB, "
              f"T=256 -> {wss(256):.0f} KB")
        print("  The tiles here still call BLAS, so this measures BLOCKING, not")
        print("  hand-written inner loops: how much you lose by cutting a big matmul")
        print("  into small ones. The loss is real and it comes from per-call")
        print("  overhead plus reduced reuse -- which is the same trade the tile-size")
        print("  choice makes inside a real GEMM, one level down.")
    return {"tiled": tiled}

Reading the implementation

Blocked matmul: partition into \(T \times T\) tiles so each tile's working set — three \(T \times T\) blocks — fits in a cache level. Each element loaded is then reused \(T\) times before eviction, so DRAM traffic drops from \(O(n^3)\) to \(O(n^3/T)\).

The working-set arithmetic is the design rule: \(3T^2 \times 8\) bytes must fit the target level. For a 32 KB L1 that is \(T \approx 36\); for a 1 MB L2, \(T \approx 209\). Real GEMMs tile at three levels simultaneously — register, L1, L2 — which is why a hand-written microkernel is a serious piece of engineering and why the measured gap to BLAS here is unsurprising.

Be honest about what this measures. The tiles still call BLAS, so this is measuring blocking overhead — how much you lose by cutting a large matmul into small ones — rather than hand-written inner loops. That loss is real (per-call overhead plus reduced reuse) and it is the same trade a real GEMM makes one level down.

What the numbers say

Output:

  C = A @ B, 512x512, 0.27 GFLOP
  implementation                      time   GFLOP/s   vs BLAS  correct
  row-at-a-time (python loop)      297.1ms       0.9     0.00x     True
  tiled, T=32                       18.8ms      14.3     0.03x     True
  tiled, T=64                        6.1ms      44.0     0.08x     True
  tiled, T=128                       2.9ms      93.4     0.18x     True
  tiled, T=256                       1.5ms     181.1     0.34x     True
  numpy (BLAS)                       0.5ms     527.4     1.00x     True
  working set per tile: T=32 -> 24 KB, T=64 -> 96 KB, T=128 -> 384 KB, T=256 -> 1536 KB
  The tiles here still call BLAS, so this measures BLOCKING, not
  hand-written inner loops: how much you lose by cutting a big matmul
  into small ones. The loss is real and it comes from per-call
  overhead plus reduced reuse -- which is the same trade the tile-size
  choice makes inside a real GEMM, one level down.

Beyond the toy

• The theoretical floor. Hong & Kung's red-blue pebble game gives \(\Omega(n^3/\sqrt{M})\) as the minimum DRAM traffic for matmul with \(M\) words of fast memory. Tiling with \(T = \sqrt{M/3}\) achieves it up to constants, so this is not a heuristic — it is optimal.
• Cache-oblivious algorithms get the same asymptotic traffic by recursive subdivision, without knowing the cache size, which matters when you cannot tune per machine.
• What a real GEMM adds: register blocking (accumulate in registers, not memory), packing (copy tiles into contiguous, aligned buffers so the inner loop streams), explicit SIMD, software prefetch, and multithreading with a NUMA-aware partition. Goto & van de Geijn's paper is the canonical description, and the gap between this block and BLAS is exactly that list.

Block 5 — Quantisation

Teaches: fewer bytes per weight IS higher arithmetic intensity

The problem. Fewer bytes per weight is higher arithmetic intensity. In a memory-bound regime that is the most direct lever available — and this block measures the honest version, including the part numpy cannot deliver.

@block(5, "Quantisation", "fewer bytes per weight IS higher arithmetic intensity")
def b5(s, show):
    if show:
        rng = np.random.default_rng(5)
        n = 2048
        W = rng.normal(0, 0.5, (n, n)).astype(np.float32)
        x = rng.normal(0, 1, n).astype(np.float32)
        scale = np.abs(W).max() / 127.0
        Wq = np.clip(np.round(W / scale), -127, 127).astype(np.int8)
        ref = W @ x
        deq = (Wq.astype(np.float32) * scale) @ x
        t32 = best(lambda: W @ x, 20)
        t8 = best(lambda: (Wq.astype(np.float32) * scale) @ x, 20)
        t8b = best(lambda: Wq.astype(np.float32) @ x, 20)
        err = np.abs(deq - ref).max() / np.abs(ref).max()
        cos = float(deq @ ref / (np.linalg.norm(deq) * np.linalg.norm(ref)))
        print(f"  {n}x{n} weight matrix, matrix-vector product (the decode shape)")
        print(f"  {'precision':<20}{'weight bytes':>14}{'time':>10}{'GB/s':>9}"
              f"{'rel error':>12}")
        for lbl, byts, tt, e in (("float32", W.nbytes, t32, 0.0),
                                 ("int8 + dequant", Wq.nbytes, t8, err)):
            print(f"  {lbl:<20}{byts/1e6:>12.1f}M{tt*1e6:>9.0f}us"
                  f"{byts/tt/1e9:>9.1f}{e:>12.2e}")
        print(f"  cosine similarity of the two outputs: {cos:.6f}")
        print(f"  4x fewer weight bytes; measured speedup {t32/t8:.2f}x, and the")
        print("  dequantisation itself costs most of what the smaller load saved.")
        print("  A real int8 kernel keeps the arithmetic in int8 and dequantises the")
        print("  ACCUMULATOR once, which is why production quantisation needs kernel")
        print("  support and not just a smaller dtype in memory. Numpy has no int8")
        print("  GEMM, so what this block honestly measures is the memory saving and")
        print("  the accuracy cost -- both real, and the speedup is the part you")
        print("  cannot get without writing the kernel.")
    return {}

Reading the implementation

Symmetric per-tensor int8 quantisation: \(s = \max|W|/127\), then \(W_q = \text{round}(W/s)\). Dequantise as \(W_q \times s\).

The measurement is deliberately honest about its limit. numpy has no int8 GEMM, so Wq.astype(float32) * scale @ x converts back to float before multiplying — which means what this block truthfully measures is the memory saving and the accuracy cost, not the speedup. The speedup requires a kernel that keeps the arithmetic in int8 and dequantises the accumulator once, and saying so is more useful than reporting a number that comes from a different mechanism.

The accuracy result is the transferable part: cosine similarity ~0.9999 at 4× fewer bytes. Distance and dot-product rankings are remarkably robust to quantisation, which is the same premise underlying P02's PQ codes.

What the numbers say

Output:

  2048x2048 weight matrix, matrix-vector product (the decode shape)
  precision             weight bytes      time     GB/s   rel error
  float32                     16.8M      121us    138.7    0.00e+00
  int8 + dequant               4.2M     1689us      2.5    8.45e-03
  cosine similarity of the two outputs: 0.999941
  4x fewer weight bytes; measured speedup 0.07x, and the
  dequantisation itself costs most of what the smaller load saved.
  A real int8 kernel keeps the arithmetic in int8 and dequantises the
  ACCUMULATOR once, which is why production quantisation needs kernel
  support and not just a smaller dtype in memory. Numpy has no int8
  GEMM, so what this block honestly measures is the memory saving and
  the accuracy cost -- both real, and the speedup is the part you
  cannot get without writing the kernel.

Beyond the toy

• Per-tensor vs per-channel. One scale for the whole matrix is simple and loses accuracy when channels have very different ranges. Per-channel (one scale per output row) is standard and nearly free.
• Outliers are the hard part in LLM quantisation: a handful of activation channels have magnitudes 100× the rest, and they dominate the scale. LLM.int8() keeps those channels in fp16; SmoothQuant migrates the difficulty from activations into weights; GPTQ and AWQ use calibration data to choose scales that minimise output error rather than weight error.
• The formats. int8 (per-channel, well understood), int4 (needs group-wise scales, ~3.5 bits effective), fp8 e4m3/e5m2 (hardware support on H100-class, keeps dynamic range), and 1.58-bit ternary schemes at the research edge.
• Quantisation is a bandwidth optimisation first. In memory-bound decode it buys nearly its full ratio; in compute-bound prefill it buys only what the tensor-core throughput ratio gives.

Block 6 — Batching and Little's Law

Teaches: the only free speedup in a memory-bound regime

The problem. The only free speedup in a memory-bound regime. Weights are read once regardless of batch size, so every additional item in the batch is nearly free until the kernel becomes compute-bound.

@block(6, "Batching and Little's Law", "the only free speedup in a memory-bound regime")
def b6(s, show):
    if show:
        rng = np.random.default_rng(6)
        n = 2048
        W = rng.normal(0, .5, (n, n)).astype(np.float32)
        print(f"  one {n}x{n} weight matrix, batch of B vectors:")
        print(f"  {'batch':>7}{'time':>10}{'per-item':>11}{'GFLOP/s':>10}"
              f"{'intensity':>11}{'weight reads':>14}")
        t1 = None
        for B in (1, 2, 8, 32, 128):
            X = rng.normal(0, 1, (n, B)).astype(np.float32)
            t = best(lambda: W @ X, 10)
            t1 = t1 or t
            f = 2 * n * n * B
            byt = W.nbytes + X.nbytes + 4 * n * B
            print(f"  {B:>7}{t*1e6:>8.0f}us{t*1e6/B:>10.1f}us"
                  f"{f/t/1e9:>10.1f}{f/byt:>11.2f}{'1':>14}")
        print("  The weight matrix is read ONCE regardless of batch size, so every")
        print("  extra request in the batch is nearly free until the kernel becomes")
        print("  compute-bound. That is why LLM serving batches aggressively and why")
        print("  batch-1 latency is the worst possible operating point: you pay the")
        print("  full 16 MB weight read to produce a single token.")
        print("  Little's Law gives the other half: L = lambda x W. To keep a batch of")
        print("  32 in flight at 20 ms per batch you need 1600 requests/second of")
        print("  arrival. Below that, the batch never fills and you are choosing")
        print("  between latency and utilisation -- which is what a scheduler's")
        print("  max-wait parameter actually configures.")
    return {}

Reading the implementation

One weight matrix, batch of \(B\) vectors. FLOPs scale as \(B\); bytes scale as \(W + Bx\), which is dominated by \(W\) for small \(B\). So intensity scales almost linearly with batch until the weight term stops dominating.

That is the entire economics of inference serving, and it is why batch-1 latency is the worst possible operating point: you pay the full weight read to produce a single result.

What the numbers say

Output:

  one 2048x2048 weight matrix, batch of B vectors:
    batch      time   per-item   GFLOP/s  intensity  weight reads
        1     123us     122.9us      68.3       0.50             1
        2     532us     265.8us      31.6       1.00             1
        8     534us      66.7us     125.8       3.97             1
       32     253us       7.9us    1060.8      15.52             1
      128     458us       3.6us    2342.5      56.89             1
  The weight matrix is read ONCE regardless of batch size, so every
  extra request in the batch is nearly free until the kernel becomes
  compute-bound. That is why LLM serving batches aggressively and why
  batch-1 latency is the worst possible operating point: you pay the
  full 16 MB weight read to produce a single token.
  Little's Law gives the other half: L = lambda x W. To keep a batch of
  32 in flight at 20 ms per batch you need 1600 requests/second of
  arrival. Below that, the batch never fills and you are choosing
  between latency and utilisation -- which is what a scheduler's
  max-wait parameter actually configures.

Beyond the toy

Little's Law supplies the other half. \(L = \lambda W\): to keep a batch of 32 in flight at 20 ms per batch you need 1,600 requests/second of arrival. Below that, the batch never fills and you are choosing between latency (dispatch a partial batch) and utilisation (wait) — which is exactly what a serving system's max_wait_ms parameter configures.

The refinements that follow, all of which exist to raise decode intensity:

• Continuous batching — admit new requests into the running batch at every decode step rather than waiting for the whole batch to finish. Removes head-of-line blocking from a long generation.
• Paged attention — allocate KV cache in fixed blocks like OS pages (P12) so memory is not over-provisioned per sequence, which raises the batch size that fits.
• Speculative decoding — a draft model proposes \(k\) tokens, the target verifies them in one pass, converting \(k\) memory-bound steps into one compute-bound one with identical output distribution.

Block 7 — The decode wall

Teaches: why generation is memory-bound and prefill is not

The problem. Why generation is memory-bound and prefill is not. This is the single most consequential arithmetic in LLM serving, and it is one division.

@block(7, "The decode wall", "why generation is memory-bound and prefill is not")
def b7(s, show):
    if show:
        print("  A transformer layer, hidden d, batch B, sequence S. Weights are")
        print("  ~12d^2 bytes in fp16; the matmuls are ~24 B S d^2 FLOPs.")
        print(f"  {'phase':<22}{'B':>4}{'S':>7}{'intensity':>12}{'bound by':>11}"
              f"{'note':>22}")
        d = 4096
        for lbl, B, S in (("prefill, 2k prompt", 1, 2048), ("decode, 1 token", 1, 1),
                          ("decode, batch 32", 32, 1), ("decode, batch 256", 256, 1)):
            flops = 24 * B * S * d * d
            byts = 12 * d * d + 4 * B * S * d
            ai = flops / byts
            note = "reads 200MB per token" if B == 1 and S == 1 else ""
            print(f"  {lbl:<22}{B:>4}{S:>7}{ai:>12.1f}"
                  f"{('COMPUTE' if ai > s['ridge32'] else 'MEMORY'):>11}{note:>22}")
        print(f"  (fp32 ridge on this machine = {s['ridge32']:.1f} FLOP/byte; on an H100")
        print("  with ~990 TFLOP/s and ~3.35 TB/s it is about 295, so the same table")
        print("  on a GPU pushes even batch-256 decode into the memory-bound column.)")
        print("  Prefill has S=2048 tokens sharing one weight read, so it is compute-")
        print("  bound and scales with FLOPs. Decode has S=1: the SAME weights are")
        print("  read to produce a single token. Batching is the only lever that")
        print("  raises decode intensity, which is the whole reason continuous")
        print("  batching, paged attention and speculative decoding exist -- all")
        print("  three are attempts to get more work per weight read. See proofs.md")
        print("  P6 for the derivation.")
    return {}

Reading the implementation

For a transformer layer with hidden size \(d\), batch \(B\), sequence \(S\): weights are ~\(12d^2\) bytes in fp16 and the matmuls are ~\(24BSd^2\) FLOPs, so

\[ I = \frac{24BSd^2}{12d^2 + 4BSd} \approx 2BS \quad \text{for } BSd \ll 3d^2 \]

Prefill has \(S = 2048\): thousands of tokens share one weight read, so intensity is high and the phase is compute-bound. Decode has \(S = 1\): the same weights are read to produce one token, intensity ≈2 FLOP/byte, and the phase is memory-bound by two orders of magnitude.

Batching is the only lever that raises decode intensity, because it is the only term in the numerator you control.

What the numbers say

Output:

  A transformer layer, hidden d, batch B, sequence S. Weights are
  ~12d^2 bytes in fp16; the matmuls are ~24 B S d^2 FLOPs.
  phase                    B      S   intensity   bound by                  note
  prefill, 2k prompt       1   2048      3510.9    COMPUTE                      
  decode, 1 token          1      1         2.0     MEMORY reads 200MB per token
  decode, batch 32        32      1        63.8    COMPUTE                      
  decode, batch 256      256      1       501.6    COMPUTE                      
  (fp32 ridge on this machine = 38.0 FLOP/byte; on an H100
  with ~990 TFLOP/s and ~3.35 TB/s it is about 295, so the same table
  on a GPU pushes even batch-256 decode into the memory-bound column.)
  Prefill has S=2048 tokens sharing one weight read, so it is compute-
  bound and scales with FLOPs. Decode has S=1: the SAME weights are
  read to produce a single token. Batching is the only lever that
  raises decode intensity, which is the whole reason continuous
  batching, paged attention and speculative decoding exist -- all
  three are attempts to get more work per weight read. See proofs.md
  P6 for the derivation.

Beyond the toy

Concrete consequences, using published H100 figures (~3.35 TB/s):

No arithmetic optimisation moves those numbers. Every real lever is a byte lever: quantise the weights, share KV heads (GQA/MQA), or amortise the read across more sequences.

And the ridge points make it worse over time: this machine's fp32 ridge is ~39, an A100's is ~156, an H100's ~295. Every accelerator generation has widened the gap between compute and bandwidth, so a kernel that was compute-bound on a V100 can be memory-bound on an H100 with no code change. That is the strongest argument in the whole track for re-measuring rather than inheriting conclusions.

The assembly

Every block above, wired together into one working system:

def assembly(s):
    print("\nSeven blocks = a performance model. Predict a workload before running it.\n")
    rng = np.random.default_rng(14)
    d, L = 1024, 6
    Ws = [rng.normal(0, .02, (d, d)).astype(np.float32) for _ in range(L)]
    def layer(x, W): return np.maximum(x @ W, 0)
    print(f"  Workload: a {L}-layer MLP, width {d}, fp32, batch B.")
    print(f"  Per batch: {L} matmuls = {2*L*d*d/1e6:.1f} MFLOP per item,")
    print(f"  weights = {L*d*d*4/1e6:.1f} MB read once per batch.\n")
    print(f"  {'batch':>6}{'intensity':>11}{'roofline says':>15}{'predicted':>12}"
          f"{'measured':>11}{'ratio':>8}{'bound':>9}")
    ratio_b1 = None
    for B in (1, 4, 16, 64, 256):
        X = rng.normal(0, 1, (B, d)).astype(np.float32)
        flops = 2 * L * B * d * d
        byts = L * d * d * 4 + 2 * B * d * 4 * L
        ai = flops / byts
        pred_rate = min(s["flops32"], ai * s["bw"])      # fp32 workload -> fp32 peak
        pred_t = flops / pred_rate
        def run():
            x = X
            for W in Ws: x = layer(x, W)
            return x
        t = best(run, 5)
        bound = "COMPUTE" if ai > s["ridge32"] else "MEMORY"
        if B == 1: ratio_b1 = t / pred_t
        print(f"  {B:>6}{ai:>11.2f}{bound:>15}{pred_t*1e6:>10.0f}us"
              f"{t*1e6:>9.0f}us{t/pred_t:>8.2f}{bound:>9}")
    print("\n  Now hold the batch at 1 and grow the weights instead, so the same")
    print("  model is priced against working sets that do and do not fit in cache:\n")
    print(f"  {'width':>7}{'layers':>8}{'weight MB':>12}{'predicted':>12}"
          f"{'measured':>11}{'ratio':>8}")
    fp_ratios = []
    for dd, LL in ((1024, 6), (2048, 6), (2048, 12), (4096, 8)):
        Ws2 = [rng.normal(0, .02, (dd, dd)).astype(np.float32) for _ in range(LL)]
        X = rng.normal(0, 1, (1, dd)).astype(np.float32)
        byts = LL * dd * dd * 4 + 2 * dd * 4 * LL
        def run2():
            v = X
            for W in Ws2: v = np.maximum(v @ W, 0)
            return v
        tt = best(run2, 5); pred = byts / s["bw"]
        fp_ratios.append(tt / pred)
        print(f"  {dd:>7}{LL:>8}{LL*dd*dd*4/1e6:>12.1f}{pred*1e6:>10.0f}us"
              f"{tt*1e6:>9.0f}us{tt/pred:>8.2f}")
    print(f"\n  Look at the 25 MB row twice. The batch table measured it at "
          f"{ratio_b1:.2f} and")
    print(f"  the footprint table at {fp_ratios[0]:.2f} -- the same weights, the same")
    print("  arithmetic, both sitting AT the bound rather than comfortably above it,")
    print(f"  and wandering by {abs(ratio_b1-fp_ratios[0]):.2f} between two runs in the same process.")
    print("  A ratio that hovers at or below 1.0 is the signature of a PARTLY cache-")
    print("  resident working set: the weights are re-read every iteration, some")
    print("  fraction survives in this machine's last-level cache, and how large that")
    print("  fraction is depends on what else touched memory first. The model charged")
    print("  the kernel for 25 MB of DRAM traffic that it only partly paid.")
    print(f"  Past 100 MB the ambiguity disappears: {fp_ratios[1]:.2f}, {fp_ratios[2]:.2f}, "
          f"{fp_ratios[3]:.2f}. The working")
    print("  set no longer fits, every byte really does come from DRAM, and the")
    print("  roofline becomes a bound the kernel reaches about half of. The cliff")
    print("  between those two regimes is the cache, measured without ever naming")
    print("  its size -- and it is the same cliff as P12's page-fault curve and")
    print("  P04's Bloom filter, one level up the hierarchy.")
    print("\n  Two different bugs produced a sub-1.0 ratio in this file: an fp64")
    print("  ceiling on an fp32 workload (block 1, a factor of 4) and a byte count")
    print("  that assumed DRAM for data sitting in cache (above). Both were invisible in the")
    print("  absolute timings and both were obvious the moment the number was divided")
    print("  by a bound it could not legally cross. That is what the model is FOR --")
    print("  not predicting runtime, but making a specific class of mistake loud.")
    print("\n  The prediction uses two numbers measured in block 1 and a FLOP count")
    print("  done on paper. No profiler, no counters. Read the batch table's gap")
    print("  where it is largest: at batch 256 we reach roughly half of peak, because")
    print("  a 1024x256 matmul is skinnier than the square one that set the ceiling")
    print("  and there is a full-size ReLU pass between every layer.")
    print("\n  This is the deliverable of hardware-aware ML: not a faster kernel, but")
    print("  the ability to say IN ADVANCE which optimisations can possibly help.")
    print("  If a kernel sits at 0.9 of its roofline, rewriting the inner loop is")
    print("  wasted work and the only remaining moves are algorithmic -- fewer bytes")
    print("  (quantisation, block 5) or more work per byte (batching, block 6).")
    print("\n  Built: measured roofline -> arithmetic intensity -> prediction vs")
    print("  measurement -> tiling -> quantisation -> batching -> the decode wall.")
    print("  Missing, on the project page: hardware counters via perf (m3), a real")
    print("  hand-written GEMM microkernel with register blocking (m5), operator")
    print("  fusion measured end to end (m7), GPU occupancy and warp scheduling")
    print("  (m9-m10), and E2 -- the experiment that finds this machine's cache")
    print("  hierarchy from a latency curve rather than from a spec sheet.")

Output:

Seven blocks = a performance model. Predict a workload before running it.

  Workload: a 6-layer MLP, width 1024, fp32, batch B.
  Per batch: 6 matmuls = 12.6 MFLOP per item,
  weights = 25.2 MB read once per batch.

   batch  intensity  roofline says   predicted   measured   ratio    bound
       1       0.50         MEMORY       460us      424us    0.92   MEMORY
       4       1.98         MEMORY       462us      501us    1.08   MEMORY
      16       7.76         MEMORY       473us      553us    1.17   MEMORY
      64      28.44         MEMORY       516us     1238us    2.40   MEMORY
     256      85.33        COMPUTE      1544us     3398us    2.20  COMPUTE

  Now hold the batch at 1 and grow the weights instead, so the same
  model is priced against working sets that do and do not fit in cache:

    width  layers   weight MB   predicted   measured   ratio
     1024       6        25.2       460us      551us    1.20
     2048       6       100.7      1836us     3711us    2.02
     2048      12       201.3      3672us     7625us    2.08
     4096       8       536.9      9788us    27413us    2.80

  Look at the 25 MB row twice. The batch table measured it at 0.92 and
  the footprint table at 1.20 -- the same weights, the same
  arithmetic, both sitting AT the bound rather than comfortably above it,
  and wandering by 0.28 between two runs in the same process.
  A ratio that hovers at or below 1.0 is the signature of a PARTLY cache-
  resident working set: the weights are re-read every iteration, some
  fraction survives in this machine's last-level cache, and how large that
  fraction is depends on what else touched memory first. The model charged
  the kernel for 25 MB of DRAM traffic that it only partly paid.
  Past 100 MB the ambiguity disappears: 2.02, 2.08, 2.80. The working
  set no longer fits, every byte really does come from DRAM, and the
  roofline becomes a bound the kernel reaches about half of. The cliff
  between those two regimes is the cache, measured without ever naming
  its size -- and it is the same cliff as P12's page-fault curve and
  P04's Bloom filter, one level up the hierarchy.

  Two different bugs produced a sub-1.0 ratio in this file: an fp64
  ceiling on an fp32 workload (block 1, a factor of 4) and a byte count
  that assumed DRAM for data sitting in cache (above). Both were invisible in the
  absolute timings and both were obvious the moment the number was divided
  by a bound it could not legally cross. That is what the model is FOR --
  not predicting runtime, but making a specific class of mistake loud.

  The prediction uses two numbers measured in block 1 and a FLOP count
  done on paper. No profiler, no counters. Read the batch table's gap
  where it is largest: at batch 256 we reach roughly half of peak, because
  a 1024x256 matmul is skinnier than the square one that set the ceiling
  and there is a full-size ReLU pass between every layer.

  This is the deliverable of hardware-aware ML: not a faster kernel, but
  the ability to say IN ADVANCE which optimisations can possibly help.
  If a kernel sits at 0.9 of its roofline, rewriting the inner loop is
  wasted work and the only remaining moves are algorithmic -- fewer bytes
  (quantisation, block 5) or more work per byte (batching, block 6).

  Built: measured roofline -> arithmetic intensity -> prediction vs
  measurement -> tiling -> quantisation -> batching -> the decode wall.
  Missing, on the project page: hardware counters via perf (m3), a real
  hand-written GEMM microkernel with register blocking (m5), operator
  fusion measured end to end (m7), GPU occupancy and warp scheduling
  (m9-m10), and E2 -- the experiment that finds this machine's cache
  hierarchy from a latency curve rather than from a spec sheet.

The design space

The roofline is the simplest useful performance model, and its value is that it is falsifiable. Richer models exist and each buys a specific kind of accuracy.

Little's Law deserves equal billing: \(L = \lambda W\). To sustain \(\lambda\) requests per second at latency \(W\), you need \(L\) in flight. With DRAM at 121 ns and a need for 50 GB/s at 64 B per line, you need \(50\text{e}9 / 64 \times 121\text{e-}9 \approx 94\) concurrent line fetches. A core with ~10--16 line-fill buffers cannot reach that alone — which is why single-threaded bandwidth is always far below peak, and why the memory-level parallelism argument in P02 matters more than any arithmetic tweak.

Two modelling bugs this page catches

Both produced a ratio below 1.0 — measured faster than the bound — which is physically impossible and therefore diagnostic.

1. Wrong ceiling (dtype). fp32 peak is ~4× fp64 peak on this machine (1937 vs 469 GFLOP/s measured), so pricing an fp32 workload against an fp64 roofline reports kernels running at twice the speed of light. There is no single roofline for a machine; there is one per dtype, and the ridge moves with it (9.4 vs 38.8 FLOP/byte here).
2. Wrong byte count (cache residency). A 25 MB working set re-read every iteration stays in last-level cache, so the model charges DRAM traffic the kernel never paid. The ratio then wanders across 1.0 from run to run. Past 100 MB the ambiguity disappears and the ratio settles near 2.

The general rule: a sub-1.0 ratio is never a fast kernel, it is a broken model — and the two ways to break it are the ceiling and the byte count.

The hardware, level by level

CPU

GPU

The mental model is occupancy and coalescing, not cores. A warp of 32 threads issues one memory transaction if their addresses fall in the same lines (coalesced) and up to 32 if they do not. Shared memory is a programmer-managed scratchpad (up to 228 KB/SM on H100) that plays the role of L1 tiling in a CPU GEMM. Tensor cores execute small matrix ops (mma.m16n8k16) and are the only way to reach quoted peak — a hand-written FMA loop tops out around 5% of it.

TPU

A 128×128 systolic MXU with compiler-scheduled data movement and no cache hierarchy to speak of. Each value read from memory is reused 128 times inside the array — the hardware expression of the reuse argument in proofs.md P15. Consequences: matmul dimensions that are not multiples of 128 waste the array, and the compiler (XLA) must know the shapes statically, which is why dynamic shapes are painful on TPU and fine on GPU.

Comparative ridges

Every accelerator generation has made the memory wall worse, because FLOP/s grows faster than bandwidth. A kernel that was compute-bound on a V100 can be memory-bound on an H100 with no code change — which is the strongest possible argument for re-measuring the model rather than inheriting conclusions.

Optimisation levers, in the order they pay

1. Fewer bytes. Quantisation (fp16 → int8 → int4) raises intensity directly. Block 5 measures the honest version: numpy has no int8 GEMM, so what is demonstrated is the memory saving and the accuracy cost, and the speedup is the part that requires writing the kernel.
2. More work per byte. Batching (P01, block 6 here) is the only lever that raises decode intensity, and it is why continuous batching, paged attention and speculative decoding all exist.
3. Better reuse. Tiling/blocking so the working set fits a cache level. Block 4 measures blocking specifically — how much you lose by cutting a big matmul into small ones — which is the same trade a real GEMM makes one level down with register blocking.
4. Fusion. Remove intermediate round trips (P13).
5. Only then, the inner loop. If a kernel sits at 0.9 of its roofline, rewriting it is wasted work, and the remaining moves are all algorithmic.

That ordering is the deliverable of hardware-aware ML: not a faster kernel, but the ability to say in advance which optimisations can possibly help.

Advanced topics

• Sparsity: 2:4 structured sparsity is supported in hardware (2× on tensor cores) because unstructured sparsity's irregular access destroys coalescing — a case where the hardware dictates the algorithm's shape.
• Communication-avoiding algorithms: the \(\Omega(n^3/\sqrt{M})\) lower bound on matmul memory traffic (Hong–Kung) is the theoretical floor tiling approaches; 2.5D and CARMA algorithms extend it to distributed memory.
• Arithmetic intensity of the whole pipeline, not one kernel — a fused attention block has very different intensity from its parts, which is exactly FlashAttention's argument.
• Amdahl and Gustafson bound the whole exercise: with 95% of runtime parallelised, the maximum speedup is 20× no matter the machine.
• Energy is increasingly the real constraint: a DRAM access costs ~100× the energy of an fp32 add, so intensity is also an energy metric.

How this connects to the rest of the track

• P01's decode wall is this model applied to a transformer.
• P13 supplies the kernels and the fusion opportunities.
• P02 is the memory-bound extreme — 0.25 FLOP/byte, latency-bound rather than bandwidth-bound.
• P12 is the same hierarchy one level down, where the miss goes to disk instead of DRAM.
• P15 uses this model to put a hardware ceiling beside every measured QPS number.

Failure modes at scale

• Benchmarking the cache rather than the memory system — the second bug above, and the reason numbers.md §14 exists.
• Peak from a spec sheet instead of a measurement; nobody reaches vendor peak, and the gap is workload-specific.
• Ignoring the dtype, as above.
• Measuring a constant-folded loop: the compiler deletes work that has no observable effect. This track hit it twice — a timed loop removed entirely, and a Rust benchmark reporting 0.0000 ms until black_box was added.
• Optimising the wrong kernel because the profile was taken at the wrong batch size; intensity, and therefore which wall you hit, changes with batch.

Primary sources

• Williams, Waterman & Patterson, Roofline: An Insightful Visual Performance Model (CACM 2009).
• Hong & Kung, I/O Complexity: The Red-Blue Pebble Game (STOC 1981).
• Little, A Proof for the Queuing Formula \(L = \lambda W\) (1961).
• Goto & van de Geijn, Anatomy of High-Performance Matrix Multiplication (TOMS 2008) — what block 4 is a shadow of.
• Jouppi et al., In-Datacenter Performance Analysis of a Tensor Processing Unit (ISCA 2017).
• Hennessy & Patterson, Computer Architecture: A Quantitative Approach (6th ed.), chapters 2 and 4.

Running it

python3 handson/h14_hardware.py            # every block, then the assembly
python3 handson/h14_hardware.py --block 3  # just block 3 and its prerequisites
python3 handson/h14_hardware.py --quiet    # the assembly only

What to do with this

Apply it to something you already run. Count the FLOPs and the bytes on paper, measure the two machine constants once, and predict the runtime before you measure it. The value is not the prediction --- it is that a ratio far from 1.0 tells you which of your assumptions is wrong, and it does so before you have spent a week optimising a kernel that was already at its bound.

Milestones, experiments, readings and exit criteria for this project: P14 — Hardware-Aware ML System.