Why is NumPy 100x faster?

Adding up ten million numbers in a plain Python loop takes about a tenth of a second. NumPy does it in about a millisecond. This lesson is about where that factor of a hundred comes from, and the answer is only partly “NumPy is written in C”. The rest is about where the numbers sit in memory and the order you visit them in.

That second part follows you everywhere Python’s speed does not. It decides whether a feature pipeline over your ranker’s scores takes one second or fifteen, why a matrix walked the wrong way round is six times slower in any language, and why PyTorch code is full of calls that rearrange a tensor before the real work starts. On a GPU, which is built to be fed huge amounts of memory per second, it matters even more.

It builds on How slow is slow?, where you measured the layers of fast memory and how long each one takes, and Where do my variables live?, where you saw that a Python list is a row of pointers to objects scattered around memory. It hands the interpreter’s share to Why is Python slow? and the memory share to Why ~17 tokens per second?.

Today

• Split “NumPy is faster” into its two real causes: the cost of Python running bookkeeping for every element, and the cost of fetching the numbers from memory. Say how big each one is on your machine.
• Explain why reading one number from memory actually brings in the 64 or 128 bytes around it, and how the chip guesses what you will read next and fetches it early.
• Predict, before running, how a Python list, an array.array, a NumPy array, a NumPy array read with gaps, and a NumPy array read in random order will compare, and say why.
• Explain the two ways to lay a grid of numbers out in memory (row after row, or column after column), and show a 6x to 30x slowdown from walking it in the wrong order, in Python compiled with Numba.
• Explain, with a picture you could draw on a whiteboard, how a CPU core does several additions at once, and why a plain loop that adds into one running total can’t.
• Carry it to the GPU: why it matters that neighbouring GPU threads read neighbouring addresses, and why .contiguous() shows up all over PyTorch code.

Prerequisites: How slow is slow? (the cache ladder, locality, latency vs bandwidth) and Where do my variables live? (a list is an array of pointers to objects scattered on the heap; a NumPy array is one block of raw numbers). This lesson does not re-teach the ladder. It uses it. All code is Python; the tight loops use Numba, which the Build section introduces. Rough time: one evening, plus half an hour on the Spark.

The Problem

You have 10 million relevance scores from the cross-encoder, one float each, and you want their sum (or their mean, or a softmax denominator; same loop). You write the obvious thing:

s = 0.0
for x in scores:
    s += x

On an M4 Pro MacBook this takes 107 ms. np.sum on the same numbers in a NumPy array takes 1.2 ms. That is the famous gap: about 90x.

Ask around and you’ll hear “NumPy is written in C”. True, and only half the story. Here are two more measurements from the same machine, same 10 million numbers, same Python version:

• Keep the Python loop exactly as written, but build the list so its float objects sit in memory in shuffled order (same values, same list length, same bytecode). The loop now takes 457 ms. Four times slower, and the interpreter is doing exactly the same work.
• Keep NumPy, but sum the numbers in random order with a[perm].sum(). It takes 42 ms. Thirty-five times slower than NumPy in order, and it’s all C.

So there are two different bills. One is paid by running Python bytecode for every element. The other is paid by the memory system, and it charges by where your bytes sit and in what order you touch them. NumPy wins on both, which is why its lead looks so big. But you can lose the memory bill in C, in NumPy, and in PyTorch on a GPU, and then “it’s written in C” saves you nothing.

This lesson takes the two bills apart and measures each one.

The Idea

Two bills

Think of summing numbers as a warehouse job. Someone has to walk to each shelf and bring the item back (the memory bill), and someone has to fill in a form for every item (the interpreter bill). NumPy is a forklift on a straight aisle: no forms, and the items are lined up in a row. A Python loop over a list fills in a form for each item, and the items are wherever there was space when they were stored.

The form first. When CPython runs for x in scores: s += x, it doesn’t add two numbers. For every element it runs a handful of bytecodebytecodeThe simple instructions CPython compiles your Python into (FOR_ITER, STORE_FAST, BINARY_OP, …) and then executes one at a time in a big loop written in C. instructions: fetch the next item from the list, store it in x, load s and x, check both of their types, find the float-add routine, add the two doubles, and allocate a brand-new float object for the result (Python floats can’t be changed in place), then drop the old s and update a few reference countreference countThe counter in every CPython object of how many references point at it. Touched on almost every operation; zero frees the object.s. That costs around 10 ns per element on a fast laptop. We’ll call this interpreter overheadinterpreter overheadThe work CPython does per operation that isn’t your actual math: dispatching bytecodes, checking types, boxing results into objects, refcounting. Usually 10-50x the cost of the arithmetic itself.. Why is Python slow? takes that apart instruction by instruction. Here we only need its size.

Now the walk. The add itself needs the 8 bytes of a double, and those 8 bytes have to come from somewhere on the ladder from How slow is slow?: L1 at about a nanosecond, RAM at about a hundred. The memory bill is decided by which rung each load hits, and that is decided almost entirely by layout and order.

Here are the numbers for one element, measured on an M4 Pro (CPython 3.14, NumPy 2.5, 10 million float64s, best of several runs). You’ll get your own in the Build.

Read the table in pairs. The first two rows run identical bytecode, so the 35 ns difference is pure memory bill. sum(list) does the same walk as the loop with no bytecode per element, so loop minus sum (about 9 ns) is roughly the interpreter bill. And when the objects are scattered, even the C loop in sum slows down by 6x. So a fair summary of the famous 90x: about 5-10x of it is the interpreter, and the rest is memory layout and what layout lets the CPU do. Those last words matter, and we’ll come back to them with SIMD.

The cache line

You met the cache linecache lineThe fixed-size block that moves between RAM and the caches as one unit. 64 bytes on x86 and on most ARM server cores, including the Spark’s. 128 bytes on Apple M-series. in the latency lesson. Here it becomes the main character, so let’s look at it closely.

The memory system never moves one byte, or one double. When your code loads 8 bytes that aren’t in cache, the hardware fetches the whole aligned line those 8 bytes sit in. On your Mac that’s 128 bytes (check: sysctl hw.cachelinesize prints 128). On the Spark and on Intel and AMD chips it’s 64.

Two habits make this pay off, and they have names.

• Spatial localitySpatial localityUsing data that sits next to data you just used. A line holds 8 or 16 doubles, so reading neighbours means most loads hit the line you already paid for.: after a[5], you read a[6]. It’s in the line already.
• Temporal localityTemporal localityUsing the same data again soon. If you come back before the line gets pushed out by other data, the second visit is a cache hit.: after reading a[5], you read it again soon, before newer data pushes the line out of the cache.

Every speedup in this lesson is one of those two, or the third trick below.

The prefetcher: a clerk who reads ahead

Spatial locality alone would still leave you paying full RAM latency once per line: 100 ns every 16 doubles is about 6 ns per double, much slower than the 0.12 ns NumPy gets. The missing piece is the hardware prefetcherhardware prefetcherCircuitry next to the caches that watches the stream of addresses that miss, spots a pattern (next line, or every Nth line), and requests lines before the program asks for them..

Picture a clerk at the warehouse door watching what you request. You ask for shelf 40, then 41, then 42. Without being told, the clerk sends runners for 43, 44, 45, 46. By the time you ask for 43, it’s on the counter. The clerk can also spot a fixed stride (40, 48, 56…). What the clerk cannot do is guess the next item in a random list, or the next object behind a pointer it hasn’t loaded yet.

So in-order access converts a latency problem into a bandwidth problem: many lines in flight at once, arriving back to back, which is exactly the “bytes in flight” idea from the latency lesson. Random access leaves you with latency. On many chips the prefetcher also stops at page boundaries (4 KB on x86 and usually on the Spark’s Linux; 16 KB on macOS), since the next page might not even be mapped.

Feel it

Every pattern below reads exactly 16 numbers. Watch how many lines it costs. Flip the line size between 64 B and 128 B, and try both Python list layouts.

Things to notice:

1. In order costs one line per 16 numbers on the Mac (one per 8 on x86). 100% useful.
2. Stride 8 is a 64-byte step. With 64 B lines that’s one number per line, 1/8 useful: you need 8x the memory traffic for the same work. With 128 B lines it’s half as bad. The same code has a different penalty on the Mac and the Spark.
3. Column is the column walk of a matrix (more on this below): one line per number, every time.
4. Random costs about one line per number and the prefetcher can’t help.
5. Python list in order is not that bad for memory! More bytes than NumPy (an 8-byte pointer plus a 32-byte object block per number), but they’re all in a forward walk. Shuffle it and it turns into the random case.

Why pointer-chasing lists lose

From Where do my variables live?: a list of floats is an array of 8-byte pointers, and each pointer leads to a separate 24-byte float object somewhere on the heap. Reading element i is two loads: the pointer, then the object it points to. The pointer array is dense and walked in order, so it’s cheap. The objects are where it gets interesting.

When objects are created in order (scores = a.tolist(), a list comprehension, appending in a loop), CPython’s small-object allocator tends to hand out neighbouring addresses one after another. Object i+1 usually sits right after object i. Following the pointers walks forward through memory, and the prefetcher is happy.

When the visiting order doesn’t match the creation order (you sorted the list, shuffled it for a training epoch, built it out of dict values, or the objects were made at different times by different code), each pointer jumps somewhere new. Now every element is a likely cache miss, and the prefetcher can’t predict the next address because it doesn’t know it until the pointer is loaded.

Why 45 ns and not the full ~100 ns RAM latency? The pointers come from the in-order array, so the CPU knows the next few object addresses early and can have a couple of misses in flight at once. That overlap is called memory-level parallelismmemory-level parallelismHaving several cache misses outstanding at the same time, so their waits overlap. Out-of-order cores can track a few dozen per core, if they can find independent loads to issue.. But the interpreter runs a long stretch of instructions between one object load and the next, which limits how far ahead the core can look, so it only overlaps a few. (That’s my best explanation for the exact number. The direction, a big slowdown from scattering alone, is the robust part.)

This is exactly the thing NumPy removes. A NumPy array has no per-element objects and no pointers: element i is at base + 8*i. The address of every element is known without loading anything. Which leads to the question of what happens when NumPy’s elements are not next to each other.

Rows, columns and strides

RAM is one long line of bytes. A matrix is two-dimensional. Something has to decide how to flatten it. C, NumPy (by default) and PyTorch use row-major orderrow-major orderStore a matrix one row after another: all of row 0, then all of row 1. Moving along a row moves to the next address. Also called C order.: row 0, then row 1, and so on. Fortran, MATLAB, Julia and many BLAS libraries use column-major ordercolumn-major orderStore a matrix one column after another. Moving down a column moves to the next address. Also called Fortran order..

NumPy generalises this with stridesstridesFor each dimension, how many bytes to jump to move one step along it. An (R, C) float64 array in row order has strides (8C, 8). The address of a[i, j] is base + istride0 + jstride1.. Every array is just a pointer, a shape, and a strides tuple. Slicing and transposing don’t copy any numbers; they make a new small array object with different strides pointing into the same buffer.

For a small matrix this doesn’t matter: the whole thing fits in cache, and every walk is a hit. For a big one it decides everything. On the M4 Pro, summing a 4096 × 4096 matrix of doubles (134 MB) with a two-line loop compiled by Numba (you’ll write it in the Build; a C version measured within a few percent of these):

Same additions, same numbers, same compiler, only the loop order swapped. And look at the last two rows: the compiler setting made the good order 5x faster and did nothing for the bad one. That’s the third trick, coming up.

NumPy is usually smart about this, but not always. For a reduction like M.T.sum(axis=1), NumPy’s iterator notices the strides and walks memory in its natural order anyway (I measured 1.36 ms vs 1.19 ms for the un-transposed version: nearly the same). But some operations must read in one order and write in another. M.T.copy() has to produce a new row-major array from a transposed view, so either its reads or its writes go down columns: 5.4 ms vs 1.2 ms for M.copy(). Every time you call np.ascontiguousarray, .copy() on a transpose, or hand a transposed array to a library that needs C order, you pay that.

You can see strides directly:

>>> a = np.zeros((1000, 10000))
>>> a.strides, a.T.strides, a[:, ::2].strides
((80000, 8), (8, 80000), (80000, 16))
>>> a.flags['C_CONTIGUOUS'], a.T.flags['C_CONTIGUOUS']
(True, False)

SIMD: doing several adds at once

The memory bill explains why a column walk is slow. It doesn’t explain why NumPy, in order, gets 0.12 ns per element while a plain compiled loop (Numba or C, they land within a few percent) over the same in-order data gets about 0.5. Both read memory perfectly. The difference is in how the arithmetic is done, and it has two parts.

Part 1: the dependency chain. The plain loop is s += a[i]. Each add needs the result of the previous add. A floating-point add takes a few clock cycles to produce its answer, so the core sits waiting for s before it can start the next one, even though it has spare adders. One cook, one pot, and every ingredient has to wait for the last one to dissolve. That’s where 0.5 ns per element comes from: it’s the latency of one add, not memory.

Fix: keep four running sums, s0 += a[i], s1 += a[i+1], s2 += a[i+2], s3 += a[i+3], and add them up at the end. Four independent chains, and the core runs them in parallel on its several adders. I measured 0.13-0.15 ns per element against 0.46-0.48 for one accumulator: about 3.4x faster from rearranging the arithmetic, no special instructions needed. This is instruction-level parallelisminstruction-level parallelismA single core executing several independent instructions at the same time. Modern cores have several adders, multipliers and load units and will overlap any instructions that don’t depend on each other..

Part 2: wide registers. Modern cores also have SIMDSIMDSingle Instruction, Multiple Data. One instruction operates on a short vector of numbers packed into one wide register, e.g. adds 2 doubles or 4 floats at once on 128-bit NEON. instructions. A vector registervector registerA wide CPU register holding several numbers side by side. 128 bits on ARM NEON (Apple M and the Spark’s cores), 256 bits with x86 AVX2, 512 bits with AVX-512. holds several numbers side by side, and one instruction adds all the lanes at once. ARM’s NEON registers are 128 bits: 2 doubles or 4 floats per instruction. Combine wide registers with several accumulators and you get the 0.07 ns per element I measured for an in-cache sum.

Why doesn’t the compiler do this for you? Because floating-point addition isn’t associative. (a + b) + c and a + (b + c) can differ in the last bits, since each add rounds. Splitting one sum into four partial sums changes the order of additions, so it can change the answer, and compilers aren’t allowed to change your answer without permission. Numba uses the same optimiser as the clang C compiler (it is called LLVM), and it follows the same rule. @njit(fastmath=True) is that permission, the same thing C programmers ask for with -ffast-math (it grants several other, riskier liberties too; don’t use it blindly in real code). For integers, where addition is exact, the compiler vectorizes the plain loop by itself: a one-accumulator int64 sum in Numba ran at 0.06 ns per element with no flag at all.

NumPy gets permission a cleaner way: its sum for floats is documented to use pairwise summation, adding numbers in blocks and then adding the block sums. That’s both more accurate than one long running sum and naturally split into independent chains the CPU can vectorize. So np.sum gets the memory bill right (one dense run) and the arithmetic bill right (many independent lanes).

SIMD helps until memory is the wall. Look at the numbers when the array is 800 MB instead of fitting in L1: the vectorized sum drops from 0.07 to 0.09 ns per element, about 89 GB/s, because now it’s waiting on RAM, not the adders. The one-accumulator sum barely changes (0.46 → 0.48), because it was never waiting on memory. That’s the roofline from Why ~17 tokens per second? in miniature: a slow kernel is compute-bound, a fast kernel on big data hits the memory roof.

And notice what SIMD needs: the numbers packed side by side in memory, so one load fills one vector register. A Python list can never feed it, because the doubles are inside separate objects. A strided array feeds it badly. Layout isn’t just about cache lines; it’s what allows the third trick at all.

Five ways, one chart

Now put all of it together.

How to read it, one gap at a time:

• list loop → sum(list) (10.7 → 1.9): remove the interpreter, keep the layout. About 5.6x. This is the interpreter bill.
• sum(list) → contiguous NumPy (1.9 → 0.12): remove the pointers and objects, so each number is 8 bytes instead of 40, and the C loop can use several accumulators and SIMD. About 16x. This is the layout bill.
• list loop → shuffled list loop (10.7 → 45.7) and sum(list) → shuffled (1.9 → 11.2): same code, scattered objects. The memory bill, with nothing else changed.
• contiguous → stride 64 B → stride 128 B (0.12 → 0.70 → 2.4): NumPy, same count of numbers, more lines per number. At stride 64 B on the Mac, two numbers still share each 128 B line; at stride 128 B, every number costs its own line. The prefetcher still predicts the stride, so this is bandwidth wasted, not latency.
• random gather (4.2): about one line per number with no prefetching, saved from disaster only by memory-level parallelism. (It also allocates an 80 MB temporary for a[perm], part of the cost.)

The same story on a GPU

Everything above gets louder on a GPU, because a GPU is a machine designed to be fed by bandwidth.

An NVIDIA GPU runs threads in groups of 32 called a warpwarpA group of 32 GPU threads that execute the same instruction at the same time, each on its own data. The basic unit the GPU schedules.. When a warp executes a load, all 32 threads issue their address at the same moment, and the memory system looks at the 32 addresses together. It moves data in 32-byte sectorsectorThe 32-byte unit an NVIDIA GPU’s caches and memory move. Four sectors make one 128-byte cache line.s. If the 32 addresses are neighbours, the warp’s whole request is served with a handful of sectors. That’s coalesced accesscoalesced accessWhen the threads of a warp read or write consecutive addresses, so the hardware can serve all of them with the minimum number of wide memory transactions..

GPU memory is fast (the Spark’s spec sheet says ~273 GB/s), but inference kernels are already memory-bound, as you saw in Why ~17 tokens per second?. An uncoalesced kernel that moves 8x the bytes runs close to 8x slower. That’s why so much GPU kernel engineering is about layout: tiling matrices into blocks that warps read row-wise, staging them in fast on-chip shared memoryshared memory (GPU)A small scratchpad on each GPU multiprocessor (tens to a couple hundred KB) that a block of threads manages by hand. Used to load a tile once, coalesced, then reuse it many times. Not the same as the OS shared memory between processes., and reusing each tile many times before going back to main memory. FlashAttention, at its heart, is temporal locality done by hand.

Where you meet this in PyTorch: .contiguous(). A PyTorch tensor is a pointer, a shape and strides, exactly like a NumPy array. x.T, x.transpose(1, 2), x.permute(0, 2, 1, 3) and many slices just swap strides. They’re free, and the result is non-contiguous. Then one of three things happens:

1. An op refuses it. .view() needs strides it can reinterpret without moving data, and fails on a transpose: view size is not compatible with input tensor's size and stride ... Use .reshape(...) instead. .reshape() quietly copies when it has to.
2. An op handles it, slowly. Elementwise ops on one tensor keep its layout and run at full speed. Mixing layouts is different: on my Mac’s CPU, x + y on two 4096 × 4096 float32 tensors took 1.1 ms, and x.T + y took 11.7 ms. One of the two inputs has to be read down columns.
3. An op copies it for you. Many kernels (custom CUDA ops, some attention implementations, anything that wants a raw pointer) call .contiguous() internally, which allocates a new tensor and does the transposing copy, often without you noticing. That copy is the M.T.copy() cost from earlier, on every call.

So when you see q = q.transpose(1, 2).contiguous() in an attention implementation, it’s someone paying the copy once, on purpose, so that the next big kernel gets coalesced reads. .contiguous() on a tensor that already is contiguous returns it unchanged, for free. And the CNNs you’ve trained have a version of this too: memory_format=torch.channels_last stores images as N, H, W, C instead of N, C, H, W, so the channels of one pixel sit next to each other, which is what convolution kernels on modern GPUs want to read together.

Build It

Everything goes in ~/workbench/caches/. You’ll write three programs, all Python: the five-ways sum in plain Python and NumPy, then the row/column walk and the dependency-chain sum with Numba. Then you’ll repeat the Numba ones on the Spark with hardware counters.

Why Numba for the last two? A plain Python loop spends 10-50 ns of bookkeeping per step, and the effects you’re after are a fraction of a nanosecond to a few nanoseconds per element, so the interpreter would bury them. NumbaNumbaA library that compiles one Python function into machine code the first time you call it. It understands a subset of Python (numbers, loops, NumPy arrays), and inside that subset loops run at roughly C speed. fixes that: put @njit above a function, and the first time you call it Numba works out the argument types and turns the function into the same kind of machine code a C compiler would produce. Every later call runs that machine code directly, with no interpreter in the way. If you met it in How slow is slow?, it’s the same tool.

mkdir -p ~/workbench/caches && cd ~/workbench/caches
python3 -m venv .venv && . .venv/bin/activate && pip install numpy numba
python -c "import numba; print(numba.__version__)"

# macOS
sysctl hw.cachelinesize hw.perflevel0.l1dcachesize hw.perflevel0.l2cachesize

# the Spark (Linux, ARM64)
lscpu -C
cat /sys/devices/system/cpu/cpu0/cache/index0/coherency_line_size

# five.py
import time, array
import numpy as np

N = 10_000_000

def bench(name, fn, reps=5):
    """Run fn() reps times; print best time in ms and ns per element, and the result."""
    ...

if __name__ == "__main__":
    rng = np.random.default_rng(0)
    a = rng.random(N)                       # contiguous float64, 80 MB
    lst = a.tolist()                        # 10M float objects, created in order
    arr = array.array("d", a.tobytes())     # packed doubles, no objects
    perm = rng.permutation(N)
    lst_shuf = ...                          # same values, objects visited in scattered order
    s16 = ...                               # a view with stride 16 elements (128 B), N elements long
    s8 = ...                                # stride 8 elements (64 B), N elements long
    # bench: python loop over lst, over lst_shuf, over arr
    #        sum(lst), sum(lst_shuf)
    #        a.sum(), s8.sum(), s16.sum(), a[perm].sum()

# rowcol.py   run: python rowcol.py 4096
import sys, time
import numpy as np
from numba import njit

@njit
def sum_rows(m):
    # TODO: sum all n*n elements of the 2-D array m, visiting row by row
    # (i outer, j inner, reading m[i, j]). Return the sum.
    ...

@njit
def sum_cols(m):
    # TODO: same sum, visiting column by column (j outer, i inner).
    ...

# Then make fastmath=True copies of both, e.g.
# sum_rows_fm = njit(fastmath=True)(sum_rows.py_func)
# (.py_func is the original Python function behind an @njit one.)

def best_ns_per_elem(f, m, reps=5):
    # TODO: call f(m) once untimed. The first call compiles the function, and
    # also makes the OS hand over the matrix's memory. Then time `reps` calls
    # with time.perf_counter_ns() and return (best time / m.size, the sum).
    ...

n = int(sys.argv[1]) if len(sys.argv) > 1 else 4096
m = (np.arange(n * n) % 7).astype(np.float64).reshape(n, n)   # row-major, like C
# TODO: print n, MB, ns/elem for rows and cols (plain and fastmath),
# the cols/rows ratio, and the sums (they must match).

@njit
def sum1(a): ...          # one accumulator: s += a[i]

@njit
def sum4(a): ...          # four accumulators s0..s3, step i by 4, combine at the end
                          # (don't forget the leftover elements when len(a) % 4 != 0)

sum1_fast = njit(fastmath=True)(sum1.py_func)   # sum1, compiler may reorder and vectorize

@njit
def repeat(f, a, reps):
    # TODO: call f(a) `reps` times and return something that depends on every result.
    ...

sudo apt install linux-tools-common linux-tools-$(uname -r)   # if perf is missing
perf list | grep -i -E "cache|l1d|ll_"                        # what this CPU offers
perf stat -e cache-references,cache-misses,L1-dcache-loads,L1-dcache-load-misses \
    python rowcol.py 4096 cols

Measure

Fill in your numbers. Best of several runs. ns per element.

Reference run on an M4 Pro (Python rows: CPython 3.14.6, NumPy 2.5.4; Numba rows: CPython 3.13.9, Numba 0.68, NumPy 2.5.4):

What to look for:

The two bills don’t simply add. Loop minus sum(list) is 8.8 ns per element for the in-order list (10.7 − 1.9), but 34.5 ns for the shuffled one (45.7 − 11.2). If the interpreter bill were a fixed tax, those would match. They don’t because layout and code interact: in the tight C loop inside sum, many object loads can be in flight at once, so the misses overlap. In the interpreter loop, each load is surrounded by so much other work that fewer overlap, and each miss is exposed more. “5x interpreter, 16x layout” is a good first story; this is the truer second one.

The row/column gap appears when the matrix stops fitting in cache. At n = 256 (512 KB) it was 1.1x for me, at 1024 (8 MB) 1.6x, at 4096 (134 MB) 5.8x, and 8192 (537 MB) stayed at 5.7x. Compare those sizes with your L2 and with Apple’s system-level cache (or the Spark’s L3). On the Spark, with 64 B lines, the column walk wastes 7/8 of each line instead of 15/16, so if anything expect a slightly smaller ratio than the Mac’s at the same size, but a slower row walk too. Measure it.

fastmath=True widens the gap. It only speeds up the row walk, because only the row walk was waiting on arithmetic. The column walk was waiting on memory all along, and no amount of SIMD helps a loop that’s waiting on RAM.

The SIMD sum hits the memory roof. In L1 the vectorized sum is 6.5x faster than sum1. From RAM it slows to about 89 GB/s on my Mac: roughly one core’s practical bandwidth. sum1 barely cares where the data is: 0.46 vs 0.48. Slow code never sees the memory wall; fast code always does.

Sorted gather is 8x faster than random gather (0.54 vs 4.2) with the same instructions, because sorted indices walk memory forward and the prefetcher can follow.

Break It

1. Make the matrix a power of two. Run rowcol.py with n = 4000, 4096 and 4100. On my Mac the column walk took 1.88, 2.78 and 1.81 ns per element. 4096 is about 50% worse than its neighbours, and the row walk didn’t move (0.48 for all three). Each column step at n = 4096 is exactly 32 KB, and a cache decides which of its slots a line may go in from a few bits of the address. Addresses exactly 32 KB apart share those bits, so they all compete for the same few slots and keep evicting each other, while the rest of the cache sits empty. These are called conflict missesconflict missesCache misses caused not by the cache being full but by too many lines mapping to the same small group of slots (a “set”). Classic with power-of-two strides.. The classic fix is padding: allocate rows of 4096 + 8 and ignore the extra. This is why some libraries pad their leading dimension, and why FFT and matrix code benchmarks at 4096 sometimes look weirdly bad.

2. Let the compiler cheat. Change repeat in chain.py so it calls f(a) but throws the answer away (return 0.0 at the end). I got 0.0000045 ns per element for the in-cache sum1, which works out to nearly two million GB/s: the compiler saw a function with no side effects whose results nobody used, and removed the calls. Make sum1 itself return nothing and the loop inside it vanishes too. Every microbenchmark has to make the compiler believe the work matters, and Numba’s compiler is as ruthless about this as a C one.

3. Shrink the problem until NumPy loses. Time sum(l) and a.sum() for 10, 100 and 1,000 elements. On my Mac: 59 ns vs 237 ns at 10 elements, about even at 100, and NumPy 5x ahead at 1,000. Every NumPy call pays a fixed few hundred nanoseconds (argument parsing, dtype dispatch, allocating the result) before it touches a single number. Then time np.array(l).sum(): it loses at every size, because converting a list to an array is itself a Python-speed loop over objects. NumPy pays off when data lives in arrays and many operations happen per conversion.

4. Defeat the prefetcher with a stride it can’t track. In five.py, sum a[::k] for k = 1, 2, 4, 8, 16, 32, 64, 512 and plot ns per element. You should see the cost rise until every element is its own line (k = 16 on the Mac, 8 on the Spark), then mostly flatten, since one line per element is one line per element. For large k, each step also crosses a 4 KB or 16 KB page, the prefetcher (on many chips) won’t cross pages, and the address translations start missing in the TLBTLBTranslation lookaside buffer: a small cache of recent virtual-to-physical address translations. Big strides touch a new page per element and miss here too.. Look for a second rise.

Explain It Back

Answer out loud, in plain words, as if an interviewer asked. Then check.

Go Deeper

• Ulrich Drepper, “What Every Programmer Should Know About Memory” (2007). The long, definitive tour of caches, lines, prefetching, TLBs and how to write code for them. Sections 3 and 6 are the core of this lesson, with real measurements.
• Randal Bryant and David O’Hallaron, Computer Systems: A Programmer’s Perspective, Chapter 6, “The Memory Hierarchy”. The textbook version, including the “memory mountain” plot (bandwidth vs stride vs size) that your Break It 4 is a slice of. The best structured explanation of locality there is.
• Igor Ostrovsky, “Gallery of Processor Cache Effects” (blog post, 2010). Seven tiny C# experiments, each showing one cache effect: line size, L1/L2 sizes, instruction-level parallelism, associativity, false sharing. Fast to read, and every one is worth reproducing.
• NumPy documentation, “The N-dimensional array (ndarray)” and “Internal memory layout of an ndarray” (numpy.org). Exactly how strides, views, C and Fortran order work. Plus the numpy.sum docs’ note on pairwise summation.
• NVIDIA, CUDA C++ Best Practices Guide, section “Coalesced Access to Global Memory”, and the CUDA C++ Programming Guide on device memory accesses. The official account of warps, sectors and what counts as coalesced, with stride diagrams.
• Tri Dao et al., “FlashAttention: Fast and Memory-Efficient Exact Attention with IO-Awareness” (2022). The ML paper where “count the bytes, tile for locality” turned a standard operation into a 2-4x faster one. Read it after this lesson and Why ~17 tokens per second?; it will read like an application of both.
• Maurice Wilkes, “Slave Memories and Dynamic Storage Allocation” (IEEE Transactions on Electronic Computers, 1965). Two pages. The original cache idea, stated so clearly it’s a little humbling.

Connections

This lesson stands on How slow is slow?, which gave you the rungs and the cache line, and Where do my variables live?, which gave you the list-of-pointers picture. Now you know what those pointers cost.

It hands the interpreter bill to Why is Python slow?, which opens up the bytecode loop that cost ~10 ns per element here. It hands the memory bill to Why ~17 tokens per second?: once your code uses every byte of every line, the only thing left is how many bytes per second the memory can deliver, and on the Spark that’s the number that sets your tokens per second. And the page-sized steps you saw in Break It 4 are the visible edge of Virtual memory and the OOM killer.