6 ms·
Cache-aware matrix multiplication – naive isn't that bad
- dmlorenzetti 13y agoProbably 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_Algebra_Software http://en.wikipedia.org/wiki/Automatically_Tuned_Linear_Alge...
- khawkins 13y agoI was actually unfamiliar with ATLAS, but I am familiar with Eigen, which I use whenever implementing matrix operations in C++: http://eigen.tuxfamily.org http://eigen.tuxfamily.org With so many available options, you should never really be implementing your own matrix multiplication.
- temp453463343 13y agoFiguring 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?
- pbsd 13y agoI'm not exactly an expert on this kind of thing, but from what I understand, Eigen and Armadillo [1] are the best performing libraries in this space. I'd probably go with Eigen, since it appears to be used in more high-profile places and thus more likely to keep being maintained. uBLAS, on the other hand, does not perform that well in comparison, but it is more likely to keep being maintained as part of Boost. Blaze [2] is supposed to be very high-performance, but it is relatively new, so there's no way to tell where the chips may fall on this one. [1] http://arma.sourceforge.net/ http://arma.sourceforge.net/ [2] https://code.google.com/p/blaze-lib/ https://code.google.com/p/blaze-lib/
- temp453463343 13y agohow about ATLAS ?
- pbsd 13y agoWell, I was comparing apples to apples (i.e., C++ expression template linear algebra libraries). Some (most?) of them [1,2] are linkable against BLAS routines from Intel MKL or ATLAS or whatever. [1] http://eigen.tuxfamily.org/dox/TopicUsingIntelMKL.html http://eigen.tuxfamily.org/dox/TopicUsingIntelMKL.html [2] http://arma.sourceforge.net/faq.html#dependencies http://arma.sourceforge.net/faq.html#dependencies
- temp453463343 13y agoThanks for that clarification
- wallbear 13y agoArmadillo is also used in pretty well known places: NASA, Siemens, Intel, Boeing, Stanford, CMU, MIT, Schlumberger, Deutsche Bank, US Navy, etc. (source: access logs for Armadillo website). It's also the basis for the MLPACK machine learning library: http://mlpack.org/ http://mlpack.org/ Additionally, there are Armadillo bindings to Python and the R language: http://sourceforge.net/projects/armanpy/ http://sourceforge.net/projects/armanpy/ http://cran.r-project.org/web/packages/RcppArmadillo/ http://cran.r-project.org/web/packages/RcppArmadillo/
- malingo 13y agoThere's a handful of BLAS implementations: http://en.wikipedia.org/wiki/Basic_Linear_Algebra_Subprograms#Implementations http://en.wikipedia.org/wiki/Basic_Linear_Algebra_Subprogram...
- betterunix 13y agoNaive 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.
- stephencanon 13y agoFor 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.
- malingo 13y agoYou don't always want a multithreaded library implementation though, if the code using the library is already heavily multithreaded.
- stephencanon 13y agoSure (though it’s worth noting that a correctly-designed high performance library should deliver good performance in that scenario as well; sadly many do not). Single-threaded time on the same system: 5.5 seconds.
- betterunix 13y agoThe 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 13y agoAs 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_t 13y agoAnd 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.
- pbsd 13y agoUnless 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.
- deleted 13y ago[deleted]
- octo_t 13y agoI'll admit that branch misprediction will also screw you over, large event loops are sometimes large enough to escape the instruction cache (especially with function inlining)
- stephencanon 13y agoFor 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 13y agoThat 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/gotoPaper.pdf http://www.cs.utexas.edu/users/pingali/CS378/2008sp/papers/g...
- lettergram 13y agoI 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.
- stephencanon 13y agoIt’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 13y agoI 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.)
- carterschonwald 13y agoDoing 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!
- deleted 13y ago[deleted]
- moomin 13y agoJust 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.
- gngeal 13y agoATLAS 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.
- moomin 13y agoActually, 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.
- gngeal 13y agoI may be mistaken, but don't ATLAS, BLAS, LINPACK etc. only operate on single and double precision floats? What if the elements you're operating on are, say, rationals or polynomials? Somehow it seemed to me that linear algebra is a somewhat more abstract concept than what Fortran numerical libraries usually tend to implement. Also, the drawback of anything written in this way is that it doesn't compose well. The individual operations may well have been optimized into oblivion, but their composition may not necessarily be what you'd like it to be. One can also say goodbye to higher-order functions. But I'm not exactly a Fortran person when it comes to code style, so I may be prejudiced against overly loopy code. (BTW, am I the only one who thinks that the Sourceforge (http://sourceforge.net/projects/math-atlas/ http://sourceforge.net/projects/math-atlas/) comments are machine-generated? I mean, comments like "An essential Windows program." or "The best program that I've ever used." are hardly what I'd expect for a numerical library.)