newbie question re linear algebra in Haskell

[Hal Daume III]

Hi Jeff,

> (1)  Create a wrapper for the standard BLAS library.  Since an optimized
> BLAS is available for various architectures with the same interface,
> this would provide a way of optimizing the library for a particular
> architecture.

I started working on this a few months back, but it's actually quite a
large endeavor.  Well, actually, BLAS isn't too bad.  It's LAPACK that's a
monster.  IIRC, I finished porting BLAS
(http://www.isi.edu/~hdaume/HBlas)...look at "documentation" to see the
layout of the libraries.

> (2) Encapsulate this (using a monad?) in a way that provides safe access
> to the imperative routines from Haskell.

On my list of things to do, but again, haven't touched it in a while.

> (3) Using this as a foundation, define a matrix algebra in a "smart" way
> which avoids redundant evaluations.  For instance, if A is an arbitrary
> invertible matrix, we know that
>         A - A = 0
>         A + 0 = A
>         A * Indentity = A
>         A * inverse A = Identity
> so if we write an expression like "X = A * inverse A + B - B" we should
> expect our library to compute the correct result, "Identity" without
> actually performing any inversion, multiplication or elementwise matrix
> addition.

This is a bit more complicated, since in general, knowing if A=A is just
as expensive as doing the addition/etc.  One way off the top of my head to
model it on top of the current HBlas library is to associate with each
matrix a unique ID and a boolean.  Then, when matrices are created from
scratch you generate a new ID and set the boolean to false.  If you invert
the matrix, you leave the ID the same, but flip the boolean.  The problem
is with something like "A + B - B" versus "B - B + A".  One of these must
actually do the addition, depending on how the associativity is
defined.  Otherwise, in the first case, if you do the addition inner-most,
you won't know "soon enough" that you don't really need to do it...

> (4) Create a Haskell module which exposes a pure functional interface,
> hiding the imperative subsystem from the user.

I don't think this will work.  I think you need the imperative interface,
otherwise you are going to have to resort to copying of the arrays to
maintain purity.  IMO this is the wrong thing to do.

> My questions are:  (a) Is this idea sound? (b) Has someone already done
> this?  (c) Is there a better way of doing efficient linear algebra in
> Haskell?  I would greatly appreciate any advice in this regard.

I hope I answered some questions.