Matrix library in Haskell

23 Apr 2002

      Hello,

        I am just discovering Haskell, so sorry if this is not the
right place to ask. I want to use it for some numerical
calculations. I need something higher level than C++ and faster than
Python or Matlab and from the initial experiments it seems that
Haskell could be the right choice. My question is: Is there any
matrix/vector library available with common operations (dot products,
matrix products, linear system solution etc.) ? I could not find any.

I am including her my first try on such a library, in case it might be
useful for somebody. However, I am not perfectly happy with the
design. In particular I would like to define MatrixClass and
VectorClass so that applying a getRow operation on a matrix instance
would yield automatically the correct instance of the VectorClass, but
I do not know how to express this interdependency, so for the moment I
have dropped the classes.

Yours,

Jan

-- 
-------------------------------------------------------------------------
Jan Kybic <kybic@ieee.org>      Robotvis, INRIA, Sophia-Antipolis, France
       or <Jan.Kybic@sophia.inria.fr>,tel. work +33 492 38 7589, fax 7845
                    http://www-sop.inria.fr/robotvis/personnel/Jan.Kybic/

-- Module implementing vectors, matrices, and operations on them
-- it uses the multiparameter class extension 
-- 
-- $Id: LinearAlgebra.hs,v 1.4 2002/04/22 11:44:41 jkybic Exp $
-- Jan Kybic, April 2002

module LinearAlgebra where

import Array
import Complex
import Observe

--------------------------------------------------------------------
-- now define the concrete vector and array implementations
-- TODO: This could be accelerated using IArray/ UArray
-- TODO: Make it a proper instance of Show

newtype ArrayVectorType a = ArrayVector (Array Int a) deriving  Show
newtype ArrayMatrixType a = ArrayMatrix (Array (Int,Int) a) deriving  Show

listToVec l = ArrayVector (array (0,uplim) $ zip [0..uplim] l) where
	  		    uplim = (length l) -1

vecToList (ArrayVector v) = [ v!i | i <- indices v  ]

dot (ArrayVector a) (ArrayVector b) = sum [ a!i * b!i | i <- indices a ]

norm2 a = dot a a 

norm a = sqrt (norm2 a)  

(!@) (ArrayVector v) i = v!i
sizeV (ArrayVector v) = rangeSize $ bounds  v
indicesV v= [0..(sizeV v - 1)]

scaleV (ArrayVector a) s = ArrayVector ( fmap (\x -> s * x) a )     

(+@) a b = listToVec [ a!@i + b!@i | i <- indicesV a ]

(-@) a b = listToVec [ a!@i - b!@i | i <- indicesV a ]

-- matrix operations

(!@@) (ArrayMatrix m) (i,j) = m!(i,j)	    

listToMat l = 
	  let { bnds=((0,0),(m,n)) ; m=length l -1 ; n= length (head l) -1} 
	    in ArrayMatrix ( array bnds 
			     [ ((i,j),(l!!i)!!j) | (i,j) <- range bnds ] )
matToList a = 
    [ [ a !@@ (i,j) | j <- range $ cbounds a ] | i <- range $ rbounds a ]

boundsM (ArrayMatrix a) = bounds a
rbounds m = let z=boundsM m in (fst(fst z),fst(snd z))
cbounds m = let z=boundsM m in (snd(fst z),snd(snd z))
nrows = rangeSize . rbounds 
ncols = rangeSize . cbounds
getRow a i = -- extract row i as vector
      listToVec [ a!@@(i,j) | j <- range $ cbounds a ] 
getCol a j = -- extract column j as vector
      listToVec [ a!@@(i,j) | i <- range $ cbounds a ] 
showMat a = "[ " ++ concat 
	       [ showVec (getRow a i) ++ "  " | 
				i <- range (rbounds a) ] ++ "]"
scaleMat (ArrayMatrix a) s = ArrayMatrix ( fmap (\x -> s * x) a )     
matMult a b =
     let bnds = transTup (rbounds a,cbounds b) in
     ArrayMatrix (array bnds [ ((i,j), (getRow a i) `dot` getCol b j) | 
		 (i,j) <- range bnds ])
idMat n = -- create an identity matrix 
       let bnds=((0,0),(n-1,n-1)) in
		ArrayMatrix (accumArray (+) 0.0 bnds 
		     [ ((i,i),1.0) | i <- [0..n-1] ])

-- cross takes two vectors of length three and calculates their cross product
-- uncomment the run-time checks if you prefer safety over speed
-- the signature is a little limiting to avoid uncertainity of the resulting
-- vector

cross a b 
--    | or [ sizeV a /=3 , sizeV b /=3 ] = 
--	error "Cross product needs length 3 vectors"
--    | otherwise 
                = listToVec [ a !@ 1 * b !@ 2 - a !@ 2 * b !@ 1,
	                    - a !@ 0 * b !@ 2 + a !@ 2 * b !@ 0,
                              a !@ 0 * b !@ 1 - a !@ 1 * b !@ 0 ]

vangle a b = acos (cosvangle a b) -- angle between two vectors

cosvangle' a b = -- the cos of an angle between two vectors
    let  mag= sqrt ( norm2 a * norm2 b )
    in if mag>0.0 then (dot a b)/mag
       else 0.0

-- the above version had numerical problems
cosvangle a b = max ( min (cosvangle' a b) 1.0 ) (-1.0) 

showVec a = -- converts a vector to a string
    "[ " ++ concat [ show (a !@ i) ++ " " | i <- [0..sizeV a-1] ]
	 ++ "]"

subvector a (i,j) = -- another vector with indices i..j
         listToVec [ a !@ k | k <- [i..j] ]

transTup ((a,b),(c,d)) = ((a,c),(b,d)) -- transpose a tuple of tuples

-- type synonymes

type DoubleVector = ArrayVectorType Double
type DoubleMatrix = ArrayMatrixType Double

-- Observable instance

instance (Observable a) => Observable (ArrayVectorType a) where
	 observer (ArrayVector a) = send "ArrayVector" 
	                            (return ArrayVector << a)

instance (Observable a) => Observable (ArrayMatrixType a) where
	 observer (ArrayMatrix a) = send "ArrayMatrix" 
	                            (return ArrayMatrix << a)

-- Functor instance

instance Functor ArrayVectorType where
	 fmap f (ArrayVector a) = ArrayVector ( fmap f a )

-- end of Linear Algebra.hs