Purely funcional LU decomposition

Rafael Gustavo da Cunha Pereira Pinto

3 Feb 2009 3 Feb '09

4:37 p.m.

Hello folks After a discussion on whether is possible to compile hmatrix in Windows, I decided to go crazy and do a LU decomposition entirely in Haskell... At first I thought it would be necessary to use a mutable or monadic version of Array, but then I figured out it is a purely interactive process. I am releasing this code fragment as LGPL. {- This program is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program. If not, see http://www.gnu.org/licenses/. -} import Data.Array.IArray type Dim=(Int,Int) lu::Array Dim Double -> (Array Dim Double,Array Dim Double) lu a = (aa l,aa u) where (l,u)=lu' a [] [] aa = accumArray (+) 0 (bounds a) lu'::(Floating e) => Array Dim e -> [(Dim,e)] -> [(Dim,e)] -> ([(Dim,e)],[(Dim,e)]) lu' a l u=if (ui==li) then ( ((ui,uj),1.0):l,((ui,uj),a!(ui,uj)):u) else (lu' an (l++ln) (u++un)) where k=li ((li,lj),(ui,uj))=bounds a lik i=(a!(i,k)/a!(k,k)) un=[((k,j),a!(k,j)) | j<-[lj..uj]] ln=((lj,lj),1.0):[((i,k),lik i) | i<-[li+1..ui]] an=array ((li+1,lj+1),(ui,uj)) [((i,j),e_an i j) | i<-[li+1..ui],j<-[lj+1..uj]] e_an i j=a!(i,j)-(lik i)*a!(k,j) Still needed: 1) Partial pivoting 2) Profiling... Lots! 3) Parallelize -- Rafael Gustavo da Cunha Pereira Pinto Electronic Engineer, MSc.

Attachments:

attachment.html (text/html — 6.6 KB)

Show replies by date

Paulo Tanimoto

3 Feb 3 Feb

8:26 p.m.

Hi Rafael, 2009/2/3 Rafael Gustavo da Cunha Pereira Pinto :

...

Hello folks

After a discussion on whether is possible to compile hmatrix in Windows, I decided to go crazy and do a LU decomposition entirely in Haskell...

At first I thought it would be necessary to use a mutable or monadic version of Array, but then I figured out it is a purely interactive process.

I am releasing this code fragment as LGPL.

Pretty cool, thanks for releasing this into the wild. I remember looking into this about a year ago. By the way, have you seen Matt's DSP library? http://haskelldsp.sourceforge.net/ He's got LU and others in there, if my memory serves me. The last release seems to be 2003, so it might be worth emailing him to see what happened and if he has plans for the future. Regards, Paulo

Rafael Gustavo da Cunha Pereira Pinto

4 Feb 4 Feb

9:49 a.m.

Matt's code is pretty comprehensive. His LU implementation is much cleaner, essentially because he used the "ijk" format, while I used the "kij". I'll take a look and e-mail him eventually. Thanks! On Tue, Feb 3, 2009 at 23:26, Paulo Tanimoto wrote:

...

Hi Rafael,

2009/2/3 Rafael Gustavo da Cunha Pereira Pinto
...
:

Hello folks

After a discussion on whether is possible to compile hmatrix in Windows, I decided to go crazy and do a LU decomposition entirely in Haskell...

At first I thought it would be necessary to use a mutable or monadic version of Array, but then I figured out it is a purely interactive process.

I am releasing this code fragment as LGPL.

Pretty cool, thanks for releasing this into the wild. I remember looking into this about a year ago. By the way, have you seen Matt's DSP library?

http://haskelldsp.sourceforge.net/

He's got LU and others in there, if my memory serves me. The last release seems to be 2003, so it might be worth emailing him to see what happened and if he has plans for the future.

Regards,

Paulo _______________________________________________ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe

-- Rafael Gustavo da Cunha Pereira Pinto Electronic Engineer, MSc.

Thomas Davie

1:42 a.m.

On 3 Feb 2009, at 22:37, Rafael Gustavo da Cunha Pereira Pinto wrote:

...

Hello folks

After a discussion on whether is possible to compile hmatrix in Windows, I decided to go crazy and do a LU decomposition entirely in Haskell...

At first I thought it would be necessary to use a mutable or monadic version of Array, but then I figured out it is a purely interactive process.

I am releasing this code fragment as LGPL.

Shinyness indeed – a quick note though, as ghc doesn't support dynamic linking of Haskell code, the above is equivalent to the GPL. Would be lovely if you packaged this up and stuck it on Hackage :) Bob

Rafael Gustavo da Cunha Pereira Pinto

4:41 a.m.

On Wed, Feb 4, 2009 at 04:42, Thomas Davie wrote:

...

Shinyness indeed – a quick note though, as ghc doesn't support dynamic linking of Haskell code, the above is equivalent to the GPL.

I always use LGPL. Anyway, I will keep it that way, as I still have hopes on dynamic linking Haskell becoming a reality

...

Would be lovely if you packaged this up and stuck it on Hackage :)

I may do that when I am finished with all the other stuff. I am building a (very primitive) circuit simulator using Haskell. It is a toy I am doing just to learn Haskell (until now I learned how to use Arrays and Parsec)...

...

Bob

-- Rafael Gustavo da Cunha Pereira Pinto Electronic Engineer, MSc.

Neal Alexander

2:15 a.m.

Array is no good man! Quad Tree matrices perform much nicer from what I've seen. I wrote some matrix stuff based on D. Stott Parker's "Randomized Gaussian elimination" papers (http://www.cs.ucla.edu/~stott/ge/). He presents some recursive block based methods of solving without pivoting. I want to upload the full library to Hackage, but it'd be nice to have some people look through it before - mainly because i never took linear algebra heh. QuickCheck seems to be holding things together though. Check these definitions (where the input matrix has been pre-transformed with the butterfly if needed -- revertRBT . invert . applyRBT): gauss (M.Scalar v) = (M.Scalar 1, M.Scalar v) gauss (M.Matrix (a,b,c,d)) = (M.Matrix (l, M.Scalar 0, x, l'), M.Matrix (u, y, M.Scalar 0, u')) where (l,u) = gauss a (l',u') = gauss $ d - (x * y) x = c * (invert u) y = (invert l) * b invert (M.Scalar v) = M.Scalar (1 / v) invert (M.Matrix (a,b,c,d)) = M.Matrix (w,x,y,z) where w = aI - (aIb * y) x = negate $ aIb * z y = negate $ z * caI z = invert $ d - (caI * b) ----------- aI = invert a caI = c * aI aIb = aI * b benchmarked against (in order of speed): 1. Mesa GL math in C called through FFI. 2. Finger tree based port. 3. Mutable array based port. 2 and 3 were scrapped right away. Heres some of the benchmark vs the C version. ------------------------------------ Matrix Inversion - 999999999x ------------------------------------ ------- haskell ------------ ContextSwitches - 224017 First level fills = 0 Second level fills = 0 ETime( 0:00:41.328 ) UTime( 0:00:40.875 ) KTime( 0:00:00.000 ) ITime( 0:00:41.484 ) ---------- C --------------- ContextSwitches - 640016 First level fills = 0 Second level fills = 0 ETime( 0:01:56.109 ) UTime( 0:01:55.359 ) KTime( 0:00:00.000 ) ITime( 0:01:54.328 ) ------------------------------------ Matrix Inversion - 9999999x ------------------------------------ ------- haskell ------------ ContextSwitches - 6204 First level fills = 0 Second level fills = 0 ETime( 0:00:00.687 ) UTime( 0:00:00.421 ) KTime( 0:00:00.000 ) ITime( 0:00:01.031 ) ---------- C --------------- ContextSwitches - 7396 First level fills = 0 Second level fills = 0 ETime( 0:00:01.171 ) UTime( 0:00:01.171 ) KTime( 0:00:00.000 ) ITime( 0:00:01.250 ) ------------------------------------- Matrix multiplication - 999999999x ------------------------------------- ------- haskell ------------ ContextSwitches - 234644 First level fills = 0 Second level fills = 0 ETime( 0:00:42.093 ) UTime( 0:00:41.609 ) KTime( 0:00:00.000 ) ITime( 0:00:42.187 ) ---------- C --------------- ContextSwitches - 228596 First level fills = 0 Second level fills = 0 ETime( 0:00:41.609 ) UTime( 0:00:41.375 ) KTime( 0:00:00.000 ) ITime( 0:00:41.515 )

Michael D. Adams

5 Feb 5 Feb

9:44 a.m.

On Wed, Feb 4, 2009 at 2:15 AM, Neal Alexander wrote:

...

Array is no good man! Quad Tree matrices perform much nicer from what I've seen.

I've co-authored a paper on that topic ("Seven at one stroke: results from a cache-oblivious paradigm for scalable matrix algorithms" http://portal.acm.org/citation.cfm?id=1178597.1178604). The current implementation is C based (to get base-case performance) but we would love to get that work ported to Haskell and a much more functional setting. Michael D. Adams

Rafael Gustavo da Cunha Pereira Pinto

10:39 a.m.

On Wed, Feb 4, 2009 at 05:15, Neal Alexander wrote:

...

Array is no good man! Quad Tree matrices perform much nicer from what I've seen.

I wrote some matrix stuff based on D. Stott Parker's "Randomized Gaussian elimination" papers (http://www.cs.ucla.edu/~stott/ge/http://www.cs.ucla.edu/%7Estott/ge/). He presents some recursive block based methods of solving without pivoting.

I want to upload the full library to Hackage, but it'd be nice to have some people look through it before - mainly because i never took linear algebra heh. QuickCheck seems to be holding things together though.

I read some of those papers and, yes, it is impressive. OTOH, I really need pivoting (some ill-conditioned matrices are expected), which they all claim to be hard to implement using quadtrees and I'm afraid that using RBTs might introduce other sources of error.

Dan Piponi

4 Feb 4 Feb

2:09 p.m.

2009/2/3 Rafael Gustavo da Cunha Pereira Pinto :

...

After a discussion on whether is possible to compile hmatrix in Windows, I decided to go crazy and do a LU decomposition entirely in Haskell... import Data.Array.IArray

...

e_an i j=a!(i,j)-(lik i)*a!(k,j)

There are three different representations of arrays in this code: arrays, lists and a functional one where f represents an array with elements (f i j). Looks like a candidate for a stream fusion type thing so the user only writes code using one representation and some RULES switch back and forth between representations behind the scenes. -- Dan

Rafael Gustavo da Cunha Pereira Pinto

3:57 p.m.

On Wed, Feb 4, 2009 at 17:09, Dan Piponi wrote:

...

2009/2/3 Rafael Gustavo da Cunha Pereira Pinto
...
: After a discussion on whether is possible to compile hmatrix in Windows, I decided to go crazy and do a LU decomposition entirely in Haskell... import Data.Array.IArray ... e_an i j=a!(i,j)-(lik i)*a!(k,j)

There are three different representations of arrays in this code: arrays, lists and a functional one where f represents an array with elements (f i j). Looks like a candidate for a stream fusion type thing so the user only writes code using one representation and some RULES switch back and forth between representations behind the scenes. -- Dan

Those different representations are derived from my (very) low Haskell handicap! :-D 1) I used lists for the intermediate L and U matrices so I wouldn't have the overhead of calling accumArray and assocs everywhere 2) I used the function that returns an array element just to shorten the line length, and to make sure I was typing it right! When I started this, I thought I would end with a couple of folds, but then I turned to the recursive version, so I could implement pivoting on matrix a' more easily. List comprehension also derived from this. I started thinking of a fold for L matrix and then it hit me that it was a list comprehension... Matt's DSP library has one of the most elegant implementations I saw, but it is harder to implement pivoting, which is really important for my application. lu :: Array (Int,Int) Double -- ^ A -> Array (Int,Int) Double -- ^ LU(A) lu a = a' where a' = array bnds [ ((i,j), luij i j) | (i,j) <- range bnds ] luij i j | i>j = (a!(i,j) - sum [ a'!(i,k) * a'!(k,j) | k <- [1 ..(j-1)] ]) / a'!(j,j) | i<=j = a!(i,j) - sum [ a'!(i,k) * a'!(k,j) | k <- [1 ..(i-1)] ] bnds = bounds a -- Rafael Gustavo da Cunha Pereira Pinto Electronic Engineer, MSc.

Dan Piponi

4:14 p.m.

On Wed, Feb 4, 2009 at 12:57 PM, Rafael Gustavo da Cunha Pereira Pinto wrote:

...

Those different representations are derived from my (very) low Haskell handicap! :-D

BTW That wasn't intended in any way as a criticism. I think it's interesting to look at real matrix code that people have written and think about what would be needed in a library to make it easier to write. -- Dan

Rafael Gustavo da Cunha Pereira Pinto

4:31 p.m.

On Wed, Feb 4, 2009 at 19:14, Dan Piponi wrote:

...

On Wed, Feb 4, 2009 at 12:57 PM, Rafael Gustavo da Cunha Pereira Pinto wrote:

...
Those different representations are derived from my (very) low Haskell handicap! :-D

BTW That wasn't intended in any way as a criticism. I think it's

Not taken!

...

interesting to look at real matrix code that people have written and think about what would be needed in a library to make it easier to write. -- Dan

What I miss most is a data structure with O(1) (amortized) direct access. One of the changes I thought today was to remove the ++ operator and create a list of lists that I would concat in the last call. -- Rafael Gustavo da Cunha Pereira Pinto Electronic Engineer, MSc.

wren ng thornton

5 Feb 5 Feb

8:46 p.m.

Rafael Gustavo da Cunha Pereira Pinto wrote:

...

What I miss most is a data structure with O(1) (amortized) direct access.

Finger trees get close, O(log(min(i,n-i))): http://hackage.haskell.org/packages/archive/containers/latest/doc/html/Data-... (And Theta(1) amortized for all dequeue operations to boot.)

...

One of the changes I thought today was to remove the ++ operator and create a list of lists that I would concat in the last call.

Don't do it manually, let an abstract type do the optimization for you. This idea[1] is often known as "difference lists", which the estimable Don Stewart has already packaged up for us: http://hackage.haskell.org/cgi-bin/hackage-scripts/package/dlist :) [1] Actually, not quite that idea. The idea instead is to represent a list as a function [a]->[a] and then compose these functions (rather than concatenating the strings). At the very end once you're done, pass in [] and get back the whole list in O(n) time, rather than the O(n^2) that commonly arises from using (++). -- Live well, ~wren

5949

Age (days ago)

5952

Last active (days ago)

List overview

Download

12 comments

7 participants

participants (7)

Dan Piponi
Michael D. Adams
Neal Alexander
Paulo Tanimoto
Rafael Gustavo da Cunha Pereira Pinto
Thomas Davie
wren ng thornton