4 — Memory & batching

Batching turns memory-bound into compute-bound

Emmanuel Pilliat — ENSAI 3A

Goals

• How long the CPU waits for memory.
• Loop order can cost more than the maths.
• Batching: more work per byte loaded, on the CPU.

Today’s enemy: the distance to your data.

▶ Script: 4_memory_and_gpu.jl (§1–3) · Notebook: Memory & GPU · 🐍 numpy (Colab)

⬅ 3 — Threads & concurrency · 5 — GPU ➡

Compilation made one core fast; threads used many. Today’s enemy is neither. One script and one notebook cover this deck AND the GPU deck.

Part A — the CPU waits for memory

The memory hierarchy

• Add two floats: < 1 ns. Fetch from RAM: ~100 ns.

Level	Size	Latency	If L1 = 1 s
registers	~100s B	~0	instant
L1	32–64 KB / core	~1 ns	1 s
L2	256 KB–1 MB	~4 ns	4 s
L3	8–32 MB	~15 ns	15 s
RAM	16–64 GB	~100 ns	~2 min

While waiting, the core could have done hundreds of operations. This machine: L1d 544 KiB total (14 instances), L2 18 MiB (9 instances), L3 24 MiB shared. Hybrid CPU (P/E/LP-E cores) — per-core L1 differs, hence “instances”. lscpu shows it; the script asks students to compare with their own.

Caches bet on locality

• Temporal: used now → used again soon.
• Spatial: used now → neighbors next.
• Memory moves in lines of 64 B = 8 Float64.

Fast code honors that bet.

Bet before you look 🎲

Sum the same array many times, from 8 KB to 64 MB.

L1 ≈ 1 ns, RAM ≈ 100 ns.

How much slower per access at 64 MB?

(~100×? Write it down.)

▶ Script: cell “1. The memory hierarchy — the cache cliff”

The cache cliff, measured

for _ in 1:passes                 # repeated_sum(a, passes)
    @inbounds @simd for i in eachindex(a)
        s += a[i]
    end
end

footprint	ns / access	lives in
8 KB	0.049	L1
4 MB	0.128	L2 (18 MiB)
64 MB	0.317	RAM (> 24 MiB L3)

Steps visible. But 6.5×, not ~100×. Why?

Absolute values are machine- and load-sensitive: run at rest and look for the JUMPS, not the exact ns. Full sweep (8 KB → 64 MB) is in the notebook/script.

Why 6.5×? Bandwidth vs latency

• Table = wait for one unpredictable access (latency).
• Our loop reads in order → hardware loads next lines early.
• We pay the flow rate (bandwidth), not the wait.

100× is real when the next address is unpredictable (linked list, random indices).

Not “the table is wrong” nor “machine noise”. The part of the hardware that loads lines early is the prefetcher. Latency and bandwidth are different resources; a streaming sum is bandwidth-bound, a linked-list walk is latency-bound.

Two walks, same maths

Sum an 8000×8000 matrix, two loop orders:

# COLUMN by column
for j in axes(M, 2)
    for i in axes(M, 1)
        s += M[i, j]
    end
end

# ROW by row
for i in axes(M, 1)
    for j in axes(M, 2)
        s += M[i, j]
    end
end

• Julia: column-major (M[i,j], M[i+1,j] adjacent).
• numpy: row-major. Your reflexes are reversed.

▶ Script: cell “2. Locality — the same sum, walked two ways”. Predict first.

11.9× for swapping two loops

8000² Float64 = 512 MB ≫ cache:

walk	time	each 64 B line…
columns	30.4 ms	fully used (8 values)
rows	361.4 ms	one value, then evicted

• Same result: 11.9× slower.
• Rule: first index in the inner loop.

Wrong order: right result, silently slow.

Prove it’s the cache 🔬

If Julia just “preferred columns”, the gap would show at any size.

Shrink the matrix. When does the gap vanish?

▶ Script: cell “Proof it IS the cache”

Prove it’s the cache — debrief

size	footprint	rows / columns
100²	78 KB — fits	0.9×
4000²	122 MB — spills	7.2×

• Gap appears when data outgrows the cache.

Not “columns are fast”: the cache decides.

Small matrix: everything is already close, so the order cannot matter. Strongest argument of the module: it rules out every explanation except the cache. Pluto: interactive slider. Script: full table (100² → 4000²).

numpy hides it, doesn’t remove it

C-order vs F-order arrays, 8000²:

numpy	ratio	why
A.sum(axis=0)	≈ 1×	numpy reorders
copy F→C / C→C	13.2× (109.5 / 8.3 ms)	can’t reorder a copy

numpy hides it; Julia shows it. Neither removes it.

On sum, someone wrote the reordering for you. A mismatched copy or a hand-written loop: the full penalty is back. Measured: sum(axis=0) C/F = 1.14× at 8000², 0.77× at 2000² — noise around 1×, and load-sensitive (0.45× seen under load). Don’t over-read either side of 1×. Demo lives in the Python notebook (python/4).

Allocating is memory traffic too

z = 2 .* x .+ y      # NEW 76 MiB array
z .= 2 .* x .+ y     # reuses z

▶ Script: cell “Allocating is memory traffic too”

n = 10⁷	time	allocated
z = …	54 ms	76 MiB
z .= …	17 ms	0 bytes

One dot: ~3×. But measure: the gain isn’t always there.

A new array is fresh memory to fetch and fill. Same for BLAS: A * B allocates, mul!(C, A, B) reuses C. At n = 2000 on this laptop the gain swings from no gain to ~3× from run to run (hybrid cores): measure. numpy gains only ~1.4× with out=: it can’t FUSE 2*x + y into one pass. For matrix × VECTOR the output is tiny and the gain vanishes. Measured at rest: mul!(y, W, x) on matrix-VECTOR was actually SLOWER than W * x (760 vs 468 µs). On the GPU, in place gains nothing in Julia either (2.10 vs 1.82 ms; CUDA.jl allocates from a pool) — whereas torch in place does gain (1.35 → 0.80 ms). The law is “allocation = memory traffic”.

Part B — batching: memory-bound → compute-bound

Arithmetic intensity: FLOP per byte

n×n, Float64	FLOP	bytes	intensity
mat × vec	2n²	8n²	0.25 — constant
mat × mat	2n³	8n²	n/4 — grows

• Low: compute waits for loads → memory-bound.
• High: each value reused → compute-bound.

Compare GFLOP/s, not seconds.

▶ Script: cell “3. Arithmetic intensity”

Work per byte loaded; a FLOP = one + or *. Don’t compare seconds: they do different work. GFLOP/s = 10⁹ FLOP per second. Intensity is per BYTE, not per value (8× apart for Float64) — be consistent. BLAS-3 (GEMM) tiling is what realizes the n/4: each tile is reused from cache.

Batching = mat×vec → mat×mat

• A layer = matrix product. W * x: each weight used once.
• B examples as columns of X: W * X reuses each weight B times.

y = W * x     # 1 example
Y = W * X     # B examples

Not a GPU trick: it raises intensity. Must pay on the CPU.

Batching on the CPU: bet 🎲

One layer W (4096² Float32), B examples as columns.

What happens to the time per example?

▶ Script: cell “Batching: one layer W, B examples”

Batching, measured on the CPU

B	total	µs / example	GFLOP/s
1	2.55 ms	2551.8	13.1
8	5.65 ms	706.0	47.5
64	11.06 ms	172.9	194.1
256	39.06 ms	152.6	219.9
1024	141.80 ms	138.5	242.3

18.4× per example, 13.1 → 242.3 GFLOP/s. No GPU.

Same CPU, far more of it used. The GPU deck pushes the same idea further. Honesty check (measured): B=1 here is a skinny GEMM (W * X[:,1:1]) at 2.55 ms; a true gemv (W * x) is 1.09 ms — the B=1 baseline is 2.3× pessimistic. Against the gemv baseline the per-example gain is still ~7.9×. Say it if asked. numpy (python/4) gives almost the same numbers: 2101 → 113.6 µs (18.5×).

Why B = 1024 barely helps any more

2551.8 → 172.9 (B=64) → 152.6 (B=256) → 138.5 µs: flattens.

• W = 64 MiB > 24 MiB L3: a fixed RAM stream per batch.
• Batching divides it by B:

\[\text{time per example} \approx \text{compute} + \frac{\text{stream } W}{B}\]

• B ≈ 256: second term tiny.

Flat = success: now compute-bound.

242 GFLOP/s is what this run reached, not the machine’s peak. The 1/B model fits: 138.5 µs ≈ pure compute per example; 2551.8 ≈ that plus the whole W stream. This is why python/4 uses n=4096: at 2000² W (16 MB) fits in L3 and the mechanism disappears.

Going further

• ▶ Script: scripts/4_memory_and_gpu.jl — §1–3 here, §4–7 GPU deck.
• ▶ Notebook: Memory & GPU, same split.
• 🐍 numpy: python/4_memory_and_gpu_python.ipynb
• ✏️ Activity: Julia · 🐍 Python (Colab)
• 🎯 Batching is not a GPU trick: ~18× on a laptop CPU.

⬅ Previous: 3 — Threads & concurrency · Next: 5 — GPU ➡