Skillversitas

Memory Layout and Performance

NumPy's speed comes from three places: contiguous typed memory the C engine can walk without per-element dispatch, operations expressed as whole-array loops in C, and a layout the CPU's prefetcher and SIMD units can exploit. Understanding where time actually goes - allocation, data movement, or compute - is what separates "it works" from "it is fast." This chapter is the performance model: how data is laid out in memory, which operations allocate, and how to measure so you optimize the right thing.

C-Order, F-Order, and Strides

import numpy as np

m = np.arange(12).reshape(3, 4)      # C-contiguous by default
print(m.flags["C_CONTIGUOUS"])       # True - row 0, then row 1, ...
print(m.strides)                      # (32, 8) - bytes to step rows, then cols

f = np.asfortranarray(m)             # Fortran order: columns contiguous
print(f.flags["F_CONTIGUOUS"])       # True
print(f.strides)                      # (8, 24) - the column step is large now

# a transpose is non-contiguous in both orders, yet still fast for most ops
t = m.T
print(t.flags["C_CONTIGUOUS"], t.flags["F_CONTIGUOUS"])  # False False

Two storage orders exist. C-order (the default) stores rows contiguously; F-order stores columns contiguously (the Fortran convention). The strides tuple is the low-level truth: the bytes to advance one index along each axis. Most NumPy operations handle either order fine, but routines that expect a specific layout - some C libraries, certain reductions - run faster when the data is already in the order they want. ascontiguousarray and asfortranarray coerce the order explicitly when that matters.

What Strides Actually Control

a = np.arange(6).reshape(2, 3)
print(a.strides)          # (24, 8) for int64: 3 elements * 8 bytes per row

# a slice changes the offset and the shape, never the data
row = a[1]
print(row.strides)         # (8,) - stepping one element is 8 bytes
print(row.data_ptr() - a.data_ptr())   # 24 - the byte offset into a

# a non-contiguous slice: the step skips bytes
sub = a[:, ::2]            # every other column
print(sub.strides)          # (24, 16) - 8 bytes skipped between columns
print(sub.flags["C_CONTIGUOUS"])  # False - the engine walks with gaps

Strides are the mechanism behind views: a strided slice keeps the same buffer but steps by more than itemsize, so the elements are not physically adjacent. The engine still walks them correctly - it just cannot use the tightest SIMD loops or hand the block to a C routine that assumes contiguity. That is the precise reason a [::2] slice is a view (strides changed) yet can be slower than a copied, compact array in a compute-heavy pipeline.

The Vectorization Rule

import time

n = 1_000_000
a = np.random.default_rng(0).normal(size=n)
b = np.random.default_rng(1).normal(size=n)

# the Python loop: per-element dispatch, the slow path
t0 = time.perf_counter()
py = [ai * bi for ai, bi in zip(a.tolist(), b.tolist())]
t_loop = time.perf_counter() - t0

# the ufunc: one C loop over contiguous memory
t0 = time.perf_counter()
vec = a * b
t_vec = time.perf_counter() - t0

print(f"loop {t_loop:.4f}s  vec {t_vec:.4f}s  speedup {t_loop/t_vec:.0f}x")

The speedup - typically 50 to 100x - is the entire argument for NumPy. The loop pays interpreter dispatch, pointer chasing, and no instruction-level parallelism for every element; the ufunc is one tight C loop the CPU can vectorize. The rule: if you can express a computation as a whole-array operation, do it that way. The Python loop is reserved for control flow the array engine cannot express - early exit, stateful logic, calls into non-NumPy code.

Allocation Is the Hidden Cost

a = np.random.default_rng(0).normal(size=100_000)
b = np.random.default_rng(1).normal(size=100_000)

# naive: each expression allocates a fresh temporary
c = a * b + a - b + 1.0            # several temporaries, several allocations

# zero-allocation form: reuse one buffer, write in place
out = np.empty_like(a)
np.multiply(a, b, out=out)
np.add(out, a, out=out)
np.subtract(out, b, out=out)
out += 1.0                          # in-place, no new buffer

In a hot loop, the temporary that every a * b + c allocates is often the cost - not the arithmetic. The out= parameter and augmented assignment (*=, +=) write into a named buffer, removing the allocation. The pattern: preallocate the result once, then chain out= calls. That is the difference between a pipeline that scales and one that thrashes the allocator.

einsum for Contracted Products

X = np.random.default_rng(0).normal(size=(1000, 50, 8))
W = np.random.default_rng(1).normal(size=(1000, 8, 12))

# the explicit matmul path
r1 = X @ W                                   # (1000, 50, 12)

# the same contraction as an einsum string - explicit about the axes
r2 = np.einsum("bai,bi->bao", X, W)          # identical result
print(np.array_equal(r1, r2))

# einsum also expresses things @ cannot, e.g. a per-batch trace
T = np.einsum("bii->b", X)                   # trace over the last two axes

einsum names the axes of every operand and the output with a string, making the contraction explicit and often faster than an equivalent reshape-plus-matmul chain (it can fuse operations and skip intermediates). Use it when the axis bookkeeping is non-trivial - batched contractions, diagonal extraction, outer products - where the string reads as the math. For the common A @ B, plain @ is clearer and equally fast.

dtype and Memory Pressure

  • float32 halves memory and can halve runtime on compute-bound, precision-tolerant work - a deliberate trade to make once, with a test, not by default.
  • int64 doubles the footprint of an index array versus int32; for large index sets, intp (or int32 when the size fits) is the leaner choice.
  • A boolean mask is one byte per element - a mask over a billion-element array is a gigabyte. Build masks in chunks when the full mask will not fit comfortably.
  • Prefer empty over zeros when you will write every element before reading it - you skip the fill loop entirely.

Measuring Before Optimizing

import time
import numpy as np

a = np.random.default_rng(0).normal(size=(2000, 2000))

def timeit_block(fn, reps=5):
    best = float("inf")
    for _ in range(reps):
        t0 = time.perf_counter()
        fn()
        best = min(best, time.perf_counter() - t0)
    return best

t_matmul = timeit_block(lambda: a @ a)
t_trans  = timeit_block(lambda: a.T)        # a view - near-zero work
print(f"matmul {t_matmul:.4f}s   transpose {t_trans:.6f}s")

The discipline: time the real expression, at the real shape and dtype, taking the best of a few reps - and only then decide where the cost is. The usual findings, in order of payoff: a Python loop that should be a ufunc; an allocation in a hot loop that should use out=; a float64 that could be float32; a non-contiguous slice feeding a layout-sensitive routine. Everything else is micro-optimization - measure it, because the "obviously slow" line is rarely the real one.

Summary

NumPy's speed is contiguous typed memory plus C loops plus a layout the CPU can exploit. The performance model: understand strides and C/F order (coercing with ascontiguousarray/asfortranarray when a routine wants a specific layout), express computations as whole-array ufuncs instead of Python loops, eliminate hot-loop allocations with out= and augmented assignment, use einsum for non-trivial contractions, and make the float32 and index-width choices deliberately. Above all, measure the real expression at the real shape before changing anything - the allocation or the loop is almost always the cost, not the arithmetic.