Overloaded matrix operations

9 Nov 2000

      Unclassified

LS,

here an haskell file that implements overloaded matrix operations. Such
overloading is not found
in Jan Skibinski's Orthogonals.lhs program. I would like to make this
code public domain
and add this file to the "hslibs", does anyone know how to go about it?

Greetings Frank Mulder

-- overloaded matrix operations
--
-- by F.W. Mulder
-- fwmulder@signaal.nl
-- Hollandse Signaalapparaten B.V.
--
module Matrix (
t,Mat(M,V,VT,S),Mtx((./),(.*),(.~),det,fromScalar,fromVec,
  fromMat,vecLength,matDims,matShow,solve,lu,ldl,zeros,
  ones,diag,trace,distL,isCovMat),
  SM,sm_find, sm_mk_spm, sm_mk_mat, sm_size,sm_rows,sm_cols)
  where

import List
import Maybe
-------------------- indexless linear algebra dense matrices
------------------
t = -1000

infixr 9 .~

data Mat = M [[Double]] | V [Double] | VT [Double] | S Double
  deriving (Eq,Show,Read)
class Mtx a where
 (./) :: a -> Double -> a
 (.*) :: a -> Double -> a
 (.~) :: a -> Integer -> a
 det :: a -> Double
 fromScalar :: a -> Double
 fromVec :: a -> [Double]
 fromMat :: a -> [[Double]]
 vecLength :: a -> Int
 matDims :: a -> (Int,Int)
 matShow :: a -> String
 solve :: (a -> (a,a)) -> a -> a -> a
 lu :: a -> (a,a)
 ldl :: a -> (a,a)
 zeros :: Int -> a
 ones :: Int -> a
 diag :: a -> a
 trace :: a -> Double
 distL :: a -> a -> Double
 isCovMat :: a -> Bool
instance Mtx Mat where
 (./) (M a) b = M [[e / b | e <- r] | r <- a]
 (./) (V a) b = V [e / b | e <- a]
 (.*) (M a) b = M [[e * b | e <- r] | r <- a]
 (.*) (V a) b = V [e * b | e <- a]
 (.~) (M a) n  | n == -1 = M (inv' a)
       | n == t  = M (transpose a)
   | otherwise = error "Matrix:inv"
 (.~) (V a) n  | n == t  = VT a
   | otherwise = error "Matrix:inv"
 (.~) (VT a) n  | n == t  = V a
   | otherwise = error "Matrix:inv"
 det (M a) = det' a
 det _ = error "Matrix:det"
 fromScalar (S a) = a
 fromVec (V a) = a
 fromVec (VT a) = a
 fromMat (M a) = a
 vecLength (V x) = length x
 matDims (M a) = (length a, length (head a))
 matDims _ = error "Matrix:matDims"
 matShow (M b) = "\n\n" ++ concat [show r ++ "\n" | r <- b] ++ ";"
 solve fac m = solve' fac m
 lu m = lu' m
 zeros n = V (take n (repeat 0.0))
 ones n = V (take n (repeat 1.0))
 diag (M a) = V (diag' a)
 trace (M a) = sum (diag' a)
 distL a b = distL' a b
 isCovMat a = isCovMat' a
instance Num a => Num [a] where
 (+) a b = zipWith (+) a b
 (-) a b = zipWith (-) a b
instance Num Mat where
 (+) (M a) (M b) = M (a + b)
 (+) (V a) (V b) = V (a + b)
 (+) (VT a) (VT b) = VT (a + b)
 (-) (M a) (M b) = M (a - b)
 (-) (V a) (V b) = V (a - b)
 (-) (VT a) (VT b) = VT (a - b)
 (*) (M a) (M b) =  M (a `xmm` b)
 (*) (M a) (V b) =  V (a `xmv` b)
 (*) (VT a) (V b) = S (inprod a b)
 (*) (V a) (VT b) = M (outprod a b)
 (*) (VT a) (M b) = VT (a `xvm` b)

diag' [r] = [head r]
diag' ([]:_) = error "Matrix:diag'"
diag' (r:rs) = head r : diag' [tail r | r <- rs]
dims' :: [[Double]] -> (Int,Int)
dims' [r] = (length r, 0)
dims' (r:rs) = (length r, length rs + 1)
isCovMat' b = n == m && all (\e -> e > 0.0) (fromVec (diag b))
  where (n,m) = matDims b
------------------- standard determinand, inverse
-----------------------------
det' [[x]]   = x
det' a@(v:_) = inprod v [cofactor a 0 j | j <- [0..n-1]]
            where n = length v

cofactor a i j = fromIntegral ((-1)^(i+j)) * det' (minor a i j)

minor a i j = [omit j v | v <- omit i a]

omit 0 (x:xs) = xs
omit _ []     = []
omit n (x:xs) = (x:omit (n-1) xs)

inv' [[x]] = [[1.0/x]]
inv' a     = [[(cofactor a j i) / d | j <- [0..n]] | i <- [0..n]]
          where n = length a - 1
                d = det' a
----------------------- vector products
---------------------------------------
inprod [] _  = 0.0
inprod (x:xs) (y:ys) = x * y + inprod xs ys
outprod a b = [[e * f | f <- b] | e <- a]
----------------------- matrix products
---------------------------------------
xmm a b = [[inprod r c | c <- transpose b] | r <- a]
xmv a b = [inprod r b | r <- a]
xvm a b = [inprod a c | c <- b]
----------------------- LU decomposition
--------------------------------------
lu' (M a) = factLU a [head a] []
factLU [[a]] u l = (M (transpose (l ++ [[1.0]])), M u)
factLU a u l = factLU a' (u ++ [head a']) (l ++ [[1.0] ++ cm])
  where (c,cs) = unzip [(head r, tail r) | r <- a]
   cm = [e / (head c) | e <- (tail c)]
   rm = head cs
 a' = tail cs --- \\ outprod cm rm

----------------------- choleski decomposition
--------------------------------
ldl' (M (r:rs)) = (M rt, M (transpose rt))
  where rt = ldl2 1 ([],r,rs)
 ldl2 j (u,r,[]) = [factCH r u j 1 []]
 ldl2 j (u,r,l)  = (ro : ldl2 (j+1) (u++[ro],head l, tail l))
  where ro = factCH r u j 1 []

factCH [] _ _ _ _ = []
factCH (a:as) u j i lj | i > j     = []
         | i == j    = (el1 : factCH as u j (i+1) (lj++[el1]))
         | otherwise = (el2 : factCH as u j (i+1) (lj++[el2]))
 where
  el1 = sqrt (a - inprod lj lj)
  el2 = (a - inprod li lj) / lii
  (li,(lii:_)) = splitAt (i-1) (u!!(i-1))

-- forward and backward substitution for lower and upper triang.
matrices -----
forwL l b = subst_ l b []
backU u b = reverse (subst_ (map reverse (reverse u)) (reverse b) [])

subst_ [] _ s = s
subst_ (l:ls) (b:bs) s = subst_ ls bs (s ++ [xi])
  where xi = (b  - inprod (init l) s) / (last l)

solve' fac a (V b) = let (M l,M u) = fac a in V (backU u (forwL l b))

distL' (V x) (M r) = inprod x (backU (transpose r) (forwL r x))

-- Matrix.hs
-- overloaded matrix operations
--
-- by F.W. Mulder
-- fwmulder@signaal.nl
-- Hollandse Signaalapparaten B.V.
--
module Matrix (	t,Mat(M,V,VT,S),Mtx((./),(.*),(.~),det,fromScalar,fromVec,
		fromMat,vecLength,matDims,matShow,solve,lu,ldl,zeros,
		ones,diag,trace,distL,isCovMat),
		SM,sm_find, sm_mk_spm, sm_mk_mat, sm_size,sm_rows,sm_cols) 
		where

import List
import Maybe
-------------------- indexless linear algebra dense matrices ------------------
t = -1000

infixr 9 .~

data Mat = M [[Double]] | V [Double] | VT [Double] | S Double
		deriving (Eq,Show,Read)
class Mtx a where
	(./) :: a -> Double -> a
	(.*) :: a -> Double -> a
	(.~) :: a -> Integer -> a
	det :: a -> Double
	fromScalar :: a -> Double
	fromVec :: a -> [Double]
	fromMat :: a -> [[Double]]
	vecLength :: a -> Int
	matDims :: a -> (Int,Int)
	matShow :: a -> String
	solve :: (a -> (a,a)) -> a -> a -> a
	lu :: a -> (a,a)
	ldl :: a -> (a,a)
	zeros :: Int -> a
	ones :: Int -> a
	diag :: a -> a
	trace :: a -> Double
	distL :: a -> a -> Double
	isCovMat :: a -> Bool
instance Mtx Mat where
	(./) (M a) b = M [[e / b | e <- r] | r <- a]
	(./) (V a) b = V [e / b | e <- a]
	(.*) (M a) b = M [[e * b | e <- r] | r <- a]
	(.*) (V a) b = V [e * b | e <- a]
	(.~) (M a) n 	| n == -1 = M (inv' a)
		    	| n == t  = M (transpose a)
			| otherwise = error "Matrix:inv"
	(.~) (V a) n 	| n == t  = VT a
			| otherwise = error "Matrix:inv"
	(.~) (VT a) n 	| n == t  = V a
			| otherwise = error "Matrix:inv"
	det (M a) = det' a
	det _ = error "Matrix:det"
	fromScalar (S a) = a
	fromVec (V a) = a
	fromVec (VT a) = a
	fromMat (M a) = a
	vecLength (V x) = length x
	matDims (M a) = (length a, length (head a))
	matDims _ = error "Matrix:matDims"
	matShow (M b) = "\n\n" ++ concat [show r ++ "\n" | r <- b] ++ ";"
	solve fac m = solve' fac m
	lu m = lu' m
	zeros n = V (take n (repeat 0.0))
	ones n = V (take n (repeat 1.0))
	diag (M a) = V (diag' a)
	trace (M a) = sum (diag' a)
	distL a b = distL' a b
	isCovMat a = isCovMat' a
instance Num a => Num [a] where
	(+) a b = zipWith (+) a b
	(-) a b = zipWith (-) a b
instance Num Mat where
	(+) (M a) (M b) = M (a + b)
	(+) (V a) (V b) = V (a + b)
	(+) (VT a) (VT b) = VT (a + b)
	(-) (M a) (M b) = M (a - b)
	(-) (V a) (V b) = V (a - b)
	(-) (VT a) (VT b) = VT (a - b)
	(*) (M a) (M b) =  M (a `xmm` b)
	(*) (M a) (V b) =  V (a `xmv` b)
	(*) (VT a) (V b) = S (inprod a b)
	(*) (V a) (VT b) = M (outprod a b)
	(*) (VT a) (M b) = VT (a `xvm` b)

diag' [r] = [head r]
diag' ([]:_) = error "Matrix:diag'"
diag' (r:rs) = head r : diag' [tail r | r <- rs]
dims' :: [[Double]] -> (Int,Int)
dims' [r] = (length r, 0)
dims' (r:rs) = (length r, length rs + 1)
isCovMat' b = n == m && all (\e -> e > 0.0) (fromVec (diag b))
  where (n,m) = matDims b
------------------- standard determinand, inverse -----------------------------
det' [[x]]   = x
det' a@(v:_) = inprod v [cofactor a 0 j | j <- [0..n-1]]
            where n = length v

cofactor a i j = fromIntegral ((-1)^(i+j)) * det' (minor a i j)

minor a i j = [omit j v | v <- omit i a]

omit 0 (x:xs) = xs
omit _ []     = []
omit n (x:xs) = (x:omit (n-1) xs)

inv' [[x]] = [[1.0/x]]
inv' a     = [[(cofactor a j i) / d | j <- [0..n]] | i <- [0..n]]
          where n = length a - 1
                d = det' a
----------------------- vector products ---------------------------------------
inprod [] _  = 0.0
inprod (x:xs) (y:ys) = x * y + inprod xs ys
outprod a b = [[e * f | f <- b] | e <- a]
----------------------- matrix products ---------------------------------------
xmm a b = [[inprod r c | c <- transpose b] | r <- a]
xmv a b = [inprod r b | r <- a]
xvm a b = [inprod a c | c <- b]
----------------------- LU decomposition --------------------------------------
lu' (M a) = factLU a [head a] []
factLU [[a]] u l = (M (transpose (l ++ [[1.0]])), M u)
factLU a u l = factLU a' (u ++ [head a']) (l ++ [[1.0] ++ cm])
  where	(c,cs) = unzip [(head r, tail r) | r <- a]
  	cm = [e / (head c) | e <- (tail c)]
  	rm = head cs
	a' = tail cs --- \\ outprod cm rm

----------------------- choleski decomposition --------------------------------
ldl' (M (r:rs)) = (M rt, M (transpose rt))
  where	rt = ldl2 1 ([],r,rs)
	ldl2 j (u,r,[]) = [factCH r u j 1 []] 
	ldl2 j (u,r,l)  = (ro : ldl2 (j+1) (u++[ro],head l, tail l))
		where	ro = factCH r u j 1 []

factCH [] _ _ _ _ = []
factCH (a:as) u j i lj | i > j     = []
		       | i == j    = (el1 : factCH as u j (i+1) (lj++[el1]))
		       | otherwise = (el2 : factCH as u j (i+1) (lj++[el2]))
	where	
		el1 = sqrt (a - inprod lj lj)
		el2 = (a - inprod li lj) / lii
		(li,(lii:_)) = splitAt (i-1) (u!!(i-1))

-- forward and backward substitution for lower and upper triang. matrices -----
forwL l b = subst_ l b []
backU u b = reverse (subst_ (map reverse (reverse u)) (reverse b) [])

subst_ [] _ s = s
subst_ (l:ls) (b:bs) s = subst_ ls bs (s ++ [xi])
  where	xi = (b  - inprod (init l) s) / (last l)

solve' fac a (V b) = let (M l,M u) = fac a in V (backU u (forwL l b))

distL' (V x) (M r) = inprod x (backU (transpose r) (forwL r x))
---------------------- sparse matrices ----------------------------------------
type SM = ((Int,Int),Double)

sm_size [[]] = (0,0)
sm_size (r:rs) = (length rs + 1, length r)

sm_find ix sm = if (cc /= Nothing) then (Just (ix,fromJust cc)) else Nothing
  where	cc = lookup ix sm 
sm_mk_spm b = ([((i,j),e) | (i,r) <- zip [1..n] b, (j,e) <- zip [1..m] r],n,m)
  where (n,m) = sm_size b
sm_mk_mat (sm,n,m) large = 
	[[mkc (sm_find (i,j) sm) | j <- [1..m]] | i <- [1..n]]
  where mkc Nothing = large
	mkc (Just (_,c)) = c
sm_rows [] = []
sm_rows [e] = [[e]]
sm_rows (e@((i1,_),_):es) = (e : rw) : sm_rows rws
  where	(rw,rws) = partition (\((i,_),_) -> i == i1) es
sm_cols [] = []
sm_cols [e] = [[e]]
sm_cols (e@((_,j1),_):es) = (e : cl) : sm_cols cls
  where	(cl,cls) = partition (\((_,j),_) -> j == j1) es

Overloaded matrix operations

f.w. mulder