Loop Tiling for Memory Hierarchy Optimization

The Problem with Large Datasets: Cache Thrashing

Let's revisit our i-k-j matrix multiplication. The memory access is sequential, which is great for spatial locality. However, consider the access to matrix B. For each outer loop iteration i, we stream through the entire B matrix. If B is 8MB (a 1024x1024 matrix of doubles) and our L2 cache is only 1MB, then by the time we start processing the next i value, the parts of B we need first have long been evicted from the cache.

This phenomenon is sometimes called cache thrashing or cache pollution: we repeatedly fill the cache with data, only to have it pushed out by new data before we can reuse it. Loop tiling is the solution.

The core idea is to break the large matrix multiplication into a series of smaller matrix multiplications on sub-matrices (tiles) that do fit in the cache.

To get a more formal definition and see how this applies to matrix multiplication, the following video explains the concept and the problems it solves.

Implementing Loop Tiling

So, how do we transform our code to implement tiling? We introduce additional outer loops that iterate over the tiles, and the original inner loops now iterate within a tile.

Let's first look at a simple 2D loop, separate from the complexity of matrix multiplication, to see the pattern clearly.

Now, let's apply this pattern to the three nested loops of matrix multiplication. This involves adding tiling loops for multiple dimensions.

The video below provides a formal introduction to this transformation and then analyzes the significant reduction in cache misses that it achieves. Given your background in mathematics, you may find the quantitative analysis of cache misses particularly interesting.

The transformation results in code that looks something like this, using a block size B:

// Tiled Matrix Multiplication
for (int i_tile = 0; i_tile < N; i_tile += B) {
    for (int k_tile = 0; k_tile < N; k_tile += B) {
        for (int j_tile = 0; j_tile < N; j_tile += B) {
            // Mini matrix multiplication on tiles of size B x B
            for (int i = i_tile; i < i_tile + B; i++) {
                for (int k = k_tile; k < k_tile + B; k++) {
                    for (int j = j_tile; j < j_tile + B; j++) {
                        C[i][j] += A[i][k] * B[k][j];
                    }
                }
            }
        }
    }
}
// Note: Boundary checks using min() are omitted for clarity but are needed in a real implementation.

Practical Considerations: Tile Size and Cache Hardware

The theory is elegant, but how does it work in practice? Two key questions arise:

1. How do we choose the optimal tile size?
2. Does our simplified model of the cache hide any real-world complexities?

1. Choosing the Tile Size

The tile size, B, is the most critical parameter.

• If B is too small, the overhead of the extra outer loops can dominate, and we don't get enough data reuse.
• If B is too large, the working set—the three B x B tiles of matrices A, B, and C—won't fit in the target cache level, and we're back to thrashing.

The working set size is approximately 3 * B * B * sizeof(double). If you are targeting your L2 cache, which might be 1MB, you would choose B such that this value is comfortably less than 1MB.

However, the absolute best way to find the optimal tile size is empirical. You test a range of values and measure the performance.

The following resource walks through this exact process.

2. The Reality of Cache Hardware: Conflict Misses

There's a fascinating and subtle issue that can arise. Our analysis assumes a "fully associative" cache, where a block of data from memory can be placed anywhere in the cache. Real-world caches are N-way set-associative, meaning a memory address can only be mapped to a small set of locations within the cache.

What happens if different parts of our data tiles all map to the same cache set? Even if the cache is large enough overall, that one set can fill up, causing data to be evicted. These are called conflict misses. This can happen with matrix operations if the matrix dimensions and tile sizes align in an unlucky way with the cache's geometry.

The same blog post that discussed finding the tile size also has an excellent, detailed explanation of this phenomenon. This is a deeper dive that connects our high-level C++ code directly to the low-level hardware behavior.

Conclusion

In this lesson, we added a crucial and powerful technique to our optimization toolbox. Where loop interchange and unrolling help fine-tune performance for data that's already in the L1 cache, loop tiling manages the flow of data across the entire memory hierarchy.

Key Takeaways:

• Loop Tiling (Blocking) is a loop transformation that partitions a large iteration space into smaller blocks, or tiles.
• The primary goal is to improve temporal locality and data reuse. By processing a small tile of data intensively, we minimize capacity misses to slower levels of memory.
• The size of the tile is chosen so that the working set of data for that tile's computation fits within a targeted level of the cache (e.g., L2 or L3).
• For matrix multiplication, this transforms one large operation into a sequence of many small matrix multiplications, drastically reducing the total data moved from main memory.
• Finding the optimal tile size is often an experimental process, and performance can be affected by low-level hardware details like cache associativity, which can cause conflict misses.

Preview of the Next Lesson

So far, our single-core optimizations have focused on making the most of the memory system and instruction pipeline. In our next lesson, we will explore a different dimension of parallelism available on a single core: data-level parallelism. We will explain the concept of SIMD (Single Instruction, Multiple Data) and learn how modern CPUs can perform a single operation (like an addition or multiplication) on multiple pieces of data simultaneously using vector registers and instructions.