In this article I want to talk about generalized matrix multiplication.

As we all know, for 2dimensional matrices, there is 2 common ways of multiplying (I will ignore per-element-multiplication for an instance):

The cross-product:

[1 2 [5 6 [(1*5 + 2*7) (1*6 + 2*8) [19 22

3 4] x 7 8] = (3*5 + 4*7) (3*6 + 4*8)] = 43 50]

And the dot-product:

[1 2 [5 6

3 4] .* 7 8] = 1*5 + 2*6 + 3*7 + 4*8 = 59

If we remember the rule about cross-multiplication, it states that we can multiply two matrices if their shared dimension is the same. Formally, we can multiply matrices:

k*n x n*l => k*n

We also know that multiplication of numbers is associative, hence there is no reason that the shared dimension can not be more or less than a dimension. In fact, if we re-examine the dot-product, we see it’s another instance of a cross-product:

1*(m*m) x (m*m)*1 => 1*1

Therefore, we can generalize this principle to enable the multiplication of matrices of M and N dimensions to become a matrix of (M-L)+(N-L) dimensions, if, for dimensions:

a_1*..*a_(m-l)*a_(m-l+1)*..*a_m

x

b_1*...*b_l*b_(l+1)*...*b_n

it holds that:

Forall k: 0 .. l-1 .

a_(m-l+k) = b_(k+1)

And then we will have result, a matrix of the dimensions:

a_1*...*a_(m-l)*...*b_(l+1)*...*b_n

If there’s several l’s for which the equations above holds true, it is obvious that there is more than one way to multiply the matrices. This is the reason why typically there’s no agreed upon way to multiply matrices. There’s too many different possibilities, the problem is underconstrained. However, if we add as constraint the result dimensions we need, we can solve this problem. Thanks to Haskell’s advanced type-system, we can do just that!

> {-# OPTIONS_GHC -fglasgow-exts -fallow-undecidable-instances -fallow-incoherent-instances -fallow-overlapping-instances #-}

> module Matrices where

First we define type level numbers:

> data One f a = One (f a) deriving (Eq, Show)

> data S n f a = S (n f a) (f a) deriving (Eq, Show)

> tPred :: S n f a -> (n f a)

> tPred = undefined

> type Two f = S One f

> type Three f = S (S One) f

> type Four f = S (S (S One)) f

Of course it is useful to be able to get the actual number encoded in the type information:

> class Number a where

> tNumber :: a f x -> Int

> tData :: a f x -> f x

> instance Number One where

> tNumber = const 1

> tData (One a) = a

> instance (Number n) => Number (S n) where

> tNumber x = 1 + tNumber (tPred x)

> tData (S n a) = a

Let us first define a generic equality typeclass. This will be reused later for more than just type-level numbers

Next, since we want to encode arbitrary dimensions, we need to be able to encode these dimensions at the type-level. Therefore, we need lists of these type-level numbers. Because we only care about lists of these type-level numbers and not more generic of other types, we can see the application of a number to another number as consing those numbers together. Of course, there must be some way to decide when here is no more numbers left. For this we introduce the marker Nil. Additionally, we need to be able to get the tail and the head of such a list of numbers. The head of such a type-level list is trivial, it’s just the type-number. As for the tail, the tail of a type-level list is the data contained within the number

> data Nil a = Nil a deriving (Eq, Show)

We want to also be able to determine the dimensions stored in such a list. For this we introduce another type-class:

> class Dimension d where

> tDimension :: d a -> [Int]

> instance Dimension Nil where

> tDimension = const []

> instance (Number na, Dimension nb) => Dimension (na nb) where

> tDimension x = tNumber x:tDimension (tData x)

Let’s test this with a simple example:

> testNumberDisplay = do

> let x = undefined :: (Two (Three (Four Nil))) ()

> print $ tNumber x

> print $ tDimension x

Now it is easy to store matrices as such constructs. Let’s add a utility function that constructs a matrix from a List:

> class ToMatrix m where

> toMatrix :: [a] -> (m a)

> toMatrixFold :: [a] -> (m a, [a])

> toMatrix xs = fst $ toMatrixFold xs

> instance ToMatrix Nil where

> toMatrixFold (x:xs) = (Nil x, xs)

> instance (ToMatrix m) => ToMatrix (One m) where

> toMatrixFold as = (One b, xs)

> where (b,xs) = toMatrixFold as

> instance (ToMatrix m, ToMatrix (n m)) => ToMatrix (S n m) where

> toMatrixFold as = (S n b, xs)

> where (b, as') = toMatrixFold as

> (n, xs) = toMatrixFold as'

It would be nice if it were possible to display these matrices of items. For this we will create a convenience class that turns a matrix into a nested list.

> class FromMatrix m b | m -> b where

> fromMatrix :: m -> b

> instance FromMatrix (Nil a) a where

> fromMatrix (Nil x) = x

> instance (FromMatrix (m a) b) => FromMatrix (One m a) [b] where

> fromMatrix (One a) = [fromMatrix a]

> instance (FromMatrix (m a) b,

> FromMatrix (n m a) [b]) => FromMatrix (S n m a) [b] where

> fromMatrix (S n a) = fromMatrix a:fromMatrix n

Now it’s time to test whether these different functionalities work as expected.

> testToMatrix = do

> let as = [1..]

> let x = toMatrix as :: (Two (Two Nil)) Int

> print $ fromMatrix $ x

> let x = toMatrix as :: (Three (Three (Three Nil))) Int

> print $ fromMatrix $ x

Finally, now that we have matrix structures that can be created and displayed easily, it is time to create a matrix multiplication algorithm. However, to keep it generic (for instance to enable the Floyd-Warshall Algorithm), we will try to keep the actually ‘addition’ and ‘multiplication’ used to remain generic. Therefore, we first introduce a semiring. The semi-ring for numbers is trivial. As for matrices, it is possible to introduce

First we define how to add matrices (treating them as higher dimensional matrices) and how to left-multiply them with constants:

> class (Num a) => Vector m a where

> (<+>) :: (m a) -> (m a) -> (m a)

> (*>) :: a -> (m a) -> (m a)

> instance (Num a) => Vector Nil a where

> (Nil a) <+> (Nil b) = Nil $ a + b

> a *> (Nil b) = Nil $ a * b

> instance (Num a, Vector m a) => Vector (One m) a where

> (One ma) <+> (One mb) = One $ ma <+> mb

> a *> (One mb) = One $ a *> mb

> instance (Num a, Vector m a, Vector (n m) a) => Vector (S n m) a where

> (S na ma) <+> (S nb mb) = S (na <+> nb) (ma <+> mb)

> a *> (S nb mb) = S (a *> nb) (a *> mb)

Next we define how to apply a dot-product of two matrices where the second one may have additional dimensions that are then treated as vectors.

> class (Vector v a) => DotMatrix m n v a | m n -> v where

> (<.>) :: m a -> n a -> v a

> instance (Vector v a) => DotMatrix Nil v v a where

> (Nil a) <.> vb = a *> vb

> instance (DotMatrix m n v a) => DotMatrix (One m) (One n) v a where

> (One ma) <.> (One mb) = ma <.> mb

> instance (DotMatrix m n v a,

> DotMatrix (o m) (o n) v a)

> => Dot Matrix (S o m) (S o n) v a where

> (S na ma) <.> (S nb mb) = (na <.> nb) <+> (ma <.> mb)

A simple test can now be performed:

> testMultMatrix = do

> let as = [1..]

> let x = toMatrix as :: (Two (Two Nil)) Int

> let y = toMatrix as :: (Two (Two (Three (Four Nil)))) Int

> let z = x <.> y :: (Three (Four Nil)) Int

> print $ fromMatrix $ x

> print $ fromMatrix $ z

Dot-product works fine, even when working on non-scalars in the second matrix. Sadly, I can not force the matrix cross-product to use the dot-product when required. As such, this is where I’m stuck. Comments are more than welcome.

The code I had tried can be seen below:

Finally, we can define the multiplication for matrices. If only this would work… sadly I do not know how to force dotmatrix of the right submatrices.

> class MultMatrix k l m a where

> (<*>) :: k a -> l a -> m a

> instance (DotMatrix (k Nil) (k m) m a) => MultMatrix (k Nil) (k m) m a where

> ma <*> mb = ma <.> mb

> instance (MultMatrix k l m a) => MultMatrix (One k) l (One m) a where

> (One ma) <*> mb = One $ ma <*> mb

> instance (MultMatrix k l m a,

> MultMatrix (o k) l (o m) a)

> => MultMatrix (S o k) l (S o m) a where

> (S na ma) <*> mb = S (na <*> mb) (ma <*> mb)

Stefan O’Rear told me that the reason it is not working is because GHC /= Prolog. More specifically, it always cuts when trying an option. Now the question, of course, is how to fix that.

Additionally, I can remove all options except for -fallow-undecidable-instances. If I remove that option, then the fundeps on FromMatrix and DotMatrix give errors. Removing them, however, makes the typechecker give more errors… So I wonder if the design is borked