What SIMD actually is: one instruction, many lanes
Scalar code processes one value per instruction: load a, load b, add them, store the result, repeat. SIMD widens that pipe. A vector register is a fixed-width bucket — 128, 256, or 512 bits — partitioned into equal lanes. Pack eight 32-bit floats into a 256-bit register and a single vaddps instruction adds eight pairs simultaneously, producing eight sums in the time one scalar add would take. The same instruction is applied to every lane in lock-step; the lanes do not interact.
The number of lanes is just the register width divided by the element size. A 512-bit register holds 512 / 32 = 16 fp32 values, or 512 / 16 = 32 bf16 values, or 512 / 8 = 64 int8 values — narrower data means more lanes and more throughput per instruction. This is data parallelism: the same operation over many independent elements. It is not threads and not cores — it is parallelism inside a single core, inside a single instruction. A CPU with 8 cores, each running 16-wide SIMD, has 8 × 16 = 128 fp32 operations in flight per clock before you even count multiple execution ports. For transformer math, whose inner loops are long runs of identical arithmetic over tensors, that multiplier is the whole game.
The instruction sets: AVX2, AVX-512, and ARM NEON
On x86, the relevant families are AVX2 and AVX-512. AVX2 uses 256-bit YMM registers — 8 fp32 lanes — and is nearly universal on desktop and server CPUs from the last decade. AVX-512 doubles the width to 512-bit ZMM registers — 16 fp32 lanes — and, just as importantly, adds mask registers (k0–k7) that let you predicate individual lanes on and off. It also expands the register file from 16 to 32 vector registers, which matters a great deal for the register tiling that fast matmul kernels depend on.
On ARM — Apple Silicon, AWS Graviton, most phones — the baseline is NEON, 128-bit registers with 4 fp32 lanes. That sounds narrow next to AVX-512, but ARM cores typically issue several NEON instructions per cycle from multiple pipelines, so real throughput is closer than the lane count alone suggests. Newer ARM adds SVE/SVE2 (Scalable Vector Extension), whose defining trick is a vector-length-agnostic programming model: you write the loop once with predication and it runs correctly whether the hardware vector is 128 or 2048 bits wide. The practical upshot is that portable transformer kernels are written against whichever ISA is present, usually behind intrinsics or a library that dispatches at runtime to the widest set the CPU actually supports.
FMA: the fused multiply-add at the heart of it all
Transformer math is dominated by one primitive: multiply two numbers and add the product to a running total. A dot product is nothing but sum += a_i * b_i repeated. The fused multiply-add (FMA) instruction does exactly this in a single operation across every lane: d = a * b + c, computed with one rounding step instead of two. That single-rounding is not just faster — it is more accurate, because the intermediate product is kept at full internal precision before the add.
Crucially, FMA counts as two floating-point operations (a multiply and an add) but issues as one instruction, so it doubles your effective FLOP rate. A 16-lane AVX-512 FMA performs 16 × 2 = 32 floating-point operations per instruction. Better still, most server cores have two FMA execution units, so with enough independent accumulators to hide the instruction latency you can retire two FMAs per clock: 16 lanes × 2 ops × 2 units = 64 FLOPs/cycle per core. This is why well-written matmul kernels keep several accumulator registers live at once — not for the parallelism of the data, which SIMD already exploits, but to keep both FMA pipelines fed despite the multi-cycle latency of each individual FMA. FMA throughput, not add or multiply throughput, is the number that sets the compute ceiling for a transformer on CPU.
Vectorizing the matmul inner loop
Every heavy transformer operation — the QKV projections, the attention output projection, the two big feed-forward matrices — is a matrix multiply, and matmul is where SIMD earns its keep. Consider C = A × B where A is [M, K] and B is [K, N]. The naive triple loop computes each output C[i,j] = Σ_k A[i,k] * B[k,j] as a scalar dot product. The vectorized version instead computes a whole strip of the output row at once: hold C[i, j:j+16] in a ZMM accumulator, and for each k broadcast the single scalar A[i,k] across all 16 lanes and FMA it against the contiguous vector B[k, j:j+16].
Written as the outer-product or broadcast formulation, the inner loop becomes a clean sequence of acc = fma(broadcast(A[i,k]), B[k, j:j+16], acc) with no horizontal reduction at all — each lane accumulates its own output column independently. Real kernels go further and register-tile: they compute, say, a 6 × 32 block of C using a dozen accumulator registers, reusing each loaded value of A and B across the whole block so that arithmetic dominates memory traffic. Layered on top is cache blocking — tiling the K and N loops so the working set fits in L1/L2 — and often a packing step that copies B into a contiguous, vector-friendly layout. This is the design at the core of GEMM libraries like oneDNN, OpenBLAS, and llama.cpp’s hand-tuned kernels.
Horizontal reductions: the awkward part
Lanes are independent by design, which is wonderful until you need to combine them. A dot product done the ‘obvious’ way — multiply two vectors elementwise, then sum the products — ends with 16 partial results sitting in 16 lanes that must be collapsed into one scalar. That collapse is a horizontal reduction, and it is comparatively expensive: SIMD hardware is optimized for vertical (lane-parallel) work, while summing across lanes requires shuffle and permute instructions that shift data between lanes, typically a log2(lanes) tree of shuffle-and-add steps.
The professional’s answer is to avoid horizontal reductions in the hot loop entirely. That is precisely why the matmul kernel above uses the broadcast formulation: each lane owns a distinct output column, so the accumulation is purely vertical and the reduction never happens inside the K loop. Where a reduction is genuinely unavoidable — the sum in softmax, the mean and variance in layernorm, an L2 norm — you keep several independent vector accumulators running through the loop to hide FMA latency, and perform a single horizontal reduction at the very end, once, over the combined accumulators. Amortized over thousands of loop iterations, that one final shuffle tree costs essentially nothing. The rule of thumb: reduce vertically for as long as you can, and pay the horizontal tax exactly once.
Vectorizing softmax
Attention needs a softmax over each row of scores: softmax(x_i) = exp(x_i) / Σ_j exp(x_j). For numerical safety it is computed in the shifted form exp(x_i - max_j x_j) so the exponentials never overflow. That gives three passes over the row, and each vectorizes cleanly. Pass one finds the maximum: load the row 16 lanes at a time and keep a running lane-wise vmax, then do a single horizontal max at the end. Pass two computes exp(x_i - max) and a running lane-wise sum, closing with one horizontal add. Pass three multiplies each element by the reciprocal of that sum.
The interesting part is exp itself, which has no hardware instruction. Fast kernels use a vectorized polynomial approximation: reduce the argument by splitting x = n·ln2 + r, compute 2^n by directly assembling the float exponent bits, and approximate exp(r) for the small remainder with a low-degree polynomial evaluated in FMA-friendly Horner form — all in vector registers, all 16 lanes at once, no branches. A common optimization fuses passes two and three, or even folds the whole softmax into the attention kernel (the online-softmax trick that FlashAttention popularized) so the score row is never fully materialized. The reductions are real but, done once at the end of each pass, they are a rounding error against the exponential work.
Vectorizing layernorm and RMSNorm
Normalization runs on every token at every layer, so its kernel matters. LayerNorm over a hidden vector x of dimension d computes mean = (1/d) Σ_i x_i, then var = (1/d) Σ_i (x_i - mean)^2, then normalizes and applies a learned scale and shift: y_i = γ_i (x_i - mean) / √(var + ε) + β_i. The two sums are horizontal reductions, but each is fed by a vertical pass: accumulate Σ x_i and Σ x_i^2 lane-wise across the vector (a single pass suffices, since var = E[x^2] - E[x]^2), then reduce once. The final normalize-and-affine step is pure vertical arithmetic: broadcast the scalar mean and inverse-standard-deviation across all lanes and FMA against the γ and β vectors.
RMSNorm, used by Llama-style models, is even friendlier: it drops the mean-centering entirely and normalizes by the root-mean-square, y_i = γ_i · x_i / √((1/d) Σ_j x_j^2 + ε). That is a single sum-of-squares reduction and one broadcast-multiply — fewer passes, one reduction, and no subtraction of the mean. The inverse square root is handled with a fast reciprocal-sqrt approximation refined by a Newton step, again vectorized across lanes. Because d (say 2048 or 4096) is a tidy multiple of the lane count, these kernels are almost all vertical FMA with a single cheap horizontal reduction at the tail — close to the ideal SIMD shape.
A worked lane and throughput example
Make the numbers concrete. Take one CPU core at 3.0 GHz with AVX-512 and two FMA units, running fp32. Per cycle it can retire:
lanes per ZMM (fp32) = 512 / 32 = 16
FLOPs per FMA per lane = 2 (one multiply + one add)
FMA units (ports) = 2
-------------------------------------------------
FLOPs per cycle = 16 x 2 x 2 = 64
peak per core = 64 x 3.0e9 = 192 GFLOP/s (fp32)
8-core peak = 192 x 8 = 1.54 TFLOP/s (fp32)Now compare against scalar code, which does one FMA lane per instruction: 1 × 2 × 2 × 3.0e9 = 12 GFLOP/s. The SIMD version is 16× faster per core — exactly the lane count, as expected, since the win comes entirely from processing 16 elements per instruction. Switch the data to bf16 and the register holds 512 / 16 = 32 lanes, doubling the ceiling again; drop to int8 with a VNNI dot-product instruction and it doubles once more. So a single core spans roughly 12 GFLOP/s scalar to ~200 GFLOP/s fp32 SIMD to well over 700 GOP/s int8 — a 60-plus× range decided almost entirely by how wide your instructions are and how narrow your data is. That spread is why choosing the right ISA and precision is the single biggest lever in CPU transformer performance.
The catch: memory bandwidth caps the speedup
Those peak FLOP numbers are a ceiling you will rarely touch, because a transformer at inference is usually memory-bandwidth-bound, not compute-bound. The deciding metric is arithmetic intensity: FLOPs performed per byte read from memory. During single-stream decode (batch size 1), each output token is a matrix-vector product against the weights, and every weight is loaded from RAM, used for exactly one multiply-add (2 FLOPs), and discarded. In fp32 that is 2 FLOPs / 4 bytes = 0.5 FLOP/byte — catastrophically low.
The roofline model makes the consequence exact. Sustainable throughput is min(peak_FLOPs, arithmetic_intensity × bandwidth). With, say, 50 GB/s of usable memory bandwidth and intensity 0.5, you are capped at 0.5 × 50e9 = 25 GFLOP/s — less than one seventh of that 192 GFLOP/s SIMD peak. The FMA units sit idle waiting for weights to arrive. This is why decode speed on CPU tracks memory bandwidth almost linearly, why quantization is the highest-leverage optimization (halving the weight bytes with int8/int4 roughly doubles decode throughput by doubling arithmetic intensity), and why prefill — a matrix-matrix product with high weight reuse and therefore high arithmetic intensity — can be compute-bound and actually approach the SIMD peak. SIMD raises the roof; bandwidth decides how much of it you can use.
Low precision on x86: AVX-512 BF16 and VNNI
If bandwidth is the wall, narrower numbers are the ladder over it, and recent x86 ships dedicated low-precision instructions. AVX-512 BF16 adds vdpbf16ps, which takes bf16 inputs and accumulates dot products into fp32 — a bf16 vector holds twice the lanes of fp32, so you halve the memory footprint of weights and activations while keeping fp32 accumulation precision, which is exactly the numerics transformers tolerate well. Bfloat16’s trick is that it keeps the full 8-bit exponent of fp32 and simply truncates the mantissa, so the dynamic range is identical and conversion is nearly free.
AVX-512 VNNI (Vector Neural Network Instructions) targets integer inference. vpdpbusd multiplies 8-bit integers and accumulates into 32-bit integers in a single instruction — what used to take a multiply, a widening, and an add now issues as one op, roughly 3× fewer instructions for an int8 dot product, on top of the 4× lane density of int8 over fp32. For a quantized small model this is a double win: fewer bytes moved (relieving the bandwidth bottleneck) and more work per instruction. The cost is quantization’s usual bookkeeping — per-channel scales, zero-points, and careful calibration to hold accuracy — but for CPU-served models the throughput gain is large enough that int8 VNNI is the default fast path in libraries like oneDNN and the quantized kernels in llama.cpp.
Beyond SIMD: Intel AMX and matrix tiles
SIMD applies one instruction to a one-dimensional vector. Intel’s AMX (Advanced Matrix Extensions), introduced on Sapphire Rapids Xeons, goes a dimension further: it exposes eight two-dimensional tile registers (up to 16 × 64 bytes each) and a TDP instruction that performs an entire tile matrix multiply — a block of C += A × B — in one operation. Where an FMA does a vector’s worth of multiply-adds, one AMX instruction does a small matrix’s worth, dramatically raising the FLOPs delivered per instruction issued and per byte fetched.
AMX is built for exactly the low-precision GEMM that transformer inference runs on: it operates on int8 and bf16 tiles with int32/fp32 accumulation, and a single core can deliver on the order of a teraflop of bf16 throughput — several times what AVX-512 FMA reaches. For prefill and for batched serving, where arithmetic intensity is high enough to feed it, AMX turns a CPU into a credible matrix engine. It does not repeal the bandwidth wall for batch-1 decode — a matrix engine still starves without operand reuse — but it is the clearest sign that CPUs are being redesigned around the shape of transformer math, and libraries increasingly dispatch to AMX tiles for the GEMM-heavy phases when the hardware is present.