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

The "semantics you're looking for" are just inconsistent with how the Haskell Num class is defined. There is no such thing as data promotion in Haskell (for very good reasons that are hard to explain but easy found out the more experience you gain with Haskell).

I'd not recommend defining a new *, nor your own operator name that has the signature a -> Matrix a -> Matrix a, it would be confusing. Rather, you should look where this operation is already defined: as user5402 has it, this is a scalar multiplication. This is obviously not just meaningful for matrices / linear operators, but already to vectors. And behold, here is the class for vector spaces!

But like suggested by Ganesh and Odomontois for Num, you'll also need to expand you data type to make proper use of this. The basic problem is that matrix multiplication is really only defined for matching covariant-contravariant dimensions, but your approach has no way to ensure this. Ideally, the type checker should infer the correct dimension, but instead of that you can have a special case, not just for the identity but general diagonal matrices.

data LinOp a = Diagonal a
             | GMatrix (Matrix a)

instance Functor LinOp where
  fmap f (Diagonal a) = Diagonal (f a)
  fmap f (GMatrix m) = GMatrix $ fmap f m

instance (Num a) => AdditiveGroup (Matrix a) where
  zeroV = Diagonal 0
  negateV = fmap negate
  (Diagonal x) ^+^ (Diagonal y) = Diagonal $ x+y
  ...

instance (Num a) => VectorSpace (Matrix a) where
  type Scalar (Matrix a) = a
  μ *^ m = fmap (μ*) m

instance (Num a) => Num (Matrix a) where
  fromInteger = Diagonal . fromInteger
  (+) = (^+^)
  ...

Answer 2

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 3

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.

Answer 4

Does this work?

scalarMult :: Num a => a -> Matrix a -> Matrix a
scalarMult c (Matrix x) = Matrix $ array (range x) [(ij,(x!ij) * c) | ij <- range x]

In general you can't implement scalar multiplication with fromInteger and ordinary matrix multiplication because you don't know what size matrix to create - it will depend on its context.

However, you could possibly use this approach if you kept track of the matrix size in the type system. Then the compiler could infer what size to promote scalars to, and it also could detect dimension mismatches for matrix operations at compile time.