How Memory Holds a Matrix

The idea in one sentence

Memory is one long line of numbered bytes, so a matrix is a rule for turning a row and a column into an address, and everything fast or slow about matrix code starts with that rule.

The line under the grid

You draw a matrix as a grid because a grid is how the algebra reads. Rows across, columns down, an entry at every crossing. The machine has none of that. It has a line of bytes, numbered from zero, and one instruction that says fetch the byte at number $q$.

So the grid has to be a promise. Somebody has to decide, once, which byte holds the entry in row $i$ and column $j$, and then keep the promise forever. That decision is called the layout, and it is the whole subject of this chapter.

The simplest promise is the one C makes. Lay row 0 down first, left to right. Then row 1 right after it. Then row 2. Nothing is skipped and nothing is repeated, so the entry in row $i$ and column $j$ is preceded by $i$ whole rows and $j$ more entries.

Write $\alpha(i, j)$ for the address of that entry, $p$ for the address of the first entry, $s$ for the size of one number in bytes, and $n$ for the number of columns. Then

Everything about the layout is in two numbers per axis. The shape says how many steps an index may take, and the stride says how many bytes one step costs. Moving one column to the right costs $s$ bytes. Moving one row down costs $n\,s$ bytes, because a row is $n$ numbers long.

Give those two multipliers names and the formula stops depending on the convention. Let $s_0$ be the stride along rows and $s_1$ the stride along columns:

A $d$-dimensional array works the same way, with one stride per axis:

That sum is the entire theory. The rest of this chapter is what happens when you choose the $s_k$ badly.

Two conventions, and who chose which

C fills a row before moving down. Fortran fills a column before moving right. Both are consistent, both cost one multiplication and one addition, and the two of them split the numerical world in half. Fortran’s promise is

with $m$ the number of rows. The names for the two habits are row-major and column-major, and Wikipedia’s row- and column-major order article keeps the roll call.

The last row is the one that costs people time. NumPy is row-major, and every serious linear algebra library underneath it is column-major, because BLAS and LAPACK grew up in Fortran. That mismatch is handled for you, and it is handled by a flag rather than by copying anything, which is a trick this chapter earns two sections from now.

NumPy will build either layout on request, and it will tell you which one you have.

import numpy as np

A = np.arange(12.0).reshape(3, 4)          # row-major, the default
F = np.asfortranarray(A)                   # the same numbers, column-major

print(A.strides)                           # (32, 8)
print(F.strides)                           # (8, 24)
print(np.array_equal(A, F))                # True

print(A.ravel(order='K'))                  # [ 0.  1.  2.  3.  4.  5.  6.  7.  8.  9. 10. 11.]
print(F.ravel(order='K'))                  # [ 0.  4.  8.  1.  5.  9.  2.  6. 10.  3.  7. 11.]

Read A.strides as the two multipliers from Figure 2. One step down a column costs 32 bytes, because a row is 4 doubles long and you have to clear it. One step right costs 8, the size of a double. In F those numbers trade places: a step right costs 24 because a column is 3 doubles long, and a step down costs 8. The order='K' flag on ravel means walk in memory order, so the two printed lists are literally the two lines of bytes.

The matrices are equal. The bytes are in a different order. Nothing else about them differs, and everything about their speed does.

Strides make views free

Once the layout is two numbers per axis, a surprising amount of array manipulation costs nothing. It is arithmetic on the strides, and the data never moves.

Transposing swaps the two strides. Row 5 of $A^{\T}$ is column 5 of $A$, so the stride that used to carry you along a row now carries you down a column, and the promise still holds for every entry.

Slicing with a step multiplies a stride. Reshaping a contiguous array recomputes the strides from the new shape. Adding a length-one axis inserts a stride nobody will ever use. All of it is bookkeeping.

import numpy as np

A = np.arange(12.0).reshape(3, 4)
print(A.T.strides)                         # (8, 32)
print(np.shares_memory(A, A.T))            # True

B = np.arange(24.0).reshape(2, 3, 4)
print(B.strides)                           # (96, 32, 8)
print(B[:, ::2, ::2].strides)              # (96, 64, 16)
print(B.transpose(2, 0, 1).strides)        # (8, 96, 32)
print(B.reshape(6, 4).strides)             # (32, 8)

The three strides of B read from the outside in. One step along the first axis skips a whole $3 \times 4$ block, which is 96 bytes. One step along the second skips a row of 4 doubles, which is 32. One step along the last is a single double. Take every other element on the last two axes and those two strides double. Nothing was copied, and np.shares_memory says so.

This is also why the transpose flag in BLAS is not a courtesy. A row-major matrix, read with the strides swapped, is a column-major matrix. The Fortran library and the C caller can agree about the bytes and disagree about which index runs fastest, and neither of them has to move anything.

Sixteen numbers arrive whether you want them or not

So far every address costs the same. That is the part memory does not honor.

Between the core and main memory sits a cache, and the cache does not deal in bytes. It deals in blocks called cache lines. Ask for one number and the machine fetches the whole aligned block around it, because fetching 128 bytes costs almost exactly what fetching 8 would.

You can read the number off the machine you are sitting at. On the Apple silicon Mac this chapter was written on, sysctl hw.cachelinesize prints 128. On most x86 desktops it prints 64. The size changes; the fact does not.

Now put the cache line next to the address formula and the whole chapter falls out. Walking along a row of a row-major matrix moves 8 bytes at a time, so 16 consecutive entries ride in on one line and 15 of them were free. Walking down a column moves $n\,s$ bytes at a time. Once $n\,s$ reaches 128, every single entry needs a line of its own, and the other 15 numbers on each line are thrown away.

That is a sixteen-fold difference in traffic for the same count of arithmetic operations, and nothing in the source code announces it. The loop looks the same either way. Only the order of the two indices tells you which one you wrote.

What a stride costs, measured

The argument above is worth checking rather than believing, and checking it is easy. Take one long array, sum every element, then sum every second element, every fourth, and so on. Each run touches fewer elements, so measure the cost per element touched. If the cache line story is right, that cost should climb while the stride is still smaller than a line, then stop climbing the moment each element needs a line of its own.

import numpy as np
from timeit import repeat

a = np.ones(1 << 25, dtype=np.float64)     # 256 MiB, far larger than any cache

for step in (1, 2, 4, 8, 16, 32, 64):
    v = a[::step]
    t = min(repeat(lambda: v.sum(), number=3, repeat=9)) / 3
    print(step * 8, "bytes  %.3f ns/element" % (t / v.size * 1e9))

#    8 bytes  0.134 ns/element
#   16 bytes  0.212 ns/element
#   32 bytes  0.387 ns/element
#   64 bytes  0.708 ns/element
#  128 bytes  2.596 ns/element
#  256 bytes  2.682 ns/element
#  512 bytes  2.581 ns/element

The knee lands on 128 bytes, which is what sysctl reported, and the flat stretch after it is the proof: a stride of 512 bytes is no worse than 128, because both throw away everything except one number per line. Between the two ends is a factor of about 19.

Force it and the difference arrives. Copying a matrix into a fresh C-order buffer is a walk along rows; copying its transpose into the same kind of buffer is a walk down columns, and no reordering can rescue both sides at once.

import numpy as np
from timeit import repeat

M = np.arange(4096 * 4096, dtype=np.float64).reshape(4096, 4096)
bench = lambda f: min(repeat(f, number=5, repeat=7)) / 5

good = bench(lambda: M.copy())             # 0.002230 s
bad = bench(lambda: M.T.copy())            # 0.073013 s
print(bad / good)                          # 32.73782816383471

The same 128 MiB in, the same 128 MiB out, thirty times the time. Repeated runs here put the factor between 27 and 38, which is a wide spread and a small point: the exact number depends on what else the machine is doing, and the order of magnitude does not. The only difference between the two lines is the order of the visits.

The number BLAS asks for that nothing else does

One more consequence, and it explains a signature that puzzles everyone who meets it.

A library rarely gets a whole matrix. It gets a window inside a larger allocation: a block of a big matrix, a slice of a batch, a panel that some blocking scheme cut out. The window has its own height and width. It does not have its own stride. Its rows are still spaced by the width of the thing it was cut from.

import numpy as np

big = np.zeros((1000, 1000))
win = big[100:110, 200:206]

print(win.shape)                           # (10, 6)
print(win.strides)                         # (8000, 8)
print(win.flags['C_CONTIGUOUS'])           # False
print(np.shares_memory(big, win))          # True

The row stride is 8000 bytes, which is 1000 doubles, which is the parent’s width. That is why every BLAS routine carries a size argument you did not expect. Here is the workhorse, the one that occupies the rest of this part of the course:

void sgemm_(const char *transa, const char *transb,
            const int *m, const int *n, const int *k,
            const float *alpha,
            const float *a, const int *lda,
            const float *b, const int *ldb,
            const float *beta,
            float *c, const int *ldc);

The shapes are m, n and k. The three extra integers lda, ldb and ldc are the leading dimensions. The reference documentation defines lda as the first dimension of a as declared in the calling program, at least max(1, m) when transa is N. Shape and stride are asked for separately because they are separate facts, exactly as in Figure 2. It is also what lets one call work on a block of a much larger matrix without copying it out first, and the next chapter leans hard on that.

Fields or columns

The last layout choice has nothing to do with matrices, and it turns up wherever data has named parts. Say every record has an $x$, a $y$ and a $z$. You can store one whole record after another, or you can store all the $x$ together, then all the $y$, then all the $z$. The first is called an array of structs, the second a struct of arrays.

import numpy as np
from timeit import repeat

n = 1 << 22
aos = np.zeros(n, dtype=[('x', 'f8'), ('y', 'f8'), ('z', 'f8')])
aos['x'] = 1.0
soa = np.ones(n, dtype=np.float64)

bench = lambda f: min(repeat(f, number=3, repeat=9)) / 3
print(aos['x'].strides, soa.strides)       # (24,) (8,)
print(bench(lambda: aos['x'].sum()))       # 0.001221
print(bench(lambda: soa.sum()))            # 0.000542

A stride of 24 bytes against a stride of 8, and a factor of 2.25 on the machine this was written on. The struct of arrays is the layout every numerical framework quietly imposes when it asks you for an array of features rather than a list of objects.

This also answers a question Lecture 3 left open. A batch of images with shape $(5, 10, 784)$ is not a nest of Python containers. It is one block of 39200 numbers carrying three strides, and every batched matrix multiply in a forward pass is a promise about them. The shapes tell you what the algebra does. The strides tell you what it costs.

Where this is going

Nothing in this chapter multiplied two matrices. Everything in it decided how fast that will go.

The next chapter takes one line, $C \mathrel{+}= AB$ at size 1024, and drives it from two and a half seconds down to under three milliseconds on one chip, without changing a single arithmetic operation. Every step is a decision about the order the memory is walked, which is to say every step is a sentence from this chapter cashed in. The row walk and the column walk of Figure 5 turn out to be two rungs of that ladder, and the swap between them is the largest single jump on it.