Monday, October 18, 2010

Fast matrix multiplication in R: Strassen's algorithm

I tried to implement the Strassen's algorithm for big matrices multiplication in R.

Here I present a pdf with some theory element, some example and a possible solution in R.
I'm not a programmer, so the function is not optimize, but it works.

I want to thank G. Grothendieck: suggested me a very nice way on StackOverFlow to create a bigger square matrix starting from small one.

This is just a first version of the function; it needs more work on it. If someone want to collaborate, I'll be very happy.
Finally if you find my code useful for your work, I'd love to be cited (ask me via e-mail how to cite me: todoslogos -at- gmail . com).


  1. Nice post.

    One modification you could make that would get rid of the need for strassen2(), strassen3(), etc. would be to check the sizes of your submatrices, and choose your multiplication function accordingly. %*% if small, strassen if big.

    Also, the times in several of your benchmarks are so small as to probably be meaningless, I think you should probably give them at least 10 seconds' worth of computation to chew on if you want to make a solid case.

    Finally, I'd be surprised if the underlying implementation of %*% wasn't already pretty well-optimized using some spiffy linear algebra libraries. Or maybe that's what the Matrix package is for, I dunno.

  2. Thanks Kenahoo for your nice comment!
    And also for your suggestion, to integrate all function into only one. Maybe I will add the upadated code, under the pdf, asap.
    About the timing computation, I'm looking for better method to assess differences between function. For example, would be great to have a function that return the number of computation done, or the number of multiplication/addition, to make a comparison between the two function.
    And yes, I was surprised too by noting that %*% is not the fastest choice. As far as I know there isn't any implementation of function for fast matrix multiplication, in R. This is the reason why I wrote this post :)

  3. I'm trying to learn the Strassen algorithm right now, so I may well be wrong, but shouldn't M1 in your example be
    29972 28340
    16254 16243

    instead of

    29972 16254
    28340 16243

    i.e., rotated 90 degrees?

    128 148
    77 78

    48 137
    161 73

    So (A11+A22)(B11+B22) should equal:
    29972 28340
    16254 16243

    I'm using an online calculator to get this result.