Cache-aware matrix multiplication – naive isn't that bad

http://functionspace.org/articles/40/Cache-aware-Matrix-Multiplication---Naive-isn-t-that-bad-

by m_class • 13 years ago

58 38 13 years ago

Probably anybody who actually implements matrix multiplication knows about ATLAS already, but this software does what the article advocates, and then some, tuning to your specific machine cache: http://en.wikipedia.org/wiki/Automatically_Tuned_Linear_Alge...

khawkins• 13 years ago

I was actually unfamiliar with ATLAS, but I am familiar with Eigen, which I use whenever implementing matrix operations in C++: http://eigen.tuxfamily.org

With so many available options, you should never really be implementing your own matrix multiplication.

temp4534633433• 13 years ago

Figuring out which libary to use is a huge PITA. I don't really understand when I should use uBLAS, mtl, or eigen or whatever else is out there.

Anyone have some insight?

malingo2• 13 years ago

There's a handful of BLAS implementations:

http://en.wikipedia.org/wiki/Basic_Linear_Algebra_Subprogram...

betterunix6• 13 years ago

Naive is going to be "that bad" at that size. I tried this on my own system, 4096x4096 matrix multiplication, with the naive algorithm, the naive algorithm with swapped loops, and Strassen's algorithm. I had my implementation of Strassen's algorithm switch to the naive algorithm (with bad locality of reference) for 64x64 matrices. The results were 530s for the textbook algorithm, 140s for the textbook algorithm with better cache access, and 94s for Strassen's algorithm.

stephencanon5• 13 years ago

For context, using a hand tuned multithreaded and cache-blocked library implementation of O(n^3) multiplication, a 4096x4096 matrix multiply requires less than two seconds on my current-gen Core i5.

malingo2• 13 years ago

You don't always want a multithreaded library implementation though, if the code using the library is already heavily multithreaded.

betterunix6• 13 years ago

The real question is this: how fast would a hand-tuned, multithreaded implementation of Strassen's algorithm be on that system? How fast would Strassen's algorithm be if it used the hand-tuned, multithreaded textbook algorithm for smaller matrices?

I am not denying that constant factor improvements are good, all I am saying is that at the sizes the article is talking about Strassen's algorithm is likely to be advantageous. For what it's worth, having my implementation of Strassen's algorithm use the textbook algorithm with better locality of reference reduces the running time to 61 seconds. Improving the running time of the textbook algorithm will also improve the running time of Strassen's algorithm.

(Edit: It is also worth mentioning that my implementation of Strassen's algorithm could be sped up with a few basic improvements. For example, right now it requires a lot of calls to malloc/free and it needlessly copies matrices. This is not even remotely production-grade software, it was just hastily thrown together.)

6cxs2hd6• 13 years ago

As a general point, it's easy to forget all the layers.

If you learn C, you feel close to the metal. I have an array. It is contiguous bytes in RAM. I've got the power!

Well maybe.

The array might be contiguous virtual memory, but backed by disjoint blocks of physical memory.

Then there's cache. In front of ... another cache.

At least with C you're closeER to the metal.

octo_t2• 13 years ago

And some of your code is in the instruction cache and some isn't - meaning a sudden context[1] switch has to read a lot of code from memory which is slow(ish)

[1] - not a thread switch (although that would be as costly), but an if-statement in a while loop changing from false to true etc.

pbsd3• 13 years ago

Unless it's a friggin' gigantic loop, instruction cache miss will usually not be the issue there; branch misprediction will.

Instruction cache misses are particularly problematic in codebases with many virtual function calls (or plain function pointers in C), which can't be easily prefetched ahead of time.

lettergram• 13 years ago

I remember learning this in one of my courses, it's actually relatively easy to create a program for matrix multiplication. Once you know the size of your cache the speed increase by multiplying just small sections of the matrix is a pretty large increase.

stephencanon5• 13 years ago

It’s generally easy to get a “good enough” matrix multiply (achieving, say, 40-50% of peak). Pressing into the 85+% of peak range usually requires quite a lot of work.

gus_massa• 13 years ago

I think this underestimate the idea of transposing one of the matrices. Perhaps it will be discussed in a future article.

  n = 4096;
  for(i=0;i<n;i++)
  {
      for(j=0;j<n;j++)
      {
           for(k=0;k<n;k++)
              C[i][j]+=A[i][k]*Bt[j][k];
      }
  }

 Transposing a single matrix is quite fast, and gfortran with -O3 does this trick under the hood, so to see the difference you must use -O0 (or perhaps -O1). (And remember that in Fortran the indices are in the other order.)

stephencanon5• 13 years ago

For anyone who wants to learn about this subject in depth (it is a surprisingly interesting subject once you push past the most basic blocking and vectorization techniques), the canonical reference on memory access patterns for truly high-performance matrix multiplication is “Anatomy of High-Performance Multiplication” by Goto and van de Geijn; it’s also a very readable paper.

mturmon• 13 years ago

That is a very detailed review based on first-principles reasoning. Thanks.

I think the second author is van de Geijn. Here's a link to the PDF: http://www.cs.utexas.edu/users/pingali/CS378/2008sp/papers/g...

carterschonwald• 13 years ago

Doing cache oblivious blocking can make a huge 1000x perf difference over naive matrix multiply. And C isn't needed for most of it.

I've some pretty cute sequential dgemm code written in mostly Haskell that's robustly within 5x of openblas dgemm. And I've not tuned it much / at all!

moomin4• 13 years ago

Just to point out the obvious: the thing to take away from this article isn't that this is the right way to do matrix multiplication (ATLAS and BLAS are), but that cache considerations can beat even complexity concerns. Bjarme Stroustroup is quite fond of demonstrating the superiority of vector even when complexity would suggest list was the right approach.

gngeal2• 13 years ago

ATLAS and BLAS are

What you probably meant was that ATLAS and BLAS are implementations of the right way, not the way itself. Because the alternative would be really, really sad.

moomin4• 13 years ago

Actually, I meant that, for most of us, using ATLAS or BLAS is superior to any implementation of our own. I'll allow exceptions for people who have unusual requirements or a phenomenal amount of skill.

betterunix6• 13 years ago

The article showed nothing of the kind. All that it shows is that cache considerations matter, not that they necessarily beat algorithmic improvements.

Really, nobody has ever claimed that constant factor improvements are irrelevant. They certainly do matter: twice as fast is twice as fast no matter what problem size you are talking about. The point of asymptotic analysis is how the running time of an algorithm will change as the input size changes. At scale, even a very bad implementation of Strassen's algorithm will be faster than a highly tuned implementation of the textbook algorithm. The constant factors matter in determining how large the input must be for one algorithm to outperform another.

"Bjarme Stroustroup is quite fond of demonstrating the superiority of vector even when complexity would suggest list was the right approach."

[citation needed]

Most of the time people who claim that vectors are outperforming lists on tasks where lists are "asymptotically better" they are actually hiding some asymptotic cost. Linked lists have poor asymptotic performance for most interesting tasks, just like vectors.

mpyne2• 13 years ago

Bjarne Stroustroup is quite fond of demonstrating the superiority of vector even when complexity would suggest list was the right approach. [citation needed]

You might want to read his Fourth Edition of The C++ Programming Language, where he basically recommends std::vector<> as the default container unless you're dealing with large amounts of data that you add into the vector without using .push_back().

In essence, std::vector<> first, switch to std::list<> if profiling shows a speedup.

Most of the time people who claim that vectors are outperforming lists on tasks where lists are "asymptotically better" they are actually hiding some asymptotic cost. Linked lists have poor asymptotic performance for most interesting tasks, just like vectors.

On the contrary, perhaps the most interesting task we ever do with a list is iterate over it and vectors kick ass there. There's also hardly any indirection or special 'short list' optimization hacks required, such as used in Qt's QList to make its performance reasonable.

moomin4• 13 years ago

Just to be clear, I said cache cOnsiderations can beat complexity, not that it will always beat complexity at the kind of scale that busts the cache.

The interesting point is that you're not necessarily comparing like with like. For the right sized matrix, one is bound by compute speed and one by memory access.

Cache-aware matrix multiplication – naive isn't that bad

Related Stories