Scalar multiplication with a new Matrix type in Haskell

Question

I've been programming in C-type and Lisp-type languages for a few decades and Haskell for a few years. Now in order to further my understanding of typeclasses and other nifty more advanced features of Haskell, such as parallelism, I've been trying to create a new data type with matrix semantics, for now piggy-backed on top of the library's Array type.

I am using the latest Haskell Platform 2013.2.0.0 with the included ghc.

newtype Matrix a = Matrix (Array (Int,Int) a)

(I understand that data would work instead of newtype, but that in this instance, you get the same semantics with better performance using newtype instead).

Some simple functions to create an identity matrix, for example, work just fine:

unityMatrix:: (Num a) => Int -> Matrix a
unityMatrix d = Matrix $ let b = ((0,0),(d-1,d-1)) in array b
    [((i,j),if i==j then 1 else 0)|(i,j)<-range b]

So does creating a basic show function:

instance Show a => Show (Matrix a) where
    show (Matrix x) =
        let
            b = bounds x
            rb = ((fst.fst) b, (fst.snd) b)
            cb = ((snd.fst) b, (snd.snd) b)
        in
            intercalate "\n" [unwords ([show (x!(r,c))|c<-range cb])|r<-range rb]

Then, for the regular arithmetic operators to work, I add this:

instance Num a => Num (Matrix a) where
    fromInteger x =
        let
            b = ((0,0),(0,0))
        in Matrix $ array b
            [((i,j),if i == j then (fromInteger x) else 0)|(i,j)<-range b]

    (Matrix x) + (Matrix y) =
        let
            b = bounds x
        in
            if b /= bounds y then error "Unmatched matrix addition" else Matrix $ array b
                [(ij,x!ij + y!ij)|ij<-range b]

    signum (Matrix x) =
        let
            b = bounds x
        in Matrix $ array b
            [(ij,signum (x!ij))|ij<-range b]

    abs (Matrix x) =
        let
            b = bounds x
        in Matrix $ array b
            [(ij,abs (x!ij))|ij<-range b]

    (Matrix x) - (Matrix y) =
        let
            b = bounds x
        in
            if b /= bounds y then error "Unmatched matrix subtraction" else Matrix $ array b
                [(ij,x!ij - y!ij)|ij<-range b]

    (Matrix x) * (Matrix y) =
        let
            b = (((fst.fst.bounds) x, (fst.snd.bounds) x),((snd.fst.bounds) y, (snd.snd.bounds) y))
            kb = ((snd.fst.bounds) x, (snd.snd.bounds) x)
        in
            if kb /= ((fst.fst.bounds) y, (fst.snd.bounds) y) then error "Unmatched matrix multiplication" else  Matrix $ array b
                [((i,j),sum [(x!(i,k)) * (y!(k,j))|k<-range kb])|(i,j)<-range b]

(My apologies if the indentation is screwed up here--the actual code is properly indented and compiles.)

So far, so good, though it is a little annoying to have to define a fromInteger function which does not really have any meaning in matrix semantics, but creating a 1x1 matrix with the value is as reasonable as anything else.

My problem is trying to get proper semantics for multiplication of a scalar (i.e., a type of the Num typeclass) with a matrix. By mathematical convention that means just element-wise multiplication with the scalar.

However, no matter what syntax I try, I don't get the right result. By default this implementation just promotes, e.g., an Int to a Matrix using fromInteger which (unless the Matrix is already 1x1) results in a "Unmatched matrix multiplication" error. I have tried all alternative syntaxes I can think of to define a alternate code for this type of multiplication without success, such as:

(Matrix x) * (y::Int) =

or

(Matrix x) * (Num y) =

or

(*):: (Num a) => Matrix -> a -> Matrix

But all of these give me various syntax errors.

How should I go about defining the scalar by matrix multiplication so that it does what you would expect it to do? I feel that enabling the non-standard pattern guard feature might help, but I'm not quite sure how to use that correctly in this context.

I realize that I could just special case the Matrix multiplication to allow multiplication of any Matrix with a 1x1 Matrix which I guess would work. But that would be (a) inelegant, (b) un-Haskell-y, (c) probably inefficient as it would require every scalar to be wrapped into a matrix before being multiplied, and (d) would permit some code (such as multiplication of an arbitrary matrix with any 1x1 matrix) to run when it should result in a n error.

I also understand that there are probably excellent Matrix implementations out there which somehow sidestep this problem. But using them would defeat my purpose to learn.

Answer 1

This is the type of the standard (*) operator:

(*) :: Num a => a -> a -> a

In other words, the two arguments to this operator have to be the same type, so what you are asking for is not possible as it stands.

I can see a couple of options:

• As suggested in the comments, define your own (*), perhaps with a type class to go with it that standard types like Integer can also be a member of. This answer might provide some inspiration.

• Add scalars to your matrix type:

  data Matrix a = Scalar a | Matrix (Array (Int,Int) a)

  (perhaps the name could now be improved - also note that now it has to be data rather than newtype. I doubt the performance difference will matter in practice.)

Answer 2

Just an attemp to extend Ganesh's answer. We could redefine Scalar matrix as Unity matrix of unindentidied size.

import Data.List (transpose)

data Matrix a = Matrix [[a]] | UnityMatrix a deriving (Show)

instance Functor Matrix where
    fmap f (Matrix x) = Matrix $ fmap (fmap f) x
    fmap f (UnityMatrix x) = UnityMatrix $ f x 

fmap2::(Num a) => (a->a->b)->Matrix a->Matrix a->Matrix b    
fmap2 f (Matrix x) (Matrix y) = Matrix $ zipWith ( zipWith f ) x y
fmap2 f m@(Matrix x) u@(UnityMatrix  y) = fmap2 f m $ expandUnity u
fmap2 f u@(UnityMatrix  y) m@(Matrix x) = fmap2 f (expandUnity u) m
fmap2 f (UnityMatrix x) (UnityMatrix y) = UnityMatrix $ f x y

expandUnity (UnityMatrix a) = Matrix [replicate i 0 ++ a : repeat 0| i <-[0..]]

instance Num a => Num (Matrix a) where
    fromInteger = UnityMatrix . fromInteger

    signum = fmap signum
    abs = fmap abs

    (+) = fmap2 (+)
    (-) = fmap2 (-)    

    (Matrix x) * (Matrix y) = Matrix [[sum $ zipWith (*) a b | b <- transpose y ]| a <- x ]
    m@(Matrix x) * (UnityMatrix y) = fmap (*y) m
    (UnityMatrix y) * m@(Matrix x) = fmap (y*) m
    (UnityMatrix x) * (UnityMatrix y) = UnityMatrix $ x * y

main = print $  3 * Matrix [[1,2,3],[4,5,6]]  + 2

This code contains lot of repeating - but this in case if any of (+) , (-) or (*) operators in base type is not commutative. Also underlying type changed to list of lists instead of matrix. But this is only to ease of idea demo.