[Haskell-cafe] Haskell, C and Matrix Multiplication

[Blake Rain]

Dear Haskellers,

I thought I'd take some time to share something of my weekend with
you  all.   Not because  of  anything new,  but  because  it is  a
feel-good example  of Haskell  being win. A  good read  for Monday
morning, perhaps.

What  feels  like  decades  ago,   I  used  to  work  in  computer
graphics.   Mostly    my   work   was   on    software   for   the
amateur/semi-professional  film   industry  -  developing  various
filters and effects for the likes of Adobe AfterEffects and so on.

I  have   not  touched   computer  graphics  from   a  programming
perspective for quite a long time  now, and I thought it was about
time that  I flexed those  old muscles. So,  I have taken  it upon
myself to write an SVG renderer for Blender.

To cut what is going to be  a very long and humdrum story short, I
needed to write some matrix arithmetic code.

Having not done  this in a long  time, I thought I'd see  if I can
still actually remember what matrix multiplication is by writing a
function to multiply two 4x4 matrices.  I thought I'd write one in
C and then one in Haskell, and see how long they took and how easy
they were to write.

Here is the one in C:

void multiplyM4 (float m1[4][4], float m2[4][4], float m3[4][4]) {
m1[0][0] = m2[0][0] * m3[0][0] + m2[0][1] * m3[1][0] + m2[0][2] *
m3[2][0] + m2[0][3] * m3[3][0];
m1[0][1] = m2[0][0] * m3[0][1] + m2[0][1] * m3[1][1] + m2[0][2] *
m3[2][1] + m2[0][3] * m3[3][1];
m1[0][2] = m2[0][0] * m3[0][2] + m2[0][1] * m3[1][2] + m2[0][2] *
m3[2][2] + m2[0][3] * m3[3][2];
m1[0][3] = m2[0][0] * m3[0][3] + m2[0][1] * m3[1][3] + m2[0][2] *
m3[2][3] + m2[0][3] * m3[3][3];

m1[1][0] = m2[1][0] * m3[0][0] + m2[1][1] * m3[1][0] + m2[1][2] *
m3[2][0] + m2[1][3] * m3[3][0];
m1[1][1] = m2[1][0] * m3[0][1] + m2[1][1] * m3[1][1] + m2[1][2] *
m3[2][1] + m2[1][3] * m3[3][1];
m1[1][2] = m2[1][0] * m3[0][2] + m2[1][1] * m3[1][2] + m2[1][2] *
m3[2][2] + m2[1][3] * m3[3][2];
m1[1][3] = m2[1][0] * m3[0][3] + m2[1][1] * m3[1][3] + m2[1][2] *
m3[2][3] + m2[1][3] * m3[3][3];

m1[2][0] = m2[2][0] * m3[0][0] + m2[2][1] * m3[1][0] + m2[2][2] *
m3[2][0] + m2[2][3] * m3[3][0];
m1[2][1] = m2[2][0] * m3[0][1] + m2[2][1] * m3[1][1] + m2[2][2] *
m3[2][1] + m2[2][3] * m3[3][1];
m1[2][2] = m2[2][0] * m3[0][2] + m2[2][1] * m3[1][2] + m2[2][2] *
m3[2][2] + m2[2][3] * m3[3][2];
m1[2][3] = m2[2][0] * m3[0][3] + m2[2][1] * m3[1][3] + m2[2][2] *
m3[2][3] + m2[2][3] * m3[3][3];

m1[3][0] = m2[3][0] * m3[0][0] + m2[3][1] * m3[1][0] + m2[3][2] *
m3[2][0] + m2[3][3] * m3[3][0];
m1[3][1] = m2[3][0] * m3[0][1] + m2[3][1] * m3[1][1] + m2[3][2] *
m3[2][1] + m2[3][3] * m3[3][1];
m1[3][2] = m2[3][0] * m3[0][2] + m2[3][1] * m3[1][2] + m2[3][2] *
m3[2][2] + m2[3][3] * m3[3][2];
m1[3][3] = m2[3][0] * m3[0][3] + m2[3][1] * m3[1][3] + m2[3][2] *
m3[2][3] + m2[3][3] * m3[3][3];
}

Urgh! This took me over 20  minutes to write out, and a further 34
minutes  of  debugging.   That's  54  minutes to  write  a  simple
function!   This included 23  invocations of  the compiler  and 18
invocations of  the test program.  I was disgusted with  myself! I
should have automated the writing of it with Haskell!

The  reason for this  was mostly  that I  kept pressing  the wrong
buttons  when writing those  many, many  array indices.   I became
blind  to the  numbers in  the  end, and  it became  one of  those
nightmare  debugging  sessions where  you  check your  printMatrix
function and  other stupid things.  I'd  like to say that  I was a
bit tired, but that's just an excuse.

Some points about the C version:

  1. It's only 21 lines of code.

  2. It  took me 54 minutes  to write and debug.   Therefore I can
     write eight of these in a standard working day. Hmm.

  3. It's written in C,  so everybody and their dog understands it
     completely. Honest.

  4. It's written  in C, so it can just be  pasted into nearly any
     other   language   (because   they're  all   just   augmented
     C). Therefore no one needs to understand it.

  5. It does not consume much space (in terms of memory).

  6. It's pretty fast (being written in C):

     Performing  100,000,000 matrix  multiples  takes on  average
     6.117s (with gcc -O2).

  7. Passing a NULL pointer anywhere will cause it to explode.

  8. If the pointers are not to arrays of 16 floats, it will trash
     memory  (being  as  it  destructively  updates  one  of  it's
     arguments).

  9. It will only operate on 4x4 float matrices.

10. How am I storing matrices again...? Column major?


So, after drinking  some coffee and having a  smoke, I started the
Haskell version:

[foldl (+) 0 $ zipWith (*) x y | x <- m1, y <- transpose m2]

Ah.

So... anyway, some points about the Haskell version:

  1  It's only  one line  of code  - that's  20 less  than  the C
     version.

  2.  It took  me just  over five  minutes to  think about  and to
     code. That's a saving of  49 minutes (1180% or something). It
     compiled straight away and worked first time. I can write too
     many of these in a day.

  3. It's laughably clear and concise. Anyone can easily tell what
     I'm doing.

  4. It  does *not*  operate in constant  space. I think.  I can't
     tell. I don't know how to tell. Maybe it does.

  5. It's very fast (being written in Haskell):

     Performing  100,000,000  matrix  multiples takes  on  average
     1.435s  (with ghc  -O2).

  6. Passing in empty lists will not crash it.

  7. It will operate on square matrices of any size (3x3, 4x4...).

  8. It will operate on a matrix of any type that is of class Num.

The Haskell  version seems to win  on all counts. 

Interestingly  enough (to some at least), other  matrix operations
also become  trivial when expressed  in Haskell:

  Determinant is  a bit fiddly,  but quite trivial. However  it is
  woefully slow at  O(n!). I need to make  use of LU decomposition
  (in which the  determinant is the sum of  the diagonal values in
  U).

  Adjoint is even more of a fiddle, but actually a one-liner after
  reusing  bits   of  the  determinant   function  (including  the
  determinant function).

  Matrix  inversion is  easy in  one  line, but  dependent on  the
  determinant and the adjoint, therefore also woefully slow.

Just another example of Haskell being awesome :)

Now over to you guys! What experiences have you had which make you
glad you use Haskell? Anything make you smile recently?

- Blake