A modern CPU core has a secret second gear. Alongside the scalar registers that hold one number at a time sit wide vector registers that hold eight, sixteen, or more numbers, and a single machine instruction can add, multiply, or fuse-multiply-add across all of them at once. This is SIMD — Single Instruction, Multiple Data — and it is the difference between a CPU that limps through a transformer at a token per second and one that streams a small language model at conversational speed. Almost everything a transformer does at inference time is arithmetic over long, regular arrays: matrix multiplications, a softmax over a row of scores, a layernorm over a hidden vector. Those are exactly the shapes SIMD was built for. This piece works through the whole picture from first principles: what a lane is, what AVX2, AVX-512, and ARM NEON give you, how the matmul inner loop, softmax, and layernorm actually vectorize, the fused multiply-add and the awkward horizontal reduction, the theoretical peak throughput and the memory-bandwidth wall that keeps you from reaching it, the low-precision instructions (BF16, VNNI, AMX) that push the ceiling higher, the alignment and tail-handling details that separate a correct kernel from a fast one, and why all of it is the bedrock of CPU-based small-model inference.
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.
Alignment: feeding the vector units cleanly
SIMD loads and stores move a whole register’s worth of bytes at once, and the address matters. An aligned load — where the address is a multiple of the register width (32 bytes for AVX2, 64 for AVX-512) — is the fast path; the data sits neatly within cache lines and the load completes in the fewest cycles. A misaligned vector load can straddle two cache lines and cost extra, and on older ISAs the wrong instruction on a misaligned address faults outright. So high-performance kernels align their buffers deliberately — allocating tensors on 64-byte boundaries with aligned allocators — and the packing step in a GEMM does double duty: it not only reorders B into a vector-friendly layout but lands it on aligned addresses.
Alignment interacts with the whole memory hierarchy. Contiguous, aligned, predictably strided access lets the hardware prefetcher run ahead and hide latency, and it keeps each cache line fully used rather than half-wasted across a boundary. This is a big reason kernels pack and tile rather than reading weights in their natural layout: the goal is not just to vectorize the arithmetic but to present memory to the SIMD units as long, aligned, sequential streams so the load ports and prefetchers never stall the FMA pipelines. On modern x86, unaligned-load penalties have shrunk, but aligned, contiguous, cache-line-friendly access is still the difference between feeding the vector units and starving them.
Tail handling: when the array does not divide evenly
Vector loops chew through data in fixed strides — 16 fp32 elements at a time on AVX-512 — but real tensors have dimensions that are not always tidy multiples of the lane count. A hidden size of 4096 divides evenly by 16; a vocabulary of 32000 or an odd sequence length may not. Whatever is left after the last full vector is the tail (or remainder), and handling it correctly is where naive kernels either read out of bounds or silently corrupt neighboring memory.
The classic fix is a scalar remainder loop: run the wide vector loop over the largest multiple of the lane count, then finish the last few elements one at a time. It is correct but adds a branch and leaves the vector units idle for the tail. AVX-512 offers a cleaner answer with its mask registers: a masked load/store activates only the valid lanes on the final iteration — compute a mask like (1 << remainder) - 1, and the hardware processes only those lanes while leaving the rest untouched, no separate loop and no out-of-bounds access. ARM SVE bakes this in as the core idiom: every loop is predicated, so the last partial iteration is handled by the same code with a narrower predicate. Tail handling sounds like a footnote, but a wrong remainder is the most common source of both crashes and subtle numerical bugs in hand-written SIMD kernels — and an unvectorized tail on a short row can quietly dominate the runtime.
Why this is the foundation of CPU small-model inference
Put it together and SIMD is the reason a small language model can run usably on a laptop or a commodity server with no GPU. A GPU wins on raw parallelism, but a modern CPU with AVX-512 or NEON, fed a well-tiled, well-aligned, low-precision kernel, delivers hundreds of GFLOP/s per core and single-digit-to-teens tokens per second on a 1–8B parameter model — enough for local assistants, on-device privacy, edge deployment, and batch back-ends where a GPU is overkill or unavailable. Projects like llama.cpp, oneDNN, and ONNX Runtime are, at their performance core, libraries of exactly the SIMD kernels described here: vectorized GEMM, fused vectorized softmax and RMSNorm, and quantized int8/int4 inner loops.
The mental model to carry away is a chain of ceilings. SIMD width sets the compute peak; FMA and multiple ports raise it; low precision (bf16, int8, VNNI, AMX) raises it further and relieves the memory wall by shrinking the bytes moved; but arithmetic intensity and memory bandwidth decide how much of that peak you can actually reach, which is why quantization so often beats micro-optimizing the math. Alignment and tail handling are the correctness-and-efficiency details that keep the vector units fed and the results right. Master those layers — lanes, FMA, reductions, roofline, precision, alignment — and you understand why the same model runs at a crawl in a naive loop and at conversational speed in a tuned CPU runtime.