SIMD, GEMM, and im2col Convolution

This is Lesson 4, Part 2 of the Deep Learning Inference Engineering series. Part 1 laid out nine kernel-optimization dimensions — loop ordering, blocking, cache reuse, SIMD, register blocking, prefetching, alignment, threading, and specialized instructions — as the toolbox a kernel author reaches for once a bottleneck has been diagnosed. This part goes deep on the two dimensions where the payoff is largest: SIMD, worked through an honest speedup calculation instead of a hand-wave, and GEMM, the single most important kernel in deep learning, together with the classic im2col trick that lets convolution ride on the same highly-tuned GEMM path as every fully-connected layer and transformer projection. SIMD and cache blocking both feed directly into that one kernel.

SIMD, Precisely

What a Vector Instruction Actually Does

A scalar instruction takes one (or two) operands and produces one result. A SIMD (Single Instruction, Multiple Data) instruction takes a vector register — a wide register packed with several scalar values laid out contiguously — and applies the same operation to every lane of that register in one instruction issue. The hardware doesn't run a loop internally; it has physically replicated ALU lanes wired to the same instruction decode, so all lanes compute in parallel within (in the common case) a single cycle of throughput.

Scalar multiply, 8 elements:
  c0 = a0 * b0        (instruction 1 → 1 result)
  c1 = a1 * b1        (instruction 2 → 1 result)
  c2 = a2 * b2        (instruction 3 → 1 result)
  ...
  c7 = a7 * b7        (instruction 8 → 1 result)
  → 8 instructions issued, 8 results produced
 
AVX2 vector multiply, 8 elements:
  [c0 c1 c2 c3 c4 c5 c6 c7] = [a0 a1 a2 a3 a4 a5 a6 a7] * [b0 b1 b2 b3 b4 b5 b6 b7]
  → 1 instruction issued, 8 results produced

The width of that register is fixed by the instruction set and the microarchitecture:

x86 AVX2: 256-bit YMM registers. A 32-bit float32 is 32 bits wide, so one YMM register holds 8 of them (256 / 32 = 8). The same register holds 4 float64 values, or 32 int8 values, since the lane count is just register-width divided by element width.
x86 AVX-512 (where available): 512-bit ZMM registers, 16 float32 lanes per instruction.
ARM NEON: 128-bit registers, 4 float32 lanes per instruction (or 2 float64, or 16 int8 — same divide-by-element-width logic).

For an edge-inference engineer, NEON's 4-wide fp32 lanes (or its dot-product and bfloat16 matrix-multiply extensions on newer ARM cores) are usually the ceiling; AVX2's 8-wide fp32 lanes are the desktop/server-class ceiling. Both are strictly about throughput — more useful work retired per instruction fetched and decoded — not about latency of any single operation.

The Naive Expectation, and Why It's Wrong

The naive mental model says: "AVX2 has 8 lanes, so my loop over N float32 elements runs 8 times faster once vectorized." That statement is true only about the arithmetic unit — never automatically true about the loop, because the arithmetic unit is rarely what a real loop is waiting on. This is exactly the compute-bound-versus-memory-bound distinction Lesson 1 built the Roofline model around, and it applies to a single vectorized loop just as much as it applies to a whole kernel.

Take one of the simplest loops that exists — a scaled vector accumulate, y[i] = a * x[i] + y[i], the classic AXPY pattern from linear algebra. Run the numbers for a concrete (illustrative, order-of-magnitude realistic) CPU core:

Clock: 3.0 GHz
2 FMA (fused-multiply-add) execution ports per core
AVX2, 8-wide float32 lanes
Sustainable single-core DRAM bandwidth: 24 GB/s

Peak compute, scalar (no SIMD):

3.0 × 10⁹ cycles/s × 2 FMA ports × 2 FLOPs/FMA = 12 GFLOP/s

Peak compute, AVX2 (8-wide):

3.0 × 10⁹ cycles/s × 2 FMA ports × 8 lanes × 2 FLOPs/FMA = 96 GFLOP/s

96 / 12 = 8× — exactly the lane width, as expected, because the compute ceiling really does scale linearly with lane count. So far the naive expectation holds. The problem is that AXPY's inner-loop body isn't compute-bound in the first place.

Per element, the loop touches three 4-byte float32 values: load x[i], load y[i], store the updated y[i] back — 12 bytes of DRAM traffic — to perform one FMA, which is 2 FLOPs (one multiply, one add). That's an arithmetic intensity of:

AI = 2 FLOPs / 12 bytes ≈ 0.167 FLOP/byte

Feed that into the Roofline formula from Lesson 1, Attainable = min(Peak_compute, AI × B_memory):

Memory ceiling = AI × B_memory = 0.167 FLOP/byte × 24 GB/s = 4.0 GFLOP/s

Now compare all three numbers side by side:

Quantity	Value
Scalar compute peak	12 GFLOP/s
AVX2 compute peak	96 GFLOP/s
Memory-bandwidth ceiling for this AI	4.0 GFLOP/s

The memory ceiling (4.0 GFLOP/s) is below both compute peaks. The scalar version was already memory-bound before a single AVX2 intrinsic was written — its 12 GFLOP/s of theoretical compute throughput was never reachable, because DRAM couldn't supply operands fast enough. Vectorizing to AVX2 raises the compute ceiling to 96 GFLOP/s, which makes the kernel more memory-bound, not less: the achievable throughput for both versions is pinned near the same ~4 GFLOP/s memory roof. The honest expected speedup from vectorizing this specific loop is close to 1×, not 8× — perhaps a modest 1.1–1.3× in practice from fewer instructions being issued and decoded (less loop-overhead and branch-prediction pressure per element), but nowhere near the lane-count multiplier.

This is the single most important caveat in this entire section, so it's worth stating as a rule rather than a one-off observation:

A SIMD instruction set raises the compute ceiling of the Roofline model. It does nothing to the memory-bandwidth ceiling. If a kernel's arithmetic intensity already puts it in the memory-bound region at scalar width, widening the vector register only makes the compute ceiling more irrelevant — the achieved speedup collapses toward 1× regardless of how many lanes the ISA offers.

Elementwise operators — activation functions, bias-add, residual-add, most normalization arithmetic — are exactly this shape: one or two loads, one store, a handful of FLOPs per element. They are the kernels where "we vectorized it and got nothing" is the expected outcome, not a bug.

Plotting this on the same style of Roofline chart Lesson 1 introduced makes the collapse visually obvious. Each width of the same core has its own ridge point (AI* = P_peak / B_memory): 0.5 FLOP/byte for the scalar core (12 / 24), 4.0 FLOP/byte for the AVX2 core (96 / 24). AXPY's 0.167 FLOP/byte sits to the left of both ridge points, meaning both versions are pinned to the same sloped memory-bandwidth line rather than either flat compute roof:

Both the scalar dot and the AVX2 dot for this kernel land at essentially the same point on the sloped memory-bandwidth line, far below either flat compute roof — which is the chart-form statement of the same conclusion reached numerically above: widening the vector register moved the ceiling, not the achieved point.

The Non-Vectorizable Tail

There is a second, smaller but very real, source of shortfall below the naive N× figure: the vector width almost never divides the problem size evenly.

Take a loop over N = 1,000,003 float32 elements on an AVX2 target (8-wide). The number of full vector iterations is:

1,000,003 / 8 = 125,000.375
→ 125,000 full 8-wide vector iterations, covering 1,000,000 elements
→ 3 leftover ("tail") elements that don't fill a register

Those 3 tail elements have to be handled somehow — either a scalar cleanup loop after the main vector loop, or a masked/predicated vector instruction (AVX-512's per-lane masking, or AVX2's maskload/maskstore) that processes a partial register while suppressing the lanes past the end of the array. Either way, that tail path:

can't benefit from the full vector width by construction (there's nothing to fill the other 5 lanes with),
often can't be branch-predicted as cleanly as the steady-state loop body, since it runs exactly once per call rather than thousands of times,
and, for small arrays relative to the vector width (a 1×1 pointwise convolution over a tiny channel count is a realistic edge-inference example), can dominate total runtime rather than being a rounding error.

For N = 1,000,003 the tail is 0.0003% of the work — genuinely negligible. For a depthwise convolution kernel iterating over C = 19 channels on a NEON target with 4-wide lanes, the tail is 19 mod 4 = 3 leftover channels out of 19 — nearly 16% of the loop trip count running at effectively scalar throughput. Small-N kernels are exactly where "SIMD gives you N× for free" breaks down hardest, and small-N kernels are exactly what edge-inference workloads are full of: small channel counts, small spatial tiles, batch size one.

Combine the two effects — memory-bandwidth ceilings that vector width doesn't touch, and remainder loops that can't fill a register — and the practical rule of thumb an inference engineer should carry is: treat the vector width as the theoretical ceiling and expect the achieved number to land well below it, until an arithmetic-intensity calculation says otherwise.

When SIMD Actually Earns Its Keep

That "until an arithmetic-intensity calculation says otherwise" clause matters, because SIMD genuinely does deliver close to full lane-width speedups — but only on kernels whose arithmetic intensity is high enough that the compute ceiling, not the memory ceiling, is what's binding. A kernel with an inner loop that reuses each loaded operand many times before evicting it — rather than touching each byte once and moving on, the way AXPY does — pushes arithmetic intensity up into the compute-bound region of the Roofline plot, where a wider vector register (or more FMA ports, or tensor cores) is exactly the lever that moves the needle.

There is exactly one kernel, more than any other in deep learning, engineered from the ground up to have that property. That's the subject of the rest of this post.

Lane Width Scales With Element Width, Not Just the ISA

One more piece of arithmetic is worth doing before leaving SIMD, because it connects directly back to Lesson 3's quantization material. AVX2's register is a fixed 256 bits, but the lane count is register width / element width, which means the lane count changes depending on what's being packed into the register, independent of the ISA itself:

256-bit AVX2 register:
  float32 (32-bit):  256 / 32 = 8 lanes
  int16   (16-bit):  256 / 16 = 16 lanes
  int8    (8-bit):   256 / 8  = 32 lanes

This is why quantizing a GEMM or convolution down to int8 doesn't just shrink memory traffic by roughly 4× relative to fp32 (Lesson 3's argument) — it also quadruples the number of elements one SIMD instruction can process, from 8 fp32 lanes to 32 int8 lanes on the same physical register, if the kernel's arithmetic intensity is high enough for that extra compute headroom to actually translate into throughput, per this section's memory-ceiling caveat. ARM's dot-product extensions push the same idea one step further at the instruction level: SDOT/UDOT accumulate four int8×int8 products directly into a single int32 lane in one instruction, rather than requiring separate multiply and add instructions per group of four — collapsing what would otherwise be several instructions into one, on top of the lane-count increase from going 8-bit. That combination — more lanes per register, plus fused multiply-accumulate-of-groups instructions — is a meaningful part of why quantized inference kernels on ARM mobile and edge cores routinely outrun a naive "4× fewer bytes" estimate.

GEMM: The Central Kernel of Deep Learning

What GEMM Is

GEMM — General Matrix Multiply — computes:

C = α·(A × B) + β·C

where A is M × K, B is K × N, and C is M × N, and α, β are scalars (usually α = 1, β = 0 for a plain matmul, or β = 1 to accumulate into an existing C, e.g. for a bias or residual). Every entry of the output is a dot product:

C[i][j] = Σₖ A[i][k] × B[k][j]     for k = 0 .. K-1

That's the entire definition. What makes GEMM worth an entire section of this series — arguably the most important section — is not the definition itself, but the fact that an enormous fraction of everything a deep neural network computes reduces exactly to this operation, and that reduction is not a coincidence or an implementation convenience. It falls directly out of what a linear (or linear-plus-nonlinearity) layer is.

Why a Fully-Connected Layer Literally Is a Matmul

A fully-connected (dense/linear) layer maps an input vector to an output vector via a learned weight matrix and bias:

y_j = Σₖ x_k × W[k][j] + b_j     for each output unit j

Stack a batch of B input rows into a matrix X of shape B × K (each row one example, K input features), and lay the layer's weights out as a K × N matrix W (K inputs, N output units). The entire batch's output, Y of shape B × N, is:

Y = X × W + b

with the bias broadcast across rows. Y = XW + b is not "analogous to" a matmul, or "computable via" a matmul — it is one, entry for entry, by the same dot-product definition given above. The FLOP count for the matmul term is 2 × B × K × N (a multiply and an add per MAC, summed over B × K × N output-times-reduction elements); the bias add is O(B × N), which for any layer with more than a handful of input features is dwarfed by the matmul term. A fully-connected layer's runtime, in other words, is GEMM's runtime, full stop.

Convolution Reduces to GEMM Too

A convolutional layer is not obviously a matmul — its inputs and outputs are indexed by spatial position, not just a flat feature index, and each output pixel only depends on a local receptive-field window rather than the entire input. But each individual output value is still a dot product: sum over the input channels and the kernel's spatial window of weight × input, exactly the same primitive operation as one row of an FC layer's output. The im2col technique, covered in full numeric detail in the next section, makes that structural fact explicit by physically rearranging the convolution's input into a matrix such that a single GEMM call produces the entire output feature map. This is why convolution belongs in the same section as fully-connected layers rather than getting treated as an unrelated primitive: after im2col, it is the same primitive, run through the same code path.

Transformers Are Mostly GEMMs Too

A transformer block's dominant cost, by FLOP count, is a handful of matmuls:

Q = X × Wq          # GEMM: (seq_len × d_model) × (d_model × d_k)
K = X × Wk          # GEMM
V = X × Wv          # GEMM
scores = Q × Kᵀ      # GEMM: (seq_len × d_k) × (d_k × seq_len)
context = softmax(scores / √d_k) × V   # softmax is elementwise/row-wise; the × V is a GEMM
output = context × Wo                   # GEMM

...followed by a feed-forward sublayer that's just two more FC layers (two more GEMMs) with a nonlinearity sandwiched between them. Softmax and layer normalization are real costs, but they're elementwise or row-reduction operations over data the surrounding GEMMs already produced — low-FLOP relative to the matmuls, in the same sense the FC layer's bias-add is low-FLOP relative to XW.

Why This Convergence Matters

Put the three cases together and the conclusion is structural, not coincidental: fully-connected layers are GEMMs by definition, convolutions become GEMMs via im2col (or an equivalent implicit formulation), and transformer attention and feed-forward sublayers are built almost entirely out of GEMMs. That means a hardware vendor or software stack that makes one operation — GEMM — extremely fast has, in one stroke, made the overwhelming majority of a modern network's FLOPs fast. This is precisely why GPU tensor cores, NPU MAC arrays, and every serious ML framework's compute backend (cuBLAS, oneDNN/MKL-DNN, XNNPACK, Arm Compute Library, BLIS, OpenBLAS) funnel almost everything through a small number of hand-tuned GEMM kernels rather than writing a bespoke kernel per operator type. Kernel engineering effort concentrated on GEMM has the highest possible leverage of any single optimization target in the entire inference stack.

Why GEMM, Specifically, Can Be Made Compute-Bound

The SIMD section above ended on the observation that SIMD only pays off on kernels with high arithmetic intensity. GEMM is the kernel where that condition is achievable by design, and it's worth being precise about why, because it's the same mechanism Lesson 2 used for tiling in general, now applied to the one operation that matters most.

An M × K × N GEMM performs 2MNK FLOPs. Naively — three nested loops, no attention to cache behavior — computing C[i][j] += A[i][k] * B[k][j] for every (i, j, k) triple re-reads A and B from memory on almost every access, because the working set doesn't fit in any cache and gets evicted before it's reused. Lesson 1 worked exactly this scenario for a 1024³ GEMM and found roughly 0.25 FLOP/byte — solidly memory-bound, wasting essentially all of that 96 GFLOP/s AVX2 compute ceiling from the previous section, because DRAM bandwidth caps the achievable throughput far below it, for exactly the same reason AXPY was capped.

Tiling (blocking) the same computation — restructuring the loop order so that a sub-block of A and a sub-block of B are loaded once and then reused many times against each other while they're still resident in a fast cache — changes the traffic-to-compute ratio dramatically. Lesson 1's tiled version of that same 1024³ GEMM reached roughly 170.67 FLOP/byte: over 680× higher arithmetic intensity than the naive version, moved from deep in the memory-bound region to comfortably compute-bound. GEMM's O(MNK) compute over O(MK + KN + MN) data is what makes this possible in the first place — unlike an elementwise op, where every byte read is used exactly once no matter how you order the loops, a GEMM's reduction dimension K gives you a genuine reuse opportunity: each element of A participates in N different output dot products, and each element of B participates in M of them. Blocking is the mechanism that actually captures that reuse before the data falls out of cache; it doesn't change the FLOP count or the mathematical result, it changes how many times each byte crosses each level of the memory hierarchy.

The Blocking Hierarchy: Registers → L1 → L2 → L3

The specific blocking scheme used by essentially every serious CPU GEMM implementation today — OpenBLAS, BLIS, Intel MKL, and the reference algorithm most of them trace back to — comes from Kazushige Goto and Robert van de Geijn's paper "Anatomy of High-Performance Matrix Multiplication," and its portable, extensible descendant, the BLIS framework (Van Zee & van de Geijn). The structure is often described as "five loops around a micro-kernel": the outer loops progressively slice A and B into blocks small enough to fit each level of the cache hierarchy, and the innermost micro-kernel does the actual arithmetic entirely out of registers.

for jc in steps of NC over N:               # L3 block of B (K x NC), packed & resident in L3
  pack B[:, jc:jc+NC]  →  B̃   (K x NC)
 
  for pc in steps of KC over K:              # shared reduction-dimension slab
    for ic in steps of MC over M:            # L2 block of A (MC x KC), packed & resident in L2
      pack A[ic:ic+MC, pc:pc+KC]  →  Ã  (MC x KC)
 
      for jr in steps of NR over NC:         # L1 micro-panel of B (KC x NR)
        for ir in steps of MR over MC:       # register micro-panel of A (MR x KC)
 
          # micro-kernel: accumulate an MR x NR tile of C entirely in registers,
          # streaming KC steps of a rank-1 update per step
          C[ir:ir+MR, jr:jr+NR] += Ã[ir:ir+MR, :] × B̃[:, jr:jr+NR]

Mapped onto cache levels, and echoing the design BLIS's own documentation describes — B's packed panel resident in L3, A's packed block resident in L2, a micro-panel of B streamed through L1, and C itself never packed at all, streamed directly between the register file and main memory a tile at a time:

Block	Shape	Lives in	What it buys
Register block	`MR × NR`	vector/FP registers	The `C` tile being accumulated; read and written zero times to any cache during the entire `KC`-step reduction — it's pinned in registers
L1 block	`KC × NR`	L1 cache	A micro-panel of `B`, reused once per `MR`-step of the register loop
L2 block	`MC × KC`	L2 cache	A packed block of `A`, reused across every `NC`-wide sweep of the L1/register loops
L3 block	`KC × NC`	L3 / last-level cache	A packed panel of `B`, reused across every `MC`-wide sweep of the L2/L1/register loops

A concrete (illustrative, order-of-magnitude) sizing pass makes the "fits in the cache" constraint tangible, using round cache sizes for a hypothetical modern desktop core — L1 = 32 KB, L2 = 512 KB, L3 = 8 MB, float32 throughout:

Register block, MR = NR = 8. The C tile is 8 × 8 × 4 bytes = 256 bytes. On an AVX2 target, one YMM register holds 8 float32 values, so one row of the tile fits exactly one register — 8 registers hold the whole accumulator, leaving headroom in a 16-register file for the A and B values streaming through the innermost loop.
L1 block, KC = 256. The micro-panel of B is KC × NR × 4 bytes = 256 × 8 × 4 = 8 KB — comfortably inside a 32 KB L1 alongside the streaming slice of A and whatever else shares the cache.
L2 block, MC = 128. The packed block of A is MC × KC × 4 bytes = 128 × 256 × 4 = 128 KB — fits inside a 512 KB L2 with room to spare for associativity slack and the L1 traffic passing through it.
L3 block, NC = 2048. The packed panel of B is KC × NC × 4 bytes = 256 × 2048 × 4 = 2 MB — fits inside an 8 MB L3.

These exact numbers vary by microarchitecture and by which BLIS configuration a library ships for a given target — the point isn't the specific constants, it's the method: pick each block dimension so the tile that has to stay resident at that level of the hierarchy actually fits, the same SRAM-budget-driven derivation Lesson 2 used to size a tile for a hypothetical NPU's on-chip memory, and the same register-blocking idea Part 1 named as one of the nine kernel-optimization dimensions, now shown as the innermost layer of a four-level hierarchy rather than a single flat trick.

The payoff of getting all four levels right simultaneously is exactly what turned 0.25 FLOP/byte into 170.67 FLOP/byte in Lesson 1's worked example: every element of the packed A block gets reused NC times (once per column of the current B panel) before it's evicted, and every element of the packed B panel gets reused MC times, instead of each element being touched once and discarded. That reuse is what raises arithmetic intensity into the compute-bound region — which is precisely the region where the previous section's lesson flips: a wider SIMD register (or more FMA ports, or a systolic MAC array) now translates almost directly into a proportional speedup, because the bottleneck genuinely is the arithmetic units, not the memory bus.

Making the Reuse Count Explicit

It's worth deriving that reuse mechanism directly, rather than only citing the before/after numbers, because it's a simple counting argument once it's laid out one level at a time.

Start at the innermost level, the register-block/L1 pairing. For one step of the KC-deep reduction, the microkernel performs a rank-1 update: it multiplies one column of the A micro-panel (MR values) against one row of the B micro-panel (NR values) and accumulates the resulting MR × NR outer product into the register-resident C tile:

C[ir:ir+MR, jr:jr+NR] += A[:, k] ⊗ B[k, :]     (outer product, one step of k)

Drawn out for a small MR = NR = 4 tile, one step of that outer product looks like this — a column of A values broadcast down, a row of B values broadcast across, and every pairwise product landing in its own register-resident C entry:

Repeat that for each of the KC values of k in the reduction slab, accumulating into the same 16 register locations every time rather than writing back to memory between steps, and the tile's final value is the correct sum over the full KC-deep reduction — with only 2 × KC × (MR + NR) total register loads for 2 × KC × MR × NR FLOPs of useful work.

Loading MR values of A and NR values of B — MR + NR loads total — produces MR × NR output entries, each requiring one multiply and one add, so 2 × MR × NR FLOPs. The ratio of FLOPs to register loads for this single step is 2 × MR × NR / (MR + NR), which grows with the tile dimensions: double MR and NR and the FLOP count quadruples while the load count only doubles. That's the entire mechanism blocking exploits, made numeric: promote every operand from "read once, used once" to "read once, used MR or NR times," and the register file's size is the only thing capping how large MR × NR can get.

The same counting argument repeats one level up. The packed MC × KC block of A sitting in L2 is not reloaded from L3 or DRAM once per micro-panel of B that sweeps past it during the NC-wide inner sweep — it is packed into contiguous, cache-resident storage exactly once, and then reused for the entire sweep. That's precisely how the illustrative 128 KB L2-resident block of A computed earlier contributes to an entire 2 MB L3-resident panel of B's worth of output columns without a single additional trip to DRAM for A. Multiply the reuse factor through all four levels — registers, L1, L2, L3 — and the aggregate effect is exactly the jump Lesson 1 measured for the 1024³ case: from 0.25 FLOP/byte (every operand reloaded from DRAM roughly every time it's touched) to 170.67 FLOP/byte (every operand loaded from DRAM once, then reused dozens to hundreds of times out of registers, L1, and L2 before it's evicted).

im2col: Turning Convolution Into GEMM

The Idea

A 2D convolution slides a small kernel over a larger input and, at every valid position, computes a dot product between the kernel and the local patch of input it currently overlaps. im2col ("image to column") makes that dot-product structure explicit and mechanical: extract every overlapping patch the kernel will ever be applied to, flatten each patch into a single row (or column) of a matrix, and stack all of those rows together. The result is a matrix whose rows are exactly the input vectors a GEMM needs, and the convolution's kernel weights, flattened the same way, become the matrix a GEMM multiplies against. One matrix multiply then produces every output position at once.

A Fully Worked Numeric Example

Take a 4 × 4 single-channel input I and a 3 × 3 kernel W, both small enough to trace by hand. Convolution here means the deep-learning convention (cross-correlation, no kernel flip), stride 1, no padding, so the output is (4 − 3 + 1) × (4 − 3 + 1) = 2 × 2.

I =                          W =
 1   2   0   1                1   0   1
 3   1   2   0                0   1   0
 0   1   3   2                1   0   1
 2   0   1   3

Step 1 — direct convolution, computed the ordinary way. For each output position, take the 3 × 3 patch of I it aligns with, multiply elementwise with W, and sum.

Output position (0,0) — patch is rows 0–2, columns 0–2:

patch =        elementwise × W =         sum:
 1  2  0        1·1  2·0  0·1             1 + 0 + 0
 3  1  2        3·0  1·1  2·0     →       + 0 + 1 + 0    = 5
 0  1  3        0·1  1·0  3·1             + 0 + 0 + 3

Output position (0,1) — patch is rows 0–2, columns 1–3:

patch =        elementwise × W =         sum:
 2  0  1        2·1  0·0  1·1             2 + 0 + 1
 1  2  0        1·0  2·1  0·0     →       + 0 + 2 + 0    = 8
 1  3  2        1·1  3·0  2·1             + 1 + 0 + 2

Output position (1,0) — patch is rows 1–3, columns 0–2:

patch =        elementwise × W =         sum:
 3  1  2        3·1  1·0  2·1             3 + 0 + 2
 0  1  3        0·0  1·1  3·0     →       + 0 + 1 + 0    = 9
 2  0  1        2·1  0·0  1·1             + 2 + 0 + 1

Output position (1,1) — patch is rows 1–3, columns 1–3:

patch =        elementwise × W =         sum:
 1  2  0        1·1  2·0  0·1             1 + 0 + 0
 1  3  2        1·0  3·1  2·0     →       + 0 + 3 + 0    = 7
 0  1  3        0·1  1·0  3·1             + 0 + 0 + 3

Direct-convolution output:

5   8
9   7

Step 2 — build the im2col matrix. Flatten each 3 × 3 patch, row-major, into one row of a matrix. There are 4 output positions, and each patch has 3 × 3 = 9 elements, so the im2col matrix is 4 × 9:

             col0 col1 col2 col3 col4 col5 col6 col7 col8
row (0,0):    1    2    0    3    1    2    0    1    3
row (0,1):    2    0    1    1    2    0    1    3    2
row (1,0):    3    1    2    0    1    3    2    0    1
row (1,1):    1    2    0    1    3    2    0    1    3

Flatten the kernel W the same way, row-major, into a 9 × 1 column vector:

w = [1, 0, 1, 0, 1, 0, 1, 0, 1]ᵀ

Step 3 — GEMM. Multiply the 4 × 9 im2col matrix by the 9 × 1 weight vector:

row (0,0) · w = 1·1 + 2·0 + 0·1 + 3·0 + 1·1 + 2·0 + 0·1 + 1·0 + 3·1 = 5
row (0,1) · w = 2·1 + 0·0 + 1·1 + 1·0 + 2·1 + 0·0 + 1·1 + 3·0 + 2·1 = 8
row (1,0) · w = 3·1 + 1·0 + 2·1 + 0·0 + 1·1 + 3·0 + 2·1 + 0·0 + 1·1 = 9
row (1,1) · w = 1·1 + 2·0 + 0·1 + 1·0 + 3·1 + 2·0 + 0·1 + 1·0 + 3·1 = 7

GEMM output, reshaped back to the 2 × 2 spatial grid:

5   8
9   7

This matches the direct-convolution result computed in Step 1, entry for entry — exactly as it must, since every row of the im2col matrix and every step of the GEMM is, mechanically, the same multiply-and-sum as the direct patch computation. Nothing about the numerical result changes; im2col is purely a data-layout transformation that lets a general-purpose, exhaustively-optimized GEMM kernel run the computation instead of a bespoke convolution loop.

Generalizing to Multiple Input Channels

The example above used one input channel to keep the arithmetic traceable by hand. A real conv layer has C_in input channels, and it's worth checking on a small example that going multi-channel doesn't introduce anything new — im2col just adds channel as another axis to flatten into the same row.

Take a 2 × 2 input with 2 channels, and a 2 × 2 kernel with 2 channels (so the kernel exactly covers the input — only one valid output position, which keeps the check small):

Channel 0 input:      Channel 1 input:      Channel 0 kernel:   Channel 1 kernel:
 1   2                  5   6                 1   0               1   1
 3   4                  7   8                 0   1               1   1

Direct convolution, summed across both channels: channel 0 contributes 1·1 + 2·0 + 3·0 + 4·1 = 5; channel 1 contributes 5·1 + 6·1 + 7·1 + 8·1 = 26. Total output = 5 + 26 = 31.

im2col, with the single output position's patch flattened channel-by-channel, row-major within each channel:

im2col row = [ 1  2  3  4 | 5  6  7  8 ]     (channel 0 patch, then channel 1 patch)
 
w           = [ 1  0  0  1 | 1  1  1  1 ]    (channel 0 kernel, then channel 1 kernel)

GEMM (a 1 × 8 row times an 8 × 1 column, this time — one output position, C_in × kernel_h × kernel_w = 2 × 2 × 2 = 8 reduction elements):

row · w = 1·1 + 2·0 + 3·0 + 4·1 + 5·1 + 6·1 + 7·1 + 8·1 = 1 + 0 + 0 + 4 + 5 + 6 + 7 + 8 = 31

31 matches the direct-convolution sum exactly. The only thing that changed versus the single-channel example is the row length — kernel_h × kernel_w × C_in instead of just kernel_h × kernel_w — and correspondingly the weight side becomes a (kernel_h × kernel_w × C_in) × C_out matrix rather than a single column once there's more than one output filter, so one GEMM call produces every output channel at every spatial position simultaneously. Nothing about the reduction itself changes shape; channels just fold into the same dot product as more terms.

The Memory Tradeoff

Look at what im2col did to the data size in the toy example: the original input I has 4 × 4 = 16 elements, but the im2col matrix has 4 × 9 = 36 elements — a 2.25× blowup, and every overlapping element got physically duplicated into multiple rows (the value 2 at position I[0][1], for instance, appears in three different im2col rows, since it falls inside three different 3 × 3 receptive fields). That duplication is the direct, mechanical cost of turning a sliding-window operation into a flat matrix: overlapping patches share input pixels, and im2col materializes each patch independently, so shared pixels get copied once per patch they participate in.

The toy example's 2.25× ratio understates the effect for realistic layer sizes, because the output there (2 × 2 = 4 positions) is tiny relative to the kernel (3 × 3 = 9 elements) — most input pixels don't get the chance to appear in many patches before the input runs out. Scale up to a size that's actually representative of an early CNN layer: a 224 × 224 × 3 input, 3 × 3 × 3 kernel, stride 1, padding 1 (so the output stays 224 × 224 spatially):

Original input:      224 × 224 × 3        = 150,528 elements  (≈ 588 KiB at fp32)
im2col matrix:        (224 × 224) rows × (3 × 3 × 3) cols
                     = 50,176 × 27
                     = 1,354,752 elements  (≈ 5.17 MiB at fp32)
 
Blowup factor ≈ 1,354,752 / 150,528 = 9.0×

That 9.0× is not a coincidence — it's exactly kernel_h × kernel_w (3 × 3). With stride-1, padding-preserving convolution, the number of output positions is approximately equal to the number of input pixels, so the ratio of im2col-matrix size to input size converges to the number of positions each interior pixel gets duplicated into, which is the kernel's spatial footprint. A 5 × 5 kernel would blow the input up by roughly 25×; a 7 × 7 kernel by roughly 49×. For large feature maps and larger kernels, that additional memory traffic — writing the unfolded matrix out to DRAM and reading it back in for the GEMM — can itself become the bottleneck, exactly the kind of self-inflicted memory-bound regime the SIMD section warned about, except now the extra bytes come from the data-layout transformation rather than the underlying math.

This is the classic tradeoff im2col makes, and it's worth stating plainly:

im2col trades memory for the ability to reuse an already highly-optimized GEMM kernel instead of writing and tuning a bespoke convolution loop. It is a strictly software-engineering decision — do the extra data movement and inherit BLIS/cuBLAS-grade GEMM performance for free, or skip the data movement and hand-roll a direct-convolution kernel that has to independently earn back everything GEMM's blocking hierarchy already gives you.

Implicit im2col: Getting the Reuse Without the Blowup

Production inference stacks rarely materialize the full im2col matrix in memory, precisely because of the blowup just computed. NVIDIA's cuDNN and CUTLASS libraries instead implement what they call implicit GEMM: the convolution is still executed as a matrix multiply — same loop structure, same blocking hierarchy, same tensor-core or MAC-array utilization as a real GEMM — but the "im2col matrix" is never written out as a physical buffer. Instead, the GEMM's inner loop computes each element's address on the fly, reading directly from the original, un-duplicated input tensor and reconstructing exactly the value the explicit im2col matrix would have held at that position, at the moment the micro-kernel needs it. The library gets im2col's central benefit — reformulating convolution as the operation every optimization dollar has already been spent on — without paying for the redundant copies of every overlapping pixel. This is exactly the compiler-level move Lesson 2 called operator fusion: avoid materializing an intermediate tensor entirely by folding its construction into the consumer's memory-access pattern, applied here to a specific, extremely high-value intermediate — the unfolded convolution matrix — rather than a generic elementwise chain.

Quick Reference

Concept	Key relationship
SIMD lane count	`register width / element width` — 8 for AVX2 fp32, 4 for NEON fp32
SIMD compute peak	scales linearly with lane count
SIMD memory ceiling	unaffected by lane count — set by `AI × B_memory`
Realistic SIMD speedup	`min(lane-count speedup, memory-ceiling speedup)` — often far below lane count for low-AI kernels
Non-vectorizable tail	`N mod lane_width` elements handled at scalar/masked throughput
GEMM FLOP count	`2 × M × N × K`
Why FC is GEMM	`Y = XW + b`, by the layer's own definition
Why conv becomes GEMM	im2col reduces each output pixel to a dot product, same primitive as an FC row
GEMM blocking hierarchy	register block (`MR×NR`) → L1 (`KC×NR`) → L2 (`MC×KC`) → L3 (`KC×NC`)
im2col memory blowup (stride 1, same-size output)	approaches `kernel_h × kernel_w`
Implicit im2col	computes im2col addresses on the fly inside the GEMM, no materialized buffer