# Use of Strassen's algorithm

Strassen's algorithm is an algorithm for matrix multiplication that is asymptotically faster than the naive one. In practice, the matrices involved have to be quite large before Strassen's algorithm becomes faster than the naive one.

## Question

Does Mathematica ever use Strassen's algorithm? If so, what is the size of the smallest (square) matrices on which it would do so?

• I believe that for numeric matrices beyond a certain threshold size, Mathematica outsources matrix operations to LAPACK/BLAS subroutines, so I think your question is logically equivalent to the question "Does LAPACK ever use Strassen's algorithm?". I don't know much about how LAPACK works internally, though, so other people can probably answer that. Commented Aug 1, 2014 at 15:55
• Reference: From the documentation "tutorial/SomeNotesOnInternalImplementation", section "Approximate Numerical Linear Algebra". Commented Aug 1, 2014 at 15:57

The answer is probably far from what you expect (other end of the spectrum, so to speak).

As noted in a comment, for numerical linear algebra Mathematica, at some level, uses library BLAS. I believe this does not use asymptotically fast matrix products for two reasons. One is that those methods are not able, as best I recall, to take advantage of data locality in the way that highly optimized level 3 BLAS does using traditional multiplication of matrices. The other is that the numerical stability of various algorithms is worsened by fast multiplication methods. (I think this may also be the case for ffts vs dfts, despite the former involving far fewer operations; has to do with reuse of correlated error I think).

Important caveat: Either or both of these reasons may have become invalid since last I had read anything on this topic. And regardless of what I wrote, level 3 BLAS might or might not be making use of fast multiplication.

There is at least one place where Strassen (7-for-8) multiplication is in fact used. It is when using 2x2 matrices that are needed for asymptotically fast integer gcd and continued fractions implementation code. Strassen in this case does not itself make the algorithm asymptotically faster, but since the bulk of cost is in large integer multiplications one gets a 1/8 speed boost on those parts of the algorithms that use 2x2 matrix products. If memory serves me correctly this gives an overall gain factor somewhere close to 1/8 (maybe 10% or so).

• Paraphrase quote: "For numerical linear algebra, Mathematica uses BLAS and I don't think BLAS uses Strassen's algorithm." Who/where should I ask to confirm this? Commented Aug 3, 2014 at 17:57
• @TysonWilliams: BLAS does use Strassen in DGEMMS*. But Daniel's description of why to avoid that algorithm for floating-point calculations is apt. Commented Aug 5, 2014 at 14:41
• @Charles That is good to know, but my original question remains. For what dimension matrices does BLAS use DGEMMS*? On page 108 of Numerical Recipes: The Art of Scientific Computing (3rd ed.), which was published in 2007, it says that BLAS uses Strassen's algorithm when the size of the matries is about 100. Daniel's answer (which he qualifies with "I believe") is in contradiction with this quote. Commented Aug 6, 2014 at 11:59
• The link mathworks.com/matlabcentral/newsreader/view_thread/242624 might be useful. Especially the response by Tim Davis (author of TAUCS). Commented Aug 6, 2014 at 17:46
• Also I would imagine that whether Strassen or Winograd or similar is used will depend on what BLAS implementation is under the hood, and how it is invoked by LAPACK. Commented Aug 6, 2014 at 17:48