On my shelf I have the book "Basic Topology" by Armstrong. After you've fought your way through 173 pages you eventually get to the section on simplicial topology and you can start playing with one of the basic tools of modern topology: homology groups . But you don't want to go through all that hassle. If you can read Haskell code then you can get there in 20 lines. But don't get too excited - this isn't meant to be an introduction to algebraic topology. Instead it's a demonstration of how incredibly expressive Haskell is. Before starting on this project I expected it to take a couple of weeks of early morning hacking and many hundreds of lines of code. It actually took a couple of hours.



So the goal is this: use the standard chain complex from simplicial homology to compute the Betti numbers.





First we need some code to compute the rank of a matrix. As this is general purpose code I'm not counting it in my 20 line budget. It took me a little while to find a recursive implementation of this suitable for Haskell. Basically, the idea is that a zero matrix has zero rank and the following operations leave the rank of a matrix unchanged:





Deleting a row of zeroes

Deleting a column of zeroes

Multiplying a row or column by a non-zero constant

Adding one row to another.





import List



rank m = let m' = removeZeroRows $ removeZeroColumns m in

if m'==[]

then 0

else let (zero,(pivot:rest)) = break ((0 /=) . head) m' in

1+rank (map (cancelZero pivot) (zero++rest))

where

removeZeroColumns x = let n = foldl1 min (map countZeros x) in

map (drop n) x

removeZeroRows = filter (or . map (0 /=))

countZeros = length . fst . break (0 /=)

cancelZero (a:as) (b:bs) = zipWith (\x y -> a*y-b*x) as bs





d x = zip (d' x) (cycle [1,-1]) where

d' [] = []

d' (x:xs) = xs : (map (x:) (d' xs))



matrix basis1 basis2 op = map (coefficients basis1 . op) basis2 where

coefficients basis x = map (flip coefficient x) basis

coefficient e = sum . lookup e

lookup k [] = fail ""

lookup k ((k',v):xs) = if k == k' then return v else lookup k xs



grade n = filter ((==(n+1)) . length)



dim c = foldl1 max (map length c)-1



dmatrix _ 0 = []

dmatrix c n = matrix (grade (n-1) c) (grade n c) d



betti c = let d = dim c

ranks = map rank $ map (dmatrix c) [0..d]

dims = map (\i -> length (grade i c)) [0..d] in

h dims ranks where

h (c:cs) (d:d':ds) = c-d-d' : h cs (d':ds)

h [c] [d] = [c-d]





One last thing is needed: if there is a non-zero value standing on its own in a row or column then deleting both the row and column that it is in reduces the matrix rank by one. Here's the complete code for the rank function:With that out of the way, here is the code to compute the homology of a simplicial complex:That's it! No hundred thousand line algebra library, or polygonal mesh geometry library. Just the standard Haskell List library.



At this point I suggest looking at something like these notes on simplicial complexes. My d function corresponds to ∂ defined at the bottom of page 1.





Let's take the example at the bottom of page 2 and the top of page 3:





example = [

[1,2,3],

[1,2],[1,3],[2,3],[2,4],[3,4],

[1],[2],[3],[4],[5]]



-1

(Note I don't have anything corresponding to F. This is always the set containing the empty set so the code is hard coded to work as if it were there.) If you play with the dmatrix routine you'll see that it gives the matrices on page 3 (transposed).



And now we can fire up ghci or hugs and get:





*Main> dmatrix example 1

[[-1,1,0,0,0],[-1,0,1,0,0],[0,-1,1,0,0],[0,-1,0,1,0],[0,0,-1,1,0]]

*Main> betti example

[2,1,0]



The example has 2 connected components, 1 1-hole, and no 2-holes. (Quick one line summary of the whole of algebraic topology: Make a hole in a piece of paper: that's a 1-hole. Hollow out the interior of a solid sphere: that's a 2-hole, etc.)



Here are some more topological spaces to play with:





point = [[0]]

line = [[0],[1],[0,1]]

ball n = sublists [0..n] \\ [[]]

sphere n = ball (n+1) \\ [[0..(n+1)]]

torus = [[1],[2],[3],[4],[5],[6],[7],[8],[9],

[1,2],[2,3],[1,3],

[5,9],[8,9],[5,8],

[4,6],[6,7],[4,7],



[1,5],[4,5],[1,4],

[2,9],[6,9],[2,6],

[3,8],[7,8],[3,7],



[1,9],[3,9],[1,8],

[4,9],[6,8],[5,7],

[1,6],[2,7],[1,7],



[1,5,9],[1,2,9],[2,3,9],[3,8,9],[1,3,8],[1,5,8],

[4,5,9],[4,6,9],[6,8,9],[6,7,8],[5,7,8],[4,5,7],

[1,4,6],[1,2,6],[2,6,7],[2,3,7],[1,3,7],[1,4,7]]





*Main> betti (sphere 3)

[1,0,0,1]

*Main> betti (sphere 0)

[2]

*Main> betti torus

[1,2,1]





sublists [] = [[]]

sublists (a:as) = let b = sublists as in b ++ map (a:) b



properSublists x = sublists x \\ [x]



chain f x = [x:xs| x'

bary x = x >>= chain (filter (not . null) . properSublists)





*Main> betti (bary torus)

[1,2,1]

*Main>



dmatrix ((bary torus)) 2

The torus is copied from page 4 of the notes. Here are some more examples:But that's not all. Pay your $19.95 right now and I'll throw in barycentric division for free:As expected, barycentric division of a simplicial complex leaves the homology unchanged:OK, I admit, it's a little slow on the last one. But it does have to perform unoptimised Gaussian elimination on fairly large matrices. (Type



And here's my laundry list of things to do with this project:





More efficiency. The vectors over the simplices are stored incredibly inefficiently.

Compute homology groups over Z, not Q. I just use dimension counting to compute the dimensions of the homology groups over a Q. I really ought to actually compute the kernels of these matrices and then find representatives of the homology groups. I also need to write code to manipulate Abelian groups rather than vector spaces. To my shame, I don't know what the algorithms to use for the latter are.

Compute cohomology groups.

Compute homology groups associated to a poset. (The barycentric subdivision code kinda does this already.) Actually, this is where I started from originally. alpheccar asked if there was a connection between certain types of Hopf algebra and the Euler characteristic. I think the answer is yes: it comes from computing the Euler characterstic of posets used to define the Hopf algebra. Actually, the whole renormalisation with Hopf algebras thing that alpheccar points out is reminiscent of a really beautiful paper I read years ago on "fat graphs" that computes the Euler characteristics of the very moduli spaces used in string theory.





euler c = sum (zipWith (*) (betti c) (cycle [1,-1]))



And the last thing on my list I can do right now. Computing the Euler characteristic, via the Betti numbers, is a one liner:Of course I showed there was an easier way to compute this earlier



How can Haskell not be the programming language that all mathematicians should learn?





Update: I just found some slides on the "fat graph" thing I mentioned above. I rate it as one of the most beautiful pieces of mathematics ever. Some time I'll devote a whole post to it.

Labels: haskell, mathematics