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.