Multi-dimensional Arrays

Memory is one-dimensional. Problems are not. The whole story of multi-dimensional arrays — the entire reason NumPy, PyTorch, BLAS, OpenCV, every game engine, and every relational database had to make hard, opinionated decisions about layout — is the friction between those two facts. Section 1 showed that an array is a contiguous run of bytes plus an address formula. This section generalizes that formula to any number of dimensions, exposes the single concept (strides) that unifies all layouts, and shows why the same algorithm can be 10× faster or slower depending on which axis you put on the inside of the loop.

What a Multi-Dimensional Array Really Is

Mathematically, an $N$-dimensional array (a tensor) of shape $(n_1, n_2, \dots, n_N)$ is a function from a Cartesian product of index ranges into an element type:

$A : \prod_{d=1}^{N} \{0, 1, \dots, n_d - 1\} \to T.$

For $N=1$ this is a vector; for $N=2$ a matrix; for $N=3$ a cube of values; for $N=4$ a stack of cubes (mini-batches in deep learning); and so on. The total element count is $\prod_{d=1}^{N} n_d$.

That is the logical picture. Physically, the elements have to be ordered into a 1D byte strip, because that is the only thing memory understands. A multi-dimensional array is therefore two things glued together:

1. A contiguous (or near-contiguous) buffer of $\prod_d n_d$ elements in memory.
2. A layout function that maps an index tuple $(i_1, \dots, i_N)$ to a flat offset into that buffer.

Different libraries pick different layout functions. Picking one is the entire content of the next two sections.

Two Ways to Lay Out 2D in 1D Memory

For a 2D array $A$ of shape $(n_{\text{rows}}, n_{\text{cols}})$, there are two canonical orderings.

Row-major (C, Python lists, NumPy default, PyTorch default, OpenCV, Pillow)

Lay out row 0 entirely, then row 1, then row 2, and so on. Within a row, columns are contiguous. The flat offset (in elements) is

$\text{off}(i, j) = i \cdot n_{\text{cols}} + j.$

Column-major (Fortran, MATLAB, R, Julia, BLAS, LAPACK, Eigen)

Lay out column 0 entirely, then column 1, and so on. Within a column, rows are contiguous. The flat offset is

$\text{off}(i, j) = j \cdot n_{\text{rows}} + i.$

Same logical matrix, same element values, completely different physical addresses. The choice is not arbitrary: it leaks into the API of every numerical library you will ever use, and it determines whether a given loop is fast or slow.

The General N-D Formula and Strides

Both layouts are special cases of a single formula. Define the stride $s_d$ of axis $d$ as “how many elements to step in the buffer when index $i_d$ increases by 1, with all other indices held fixed”. Then

$\text{off}(i_1, i_2, \dots, i_N) = \sum_{d=1}^{N} i_d \cdot s_d.$

For row-major (C-order) on a tensor of shape $(n_1, \dots, n_N)$, the strides are determined by the recursion: the last axis varies fastest, so its stride is 1; each earlier axis strides over a full slab of all later axes. Concretely,

$s_N = 1, \quad s_d = s_{d+1} \cdot n_{d+1} = \prod_{k=d+1}^{N} n_k.$

For a $(n_0, n_1, n_2)$ 3D tensor in row-major order, the strides are $(n_1 n_2, \; n_2, \; 1)$. Walking $i_2$ by 1 advances by 1 element. Walking $i_1$ by 1 advances by $n_2$ elements (a full row). Walking $i_0$ by 1 advances by $n_1 n_2$ elements (a full slab). For column-major (F-order), the formula mirrors with the first axis varying fastest:

$s_1 = 1, \quad s_d = \prod_{k=1}^{d-1} n_k.$

True 2D Arrays vs Jagged Arrays of Pointers

Not everything advertised as a “2D array” is laid out as one contiguous block. There is a critical distinction between a true 2D array — where all $m \cdot n$ elements live in one contiguous run — and a jagged array (or array of arrays), which is a 1D array of pointers, each pointing to an independently allocated row.

The two-load form is the killer. When A[i][j] requires loading A[i] first to discover the row pointer and then dereferencing it, the second load cannot start until the first finishes (a data-dependent load), the prefetcher has nothing to predict (the row pointers can point anywhere), and each row may live in a different page of the heap. For a 1024×1024 matrix, the difference between a contiguous float A[1024][1024] and a float **A with 1024 separate mallocs is routinely 3–10× on a tight numerical loop.

The Loop-Order Rule and Cache Reality

The address formula tells you indexing is $O(1)$. The cache hierarchy tells you that the order in which you ask for those $O(1)$ elements determines whether you run at L1 speed or DRAM speed — a roughly 100× gap. The rule is simple and exception-free:

Iterate the dimension whose stride is 1 in the innermost loop. On row-major data, that is the last index. On column-major data, the first.

Why? Because a stride-1 access pattern reads consecutive bytes. A 64-byte cache line then carries 16 float32s or 8 float64s of useful work for the price of one DRAM round trip, and the hardware prefetcher streams the next line into L1 before you ask for it. A stride-$n_{\text{cols}}$ access pattern, by contrast, jumps a full row each iteration, touches a fresh cache line every time, and gives the prefetcher nothing to grab.

1# Row-major data (NumPy default). Inner loop varies LAST index → stride-1 → fast.
2for i in range(m):
3    for j in range(n):
4        out += A[i, j]                # consecutive 8-byte reads
5
6# Same logical work, but inner loop varies FIRST index → stride-n → cache miss every step.
7for j in range(n):
8    for i in range(m):
9        out += A[i, j]                # 8*n-byte jumps; prefetcher gives up

On a 4096×4096 float64 matrix (128 MB — well outside any L2 or L3) the second loop is typically 5–15× slower than the first on the same CPU, with the same FLOP count, the same arithmetic operations, and the same correctness. Nothing changed except the order in which we asked memory for our own data.

Interactive: ND-Tensor Stride Inspector

The widget below makes the abstract formulas concrete. Pick a 2D or 3D shape, drag the index sliders, and watch three things at once: (1) the logical grid view, (2) the physical 1D strip in memory order, and (3) the offset formula expanded with your specific numbers. In 2D mode, toggle row-major / column-major — the data does not move, but the strip's ordering does. In 3D mode, flip the transpose switch to swap the strides of axis 0 and axis 2; the same buffer now appears in a different visual order, which is exactly how a zero-copy transpose works in PyTorch.

Two patterns to look for. First, when you change the layout in 2D mode, observe that the same logical cell [i, j] lands at completely different positions in the strip. Second, the transpose toggle in 3D mode flips strides without ever rewriting the buffer; that is why x.transpose(0, 2), x.permute(...), and x[::-1] all run in $O(1)$ regardless of how big the tensor is.

Case Study: Matrix Multiplication and Loop Order

Matrix multiplication $C = AB$ with $A \in \mathbb{R}^{m \times k}$, $B \in \mathbb{R}^{k \times n}$, $C \in \mathbb{R}^{m \times n}$ is the textbook playground for loop-order effects. The naive triple loop has six possible orderings of the three indices (i, j, k); they all produce identical answers but exhibit dramatically different cache behavior on row-major storage.

On a 1024×1024 float32 multiply, the i k j order routinely beats i j k by 3–5× in plain C with optimization on, simply because the innermost loop now streams contiguously through both $B$ and $C$. This is also why production BLAS (OpenBLAS, MKL, BLIS) does not use any of these orderings directly. They use tiled or blocked matrix multiply: chop $A$, $B$, and $C$ into small panels that fit in L1 / L2, multiply panel-by-panel, and reuse each panel many times before evicting it. The result is the difference between 2 GFLOPS (naive) and 50–200 GFLOPS (BLAS) on the same CPU.

Higher Dimensions: Images, Video, Mini-Batches

Three- and four-dimensional arrays are the workhorses of perception, graphics, and deep learning. The shape conventions are not historical accidents — they are deliberate choices about which axis is stride-1, because that determines which kernels run fast.

Code: A Hand-Built Array2D in Python

To make the formulas concrete, here is a 2D array implemented from scratch on top of a flat Python list. Notice how every read and write reduces to the row-major offset formula.

Array2D(4, 5) → self.data is a list of 20 zeros

A[2, 3] → 2 * 5 + 3 = 13

self.data = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

After this line, position 13 of the flat buffer holds 7.

[0,0,0,0,0, 0,0,0,0,0, 0,0,0,7,0, 0,0,0,0,0]

13

1class Array2D:
2    def __init__(self, n_rows, n_cols, fill=0):
3        self.n_rows = n_rows
4        self.n_cols = n_cols
5        self.data = [fill] * (n_rows * n_cols)   # one flat 1-D buffer
6
7    def _flat(self, i, j):
8        if not (0 <= i < self.n_rows and 0 <= j < self.n_cols):
9            raise IndexError(f"({i}, {j}) out of bounds for ({self.n_rows}, {self.n_cols})")
10        return i * self.n_cols + j               # row-major offset
11
12    def __getitem__(self, key):
13        i, j = key
14        return self.data[self._flat(i, j)]
15
16    def __setitem__(self, key, value):
17        i, j = key
18        self.data[self._flat(i, j)] = value
19
20
21A = Array2D(4, 5)
22A[2, 3] = 7
23print(A.data)            # [0,0,0,0,0, 0,0,0,0,0, 0,0,0,7,0, 0,0,0,0,0]
24print(A._flat(2, 3))     # 13

Swapping the offset formula to j * self.n_rows + i would convert this class to column-major. The data buffer would not change shape; only the mapping function would. Layout is software, not data.

Code: Strides, Views, and Transpose-Without-Copy (NumPy / PyTorch)

Real libraries do not write a fresh class for every shape. They store one buffer plus a shape tuple and a strides tuple, and let the strides tuple do all the work. Here is how to read it out and watch a transpose happen with zero copy.

shape = (2, 3, 4), 96 bytes total, values 0..23

A.shape = (2, 3, 4)

A.strides = (48, 16, 4)

B.shape = (4, 3, 2)

B.strides = (4, 16, 48)

share memory? True

i=1, j=2, k=3

addr_A = base + 92

addr_B = base + 92  (matches addr_A)

values: 23 == 23

1import numpy as np
2
3A = np.arange(24, dtype=np.int32).reshape(2, 3, 4)
4print("A.shape  =", A.shape)
5print("A.strides=", A.strides)        # bytes per axis step
6print("A.data   =", hex(A.ctypes.data))
7
8B = A.transpose(2, 1, 0)              # zero-copy: just permutes shape and strides
9print("B.shape  =", B.shape)
10print("B.strides=", B.strides)
11print("share memory?", np.shares_memory(A, B))
12
13# Same element, two index tuples, same byte address
14i, j, k = 1, 2, 3
15addr_A = A.ctypes.data + i*A.strides[0] + j*A.strides[1] + k*A.strides[2]
16addr_B = B.ctypes.data + k*B.strides[0] + j*B.strides[1] + i*B.strides[2]
17print("A[1,2,3] addr:", hex(addr_A))
18print("B[3,2,1] addr:", hex(addr_B))
19print("values:", A[i, j, k], "==", B[k, j, i])

Identical idea in PyTorch: x.permute(2, 1, 0), x.transpose(0, 2), and x.t() all return non-contiguous views over the same storage. When a downstream op (like view or reshape) needs a contiguous buffer, you must call .contiguous() — which is the moment the bytes actually get copied. Until then, all the rearranging happens in stride metadata.

The same idea in C++ / Eigen

1#include <Eigen/Dense>
2using Eigen::MatrixXd;
3
4MatrixXd A(1000, 1000);
5A.setRandom();
6
7// Eigen is column-major by default. A.transpose() returns an expression
8// template that swaps strides — no allocation, no memcpy.
9auto At = A.transpose();              // O(1)
10
11double s = (A * At).trace();          // multiply uses BOTH layouts efficiently

Eigen, BLAS, and LAPACK all expose a “leading dimension” (lda) parameter on their matrix kernels. That parameter is the stride between consecutive columns (or rows). It exists precisely so that submatrix views — which have a different leading dimension than the parent matrix — can be passed in without copying.

Where This Shows Up in the Real World

Pitfalls and Gotchas

1. Mixing row-major and column-major libraries. Passing a NumPy row-major matrix into a Fortran-style BLAS routine without setting the “transposed” flag silently solves the transposed problem. Symptoms: numerically plausible but wrong results, or shape mismatches that only appear for non-square matrices.
2. “2D arrays” that are arrays of pointers. Java's int[][], C's int **A, Python's list-of-lists. The rows are not contiguous with each other and may not even be the same length. For numerical work, flatten to a single buffer.
3. Confusing logical shape with memory layout. A tensor of shape $(2, 3, 4)$ can be stored in 6 different stride orderings (3! axis permutations), all containing the same elements. Two tensors with the same shape can have wildly different cache behavior if their strides differ.
4. Wrong loop order on multidimensional data. The single most common cache-induced slowdown. On row-major data, the rightmost index varies in the innermost loop. On column-major data, the leftmost index does. Get it backwards and you pay 5–15× for the same computation.
5. Forgetting .contiguous() before view / reshape. A non-contiguous tensor (e.g., the output of transpose) cannot be reshaped without copying. PyTorch will throw “view size is not compatible with input tensor's size and stride”. Calling .contiguous() forces the copy explicitly so you know it is happening.
6. Strided views aliasing the original. x[::2] shares storage with x. Writing to it mutates x. If you need independence, call .copy() in NumPy or .clone() in PyTorch.
7. NCHW vs NHWC mismatches in deep learning pipelines. A model trained in PyTorch (NCHW) will silently transpose every input if the data loader hands it NHWC. The correctness is preserved by the transpose, but throughput collapses because every batch now allocates a new contiguous buffer.
8. Padding / row stride differences in image libraries. OpenCV often pads each row to a multiple of 4 or 16 bytes for SIMD alignment. The image's step (row stride) is therefore not always width * channels. Iterating (width * channels) instead of step reads garbage from the padding.

Hold on to one image. A multi-dimensional array is a single contiguous buffer plus a stride vector. Indexing computes $\sum_d i_d \cdot s_d$, one multiply-add per axis, and that is it. Row-major and column-major are just two specific stride choices. Transpose, slice, permute, broadcast, and reshape are stride manipulations — the data does not move. The cache hierarchy turns this picture into a performance contract: iterate the stride-1 axis on the inside, or pay 10×. Every fast numerical library — NumPy, PyTorch, BLAS, cuDNN, OpenCV, Parquet readers, game engines — is, at its core, a careful negotiation between this contract and the shape of the problem.