newbie question re linear algebra in Haskell
Hal Daume III
hdaume@ISI.EDU
Mon, 18 Nov 2002 13:39:05 -0800 (PST)
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.