Matrix library in Haskell

Jan Kybic kybic@ieee.org
23 Apr 2002 13:33:43 +0200

```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

```