Scalar Multiplication with a New Type of Matrix in Haskell

Question

I have been programming in languages ​​like C and Lisp for several decades and Haskell for several years. Now, to learn more about Haskell's classes and other great more advanced features, such as parallelism, I tried to create a new data type with matrix semantics, currently supported by a copy on top of the Array library type.

I am using the latest Haskell 2013.2.0.0 platform with ghc enabled.

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

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

Some simple functions for creating an identical matrix, for example, work very well:

 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] 

Also creates the 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 regular arithmetic operators to work, I add the following:

 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] 

(I apologize if the indentation is screwed here - the actual code is correctly set back and compiled.)

So far, it’s so good, although it’s a little annoying to define the fromInteger function, which actually does not make any sense in the semantics of the matrix, but creating a 1x1 matrix with a value is just as reasonable as everything else.

My problem is trying to get the correct semantics for multiplying a scalar (i.e. type of Num class) with a matrix. By mathematical convention, which simply means multiplication by type with a scalar.

However, no matter what syntax I try, I do not get the correct result. By default, this implementation only promotes, for example, Int to the matrix, using fromInteger, which (if the matrix is ​​already 1x1) results in the error "Unmatched matrix multiplication." I tried all the alternative syntaxes that I can think of in order to define alternative code for this type of multiplication without success, for example:

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

or

 (Matrix x) * (Num y) = 

or

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

But they all give me various syntax errors.

How do I determine the definition of a scalar by multiplying the matrix so that it does what you expect from it? I believe that enabling the custom security template feature may help, but I'm not quite sure how to use this in this context.

I understand that I could, by a special case, multiply the matrix to allow the multiplication of any matrix with a 1x1 matrix, which I think will work. But this would be: a) inelegant, (b) un-Haskell-y, (c) probably inefficient, since it would take each scalar to be wrapped in a matrix before multiplication and (d) would allow some code (e.g. multiplication of an arbitrary matrix with any 1x1 matrix) to start when it should lead to an error.

I also understand that there are probably excellent Matrix implementations that somehow circumvented this issue. But using them could defeat my goal to learn.

Answer 1

The semantics you are looking for simply contradicts the way the Haskell Num class is defined. There is no such thing as data promotion in Haskell (for very good reasons, which are difficult to explain but easy to find out the more experience you get with Haskell).

I would not recommend defining a new * , as well as the own name of the operator, which has the signature a -> Matrix a -> Matrix a , this will be confused. Rather, you should see where this operation is already defined: like user 5402, this is scalar multiplication. This obviously does not make sense for matrices / linear operators, but already for vectors. And now, here is the class for vector spaces !

But, as suggested by Ganesh and Odomontois for Num , you will also need to expand your data type in order to use this correctly. The main problem is that matrix multiplication is really only defined for matching covariant-contravariant measurements, but your approach has no way to provide this. Ideally, the checking type should output the correct size, but instead you can have a special case, not just for identical, but common diagonal matrices.

 data LinOp a = Diagonal a | GMatrix (Matrix a) instance Functor LinOp where fmap f (Diagonal a) = Diagonal (fa) fmap f (GMatrix m) = GMatrix $ fmap fm 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 standard operator (*) :

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

In other words, the two arguments to this operator must be of the same type, so what you are asking for is not possible because it is worth it.

I see a couple of options:

• As pointed out in the comments, define your own (*) , possibly with a type class, to go with it, which standard types like Integer can also join. This answer may give some inspiration.

• Add scalars to your matrix type:

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

  (Perhaps the name can now be improved - also note that now it should be data , not newtype . I doubt that the difference in performance will matter in practice.)

Answer 3

Just to expand the answer of Ganesh . We could redefine the Scalar matrix as a Unity unindentidied size matrix.

 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 $ fx fmap2::(Num a) => (a->a->b)->Matrix a->Matrix a->Matrix b fmap2 f (Matrix x) (Matrix y) = Matrix $ zipWith ( zipWith f ) xy fmap2 f m@ (Matrix x) u@ (UnityMatrix y) = fmap2 fm $ expandUnity u fmap2 f u@ (UnityMatrix y) m@ (Matrix x) = fmap2 f (expandUnity u) m fmap2 f (UnityMatrix x) (UnityMatrix y) = UnityMatrix $ fxy 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 (*) ab | 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 many repetitions - but this is the case if any of the (+), (-) or (*) operators in the base type is not commutative. Also, the base type is changed to a list of lists instead of a matrix. But this is only to facilitate the demonstration of ideas.

Answer 4

It works?

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

In the general case, you cannot implement scalar multiplication using fromInteger and ordinary matrix multiplication, because you do not know which size matrix to create - this will depend on its context.

However, you could use this approach if you were tracking the size of a matrix in a type system. Then the compiler could determine what size to advance the scalars, and it could also detect measurement mismatches for matrix operations at compile time.