The Cost of Matrix Multiplication — Where Almost All of AI's Compute Goes

Why GEMM is everything: the anatomy of O(n³), memory bandwidth and arithmetic intensity, what actually makes a GPU fast, and the intuition behind tiling — ending with you able to estimate a model's training and inference FLOPs yourself.

The metaphor: the machine and the trucks

A factory installs a ferociously fast machine tool: it can cut a thousand parts per second. The trucks delivering raw stock, however, arrive with ten parts per second. This factory's output is set by the number of trucks, not by the machine. The machine stands idle 99% of the time.

Modern computers are that factory. Arithmetic units are ferociously fast; the bandwidth that carries data out of memory has not kept pace. Counting how many operations you perform therefore does not predict runtime. You have to count how many bytes you move at the same time.

And in AI, the overwhelming majority of that time goes into a single operation: matrix multiplication.

Why matrix multiplication is everything

Take a Transformer apart and what falls out is almost entirely matrix products: the projections to Q, K and V, the attention scores, the weighted sum over values, the output projection, and the two linear layers of the FFN. Convolutions, too, are typically lowered to a matrix product via im2col or an equivalent rearrangement.

That is why the library world gave this one operation a dedicated name — GEMM, general matrix multiply — and has been optimizing it for decades. Strip the marketing off "making AI faster" and a large fraction of it is "making GEMM faster".

The anatomy of $O(n^3)$

Multiply an $m \times k$ matrix $A$ by a $k \times n$ matrix $B$ to get an $m \times n$ matrix $C$. One output cell is defined as

Read in words: "row $i$ of $A$ against column $j$ of $B$, multiplied elementwise and summed" — a single dot product. There are $mn$ cells, each costing $k$ multiplications and roughly $k$ additions. Counting a multiply and an add as one operation each,

$$\text{FLOPs} = 2\,m\,k\,n$$

which says nothing more exotic than: cells in the output ($mn$), times multiply-adds needed to fill one cell ($k$), times two operations per multiply-add.

That is $2n^3$ for square matrices. Double $n$ and the work grows eightfold. That is the content of $O(n^3)$, and it is also what people mean by a model's "compute".

Numbers give it scale. A square product at $n = 4096$ costs about $1.4 \times 10^{11}$ FLOPs — over a hundred billion operations for one multiplication. A Transformer stacks dozens of layers of those and repeats them per token.

It is worth feeling, with your hands, that the atom of all this is a dot product.

That covers the operation count. But runtime is not determined by $2n^3$ alone — the trucks from the opening are about to matter.

Arithmetic intensity: how much work per byte delivered

Here is the factory, written as a formula. Define the arithmetic intensity:

Read in words: the numerator is how much work you did, the denominator is how many bytes you paid to move in order to do it. Think of $I$ as fuel economy for data. It is the ratio of work done to data delivered — how many times you managed to reuse each byte you paid to fetch. An operation with low $I$ stalls on bandwidth no matter how fast the arithmetic units are (memory bound). With high $I$ you get to use the units at full tilt (compute bound). Achievable performance is capped by

$$\text{performance} \le \min\big(\text{peak compute},\ \ I \times \text{memory bandwidth}\big)$$

which says your real speed is pinned by whichever is lower: what the arithmetic units are capable of, or what the bandwidth can keep them fed with. The first term is the machine's ceiling, the second is "bytes per second times work per byte". Plotting those two lines gives you the roofline model.

Run the numbers for GEMM. For square $n \times n$ matrices you move three of them, $3n^2$ elements, which in fp32 is $12n^2$ bytes. The work is $2n^3$, so

$$I = \frac{2n^3}{12n^2} = \frac{n}{6}$$

Put in words, the arithmetic intensity of a matrix product is nothing but its side length divided by six. Intensity grows linearly with $n$: the bigger the matrix, the more times each fetched byte gets reused. That is the distinguishing property of matrix multiplication. At $n = 4096$, $I \approx 680$ — comfortably compute bound.

Contrast a matrix-vector product (GEMV): $2n^2$ FLOPs against roughly $4n^2$ bytes, so $I \approx 0.5$, fixed. Each element of the matrix is used exactly once, leaving nothing to reuse. This is why single-stream LLM token generation is slow and batching helps so dramatically: batching converts GEMV into GEMM and raises $I$.

What actually makes a GPU fast

Two things: sheer count of arithmetic units, and instructions cut to the shape of a matrix product.

A CPU is built to finish one instruction as quickly as possible, spending enormous silicon area on branch prediction and speculation. A GPU assumes the opposite workload — the same operation applied across a flood of data — so it strips control logic and packs the area with arithmetic units instead. Recent accelerators go further, with dedicated units that perform a small matrix multiply-accumulate as a single instruction. It is a deliberate bet on operations with high $I$.

Which is where intensity bites back. Multiply your arithmetic units by a thousand and any low-$I$ operation still stalls on bandwidth. A GPU is fast because matrix multiplication has high arithmetic intensity, not because computation in general got faster. That asymmetry drives real implementation decisions.

Tiling: cut the matrix into cache-sized blocks

$I = n/6$ assumed the ideal: $A$, $B$ and $C$ all resident in fast memory. In reality a $4096 \times 4096$ matrix does not fit in cache. A naive triple loop walks a column of $B$ on every pass of the inner loop, going out to main memory each time, and the effective intensity collapses to around 0.25.

The fix is tiling (blocking). Cut the matrices into $T \times T$ blocks so that one block each of $A$, $B$ and $C$ live in fast memory simultaneously. Once they are resident, that block performs $T^3$ operations on $3T^2$ elements of data. Across the whole matrix this gives

$$I \approx \frac{T}{6}$$

which says the number governing your fuel economy is no longer the side of the whole matrix, $n$, but the side of the block you managed to keep resident, $T$. Bigger $T$ is better, bounded by the capacity of cache — or shared memory on a GPU. What BLAS implementations have polished for decades is, at bottom, exactly this partitioning, applied at several levels down to registers. The algorithm stays $O(n^3)$ while data movement alone changes runtime by factors of several to an order of magnitude — the most blatant instance of the "outside the complexity class" effects from the complexity article.

There are also ways to attack the $O(n^3)$ itself. Strassen's algorithm uses divide and conquer to reach $O(n^{2.807})$ and does get used for very large matrices, at some cost in numerical stability. The theoretical exponent has been pushed below 2.4, but with constants so vast that those are "galactic algorithms" — provably better, never run.

Estimating FLOPs yourself

This is the destination. First, a tool.

def gemm_cost(m, k, n, bytes_per_elem=4):
    flops = 2 * m * k * n                          # one multiply-add = 2 FLOP
    moved = (m * k + k * n + m * n) * bytes_per_elem
    return flops, flops / moved                    # (FLOPs, arithmetic intensity)

print(gemm_cost(4096, 4096, 4096))   # → intensity 682.7 FLOP/byte
print(gemm_cost(4096, 4096, 1))      # → intensity   0.5 FLOP/byte

Same matrix $A$, and the intensity differs by more than a factor of a thousand depending on whether the other operand is a matrix or a vector. The value of batching falls straight out of the arithmetic.

For a whole model, estimate from the parameter count $N$. Each parameter contributes a multiply and an add, so a forward pass costs about $2N$ FLOPs per token. Training adds a backward pass that computes gradients with respect to both inputs and weights, costing roughly twice the forward, for a total of about $6N$ FLOPs per token (backpropagation). Training on $D$ tokens therefore costs

which says, in words: about six operations for every pairing of one parameter with one token. The bill for a training run is essentially model size $N$ multiplied by how much text $D$ you push through it. This approximation ignores attention's $s^2$ term and normalization, and holds while $d_{\text{model}} > s/12$, i.e. while the sequence length $s$ stays below about $12\,d_{\text{model}}$. At long context the attention term stops being negligible — remember that caveat.

Put numbers in. With $N = 7 \times 10^9$ and $D = 10^{12}$ tokens, $C \approx 4.2 \times 10^{22}$ FLOPs. On a setup sustaining $10^{15}$ FLOP/s — one PFLOP/s — that is $4.2 \times 10^{7}$ seconds, about 1.3 years. One line of arithmetic explains why training runs are spread across hundreds or thousands of accelerators.

Four things that matter in practice

1. Suspect the batch size When inference is slow, look at $I$ first. Batch-1 generation is GEMV and memory bound, so adding compute does nothing. Batching is the only lever that works.

2. Track effective FLOPs over peak FLOPs Back out effective performance from measured time and divide by peak. A low ratio tells you the units are idling, and that the problem is data movement rather than arithmetic.

3. Lowering precision is buying bandwidth Moving from fp32 to fp16 or int8 does not only speed up arithmetic — it halves or quarters the bytes moved, which raises $I$. That is why quantization pays off most in memory-bound regimes.

4. Do not write your own GEMM Tiling, prefetching and register allocation carry decades of accumulated tuning inside existing BLAS implementations. What you should be shaping is the form in which you call GEMM.

Summary

• AI compute reduces almost entirely to GEMM: $2mkn$ FLOPs for $m \times k$ by $k \times n$, or $2n^3$ when square
• Operation count alone does not set runtime. Arithmetic intensity $I$ = FLOPs / bytes moved decides memory bound versus compute bound
• GEMM has $I \approx n/6$ and improves with size; GEMV is stuck at $I \approx 0.5$. That gap is why batching works
• Tiling raises effective intensity toward $I \approx T/6$, bounded by cache. Runtime moves by an order of magnitude with the complexity unchanged
• Rules of thumb: $2N$ FLOPs per token for inference, $6ND$ for training — enough to size a run yourself

Next: how this computation gets split across many units — the limits of parallelism and the communication cost of distributed training.