Matmul Arithmetic Intensity

What arithmetic intensity is

Every kernel does two things a machine charges for: it computes, and it moves data. Arithmetic intensity is the ratio of the first to the second — the number of floating-point operations performed per byte transferred between the compute units and memory:

I = (FLOPs performed) / (bytes moved)   [FLOP/byte]

The ratio tells you which resource is the bottleneck. A machine has a peak compute rate (FLOP/s) and a peak memory bandwidth (bytes/s). If a kernel supplies operations faster than memory can feed operands, the math units sit idle waiting — the kernel is memory-bound. If operands arrive faster than the math units can consume them, the bus is idle and the kernel is compute-bound. Intensity is the coordinate that places a kernel on one side or the other, and matmul is the canonical example because its intensity is a knob you can turn with the shapes.

Counting the FLOPs: 2MNK

Take the standard product C = A · B with A: [M, K], B: [K, N], and C: [M, N]. Every output element C[i,j] is a dot product of one row of A with one column of B, both of length K:

C[i,j] = Σ_k A[i,k] * B[k,j]     (k = 1..K)

That inner sum costs K multiplies and K adds — 2K FLOPs by the universal convention that a multiply and an add are one op each. There are M × N output elements, so the whole product costs

FLOPs = 2 · M · N · K

This 2MNK is exact and worth memorizing: it is the numerator of every intensity calculation. A fused multiply-add executes the multiply and add together, but it still counts as two FLOPs — the hardware just does them in one instruction.

Counting the bytes: what actually moves

The denominator is subtler because it depends on how much data the kernel is forced to fetch from slow memory. In the ideal case each matrix is read (or written) exactly once: A contributes MK elements, B contributes KN, and C contributes MN written back. At p bytes per element the ideal traffic is

bytes_ideal = p · (M·K + K·N + M·N)

with p = 2 for fp16/bf16, p = 4 for fp32. This is a lower bound on traffic and therefore an upper bound on intensity: it assumes every value fetched is reused for all the arithmetic that needs it before being evicted. A naive triple loop does far worse — it re-reads B from memory for every row of A — so its real byte count is much larger and its real intensity much lower. The gap between ideal and naive is exactly the reuse that tiling exists to capture.

The intensity formula and its shape

Divide the two counts and the intensity of a matmul is

I = 2MNK / [ p · (MK + KN + MN) ]   FLOP/byte

The cubic numerator against the quadratic denominator is the whole story. For a square problem M = N = K = n it collapses to

I = 2n³ / (p · 3n²) = 2n / (3p)

Intensity grows linearly with n. Double the matrix and you double the arithmetic squeezed out of each byte. That is the property no other common kernel has: elementwise ops, normalization, and activations are stuck at a fixed, low intensity no matter how big the tensor, because each element is touched a constant number of times. Matmul is the one primitive whose intensity you can scale up simply by making it bigger, which is why hardware and models are both built around it.

Why intensity rises with size: reuse

The linear growth is a reuse argument in disguise. In C = A · B, each element of A participates in N different output dot products (once per column of C), and each element of B participates in M of them. So a single fetched value of A does 2N FLOPs of useful work over its lifetime, and a value of B does 2M. The bytes you pay for a value are fixed; the arithmetic it enables scales with the other dimension.

This is why the ‘fat’ dimensions matter. If M, N, and K are all large, every operand is reused hundreds or thousands of times and the fetch cost amortizes away. If one dimension is tiny, the operands along it are barely reused and the kernel drowns in traffic. Arithmetic intensity is just a bookkeeping of average reuse: high reuse means high intensity means compute-bound.

The roofline ridge point

Whether a given intensity is ‘enough’ depends on the machine. The roofline model draws attainable performance as min(peak_FLOPs, I × bandwidth): below a threshold intensity you are on the sloped bandwidth roof, above it you are on the flat compute roof. The crossover — the ridge point — is the machine balance:

I_ridge = peak_FLOPs / peak_bandwidth

For an H100 at roughly 990 TFLOP/s bf16 and about 3.35 TB/s of HBM bandwidth, I_ridge ≈ 990e12 / 3.35e12 ≈ 295 FLOP/byte. A kernel needs intensity above ~300 just to reach the compute roof on that chip. Setting the square-matmul intensity 2n/(3p) equal to 295 with p = 2 gives n ≈ 885: below that a perfectly tiled square GEMM cannot saturate the FLOP/s no matter how good the kernel, because there simply is not enough reuse in the problem. Bigger machines have higher ridge points, which is why matmuls keep needing to grow to stay compute-bound.

A worked example: a 4096 cube

Take M = N = K = 4096 in bf16 (p = 2). The arithmetic is

FLOPs = 2 · 4096³ ≈ 1.37e11  (137 GFLOP)

and the ideal traffic is

bytes = 2 · 3 · 4096² ≈ 1.01e8  (~100 MB)

so I ≈ 1.37e11 / 1.01e8 ≈ 1365 FLOP/byte — exactly 2n/(3p) = 4096/3. That is well above the ~295 ridge point, so this GEMM is firmly compute-bound: on an H100 the lower bound on runtime is set by compute, 137e9 / 990e12 ≈ 138 µs, not by moving those 100 MB, which would take only 100e6 / 3.35e12 ≈ 30 µs. The math units are the binding constraint, which is the regime you want: it means the hardware’s headline FLOP/s number is actually reachable.

GEMV: why it is stuck near intensity one

Now collapse N to 1 — a matrix times a vector, y = A · x with A: [M, K], x: [K], y: [M]. The FLOPs are 2MK, and the bytes are dominated by reading the matrix once, p · MK (the vectors are negligible when M, K are large). So

I ≈ 2MK / (p · MK) = 2 / p

That is 1.0 FLOP/byte in fp16 and 0.5 in fp32 — independent of size. The reason is stark: with N = 1 every matrix element is used exactly once, so there is no reuse to amortize the fetch. You pay for the whole matrix and get two FLOPs per element in return, forever. GEMV is hopelessly memory-bound on any modern chip whose ridge point is in the hundreds, which is precisely the situation in autoregressive decoding: generating one token at a time makes each weight matrix multiply a GEMV, so token-by-token inference is bandwidth-bound and the FLOP/s of the accelerator barely matter.

Batching turns GEMV back into GEMM

The escape from intensity one is to stop doing one vector at a time. Stack B vectors into a matrix X: [K, B] and the matrix-vector product becomes a matrix-matrix product Y = A · X with N = B. Now each element of the weight matrix A is reused B times, once per column, and the intensity climbs toward

I ≈ 2MKB / [ p · (MK + KB + MB) ] → grows with B

When MK dominates (the usual case for a large weight matrix), the matrix is read once but drives B times more arithmetic, so intensity scales roughly with the batch until other terms catch up. This is the entire quantitative case for batching in inference: it does not reduce the work, it raises reuse, moving the weight-load cost from being paid per token to being shared across a batch. It is why throughput-oriented serving batches aggressively while latency-bound single-stream decoding stays memory-bound.

Tiling: recovering the reuse a naive loop wastes

The ideal byte count assumed each operand is fetched once, but a naive triple loop re-streams B from memory for every row of A, inflating traffic by a factor of M and crushing the achieved intensity far below the theoretical 2n/(3p). Tiling closes that gap. Block the output into tiles of Mt × Nt, hold the accumulating C tile in registers, and stream matching strips of A and B through fast cache:

for each C-tile [Mt x Nt]:
    keep accumulators in registers
    for k-strip: load A[Mt x Kt], B[Kt x Nt] into cache, FMA

Inside the tile, each loaded A value is reused Nt times and each B value Mt times, so the effective reuse — and therefore the achieved intensity — scales with the tile dimensions. Register and cache capacity cap how large the tile can be, which in turn caps achievable intensity. Good matmul kernels are, at their core, an exercise in choosing tile sizes that push realized reuse as close to the 2n/(3p) ceiling as the memory hierarchy allows.

Precision, and why bytes-per-element matters

The factor p sits in the denominator of every intensity formula, so narrowing the datatype directly raises intensity. Going from fp32 (p = 4) to bf16 (p = 2) doubles the FLOP/byte for the same shapes, and int8 or fp8 (p = 1) doubles it again. That is a second lever alongside shape: for a memory-bound GEMV, halving the bytes per weight halves the traffic and roughly doubles achievable throughput, which is a large part of why weight quantization is so effective for single-stream decoding.

The subtlety is that lowering precision usually raises the hardware’s peak FLOP/s too, which pushes the ridge point up. So quantization helps most in the memory-bound regime, where you are on the bandwidth roof and cutting bytes is a direct win; in the compute-bound regime the benefit comes instead from the faster low-precision math units, not from the intensity change.

Implications for CPU and small-model inference

On a CPU the ridge point is lower than on a GPU — tens of GFLOP/s per core against maybe tens of GB/s per socket still lands the balance in the low tens of FLOP/byte — but the logic is identical, and the memory hierarchy is deeper and less forgiving. Small-model inference lives mostly in the GEMV regime because it runs batch-one, so it is bandwidth-bound: the win comes from quantizing weights to shrink p, keeping the working set inside the L2/L3 cache, and blocking so a weight tile loaded once serves many FLOPs before eviction.

The common pitfall is to chase peak FLOP/s numbers that a batch-one workload can never reach, or to compute intensity with the ideal byte count and then be surprised the kernel runs slower — because the realized traffic, set by cache behavior and tiling, is what the hardware actually experiences. Estimate with the ideal ratio to know the ceiling; profile the realized ratio to know how far below it you are.