Hopf Algebra = Group + Monad

I have two goals in this post. Firstly I want to revisit the vector space monad I've used before. This time I want to show how it provides a great domain specific language for working with multilinear functions on vector spaces. Secondly I want to apply this monad to talk about Hopf algebras. In particular I want to address the issue of why Hopf algebras and groups are actually the same thing. In fact, I want to push this idea to the point where I can use exactly the same code to perform operations in both groups and Hopf algebras. By parameterising the definition of groups in the right way, Hopf algebras appear as if by magic simply by tweaking one parameter. What this means is that if you understand groups and vector spaces, and are comfortable with monads in Haskell, then despite their notoriety as something whose definition is a little opaque, you'll have some idea of what a Hopf algebra is well before the end this post.

Of course there's the question of why anyone would want to know about Hopf algebras. They appear in combinatorics, topology and knot theory, and theoretical physics among other subjects. In fact, there's a review paper on their ubiquity. As topics in mathematics go, this is a hot one.

First some administrivia, because this is, as always, a literate Haskell post.

> {-# LANGUAGE MultiParamTypeClasses,FlexibleInstances,FunctionalDependencies,GeneralizedNewtypeDeriving #-}

> module Main where

> import qualified Data.Map as M
> import Control.Monad
> import Test.QuickCheck

> infixl 5 .+
> infixl 6 .*

> swap (x,y) = (y,x)

Vector Spaces

We'll start with a type class for vector spaces v over a base field k:

> class Num k => VectorSpace k v | v -> k where
>     zero :: v
>     (.+) :: v -> v -> v
>     (.*) :: k -> v -> v
>     (.-) :: v -> v -> v
>     v1 .- v2 = v1 .+ ((-1).*v2)

.+ is vector space addition and .* is scalar multiplication by elements of the base field.

And now I need a concrete instance. I'll be working with vector spaces equipped with a basis. I'll use the obvious representation: a list of pairs of coefficients and basis elements. So 2x-3y is represented as V [(2,x),(3,y)]. Sometimes it will be clearer if I write 2e_x-3e_y. x and y aren't really the basis elements but labels. The e_x and e_y are the actual basis elements. But it's often more readable to write x and y.

> data V k a = V { unV :: [(k,a)] } deriving (Show)

The problem with this representation is that basis elements might appear in any order, the same basis element might appear multiple times, and some basis elements might appear with redundant zero coefficients. We need a suitable definition of equality. (I'm sidestepping the usual issues with Eq and Monad.)

> reduce x = filter ((/=0) . fst) $ fmap swap $ M.toList $ M.fromListWith (+) $ fmap swap $ x

> instance (Ord a,Num k) => Eq (V k a) where
>   V x==V y = reduce x==reduce y

> instance (Ord a,Num k,Ord k) => Ord (V k a) where
>   compare (V x) (V y) = compare (reduce x) (reduce y)

V is a functor. fmap f just applies f to the basis elements.

> instance Num k => Functor (V k) where
>     fmap f (V as) = V $ map (\(k,a) -> (k,f a)) as

Given any type a, V k a is the type of vectors over k generated by elements of a. So V k (V k a) is the type of vectors whose basis elements are themselves labeled by vectors in V k a. We can think of an element of V k (V k a) looking something like 2(2x+3y)-5(3x-y). It's tempting to say this is equal to what we get when we multiply out, but where I wrote '2x+3y' I really meant e_2x+3y. Nonetheless, we can define a function that does the multiplying out. I'll call this function join. So we get a function join :: V k (V k a) -> V k a. In fact, we get a monad.

> instance Num k => Monad (V k) where
>     return a = V [(1,a)]
>     x >>= f = join (fmap f x)
>         where join x = V $ concat $ fmap (uncurry scale) $ unV $ fmap unV x
>               scale k1 as = map (\(k2,a) -> (k1*k2,a)) as

What's really going on here is that the "vector space over basis" functor is a monad in the category Set, but in the category of Haskell types and functions we're stuck with lists.

V k is not just a monad, it's also a MonadPlus.

> instance Num r => MonadPlus (V r) where
>     mzero = V []
>     mplus (V x) (V y) = V (x++y)

So what's the use of making our vector space into a monad? If m is a monad, then any function f:a→mb can be lifted to a function f':ma→mb. f maps the basis elements in a to elements of mb. It's a standard fact about vector spaces that a linear function is completely defined by its value on a basis. f' is the corresponding linear function. So we only need to implement linear functions for basis elements and the function can be automatically lifted to a linear function on the whole vector space. An arrow of the form a→mb is a Kleisli arrow and what I've just argued is that the Kleisli category of the V k monad is the category of vector spaces with basis.

This gives us lots of advantages. For example, we can't accidentally implement a function that isn't linear this way. If we implement the action on a basis, the lift can't help but be linear. Also, working with basis elements rather than the vector space gives us greater flexibility. For example, it'll be easy to form tensor products of vector spaces by directly working with their bases.

Here's the instantiation of the vector space type:

> instance (Num k,Ord a) => VectorSpace k (V k a) where
>     zero = V []
>     V x .+ V y = V (x ++ y)
>     (.*) k = (>>= (\a -> V [(k,a)]))

> e = return :: Num k => a -> V k a
> coefficient b (V bs) = maybe 0 id (lookup b (map swap (reduce bs)))

I could have implemented .* using scale but this definition is more in the spirit of what's to come.

Groups

Now we switch to groups. Here's a first attempt at writing a group type class

> class Group1 a where
>    unit1 :: a
>    mult1 :: a -> a -> a
>    inverse1 :: a -> a

> newtype Z = Z Int deriving (Eq,Ord,Show,Arbitrary,Num)

> instance Group1 Z where
>    unit1 = Z 0
>    mult1 = (+)
>    inverse1 = negate

Now we can try testing some of the group laws. For example

> test_inverse1 a = mult1 a (inverse1 a)
> ex1 = quickCheck (\a -> test_inverse1 a==(unit1::Z))

If we run ex1 we gain some confidence that we have a group. If we want to be more categorial then we need to express as much as possible as compositions of functions. So let's break down test_inverse1 as a composition of elementary operations.

> diag a = (a,a)
> both f g (a,b) = (f a,g b)
> test_inverse2 = uncurry mult1 . (id `both` inverse1) . diag
> ex2 = quickCheck (\a -> test_inverse2 a==(unit1::Z))

There's a nice symmetry in test_inverse2. If we make diag part of the definition of a group that we can make this highly symmetric definition of a group:

> class Group2 a where
>     unit2 :: () -> a
>     counit2 :: a -> ()
>     mult2 :: (a,a) -> a
>     comult2 :: a -> (a,a)
>     anti2 :: a -> a

I've even made a (useless looking) symmetric partner for unit1 and I've renamed the inverse operation as anti, short for antipode. Let's make Z an instance of this class.

> instance Group2 Z where
>   unit2 () = Z 0
>   counit2 _ = ()
>   mult2 = uncurry (+)
>   comult2 a = diag a
>   anti2 = negate

I hope I haven't done anything too tricky yet. This is just the familiar definition of a group except I've made a trivial 'counit' and the diagonal function part of the definition. Now comes the most important step in this post: I'll tweak this definition so that all of the associated functions are now Kleisli arrows:

> class Group m a where
>     unit :: () -> m a
>     counit :: a -> m ()
>     mult :: (a,a) -> m a
>     comult :: a -> m (a,a)
>     anti :: a -> m a

We can now straightforwardly reimplement Z as an instance of this class.

> instance Monad m => Group m Z where
>     unit _ = return (Z 0)
>     counit _ = return ()
>     mult = return . uncurry (+)
>     comult = return . diag
>     anti = return . negate

Apart from writing it monad-agnostically, it's just like before. We can exactly recover what we had before with the help of the identitiy monad.

> newtype Identity a = I a deriving (Eq,Ord,Show)

> instance Monad Identity where
>     return x = I x
>     I x >>= f = f x

Now we can rewrite the group laws. Before checking to see if aa^-1=1 we'll make sure aa^-1=a^-1a.

> test_antipode1 a = do
>     (u,t) <- comult a
>     u' <- anti u
>     mult (u',t)

> test_antipode2 a = do
>     (u,t) <- comult a
>     t' <- anti t
>     mult (u,t')

> ex3 = quickCheck (\a -> lhs a==rhs a) where
>   lhs = test_antipode1 :: Z -> Identity Z
>   rhs = test_antipode2 :: Z -> Identity Z

For aa^-1 we need to compare with the identity. We make an identity by applying unit to (). We could just provide a () but we get something more symmetrical if we construct it using the counit.

> test_antipode3 a = do
>   x <- counit a
>   unit x

> ex4 = quickCheck (\a -> lhs a==rhs a) where
>   lhs = test_antipode1 :: Z -> Identity Z
>   rhs = test_antipode3 :: Z -> Identity Z

We can test associativity easily

> test_assoc1 (a,b,c) = do
>     ab <- mult (a,b)
>     mult (ab,c)

> test_assoc2 (a,b,c) = do
>     bc <- mult (b,c)
>     mult (a,bc)

> ex5 = quickCheck (\a -> lhs a==rhs a) where
>   lhs = test_assoc1 :: (Z,Z,Z) -> Identity Z
>   rhs = test_assoc2 :: (Z,Z,Z) -> Identity Z

In the interest of symmetry, there ought to be a coassociativity law for comultiplication. Here it is:

> test_coassoc1 x = do
>     (u,v) <- comult x
>     (s,t) <- comult v
>     return (u,s,t)

> test_coassoc2 x = do
>     (u,v) <- comult x
>     (s,t) <- comult u
>     return (s,t,v)

> ex6 = quickCheck (\a -> lhs a==rhs a) where
>   lhs = test_coassoc1 :: Z -> Identity (Z,Z,Z)
>   rhs = test_coassoc2 :: Z -> Identity (Z,Z,Z)

One issue is that we've introduced comult without making much demand on it. Here's a law that is obviously true in any group that helps to pin down comult. It's simply evaluating (p,q) →(pq,pq) in two different ways.

> test_multcomult1 (u,v) = do
>     (p,q) <- comult u
>     (r,s) <- comult v
>     u' <- mult (p,r)
>     v' <- mult (q,s)
>     return (u',v')

> test_multcomult2 (u,v) = do
>     w <- mult (u,v)
>     comult w

> ex7 = quickCheck (\a -> lhs a==rhs a) where
>   lhs = test_multcomult1 :: (Z,Z) -> Identity (Z,Z)
>   rhs = test_multcomult2 :: (Z,Z) -> Identity (Z,Z)

Out of laziness I'm not going to test every group law, but shortly I'll show you where to find all of the laws, including the new laws for comult.

Hopf Algebras

Identity Z is just the usual group of integers with addition. But the time has come to switch monads. What's H=V Float Z? It's clearly a vector space with basis elements labelled with integers. It's also an instance of Group (V Float) meaning it comes equipped with implementations of functions like unit and comult. But note this: counit is no longer trivial. We have counit :: Z -> V Float (). In other words counit lifts to a linear map from H to V Float (), which is isomorphic to Float. In fact, it sums the coefficients if the basis elements. unit now embeds the ground field in our new structure. But what's mult? Well V Float (Z,Z) is H⊗H. So mult lifts to a bilinear map from H×H to H. We have a vector space with a bilinear multiplication law. Is this law associative? We could write code to check. But here's the really cool bit: we can just reuse test_assocn from above.

> ex8 = quickCheck (\a -> lhs a==rhs a) where
>   lhs = test_assoc1 :: (Z,Z,Z) -> V Float Z
>   rhs = test_assoc2 :: (Z,Z,Z) -> V Float Z

By switching monad, we automatically get a new algebraic structure and new 'semantics' for the group laws. Unpacking ex8 reveals that it gives the associativity law for an algebra. Unpacking further group laws in this context gives bialgebra laws.

> ex9 = quickCheck (\a -> lhs a==rhs a) where
>   lhs = test_multcomult1 :: (Z,Z) -> V Float (Z,Z)
>   rhs = test_multcomult2 :: (Z,Z) -> V Float (Z,Z)

And 'monadifying' the laws about group inverses give rise to the Hopf algebra laws for the antipode:

> ex10 = quickCheck (\a -> lhs a==rhs a) where
>   lhs = test_antipode1 :: Z -> V Float Z
>   rhs = test_antipode2 :: Z -> V Float Z

In other words, a Hopf algebra over k is nothing other than a group in the generalised sense above using the monad V k. The underlying object and all of the Hopf algebra laws just follow automatically as soon as you switch monads. Groups and Hopf algebras are the same thing.

But we still have this example Hopf algebra V Float Z to think about. It's none other than the group algebra of Z. It's essentially the algebra of Laurent polynomials. Clearly we could repeat this group algebra construction for any group and Z was just an example.

Examples

Now I have this machinery for working with Hopf algebras I might as well construct some examples that aren't simply generated by some underlying group. I found some pretty combinatorial examples in this paper.

Firstly I'd like to define some combinatorial operations on lists: cat, chop, shuffle and excise.

cat simply joins a pair of lists.

> cat x = return $ uncurry (++) x

chop computes a kind of inverse of cat. This inverse is multivalued so depending on which monad is used it can return all of the possible pairs that could have been catted to make a particular list.

> chop :: (MonadPlus m, Functor m) => [a] -> m ([a], [a])
> chop [] = return ([],[])
> chop (a:as) = return ([],a:as) `mplus` fmap ((a:) `both` id) (chop as)

shuffle finds all the ways that two lists can be riffle shuffled together.

> shuffle ([],b) = return b
> shuffle (a,[]) = return a
> shuffle (a:as,b:bs) =
>     fmap (a:) (shuffle (as,b:bs)) `mplus`
>     fmap (b:) (shuffle (a:as,bs))

excise is the multivalued inverse of shuffle in that it finds all of the ways pairs of lists could have been riffle shuffled together to give a particular list.

> excise [] = return (mzero,mzero)
> excise (a:as) = do
>     (x,y) <- excise as
>     return (a:x,y) `mplus` return (x,a:y)

We can now define a pair of Hopf algebras as follows.

> newtype Hopf1 = H1 [Int] deriving (Ord,Eq,Show,Arbitrary)

> instance Group (V Float) Hopf1 where
>     unit _ = return (H1 [])
>     counit (H1 []) = return ()
>     counit _ = zero
>     mult (H1 x,H1 y) = fmap H1 $ cat (x,y)
>     comult (H1 x) = fmap (\(x,y) -> (H1 x,H1 y)) $ excise x
>     anti (H1 x) = (-1)^(length x) .* return (H1 $ reverse x)

> newtype Hopf2 = H2 [Int] deriving (Ord,Eq,Show,Arbitrary)

> instance Group (V Float) Hopf2 where
>     unit _ = return (H2 [])
>     counit (H2 []) = return ()
>     counit _ = zero
>     mult (H2 x,H2 y) = fmap H2 $ shuffle (x,y)
>     comult (H2 x) = fmap (\(x,y) -> (H2 x,H2 y)) $ chop x
>     anti (H2 x) = (-1)^(length x) .* return (H2 $ reverse x)

We can now test some of the Hopf algebra laws

> ex11 = quickCheck (\a -> lhs a==rhs a) where
>   lhs = test_assoc1 :: (Hopf1,Hopf1,Hopf1) -> V Float Hopf1
>   rhs = test_assoc2 :: (Hopf1,Hopf1,Hopf1) -> V Float Hopf1

Unfortunately my code is a inefficient, and the shuffles cause combinatorial explosions, so it's probably best not to run these examples. This post is about expressing algebraic structures rather than making practical computations in them.

> ex12 = quickCheck (\a -> lhs a==rhs a) where
>   lhs = test_coassoc1 :: Hopf1 -> V Float (Hopf1,Hopf1,Hopf1)
>   rhs = test_coassoc2 :: Hopf1 -> V Float (Hopf1,Hopf1,Hopf1)

> ex13 = quickCheck (\a -> lhs a==rhs a) where
>   lhs = test_antipode1 :: Hopf1 -> V Float Hopf1
>   rhs = test_antipode2 :: Hopf1 -> V Float Hopf1

> ex14 = quickCheck (\a -> lhs a==rhs a) where
>   lhs = test_multcomult1 :: (Hopf1,Hopf1) -> V Float (Hopf1,Hopf1)
>   rhs = test_multcomult2 :: (Hopf1,Hopf1) -> V Float (Hopf1,Hopf1)

And of course you can try the same with Hopf2.

Sweedler Notation

There's one last important thing to come out of these. We have a powerful notation for working with tensor products of vector spaces but the notation is so easy to use that it's hidden the fact that there's complex stuff going on under the hood. When we're dealing with the Identity monad then a line like (p,q) <- comult u is pretty straightforward to interpret as an assignment. But in the vector space monad it means something very different. Nonetheless, we can often reason (with care) as if it were just an assignment. This is what Sweedler Notation is about. But it kind of just appears automatically here.

Summary

What they say is right. Hopf algebras really are groups in disguise. Just take all of the functions that make up the structure of a group and convert each one to a Kleisli arrow for the vector space monad. You end up with Hopf algebras.

I'm amazed at the way I could express this level of abstraction in Haskell. I don't know of any computer algebra package that comes near to this level of expressivity. Although algebra packages often come with amazing implementations of algebraic algorithms, the programming languages they come with are generally completely uninteresting in themselves. I'm looking forward to the day when computational algebraists and programming language theorists get talking to each other.

By the way, this kind of manipulation was what monads were invented for. It's nice to see the two different applications work so well together.

The Fibonacci Numbers, Coalgebraicaly

> {-# LANGUAGE TypeSynonymInstances #-}

What kind of Haskell programmer can I call myself if I haven't written a fiendishly complex piece of code to inefficiently compute the Fibonacci numbers. Unfortunately, what follows isn't particularly fiendish and isn't wholly inefficient, but it'll have to do.

My real motivation is this: I've mentioned before that one difference between mathematics and computer science is that mathematics has more of a bias towards algebraic rather than coalgebraic structures. But coalgebras do appear in mathematics. For example, every Hopf Algebra is a coalgebra (and an algebra) and Hopf algebras play a prominent role due to the amazing discovery that they have a nice connection with combinatorics, especially helping to organise the combinatorics that come up in renormalisation.

The first thing to note is that the computer science and mathematics definitions of coalgebra are closely related but different. An algebra is essentially a vector space over some field K with a bilinear associative product and multiplicative identity. If our algebra is V, then we have a multiplication map V×V -> V that is linear in each of its arguments. That's the same as saying that multiplication is a linear map V⊗V->V. In a coalgebra we have the converse, a linear comultiplication map Δ:V→V⊗V. Similarly, a unit in an algebra gives a way to map your ground field into the algebra, you map x to x.1 so a counit gives a map from V down to K. The full definition is on wikipedia.

The Fibonacci sequence can be thought of as a function N→R, where

f(0)=1
f(1)=1
f(n)=f(n-1)+f(n-2) for n>=2

> f 0 = 1
> f 1 = 1
> f i = f (i-1)+f (i-2)

Many functions satisfy sum formulae like e^x+y=e^xe^y and tan(x+y)=(tan(x)+tan(y))/(1-tan x tan y). f also has one. What's more, the formula is essentially the definition of a coproduct on a coalgebra.

Define a Fibonacci-like sequence to be a sequence s satisfying just the last part of the definition of f:

s(n)=s(n-1)+s(n-2) for n>=2

Because I have chosen maps to R, Fibonacci-like sequences form a vector space over R.

Let S be the space of Fibonacci-like sequences. We're going to try to make this into a coalgebra. It's easy to come up with a candidate for a counit: ε(f) = f(0).

We can express this in Haskell:

> class VectorSpace v where
>    zero :: v
>    (.*) :: Float -> v -> v
>    (.+) :: v -> v -> v
>    (.-) :: v -> v -> v
>    v .- w = v .+ ((-1) .* w)

> type Sequence = Int -> Float

> instance VectorSpace Sequence where
>    zero _ = 0
>    a .* v = \i -> a*v i
>    v .+ w = \i -> v i+w i

> epsilon f = f 0

We now need to consider S⊗S. Tensor products of infinite dimensional spaces can be a bit delicate. But S is a two-dimensional space because it is completely specified by its first two terms. This means that S⊗S is the obvious thing: a 2D grid of values where each row and each column is an element of S.

> type Sequence2 = Int -> Int -> Float

> instance VectorSpace Sequence2 where
>    zero _ _ = 0
>    a .* v = \i j -> a*v i j
>    v .+ w = \i j -> v i j+w i j

The tensor product operation is straightforward:

> (<*>) :: Sequence -> Sequence -> Sequence2
> v <*> w = \i j -> v i*w j

Let's define the operator D so that (Df)(i)=f(i+1). So D is a kind of shift operator. Dⁿ shifts n times and we can form polynomials in D. Note that this space of polynomials, P, in D is also two-dimensional because D²=1+D by definition of Fibonacci-like.

> d :: Sequence -> Sequence
> d v i = v (i+1)

We can define a pairing between P and S by <f(D),s> = ε(f(D)s). This makes P dual to S. P has an obvious multiplication inherited from ordinary multiplication which we can denote by the map m:P×P→P. This suggests the idea of forming a kind of dual to multiplication, Δ:S→S⊗S, defined by

<f(D)g(D),s>=<f(D)⊗g(D),Δ(s)>

The fact that we have defined things through this duality ensures that (S,ε,Δ) is a coalgebra. (Note that the pairing extends to the tensor products via <f(D)⊗g(D),s⊗t>=<f(D),s><g(D),t>

By definition <D^m⊗Dⁿ,Δ(s)>=<D^m+n,s>. So Δ(s) is simply the 2D grid made by arranging s along the edges and filling in with constant values along the diagonals.

> delta :: Sequence -> Sequence2
> delta v = \i j -> v (i+j)

Give Δ(s) we can extract s back out from it in many different ways. To get the ith element we need to pick j and k such that j+k=i and grab the ith element of the kth row. One choice is to pick k such that it is half of i rounded down:

> extract :: Sequence2 -> Sequence
> extract f i = let k = i `div` 2 in f k (i-k)

Now given a Fibonacci-like sequence s, we can think of the top left 2×2 block of Δ(s) as a matrix. Ie. s(i+j) for i=0,1;j=0,1 forms a matrix. It's clearly invertible by inspection. But there's a deeper reason why it's invertible. If there were a linear dependence between the rows of this matrix it would be a relationship between the shifts of the sequence and hence would give rise to a simpler definition of the Fibonacci sequence. For more general recursive sequences, if there is an nth order recurrence relation for the sequence, and no simpler relation, then the top left (n-1)⊗(n-1) block is full rank and invertible.

Write t_ij for the inverse of this matrix. Then

Σ_jf(i+j)t_jk=δ_ik for i,j,k=0,1.

Now here's a neat bit: any Fibonacci-like sequence is defined by its first n-1 elements. Any Fibonacci-like sequence is simply a linear combination of f and Df. So the above is actually true for any Fibonacci-like sequence. And as any shift of the Fibonacci sequence is itself a Fibonacci-like sequence, the above is true for any i,k>0.

Another matrix multiplication gives

Σ_j,kf(i+j)t_jkf(k+l)=f(i+l) for j,k=0,1;i,l≥0

That's our Fibonacci summation formula. It tells us how to get f(i+l) in terms of f(i+j) and f(k+l) for j,k=0,1.

Using the definitions above we get

Σ_j,k<Dⁱ⊗D^l,t_jkD^jf⊗D^kf=e<D^i+l,f>

So t_jkD^jf⊗D^kf satisfies the definition of Δ above, and so must in fact be the same thing. So we have:

Δ(f)=2s⊗f-f⊗Df-Df⊗f+Df⊗Df

(Where I actually computed the elements of the inverse matrix t.)

We can now use this as the basis for a recursive algorithm. We need the bootstrap function to get the recursion started:

> bootstrap :: [Float] -> Sequence -> Sequence
> bootstrap xs v i = if i<length xs then xs!!i else v i

> f' :: Sequence
> f' = bootstrap [1,1,2,3] $ extract $ 2 .* (f' <*> f') .- (f' <*> (d f')) .- ((d f') <*> f') .+ ((d f') <*> (d f'))

It's faster than the usual recursion. (Though that's not saying much.)

I'll leave it to you to extract out what's really going on above, and to generalise it to give sum formulae for other recurrence relations.

Note that the goal here isn't to write a fast algorithm. My main motivation is to view the sum formula as an identity about comultiplication and get some intuition about comultipliation so I can read some Hopf algebra papers.

> let's_go_slow = map f [0..]
> let's_go_faster = map f' [0..]

All of this, apart from the algorithm itself, is an unpacking and simplification of section 1 of Bialgebras of Recursive Sequences and Combinatorial Identities by Futia et al.

PS Tell me if you have trouble reading the above. I wrote a new pseudo-LaTeX to blogger-HTML converter and it might not actually be any good.

MSFP 2008

I was invited by Conor McBride to give a talk MSFP 2008 in Reykjavik. I've put the slides online here. While rehearsing I wrote a 'script' that I've put here. That's not actually what I said but I felt I needed to write something out beforehand just so I could make sure I had an appropriate amount of material. Apologies for the silly chatty style in a written document but that's the sort of thing that came to mind while sitting in a hotel room. I got to meet in person a bunch of people whose papers I'd read including Thorsten Altenkirch, Tarmo Uustalu and Jeremy Gibbons.

I had an opportunity to explore around Reykjavik so I also put some photos online, including a terrible one of Conor and Tarmo Uustalu (sorry!). By curious coincidence I bumped into, and got to hang out with, an old work colleague as well as a friend of the person who sits one office down from me at work, and Eben Moglen. Small world.

A blessed man's formula for holey containers

I love the way derivatives of types tell you about holes in containers. It works the other way too and holes can give insight into derivatives.

Suppose S and T are containers so that S(X) and T(X) are containers of elements of type X. Then S(T(X)) is an S-container of T-containers of X's. If we draw S as a square (rectangle actually) and T as a triangle then we can draw a picture of an example of such a thing:


Taking the nth derivative of a type gives the type with (an ordered sequence of) n holes. Here's that previous container with 3 holes made in it:

As the holes have an ordering I've numbered them from 1 to 3.

An S-container of T-containers with holes is essentially an S-container containing both ordinary T-containers, and some T-containers with holes. If we excise the T-containers with holes, we're left with an S-container containing just T-containers, and some holes. Here's a picture of that:

But there are a couple of problems with that. Once we've excised the holey sub T-containers we don't know which S-holes to plug them back into and we don't know how to reconstruct the original numbering of the holes. We need to keep a tiny bit more information. That's the set I wrote down: {{1,3},{2}}. Call each element of this set a block. I've written the elements of the blocks in ascending order and I've written the blocks in ascending order of their lowest elements. The first block corresponds to hole 1 in the S-container and the second block corresponds to hole 2. Similarly, we write the elements of the blocks into the holes in the T-containers. And that allows us to reconstruct the original ST-container.

Think about the general case. We make n holes in an ST-container. So up to n of the T-containers acquire holes. Let's say it's m holes. We can think of the S-container having m holes and these holes being filled by T-containers with b₁, b₂,...,b_m holes where the sum of the b_i is n. We're essentially just partitioning the original n holes {1,2,...,n} into m sets. So each ST-container with n holes gives a partition of {1,...,n}, an S-container with n holes, and m T-containers where the ith has b_i holes. Writing containers with n holes as nth derivatives we get

dⁿ(F(G(X))/dXⁿ = sum over each partition P of {1,...,n} of d^|P|F/dG^|P| times the product over each block B in P of d^|B|F/dX^|B|

Note that the equality above isn't just a numerical equality, it's an isomorphism of types. In fact, it's the type version of the 'combinatorial form' of the Faà di Bruno formula. Although that wikipedia page describes the formula as 'forbidding', if I've done my job right then I think the above picture make it seem almost trivial. I find it much easier to think of this version of the chain rule in terms of holes.

And the explanation for the title of this post: in 1988 Pope John Paul II beatified Faà di Bruno.

Categories of polynomials and comonadic plumbing

Suppose you have a Haskell program and you want to introduce a new global constant into your program. There are at least two approaches you could take:

1. Simply introduce a new global constant. You could name it x and write something like x=1.23456 and refer to x throughout your code. This has the advantage of being easy to implement.
   
2. Write all of your code in monadic style and make use of the reader monad. This is intrusive in the sense that you may have to make many changes to your code to support it. But it has the advantage that all of your functions now explicitly become functions of your global constant.
   

Now I’m going to roughly sketch a more categorical view of both of these approaches. So let’s restrict ourselves to the subset of Haskell that corresponds to typed lambda calculus without general recursion so that we know all of our functions will be total and correspond to the mathematical notion of a function. Then all of our functions become arrows in the category that we’ll call Hask.

Firstly consider approach (1). Suppose we want to introduce a new constant, x, of type A. Category theory talks about arrows rather than elements of objects, so instead of introducing x of type A, introduce the function x:1->A where 1 is the terminal object in Hask, normally called (). An element of A is the same thing as an element of 1->A, but in the latter case we have an arrow in the category Hask.

Before continuing, let me digress to talk about polynomials. Suppose we have a ring (with an identity) R. We define R[x], where x is an indeterminate, to be the ring of polynomials in x. Another way to describe that is to say that R[x] is the smallest ring containing R and an indeterminate x, that makes no assumptions about x other than those required to make R[x] a ring. For example we know that (1+x)(1-x)=1-x², because that must hold in any ring. Given a polynomial p in R[x] we can think of it as a function f_p from R to R. f_p(a) is the value we get when substituting the value of a for x in p. So a polynomial in R[x] is the same as a function from R to R that can be written in terms of elements of R, multiplication and addition.

We can do the same with category theory. Given a category A we can ask for the smallest category extending A and containing an indeterminate arrow x:1 -> A. Just as with polynomials we have to allow all possible arrows that can be made by composing arrows of A with x. The resulting expressions for arrows will contain x as a free variable, just like the way x appears in polynomials. In fact, by analogy we can call the resulting category, A[x], the category of polynomials in x:1->A. In the special case A=Hask, you can see that Hask[x] is the category of Haskell functions extended by a new constant of type x:1->A but assuming no equations other than those necessary to make Hask[x] a category. Just as an arrow in Hask is a Haskell function, an arrow in Hask[x] is a Haskell function making use of an as yet undefined constant x.

(I've glossed over some subtleties. Just as we need a suitable equivalence relation to ensure that (1+x)(1-x)=1-x² in R[x], we need suitable equivalence relations in our category. I'll be showing you where to find the missing details later.)

Here's the implementation of a function, h, making use of a constant x:

(Note that I'll be using Edward Kmett's category-extras shortly so I need some imports)

> import Control.Monad.Reader
> import Control.Comonad
> import Control.Comonad.Reader

> x = 1.23456

> f a = 2*a+x
> g a = x*a
> h a = f (g a)

> test1 = h 2

Now consider the second approach. The easiest thing is to just give an implementation of the above using the reader monad:

> f' a = do
>    x <- ask
>    return $ 2*a+x

> g' a = do
>    x <- ask
>    return $ x*a

> h' a = return a >>= g' >>= f'

> test2 = runReader (h' 2) 1.23456

Note how, as is typical in monadic code, I have to plumb f' and g' together using >>= so that 1.23456 is passed through f' and g'. Previously I've described another way to think about the composition of monadic functions. Using >>= we can compose functions of type a->m b and b->m c to make a function of type a->m c. The result is that given a monad we can form the Kleisli category of the monad. The objects are the same as in Hask, but an arrow from a->b in the Kleisli category is an arrow of type a->m b in Hask. It's not hard to show this satisfies all of the axioms of a category. When we program in the reader monad it's a bit like we've stopped using Hask and switched to the Kleisli category of the reader monad. It's not quite like that because we used functions like +. But in theory we could use lifted versions of those functions too, and then we'd be programming by composing things in the Kleisli category. If we call the reader monad R then we can call the corresponding Kleisli category Hask^R. (Strictly speaking that R needs a subscript telling is the type of the value we intend to ask for.)

So here's the important point: Hask[x] is the same category as Hask^R. In both cases the arrows are things, which when supplied a value of the right type (like 1.23456), give arrows in Hask from their head object to their tail object.

But there's another way to do this. We can use the reader comonad:

> f'' a = 2*extract a+askC a
> g'' a = extract a*askC a
> h'' a = a =>> g'' =>> f''

> test3 = runCoreader (h'' (Coreader 1.23456 2))

In a similar way, we're dealing with arrows of the form wa -> b and we can compose them using =>>. These arrows form the coKleisli category of the reader comonad, S, which we can write Hask_S. So we must have

Hask[x]≅Hask^R≅Hask_S.

Now some back story. Over 20 years ago I was intrigued by the idea that logic might form a category with logical ‘and’ and ‘or’ forming a product and coproduct. I came across the book Introduction to Higher Order Categorical Logic by Lambek and Scott for ₤30.00. That’s ₤60.00 at today's prices, or about $120.00. On a student grant? What was I thinking? And as it bore no relation to anything I was studying at the time, I barely understood a word of it. I was probably fairly applied at that point doing courses in stuff like solid state physics and electromagnetism as well as a bit of topology and algebra. I doubt I'd heard of lambda calculus though I could program in BASIC and APL. So there it sat on my bookshelf for 22 years. Periodically I’d look at it, realise that I still didn’t understand enough of the prerequisites, and put it back on the shelf. And then a month or so ago I picked it up again and realised that the first third or so of it could be interpreted as being about almost trivial Haskell programs. For example, on page 62 was

Proposition 7.1
The category A[x] of all polynomials in the indeterminate x:1->A over the cartesian or cartesian closed category A is isomorphic to the Kleisli category A_A=A_SA of the cotriple (S_A,&epsilon_A,δ_A).

The language is a little different. Lambek and Scott used the term cotriple instead of comonad and Kleisli category where I’d say coKleisli category. δ and ε are cojoin and coreturn. And Lambek and Scott's theorem applies to any cartesian closed category. But after staring at this claim for a while it dawned on me that all it was really saying was this: here are two ways to introduce new constants into a category. But there’s no way I would have seen that without having practical experience of programming with monads. Learning Haskell has finally paid off. It’s given me enough intuition about category theory for me to get some return on my ₤30.00 investment paid to Heffers all those years ago. I expected to take this book to my deathbed, never having read it.

Anyway, for the details I left out above, especially the correct equivalence relation on Hask[x], you'll just have to read the book yourself.

Also, note the similarity to the deduction theorem. This theorem says that if we can prove B, assuming A, then we can deduce A implies B without making any assumptions. It unifies two way to introduce a proposition A, either as a hypothesis, or as an antecedent in an implication. In fact, the above theorem is just a categorical version of the deduction theorem.

Also note the connection with writing pointfree code. In fact, the pointfree lambdabot plugin makes use good use of the reader monad to eliminate named parameters from functions.

I’m amazed by seeing a book from 1986 that describes how to use a comonad to plumb a value through some code. As far as I know, this predates the explicit use of the reader monad in a program, Wadler and Moggi’s papers on monads, and certainly Haskell. Of course monads and comonads existed in category theory well before this date, but not, as far as I know, for plumbing computer programs. I’d love to hear from anyone who knows more about the history these ideas.

Drugs, Kate Moss, and Intuitionistic Logic

Before make the post I promised I'd make I thought I'd make a digression to point out a connection between Intuitionistic Logic and a recent news story.

Intuitionistic logic is what we get when we take ordinary everyday classical logic and drop the Law of the Excluded Middle (LEM), in other words, we drop the law that says that for any proposition, either it or its negation is true. This is such an ingrained notion that it's hard to imagine giving it up. For example, it seems obvious that either it's raining or it's not raining. But there are good reasons for not taking it for granted in mathematics. The main issue is that sometimes when we use LEM we find that we can prove the existence of something, but have no way of constructing it. Mathematics is full of such proofs. The best known example is probably the proof that an irrational number raised to the power of an irrational number can be rational. Some mathematicians (and many computer scientists) only like constructive proofs, ie. proofs that actually exhibit the thing whose existence is being proved. Using intuitionistic logic is a good way to always force this to be true. So, for example, in classical logic, you might find you have proved "A or B", but not have a proof of A or a proof of B. But in intuitionistic logic, if you prove "A or B", then you must have proved A or you must have proved B.

Now (British) criminal law is more lax than mathematical logic when it comes to proof. You only need to prove something beyond reasonable doubt, rather than providing a rigorous derivation from axioms. (Though admittedly mathematicians rarely do this in practice.) But criminal law does have one place where it has higher standards of proof than classical logic: you can't necessarily convict someone of "A or B" unless you have a proof of A or a proof of B. A proof of "A or B" will not do.

It seems that there is enough evidence to show that Kate Moss recently used a controlled substance of class A or class B. Unfortunately, the law requires either a proof that she had been using a class A drug, or a proof that she had been using a class B drug. A proof that she had been using one or the other will not do. And hence Kate Moss cannot be prosecuted. In this regard, the law is Intuitionistic. If you're going to have a legal system that doesn't recognise LEM then you really need to carve up crime-space as a semilattice so you can charge people with the join of two crimes :-)

Life in a Lazy Universe

I've seen lots of articles on simulated reality recently. Supposing for a bit that this is a sensible hypothesis to entertain and not just a science fictional conceit, what might be the consequences of supposing that such a 'reality' is running on a machine with lazy evaluation?

If your task were to simulate a universe well enough for its inhabitants to be convinced it was seamless, then there is an obvious optimisation that could be made. You could simulate just those parts that could be perceived by the inhabitants. But the catch is that if an inhabitant were to explore a new region, the simulation would be required to fill in that region. Just creating a blank new area wouldn't do, it's current state would need to be consistent with having had a plausible history and so the simulation would be required to fill in not just its current state, but also its past. This is precisely what is provided by lazy evaluation - demand driven evaluation that takes place only when results are observed. It seems natural that such a simulation should make use of laziness.

But such optimisation doesn't need to be limited to lazily computing 'new regions' where by 'region' we mean a spatially localised volume. We could also imagine implementing level of detail. If we don't look closely at an object in our grasp, there's no need to compute every single detail of it. We need only provide detail at the resolution at which it can be perceived. We'd like this built into every aspect of the simulation so that anything within it is computed only to a desired accuracy. I've spoken a few times about such a data structure in posts like this. This is precisely what the real numbers give us. Computationally, a real number is an object which when given a degree of accuracy, returns a rational number (or some other finitely manipulatable structure) that represents the real to this accuracy. The wikipedia article suggests that 'the use of continua in physics constitutes a possible argument against the simulation of a physical universe'. This is diametrically opposed to what I might argue: in fact the presence of continua suggests the existence of an efficient demand driven simulation with level of detail. The problem is that the people who have been thinking about these issues have been thinking in terms of traditional imperative style programming where a real number is typically represented to some fixed finite precision. But in fact there are infinitely many real numbers that can be represented exactly if we think in terms of lazy data structures.

But this is all idle speculation if we don't make predictions. So here is one. In a computer simulated universe, all physics must obviously be computable. But all computable functions are continuous. So in a simulated universe we'd expect the physical laws involving real numbers to make use of them only in a continuous way. Sound familiar?

---

The next section is an aside and I need to give a spoiler alert as I'll be mentioning some aspects of Greg Egan's Permutation City.

What exactly what is meant by 'demand driven' above? If we have implemented a simulated universe I expect we would be interested in looking in on it from time to time. So this is what is usually meant by 'demand'. Whenever the program performs I/O to show us what is going on, it would trigger the evaluation of the lazy thunks that had been sitting around. But consider how the universe might appear to its inhabitants. Whether or not we look in on it us surely irrelevant to the inhabitants, assuming we just look and don't try to communicate. But if I/O isn't performed, functional programmers don't make much of a distinction between an unevaluated thunk and an evaluated one. They are observationally equivalent. So do we in fact need to run the simulation for its inhabitants to consider themselves to exist? Anyway, I won't say anything more about this because Permutation City deals with this issue at length. I'm just rephrasing it as a lazy evaluation issue.

My next post will be on life in a monadic universe.

---

Update: I'm going to add some detail to the above. Suppose we represent real numbers as infinite sequences of digits. Suppose also that we have some module in a simulation that, let's say, resolves collisions between particles. You might ask it "what is the velocity, to 20 digits, of particle A after the collision". It must then look at all of the relevant inputs, decide how accurately it needs them, and then propagate new requests. For example it might work like this "if I need the velocity of A to 20 digits I need the mass of particle B to 25 digits and the velocity of particle C to 26 digits..." and so on. The idea of demand driven modules that produce outputs to a given accuracy by making requests of their inputs to a given accuracy exactly parallels the subject of analysis in mathematics. There you ask questions like "if I need to know f(x,y) with an error of no more than δ, what ε of error can I tolerate in x and y?". Functions that allow this kind of analysis are precisely the continuous ones. In other words, a simulated universe that is built from continuous functions lends itself to a demand driven implementation. One might argue that the whole of the fields of analysis (and topology) are the study of functions that can be evaluated in this way. And of course, in our universe, the functions we study in physics are typically the continuous functions. So contrary to the claim that the existence of continua argues against the idea that the universe is a simulation, I'd like to point out that they might actually make it more convenient to simulate our universe. In fact, we use this fact all the time in ordinary everyday physics simulations because we know we only need to work to a certain level of accuracy to get accurate results.

The claim that this universe can be simulated, and the claim that this universe is a simulation, are philosophical claims that I don't like to waste time on. But the claim that the existence of continua blocks the notion that the universe is a simulation is a fairly precise statement of logic, mathematics and physics. I claim it is wrong, and that's the only claim I'm making here.

The Interchange Law

I find it amazing the way you can take a bunch of abstract nonsense and turn it into concrete Haskell code for doing something mundane. Having down-to-earth examples really does help with getting your head around the abstractions.

Anyway, on page 43 of Categories for the Working Mathematician (2nd Ed.) is the 'interchange law' for horizontal and vertical composition of natural transformations. Not only can we find a nice mundane example to help thinking about it, we can even test it using Quickcheck.

For convenience we can define a natural transformation type:

> {-# OPTIONS -fglasgow-exts -fno-warn-missing-methods #-}

> import Test.QuickCheck
> import Control.Monad.Writer

> type Natural f g = forall a.f a -> g a

Where f and g are intended to be instances of Functor.

If we have a natural transformation from f to g, and another from g to h, then we can compose them using the ordinary Haskell composition operator (.). In Haskell, functors are frequently containers, so we can use containers as a guiding example. If we know how to map a bag of stuff to a sack of stuff, and we know how to map a sack of stuff to a box of stuff, then we just perform the two operations in sequence and we can map a bag of stuff to a box of stuff. This is known as vertical composition.

But there is also another way to compose natural transformations known as horizontal composition. Again we can think about containers. Suppose we know how to map a bag of stuff to a sack of stuff, and we know how to map a box of stuff to a crate of stuff, then we also know how to map a bag of boxes of stuff to a sack of crates of stuff. There are two obvious ways to do this: we could convert the bag of boxes of stuff to a sack of boxes, and then convert each box to a crate. Or we could convert the boxes to crates first, and then map the resulting bag of crates to a sack of crates. We can define two binary operators `o` and `o'` to perform each of these tasks:

> o,o' :: (Functor f,Functor f',Functor g,Functor g') => Natural f' g' ->
>          Natural f g -> (forall c.f' (f c) -> g' (g c))
> o s t x = s (fmap t x)
> o' s t x = fmap t (s x)

Intuitively we'd expect these things to be equal, and in fact they always are. (Exercise: write a quickcheck for this based on the code below.)

Now we need some functorial containers to play with:

> data Pair x = Pair x x deriving (Eq,Show)
> instance Functor Pair where
>     fmap f (Pair a b) = Pair (f a) (f b)

> newtype Id x = Id x deriving (Eq,Show)
> instance Functor Id where
>    fmap f (Id x) = Id (f x)

Now we can define a bunch of natural transformations mapping between these containers and some others:

> alpha :: Natural Pair []
> alpha (Pair x y) = [x,y]

> beta :: Natural [] Maybe
> beta [] = Nothing
> beta (x:xs) = Just x

> alpha' :: Natural ((,) a) Id
> alpha' (a,x) = Id x

> beta' :: Natural Id (Either b)
> beta' (Id x) = Right x

On page 43 is the interchange law. Superficially it looks a lot like abiding. For any α, β, α' and β' with compatible types, we have

(β . α) o (β' . α') = (β o β') . (α o α')

This is the identity I want to check for the special case of the types I've chosen above. So we can define the left and right hand sides:

> lhs = (beta . alpha) `o` (beta' . alpha')
> rhs = (beta `o` beta') . (alpha `o` alpha')

> type From = Pair (Float,Integer)
> type To = Maybe (Either String Integer)

And here we go:

> test1 = quickCheck (\x -> lhs (x :: From) == (rhs x :: To))

Just type test1 in ghci to hear the sweet sound of 100 tests passing.

Anyway, we're just a hair's breadth away from defining 2-categories now. But I'll leave that for another day.

Thinking about it, I gave slogans for the previous theorems so why don't I give one for the interchange law. Take a deep breath. If you know how to convert a bowl of stuff into a box of stuff, and a box into a bag, and an urn into a crate, and a crate into a sack, then there are two equivalent ways to convert a bowl of urns into a bag of sacks: we can either construct methods going from bowl to bag and from urn to sack and combine them to go from bowl of urns to a bag of sacks, OR, we can go from a bowl of urns to a box of crates to a bag of sacks. I'm sure it's all a lot clearer now. ;-)

Well that's enough abstract nonsense. Time to get away from it all and watch some youtube videos. Eh? Oh well. Just wish me luck for tomorrow when I try to beat my fastest time in Bay to Breakers. Maybe some other bayfpers will be there too.

---

> instance Arbitrary a => Arbitrary (Id a) where
>     arbitrary = liftM Id arbitrary

> instance Arbitrary a => Arbitrary (Pair a) where
>     arbitrary = liftM2 Pair arbitrary arbitrary

Desingularisation and its applications

> {-# OPTIONS -fno-warn-missing-methods #-}

Here's a neat trick you can carry out with the Dual class that I've talked about many times before. Suppose you define a function like

> sinc x = sin x/x

Clearly this function has a singularity at x=0 where the division isn't defined. But this is an example of a removable singularity because it is possible to replace this function with another that is continuous (in fact, analytic) everywhere and is well defined at x=0. We could do this simply by defining:

> sinc' 0 = 1
> sinc' x = sin x/x

but it'd be cool if we could automate the process. Here's a function that does exactly that:

> desingularise f x = re (f (D x 1))

As if by magic, we can now compute

> ex1 = desingularise sinc 0

We can try other functions too:

> ex2 = desingularise (\x -> x/(exp x-1)) 0

So how does it work? Well I've modified my usual definition of (/) so that it is able to divide one infinitesimal number by another:

> instance Fractional a => Fractional (Dual a) where
>   D 0 a'/D 0 b' = D (a'/b') 0
>   D a a'/D b b' = D (a/b) ((b*a'-a*b')/(b*b))

desingularise 'perturbs' the value of its argument infinitesimally so that when we pass in an argument of zero we end up evaluating our function not at zero, but somewhere infinitesimally close. So our division by zero becomes division of infinitesimals and the above definition comes into play. The first clause of the definition can also be viewed as an application of l'Hôpital's rule.

But it doesn't solve all of our problems. For example

> ex3 = desingularise (\x -> sin x^2/x^2) 0

doesn't quite work. You need two applications of l'Hôpital's rule. We can implement this easily with

> ex4 = desingularise (desingularise (\x -> sin x^2/x^2)) 0

That's fine if you know how bad your singularity is in advance. If you don't, you can use power series instead. None of these methods will work for this function however:

> horrible x = exp (-1/x^2)

Anyway, the fact that it doesn't work for all functions doesn't detract from the fact that it's a pretty approach from solving the problem when you have a fairly well behaved singularity.

I've not actually seen anyone else use dual numbers to smooth over singularities in this way, but a few years ago I did come across something very closely related that apparently caused a little controversy in the computational geometry world.

Suppose you want to determine whether or not a point lies within a possibly non-convex and arbitrarily complicated polygon. One algorithm is this: follow a line in any direction from the point to infinity (though just outside the bounding box of the polygon would do) and count how many times you cross the boundary of the polygon. If you count an odd number of crossings you're in the polygon, otherwise you're outside. It's beautifully simple and if nobody had thought of it before it'd be worthy of publication.

But like the methods described in many computational geometry papers I've read, this one often fails in practice. There are at least two sources of difficulty. Firstly, computational geometry algorithms are subject to numerical errors. You often find yourself having to decide whether or not a point is on one side or another of a line, say, and getting inconsistent results because two different ways of computing the same thing result in different decisions being made because of rounding. That can be fixed by using something like arbitrary precision rational arithmetic. But even if you do this, there is still a problem.

When you follow a line from the point you're testing you may find that it doesn't simply cross one of the sides of the polygon, it may travel precisely through a vertex. If your code isn't written carefully you may accidentally count that as zero, one or two crossings. It's easy to write sloppy code that gets it right most of the time, but it's hard to work out every possible case where some kind of degeneracy might happen. What happens if the outward ray is actually collinear with a polygon face? Do vertex intersections and collinearities cover everything?

An ugly kludge is this: at the first sign of trouble, restart the algorithm but jiggle all the inputs randomly a bit. For a large class of algorithms this will work tolerably well, depending on your definition of tolerable. (Compare with Sard's theorem in Morse theory.) But it really is ugly, especially if you've gone through the hassle of using arbitrary precision arithmetic. But maybe you can see what's coming...

Back in '94 Edelsbrunner and Mücke came up with the idea of deforming the geometry infinitesimally, giving their method the utterly bizarre name "Simulation of Simplicity". I think it took people a while to realise that they were actually doing the same thing as people doing automatic differentiation and so the whole notion of an infinitesimal perturbation left people wondering what in fact they were doing. More seriously was this: writing robust computational geometry code that catches all of the degenerate cases is hard. How is it that such a simple algorithm could be slotted into a sloppily written polymorphic algorithm and turn into a watertight one? Was everyone in computational geometry going to be put out of work?

Luckily for the computational geometers it still takes quite a bit of work to get the details exactly right, and even when you do, the result isn't necessarily an efficient solution to the problem. I know that when I had a messy computational geometry problem to solve when I was working on MI:2 I just wrote code that turned a blind eye to degeneracies and aggressively pruned away anything that looked suspicious: like near zero length lines and near zero area faces. If a stray pixel poked through somewhere that it shouldn't, and tweaking the tolerance didn't fix it, an artist found a way to hide it. Sometimes beauty and elegance come second place to getting the job done on time.

---

> data Dual a = D a a deriving (Show,Eq)

> re (D a _) = a
> im (D _ b) = b

> instance Num a => Num (Dual a) where
>   fromInteger i = D (fromInteger i) 0
>   (D a a')+(D b b') = D (a+b) (a'+b')
>   (D a a')-(D b b') = D (a-b) (a'-b')
>   (D a a')*(D b b') = D (a*b) (a*b'+a'*b) 

> instance Floating a => Floating (Dual a) where
>   cos (D a a') = D (cos a) (-sin a*a')
>   sin (D a a') = D (sin a) (cos a*a')
>   exp (D a a') = D (exp a) (exp a*a')

You Could Have Defined Natural Transformations

I'll assume you know what a functor is (and optionally a bit of Haskell):

Consider these two functions

f:Z x Z -> Z    f(x,y) = x+y
g:Z x Z -> Z    f(x,y) = x

There's a fundamental difference between them. In some sense f "makes use" of the fact that we're applying functions to integers, but even though g is acting on integers, it's not using the fact that we're acting on integers. But this is a vague idea. How can we make it precise?

Consider the same in the context of computer programs. This applies in many programming languages, but I'll use Haskell here

> f (x,y) = x+y
> g (x,y) = x

When we ask ghci for the types of these it tells us that f is of type (Num t) => (t, t) -> t but that g is of type (t, t1) -> t. ghci has spotted that f is making use of an interface exposed by integer types (the type class Num) but that g makes no assumptions at all about its arguments. So in a programming language (with static types) we can make the notion above precise by having a type signature that explicitly says something like "I'm not going to make any assumptions about my arguments". But how can we say something like this in the language of set theory, say? Set theory doesn't come with a mechanism for making such promises. (Well, ZF out of the box doesn't, anyway.)

So consider a set-theoretical function like this

h:Z x Z -> Z x Z x Z
h(a,b) = (b,a,a)

What does it mean to say that h doesn't use specific facts about its arguments? It clearly doesn't mean that h doesn't depend on its arguments because the right hand side clearly depends on a and b. We need to say that the right hand side depends on the arguments, but that it somehow doesn't make specific use of them.

Consider an example:

h(1,2) = (2,1,1).

Part of what we mean by not using specific knowledge is that if we were to replace one argument by another value we expect to perform the same substitution on the right. For example, replacing 1 by 3 should give:

h(3,2) = (2,3,3)

So for any 'substitution' function k:Z -> Z we want

h(k(1),k(2)) = (k(2),k(1),k(1))

More generally, we want:

h(k(a),k(b)) = (k(b),k(a),k(a))

Note that we've made use of the functions (a,b) -> (k(a),k(b)) and (a,b,c) -> (k(a),k(b),k(c)) to implement the idea of applying the same substitution to each of the elements. We might want to work with a wide variety of mathematical 'containers' like tuples and sequences as well as more abstract objects. Can we generalise to all of these in one go?

The function that coverts a set X to X x X, or X x X x X forms part of a functor. The other part maps functions between sets to functions between containers. So we can work with the usual ordered pair functor L with LX = X x X and Lf(x,y) = (f(x),f(y)). Similarly we can use MX = X x X x X and Mf(x,y,z) = (f(x),f(y),f(z)). We can stop talking about 'containers' now and talk about functors instead. We can rewrite our rule above as

h(Lk(a,b)) = Mk(h(a,b)).

I've been working with h acting on containers of integers, but if h truly doesn't rely on specific properties of its arguments then there is a version of h for any set X, ie. for any X we have a function

h_X : X x X -> X x X x X
h_X(a,b) = (b,a,a)

and when we use a function k to substitute values in the arguments there's no reason to substitute values from the same set. For example, if we let A be the set of the 26 letters of the alphabet (using your favourite way to encode letters as sets), then we could let k be the function that maps n to the nth letter (modulo 26) and still expect

h_A(Lk(1,2)) = Mk(h_Z(1,2))

to hold. For example

h_A('a','b') = ('b','a','a')

More generally, h is a family of functions such that for any k:X -> Y we have:

h_Y o Lk = Mk o h_X

Note that there's nothing set-specific about this, it's a statement that could hold in any category.

We can take this as a definition. A natural transformation between functors L,M:C -> D is a family of morphisms h_X, one for each object X in C, such that for any morphism k:X -> Y in C we have:

h_Y o Lk = Mk o h_X

And there we have it. Starting with the idea of a function that doesn't 'peek inside' its arguments, we're led inexorably to the idea of substitutability of arguments and from there to a categorical definition of a natural transformation.

As the definition has been generalised from Set to any category you can have fun unpacking the definition again in other categories. For Set-like categories with functions as morphisms you'll often find that natural transformations are still functions that don't make too many assumptions about their arguments, but in each case the rule about what exactly constitutes 'too many assumptions' will be different. The category of vector spaces is a good place to start.

One last point: Category theory was more or less invented for natural transformations. And it's curious how it fits so nicely with the idea of abstraction in computer science, ie. the idea of a function that is limited to using a particular interface on an object.

One interesting result relating to natural transformations is the Yoneda lemma. Note also that the Theorems for Free are also largely about similar kinds of substitutions on function arguments.

Infinitesimal rotations and Lie algebras

A little while back I tried to give a rough idea of what Lie algebras were about. What I want to do now is show how with a bit of Haskell code we can directly get our hands on one. I'm going to assume you know how to make matrices for rotations, and multiply them, and some knowledge from previous posts, for which I'll give links.

Firstly a bit of Haskell administration:

> {-# OPTIONS -fno-warn-missing-methods #-}

Now we need some quick and dirty matrix code:

> data Matrix a = M [[a]] deriving (Eq,Show)

> instance Functor Matrix where
>   fmap f (M a) = M $ map (map f) a

> instance Num a => Num (Matrix a) where
>   M a * M b = M $ mult a b

A Lie Group

What I'm going to do is start by constructing elements of the group of 3D rotations, otherwise known as SO(3), and show how there's another algebraic structure hidden inside it. So let's make some rotation matrices:

> rx theta = M $ [[1,0,0],
>                 [0,cos theta,-sin theta],
>                 [0,sin theta,cos theta]]

> ry theta = M $ [[cos theta,0,sin theta],
>                 [0,1,0],
>                 [-sin theta,0,cos theta]]

> rz theta = M $ [[cos theta,-sin theta,0],
>                 [sin theta,cos theta,0],
>                 [0,0,1]]

These are the three rotations around the x-, y- and z-axes. It's traditional to build arbitrary rotations through the use of Euler angles:

> euler [a,b,c] = rx a*ry b*rz c

The 3D rotations form an example of a Lie group. (A Lie group is essentially a group where the operations like multiplication are differentiable.)

Any 3D rotation can be constructed from a single application of euler. But notice how there's a bit of ugliness in this function. I've made an arbitrary choice of which order to apply the rotations in. I could have defined:

> euler' [a,b,c] = rz c*ry b*rx a

And it's easy to show that euler≠euler'. This is because rotations don't commute. In other words, for rotations, a*b≠b*a.

We can measure the non-commutativity. Remember that any rotation has an inverse. For example rx theta*rx (-theta) gives us the identity matrix because one of these two rotations 'undoes' the other. Given any two rotations we can construct what is known as their commutator:

> commutator a b = inverse a*inverse b*a*b

The idea is that we first perform b, then a, then undo b and then undo a. If a and b commute then this expression can be rearranged to inverse a*inverse b*b*a and then the inverses cancel leaving us with the identity matrix. If they don't commute then we end up with a non-identity matrix. So the comuutator measures the extent to which matrices don't commute.

As I'm feeling lazy, I don't feel like writing inverse. Instead, as I'm only going to work with rotations, I'll use the fact that the inverse of a rotation matrix is the transpose and define:

> inverse = transpose

Try playing with expressions like commutator (rx 1) (ry 2). Note how the numbers quickly get messy. Try to write down closed form expressions for applications of euler and you'll see how complex things can get.

A Lie Algebra

But there's another way to approach the rotation group - through 'infinitesimal' rotations. In my earlier article I just talked about infinitesimal group operations in a hand-wavey way. Now I'm going to make the notion more rigorous. We just need to introduce an infinitesimal number, d, whose square is zero. I've talked about this a lot before so I'm borrowing my earlier code and defining:

> d :: Num a => Dual a
> d = D 0 1

If you try it you'll see that d*d is zero.

Now we can try making infinitesimal rotations:

> rot1 = euler [d,2*d,3*d]

Note how when we evaluate this we get 'nice' numbers. No need to worry about those trig functions any more. And if you look closely at rot1 you'll see that it's essentially the identity matrix plus an infinitesimal part. We can pull the infinitesimal part out using fmap im rot1. You may be able to guess how to build it from the arguments to euler. But first, try evaluating rot2:

> rot2 = euler' [d,2*d,3*d]

It's the same! Try other infiitesimal rotations. When dealing with infinitesimal rotations it doesn't matter what order you apply them in, you get the same result. Working with infinitesimal rotations is looking much easier than working with full=size rotations. In fact, it gets better. Try defining

> rot3 = euler [5*d,-d,2*d]

Now look at fmap im (rot1*rot3) and compare with fmap im rot1 and fmap im rot3. We can multiply infinitesimal rotations simply by adding their infinitesimal parts. In fact, we can define

> star [a,b,c] = M $ [[0, -c, b],
>                     [c,0,-a],
>                     [-b,a,0]]

So we have:

fmap im (euler [a*d,b*d,c*d]) == star [a,b,c]

and

fmap im (euler [a*d,b*d,c*d]*euler [u*d,v*d,w*d]) == fmap im (euler [(a+u)*d,(b+v)*d,(c+w)*d])

Not a single trig expression anywhere!

So we have a simplified way of viewing rotations by looking at infinitesimal rotations. A triple [a,b,c] can be thought of as representing an infinitesimal rotation through star and instead of multiplying matrices we just add the triples elementwise. These triples, together with the binary operation of addition form an example of a Lie algebra. But there's a piece missing. We have lost all information about the non-commutativity of the rotations. It's one thing to simplify, but it's another to lose an important feature of what you're looking at.

The problem is that d is 'too small'. We need an infinitesimal that doesn't go to zero the moment you square it, but is still, in some sense, infinitesimally small. We could rewrite the Dual type. But there's a trick. Define:

> e :: Num a => Dual (Dual a)
> e = D (D 0 1) (D 1 0)

(If you think of Dual as forming a tensor product as I described here then e=1⊗d+d⊗1.)

You can check that e^2 is non-zero but e^3 is zero.

Now when we compute a commutator we get something a little different:

> comm1 = commutator (euler [e,0,0]) (euler [0,e,0])

fmap re comm1 is essentially the identity as before. But if we look at fmap im comm2 there's a non-zero infinitesimal piece which is different from what we had when we worked with d. This infinitesimal piece is in fact proportional to e^2. As (im . im) (e^2) is a half, we can extract the coefficient of e^2 from comm1 using

> im2 x = im (im x)/2

In fact, we have fmap im2 comm1 == star [0,0,1]. So by choosing infinitesimals that aren't too small, we haven't lost information about non-commutativity. In fact, a bit of experimentation may convince you that with:

> shrink u = fmap (e*) u

we have:

fmap im2 (commutator (euler (shrink u)) (euler (shrink v))) == star (u `cross` v)

So let's step back a bit and recapitulate all this in something approaching English:

Given a tiny rotation we can represent it as three Euler angles a, b, c, all of which are tiny. We can think of a, b and c as forming a vector [a,b,c]. When we do this, apart from an even smaller error, multiplication of rotations becomes ordinary addition of vectors and the order of rotations isn't significant. But if we choose not to ignore this small error we see that a rotation represented by u and a rotation represented by v don't quite commute and the order does matter. The size of this error is measured by the cross product of u and v. This is intuitively plausible, we'd expect that rotations defined by vectors in a similar direction would be closer to commuting, and this is reflected in the fact that the cross product is zero for parallel vectors.

So now I can fill in the missing piece from the description of the Lie algebra I gave above. The Lie algebra so(3) consists of the 3d vectors (representing infinitesimal rotations), addition of vectors (representing multiplcation of infinitesimal rotations) and the cross product (measuring the amount by which small rotations fail to commute). This picture doesn't just apply to the rotations, a similar one applies for any Lie group. We get the same pattern of a simplified form of multiplication and a binary operation that measures non-commutativity.

But you might still ask if there's something missing. What if we defined f so that f^4 is zero but f^3 isn't? Would we extract more information about the non-commutativity of rotations? Well the interesting fact, which I won't prove here, is that it doesn't. Almost everything you need to know about a Lie group can be extracted from its Lie algebra. In fact, from a group's Lie algebra you can almost recover the original group. (In fact, what you recover is its universal cover.) But notice how there are no trig formulae involved when talking about Lie algebras. Lie algebras give a really nice way to study Lie groups without getting your hands too dirty. But this isn't the only reason to study Lie algebras, many physical properties arising from symmetries are more naturally studied through the Lie algebra because Lie algebras arise from Lie groups as soon as you start differentiating things. In particular, Lie algebras play a major role in Noether's theorem, one of the cornerstones of modern theoretical physics.

In summary then, I hope I've given some flavour of Lie algebras. Using infinitesimals you can get your hands on them directly without the use of a large amount of mathematical machinery. There has been some difficult stuff here, but I'm hoping that the freedom you now have at the Haskell prompt to play with the things I've been talking about will make up for the inadequacies of my explanations.

One last word: in principle we could do the same with E8. But we'd need 248 by 248 matrices, and the equivalent of euler would need 248 parameters. (That's a coincidence, in general the number of parameters needed to define an element of a Lie group isn't equal to the dimension of the matrix, but it is for SO(3) and E8).

---

Appendix (the bits of code I left out above)

Defining the infinitesimals:

> data Dual a = D a a deriving (Show,Eq)

Extract the 'full-size' and 'infinitesimal' parts of a number:

> re (D a _) = a
> im (D _ b) = b

> instance Num a => Num (Dual a) where
>   fromInteger i = D (fromInteger i) 0
>   (D a a')+(D b b') = D (a+b) (a'+b')
>   (D a a')-(D b b') = D (a-b) (a'-b')
>   (D a a')*(D b b') = D (a*b) (a*b'+a'*b) 

> instance Floating a => Floating (Dual a) where
>   cos (D a a') = D (cos a) (-sin a*a')
>   sin (D a a') = D (sin a) (cos a*a')

> instance Fractional a => Fractional (Dual a)

Some useful matrix and vector operations:

> mult a ([]:_) = map (const []) a
> mult a b = zipWith (:) (map (dot (map head b)) a) (mult a (map tail b))

> transpose (M a) = M $ transpose' a where
>   transpose' [] = repeat [] 
>   transpose' (xs : xss) = zipWith (:) xs (transpose' xss) 

Some simple vector operations

> dot a b = foldr (+) 0 $ zipWith (*) a b
> cross [a,b,c] [d,e,f] = [b*f-c*e,c*d-a*f,a*e-b*d]

Negative Probabilities

I'm always interested in new ways to present the ideas of quantum mechanics to a wider audience. Unfortunately the conventional description of quantum mechanics requires knowledge of complex numbers and vector spaces making it difficult to avoid falling back on misleading analogies and metaphors. But there's another approach to quantum mechanics, due to Feynman, that I have never seen mentioned in the popular science literature. It assumes only basic knowledge of probability theory and negative numbers to get a foot on the ladder.

We start with tweaking probability theory a bit. One of the axioms of probability theory says that all probabilities must lie in the range zero to one. However, we could imagine relaxing this rule even though on the face of it it seems meaningless. For example, suppose we have a coin that has a 3/2 chance of landing heads and a -1/2 chance of landing tails. We can still reason that the chance of getting two heads in a row is 3/2×3/2=9/4 by the usual multiplication rule. But obviously no situation like this could ever arise in the real world because after tossing such a coin 10 times we'd expect to see -5 tails on average.

But what if we could contrive a system with some kind of internal state governed by negative probabilities even though we couldn't observe it directly? So consider this case: a machine produces boxes with (ordererd) pairs of bits in them, each bit viewable through its own door. Let's suppose the probability of each possible combination of two bits is given by the following table:

First bit	Second bit	Probability
0	0	-1/2
0	1	1/2
1	0	1/2
1	1	1/2

Obviously if we were able to look through both doors we'd end up with the meaningless prediction that we'd expect to see 00 a negative number of times. But suppose that things are arranged so that we can only look through one door. Maybe the boxes self-destruct if one or other door is opened but you still get enough time to see what was behind the door. Now what happens?

If you look through the first door the probability of seeing 1 is P(10)+P(11)=1. We get the same result if we look through the second door. We only get probabilities in the range zero to one. As long as we're restricted to one door we get meaningful results.

If we were to perform this experiment repeatedly with different runs of the machine, each time picking a random door to look through, we'd eventually become very confident that every box contained 11. After all, if we freely choose which door to look through, and we always see 1, there's no place 0's could be 'hiding'.

But now suppose a new feature is added to the box that allows us to compare the two bits to see if they are equal. It reveals nothing about what the bits are, just their state of equality. And of course, after telling us, it self-destructs. We now find that the probability of the two bits being different is P(01)+P(10)=1. So if we randomly chose one of the three possible observations each time the machine produced the box we'd quickly run into the strange situation that the two bits both appear to be 1, and yet are different. But note that although the situation is weird, it's not meaningless. As long as we never get to see both bits at the same time we never directly observe a paradox. If we met such boxes in the real world we'd be forced to conclude that maybe the boxes knew which bit you were going to look at and changed value as a result, or that maybe you didn't have the free will to choose door that you thought you had, or maybe, even more radically, you'd conclude that the bits generated by the machine were described by negative probabilities.

That's all very well, but obviously the world doesn't really work like this and we never see boxes like this. Except that actually it does! The EPR experiment has many similarities to the scenario I described above. The numbers aren't quite the same, and we're not talking about bits in boxes, but we do end up with a scenario involving observations of bits that simply don't add up. If we do try to explain what's going on using probability theory, we either conclude there's something weird about our assumptions of locality or causality or we end up assigning negative probabilities to the internal states of our systems. In fact, you can read the details in an article by David Schneider. Being forced to conclude that we have negative probabilities in a physical system is usually taken as a sign that we have a contradiction. In the case of Bell's theorem it shows that we can't interpret what we see in terms of probability theory and hence that the weirdness of quantum mechanics can't be explained in terms of some hidden random variable that we can't see. QM simply doesn't obey the rules you'd expect of hidden variables we can't see.

But in a paper called Negative Probability, Feynman tried taking the idea of negative probabilities seriously. He showed that you could reformulate quantum mechanics completely in terms of them so that you no longer needed to think in terms of the complex number valued 'amplitudes' that physicists normally use. This means the above isn't just an analogy, it's actually a formal system within which you can do QM, although I haven't touched on the bit that refers to the dynamics of quantum systems. So if you can get your head around the ideas I've talked about above you're well on your way to understanding some reasons why quantum mechanics seems so paradoxical.

At this point you may be wondering how nature contrives to hide these negative probabilities from direct observation. Her trick is that making one kind of observation disturbs up the state of what you've observed so that you can't make the other kind of observation on a pristine state. You have to pick one kind of observation or the other. Electrons and photons really are a lot like the boxes I just described.

So why don't physicists use this formulation? Despite the fact that negative numbers seem simpler to most people than imaginary numbers, the negative number formulation of QM is much more complicated. What's more, because it makes exactly the same predictions as regular QM there's no compelling reason to switch to it. And anyway, it's not as if directly observing negative probabilities is any more intuitive or meaningful than imaginary ones. Once you've introduced negative ones, you might as well go all the way!

This all ties in with what I said a while back. The important thing about QM is that having two ways to do something can make it less likely to happen, not more.

For a different perspective this is an interesting comment.

Footnote: We can embed QM in negative probability theory. But can we do the converse? Can every negative probability distribution be physically realised in a quantum system? I've a hunch the answer is obvious but I'm too stupid to see it.

---