[Haskell-cafe] matrix computations based on the GSL

[Alberto Ruiz]

Hello! I have a few doubts concerning the LinearAlgebra library...

On Friday 08 July 2005 11:29, Henning Thielemann wrote:
>I would remove the 'matrix' portions of the function names, since this
>information is already contained in the module name. E.g. after importing
>LinearAlgebra.Matrix as Matrix, Matrix.matrix look strange, but
>Matrix.fromList says everything.

I have changed the function names as suggested. This new style is clearly 
better, allowing Vector.add, Matrix.add, Vector.Complex.add, 
Matrix.Complex.add, etc.

>If the Matrix type would be parametrised (...)

Now we can have Vector.T a and Matrix.T a for any storable a (although at this 
point most functions are only defined for Double). For example:

eig :: Matrix.T Double -> ([Double], Matrix.T Double)

>I would also prefer a vector of complex numbers for the FFT function
>instead of a nx2 matrix.

This comes from the old version in which vectors and matrices were the same 
type :-) Now we have:
 
fft1D  :: Vector.T (Complex Double) -> Vector.T (Complex Double)
ifft1D :: Vector.T (Complex Double) -> Vector.T (Complex Double)
fft2D :: Matrix (Complex Double) -> Matrix (Complex Double)
etc.


But now I have two problems:

1) If I define

foo :: Vector.T a -> Matrix.T a

Haddock (version 0.6) shows just this:

foo :: T a -> T a
 
I don't know how to tell haddock that I want the full names in the signatures.


2) When setting up the "main" module which reexports all the modules exposed 
to the user I cannot write this:

---------------------------------
module LinearAlgebra ( 
    module Vector,
    module Matrix,
    module Matrix.Complex,
    module LinearAlgebra.Algorithms
    --etc.
    ) 
where
import LinearAlgebra.Vector as Vector
import LinearAlgebra.Matrix as Matrix
import LinearAlgebra.Matrix.Complex as Matrix.Complex
import LinearAlgebra.Algorithms
--etc.
-------------------------------------

since we get lots of conflicting exports for T, add, scale, etc. 

I we use qualified imports:

-----------------------------------------------------------
module LinearAlgebra ( 
    module Vector,
    module Matrix,
    module Matrix.Complex,
    module LinearAlgebra.Algorithms
    --etc.
    ) 
where
import qualified LinearAlgebra.Vector as Vector
import qualified LinearAlgebra.Matrix as Matrix
import qualified LinearAlgebra.Matrix.Complex as Matrix.Complex
import LinearAlgebra.Algorithms
-----------------------------------------------------------

then all the names in Vector, Matrix and Matrix.Complex must be always 
qualified, which is not very practical. I am afraid that I am doing something 
wrong... 

Another possibility is using type classes, since there is a lot of types 
addmiting add, sub, (ew)mul, scale, etc. (e.g., Vector Int, Vector Double, 
Matrix (Complex Float), etc.) We could also define type classes for higher 
level linear algebra algorithms working with different base types in the 
Matrices.

Any suggestion?     
    

And some final comments.

>If there are no efficiency concerns, I would drop element-wise operations
>and prefer a matrix-map and a matrix-zipWith. If these operations shall
>remain I would somehow point to their element-wise operation in the name.

There is about 5x speed gain if we map in the C side. The "optimized" floating 
map functions could be moved to a separate module.
 
>If the Matrix type would be parametrised then Matrix.fromBlocks could use
>a more natural indexing.

>Matrix.fromBlocks :: [[Matrix a b]] -> Matrix (Int,a) (Int,b)

>The Int of the index pair would store the block number and the type a
>index would store the index within the block. Hm, but it would not be
>possible to store the sizes of the sub-matrices.  It would be possible
>only if the sub-matrices have the same size.  Maybe it is better to allow
>matrices of matrices and to be able to apply a matrix-matrix
>multiplication which in turn employs a matrix-matrix multiplication of its
>elements. But the matrix decompositions will fail on this structure - but
>possibly those which fail aren't appropriate for block matrices.

We can have a Matrix (Matrix Double) if we define the instance Storable 
(Matrix a). Although a "block matrix" is perhaps better represented at a 
higher level (storing only one copy of repeated blocks, etc.). 

>It would be nice to have functions which don't simply write the formatted
>vectors and matrices to stdout but which return them as strings. E.g.

>Matrix.format :: Int -> Matrix.T -> String
>Vector.format :: Int -> Vector.T -> String
>
>(or Matrix.toString)
>
>The results can be reused in other contexts and the 'display' routines can
>use their results.

Done.

--Alberto