Memory layout
A two-dimensional array is stored as one flat run of memory. The layout decides which neighbours are actually adjacent, and that decides how well the CPU cache helps you.
Two orders
In row-major order, used by C and by NumPy by default, each row is contiguous. In column-major order, used by Fortran, MATLAB and Julia, each column is contiguous. For an array with m rows and n columns:
| Order | Offset of element (i, j) |
|---|---|
| Row-major | i * n + j |
| Column-major | i + j * m |
Strides
NumPy generalises this with strides: the number of bytes to move for one step along each axis.
import numpy as np
a = np.zeros((3, 4), dtype=np.float64)
a.strides # (32, 8)
np.asfortranarray(a).strides # (8, 24)
a.T.strides # (8, 32)The last line shows why transpose is cheap. It copies nothing and only swaps the strides, so a.T is a view that happens to be laid out column-first.
Why loop order matters
Caches fetch whole lines of memory, commonly 64 bytes, so reading neighbouring addresses is far cheaper than jumping around. Summing a row-major matrix with the outer loop over rows walks memory in a straight line. Swapping the loops strides across a full row on every step, and on a matrix too large for the cache that can be several times slower while doing identical arithmetic.
The same effect appears at a higher level in NumPy. Reducing along the fast axis touches consecutive values, so when an array is built once and reduced many times, it pays to store it in the order that matches how it is read.
Things to watch
A view such as a[:, 0] is not contiguous, because its stride is a full row. Passing it to native code usually forces a copy. Checking a.flags["C_CONTIGUOUS"] shows whether that will happen, and reshaping a non-contiguous view may also copy rather than return a view.