Arrays vs Lists and Specialization

[Matthew Donadio]

Hi all,

I'm sorry if this topic has been rehashed a lot, but I poked around in
the mailing list archive, and didn't find what I was looking for.

I currently have some free time on my hands, and have been implementing
some digital signal processing and spectral/frequency estimation
algorithms along with the needed matrix routines in Haskell.  

For those unfamiliar with this field, most algorthms are defined in
textbooks in terms of indexing through discrete sequences.  For example
the implementation of cross-correlation

> rxy x y k | k >= 0 = sum [ x!(i+k) * (conjugate (y!i)) | i <- [0..(n-1-k)] ]
>           | k <  0 = conjugate (rxy y x (-k))
>                      where n = snd (bounds x) + 1

looks very similar to the textbook definition if one uses arrays.  A
definition using lists would probably use drop, map, and zipWith, and
look nothing like the definitions found in the standard texts.

Many spectral estimation routines are defined in terms of special
matrices (ie, Toeplitz, etc).  Arrays defined recursively by list
comprehensions make it easy to implement algorithms like Levinson-Durbin
recursion, and they look very similar to the mathematical definitions:

> levinson r p = (array (1,p) [ (k, a!(p,k)) | k <- [1..p] ], realPart (rho!p))
>     where a   = array ((1,1),(p,p)) [ ((k,i), ak k i) | k <- [1..p], i <- [1..k] ]
>	  rho = array (1,p) [ (k, rhok k) | k <- [1..p] ]
>	  ak 1 1             = -r!1 / r!0
>	  ak k i | k==i      = -(r!k + sum [ a!(k-1,l) * r!(k-l) | l <- [1..(k-1)] ]) / rho!(k-1)
>		 | otherwise = a!(k-1,i) + a!(k,k) * (conjugate (a!(k-1,k-i)))
>	  rhok 1 = (1 - (abs (a!(1,1)))^2) * r!0
>	  rhok k = (1 - (abs (a!(k,k)))^2) * rho!(k-1)

This could be rewritten for lists, but would probably need to be defined
in terms of an aux. recursive function, which destroys the simplicity of
the above definition.

OK, my question then has to do with the efficiency of lists versus
arrays.  Do the latest compilers handle handle arrays efficiently, or
are lists really the way to go?  If there is a performace difference, is
it generally big enough to warrant rewriting algorithms?

A related question is how is specilization handled in arrays with lazy
evaluation.  In the definition of levinson above, the a array is defined
in terms of the ak function.  By doing this, you save some horizontal
space, but it also unburdens the programmer from tracking the recursive
dependencies.  a!(k,k) is needed before a!(i,j) can be calculated, but
lazy evaluation takes care of this.  If the above function is
specialized for r::Array Int (Complex Double) and p::Int, would I be
correct to say that the final value of the function would be unboxed,
but all intermediate values wouldn't?  Now, in some cases, a user may
need all of the model orders from 1..p.  This is handled easilly enough
by just changing the first line. to 

> levinson r p = (a, fmap realPart rho)

Would the a matrix in the tuple be unboxed with specilization?

If anyone is interesting in what I have put together, I will be making
everything public sometime next week.  I have a lot of algorithms
implemented, but I need to clean up the documentation a bit (well, a
lot).

Thanks.

-- 
Matthew Donadio (m.p.donadio@ieee.org)