Re: [Haskell-cafe] matrix computations based on the GSL

13 Jul 2005

      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

Re: [Haskell-cafe] matrix computations based on the GSL

Alberto Ruiz