newbie question re linear algebra in Haskell
I'm a complete newcomer to Haskell, having learned about only recently. I'm intrigued by the possibility of in using it for numerical applications, specifically linear algebra. I understand that (at least in its present state) Haskell 98 isn't competitive with imperative languages when it comes to primitive matrix-vector operations, which often rely on destructive updating. It strikes me that one approach that takes advantage of the strengths of both paradigms would be create an imperative subsystem to handle primitive operations, then create a functional matrix algebra layer on top of it. Based on my very limited understanding of Haskell, I envision a linear algebra library which would work something like this: (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. (2) Encapsulate this (using a monad?) in a way that provides safe access to the imperative routines from Haskell. (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. (4) Create a Haskell module which exposes a pure functional interface, hiding the imperative subsystem from the user. 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. --Jeff Austin
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.
I'm a complete newcomer to Haskell, having learned about only recently. I'm intrigued by the possibility of in using it for numerical applications, specifically linear algebra. I understand that (at least in its present state) Haskell 98 isn't competitive with imperative languages when it comes to primitive matrix-vector operations, which often rely on destructive updating. It strikes me that one approach that takes advantage of the strengths of both paradigms would be create an imperative subsystem to handle primitive operations, then create a functional matrix algebra layer on top of it. [..]
One thing that comes to mind is Barry Jay's FISh language:
http://www-staff.it.uts.edu.au/~cbj/Publications/shapes.html#Array_Programmi...
This compiles code in a functional language with arrays down to C, by using "shape inference" to fix the size of all the arrays.
I believe FFTW (the Fastest Fourier Transform in the West) similarly uses a functional programming language to generate imperative (C) code.
--KW 8-)
--
Keith Wansbrough
--- Keith Wansbrough
I'm a complete newcomer to Haskell, having learned about only recently. I'm intrigued by the possibility of in using it for numerical applications, specifically linear algebra. I understand that (at least in its present state) Haskell 98 isn't competitive with imperative languages when it comes to primitive matrix-vector operations, which often rely on destructive updating. It strikes me that one approach that takes advantage of the strengths of both paradigms would be create an imperative subsystem to handle primitive operations, then create a functional matrix algebra layer on top of it. [..]
One thing that comes to mind is Barry Jay's FISh language:
http://www-staff.it.uts.edu.au/~cbj/Publications/shapes.html#Array_Programmi...
This compiles code in a functional language with arrays down to C, by using "shape inference" to fix the size of all the arrays.
I believe FFTW (the Fastest Fourier Transform in the West) similarly uses a functional programming language to generate imperative (C) code.
Mr. Austin migfht also want to look at some dated modules at ftp://ftp.cs.chalmers.se/pub/haskell/library/bevan/ and perhaps at Matrix Inversion using Quadtrees Implemented in Gofer (1995) Jeremy D. Frens & David S. Wise http://citeseer.nj.nec.com/frens95matrix.html Auto-Blocking Matrix-Multiplication or Tracking BLAS3 Performance from Source Code (1997) Jeremy D. Frens & David S. Wise http://citeseer.nj.nec.com/frens97autoblocking.html and
From Fast Exponentiation to Square Matrices: An Adventure in Types Chris Okasaki http://www.eecs.usma.edu/Personnel/okasaki/pubs.html#icfp99
Chris ===== Christopher Milton cmiltonperl@yahoo.com __________________________________________________ Do you Yahoo!? Yahoo! Web Hosting - Let the expert host your site http://webhosting.yahoo.com
participants (4)
-
Christopher Milton -
Hal Daume III -
Jeff Austin -
Keith Wansbrough