Memory Layout and Strides

Core Principle: Multidimensional indexing translates coordinate tuples into 1D buffer byte offsets via strides, governing memory layout efficiency, transposition costs, and CPU cache performance.

The mechanism

(32, 8)

To move one step along axis 1 (across a row), advance 8 bytes: one int64. To move one step along axis 0 (down a column), advance 32 bytes: four elements, i.e. a whole row. The address of a[i, j] is base + 32*i + 8*j.

That is the entire model. Every "free" operation in NumPy is a strides edit.

C order vs Fortran order

array([[ 0,  1,  2,  3],
       [ 4,  5,  6,  7],
       [ 8,  9, 10, 11]])

NumPy defaults to C order (row-major), like C and Python. R, MATLAB, Fortran and BLAS internals use column-major. When you receive a matrix from R via rpy2 or read a MATLAB .mat file, expect Fortran order.

Transpose is free: and that is the whole trick

False

a.T swaps the strides tuple. Nothing is copied, at any size. But the result is no longer C-contiguous, and that has consequences downstream.

Why non-contiguity costs you

Modern CPUs read memory in cache lines (~64 bytes). Walking memory sequentially means each cache line fetch serves ~8 float64 values. Walking with a large stride means each fetch serves one value and the rest of the line is discarded.

array([0., 0., 0., ..., 0., 0., 0.], shape=(10000,))

Same number of arithmetic operations, several times the runtime, entirely because of memory access order. For a 20,000 × 500 expression matrix this is usually irrelevant; for a 50,000 × 50,000 matrix it dominates.

The practical rule: if you will repeatedly reduce along a particular axis, store the array so that axis is contiguous.

expr_by_gene = np.ascontiguousarray(expr)        # fast expr.mean(axis=1)
expr_by_sample = np.asfortranarray(expr)         # fast expr.mean(axis=0)

ascontiguousarray costs one full copy, so it pays off only if you reduce many times.

Broadcasting is a zero stride

(8, 0)

A stride of 0 means "do not move": every step along that axis re-reads the same element. That is how Broadcasting avoids allocating. The result is read-only, because writing through a zero stride would be ambiguous.

Sliding windows for free

(8, 8)

An (8, 3) array occupying 10 elements of memory. This makes rolling computations essentially free until you reduce over them:

gc_100bp = sliding_window_view(is_gc, 100).mean(axis=1)
rolling_depth = sliding_window_view(coverage, 50).mean(axis=1)

The older as_strided does the same thing with no safety checks and will happily read past the end of your buffer, segfaulting or returning garbage. Use sliding_window_view; avoid as_strided.

Diagnostics

a.flags                     # full contiguity/ownership report
a.strides
a.base                      # what it views, if anything
a.nbytes                    # logical size: misleading for broadcast views
a.base.nbytes               # actual memory
np.shares_memory(a, b)

When this actually matters in bioinformatics

• Large matrices (single-cell: 100k cells × 30k genes; whole-genome coverage arrays). Axis order can be a 5× runtime difference.
• Passing arrays to compiled libraries. BLAS/LAPACK, and many C extensions, require contiguous input and will copy silently if you give them a transposed view. If a linear algebra call is unexpectedly slow, check contiguity.
• Memory-mapped files. np.memmap lets you slice a file larger than RAM. Only the touched pages load, and access order determines whether that is fast or thrashing.

  arr = np.memmap("big.dat", dtype=np.float32, mode="r", shape=(1_000_000, 100))
  chunk = arr[0:1000]         # only these pages are read
  

• Interop with R/MATLAB output: Fortran order.

Common mistakes

• Assuming .T is expensive. It is free. What may be expensive is the operation after it.
• Assuming a view is contiguous. a[::2] is a view with stride 16, not contiguous.
• Using as_strided casually. It is a loaded gun.
• Reducing along the wrong axis repeatedly on a very large array and blaming NumPy.
• Calling np.ascontiguousarray in a hot loop: it copies every time.
• Confusing nbytes with real memory on a broadcast view.

Going further

This is the boundary where NumPy stops and systems programming begins. If you find yourself here regularly, the next steps are numba (JIT with explicit control over loop order), zarr/dask (chunked out-of-core arrays), and scipy.sparse (for the >90%-zero matrices that single-cell genomics produces).

See also

ndarray · Views vs Copies · Broadcasting · Reshaping and Stacking · Vectorization · NumPy File IO

Exercise

Write layout_report(matrix, step), which compares a C-ordered integer matrix with two arrays derived from it. Return ["bad step"] when step is below one.

Otherwise return five lines: c=<bool> f=<bool> for the matrix's own contiguity flags; t_c=<bool> t_f=<bool> for its transpose, which reuses the same buffer with the strides swapped; sub_c=<bool> ratio=<int> for the strided column selection matrix[:, ::step], where the ratio is that view's column stride divided by the matrix's own column stride; ravel_shares=<bool>, whether sub.ravel() still aliases the matrix, which it can only do when the selection is contiguous; and itemsize=<int> row_stride_elems=<int>, the row stride expressed in elements rather than bytes. A matrix with a single column is both C- and F-contiguous at once, so do not assume the two flags are opposites.