12 ms·

How to write efficient matrix multiplication

400 points by Mtz1974 8y ago

romwell 8y ago
Sorry, but this is not at all how to write efficient matrix multiplication. The correct answer is: don't. Use BLAS (yes, written in Fortran) + OpenMP (if needed), and let the smart people handle this. (Unless you are one of those smart people - in which case, you should improve BLAS, which the rest of us mortals will use). The longer answer is: matrix multiplication is a highly non-trivial algorithmic problem[1], which is, in fact, still open. The asymptotically fast algorithms (such as Strassen's[2]) outperform the naive algorithm (i.e. by definition, such is in the article) when matrices get large. After skimming the article it is still unclear to me how the optimized naive algorithm fares against Strassen's, for example. The final bit of advice this article is offers is reading the paper on... implementation of BLAS (that's what Goto describes there). And so that's basically the case where avoiding Goto is considered harmful. [1]https://en.wikipedia.org/wiki/Matrix_multiplication_algorithm#Sub-cubic_algorithms https://en.wikipedia.org/wiki/Matrix_multiplication_algorith... [2]https://en.wikipedia.org/wiki/Strassen_algorithm https://en.wikipedia.org/wiki/Strassen_algorithm [3]https://en.wikiquote.org/wiki/Donald_Knuth https://en.wikiquote.org/wiki/Donald_Knuth
- forrrealman 8y ago
  Why do people keep spreading this nonsense? 1. All BLAS implementations I know of don't use Strassen's algorithm. 2. This is such a Stack Overflow answer, please make it stop. Q: "How do I do a thing?" A: "I have unilaterally decided your question is invalid and you should not do a thing." That's really useful! To the author: hack away. Goto can probably write better ASM but tutorials like this are very helpful for people who'd just like to read some interesting case studies.
  - beevai142 8y ago
    But the stackoverflow answer is right, and accounts for the fact that in the majority of cases the person who is asking a confused question is indeed confused. If what you are trying to do is to get best performance, make use of work that other people already did, rather than wasting your time on a solved problem. If you want to learn about how efficient matrix multiplication can be implemented, that is a different problem.
    - forrrealman 8y ago
      Ahhh make it stop! By the way "make use of work that other people already did" might be nice for getting something built fast but it may not be the best thing. Any issue that involves writing fast code involves tradeoffs. If you sit down to write something new you may have different views on the tradeoffs than whoever wrote "the fastest" one. Life ain't so black and white. In fact even these 'best of field' products tend to be ugly inside and unoptimized in places. (source: I develop linear algebra libraries) Bonus: for low level ASM math, every "solved" problem (which by the way it wasn't) becomes unsolved the second Intel or AMD or whoever releases a new chip or coprocessor.
      - throwawaymath 8y ago
        > (source: I develop linear algebra libraries) This is something I'm interested in contributing to. Can you name a few libraries (especially ones implementing new and interesting work) that would welcome open source contributors? Alternatively you can just contact me (via my profile info) if you're working on something in particular but would rather not be identified publicly.
        forrrealman 8y ago
        I don't see any contact info in your profile, but I would be happy to follow up via email.
        throwawaymath 8y ago
        Whoops, sorry about that. Fixed :)
        deleted 8y ago
        [deleted]
    - tzs 8y ago
      The problem with Stackoverflow answers of the form "You asked about doing X, but it looks like you are trying to do X to achieve Y, and X is not a good way to achieve Y. You want to do Z. Here's how to do Z" is that later when people who actually do need to do X are looking for help, their searches keep turning up these damned answers that don't tell how to do X. SO needs a flag on questions that is set if there is an answer to the actual question as asked and clear otherwise, and that can be used as a search filter.
- lvlnr 8y ago
  Could it, perhaps, be a tutorial that imparts some useful knowledge instead of being a suggestion that everyone should roll their own dense linear algebra? I mean, most tutorials shouldn't really exist - why would you code anything, when there is a library for it?
  - bo1024 8y ago
    The problem with this page (for me) is that it's hard to evaluate as a tutorial. It reads like describing an idea the author came up with, not how state of the art implementations work.
    - iainmerrick 8y ago
      It ends with citing and recommending a paper by another person; that doesn’t sound like “an idea the author came up with”! It’s the author trying to describe the state of the art, as they see it.
- acqq 8y ago
  Note what the article author writes: "On my machine the code runs at around 100 Gflops, which is not far from the theoretical peak, and similar to Intel's optimized implementation." So he knows that, for example, Intel provides/sells already optimized implementation. And he knows that there is GotoBLAS (1), because as you noted, he cites the paper about it. But he explains the implementation present in Glow, "a machine learning compiler and execution engine for various hardware targets." (2) Now if you want to prove to the authors of Glow why they should have "just used BLAS" I would like to read this discussion (or even better, see a proof of concept code with benchmarks etc). But I can imagine that in their case there were some reasons why not to use it, as they were aware that it exists. However it is true that the readers should be aware of the background assumptions and to be aware of the existence of the libraries. 1) https://en.wikipedia.org/wiki/GotoBLAS https://en.wikipedia.org/wiki/GotoBLAS 2) https://github.com/pytorch/glow https://github.com/pytorch/glow
- danmg 8y ago
  A hand rolled routine where you benchmarked different indexing orders, strides, interleaving on a particular sized matrix on fixed sizes of caches has the potential to outperform blas. In general, use BLAS. Those sub-cubic algorithms require the matrix to be very-very-very large to recoup the extra constant costs associated with them.
  - sampo 8y ago
    > Those sub-cubic algorithms require the matrix to be very-very-very large to recoup the extra constant costs associated with them. Wikipedia disagrees: "In practice, Strassen's algorithm can be implemented to attain better performance than conventional multiplication even for small matrices, for matrices that are not at all square, and without requiring workspace beyond buffers that are already needed for a high-performance conventional multiplication." https://en.wikipedia.org/wiki/Strassen_algorithm https://en.wikipedia.org/wiki/Strassen_algorithm
    - danmg 8y ago
      That's not what I learned in grad school when I took a class on this.
      - sdenton4 8y ago
        I used to use Atlas quite a bit, which effectively benchmarks your CPU and memory configuration at installation time to empirically determine where to switch from 'small n' algorithms to asymptomatically optimal algorithms, like Strassen. Of course, if you only some resources into optimizing the cubic time version, Strassen will look slow by comparison. Better we take a serious look at both and empirically figure out where the trade off point is...
        danmg 8y ago
        Yes atlas is good. For small, dense n you're not going to beat the cubic algorithm. This is the most common case in a lot of applications. For sparse n, use a sparse version of the cubic algorithm. For large dense n, those sub-cubic algorithms might be better. Here n would have to be so large that you'd have to factor in the cost of doing it out of core.
- andrepd 8y ago
  I hardly think the purpose of this post was to show "hey, here's what you should do if you want to multiply matrices in production: write this from scratch". More like "learn how it is done, to make a seemingly simple algorithm, for a crucial problem, fast on modern processors".
- syllogism 8y ago
  Oh, interesting point! Yes let's just use OpenBLAS, no need for anyone but one or two people to understand these things. Sure. Okay so let's have a look at OpenBLAS's Skylake kernel. Maybe we can learn some deep wisdom this author hasn't considered! Go and have a look for it, I'll wait here. Ah, what's that? THERE IS STILL NO SKYLAKE KERNEL IN OPENBLAS? Interesting. In fact if you do "sudo apt-get install libopenblas" on Ubuntu 16.04 (the one every deep learning tutorial will be recommending), you'll find you get version 2.18 of OpenBLAS. You'll find that this version of the software is quite old now, and it has a bug that causes it to fail to detect your CPU correctly. You'll find that instead of using any AVX instructions, it'll fall-back to the Prescott kernel. You'll also find this is really fucking slow. In summary: a) You have no idea what you're talking about; b) Even if you were right on these points of detail, your point of view is still terrible c) Please avoid being so beligerantly wrong in future
  - madez 8y ago
    Let me add that there exists a similar echo chamber when it comes to questions regarding cryptography. Way too often questions about how to do this or that in cryptography are put down by a "you don't" or the generic "don't roll your own crypto", which is hostile against curiosity, learning, and knowledge and on top of that unfriendly towards the one asking. Often when I read it, I feel the hacker ethos is being insulted. While I think there are reasonable arguments for the security of (certain) self-made crypto, taking that security aspect completely aside, we need new people in every field, and cryptography is no exception from that. We should embrace newcomers, not tell me to not do what would teach them so much. Admittingly, one can easily mess up, but then tell how and what to take care of instead of just telling them not to try anything at all.
    - throwawaymath 8y ago
      I feel like the admonition, "don't roll your own crypto" refers to production code. If you want to develop your own crypto for "curiosity, learning and knowledge" I don't think many people in the security community will really bother you about it. Just don't use it for securing anything in production unless it's been thoroughly audited and you have a reason not to use something that already exists. It seems pretty obvious why you shouldn't design or implement novel cryptography in production without good reason and significant expertise. It also doesn't seem like the people saying this need to spell out that this doesn't preclude learning exercises or legitimate research which won't endanger any real world information. So what's the controversy?
- deleted 8y ago
  [deleted]
- vmarsy 8y ago
  > The asymptotically fast algorithms (such as Strassen's[2]) outperform the naive algorithm (i.e. by definition, such is in the article) when matrices get large. This is true only if approaching the matrix multiplication problem purely from a theoretical perspective. But the reality is that computers have a specific architecture designed to be as fast as possible, while still respecting the laws of physics. The libraries such as BLAS deal with those architectures, and it happens that there, an O(n³) triple loop with the optimizations described in the Goto paper is what gives you the fastest performance. Maybe if we were dealing with 1b*1b matrices, Strassen's like algorithms would start to shine, but with those sizes you already exceed the maximum space on a computer, in which case you then start looking at parallel algorithms on clusters/supercomputers where the algorithms become very different (see things such as SUMMA: Scalable Universal Matrix Multiplication Algorithm). In those cases, the computation happening on each processor would still be the fast O(n³) BLAS code.
vecplane 8y ago
Somewhat related: check out gl-matrix for some fast matrix math in JavaScript - https://github.com/toji/gl-matrix https://github.com/toji/gl-matrix
- svantana 8y ago
  This looks like it's only for 3x3 and 4x4 matrices, e.g. what you usually use in 3d graphics. Or is it possible to do arbitrary matrix multiplication with this lib?
  - Lerc 8y ago
    It does not appear to be the case. On the other hand, optimal speed and JavaScript are often different ball parks. I threw together this thing for when I need some matrix stuff in a little script. Not fast, but flexible. var vop = op=>( (a,b)=>( a.map((v,i)=>op(v,b[i])) ) ); var vdiff = vop((a,b)=>a-b); var vadd = vop((a,b)=>a+b); var vdot= (a,b)=>a.reduce( (ac,av,i)=>ac+=av*b[i],0); var vlength = a=>Math.sqrt(vdot(a,a)); var vscale = (a,b)=>a.map(v=>v*b); var vdistance = (a,b)=>vlength(vdiff(a,b)); var vnormalised = (a,b=1)=>vscale(a,b/vlength(a)); var project = (point, matrix) => matrix.map( p=>vdot(p,[...point,1])); var transpose = matrix=> ( matrix.reduce(($, row) => row.map((_, i) => [...($[i] || []), row[i]]), []) ); var multiply = (a,b,...rest)=>(!b)?a:multiply(a.map( (p=>transpose(b).map(q=>vdot(p,q)))),...rest); I look forward to the day when the sufficiently-smart JIT implements the optimal Matrix Multiplication from that. :-)
    - svantana 8y ago
      Thanks! I find it weird that performance-focused javascript projects don't have live benchmarks up so you can judge for yourself. Personally I've found that a semi-clever matrix multiplier in c++ compiled with emscripten can give about 1.5 gflops (singlethreaded on a 2.5GHz i7 macbook pro), and existing js libs like sushi are a lot slower.
- megalodon 8y ago
  For anyone working with linear algebra in Node.js, I wrote C++ bindings to BLAS a couple of years ago [0]. I also started working on bindings to LAPACK, but never had the time to finish them [1]. [0] https://github.com/mateogianolio/nblas https://github.com/mateogianolio/nblas [1] https://github.com/mateogianolio/nlapack https://github.com/mateogianolio/nlapack
gtycomb 8y ago
Somewhere in the conversation Les Valiant's paper "General Context-free recognition in less than Cubic time" is worth mentioning. This has informative and clear discussion on the complexity of matrix-multiplication in programming areas such as grammar and parsing.
svantana 8y ago
I'm a big fan of writing simple, portable code that is structured enough to let the compiler optimize it to perfection. If you do it right, a naive A' * B multiplier can be within 20% of BLAS performance, at least for non-huge matrices. Now, when you need that code on some embedded platform or for emscripten or whatever, it'll just work and be really tiny too.
- PeCaN 8y ago
  Unfortunately, matrix multiplication is the kind of thing where you often need it to be fast and scale to huge matrices. I like simple code as much as the next guy, but occasionally you do need to write fast code.
- ummonk 8y ago
  It's the huge matrices where performance matters, and a naive multiplier for that will use memory rather inefficiently.
  - svantana 8y ago
    > It's the huge matrices where performance matters I wouldn't say that, just look at CGI where a lot of effort is spent making fast 4x4 matrix multiplies. But yeah, a general matrix mult needs to scale, but a lot of the time you know what size you will run and you can optimize for that.
- ziotom78 8y ago
  Except that, when you are writing scientific code to solve a large linear problem, like a weather model or asupernova simulation, the difference between some naïve code and a heavily optimized one can be ten hours against ten days...
  - carlmr 8y ago
    If it's within 20% of BLAS performance it shouldn't be the difference between ten hours and ten days. That must be something else then.
nickodell 8y ago
>The first thing that we need to do to accelerate our program is to improve the memory utilization of matrix B. Why not put matrix B in column-major order?
- costrouc 8y ago
  Matrix multiplication involves moving across one matrix in column order and the other matrix in row order. So it turns out that both row or column ordering make no difference. I think that matrix multiplication is one of the best examples of a deceptivly simple problem. It shows how far ALL code is from peak performance. We can only strive to be within an order or two from peak performance.
  - fdej 8y ago
    This trick does work. If the matrices are in row-major order, you transpose B in memory and then compute A * (B^T)^T. This multiplication reads both matrices in row order. However, while this does improve performance over the naive algorithm, it's still not as good as a tiling algorithm.
    - shaklee3 8y ago
      I've found where that causes problems is when on of the matrix dimensions is not a multiple of the cache line size. It's common on gpus to use more elements then there are in the dimension. Nvidia calls this the leading dimension, and it must be greater than or equal to the Matrix dimension. If this is the case, the transpose trick doesn't quite work anymore.
  - hedora 8y ago
    I took a different lesson from learning to implement matrix multiply: you can get a small amount of code within a few percent of peak performance if you try. So, the trick is to keep the 99% of the code that does 10% of the work within an order of magnitude of peak performance, and then get the rest within a factor of two. That gets the system within a factor of three, which is still embarrassing, but that’s what profiling tools are for. :-)
whos3424 8y ago
yo
whos3424 8y ago
please look at http://halide-lang.org/ http://halide-lang.org/ and especially https://www.youtube.com/watch?v=3uiEyEKji0M https://www.youtube.com/watch?v=3uiEyEKji0M they seem to open new ground in the area you are prusuing edit: if you know any other groups/ppl who practice a similar approach send them my way anyhow, it seems like the hand written assembly vs compiled code goes in to a new skirmish the age where the whole spectrum from generalcpu > asic becomes available to programmers all around bottle necks of parallelism, locality and redundancy seems to be the guiding lights hand written we'll be probably lead by the magnificent http://www.greenarraychips.com/home/documents/pub/AP001-MD5.html http://www.greenarraychips.com/home/documents/pub/AP001-MD5.... compiler generated by halide and the like and once we got the ball rolling and more public knowledge is available probably the locked down techniques of the major supercomputers will start to emerge ( i just assume pps from the supercomputers area have been dealing with this stuff for the last couple of decades ) can i spit some more, dear hacker news? focus, you gotta focus, forget about the compiler generating good code, it's all about the ability for us humans to explore the domain even in halide, they mention the point that you still need to sit down and code different approaches, but their tool makes it a breeze compared to exploring those different approaches using hand written assembly
matachuan 8y ago
tiling and vectorizing are well known.
deleted 8y ago
[deleted]
nautilus12 8y ago
How would this change if you were dealing with more sparse matrices? Im used to using spark mllib distributed to do matrix multiplication on large sparse coordinate matrices and it works fine for me but i dont think i need the level of optimality this approach would demand.
- knlji 8y ago
  Doesn't the multiplication of two sparse matrices increase fill-in? That cannot be a good idea in general.
  - stephencanon 8y ago
    Absolutely, but if that’s what you need to compute, that’s what you need to compute. Sometimes it is.
- hadibrais 8y ago
  Sparse matrix multiplication is generally considered to be a different problem. Sparse matrices need to be represented differently (See https://en.wikipedia.org/wiki/Sparse_matrix https://en.wikipedia.org/wiki/Sparse_matrix for an introduction). There has been recently a lot of research on how to do sparse matrix multiplication efficiently on manycore architectures. See this 2018 paper for example: https://arxiv.org/pdf/1804.01698.pdf https://arxiv.org/pdf/1804.01698.pdf. Note that sparse matrix-vector multiplication is also a different problem.
  - nautilus12 8y ago
    Thats fascinating, i wonder how much coordinate matrix multiplication in spark takes into account factors like choelesky decomposition. Noticed most of those libraries are written in C++ or fortran, is JVM based approaches generally discouraged in this realm? It fits well with our current stack based mostly around spark
    - joshuamorton 8y ago
      Not quite the same, but similar, see the flowcharts at the end of this [1]. [1]: https://www.mathworks.com/help/matlab/ref/mldivide.html https://www.mathworks.com/help/matlab/ref/mldivide.html
taeric 8y ago
I was really expecting this to just be "import library". That said, very well done article. Would be interesting to see this style expand to the more ambitious methods in https://en.wikipedia.org/wiki/Matrix_multiplication_algorithm#Sub-cubic_algorithms https://en.wikipedia.org/wiki/Matrix_multiplication_algorith...
vmarsy 8y ago
Writing the fastest matrix multiplication you can is pretty fun. It probably won't be as fast a BLAS or MKL, but if you get close to it, it's a rewarding experience. > see what else can be done I recommend reading this paper by Kazushige Goto After doing the easiest change of loop reordering for already significant perf improvements, the paper linked ("Anatomy on a high performance matrix multiplication") is really excellent. It will walk you through all the optimizations like tiling for reducing L1, L2, L3, and TLB cache misses, and leverage vectorization. Then for squeezing out even more perf I remember there's the things like loop prefetching, loop unrolling, which you would expect the compiler with the best optimizations flags even targeted for the native architecture to take care of automatically. But you realize that it's not necessarily true.
- comicjk 8y ago
  If your particular matrix problem has special properties that aren't covered by blas, you can wring out an apparently impossible improvement. I once found a matrix problem (periodic cubic splines) which turned out to have a solution in O(N) - the matrix columns were related such that it could be replaced by a vector.
  - yonkshi 8y ago
    I can attest this too, I was recently working with some massive and sparse binary matrix, with only 0.1% elements have values, it was magnitudes faster to index and sum the elements than to do a dot product. Not even a Tesla GPU was close to the performance of a custom matrix operation we wrote.
    - repsilat 8y ago
      A lot of matrix routines will have sparse matrix support. Though often the line between wanting a sparse matrix and wanting some kind of adjacency-list graph structure can be pretty fine... Also, only tangentially related, but if you can change your problem formulation to make the matrix of interest sparse(r), that can be a huge win. Two things that have helped me in the past: - Rows or columns that have mostly the same (nonzero) coefficient throughout. You can normally do some simple substitution to turn them into zeros. - Rows or columns that are naturally very similar to each other. You can often replace `y` with `y-x` or replace `x[i]` with `x[i]-x[i+1]` to turn one (or `n-1`) columns "mostly empty."
tehsauce 8y ago
A really nice video with animations of various scheduling approaches for these types of linear algebra algorithms: https://youtu.be/3uiEyEKji0M https://youtu.be/3uiEyEKji0M
fibo 8y ago
I wrote a JavaScript implementation with the naive method, but it also checks that matrix are compatibile and generalize the underlying ring operation (mul and add) so you could use it to multiply complex matrices, or bigint matrices, or every other ring you define. http://g14n.info/matrix-multiplication/ http://g14n.info/matrix-multiplication/
aqme28 8y ago
I'm surprised to see no discussion of faster-than-naive algorithms like Strassen[1] or Coppersmith-Winograd[2], which very counterintuitively run faster than N^3. Note: Only helpful for specific or very large matrices. [1]: https://en.wikipedia.org/wiki/Strassen_algorithm https://en.wikipedia.org/wiki/Strassen_algorithm [2]: https://en.wikipedia.org/wiki/Coppersmith%E2%80%93Winograd_algorithm https://en.wikipedia.org/wiki/Coppersmith%E2%80%93Winograd_a...
- nwallin 8y ago
  The problem with those algorithms is that they define "work" as "how many multiplications do I have to do?" Addition and memory accesses are defined by the analysis to be free. And you end up doing 4-5 times as many floating point additions as the naive algorithm, and the memory accesses are decidedly non-cache friendly. And modern CPUs have two multipcation units but only one addition unit, so multiplications are actually faster than addition. Plus, there are ergonomic issues with requiring square powers-of-two matrices. Matrices have to be truly spectacularly large to benefit from Strassens, and applications of the asymptotically faster algorithms don't really exist. The matrices have to be so large that they're computationally too expensive regardless of the asymptotic speedup. You're talking about dispatching work units to a cluster at that point. Most importantly though, the discussion is... boring. It's just "yup, we do Strassens, then we quickly reach the point where we go back to the naive algorithm, and things are interesting again."
  - repsilat 8y ago
    > Matrices have to be truly spectacularly large to benefit from Strassens I don't have personal experience with other mmult methods, but Wikipedia says Strassen's algorithm is at least beneficial at matrix sizes that fit in cache. I have had occasion to deal with matrices of large dimension (though usually sparse, and usually not square) and those numbers don't seem crazy at all. Are they that far off-base? (Wikipedia does say that Coppersmith-Winograd tends to never be a useful improvement.) In any case, perhaps a more friendly response might be to point out that Strassen's algorithm tends to "bail out" into regular old matrix multiplication for small `n` anyway, so this still helps provide a speedup.
fabrice_d 8y ago
Related: https://news.ycombinator.com/item?id=17058189 https://news.ycombinator.com/item?id=17058189
phlip9 8y ago
A similar article on how to optimize SGEMM on GPUs using OpenCL - https://cnugteren.github.io/tutorial/pages/page1.html https://cnugteren.github.io/tutorial/pages/page1.html
gnufx 8y ago
This thread is old, but for the sake of archives: BLIS actually tells you how to write a fast production large-matrix GEMM, and the papers linked from https://github.com/flame/blis https://github.com/flame/blis would be a better reference than the Goto and van de Geijn. For small matrices see the references from https://github.com/hfp/libxsmm https://github.com/hfp/libxsmm but you're not likely to be re-implementing that unless you're serious about a version on non-x86_64. See also https://news.ycombinator.com/item?id=17061947 https://news.ycombinator.com/item?id=17061947
- gnufx 8y ago
  I could also have referenced the original(?) take on this: https://github.com/flame/how-to-optimize-gemm https://github.com/flame/how-to-optimize-gemm

400 points by Mtz1974 8y ago

Sorry, but this is not at all how to write efficient matrix multiplication. The correct answer is: don't. Use BLAS (yes, written in Fortran) + OpenMP (if needed), and let the smart people handle this. (Unless you are one of those smart people - in which case, you should improve BLAS, which the rest of us mortals will use). The longer answer is: matrix multiplication is a highly non-trivial algorithmic problem[1], which is, in fact, still open. The asymptotically fast algorithms (such as Strassen's[2]) outperform the naive algorithm (i.e. by definition, such is in the article) when matrices get large. After skimming the article it is still unclear to me how the optimized naive algorithm fares against Strassen's, for example. The final bit of advice this article is offers is reading the paper on... implementation of BLAS (that's what Goto describes there). And so that's basically the case where avoiding Goto is considered harmful. [1]https://en.wikipedia.org/wiki/Matrix_multiplication_algorithm#Sub-cubic_algorithms https://en.wikipedia.org/wiki/Matrix_multiplication_algorith... [2]https://en.wikipedia.org/wiki/Strassen_algorithm https://en.wikipedia.org/wiki/Strassen_algorithm [3]https://en.wikiquote.org/wiki/Donald_Knuth https://en.wikiquote.org/wiki/Donald_Knuth

forrrealman 8y ago
Why do people keep spreading this nonsense? 1. All BLAS implementations I know of don't use Strassen's algorithm. 2. This is such a Stack Overflow answer, please make it stop. Q: "How do I do a thing?" A: "I have unilaterally decided your question is invalid and you should not do a thing." That's really useful! To the author: hack away. Goto can probably write better ASM but tutorials like this are very helpful for people who'd just like to read some interesting case studies.
- beevai142 8y ago
  But the stackoverflow answer is right, and accounts for the fact that in the majority of cases the person who is asking a confused question is indeed confused. If what you are trying to do is to get best performance, make use of work that other people already did, rather than wasting your time on a solved problem. If you want to learn about how efficient matrix multiplication can be implemented, that is a different problem.
  - forrrealman 8y ago
    Ahhh make it stop! By the way "make use of work that other people already did" might be nice for getting something built fast but it may not be the best thing. Any issue that involves writing fast code involves tradeoffs. If you sit down to write something new you may have different views on the tradeoffs than whoever wrote "the fastest" one. Life ain't so black and white. In fact even these 'best of field' products tend to be ugly inside and unoptimized in places. (source: I develop linear algebra libraries) Bonus: for low level ASM math, every "solved" problem (which by the way it wasn't) becomes unsolved the second Intel or AMD or whoever releases a new chip or coprocessor.
    - throwawaymath 8y ago
      > (source: I develop linear algebra libraries) This is something I'm interested in contributing to. Can you name a few libraries (especially ones implementing new and interesting work) that would welcome open source contributors? Alternatively you can just contact me (via my profile info) if you're working on something in particular but would rather not be identified publicly.
      - forrrealman 8y ago
        I don't see any contact info in your profile, but I would be happy to follow up via email.
        throwawaymath 8y ago
        Whoops, sorry about that. Fixed :)
        deleted 8y ago
        [deleted]
  - tzs 8y ago
    The problem with Stackoverflow answers of the form "You asked about doing X, but it looks like you are trying to do X to achieve Y, and X is not a good way to achieve Y. You want to do Z. Here's how to do Z" is that later when people who actually do need to do X are looking for help, their searches keep turning up these damned answers that don't tell how to do X. SO needs a flag on questions that is set if there is an answer to the actual question as asked and clear otherwise, and that can be used as a search filter.
lvlnr 8y ago
Could it, perhaps, be a tutorial that imparts some useful knowledge instead of being a suggestion that everyone should roll their own dense linear algebra? I mean, most tutorials shouldn't really exist - why would you code anything, when there is a library for it?
- bo1024 8y ago
  The problem with this page (for me) is that it's hard to evaluate as a tutorial. It reads like describing an idea the author came up with, not how state of the art implementations work.
  - iainmerrick 8y ago
    It ends with citing and recommending a paper by another person; that doesn’t sound like “an idea the author came up with”! It’s the author trying to describe the state of the art, as they see it.
acqq 8y ago
Note what the article author writes: "On my machine the code runs at around 100 Gflops, which is not far from the theoretical peak, and similar to Intel's optimized implementation." So he knows that, for example, Intel provides/sells already optimized implementation. And he knows that there is GotoBLAS (1), because as you noted, he cites the paper about it. But he explains the implementation present in Glow, "a machine learning compiler and execution engine for various hardware targets." (2) Now if you want to prove to the authors of Glow why they should have "just used BLAS" I would like to read this discussion (or even better, see a proof of concept code with benchmarks etc). But I can imagine that in their case there were some reasons why not to use it, as they were aware that it exists. However it is true that the readers should be aware of the background assumptions and to be aware of the existence of the libraries. 1) https://en.wikipedia.org/wiki/GotoBLAS https://en.wikipedia.org/wiki/GotoBLAS 2) https://github.com/pytorch/glow https://github.com/pytorch/glow
danmg 8y ago
A hand rolled routine where you benchmarked different indexing orders, strides, interleaving on a particular sized matrix on fixed sizes of caches has the potential to outperform blas. In general, use BLAS. Those sub-cubic algorithms require the matrix to be very-very-very large to recoup the extra constant costs associated with them.
- sampo 8y ago
  > Those sub-cubic algorithms require the matrix to be very-very-very large to recoup the extra constant costs associated with them. Wikipedia disagrees: "In practice, Strassen's algorithm can be implemented to attain better performance than conventional multiplication even for small matrices, for matrices that are not at all square, and without requiring workspace beyond buffers that are already needed for a high-performance conventional multiplication." https://en.wikipedia.org/wiki/Strassen_algorithm https://en.wikipedia.org/wiki/Strassen_algorithm
  - danmg 8y ago
    That's not what I learned in grad school when I took a class on this.
    - sdenton4 8y ago
      I used to use Atlas quite a bit, which effectively benchmarks your CPU and memory configuration at installation time to empirically determine where to switch from 'small n' algorithms to asymptomatically optimal algorithms, like Strassen. Of course, if you only some resources into optimizing the cubic time version, Strassen will look slow by comparison. Better we take a serious look at both and empirically figure out where the trade off point is...
      - danmg 8y ago
        Yes atlas is good. For small, dense n you're not going to beat the cubic algorithm. This is the most common case in a lot of applications. For sparse n, use a sparse version of the cubic algorithm. For large dense n, those sub-cubic algorithms might be better. Here n would have to be so large that you'd have to factor in the cost of doing it out of core.
andrepd 8y ago
I hardly think the purpose of this post was to show "hey, here's what you should do if you want to multiply matrices in production: write this from scratch". More like "learn how it is done, to make a seemingly simple algorithm, for a crucial problem, fast on modern processors".
syllogism 8y ago
Oh, interesting point! Yes let's just use OpenBLAS, no need for anyone but one or two people to understand these things. Sure. Okay so let's have a look at OpenBLAS's Skylake kernel. Maybe we can learn some deep wisdom this author hasn't considered! Go and have a look for it, I'll wait here. Ah, what's that? THERE IS STILL NO SKYLAKE KERNEL IN OPENBLAS? Interesting. In fact if you do "sudo apt-get install libopenblas" on Ubuntu 16.04 (the one every deep learning tutorial will be recommending), you'll find you get version 2.18 of OpenBLAS. You'll find that this version of the software is quite old now, and it has a bug that causes it to fail to detect your CPU correctly. You'll find that instead of using any AVX instructions, it'll fall-back to the Prescott kernel. You'll also find this is really fucking slow. In summary: a) You have no idea what you're talking about; b) Even if you were right on these points of detail, your point of view is still terrible c) Please avoid being so beligerantly wrong in future
- madez 8y ago
  Let me add that there exists a similar echo chamber when it comes to questions regarding cryptography. Way too often questions about how to do this or that in cryptography are put down by a "you don't" or the generic "don't roll your own crypto", which is hostile against curiosity, learning, and knowledge and on top of that unfriendly towards the one asking. Often when I read it, I feel the hacker ethos is being insulted. While I think there are reasonable arguments for the security of (certain) self-made crypto, taking that security aspect completely aside, we need new people in every field, and cryptography is no exception from that. We should embrace newcomers, not tell me to not do what would teach them so much. Admittingly, one can easily mess up, but then tell how and what to take care of instead of just telling them not to try anything at all.
  - throwawaymath 8y ago
    I feel like the admonition, "don't roll your own crypto" refers to production code. If you want to develop your own crypto for "curiosity, learning and knowledge" I don't think many people in the security community will really bother you about it. Just don't use it for securing anything in production unless it's been thoroughly audited and you have a reason not to use something that already exists. It seems pretty obvious why you shouldn't design or implement novel cryptography in production without good reason and significant expertise. It also doesn't seem like the people saying this need to spell out that this doesn't preclude learning exercises or legitimate research which won't endanger any real world information. So what's the controversy?
deleted 8y ago
[deleted]
vmarsy 8y ago
> The asymptotically fast algorithms (such as Strassen's[2]) outperform the naive algorithm (i.e. by definition, such is in the article) when matrices get large. This is true only if approaching the matrix multiplication problem purely from a theoretical perspective. But the reality is that computers have a specific architecture designed to be as fast as possible, while still respecting the laws of physics. The libraries such as BLAS deal with those architectures, and it happens that there, an O(n³) triple loop with the optimizations described in the Goto paper is what gives you the fastest performance. Maybe if we were dealing with 1b*1b matrices, Strassen's like algorithms would start to shine, but with those sizes you already exceed the maximum space on a computer, in which case you then start looking at parallel algorithms on clusters/supercomputers where the algorithms become very different (see things such as SUMMA: Scalable Universal Matrix Multiplication Algorithm). In those cases, the computation happening on each processor would still be the fast O(n³) BLAS code.

vecplane 8y ago

Somewhat related: check out gl-matrix for some fast matrix math in JavaScript - https://github.com/toji/gl-matrix https://github.com/toji/gl-matrix

svantana 8y ago
This looks like it's only for 3x3 and 4x4 matrices, e.g. what you usually use in 3d graphics. Or is it possible to do arbitrary matrix multiplication with this lib?
- Lerc 8y ago
  It does not appear to be the case. On the other hand, optimal speed and JavaScript are often different ball parks. I threw together this thing for when I need some matrix stuff in a little script. Not fast, but flexible. var vop = op=>( (a,b)=>( a.map((v,i)=>op(v,b[i])) ) ); var vdiff = vop((a,b)=>a-b); var vadd = vop((a,b)=>a+b); var vdot= (a,b)=>a.reduce( (ac,av,i)=>ac+=av*b[i],0); var vlength = a=>Math.sqrt(vdot(a,a)); var vscale = (a,b)=>a.map(v=>v*b); var vdistance = (a,b)=>vlength(vdiff(a,b)); var vnormalised = (a,b=1)=>vscale(a,b/vlength(a)); var project = (point, matrix) => matrix.map( p=>vdot(p,[...point,1])); var transpose = matrix=> ( matrix.reduce(($, row) => row.map((_, i) => [...($[i] || []), row[i]]), []) ); var multiply = (a,b,...rest)=>(!b)?a:multiply(a.map( (p=>transpose(b).map(q=>vdot(p,q)))),...rest); I look forward to the day when the sufficiently-smart JIT implements the optimal Matrix Multiplication from that. :-)
  - svantana 8y ago
    Thanks! I find it weird that performance-focused javascript projects don't have live benchmarks up so you can judge for yourself. Personally I've found that a semi-clever matrix multiplier in c++ compiled with emscripten can give about 1.5 gflops (singlethreaded on a 2.5GHz i7 macbook pro), and existing js libs like sushi are a lot slower.
megalodon 8y ago
For anyone working with linear algebra in Node.js, I wrote C++ bindings to BLAS a couple of years ago [0]. I also started working on bindings to LAPACK, but never had the time to finish them [1]. [0] https://github.com/mateogianolio/nblas https://github.com/mateogianolio/nblas [1] https://github.com/mateogianolio/nlapack https://github.com/mateogianolio/nlapack

gtycomb 8y ago

Somewhere in the conversation Les Valiant's paper "General Context-free recognition in less than Cubic time" is worth mentioning. This has informative and clear discussion on the complexity of matrix-multiplication in programming areas such as grammar and parsing.

svantana 8y ago

I'm a big fan of writing simple, portable code that is structured enough to let the compiler optimize it to perfection. If you do it right, a naive A' * B multiplier can be within 20% of BLAS performance, at least for non-huge matrices. Now, when you need that code on some embedded platform or for emscripten or whatever, it'll just work and be really tiny too.

PeCaN 8y ago
Unfortunately, matrix multiplication is the kind of thing where you often need it to be fast and scale to huge matrices. I like simple code as much as the next guy, but occasionally you do need to write fast code.
ummonk 8y ago
It's the huge matrices where performance matters, and a naive multiplier for that will use memory rather inefficiently.
- svantana 8y ago
  > It's the huge matrices where performance matters I wouldn't say that, just look at CGI where a lot of effort is spent making fast 4x4 matrix multiplies. But yeah, a general matrix mult needs to scale, but a lot of the time you know what size you will run and you can optimize for that.
ziotom78 8y ago
Except that, when you are writing scientific code to solve a large linear problem, like a weather model or asupernova simulation, the difference between some naïve code and a heavily optimized one can be ten hours against ten days...
- carlmr 8y ago
  If it's within 20% of BLAS performance it shouldn't be the difference between ten hours and ten days. That must be something else then.

nickodell 8y ago

>The first thing that we need to do to accelerate our program is to improve the memory utilization of matrix B. Why not put matrix B in column-major order?

costrouc 8y ago
Matrix multiplication involves moving across one matrix in column order and the other matrix in row order. So it turns out that both row or column ordering make no difference. I think that matrix multiplication is one of the best examples of a deceptivly simple problem. It shows how far ALL code is from peak performance. We can only strive to be within an order or two from peak performance.
- fdej 8y ago
  This trick does work. If the matrices are in row-major order, you transpose B in memory and then compute A * (B^T)^T. This multiplication reads both matrices in row order. However, while this does improve performance over the naive algorithm, it's still not as good as a tiling algorithm.
  - shaklee3 8y ago
    I've found where that causes problems is when on of the matrix dimensions is not a multiple of the cache line size. It's common on gpus to use more elements then there are in the dimension. Nvidia calls this the leading dimension, and it must be greater than or equal to the Matrix dimension. If this is the case, the transpose trick doesn't quite work anymore.
- hedora 8y ago
  I took a different lesson from learning to implement matrix multiply: you can get a small amount of code within a few percent of peak performance if you try. So, the trick is to keep the 99% of the code that does 10% of the work within an order of magnitude of peak performance, and then get the rest within a factor of two. That gets the system within a factor of three, which is still embarrassing, but that’s what profiling tools are for. :-)

whos3424 8y ago

please look at http://halide-lang.org/ http://halide-lang.org/ and especially https://www.youtube.com/watch?v=3uiEyEKji0M https://www.youtube.com/watch?v=3uiEyEKji0M they seem to open new ground in the area you are prusuing edit: if you know any other groups/ppl who practice a similar approach send them my way anyhow, it seems like the hand written assembly vs compiled code goes in to a new skirmish the age where the whole spectrum from generalcpu > asic becomes available to programmers all around bottle necks of parallelism, locality and redundancy seems to be the guiding lights hand written we'll be probably lead by the magnificent http://www.greenarraychips.com/home/documents/pub/AP001-MD5.html http://www.greenarraychips.com/home/documents/pub/AP001-MD5.... compiler generated by halide and the like and once we got the ball rolling and more public knowledge is available probably the locked down techniques of the major supercomputers will start to emerge ( i just assume pps from the supercomputers area have been dealing with this stuff for the last couple of decades ) can i spit some more, dear hacker news? focus, you gotta focus, forget about the compiler generating good code, it's all about the ability for us humans to explore the domain even in halide, they mention the point that you still need to sit down and code different approaches, but their tool makes it a breeze compared to exploring those different approaches using hand written assembly

matachuan 8y ago

tiling and vectorizing are well known.

deleted 8y ago

[deleted]

nautilus12 8y ago

How would this change if you were dealing with more sparse matrices? Im used to using spark mllib distributed to do matrix multiplication on large sparse coordinate matrices and it works fine for me but i dont think i need the level of optimality this approach would demand.

knlji 8y ago
Doesn't the multiplication of two sparse matrices increase fill-in? That cannot be a good idea in general.
- stephencanon 8y ago
  Absolutely, but if that’s what you need to compute, that’s what you need to compute. Sometimes it is.
hadibrais 8y ago
Sparse matrix multiplication is generally considered to be a different problem. Sparse matrices need to be represented differently (See https://en.wikipedia.org/wiki/Sparse_matrix https://en.wikipedia.org/wiki/Sparse_matrix for an introduction). There has been recently a lot of research on how to do sparse matrix multiplication efficiently on manycore architectures. See this 2018 paper for example: https://arxiv.org/pdf/1804.01698.pdf https://arxiv.org/pdf/1804.01698.pdf. Note that sparse matrix-vector multiplication is also a different problem.
- nautilus12 8y ago
  Thats fascinating, i wonder how much coordinate matrix multiplication in spark takes into account factors like choelesky decomposition. Noticed most of those libraries are written in C++ or fortran, is JVM based approaches generally discouraged in this realm? It fits well with our current stack based mostly around spark
  - joshuamorton 8y ago
    Not quite the same, but similar, see the flowcharts at the end of this [1]. [1]: https://www.mathworks.com/help/matlab/ref/mldivide.html https://www.mathworks.com/help/matlab/ref/mldivide.html

taeric 8y ago

I was really expecting this to just be "import library". That said, very well done article. Would be interesting to see this style expand to the more ambitious methods in https://en.wikipedia.org/wiki/Matrix_multiplication_algorithm#Sub-cubic_algorithms https://en.wikipedia.org/wiki/Matrix_multiplication_algorith...

vmarsy 8y ago

Writing the fastest matrix multiplication you can is pretty fun. It probably won't be as fast a BLAS or MKL, but if you get close to it, it's a rewarding experience. > see what else can be done I recommend reading this paper by Kazushige Goto After doing the easiest change of loop reordering for already significant perf improvements, the paper linked ("Anatomy on a high performance matrix multiplication") is really excellent. It will walk you through all the optimizations like tiling for reducing L1, L2, L3, and TLB cache misses, and leverage vectorization. Then for squeezing out even more perf I remember there's the things like loop prefetching, loop unrolling, which you would expect the compiler with the best optimizations flags even targeted for the native architecture to take care of automatically. But you realize that it's not necessarily true.

comicjk 8y ago
If your particular matrix problem has special properties that aren't covered by blas, you can wring out an apparently impossible improvement. I once found a matrix problem (periodic cubic splines) which turned out to have a solution in O(N) - the matrix columns were related such that it could be replaced by a vector.
- yonkshi 8y ago
  I can attest this too, I was recently working with some massive and sparse binary matrix, with only 0.1% elements have values, it was magnitudes faster to index and sum the elements than to do a dot product. Not even a Tesla GPU was close to the performance of a custom matrix operation we wrote.
  - repsilat 8y ago
    A lot of matrix routines will have sparse matrix support. Though often the line between wanting a sparse matrix and wanting some kind of adjacency-list graph structure can be pretty fine... Also, only tangentially related, but if you can change your problem formulation to make the matrix of interest sparse(r), that can be a huge win. Two things that have helped me in the past: - Rows or columns that have mostly the same (nonzero) coefficient throughout. You can normally do some simple substitution to turn them into zeros. - Rows or columns that are naturally very similar to each other. You can often replace `y` with `y-x` or replace `x[i]` with `x[i]-x[i+1]` to turn one (or `n-1`) columns "mostly empty."

tehsauce 8y ago

A really nice video with animations of various scheduling approaches for these types of linear algebra algorithms: https://youtu.be/3uiEyEKji0M https://youtu.be/3uiEyEKji0M

fibo 8y ago

I wrote a JavaScript implementation with the naive method, but it also checks that matrix are compatibile and generalize the underlying ring operation (mul and add) so you could use it to multiply complex matrices, or bigint matrices, or every other ring you define. http://g14n.info/matrix-multiplication/ http://g14n.info/matrix-multiplication/

aqme28 8y ago

I'm surprised to see no discussion of faster-than-naive algorithms like Strassen[1] or Coppersmith-Winograd[2], which very counterintuitively run faster than N^3. Note: Only helpful for specific or very large matrices. [1]: https://en.wikipedia.org/wiki/Strassen_algorithm https://en.wikipedia.org/wiki/Strassen_algorithm [2]: https://en.wikipedia.org/wiki/Coppersmith%E2%80%93Winograd_algorithm https://en.wikipedia.org/wiki/Coppersmith%E2%80%93Winograd_a...

nwallin 8y ago
The problem with those algorithms is that they define "work" as "how many multiplications do I have to do?" Addition and memory accesses are defined by the analysis to be free. And you end up doing 4-5 times as many floating point additions as the naive algorithm, and the memory accesses are decidedly non-cache friendly. And modern CPUs have two multipcation units but only one addition unit, so multiplications are actually faster than addition. Plus, there are ergonomic issues with requiring square powers-of-two matrices. Matrices have to be truly spectacularly large to benefit from Strassens, and applications of the asymptotically faster algorithms don't really exist. The matrices have to be so large that they're computationally too expensive regardless of the asymptotic speedup. You're talking about dispatching work units to a cluster at that point. Most importantly though, the discussion is... boring. It's just "yup, we do Strassens, then we quickly reach the point where we go back to the naive algorithm, and things are interesting again."
- repsilat 8y ago
  > Matrices have to be truly spectacularly large to benefit from Strassens I don't have personal experience with other mmult methods, but Wikipedia says Strassen's algorithm is at least beneficial at matrix sizes that fit in cache. I have had occasion to deal with matrices of large dimension (though usually sparse, and usually not square) and those numbers don't seem crazy at all. Are they that far off-base? (Wikipedia does say that Coppersmith-Winograd tends to never be a useful improvement.) In any case, perhaps a more friendly response might be to point out that Strassen's algorithm tends to "bail out" into regular old matrix multiplication for small `n` anyway, so this still helps provide a speedup.

fabrice_d 8y ago

Related: https://news.ycombinator.com/item?id=17058189 https://news.ycombinator.com/item?id=17058189

phlip9 8y ago

A similar article on how to optimize SGEMM on GPUs using OpenCL - https://cnugteren.github.io/tutorial/pages/page1.html https://cnugteren.github.io/tutorial/pages/page1.html

gnufx 8y ago

This thread is old, but for the sake of archives: BLIS actually tells you how to write a fast production large-matrix GEMM, and the papers linked from https://github.com/flame/blis https://github.com/flame/blis would be a better reference than the Goto and van de Geijn. For small matrices see the references from https://github.com/hfp/libxsmm https://github.com/hfp/libxsmm but you're not likely to be re-implementing that unless you're serious about a version on non-x86_64. See also https://news.ycombinator.com/item?id=17061947 https://news.ycombinator.com/item?id=17061947

gnufx 8y ago
I could also have referenced the original(?) take on this: https://github.com/flame/how-to-optimize-gemm https://github.com/flame/how-to-optimize-gemm