LUNAROPS · OPERATIONAL UPLINK 100% UPTIME 1,247d POSTS 893 JEFF.MOON@LUNAROPS.DEV UTC --:--:--

Why GPUs Beat CPUs at Matrix Math

20 min read The Physics of Computing
gpucudasimtparallel-computingtensor-coreshardware
Contents

A single GPU core is not impressive. It runs at a lower clock than a desktop CPU core, it has a feeble branch predictor, it cannot reorder instructions to any meaningful degree, and on its own it would lose a footrace to a Pentium from two decades ago. The reason a modern GPU obliterates a CPU at matrix multiplication has nothing to do with any individual core being good. It is that the GPU is a throughput machine built around a single bet: that you have so much independent arithmetic to do that the right strategy is to launch tens of thousands of threads, run them in lockstep bundles, and use that enormous pool of pending work to paper over the fact that memory is slow. The CPU makes the opposite bet. It assumes your work is serial, branchy, and latency-sensitive, so it spends its transistor budget on caches, out-of-order execution, and prediction to make one stream of instructions finish as fast as physically possible. Matrix math is the canonical workload where the GPU’s bet pays off and the CPU’s does not, and understanding exactly why requires looking at how the GPU schedules, where it keeps data, and what it added in silicon specifically to make general matrix multiply (GEMM) cheaper.

The arithmetic is embarrassingly parallel, and that changes everything

Start with the workload itself, because the hardware is shaped to fit it. A matrix multiply C = A x B for an M x K times K x N problem computes M x N output elements, and each output element is an independent dot product of length K. There are no dependencies between output elements. Element C[0][0] does not need to know anything about C[5][7]. That property, full independence across the output, is what the literature politely calls “embarrassingly parallel,” and it is the single most important fact about why GPUs win here.

The total work is 2 x M x N x K floating-point operations (one multiply and one add per term, K terms per element). For a 4096-cubed multiply that is roughly 137 GFLOP of arithmetic. The CPU’s entire architectural philosophy, finishing one instruction stream quickly, is irrelevant when you have 16 million independent dot products to compute. What you want is not a fast core. You want many ALUs all doing fused multiply-adds at once, and you want to keep them fed. The GPU is the machine that took this observation to its logical extreme.

A fused multiply-add, or FMA, is the atom of this entire discussion. It computes d = a x b + c in a single instruction with a single rounding step, and it is the workhorse of dot products. Both CPUs and GPUs have FMA units. The difference is how many, and how they are scheduled.

SIMD versus SIMT: the distinction people get wrong

Here is where careful language matters, because “GPUs are just wide SIMD” is a half-truth that obscures the actual design.

A CPU goes wide using SIMD: Single Instruction, Multiple Data. An AVX-512 instruction operates on a 512-bit vector register, which is sixteen FP32 lanes, and a single instruction multiplies all sixteen in one shot. The programmer (or compiler) is responsible for packing data into those vector registers and for handling the leftover elements when the array length is not a multiple of sixteen. The vector width is part of the instruction set and is exposed directly in the assembly. There is one program counter for those sixteen lanes, and they are explicitly a vector.

A GPU goes wide using SIMT: Single Instruction, Multiple Threads, NVIDIA’s term for its execution model. You write a scalar program, a kernel, as if for one thread. You then launch a grid of thousands or millions of these threads. The hardware groups threads into bundles of 32 called a warp, and the 32 threads of a warp share one program counter and one instruction-fetch/decode unit, executing the same instruction in lockstep on their own private data. From the silicon’s point of view a warp behaves a lot like a 32-lane SIMD unit. From the programmer’s point of view, each thread looks scalar and independent, with its own registers and its own index.

That difference is not cosmetic. It changes how divergence is handled, and divergence is the crux.

Property SIMD (CPU, e.g. AVX-512) SIMT (GPU, CUDA warps)
Unit of execution One vector register, fixed lanes One warp = 32 threads
Programming model Explicit vectors, lane-packing Scalar per-thread kernel
Width visible in ISA? Yes (128/256/512-bit) No (warp size implicit)
Branch divergence Programmer uses masks/blends Hardware masks lanes, executes both sides serially
Per-lane addressing Awkward (gather/scatter) Natural (each thread computes its own address)
Threads per scheduler 1 (the vector lane group) Many warps oversubscribed, swapped each cycle

When a CPU SIMD program hits a branch where some lanes go one way and others go the other, the programmer or compiler has to compute both sides and blend the results with a mask. When a GPU warp hits a data-dependent branch where some of its 32 threads take the if and others take the else, the hardware does something similar but automatically: it executes the if path with the else threads masked off (idle), then executes the else path with the if threads masked off. This is warp divergence, and it is the single biggest performance trap on a GPU. Both paths run; you pay for both. If all 32 threads agree on the branch, the warp is “convergent” and there is no penalty. Matrix multiply is wonderful precisely because the inner loop is branch-free: every thread does the same FMAs over the same loop bounds, so warps never diverge.

The honest framing is this. SIMT is a programming-model convenience layered on top of hardware that is, underneath, executing 32-wide lockstep vectors. It buys you per-thread addressing and per-thread control flow at the cost of running both sides of any divergent branch. For regular, dense, branch-free work like GEMM, you get the ergonomics of scalar code and the throughput of wide vectors, which is the best of both worlds. For branchy irregular work you get the throughput of wide vectors and the misery of running every path, which is the worst.

The execution hierarchy: grids, blocks, warps, and the SM

To see why oversubscription works, you have to look at the actual scheduling hierarchy. A CUDA kernel launch creates a grid of thread blocks. Each block is a group of up to 1024 threads that are guaranteed to run on the same Streaming Multiprocessor (SM) and can communicate through fast on-chip shared memory and synchronize with barriers. The block is the unit of resource allocation and cooperation. Within a block the hardware chops threads into warps of 32. The SM is the actual physical processor: a modern data-center GPU has on the order of 100 to 150 SMs, and each SM contains multiple warp schedulers, a large register file, shared memory/L1, and the FP32, FP64, INT, and tensor execution units.

                       GPU
   +-----------------------------------------------+
   |  SM 0      SM 1      SM 2     ...    SM N       |
   |  +------------------------------------------+   |
   |  |              Streaming Multiprocessor     |  |
   |  |  +-----------+  +-----------+              |  |
   |  |  | Warp Sched|  | Warp Sched|  (x4)        |  |
   |  |  +-----------+  +-----------+              |  |
   |  |   |  |  |  |     |  |  |  |                 |  |
   |  |  [warp][warp][warp][warp]...  <- many      |  |
   |  |   resident warps, swapped each cycle       |  |
   |  |                                            |  |
   |  |  Register File (e.g. 64K x 32-bit regs)    |  |
   |  |  Shared Memory / L1  (~100-228 KB)         |  |
   |  |  FP32 ALUs   Tensor Cores   LD/ST units    |  |
   |  +------------------------------------------+   |
   +-----------------------------------------------+
              |                         |
          L2 Cache (tens of MB)         |
              |                         |
          HBM / Global Memory (~1-8 TB/s, ~400+ ns latency)

The warp scheduler is the heart of the design. Each cycle, a scheduler looks at the pool of warps resident on its SM and picks one that is ready to issue an instruction, meaning that warp’s next instruction has all its operands available. It issues that warp’s instruction to the execution units. Crucially, picking a ready warp is nearly free: warp context (registers, program counter) lives in hardware the whole time the warp is resident, so a context switch between warps costs essentially nothing, unlike a CPU thread switch that saves and restores state to memory. This is zero-overhead warp switching, and it is the mechanism that hides latency.

Latency hiding: the entire trick in one idea

A CPU hides memory latency with caches and out-of-order execution. When a load misses to DRAM, costing hundreds of cycles, the CPU’s out-of-order engine tries to find independent later instructions to execute in the shadow of the miss, and a big speculation window keeps a few dozen instructions in flight. It works, but it is expensive in transistors and limited in depth.

A GPU hides memory latency with threads. It does not try to make any single warp fast. When a warp issues a global memory load and then needs the result, that warp stalls, possibly for 400-plus nanoseconds waiting on HBM. The scheduler simply ignores that stalled warp and issues from a different ready warp instead. If there are enough warps resident, there is always at least one ready to run while others wait on memory, and the execution units never go idle. The latency of any individual load is never reduced; it is merely overlapped with useful work from other warps. This is the whole game. The GPU trades single-thread latency, which it does not care about, for aggregate throughput, which it cares about exclusively.

This is why GPUs deliberately oversubscribe. To fully hide a 400 ns memory latency on an SM that can issue many FMAs per nanosecond, you need a deep pool of independent warps in flight. The ratio of resident warps to the maximum the SM can hold is called occupancy, and it is the central tuning knob of GPU performance. Low occupancy means too few warps to cover the stalls, the schedulers run dry, and the ALUs sit idle waiting on memory. The roofline mental model captures the consequence: a kernel is either compute-bound (limited by FLOP/s) or memory-bound (limited by bytes/s from HBM), and the GPU’s job is to stay on the compute roof, which requires enough parallel work to hide the memory traffic. GEMM, with its high arithmetic intensity, is the poster child for staying compute-bound.

Why register pressure and occupancy dominate

Here is the trade-off that surprises people coming from CPU-land. The thing that limits how many warps you can keep resident is not usually the number of execution units. It is on-chip storage, specifically the register file and shared memory, which are partitioned among resident threads.

Each SM has a fixed register file, for example 65,536 32-bit registers. Every thread that is resident consumes some of those registers for its private variables, and that allocation is fixed for the lifetime of the block. If your kernel compiles to use 32 registers per thread, an SM can host 65536 / 32 = 2048 threads, which is 64 warps, plenty for high occupancy. If your kernel is register-hungry and uses 128 registers per thread, the SM can only host 65536 / 128 = 512 threads, 16 warps, and your occupancy drops to a quarter. Fewer warps means less latency hiding, which can starve the ALUs. This is register pressure: the compiler spilling complex kernels into more registers directly shrinks the number of warps you can run, even though no execution unit is the bottleneck.

The same applies to shared memory. A block that allocates a lot of shared memory for its tile of the matrix limits how many blocks can co-reside on an SM, because shared memory is partitioned too. The art of writing a fast GEMM kernel is balancing tile size (bigger tiles reuse more data and reduce memory traffic) against the register and shared-memory pressure those tiles create (bigger tiles can crater occupancy). There is no free lunch; you tune for the sweet spot where you have enough occupancy to hide memory and enough on-chip reuse to stay compute-bound.

Resource per SM (representative) Approx capacity What it limits
Registers 65,536 x 32-bit Threads/warps resident (register pressure)
Shared memory / L1 ~100-228 KB Blocks resident, tile size
Max resident warps ~48-64 Ceiling on latency hiding
Warp schedulers ~4 Instructions issued per cycle

The non-intuitive lesson: making a GPU kernel “smarter” by adding more per-thread state often makes it slower, because it sacrifices the parallelism that the whole architecture depends on. CPUs reward sophisticated single-thread code; GPUs punish it.

The memory hierarchy and coalesced access

Throughput is meaningless if you cannot feed the ALUs. The GPU memory hierarchy is explicitly tiered, and unlike a CPU’s mostly transparent caches, some of these tiers are programmer-managed.

   Registers      <- per-thread, ~1 cycle, fastest, scarce
       |
   Shared Memory  <- per-block, on-chip, ~20-30 cycle, programmer-managed
       |
   L1 / L2 Cache  <- hardware-managed
       |
   Global (HBM)   <- per-GPU, huge, ~400+ cycle latency, ~1-8 TB/s
Memory space Scope Latency (rough) Bandwidth Managed by
Registers Single thread ~1 cycle Enormous Compiler
Shared memory Thread block ~20-30 cycles Very high Programmer
L2 cache Whole GPU ~200 cycles High Hardware
Global / HBM Whole GPU ~400-600 cycles ~1-8 TB/s Hardware

The performance difference between a naive and a fast GEMM is almost entirely a memory-hierarchy story, not a math story. A naive kernel where each thread reads its inputs straight from global memory recomputes the same loads over and over and is hopelessly memory-bound. A fast kernel uses shared memory as a software-managed scratchpad: a block cooperatively loads a tile of A and a tile of B from global memory into shared memory once, synchronizes, and then every thread in the block reuses those on-chip tiles many times for its FMAs. This is “blocking” or “tiling,” the same idea as cache blocking on a CPU, except here you control the fast memory by hand. Reuse is what converts a memory-bound problem into a compute-bound one and lets you climb to the compute roof.

The other half is coalescing. Because a warp issues one memory instruction for all 32 threads, the addresses those 32 threads request matter enormously. If thread i accesses element base + i, the 32 accesses fall in a contiguous, aligned block and the hardware services them in one (or a few) memory transactions: this is a coalesced access, and it achieves full bandwidth. If the 32 threads access scattered addresses, say with a large stride or a random gather, the hardware must issue many separate transactions, each fetching a wide line of which only a few bytes are used, and effective bandwidth collapses, sometimes by an order of magnitude. Writing fast GPU code is in large part the discipline of arranging data so that consecutive threads touch consecutive addresses. For matrix math this dictates whether you store and traverse matrices row-major or column-major and how you index tiles, and it is exactly the kind of structured access that dense linear algebra makes easy and that irregular pointer-chasing makes impossible.

Here is the shape of a tiled multiply, simplified, to make the reuse and coalescing concrete:

#define TILE 16
__global__ void sgemm(const float* A, const float* B, float* C,
                       int M, int N, int K) {
    __shared__ float As[TILE][TILE];
    __shared__ float Bs[TILE][TILE];

    int row = blockIdx.y * TILE + threadIdx.y;
    int col = blockIdx.x * TILE + threadIdx.x;
    float acc = 0.0f;

    for (int t = 0; t < K; t += TILE) {
        // Cooperative, coalesced load of one tile into shared memory.
        As[threadIdx.y][threadIdx.x] = A[row * K + (t + threadIdx.x)];
        Bs[threadIdx.y][threadIdx.x] = B[(t + threadIdx.y) * N + col];
        __syncthreads();

        // Reuse the on-chip tile: TILE FMAs per loaded element.
        for (int k = 0; k < TILE; ++k)
            acc = fmaf(As[threadIdx.y][k], Bs[k][threadIdx.x], acc);
        __syncthreads();
    }
    if (row < M && col < N) C[row * N + col] = acc;
}

Every thread loads two values per tile but performs TILE FMAs against shared memory, so the global-memory traffic per FMA drops by a factor of TILE. The __syncthreads() barriers ensure the tile is fully loaded before anyone reads it. Production kernels go far beyond this with register-level blocking, double buffering, and vectorized loads, but the principle is identical: load once into fast memory, reuse many times, keep accesses coalesced. If you want the broader systems picture of how these chips get deployed and fed at scale, see GPU infrastructure for ML, and for the bandwidth side of the wall, HBM4 and the memory wall for AI.

Tensor cores: what consumer FP32 ALUs do not have

Everything above describes how a GPU’s general FP32 ALUs do matrix math, and they are already very good at it. Tensor cores are a separate, specialized class of execution unit that NVIDIA introduced with the Volta architecture, and they change the economics again.

A standard FP32 ALU executes one FMA per thread per instruction: two scalar inputs, one scalar output. To compute a matrix multiply you issue an enormous number of these scalar FMAs. A tensor core does not operate on scalars. It executes a small matrix-multiply-accumulate (MMA) as a single hardware operation, computing D = A x B + C where A, B, C, and D are small matrix tiles (for example 16x16 fragments), in one instruction. Instead of one FMA, a single tensor-core MMA performs hundreds of multiply-adds at once. The hardware is physically a systolic-style array of multiply-accumulate units wired together so that a tile of inputs streams through and a tile of partial sums accumulates, with vastly less instruction-issue and register-file overhead per FLOP than scalar FMAs incur. That overhead reduction is the entire point: at the scale of GEMM, the cost of fetching, decoding, and issuing scalar instructions and shuttling operands through the register file becomes a real tax, and the MMA instruction amortizes it across a whole tile.

The catch, and the honesty, is in the data types. Tensor cores deliver their headline throughput on reduced-precision inputs, FP16, BF16, FP8, and INT8, typically accumulating in FP32 to preserve accuracy. They are built for the precision profile of deep learning, where the inputs tolerate low precision but the running sum needs more. That is why tensor cores transformed AI training and inference: those workloads are dominated by GEMM and are precision-tolerant on the multiply inputs. They are not a free speedup for everything. If your problem genuinely requires full FP32 or FP64 precision on the multiplicands, classic HPC and scientific simulation, you fall back to the standard ALUs and get far less of the advertised peak (recent architectures added some tensor support for higher precision, but the biggest multiples are on the low-precision paths). So the precise statement is: tensor cores added a tile-granular MMA execution unit that the consumer-style scalar FP32 ALU pipeline does not have, and they win biggest exactly where the math is dense, the precision is relaxed, and the problem is GEMM-shaped. For the cases where you cannot or will not write to CUDA and tensor-core intrinsics directly, the portability landscape is its own topic, covered in GPU programming without CUDA.

Where CPUs still win, honestly

It would be dishonest to leave the impression that GPUs are simply better processors. They are better at one specific shape of problem and worse at most others. The architectural choices that make a GPU dominate at GEMM make it bad at everything GEMM is not.

Branchy, control-heavy code. A parser, a compiler, a database query planner, a chess engine’s search, anything dominated by unpredictable data-dependent branches, is poison for SIMT. Warp divergence means the GPU executes every taken path serially with most lanes idle. A CPU with a good branch predictor that is right 95-plus percent of the time runs straight through. The GPU’s lockstep warp is a liability the moment threads disagree.

Serial and latency-critical work. If you have one task and you need its answer as fast as possible, the CPU wins decisively. A GPU’s strength is throughput across many independent tasks, not the latency of any single one. The whole latency-hiding scheme only helps when there is a deep pool of other work to overlap; with a single thread, the GPU just eats the full memory latency with nothing to hide it behind. CPU clocks are higher, single-thread IPC is higher, and the out-of-order engine is built precisely to minimize the time to one answer.

Cache-friendly irregular work and pointer chasing. Traversing a linked list, walking a tree, or hopping around a graph generates dependent, scattered memory accesses. These cannot be coalesced (consecutive threads do not touch consecutive addresses), they cannot easily be vectorized, and each access depends on the previous one so there is little parallelism to hide latency. A CPU’s large multi-level caches and prefetchers handle this gracefully; a GPU’s bet on coalesced bulk bandwidth is exactly the wrong bet here.

Small problems and overhead. Launching a kernel, moving data across PCIe (or even within a unified memory space), and synchronizing all carry fixed costs. For a small matrix, the overhead dwarfs the compute, and the CPU finishes before the GPU has even started. GPUs amortize their overhead only at scale.

Workload trait Favors Why
Dense, regular, independent FLOPs GPU Saturates many ALUs, hides memory with warps
Unpredictable branches CPU No warp divergence penalty; real branch prediction
Single serial task, low latency CPU High clock, OoO, no need for parallel work to hide stalls
Pointer chasing / irregular access CPU Caches and prefetch; GPU cannot coalesce
Precision-tolerant GEMM (ML) GPU + tensor cores MMA throughput on FP16/BF16/FP8
Strict FP64 with little parallelism CPU (or GPU FP64 units) Limited tensor-core benefit, overhead may dominate

The two designs are not competitors so much as specialists, which is why modern systems pair them: the CPU runs the irregular control logic, orchestrates, and handles latency-sensitive paths, while it hands the GPU the big regular numeric kernels. That heterogeneous split, and the way different instruction-set philosophies coexist, echoes the broader point that there is rarely one architecture to rule them all, a theme explored in the RISC vs CISC war ended in a tie.

Verdict

The GPU does not beat the CPU at matrix math by having faster cores. It wins by refusing to care about single-thread speed at all. It commits its entire transistor budget to throughput: thousands of modest ALUs, a SIMT model that runs threads 32 to a warp under one program counter, warp schedulers that swap among an oversubscribed pool of warps at zero cost to hide hundreds of cycles of memory latency, and a memory hierarchy with a programmer-managed shared-memory tier that, combined with coalesced access, keeps those ALUs fed. Occupancy and register pressure, not raw FLOP counts, are what determine whether the design actually delivers, because the whole machine depends on having enough parallel work resident to cover its slow memory. Tensor cores then take the dominant operation, the matrix-multiply-accumulate, and bake it into dedicated tile-granular hardware that crushes the instruction-issue overhead of scalar FMAs, at the price of leaning on reduced precision.

Matrix multiplication is the perfect workload for every one of those choices: massively parallel, branch-free, regular in its memory access, arithmetically intense, and (for deep learning) tolerant of low precision. That alignment, not any general superiority, is why GPUs win here. Hand a GPU a branchy parser or a pointer-chasing graph traversal and it loses badly, because the same lockstep, latency-hiding, bandwidth-betting design that makes GEMM fly is precisely wrong for irregular serial work. The right mental model is not “GPU beats CPU.” It is “GPU and CPU are specialists, and matrix math happens to be the GPU’s home turf.” Build your systems to play each to its strength, and let the packaging and interconnect trends in chiplets and advanced packaging keep blurring the line between them.

Sources

Comments