Matrix multiplication and how tensor accelerators work

Strip away the names and an LLM is mostly one operation, repeated billions of times: multiply two matrices. Attention projections, MLPs, the output head are all matmuls. Every accelerator built for AI is, at heart, a machine for doing that one thing fast. Understanding the cost of a matmul explains most of what you see in inference and training performance.

What a matmul costs

Multiplying an m × k matrix by a k × n matrix produces m × n outputs. Each output is a dot product of length k, k multiplies and k adds:

FLOPs = 2 · m · n · k
bytes ≥ (m·k + k·n + m·n) · bytes_per_element      read A and B, write C

The ratio of the two is the arithmetic intensity, FLOPs per byte moved. It decides everything:

• A large square matmul, 4096 × 4096 × 4096 in FP16: 137 GFLOP over about 100 MB of data, about 1,365 FLOPs per byte. On an RTX 4090 that’s about 0.83 ms of compute against 0.1 ms of memory traffic. It’s compute-bound.
• The same matrix times one vector, which is decode at batch size 1: 33.5 MFLOP over 33.5 MB of weights, 1 FLOP per byte. That’s about 33 µs of memory traffic against 0.2 µs of compute. It’s memory-bound by a factor of about 160.

Same weights, same hardware, two completely different bottlenecks. That’s the whole story of prefill versus decode (the decode note has the consequences).

The roofline

The roofline model (Williams, Waterman and Patterson, 2009) turns that into one picture. Attainable throughput is the lower of two ceilings:

attainable FLOP/s = min( peak FLOP/s ,  arithmetic intensity × memory bandwidth )

The ridge point, where the two ceilings meet, is peak / bandwidth. Below it you’re bandwidth-bound; above it, compute-bound.

GPU	Dense 16-bit tensor peak	Memory bandwidth	Ridge point
RTX 4090	~165 TFLOPS (FP32 accumulate)	1,008 GB/s	~164 FLOPs per byte
H100 SXM	~989 TFLOPS (BF16)	3.35 TB/s	~295 FLOPs per byte

Single-stream decode sits around 1–4 FLOPs per byte, two orders of magnitude below either ridge. That’s why batching raises throughput almost for free: B sequences multiply the intensity by B (the batching note has the numbers). And it’s why newer GPUs, whose compute grows faster than their bandwidth, push the ridge further out of reach of decode.

Tiling: reuse what you’ve loaded

A naive kernel reads every element of A and B from memory once per output it contributes to, which is far too much traffic. Real kernels split the output into tiles. Each tile’s slices of A and B are loaded once into fast on-chip memory (shared memory and registers) and reused for every multiply in the tile.

For a T × T output tile, the kernel loads about 2·T·k elements and performs 2·T²·k FLOPs, so the intensity grows with T. Bigger tiles mean more reuse, up to the limit of on-chip memory. Most of the art of a fast GEMM kernel is choosing tile shapes for a specific chip, and overlapping loading the next tile with computing the current one.

Tensor cores: matrix math as one instruction

A conventional GPU core does one fused multiply-add per lane per cycle. NVIDIA’s tensor cores, introduced with Volta in 2017, do a small matrix product per instruction instead. On Volta, each tensor core computed D = A × B + C on 4 × 4 tiles per clock, with FP16 inputs and FP32 accumulation. Later generations widened the tiles and added formats: TF32 and BF16 on Ampere, FP8 on Hopper, FP4 on Blackwell. Hopper also added warpgroup-level instructions and a hardware engine (TMA) for moving tiles into shared memory asynchronously.

The point isn’t the specific shapes. Dedicating silicon to the dense inner loop of matmul delivers far more FLOPs per watt and per square millimetre than general-purpose lanes. On every modern AI GPU, most of the advertised FLOPs are tensor-core FLOPs, and a kernel that doesn’t use them leaves most of the chip idle.

Systolic arrays: the TPU approach

Google’s first TPU (Jouppi et al., 2017) took the idea further with a systolic array: a 256 × 256 grid of 65,536 8-bit multiply-accumulate units. Weights are held stationary in the cells, activations flow in from one side, and partial sums flow through the grid, moving one cell per clock “like blood pumped by a heart”, which is where the name comes from. Data is reused as it moves between neighbouring cells instead of going back to memory. Later TPUs, and many inference chips, use the same principle at different scales.

Precision keeps shrinking

Every halving of element size halves the bytes moved, doubles the elements per cache line, and roughly doubles the tensor throughput on hardware that supports it. The history of AI number formats follows:

• FP32 → FP16: half the bytes, but FP16’s narrow exponent range overflows easily in training.
• BF16: FP32’s 8-bit exponent with a short mantissa. Same range as FP32 at half the size, now the training default.
• FP8 (E4M3 / E5M2): used for training at scale. DeepSeek-V3 trained mostly in FP8 (the DeepSeek note covers it).
• 4-bit weights for inference: usually stored in a block format and dequantised to 16-bit in registers just before the tensor-core instruction. That’s exactly the fused-dequant pattern in the quantisation and attention notes.

Rules I use

1. Compute the arithmetic intensity first. It tells you whether to optimise compute or memory before you write a line of kernel code.
2. Know your GPU’s ridge point. Below it, fewer bytes per FLOP (quantisation, batching, fusion) is the only lever that works.
3. Use the tensor cores. If a profile shows the matmuls on the regular FP32 pipes, fix the data type or the shapes first.

Sources: Williams, Waterman and Patterson, Roofline (CACM 2009). NVIDIA’s Volta, Ampere and Hopper architecture whitepapers. Jouppi et al., In-Datacenter Performance Analysis of a Tensor Processing Unit (ISCA 2017). Peak figures are vendor specifications; intensities and times are computed from the formulas above.