# 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. –  DumpsterDoofus Aug 1 '14 at 15:55
Reference: From the documentation "tutorial/SomeNotesOnInternalImplementation", section "Approximate Numerical Linear Algebra". –  DumpsterDoofus Aug 1 '14 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? –  Tyson Williams Aug 3 '14 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. –  Charles Aug 5 '14 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. –  Tyson Williams Aug 6 '14 at 11:59
I am looking for a definitive answer...not "I believe" or "typically used for $N \ge 100$". Someone familiar with the BLAS source code should be able to tell us a specific scenario, which preferably is based solely on a threshold of the matrices' dimensions, when DGEMMS* is called. –  Tyson Williams Aug 6 '14 at 12:04
The link mathworks.com/matlabcentral/newsreader/view_thread/242624 might be useful. Especially the response by Tim Davis (author of TAUCS). –  Daniel Lichtblau Aug 6 '14 at 17:46