Emil Hedevang | 11 Dec 16:14 2012
Picon

Large, numerical calculations in Haskell

Hi Haskell Cafe,

I need to perform very large numerical computations requiring tens of GB of memory. The computations consist essentially of generation of random numbers and discrete convolutions of large arrays of random numbers with somewhat smaller kernel arrays. Should I use Haskell and call some C libraries when necessary? Should I just write everything in C? Or should I use e.g. Chapel (chapel.cray.com)?

To be more specific, a very large array X is filled with random numbers drawn independently from some distribution. Then, the discrete convolution Y = K * X of X with a somewhat smaller kernel array K is calculated.

Typical sizes of X and K could be 1000 x 1000 x 1000 and 100 x 100 x 100, respectively. The storage of X alone requires 8 GB of memory. Luckily I have recently been blessed with a well-equipped workstation with 128 GB RAM and two 8-core Xeon E5-2650 running at 2.0 GHz.

I see immediate opportunities for parallelism. Firstly, parallel random number generators could be employed to populate the array X (e.g. SPRNG). Secondly, since the size of the kernel K is much smaller than the size of the array X, the discrete convolution can split into several smaller convolutions (see e.g. Kim and Strintzis (1980)). Finally, the discrete convolutions can be computed efficiently using Fourier methods (using e.g. FFTW).

It seems fairly straightforward to implement it all in C, perhaps with some OpenMP pragmas scattered around in the code, since the heavy numerics would be performed in preexisting libraries (SPRNG and FFTW),

I would like to use Haskell, since I like Haskell and find it tremendously fascinating, but I will not use it for the sake of just using it. I am willing to pay a small price in terms of performance, say, twice the run time, if I in turn will get more convenience and less suffering while coding. Therefore, would you claim it to be feasible to use Haskell? I suppose it very much depends on the quality of the Haskell FFI bindings to parallel random number generators and fast Fourier transforms. Unfortunately, I find it difficult to assess these qualities.

Best regards,
Emil Hedevang
_______________________________________________
Haskell-Cafe mailing list
Haskell-Cafe <at> haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe
Dominic Steinitz | 12 Dec 10:15 2012

Re: Large, numerical calculations in Haskell

Emil Hedevang <emilhedevang <at> gmail.com> writes:

> 
> 
> Hi Haskell Cafe,
> 
> I need to perform very large numerical computations requiring tens of GB of 
memory. The computations consist essentially of generation of random numbers 
and discrete convolutions of large arrays of random numbers with somewhat 
smaller kernel arrays. Should I use Haskell and call some C libraries when 
necessary? Should I just write everything in C? Or should I use e.g. Chapel 
(chapel.cray.com)?
> 
> 
Hi Emil,

The first place I would look would be repa http://repa.ouroborus.net/. IIRC it 
supports discrete convolutions and repa folks seem quite responsive.

Dominic.

PS I am planning to use repa to solve PDEs and address, at least 
partially, "the curse of dimensionality" with it.
Emil Hedevang | 12 Dec 11:03 2012
Picon

Re: Large, numerical calculations in Haskell

Hi A.M. and Dominic,


I have been following the development of Repa for some time, mostly reading blogs and papers about it, not actually writing any code. It's some quite cool stuff the developers are doing.

A discrete convolution is a stencil computation. There are (to my knowledge) essentially two approaches to calculate a discrete convolutions Y = K * X:
(1) from the definition, and
(2) using discrete Fourier transforms

Approach no. 1 runs in O(nK nX) where nK and nX are the number of elements in the arrays K and X, respectively, whereas approach no. 2 runs in O(nX log nX). If nK is very, very small, then approach no. 1 is fastest, otherwise approach no. 2 is fastest. If nK is somewhat smaller than nX, then some tricks can be applied to obtain something a bit faster than approach no. 2.

To my knowledge, Repa applies approach no. 1 in its stencil computations and discrete convolutions. Therefore it cannot be applied in my case, since I have nK = 10^6 and nX = 10^9 (or something similar), and so I must use approach no. 2 with discrete Fourier transforms (DFT). In Repa, the DFT for arbitrary array sizes is implemented using the slow O(n^2) algorithm, i.e., implemented using the definition of the DFT. If the size of the array is a power of two, then the DFT is implemented using the fast Fourier transform (FFT). Unfortunately, the Repa documentation states that the Repa FFT is approximately 50 times slower than the FFT from FFTW. (And the FFT from FFTW is even capable of calculating fast DFTs of arbitrary size). This essentially forces my to use FFTW for the actual number crunching.

I suppose it would be a very interesting research project to implement FFTW entirely in Haskell. Unfortunately, I currently do not have the time to embark on such a quest.

Regarding the generation of random numbers, Repa implements what they call a minimal standard Lehmer generator. I am not sure if that random number generator is of high enough qualify for my needs.

Cheers,
Emil


On 12 December 2012 10:15, Dominic Steinitz <dominic <at> steinitz.org> wrote:
Emil Hedevang <emilhedevang <at> gmail.com> writes:

>
>
> Hi Haskell Cafe,
>
> I need to perform very large numerical computations requiring tens of GB of
memory. The computations consist essentially of generation of random numbers
and discrete convolutions of large arrays of random numbers with somewhat
smaller kernel arrays. Should I use Haskell and call some C libraries when
necessary? Should I just write everything in C? Or should I use e.g. Chapel
(chapel.cray.com)?
>
>
Hi Emil,

The first place I would look would be repa http://repa.ouroborus.net/. IIRC it
supports discrete convolutions and repa folks seem quite responsive.

Dominic.

PS I am planning to use repa to solve PDEs and address, at least
partially, "the curse of dimensionality" with it.




_______________________________________________
Haskell-Cafe mailing list
Haskell-Cafe <at> haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe

_______________________________________________
Haskell-Cafe mailing list
Haskell-Cafe <at> haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe
Carter Schonwald | 12 Dec 20:36 2012
Picon

Re: Large, numerical calculations in Haskell

Emil, there is an actively maintained set of FFTW bindings on hackage, would that help you?


cheers
-Carter


On Wed, Dec 12, 2012 at 5:03 AM, Emil Hedevang <emilhedevang <at> gmail.com> wrote:
Hi A.M. and Dominic,

I have been following the development of Repa for some time, mostly reading blogs and papers about it, not actually writing any code. It's some quite cool stuff the developers are doing.

A discrete convolution is a stencil computation. There are (to my knowledge) essentially two approaches to calculate a discrete convolutions Y = K * X:
(1) from the definition, and
(2) using discrete Fourier transforms

Approach no. 1 runs in O(nK nX) where nK and nX are the number of elements in the arrays K and X, respectively, whereas approach no. 2 runs in O(nX log nX). If nK is very, very small, then approach no. 1 is fastest, otherwise approach no. 2 is fastest. If nK is somewhat smaller than nX, then some tricks can be applied to obtain something a bit faster than approach no. 2.

To my knowledge, Repa applies approach no. 1 in its stencil computations and discrete convolutions. Therefore it cannot be applied in my case, since I have nK = 10^6 and nX = 10^9 (or something similar), and so I must use approach no. 2 with discrete Fourier transforms (DFT). In Repa, the DFT for arbitrary array sizes is implemented using the slow O(n^2) algorithm, i.e., implemented using the definition of the DFT. If the size of the array is a power of two, then the DFT is implemented using the fast Fourier transform (FFT). Unfortunately, the Repa documentation states that the Repa FFT is approximately 50 times slower than the FFT from FFTW. (And the FFT from FFTW is even capable of calculating fast DFTs of arbitrary size). This essentially forces my to use FFTW for the actual number crunching.

I suppose it would be a very interesting research project to implement FFTW entirely in Haskell. Unfortunately, I currently do not have the time to embark on such a quest.

Regarding the generation of random numbers, Repa implements what they call a minimal standard Lehmer generator. I am not sure if that random number generator is of high enough qualify for my needs.

Cheers,
Emil


On 12 December 2012 10:15, Dominic Steinitz <dominic <at> steinitz.org> wrote:
Emil Hedevang <emilhedevang <at> gmail.com> writes:

>
>
> Hi Haskell Cafe,
>
> I need to perform very large numerical computations requiring tens of GB of
memory. The computations consist essentially of generation of random numbers
and discrete convolutions of large arrays of random numbers with somewhat
smaller kernel arrays. Should I use Haskell and call some C libraries when
necessary? Should I just write everything in C? Or should I use e.g. Chapel
(chapel.cray.com)?
>
>
Hi Emil,

The first place I would look would be repa http://repa.ouroborus.net/. IIRC it
supports discrete convolutions and repa folks seem quite responsive.

Dominic.

PS I am planning to use repa to solve PDEs and address, at least
partially, "the curse of dimensionality" with it.




_______________________________________________
Haskell-Cafe mailing list
Haskell-Cafe <at> haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe


_______________________________________________
Haskell-Cafe mailing list
Haskell-Cafe <at> haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe


_______________________________________________
Haskell-Cafe mailing list
Haskell-Cafe <at> haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe

Gmane