A Glitch In The Matrix.

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

Where r and c are 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

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

we get


 

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

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.

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:

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

and of course subtraction is very similar.

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

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

 

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:

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.

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:

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

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…

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!

 

Leave a Reply

Your email address will not be published. Required fields are marked *

ˆ Back To Top