I recently wrote this post. It contained a lot of rather obtuse mathematics, just to introduce what was basically the example problem for the article. That’s because it was a real honest-to-goodness problem I was solving and writing papers on for algebra journals… admittedly, not the best choice for a blog aimed partly at non-mathematicians. Here, I do something similar, but with a much simpler and more general purpose problem.
Overview
We’ll be playing with converting matrices to an upper triangular form, essentially using Gaussian elimination. As a reminder, here’s what that means (I’m simplifying a little):
Goal: I have a matrix, and I’d like it to be upper triangular. (Upper triangular means that everything below the main diagonal is zero.)
Rules: I’ll allow myself to do these things to the matrix: swap any two rows or columns, or add any multiple of a row or column to another one.
Strategy: I’ll look for a non-zero element in the first column. If there is one, I’ll swap rows to move it to the first row of the matrix. Then I’ll add multiples of the first row to all of the other rows, until the entire rest of the first column is equal to zero. Then (here’s where we get a little tricky) I can just do the same thing on the rest of the matrix, ignoring the first row and column. In other words, I’ve reduced the problem to a smaller version of itself. (Yep, it’s a recursive algorithm.)
The Trick Up Our Sleeve: Instead of writing our function to operate directly on some representation of matrices, we’ll make it work on a type class. This will let us play all sorts of cool tricks.
Preliminary Stuff
I’d like to be able to declare instance for my matrices without jumping through newtype hoops, so I’ll start with a language extension.
{-# LANGUAGE TypeSynonymInstances #-}
Imports:
import Data.Maybe
Next, I need a few easy utility functions on lists. There’s nothing terribly interesting here; just a swap function, and a function to apply a function to the nth element of a list.
swap :: Int -> Int -> [a] -> [a]
swap i j xs | i == j = xs
| i > j = swap j i xs
| otherwise = swap' i j xs
where swap' 0 j (x:xs) = let (b,xs') = swap'' x (j-1) xs in b : xs'
swap' i j (x:xs) = x : swap' (i-1) (j-1) xs
swap'' a 0 (x:xs) = (x, a:xs)
swap'' a j (x:xs) = let (b,xs') = swap'' a (j-1) xs in (b, x:xs')
modifynth :: Int -> (a -> a) -> [a] -> [a]
modifynth _ _ [] = []
modifynth 0 f (x:xs) = f x : xs
modifynth n f (x:xs) = x : modifynth (n-1) f xs
Finally, I need a type to represent matrices. Since this is just toy code where I don’t need really high performance, a list of lists will do just fine.
type Matrix = [[Double]]
All done. On to the interesting stuff.
Building a Type Class
I already mentioned that I don’t want to operate directly on the representation of a matrix as a list of lists. Instead, I’ll declare a type class capturing all of the operations that I’d like to be able to perform.
class Eliminable a where
(@@) :: a -> (Int,Int) -> Double
size :: a -> Int
swapRows :: Int -> Int -> a -> a
swapCols :: Int -> Int -> a -> a
addRow :: Double -> Int -> Int -> a -> a
addCol :: Double -> Int -> Int -> a -> a
I’ve reserved the operator @@ to examine an entry of a matrix, size to give me its size, and then included functions to swap rows, swap columns, add a multiple of one row to another, and add a multiple of one column to another. Now I just need an implementation for the concrete matrix type I declared earlier.
instance Eliminable Matrix where
m @@ (i,j) = m !! i !! j
size m = length m
swapRows p q m = swap p q m
swapCols p q m = map (swap p q) m
addRow k p q m = modifynth q (zipWith comb (m!!p)) m
where comb a b = k*a + b
addCol k p q m = map (\row -> modifynth q (comb (row!!p)) row) m
where comb a b = k*a + b
Done.
Programming With Our Type Class
Recall that the substantial portion of the elimination algorithm earlier was to zero out most of the first column of the matrix, leaving only the top element possibly non-zero. We’re now in a position to implement this piece of the algorithm. It’s not all that tricky.
zeroCol :: Eliminable a => a -> a
zeroCol m = let clearRow j m' = addRow (-(m' @@ (j,0) / m' @@ (0,0))) 0 j m'
clearCol m' = foldr clearRow m' [1 .. size m - 1]
in case listToMaybe [ i | i <- [0 .. size m - 1],
m @@ (i,0) /= 0 ]
of Nothing -> m
Just 0 -> clearCol m
Just i -> clearCol (swapRows 0 i m)
This function looks only at the first column of the matrix, and clears it out my moving a non-zero element to the top, and then adding the right multiple of that column to all those below it. The important thing to notice is that this was implemented for any arbitrary instance of the type class I called “Eliminable”. This will be incredibly useful in the next few steps.
Using the Type Class
The next part of the algorithm is to ignore the first row and column, and perform the same operation on the submatrix obtained by deleting them. It’s actually a bit unclear how to implement this. We have a few options:
- Modify zeroCol above, to have it take a parameter representing the current column, and do everything relative to the current column. This is pretty messy. It actually might not be too messy in this case, but if the algorithm I were implementing were a little less trivial to begin with, It might definitely be quite messy.
- Actually perform the elimination on a separate matrix, and then somehow graft the first row and column from this matrix onto that one. Again, this could get pretty messy in general.
- Change the representation of the submatrices.
I’ll choose the third. Luckily, this isn’t too tough, since we have a type class. I’ll just define a newtype, and a new instance, to encapsulate the idea of a matrix with the first row and column deleted.
newtype SubMatrix a = SubMatrix { unwrap :: a }
instance Eliminable a => Eliminable (SubMatrix a) where
(SubMatrix m) @@ (i,j) = m @@ (i+1,j+1)
size (SubMatrix m) = size m - 1
swapRows p q (SubMatrix m) = SubMatrix (swapRows (p+1) (q+1) m)
swapCols p q (SubMatrix m) = SubMatrix (swapCols (p+1) (q+1) m)
addRow k p q (SubMatrix m) = SubMatrix (addRow k (p+1) (q+1) m)
addCol k p q (SubMatrix m) = SubMatrix (addCol k (p+1) (q+1) m)
Using this new instance, I can easily complete the elimination algorithm.
eliminate :: Eliminable a => a -> a
eliminate m | size m <= 1 = m
| otherwise = unwrap . eliminate . SubMatrix . zeroCol $ m
Yep, that’s all there is to it, and we have a working elimination algorithm.
Type Class Games
Suppose, now, that I want a lower triangular matrix. It might initially seem that I’m out of luck; I need to write all this code again. That turns out not to be the case, though. If I just teach the existing code how to operate on the transpose of a matrix instead of the matrix I’ve given it, then all is well! Here goes.
newtype Transposed a = Transposed { untranspose :: a }
instance Eliminable a => Eliminable (Transposed a) where
(Transposed m) @@ (i,j) = m @@ (j,i)
size (Transposed m) = size m
swapRows p q (Transposed m) = Transposed (swapCols p q m)
swapCols p q (Transposed m) = Transposed (swapRows p q m)
addRow k p q (Transposed m) = Transposed (addCol k p q m)
addCol k p q (Transposed m) = Transposed (addRow k p q m)
To implement the lower triangular conversion, now, is simple.
lowerTriang :: Eliminable a => a -> a
lowerTriang = untranspose . eliminate . Transposed
Any number of changes to the operation we’re trying to perform can often be expressed by simply substituting a different representation for the type on which we’re performing the operation. (Thinking about this fact can actually get pretty deep.)
Side Calculations
There’s a fairly common problem that many people run into when moving from an imperative language to a functional one. This can apply to learning functional programming, converting existing imperative code, or even just translating the concepts in one’s mind when talking to someone who thinks imperatively. The problem goes something like this: you have some code that performs some computation, and now you want to change the code to add some new concept to the existing computation. Often, the new idea you’re trying to add could be performed trivially in an imperative language, by adding print statements to some function seven layers in, or by keeping track of some value in a global variable, or in some field of some object. In the functional setting, these aren’t available to you.
The minimalist answer is simply to add all the plumbing code; new parameters and return values, etc. to every function in the entire call tree. To say the least, this is unappealing! To a new Haskell programmer, the obvious answer often looks like monads. However, again, the entire call tree has to be rewritten in a monadic style, and besides, this is a tad like using a rocket launcher to rid the house of mice.
The solution I propose here is that many times, it’s sufficient to use a type class. Here’s an example.
Problem: Calculate the determinant of a matrix efficiently.
Determinants can be calculated in a lot of different ways, but one of the most common uses elimination. The interesting fact here is that once you’ve got a triangular matrix (lower or upper; doesn’t matter), then its determinant is just the product of its diagonal elements. Furthermore, we know precisely what happens to the determinant when you swap rows or columns (it flips sign, but the magnitude stays the same), or when you add a multiple of one row or column to the other (it stays the same). So a (very fast) way to calculate a determinant is to perform elimination, but also keep track, at each step, of what you’ve done to the determinant so far.
So now we need, not merely a matrix, but a pair consisting of a matrix and some side information – namely, which change we’ve made so far to the determinant.
data WithDeterminant a = WithDeterminant Double a
instance Eliminable a => Eliminable (WithDeterminant a) where
(WithDeterminant _ m) @@ (i,j) = m @@ (i,j)
size (WithDeterminant _ m) = size m
swapRows p q (WithDeterminant d m) = WithDeterminant (-d) (swapRows p q m)
swapCols p q (WithDeterminant d m) = WithDeterminant (-d) (swapCols p q m)
addRow k p q (WithDeterminant d m) = WithDeterminant d (addRow k p q m)
addCol k p q (WithDeterminant d m) = WithDeterminant d (addCol k p q m)
As before, once we’ve defined the appropriate instance, the implementation is actually quite easy.
diags :: Matrix -> [Double]
diags [] = []
diags (r:rs) = head r : diags (map tail rs)
determinant :: Matrix -> Double
determinant m = let WithDeterminant d m' = eliminate (WithDeterminant 1 m)
in d * product (diags m')
The resulting determinant function actually performs quite well. For example, calculation of the determinant of a 100 by 100 matrix is done in 1.61 seconds, fully interpreted in GHCi. I didn’t bother compiling with optimization to see how well that does, nor replacing the inefficient list-of-lists representation of matrices with one based on contiguous memory or arrays. (Edit: Compiled and optimized with GHC, but still using the list-of-list representation, the time is around a third of a second.)
(It’s worth pointing out that automatic differentiation is another very impressive example of this same technique, except that it uses the standard numeric type classes instead of a custom type class.)
Conclusion
The point of this article is that plain old type classes in Haskell can be used to make your purely functional code very flexible and versatile. By defining a type class to capture the concept of a set of related operations, I was able to achieve:
- Choice of data structures. Had I wanted to use a contiguous array instead of a list of lists, I could have easily done so.
- Easier programming. For example, operating on a submatrix of the original matrix became much easier.
- Flexible code. I was able to get lower triangular matrices, too, without rewriting the code.
- Better composability. I easily reused by upper triangular matrix calculation to find determinants, even though additional calculations had to be traced through the original.
- Separation of concerns. When I started, I never even dreamed that I might need to trace determinant calculations through the process. That got added later on, in it’s own separate bit of code. If someone else wanted different plumbing… say, for logging, or precondition checking, or estimating the possible rounding error, all they’d need to do is define a new instance of the type class.
None of this is new, of course. But type classes are definitely an underused language feature by many Haskell programmers.
, a function 
is called a left inverse of f in case 
. In other words, we want g to undo whatever f did to its parameter. A function has a left inverse if and only if it is injective, a.k.a. one-to-one.
is constructively valid. But it isn’t even classically valid, so that’s hopeless. Of course, if the domain of f is enumerable, I can in principle simply try all inputs to f in parallel until I find one that works… this is not terribly satisfying, however; it’s completely impractical for most purposes.
