4 ms·

It's typically still the O(N^3) method that's implemented in OpenBLAS, BLIS, cuBLAS, MKL, etc. There are some high performance implementations of Strassen's Alg

It's typically still the O(N^3) method that's implemented in OpenBLAS, BLIS, cuBLAS, MKL, etc. There are some high performance implementations of Strassen's Algorithm with a fixed level of recursion that start to perform much better as the matrix size gets large (see the work from UT Austin's BLIS). From my understanding for the more complex algorithms the hidden constant simply grows too large to be worth it on conventional hardware. As another poster commented, these are O(N^3) in work, not necessarily the wall clock time you see due to parallel speedup and cache effects, but will scale as O(N^3) once you're at a size past all of these. These implementations are highly tuned and optimized for hardware, much much more so than the naive 3 loop implementation. The predictable and 'easy' access patterns make the simple algorithm hard to beat.

kkylin 6y ago
Another issue of both practical and theoretical importance is whether these asymptoticaly-faster algorithms are numerically stable. Skimming the article very quickly (and not having looked at any original sources), this doesn't seem to be addressed (yet).
- jmblpati 6y ago
  Strassen is reasonably numerically stable (although not as good as the naive algorithm), and everything beyond Strassen is extremely unstable. This is due to a trick that many of these algorithms employ: you technically only need to do multiplication of 0-1 matrices to fully solve the matrix multiplication problem in theory, and so having an algorithm which computes the entries of the matrix with additive error 0.1 is sufficient to exactly solve the problem (just round entries to the nearest integer). As you can imagine, this means your algorithm can only give O(1) bits of precision unless you employ this reduction to the 0-1 case first. To understand why this even happens, let's say we want to compute the expression A*B + B*A for matrices A and B. One way we can do this is to compute the products A*B and B*A naively: two multiplications are needed. A trickier way to do this is to introduce a parameter x: we will instead compute A*A and (x*A + B/x)*(x*A + B/x) = x^2 A*A + (A*B + B*A) + x^-2 B*B. Thus (x*A + B/x)*(x*A + B/x) - x^2 A*A = (A*B + B*A) + x^-2 B*B: if x is large enough the second term vanishes and we can employ the rounding trick from before. Now this still needed two multiplications, but here one of our multiplications was A*A. If later we needed to compute A*C + C*A in the algorithm, we could then do that in only 1 additional matrix multiplication by repeating the trick. A more sophisticated version of this algorithm underlies all known approaches for matrix multiplication beyond w << 2.8.
  - kkylin 6y ago
    Thank you for the very nice explanation!
- oscardssmith 6y ago
  At least for strassen, the result is slightly less stable (increases the condition number by a constant factor). I think higham shows that any sub-cubic algorithm must have done conditioning issues, but that in practice they're not that bad for most matrices.
taxemeEvasion 6y ago
I think its also important to mention that these results are for the multiplication of two general dense matrices and full floating point accuracy. If your matrix is structured (sparse, sparse in the Fourier domain, low rank, H, HSS, etc) you can usually exploit that structure to break up the multiplication and reduce the work complexity dramatically. (These rely on building blocks of the dense general matrix results)

kkylin 6y ago

Another issue of both practical and theoretical importance is whether these asymptoticaly-faster algorithms are numerically stable. Skimming the article very quickly (and not having looked at any original sources), this doesn't seem to be addressed (yet).

jmblpati 6y ago
Strassen is reasonably numerically stable (although not as good as the naive algorithm), and everything beyond Strassen is extremely unstable. This is due to a trick that many of these algorithms employ: you technically only need to do multiplication of 0-1 matrices to fully solve the matrix multiplication problem in theory, and so having an algorithm which computes the entries of the matrix with additive error 0.1 is sufficient to exactly solve the problem (just round entries to the nearest integer). As you can imagine, this means your algorithm can only give O(1) bits of precision unless you employ this reduction to the 0-1 case first. To understand why this even happens, let's say we want to compute the expression A*B + B*A for matrices A and B. One way we can do this is to compute the products A*B and B*A naively: two multiplications are needed. A trickier way to do this is to introduce a parameter x: we will instead compute A*A and (x*A + B/x)*(x*A + B/x) = x^2 A*A + (A*B + B*A) + x^-2 B*B. Thus (x*A + B/x)*(x*A + B/x) - x^2 A*A = (A*B + B*A) + x^-2 B*B: if x is large enough the second term vanishes and we can employ the rounding trick from before. Now this still needed two multiplications, but here one of our multiplications was A*A. If later we needed to compute A*C + C*A in the algorithm, we could then do that in only 1 additional matrix multiplication by repeating the trick. A more sophisticated version of this algorithm underlies all known approaches for matrix multiplication beyond w << 2.8.
- kkylin 6y ago
  Thank you for the very nice explanation!
oscardssmith 6y ago
At least for strassen, the result is slightly less stable (increases the condition number by a constant factor). I think higham shows that any sub-cubic algorithm must have done conditioning issues, but that in practice they're not that bad for most matrices.

taxemeEvasion 6y ago

I think its also important to mention that these results are for the multiplication of two general dense matrices and full floating point accuracy. If your matrix is structured (sparse, sparse in the Fourier domain, low rank, H, HSS, etc) you can usually exploit that structure to break up the multiplication and reduce the work complexity dramatically. (These rely on building blocks of the dense general matrix results)