This is the first in a series of posts that will eventually get round to implementing a Kalman Filter (or filters) in Haskell. It’s not my intention to explain Kalman filtering in any detail – just sufficiently to complelent the Haskell code. There’s a plethora of tutorials and videos about Kalman Filtering and one of the better series is here.

As will be seen, the computations involved in Kalman Filtering are effectively expressed using matrices. The Haskell ecosystem has a number of Matrix libraries but I’ll learn more by writing my own library even if its scope is limited to the immediate requirement of Kalman Filter computations.

When thinking about implementing basic functions for matrices I wondered how – or if – the number of rows and columns could be encoded in the type rather than as a parameter in the matrix itself. So what follows is my limited exploration of type jiggery in Haskell.

We can define a Matrix type

1 |
newtype Matrix (r :: Nat) (c :: Nat) a = Matrix [[a]] deriving Show |

Where *r* and *c* or so called Phantom Type. i.e. Types that appear in the type constructor but not in the data constructor. It’s a form of ‘tag’ on the data type that exists at compile time but is lost at runtime. Having the tags a compile time prevents errors being introduced in the first place rather than having to handle them at runtime.

### Matrix Multiplication

So what sort of errors can this phantom type help catch at compile time? Matrix multiplication is an example and was the motivation behind these experiments.

Multiplication of two matrices is only defined if the number of columns in the first matrix is equal to the number of rows in the second.

**M1 = [m, n]**

**M2 =[p, q]**

then for multiplication to be valid ** n == p** and the result is

**M3 = [m, q]**And the corresponding function signature is

1 |
mul :: Num a => Matrix n m a -> Matrix m p a -> Matrix n p a |

and at the moment it doesn’t matter what the implementation is – a compile time error will occur if the *mul* function is called without the correct types.

e.g. If we try to compile this snippet

1 2 3 4 5 6 7 8 |
mul :: Num a => Matrix n m a -> Matrix m p a -> Matrix n p a mul = undefined mu :: Matrix 3 3 Int mu = Matrix [[1,2,2],[1,1,2], [4,5,6]] mv :: Matrix 4 2 Int mv = Matrix [[1,2],[3,4], [1,1], [4,4]] |

we get

1 2 3 4 5 6 7 8 9 10 11 |
Matrix.hs:89:15: error: • Couldn't match type ‘4’ with ‘3’ Expected type: Matrix 3 2 Int Actual type: Matrix 4 2 Int • In the second argument of ‘mul’, namely ‘mv’ In the expression: mul mu mv In an equation for ‘prod’: prod = mul mu mv | 89 | prod = mul mu mv | ^^ Failed, no modules loaded. |

The implementation of matrix multiplication is quite a nice algorithm to implement in Haskell and so I’ll provide two!

1 2 3 4 5 6 7 8 |
mul :: Num a => Matrix n m a -> Matrix m p a -> Matrix n p a mul (Matrix xs) (Matrix ys ) = Matrix $ with xs $ \rowi -> with ys' $ \coli -> sum $ zipWith (*) rowi coli where ys'= L.transpose ys with :: [a] -> (a -> b) -> [b] with = flip fmap |

here *with* – which is just *fmap* with its arguments swapped round – is used to ‘offer up’ a row to each column in turn, multiply corresponding entries and then add them all up.

An alternative – that doesn’t explicitly use map is below.

1 2 3 4 5 6 |
mul' :: Num a => Matrix n m a -> Matrix m p a -> Matrix n p a mul' (Matrix xs) (Matrix ys ) = Matrix $ L.transpose . map (applyRow xs) . L.transpose $ ys applyRow :: Num a => [[a]] -> [a] -> [a] applyRow [ys] xs = [sum $ zipWith (*) xs ys] applyRow (ys:yss) xs = sum (zipWith (*) xs ys) : applyRow yss xs |

### Matrix Addition and Subtraction

Adding two matrices requires that both matrices have the same number of rows and columns. So this can be constrained by the function signature:

1 |
add :: Num a => Matrix n m a -> Matrix n m a -> Matrix n m a |

and an implementation of *add* just needs to add corresponding elements in a ‘*zipWith zipWith…*‘ style like this.

1 2 |
add :: Num a => Matrix n m a -> Matrix n m a -> Matrix n m a add (Matrix xs) (Matrix ys ) = Matrix $ zipWith (zipWith (+)) xs ys |

and of course subtraction is very similar.

1 2 |
sub :: Num a => Matrix n m a -> Matrix n m a -> Matrix n m a sub (Matrix xs) (Matrix ys ) = Matrix $ zipWith (zipWith (-)) xs ys |

Clearly *add* and *sub* are very similar, they both embed a binary operation and this can be abstracted into a separate function like this

1 2 3 |
-- A binary operation on square matrices. sqOp :: Num a => (a -> a -> a) -> Matrix n m a -> Matrix n m a -> Matrix n m a sqOp f (Matrix xs) (Matrix ys ) = Matrix $ zipWith (zipWith f) xs ys |

and now *add* and *sub* can now be written in a much simpler way:

1 2 3 4 5 6 |
-- add' :: Num a => Matrix n m a -> Matrix n m a -> Matrix n m a add' = sqOp (+) sub' :: Num a => Matrix n m a -> Matrix n m a -> Matrix n m a sub' = sqOp (-) |

### Matrix Division

The division considered here is division of corresponding elements between two square matrices of the same size.

The functions for addition and subtraction have had a type constraint of *(Num a*) – however the type *Num* does not support a division operation as can be seen from GHCi:

1 2 3 4 5 6 7 8 9 10 11 |
-- λ-> :i Num class Num a where (+) :: a -> a -> a (-) :: a -> a -> a (*) :: a -> a -> a negate :: a -> a abs :: a -> a signum :: a -> a fromInteger :: Integer -> a {-# MINIMAL (+), (*), abs, signum, fromInteger, (negate | (-)) #-} |

So the problem here how does the division operation get resolved from just a *Num* type? In other words which division operation is actually used? Is it integer type division or the division of double types? A neat solution to this is to create a custom type class that resolves the actual division operation and use this new type class as a constraint on the division operation. Like this.

1 2 3 4 5 6 7 8 9 10 11 |
-- class DivSupported a where divOp :: a -> a -> a instance DivSupported Int where divOp = div instance DivSupported Integer where divOp = div instance DivSupported Double where divOp = (/) instance DivSupported Float where divOp = (/) divEls :: (DivSupported a, Num a) => Matrix n m a -> Matrix n m a -> Matrix n m a divEls = sqOp divOp |

and division now just becomes another use of the *sqOp* function. 🙂

### A Glitch In the Matrix

As I was writing the Haskell code for matrices using types to encode constraints I kept feeling something was missing. And it was that defining a type signature such as
ma :: Matrix 5 3 Double is all well and good there’s nothing to ensure that the matrix does actually have 5 rows and 3 columns. And the issue came again when I looked at making a *Monoid* instance for the Matrix type. As a reminder, a *Monoid* is a *Semigroup* with an identity element and a *SemiGroup* is a binary operation defined over a set (i.e a type). I often find it instructive and enjoyable to define Monoids and Semigroups etc. on types and Matrix is no exception.

Using GHCi to get a reminder of a *SemiGroup*:

1 2 3 4 5 6 7 |
-- λ-> :i Semigroup class Semigroup a where (<>) :: a -> a -> a GHC.Base.sconcat :: GHC.Base.NonEmpty a -> a GHC.Base.stimes :: Integral b => b -> a -> a {-# MINIMAL (<>) #-} |

we see that the operation <> needs defining. I decided to use matrix multiplication (addition would be just as ‘good’). Clearly only multiplication on square matrices could form a *SemiGroup* so the *SemiGroup* instance for square matrices is

1 2 3 |
-- instance (Num a ) => Semigroup ( Matrix (n :: Nat) (n :: Nat) a) where (<>) = mul |

restricting it to *n x n* matrices. A *Monoid* is a *SemiGroup* with an identity element and for a *[n, n]* matrix the identity element is the matrix with 1s down the diagonal and 0s elsewhere. The problem with this is that for square matrix *[n, n]* the *n* is encoded in the type and I don’t see how to ‘extract’ it – or if indeed it is possible to do so. i.e. What I’d like to write for a matrix *Monoid* would be…

1 2 3 4 5 6 7 8 9 10 |
-- instance (Num a ) => Monoid ( Matrix (n :: Nat) (n :: Nat) a) where mempty = Matrix $ identity n -- Creates array of arrays with 1 in the main diagonal and 0 elsewhere. identity :: (Num a ) => Int -> [[a]] identity n | n <= 0 = [[]] | otherwise = take n (fmap (take n) rc) where rc = (1 : repeat 0) : fmap (0 :) rc |

but, not unreasonably, the value for *n* in the *Monoid* is not available. At this point I reached the limit of my ‘type voodo’ skills in Haskell and will revisit the matrix and also have the row/column sizes as part of the data constructor as well as the type constructor. That way there’s the type safety at compile time and the flexibility of having the row/column sizes available at runtime.

The code is here on GitHub.

And thanks for reading!