pythonforbio.
[WASM idle]
NumPy08/18

Memory Layout and Strides

Starting sandbox…
Beginnerlesson2 tests

Memory Layout and Strides

In one line: an array is a flat buffer plus a rule for converting an index tuple into a byte offset — that rule is the strides, and it explains transposes, views and mysterious slowdowns.

The mechanism

a = np.arange(12, dtype=np.int64).reshape(3, 4)
a.strides        # (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

a.flags.c_contiguous       # True: last axis is contiguous in memory (row-major)
a.flags.f_contiguous       # True: first axis is contiguous (column-major)
np.ascontiguousarray(a)    # force C order
np.asfortranarray(a)       # force F order

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

big = np.zeros((100_000, 1000))     # 800 MB
big.T                                # instant
big.T.strides                        # strides swapped; not one byte moved
big.T.flags.c_contiguous             # 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.

a = np.zeros((10_000, 10_000))
a.sum(axis=1)      # walks along rows — contiguous — fast
a.sum(axis=0)      # walks down columns — strided — several times slower

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

v = np.arange(3)[:, None]              # (3, 1)
b = np.broadcast_to(v, (3, 1000))
b.strides                               # (8, 0)  ← zero!

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

from numpy.lib.stride_tricks import sliding_window_view

x = np.arange(10)
w = sliding_window_view(x, 3)     # shape (8, 3) — a VIEW
w.strides                          # (8, 8) — overlapping windows, one buffer

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

Test cases · 2

#viainputexpected stdout
1stdin1 2 3 4 5 6c_contig=True t_c_contig=False t_f_contig=True
2stdin1 2 3 4c_contig=True t_c_contig=False t_f_contig=True

Hints · 2

01Hint

Transposing does not move data — it swaps the strides.

02Hint

arr.flags['C_CONTIGUOUS'] tells you whether rows are contiguous.

starter

No output yet — run the code to populate this drawer.