Cache-friendly data

A photo has to be rotated by 90°. The program is two loops: take pixel (y, x), put it into the rotated image at (x, h − 1 − y). However you rearrange the loops, there are as many copies as pixels. Yet the time differs several times over, because memory comes into the cache in lines of 64 bytes, and what matters is how many times each line has to be fetched again. Below is the rotation on a model cache of 4 KB: the image is read on the left and written on the right, and what is in the cache right now is highlighted.

With a row-by-row traversal reading is sequential: each memory line of the source is loaded once. But writing goes down a column of the result: each next pixel goes into a new memory line, and by the time the program comes back to it with the next pixel, the cache has long since thrown it out. Switch the layer to “how many times loaded”: the right image is red, each of its memory lines loaded 16 times, once per pixel. There are 34,000 misses in all, against a minimum of 4,000. Tiles of 16 × 16 bring them down to 5,792.

Now choose a height of 256. The image is only six pixels taller, but tiles stop helping: the misses are back to 34,816, as with the row-by-row traversal. Add padding to the rows, and the misses are exactly 4,096, the minimum. How this happens we will see below, but first the general rule.

The idea

The cache helps a program that touches memory sequentially, comes back to data while they are still in the cache, and keeps together what is needed together. Each of the three rules has its technique: the traversal order, tiling, data layout. None of them changes big O: the number of operations is the same. What changes is the number of trips to memory, and a trip to memory costs a hundred times more than an operation.

How the cache, its lines and sets are built is told in the course chapter “Near and far”, and how much each level of memory costs in the article on latency numbers. Here: what to do about it in code.

Tiling

Why is a row-by-row traversal bad for rotation? For a memory line of the result to be reused, it has to survive in the cache until the next row of the source. Meanwhile the program manages to write one pixel into every row of the result, as many as the image is wide. For a photo 2000 pixels wide that is 2000 memory lines, 128 KB: more than the first-level cache.

Tiles change the count. Take a 32 × 32 square of source pixels and rotate it whole before going on. A row of the square is 32 pixels of 4 bytes, two cache lines, so the source square takes 64 cache lines and the result square as many: 8 KB in all, while a first-level cache is usually 32–48 KB. Until the square is finished, none of its lines is evicted, and each is loaded once. The tile size is chosen so that both squares fit in the cache with room to spare; for 4-byte pixels that is 16–32, sometimes 64. The same technique underlies fast matrix multiplication: linear algebra libraries multiply in blocks sized to the caches.

There is also a way that does not need to know the cache size. Cut the image in half along its long side and rotate each half, cutting it further until the pieces become small. At some level of division the pieces fit in the cache, whatever its size, and from there everything works as with tiles. Such algorithms are called cache-oblivious; they were proposed by Matteo Frigo, Charles Leiserson, Harald Prokop and Sridhar Ramachandran in 1999.

The trap: a row 4 KB long

Now the height of 256. A cache divides memory lines into sets, and a line can go only into its own set, whose number is taken from the address. In the model there are 16 sets of 4 slots, and addresses 1 KB apart fall into the same set. A row of the rotated image is h pixels of 4 bytes. At a height of 250 that is 1000 bytes, and the rows of the result spread over different sets. At 256 it is exactly 1024 bytes, and all 16 rows of the result within a tile fall into one set with four slots. They evict one another, and every write is a miss.

A real processor is the same, only scaled up: a first-level cache usually has 64 sets, and addresses every 4 KB fall into the same set. The dangerous case is a row whose length is a multiple of 4 KB, that is, an image whose height is a multiple of 1024 pixels: 2048, 3072, 4096. Here is the rotation of a real image 2000 pixels wide at heights around 3072:

This measurement was made on the author’s computer: the spike depends on how the cache of a particular processor is built, and on the site’s server it does not show. At a height of 3072 both the row-by-row and the tiled traversal in Go are several times slower than at neighbouring heights: a spike sticks out above the chart. The cure is padding: a row of the result is made slightly longer in memory than needed, and the unused pixels at the end are simply skipped. Graphics libraries store images this way, with the row length in memory (the stride) kept separately from the width. The padding rule: a row should take an odd number of cache lines. Just adding 16 pixels is not enough: at a height of 3056 that gives exactly 3072.

numpy has no spike: rotating with rot90 and copying goes through these heights evenly. And hand-written tiles in Python are slower than the ready function, since each tile costs several calls from the interpreter. For numpy the conclusion is simple: do not reinvent, use what is there.

Layout: array of structs and struct of arrays

The third rule is to keep together what is needed together. A particle model: each particle has coordinates, velocities, a mass and a number, eight fields of 8 bytes, 64 bytes, exactly one cache line. We need the sum of the masses. If the particles lie as an array of structs, each cache line brings one mass and seven unneeded fields. If the masses lie in a separate array, one line brings eight masses.

With a thousand particles they all fit in the cache, and the difference is small. From hundreds of thousands on, the array of structs falls behind: 7–11 times in Go, and 7–16 times in numpy, where a structured array is compared with an ordinary one. This is how columnar databases and formats like Parquet work: an analytical query reads two columns out of forty, and storing a table by columns pays off better than by rows. The opposite case happens too: if code almost always takes all fields of an object at once, an array of structs is better, since one cache line brings everything needed.

Worse than either layout are objects scattered over memory and linked by pointers: every step is a chase after an address, as in the measurement from the article on latencies. That is why an array is almost always walked faster than a linked list, and a B-tree, with many keys in a row in each node, often overtakes a binary search tree even in memory.

Code

The model from the scene: the same cache, the same orders, the same numbers.

The row-by-row and column-by-column traversals cost the same, 8.5 times the minimum: in one case the write misses every time, in the other the read. Tiles of 8 × 8 are worse than tiles of 16 × 16: a cache line holds 16 pixels, and a tile 8 wide uses half the line and has to come back for the other half. At a height of 256 no tiles help, while padding of 16 pixels gives exactly the minimum: each memory line is loaded once.

And this is the real rotation from the measurement: tiles and a padded result row.

What it costs

Big O does not change: rotating an image of N pixels is N copies in any order. To see the difference you need another measure, the number of block transfers between levels of memory. Such a model, the external memory model, was proposed by Alok Aggarwal and Jeffrey Vitter in 1988: fast memory is small, holding M elements, exchanges with slow memory go in blocks of B elements, and only the exchanges are counted. Walking an array of N elements costs N/B exchanges. Rotation row by row costs almost N: every write goes to a new block. Rotation in tiles costs about 2N/B, if an image tile fits in fast memory. The same model explains external sorting and B-trees, only there memory and disk take the place of cache and memory.

The price of these techniques is code complexity and machine dependence: the tile size and the right padding are tuned to the cache, and on another processor the optimum shifts. So the order is this: first find the hot loop by measuring, then estimate how often it goes to memory, and only then rearrange loops.

Problems

The problem is the same cache model, but on a matrix and with real parameters: row-by-row, column-by-column and tiled traversals, and matrices whose rows are exactly 2 KB.