Deep Implementation

Benchmark ops:byte in practice

Back on day 6, I worked through the roofline model on paper. Today I measured real ops:byte ratios for common inference kernels and watched the theory hold up, break down, and teach me things the textbook never mentioned.

On day 6, I worked through the roofline model: ops:byte ratio as a hardware constant, arithmetic intensity as the software side, and the chart that shows which resource is your bottleneck. The theory was clean. Today I wanted to see whether it actually holds up when I measure real GPU kernels.

The short answer: it does, mostly. But "mostly" is the interesting part.

The setup

I wrote a small benchmarking harness in PyTorch that measures wall-clock time for isolated operations at various sizes, then computes the effective FLOP/s and effective bandwidth. From those two numbers, I can figure out whether a kernel is compute-bound or memory-bound, and compare the result against what the roofline predicts.

import torch
import time

def bench_matmul(M, N, K, dtype=torch.float16, warmup=10, iters=100):
    A = torch.randn(M, K, dtype=dtype, device="cuda")
    B = torch.randn(K, N, dtype=dtype, device="cuda")
    # warmup
    for _ in range(warmup):
        C = A @ B
    torch.cuda.synchronize()

    start = time.perf_counter()
    for _ in range(iters):
        C = A @ B
    torch.cuda.synchronize()
    elapsed = (time.perf_counter() - start) / iters

    flops = 2 * M * N * K  # standard matmul FLOP count
    bytes_moved = (M * K + K * N + M * N) * A.element_size()
    arithmetic_intensity = flops / bytes_moved

    return {
        "elapsed_ms": elapsed * 1000,
        "tflops": flops / elapsed / 1e12,
        "bw_gb_s": bytes_moved / elapsed / 1e9,
        "arithmetic_intensity": arithmetic_intensity,
    }

Nothing fancy. The key insight is that you need both the FLOP count and the memory traffic estimate to place a kernel on the roofline. Without both, you are guessing.

Matmul: square matrices climb the roofline

I swept square matrix sizes from 128 to 8192 on an A100 80GB (FP16). The results were textbook:

The transition from memory-bound to compute-bound happens right around where the roofline says it should: when arithmetic intensity crosses the A100's ops:byte ratio of roughly 160 (for FP16 with Tensor Cores at 312 TFLOP/s and 2 TB/s bandwidth).

Decode-shaped matmul: the reality check

Square matmuls are the easy case. The interesting one is the decode-shaped matmul: a tall-skinny matrix times the weight matrix. In autoregressive decode, each step multiplies a (batch_size, hidden_dim) vector against the (hidden_dim, hidden_dim) weight. With batch size 1, that is a vector-matrix multiply.

results = bench_matmul(M=1, N=4096, K=4096, dtype=torch.float16)

Arithmetic intensity for this shape: 2 * 1 * 4096 * 4096 / ((1*4096 + 4096*4096 + 1*4096) * 2) = about 1.0 FLOP/byte. One. The kernel is loading 32 MB of weights to do 33 million FLOP. It is deeply, hopelessly memory-bound.

Measured bandwidth: about 1.6 TB/s, which is 80% of the A100's peak bandwidth. The kernel is doing its job. It is saturating the memory bus. But at batch=1, it uses less than 1% of available compute. All those Tensor Cores, idle.

Batching pulls you up the roofline

This is where the practical lesson lives. I swept batch sizes from 1 to 256:

for bs in [1, 2, 4, 8, 16, 32, 64, 128, 256]:
    r = bench_matmul(M=bs, N=4096, K=4096)
    print(f"bs={bs:4d}  AI={r['arithmetic_intensity']:.1f}  "
          f"TFLOP/s={r['tflops']:.1f}  BW={r['bw_gb_s']:.0f} GB/s")

The transition is smooth and predictable:

At batch size 256, the same matmul that was wasting 99% of compute at batch=1 is now compute-bound and achieving 250+ TFLOP/s. This is why continuous batching matters so much for serving: it is not just about throughput, it is about turning a memory-bound problem into a compute-bound one.

The practical takeaway

If you are running a decode workload at low batch sizes and wondering why your expensive GPU seems slow, it is not slow. It is starving. The weights are the bottleneck, not the math. Batch more, quantize the weights, or accept that you are paying for compute you cannot use.

Where the theory breaks down

Two things surprised me. First, the naive FLOP/byte calculation assumes memory traffic equals the sum of input and output tensors. In reality, there are L2 cache effects. Small matmuls can keep operands in L2 (40 MB on A100), which means effective bandwidth is much higher than HBM bandwidth. My 128x128 matmul was faster than the roofline predicted because it was hitting L2, not HBM.

Second, very tall-skinny shapes (like bs=1, K=4096) do not always hit peak bandwidth because the memory access pattern is less efficient. Memory controllers are optimized for large, contiguous reads. A vector-matrix multiply has inherently fragmented access patterns compared to a large square matmul. I saw 80% of peak bandwidth at bs=1 but 95% at bs=8, even though both are memory-bound.

Elementwise ops: bandwidth-bound by definition

I also measured elementwise operations (ReLU, GELU, layer norm) for completeness. These always have arithmetic intensity close to 1: you read each element, do a small amount of math, and write the result. They live firmly on the memory-bound side of the roofline, always.

This is why kernel fusion matters. If you can fuse a GELU activation into the preceding matmul, you avoid a round trip to HBM. The GELU itself is cheap. The memory traffic to do it separately is not.

What I learned

The roofline is not just a theoretical model. It is a diagnostic tool. Measure your arithmetic intensity, compare it to the hardware's ops:byte ratio, and you know exactly where to focus optimization effort. No guessing required.

The next step is to stop measuring synthetic kernels and start profiling real model inference end-to-end. That is exactly what day 43 is about: using torch.profiler to see where time actually goes inside a real forward pass.