Bird Tracks Through Math Land: Basic Matrix Ops

Math scares me. I hate to say it but it's true. I write software for super smart researchers/scientists who are solving all kinds of complex mathematical problems in computational biology and this discomfort with math has become a limiting factor in how much I can contribute to my group. So, it's time for me to fix the problem.

I don't learn well from just reading books on math. I have tried it. I always end up feeling like I've learned a few rules but I don't really "get" it. I need to have a sense for what the small pieces are and how they all fit together to form a larger system of high level abstractions/components. That's how I build software and it's the only way I am comfortable learning technical concepts. That's what got me thinking about "Bird tracks through math land". The best way for me to learn these concepts is to implement them! If you build something you have to understand it right?

I am going to do this is through a series of Literate Haskell posts. The idea is to take it slow, and explain things as fully as possible. In the end I would like to have a series of math tutorials/notes that also happen to be working code (there is nothing like a compiler to keep you honest!). My hope is that these notes will be almost completely self contained and should not require much more math background than what you get in high school. You will have to be able to read the Haskell source code though (I very highly recommend Learn You a Haskell for anyone new to Haskell).

A bit of background for anyone who is not familiar with Literate Haskell:

• Haskell is a purely functional programming language. That means that the output of a function is completely determined by it's inputs (no cheating with hidden state). This aligns closely with the mathematical definition of a function, but is very different from main stream languages like C, Java ...
• Haskell's syntax is beautiful and extremely concise. It is nearly as expressive as writing pseudocode which makes it perfect for this kind of project.
• You're reading Literate Haskell source right now... Seriously! I wouldn't lie to you. Literate Haskell treats every line that does not start with a '>' (AKA bird track) as a comment. Quoting Dr. Knuth, who conceived Literate Programming "The main idea is to regard a program as a communication to human beings rather than as a set of instructions to a computer." Which is exactly what I am trying to achieve.

The road map for "Bird tracks through math land" is still up in the air, but for now I plan on starting with Linear Algebra and working my way toward its applications in probability and statistics. I am also very interested in Bayesian Networks so I think we'll head down that path later.

OK, enough blah blah. Lets get to it. We'll start by defining the basic matrix operations.

> module Math.BTTML.BasicMatrixOps (
>           dotProd,
>           matAdd,
>           matSub,
>           matMult,
>           eqZipWith,
>           matZipWith) where
>
> import Data.List

Let's define our own version of Haskell's zipWith function. eqZipWith applies the function f to each corresponding pair of elements in the given lists. It is just like Haskell's version, except that eqZipWith requires that the list lengths match.

> eqZipWith :: (a -> b -> c) -> [a] -> [b] -> [c]
> eqZipWith f (x:xs) (y:ys) = (f x y):(eqZipWith f xs ys)
> eqZipWith f [] [] = []
> eqZipWith _ _ _ = error "List sizes should match"

Now let's create a matrix zip function. It applies a function f to each pair of corresponding elements from the two given matrices.

> matZipWith :: (a -> b -> c) -> [[a]] -> [[b]] -> [[c]]
> matZipWith f matrix1 matrix2 = eqZipWith (eqZipWith f) matrix1 matrix2

Subtract matrix1 from matrix2. This is an element-by-element operation.

> matSub :: (Num a) => [[a]] -> [[a]] -> [[a]]
> matSub = matZipWith (-)

Matrix addition works pretty much the same way

> matAdd :: (Num a) => [[a]] -> [[a]] -> [[a]]
> matAdd = matZipWith (+)

The dot product (AKA scalar product) of two vectors is the sum of the product of each pair of elements. The vectors have to be the same length.

> dotProd :: (Num a) => [a] -> [a] -> a
> dotProd vector1 vector2 = sum $ eqZipWith (*) vector1 vector2

The wikipedia page on matrix multiplication is great: http://en.wikipedia.org/wiki/Matrix_multiplication

The most important points are:

• The row count of the resulting matrix is determined by the row count of matrix1
• The column count of the resulting matrix is determined by the column count of matrix2
• element[i, j] of the resulting matrix is calculated by taking the dot product of the ith row vector from matrix1 and the jth column from matrix2. Also note that because these two vectors must be the same length, the column count of matrix1 must equal the row count of matrix2.

In order to simplify the matrix multiplication we will transpose matrix2 (moves element[i,j] to element[j,i]) which makes it trivial to pull out column vectors

> matMult :: (Num a) => [[a]] -> [[a]] -> [[a]]
> matMult matrix1 matrix2 =
>   [[row1 `dotProd` col2 | col2 <- transpose matrix2] | row1 <- matrix1]

Now here are a couple of tests:

> testMat1 =
>   [[2, 0, -1, 1],
>    [1, 2, 0, 1]]

> testMat2 =
>   [[1, 5, -7],
>    [1, 1, 0],
>    [0, -1, 1],
>    [2, 0, 0]]

> testAdd = testMat1 `matAdd` testMat1
> testSub = testAdd `matSub` testMat1
> testMult = testMat1 `matMult` testMat2

Next time we'll solve linear equations with gaussian elimination

Updates: Applied CSS colors and made some code changes suggested by Conal