Haskell for Maths

New release of HaskellForMaths

2011-11-12T21:24:00.001Z

I've just uploaded a new version v0.4.1 of HaskellForMaths, containing three new modules and a couple of other improvements. The additions are as follows:

Math.Algebras.Quaternions

This module was already present: it defines the quaternion algebra on the basis {1,i,j,k}, where multiplication is defined by:
i^2 = j^2 = k^2 = ijk = -1

This is enough information to figure out the full multiplication table. For example:
ijk = -1
=> (ijk)k = -k
=> ij(kk) = -k (associativity of multiplication)
=> ij = k
It turns out that the basis elements i,j,k anti-commute in pairs, eg ij = -ji, etc.

In this release I've added a couple of new things.

First, the quaternions are a division algebra, so I've added a Fractional instance.

Specifically, we can define a conjugation operation on the quaternions (similar to complex conjugation) via
conj (w+xi+yj+zk) = w-xi-yj-zk
Then we can define a quadratic norm via
sqnorm q = q * conj q = w^2+x^2+y^2+z^2
Since the norm is always a scalar, we can define a multiplicative inverse by
q^-1 = conj q / sqnorm q

For example:

$ cabal install HaskellForMaths
$ ghci
> :m Math.Algebras.Quaternions
> (2*i+3*j)^-1 :: Quaternion Q
-2/13i-3/13j

(If you leave out the type annotation, you'll be working in Quaternion Double.)

Second, the quaternions have an interesting role in 3- and 4-dimensional geometry.

Given any non-zero quaternion q, the map x -> q^-1 x q turns out to be a rotation of the 3-dimensional space spanned by {i,j,k}. To multiply rotations together (ie do one then another), just multiply the quaternions. This turns out to be a better way to represent rotations than 3*3 matrices:
- It's more compact - four scalars rather than nine
- They're faster to multiply - 16 scalar multiplications versus 27
- It's more robust against rounding error - whatever quaternion you end up with will still represent a rotation, whereas a sequence of matrix multiplications of rotations might not be quite a rotation any more, due to rounding error.

If you're curious, the function reprSO3 converts a quaternion to the corresponding 3*3 matrix:

> reprSO3 (1+2*i) :: [[Q]]
[[1,0,0],[0,-3/5,-4/5],[0,4/5,-3/5]]

(Exercise: Figure out why we got this matrix.)

Quaternions can also be used to represent rotations of 4-dimensional space - see the documentation.

Math.Algebras.Octonions

This is a new module, providing an implementation of the 8-dimensional non-associative division algebra of octonions. I follow Conway's notation [1], so the octonions have basis {1,e0,e1,e2,e3,e4,e5,e6}, with multiplication defined by:
e_i * e_i = -1, for i in [0..6]
e_i+1 * e_i+2 = e_i+4, where the indices are taken modulo 7.

The octonions are not associative, but they are an inverse loop, so they satisfy x^-1(xy) = y = (yx)x^-1. This is enough to enable us to deduce the full multiplication table from the relations above.

Like the quaternions, the octonions have conjugation and a norm, and multiplicative inverses:

> :l Math.Algebras.Octonions
> (2+i0+2*i3)^-1
2/9-1/9i0-2/9i3

The octonions are an exceptional object in mathematics: there's nothing else quite like them. They can be used to construct various other exceptional objects, such as the root lattice E8, or the Lie group G2. Hopefully I'll be able to cover some of that stuff in a future installment.

[1] Conway and Smith, On Quaternions and Octonions

Math.NumberTheory.Prime

The main function in this module is isPrime :: Integer -> Bool, which tells you whether a number is prime or not. It's implemented using the Miller-Rabin test.

The basic idea of the test is:
- If p is prime, then Zp is a field
- In a field, the equation x^2 = 1 has only two solutions, 1 and -1
- Given an arbitrary b coprime to p, we know from Fermat's little theorem that b^(p-1) = 1 (mod p)
- So if p-1 = q * 2^s, with q odd, then either b^q = 1 (mod p), or there is some r, 0 <= r < s with b^(q*2^r) = -1 (mod p)

The idea of the algorithm is to try to show that p isn't prime by trying to find a b where the above is not true. We take several different values of b at random, and repeatedly square b^q, to see whether we get -1 or not.

The advantage of the Miller-Rabin test, as compared to trial division say, is that it has a fast running time even for very large numbers. For example:

> :m Math.NumberTheory.Prime
> :set +s
> nextPrime $ 10^100
10000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000267
(0.09 secs, 14904632 bytes)

The potential disadvantage of the Miller-Rabin test is that it is probabilistic: There is a very small chance (1 in 10^15 in this implementation) that it could just fail to hit on a b which disproves n's primeness, so that it would say n is prime when it isn't. In practice, at those odds it's not worth worrying about.

Math.NumberTheory.Factor

The main function in this module is pfactors :: Integer -> [Integer], which returns the prime factors of a number (with multiplicity). It uses trial division to try to find prime factors less than 10000. After that, it uses the elliptic curve method to try to split what remains. The elliptic curve method relies on some quite advanced maths, but the basic idea is this:
- If p is a prime, then Zp is a field
- Given a field, we can do "arithmetic on elliptic curves" over the field.
- So to factor n, pretend that n is prime, try doing arithmetic on elliptic curves, and wait till something goes wrong.
- It turns out that if we look at what went wrong, we can figure out a non-trivial factor of n

Here it is in action:

> :m Math.NumberTheory.Factor
> pfactors $ 10^30+5
[3,5,4723,1399606163,10085210079364883]
(0.55 secs, 210033624 bytes)
> pfactors $ 10^30+6
[2,7,3919,758405810021,24032284101871]
(2.31 secs, 883504748 bytes)

I love the way it can crunch through 12-digit prime factors with relative ease.

Commutative Algebra and Algebraic Geometry

2011-09-18T19:22:00.000+01:00

Last time we saw how to create variables for use in polynomial arithmetic. This time I want to look at some of the things we can do next.

First, let's define the variables we are going to use:

> :l Math.CommutativeAlgebra.GroebnerBasis
> let [t,u,v,x,y,z,x',y',z'] = map glexvar ["t","u","v","x","y","z","x'","y'","z'"]

So now we can do arithmetic in the polynomial ring Q[t,u,v,x,y,z,x',y',z']. For example:

> (x+y)^2
x^2+2xy+y^2

The branch of mathematics dealing with the theory of polynomial rings is called commutative algebra, and it was "invented" mainly in support of algebraic geometry. Algebraic geometry is roughly the study of the curves, surfaces, etc that arise as the solution sets of polynomial equations. For example, the solution-set of the equation x^2+y^2=1 is the unit circle.

If we are given any polynomial equation f = g, then we can rewrite it more conveniently as f-g = 0. In other words, we only need to track individual polynomials, rather than pairs of polynomials. Call the solution set of f = 0 the zero-set of f.

Sometimes we're interested in the intersection of two or more curves, surfaces, etc. For example, it is well known that the hyperbola, parabola and ellipse all arise as "conic sections" - that is, as the intersection of a cone with a plane. So define the zero-set of a collection (or system) of polynomials to be the set of points which are zeros of all the polynomials simultaneously. For example, the zero-set of the system [x^2+y^2-z^2, z-1] is the unit circle x^2+y^2=1 situated on the plane z=1.

Okay, so how can commutative algebra help us to investigate curves and surfaces? Is there a way for us to "do geometry by doing algebra"? Well, first, what does "doing geometry" consist of? Well, at least some of the following:
- Looking at the shapes of curves and surfaces
- Looking at intersections, unions, differences and products of curves and surfaces
- Looking at when curves or surfaces can be mapped into or onto other curves or surfaces
- Looking at when two different curves or surfaces are equivalent, in some sense (for example, topologically equivalent)

(That phrase "curves and surfaces" is not only clumsy but also inaccurate, so from now on I'll use the proper term, "variety", for the zero-set of a system of polynomials, whether it's a set of isolated points, a curve, a surface, some higher dimensional thing, or a combination of some of the preceding.)

So can we do all those things using algebra? Well, let's have a go.

Let's start by looking at intersections and unions of varieties (remember, that's just the fancy name for curves, surfaces, etc.).

Well, we've already seen how to do intersections. If a variety V1 is defined by a system of polynomials [f1...fm], and a variety V2 is defined by [g1...gn], then their intersection is defined by the system [f1...fm,g1...gn] - the zero-set of both sets of polynomials simultaneously. We'll call this the "sum" of the systems of polynomials. (Note to the cognoscenti: yes, I'm really talking about ideals here.)

sumI fs gs = gb (fs ++ gs)

Don't worry too much about what that "gb" (Groebner basis) call is doing. Let's just say that it's choosing the best way to represent the system of polynomials. For example:

> sumI [x^2+y^2-z^2] [z-1]
[x^2+y^2-1,z-1]

Notice how the gb call has caused the first polynomial to be simplified slightly. The same variety might arise as the zero-set of many different systems of polynomials. That's something that we're going to need to look into - but later.

Okay, so what about unions of varieties. So suppose V1 is defined by [f1...fm], V2 is defined by [g1...gn]. A point in their union is in either V1 or V2, so it is in the zero-set of either [f1...fm] or [g1...gn]. So how about multiplying the polynomials together in pairs. That is, let's look at the system [fi*gj | fi <- fs, gj <- gs]. Call the zero-set of this system V. Then clearly, any point in either V1 or V2 is in V, since we then know that either all the fs or all the gs vanish at that point, and hence so do all the products. Conversely, suppose that some point is not in the union of V1 and V2. Then there must exist some fi, and some gj, which are non-zero at that point. Hence there is an fi*gj which is non-zero, so the point is not in V.

This construction is called, naturally enough, the product of the systems of polynomials.

productI fs gs = gb [f * g | f <- fs, g <- gs]

> productI [x^2+y^2-z^2] [z-1]
[x^2z+y^2z-z^3-x^2-y^2+z^2]

Just in case you're still a little unsure, let's confirm that a few arbitrary points in the union are in the zero-set of this polynomial:

> eval (x^2*z+y^2*z-z^3-x^2-y^2+z^2) [(x,100),(y,-100),(z,1)]
0
> eval (x^2*z+y^2*z-z^3-x^2-y^2+z^2) [(x,3),(y,4),(z,5)]
0

The first expression evaluates the polynomial at the point (100,-100,1), an arbitrary point on the plane z=1. The second evaluates at (3,4,5), an arbitrary point on the cone x^2+y^2=z^2. Both points are in the zero-set of our product polynomial.

Since we're in the neighbourhood, let's have a look at the other conic sections. First, let's rotate our coordinate system by 45 degrees, using the substitution x'=x+z, z'=z-x. (Okay, so this also scales - to save us having to handle a sqrt 2 factor.)

> let cone' = subst (x^2+y^2-z^2) [(x,(x'-z')/2),(y,y'),(z,(x'+z')/2)]
> cone'
-x'z'+y'^2

In these coordinates, the intersection of the cone with the plane z'=1 is the parabola x'=y'^2:

> sumI [cone'] [z'-1]
[y'^2-x',z'-1]

Alternatively, the intersection with the plane y'=1 is the hyperbola x'z'=1:

> sumI [cone'] [y'-1]
[x'z'-1,y'-1]

Okay, so we've made a start on seeing how to do geometry by doing algebra, by looking at unions and intersections of varieties. There's still plenty more to do. We mustn't forget that we have some unfinished business: we need to understand when different polynomial systems can define the same variety, and in what sense the gb (Groebner basis) function finds the "best" representation. That will have to wait for another time.

Incidentally, for the eval and subst functions that I used above, you will need to take the new release HaskellForMaths v0.4.0. In this release I also removed the older commutative algebra modules, so I revved the minor version number.

Commutative Algebra in Haskell, part 1

2011-09-04T21:11:00.000+01:00

Once again, it's been a little while since my last post, and once again, my excuse is partly that I've been too busy writing code.

I've just uploaded a new release, HaskellForMaths 0.3.4, which contains the following new modules:

Math.Core.Utils - this is a collection of utility functions used throughout the rest of the library. I've belatedly decided that it's better to put them all in one place rather than scattered here and there throughout other modules.

Math.Core.Field - this provides new, more efficient implementations of several finite fields. There already were implementations of these finite fields, in the Math.Algebra.Field.Base and ...Extension modules, as discussed here and here. However, that code was written to make the maths clear, rather than for speed. This new module is about speed. For the prime power fields in particular (eg F4, F8, F9), these implementations are significantly faster.

Math.Combinatorics.Matroid - Matroids are a kind of combinatorial abstraction of the concept of linear independence. They're something that I heard about years ago - both of my favourite combinatorics books have brief introductions - but I never bothered to follow up. Well anyway, so something finally piqued my curiosity, and I got Oxley's Matroid Theory. It turned out to be really interesting stuff, and this module is pretty much a translation of a large part of that book into Haskell code, written as I taught myself all about matroids.

Math.CommutativeAlgebra.Polynomial - Although I hadn't yet got around to discussing them in the blog, HaskellForMaths has always had modules for working with multivariate polynomials, namely Math.Algebra.Commutative.Monomial and ...MPoly. However, these were some of the earliest code I wrote, before my more recent free vector space and algebra code. So I saw an opportunity to simplify and improve this code, by building it on top of the free vector space code. Also, I'm trying to rationalise the module naming convention in HaskellForMaths, to more closely follow the categories used in arxiv.org or mathoverflow.net . In the long run, I expect this module to supercede the older modules.

Math.CommutativeAlgebra.GroebnerBasis - Again, there was already code for Groebner bases in Math.Algebra.Commutative.GBasis. This is pretty much the same code, ported to the new polynomial implementation, but I've also begun to build on this, with code to find the sum, product, intersection, and quotient of ideals.

So the matroid code was just new code that I wrote while teaching myself some new maths. But most of the other code comes from an ambition to organise and simplify the HaskellForMaths library. I've also been trying to improve the documentation.

My ultimate ambition is to get more people using the library. To do that, the structure of the library needs to be clearer, the documentation needs to be better, and I need to explain how to use it. So I thought I'd start by explaining how to use the new commutative algebra modules.

(So this is a bit of a digression from the series on quantum algebra that I've been doing the last few months. However, in terms of the cumulative nature of maths, it's probably better to do this first.)

Okay, so suppose we want to do some polynomial arithmetic. Well, first we need to create some variables to work with. How do we do that?

First, decide on a monomial ordering - that is, we need to decide in what order monomials are to be listed within a polynomial. For the moment, let's use "graded lexicographic" or Glex order. This says that you should put monomials of higher degree before those of lower degree (eg y^3 before x^2), and if two monomials have the same degree, you should use lexicographic (dictionary) order (eg xyz before y^3).

Next, decide on a field to work over. Most often, we'll want to work over Q, the rationals.

Then, our variables themselves can be of any Haskell type - but there are usually only two sensible choices:

The easiest way is to use String as the type for our variables.

Then we could make some variables like this:

> :l Math.CommutativeAlgebra.Polynomial
> let [x,y,z] = map glexvar ["x","y","z"]

And then we can do polynomial arithmetic:

> (x+y+z)^3
x^3+3x^2y+3x^2z+3xy^2+6xyz+3xz^2+y^3+3y^2z+3yz^2+z^3

If we want to use any other field besides Q, then we will have to use a type annotation to tell the compiler which field we're working over:

> let [x,y,z] = map var ["x","y","z"] :: [GlexPoly F3 String]
> (x+y+z)^3
x^3+y^3+z^3

The alternative to using String for our variables is to define our own type. For example



data Var = X | Y | Z | W deriving (Eq,Ord)



instance Show Var where

    show X = "x"

    show Y = "y"

    show Z = "z"

    show W = "w"



[x,y,z,w] = map glexvar [X,Y,Z,W]

So there you have it - now you can do polynomial arithmetic in Haskell.

So how does it work?

Well, fundamentally, k-polynomials are a free k-vector space on the basis of monomials. So we define a type to implement monomials:



data MonImpl v = M Int [(v,Int)] deriving (Eq)

-- The initial Int is the degree of the monomial. Storing it speeds up equality tests and comparisons



instance Show v => Show (MonImpl v) where

    show (M _ []) = "1"

    show (M _ xis) = concatMap (\(x,i) -> if i==1 then showVar x else showVar x ++ "^" ++ show i) xis

        where showVar x = filter ( /= '"' ) (show x) -- in case v == String

Notice that our monomial implementation is polymorphic in v, the type of the variables.

Next, monomials form a monoid, so we make them an instance of Mon (the HaskellForMaths class for monoids):



instance (Ord v) => Mon (MonImpl v) where

    munit = M 0 []

    mmult (M si xis) (M sj yjs) = M (si+sj) $ addmerge xis yjs

In principle, all we need to do now is define an Ord instance, and then an Algebra instance, using the monoid algebra construction.

However, for reasons that will become clear in future postings, we want to be able to work with various different orderings on monomials, such as Lex, Glex, or Grevlex. So we provide various newtype wrappers around this basic monomial implementation. Here's the code for the Glex ordering that we used above:



newtype Glex v = Glex (MonImpl v) deriving (Eq, Mon) -- GeneralizedNewtypeDeriving



instance Show v => Show (Glex v) where

    show (Glex m) = show m



instance Ord v => Ord (Glex v) where

    compare (Glex (M si xis)) (Glex (M sj yjs)) =

        compare (-si, [(x,-i) | (x,i) <- xis]) (-sj, [(y,-j) | (y,j) <- yjs])



type GlexPoly k v = Vect k (Glex v)



glexvar :: v -> GlexPoly Q v

glexvar v = return $ Glex $ M 1 [(v,1)]



instance (Num k, Ord v, Show v) => Algebra k (Glex v) where

    unit x = x *> return munit

    mult xy = nf $ fmap (\(a,b) -> a `mmult` b) xy

We also have similar newtypes for Lex and Grevlex orderings, which I'll discuss another time.

And that's pretty much it. Now that we have an instance of Algebra k (Glex v), we get a Num instance for free, so we get +, -, *, and fromInteger. That means we can enter expressions like the following:

> (2*x^2-y*z)^2
4x^4-4x^2yz+y^2z^2

Note that division is not supported: you can't write x/y, for example. However, as a convenience, I have defined a partial instance of Fractional, which does let you divide by scalars. That means that it's okay to write x/2, for example.

Next time, some more things you can do with commutative algebra.

The Tensor Algebra Monad

2011-07-10T20:09:00.001+01:00

It's been a little while since my last post. That's partly because I've been busy writing new code. I've put up a new release, HaskellForMaths 0.3.3, which contains three new modules:
- Math.Combinatorics.Digraph
- Math.Combinatorics.Poset
- Math.Combinatorics.IncidenceAlgebra

I'll go through their contents at some point, but this time I want to talk about the tensor algebra.

So recall that previously we defined the free vector space over a type, tensor products, algebras and coalgebras in Haskell code.

In HaskellForMaths, we always work with the free vector space over a type: that means, we take some type b as a basis, and form k-linear combinations of elements of b. This construction is represented by the type Vect k b.

Given two vector spaces A = Vect k a, B = Vect k b, we can form their tensor product A⊗B = Vect k (Tensor a b). So Tensor is a type constructor on basis types, which takes basis types a, b for vector spaces A, B, and returns a basis type for the tensor product A⊗B.

We also defined a type constructor DSum, which returns a basis type for the direct sum A⊕B.

Now, we saw that tensor product is a monoid (at the type level, up to isomorphism):
- it is associative: (A⊗B)⊗C is isomorphic to A⊗(B⊗C)
- it has a unit: the field k itself is an identity for tensor product, in the sense that k⊗A is isomorphic to A, is isomorphic to A⊗k

Given some specific vector space V, we can consider the tensor powers of V:
k, V, V⊗V, V⊗V⊗V, ...
(We can omit brackets in V⊗V⊗V because tensor product is associative.)

And indeed we can form their direct sum:
T(V) = k ⊕ V ⊕ V⊗V ⊕ V⊗V⊗V ⊕ ...
(where an element of T(V) is understood to be a finite sum of elements of the tensor powers.)

This is a vector space, since tensor products and direct sums are vector spaces. If V has a basis e1,e2,e3,..., then a typical element of T(V) might be something like 3 + 5e2 + 2e1⊗e3⊗e1.

Now the interesting thing is that T(V) can be given the structure of an algebra, as follows:
- for the unit, we use the injection of k into the first direct summand
- for the mult, we use tensor product

For example, we would have
e2 * (2 + 3e1 + e4⊗e2) = 2e2 + 3e2⊗e1 + e2⊗e4⊗e2

With this algebra structure, T(V) is called the tensor algebra.

So how should we represent the tensor algebra in HaskellForMaths? Suppose that V is the free vector space Vect k a over some basis type a. (Recall also that the field k itself can be represented as the free vector space Vect k () over the unit type.) Can we use the DSum and Tensor type constructors to build the tensor algebra? Something like:
Vect k (DSum () (DSum a (DSum (Tensor a a) (DSum ...))))

Hmm, that's not going to work - we can't build the whole of what we want that way. (Unless some type system wizard knows otherwise?) So instead of representing the direct sum and tensor product at the type level, we're going to have to do it at the value level. Here's the definition:

data TensorAlgebra a = TA Int [a] deriving (Eq,Ord)

Given the free vector space V = Vect k a over basis type a, then TensorAlgebra a is the basis type for the tensor algebra over a, so that Vect k (TensorAlgebra a) is the tensor algebra T(V). The Int in TA Int [a] tells us which direct summand we're in (ie which tensor power), and the [a] tells us the tensor multiplicands. So for example, e2⊗e1⊗e4 would be represented as TA 3 [e2,e1,e4]. Then Vect k (TensorAlgebra a) consists of k-linear combinations of these basis elements, so it is the vector space T(V) that we are after.

Here's a Show instance:

instance Show a => Show (TensorAlgebra a) where
    show (TA _ []) = "1"
    show (TA _ xs) = filter (/= '"') $ concat $ L.intersperse "*" $ map show xs

It will be helpful to have a vector space basis to work with, so here's one that we used previously:

newtype EBasis = E Int deriving (Eq,Ord)

instance Show EBasis where show (E i) = "e" ++ show i

Then, for example, our Show instance gives us:

> :l Math.Algebras.TensorAlgebra
> return (TA 0 []) <+> return (TA 2 [E 1, E 3])
1+e1*e3

(Recall that the free vector space is a monad, hence our use of return to put a basis element into the vector space.)

So note that in the show output, the "*" is meant to represent tensor product, so this is really 1+e1⊗e3. You can actually get Haskell to output the tensor product symbol - just replace "*" by "\x2297" in the definition of show - however I found that it didn't look too good in the Mac OS X terminal, and I wasn't sure it would work on all OSes.

Ok, how about an Algebra instance? Well, TensorAlgebra a is basically just a slightly frilly version of [a], so it's a monoid, and we can use the monoid algebra construction:

instance Mon (TensorAlgebra a) where
    munit = TA 0 []
    mmult (TA i xs) (TA j ys) = TA (i+j) (xs++ys)

instance (Num k, Ord a) => Algebra k (TensorAlgebra a) where
    unit x = x *> return munit
    mult = nf . fmap (\(a,b) -> a `mmult` b)

So now we can do arithmetic in the tensor algebra:

> let e_ i = return (TA 1 [E i]) :: Vect Q (TensorAlgebra EBasis)
> let e1 = e_ 1; e2 = e_ 2; e3 = e_ 3; e4 = e_ 4
> (e1+e2) * (1+e3*e4)
e1+e2+e1*e3*e4+e2*e3*e4

We've got into the habit of using QuickCheck to check algebraic properties. Let's just check that the tensor algebra, as we've defined it, is an algebra:

instance Arbitrary b => Arbitrary (TensorAlgebra b) where
    arbitrary = do xs <- listOf arbitrary :: Gen [b] -- ScopedTypeVariables
                   let d = length xs
                   return (TA d xs)

prop_Algebra_TensorAlgebra (k,x,y,z) = prop_Algebra (k,x,y,z)
    where types = (k,x,y,z) :: ( Q, Vect Q (TensorAlgebra EBasis), Vect Q (TensorAlgebra EBasis), Vect Q (TensorAlgebra EBasis) )

> quickCheck prop_Algebra_TensorAlgebra
+++ OK, passed 100 tests.

Ok, so what's so special about the tensor algebra? Well, it has a rather nice universal property:
Suppose A = Vect k a, B = Vect k b are vector spaces, and we have a linear map f : A -> B. Suppose that B is also an algebra. Then we can "lift" f to an algebra morphism f' : T(A) -> B, such that the following diagram commutes.

In other words, f' agrees with f on the copy of A within T(A): f = f' . i

Ah, but hold on, I didn't say what an algebra morphism is. Well, it's just the usual: a function which "commutes" with the algebra structure. Specifically, it's a linear map (so that it commutes with the vector space structure), which makes the following diagrams commute:

So how does this universal property work then? Here's the code:

injectTA :: Num k => Vect k a -> Vect k (TensorAlgebra a)
injectTA = fmap (\a -> TA 1 [a])

liftTA :: (Num k, Ord b, Show b, Algebra k b) =>
     (Vect k a -> Vect k b) -> Vect k (TensorAlgebra a) -> Vect k b
liftTA f = linear (\(TA _ xs) -> product [f (return x) | x <- xs])

In other words, any tensor product u⊗v⊗... is sent to f(u)*f(v)*...

Let's look at an example. Recall that the quaternion algebra H has the basis {1,i,j,k}, with i^2 = j^2 = k^2 = ijk = -1.

> let f = linear (\(E n) -> case n of 1 -> 1+i; 2 -> 1-i; 3 -> j+k; 4 -> j-k; _ -> zerov)
> let f' = liftTA f
> e1*e2
e1*e2
> f' (e1*e2)
2

Recall that we usually define a linear map by linear extension from its action on a basis - that's what the "linear" is doing in the definition of f. It's fairly clear what f' is doing: it's basically just variable substitution. That is, we can consider the basis elements ei as variables, and the tensor algebra as the algebra of non-commutative polynomials in the ei. Then the linear map f assigns a substitution to each basis element, and f' just substitutes and multiplies out in the target algebra. In this case, we have:
e1⊗e2 -> (1+i)*(1-i) = 1-i^2 = 2

We can use QuickCheck to verify that liftTA f is indeed the algebra morphism required by the universal property. Here's a QuickCheck property for an algebra morphism. (We don't bother to check that f is a linear map, since it's almost always clear from the definition. If in doubt, we can test that separately.)

prop_AlgebraMorphism f (k,x,y) =
    (f . unit) k == unit k &&
    (f . mult) (x `te` y) == (mult . (f `tf` f)) (x `te` y)

This is just a transcription of the diagrams into code.

In order to test the universal property, we have to check that liftTA f is an algebra morphism, and that it agrees with f on (the copy of) V (in T(V)):

prop_TensorAlgebra_UniversalProperty (fmatrix,(k,x,y),z) =
    prop_AlgebraMorphism f' (k,x,y) &&
    (f' . injectTA) z == f z
    where f = linfun fmatrix
          f' = liftTA f
          types = (fmatrix,(k,x,y),z) :: (LinFun Q EBasis HBasis,
                                         (Q,Vect Q (TensorAlgebra EBasis), Vect Q (TensorAlgebra EBasis)),
                                         Vect Q EBasis)

So the key to this code is the parameter fmatrix, which is an arbitrary (sparse) matrix from Q^n to H, the quaternions, from which we build an arbitrary linear function. Note that the universal property of course implies that we can choose any algebra as the target for f - I just chose the quaternions because they're familiar.

> quickCheck prop_TensorAlgebra_UniversalProperty
+++ OK, passed 100 tests.

With this construction, tensor algebra is in fact a functor from k-Vect to k-Alg. The action on objects is V -> T(V), Vect k a -> Vect k (TensorAlgebra a). But a functor also acts on the arrows of the source category.

How do we get an action on arrows? Well, we can use the universal property to construct one. If we have an arrow f: A -> B, then (injectTA . f) is an arrow A -> T(B). Then we use the universal property to lift to an arrow f': T(A) -> T(B).

Here's the code:

fmapTA :: (Num k, Ord b, Show b) =>
    (Vect k a -> Vect k b) -> Vect k (TensorAlgebra a) -> Vect k (TensorAlgebra b)
fmapTA f = liftTA (injectTA . f)

For example:

newtype ABasis = A Int deriving (Eq,Ord,Show)
newtype BBasis = B Int deriving (Eq,Ord,Show)

> let f = linear (\(A i) -> case i of 1 -> return (B 1) <+> return (B 2);
                                      2 -> return (B 3) <+> return (B 4);
                                      _ -> zerov :: Vect Q BBasis)
> let f' = fmapTA f
> return (TA 2 [A 1, A 2]) :: Vect Q (TensorAlgebra ABasis)
A 1*A 2
> f' it
B 1*B 3+B 1*B 4+B 2*B 3+B 2*B 4

So this is variable substitution again. In this case, as f is just a linear map between vector spaces, we can think of it as something like a change of basis of the underlying space. Then f' shows us how the (non-commutative) polynomials defined over the space transform under the change of basis.

Let's just verify that this is a functor. We have to show:
- That fmapTA f is an algebra morphism (ie it is an arrow in k-Alg)
- That fmapTA commutes with the category structure, ie fmapTA id = id, and fmapTA (g . f) = fmapTA g . fmapTA f.

Here's a QuickCheck property:

prop_Functor_Vect_TensorAlgebra (f,g,k,x,y) =
    prop_AlgebraMorphism (fmapTA f') (k,x,y) &&
    (fmapTA id) x == id x &&
    fmapTA (g' . f') x == (fmapTA g' . fmapTA f') x
    where f' = linfun f
          g' = linfun g
          types = (f,g,k,x,y) :: (LinFun Q ABasis BBasis, LinFun Q BBasis CBasis,
                                  Q, Vect Q (TensorAlgebra ABasis), Vect Q (TensorAlgebra ABasis) )

> quickCheck prop_Functor_Vect_TensorAlgebra
+++ OK, passed 100 tests.

So can we declare a Functor instance? Well no, actually. Haskell only allows us to declare type constructors as Functor instances, whereas what we would want to do is declare the type function (\Vect k a -> Vect k (TensorAlgebra a)) as a Functor, which isn't allowed.

Ok, so we have a functor T: k-Vect -> k-Alg, the tensor algebra functor. We also have a forgetful functor going the other way, k-Alg -> k-Vect, which consists in taking an algebra, and simply forgetting that it is an algebra, and seeing only the vector space structure. (As it does at least remember the vector space structure, perhaps we should call this a semi-forgetful, or merely absent-minded functor.)

The cognoscenti will no doubt have seen what is coming next: we have an adjunction, and hence a monad.

How so? Well, it's obvious from its type signature that injectTA is return. For (>>=) / bind, we can define the following:

bindTA :: (Num k, Ord b, Show b) =>
    Vect k (TensorAlgebra a) -> (Vect k a -> Vect k (TensorAlgebra b)) -> Vect k (TensorAlgebra b)
bindTA = flip liftTA

Note that in addition to flipping the arguments, bindTA also imposes a more restrictive signature than liftTA: the target algebra is constrained to be a tensor algebra.

> let f = linear (\(A i) -> case i of 1 -> return (TA 2 [B 1, B 2]);
                                      2 -> return (TA 1 [B 3]) + return (TA 1 [B 4]);
                                      _ -> zerov :: Vect Q (TensorAlgebra BBasis))
> return (TA 2 [A 1, A 2]) `bindTA` f
B 1*B 2*B 3+B 1*B 2*B 4

So the effect of bind is to feed a non-commutative polynomial through a variable substitution.

Monads are meant to satisfy the following monad laws:
- "Left identity": return a >>= f == f a
- "Right identity": m >>= return == m
- "Associativity": (m >>= f) >>= g == m >>= (\x -> f x >>= g)

As usual, we write a QuickCheck property:

prop_Monad_Vect_TensorAlgebra (fmatrix,gmatrix,a,ta)=
    injectTA a `bindTA` f == f a                                  && -- left identity
    ta `bindTA` injectTA == ta                                    && -- right identity
    (ta `bindTA` f) `bindTA` g == ta `bindTA` (\a -> f a `bindTA` g) -- associativity
    where f = linfun fmatrix
          g = linfun gmatrix
          types = (fmatrix,gmatrix,a,ta) :: (LinFun Q ABasis (TensorAlgebra BBasis),
                                             LinFun Q BBasis (TensorAlgebra CBasis),
                                             Vect Q ABasis, Vect Q (TensorAlgebra ABasis) )

> quickCheck prop_Monad_Vect_TensorAlgebra
+++ OK, passed 100 tests.

Once again, we can't actually declare a Monad instance, because our type function (\Vect k a -> Vect k (TensorAlgebra a)) is not a type constructor.

So, we have a functor, and indeed a monad, T: k-Vect -> k-Alg. Now recall that the free vector space construction (\a -> Vect k a) was itself a functor, indeed a monad, from Set -> k-Vect. What happens if we compose these two functors? Why then of course we get a functor, and a monad, from Set -> k-Alg. In Haskell terms, this is a functor a -> Vect k (TensorAlgebra a).

What does this functor look like? Well, relative to a, Vect k (TensorAlgebra a) is the free algebra on a, consisting of all expressions in which the elements of k and the elements of a are combined using (commutative) addition and (non-commutative) multiplication. In other words, the elements of a can be thought of as variable symbols, and Vect k (TensorAlgebra a) as the algebra of non-commutative polynomials in these variables.

Here's the code:

injectTA' :: Num k => a -> Vect k (TensorAlgebra a)
injectTA' = injectTA . return

liftTA' :: (Num k, Ord b, Show b, Algebra k b) =>
     (a -> Vect k b) -> Vect k (TensorAlgebra a) -> Vect k b
liftTA' = liftTA . linear

fmapTA' :: (Num k, Ord b, Show b) =>
    (a -> b) -> Vect k (TensorAlgebra a) -> Vect k (TensorAlgebra b)
fmapTA' = fmapTA . fmap

bindTA' :: (Num k, Ord b, Show b) =>
    Vect k (TensorAlgebra a) -> (a -> Vect k (TensorAlgebra b)) -> Vect k (TensorAlgebra b)
bindTA' = flip liftTA'

The only one of these which might require a little explanation is liftTA'. This works by applying a universal property twice, as shown by the diagram below: first, the universal property of free vector spaces is used to lift a function a -> Vect k (TensorAlgebra b) to a function Vect k a -> Vect k (TensorAlgebra b); then the universal property of the tensor algebra is used to lift that to a function Vect k (TensorAlgebra a) -> Vect k (TensorAlgebra b).

Here's an example, which shows that in the free algebra as in the tensor algebra, bind corresponds to variable substitution:

> let [t,x,y,z] = map injectTA' ["t","x","y","z"] :: [Vect Q (TensorAlgebra String)]
> let f "x" = 1-t^2; f "y" = 2*t; f "z" = 1+t^2
> (x^2+y^2-z^2) `bindTA'` f
0

What is a Coalgebra?

2011-04-23T21:48:00.000+01:00

Last time we saw how to define an algebra structure on a vector space, in terms of category theory. I think perhaps some readers wondered what we gained by using category theory. The answer may be: not much, yet. However, in due course, we would like to understand the connection between quantum algebra and knot theory, and for that, category theory is essential.

This week, I want to look at coalgebras. We already saw, in the case of products and coproducts, how given a structure in category theory, you can define a dual structure by reversing the directions of all the arrows. So it is with algebras and coalgebras.

Recall that an algebra consisted of a k-vector space A together with linear functions
unit :: k -> A
mult :: A⊗A -> A
satisfying two commutative diagrams, associativity and unit.

Well, a coalgebra consists of a k-vector space C together with two linear functions:
counit :: C -> k
comult :: C -> C⊗C
satisfying the following two commutative diagrams:

Coassociativity:

This diagram is actually shorthand for the following diagram:

The isomorphisms at the top are the assocL and assocR isomorphisms that we defined here.

Counit:

These are just the associativity and unit diagrams, but with arrows reversed and relabeled.

Recall that when we say that a diagram commutes, we mean that it doesn't matter which way you follow the arrows, you end up with the same result. In other words, what these diagrams are saying is:

comult⊗id . comult == id⊗comult . comult  (Coassoc)
counit⊗id . comult == unitInL             (Left counit)
id⊗comult . comult == unitInR             (Right counit)

(where unitInL, unitInR are the isomorphisms that we defined in here.)

In HaskellForMaths, recall that we work with free vector spaces over a basis type, so the definition comes out slightly different:

module Math.Algebras.Structures where

...

class Coalgebra k c where
    counit :: Vect k c -> k
    comult :: Vect k c -> Vect k (Tensor c c)

What this definition really says is that c is the basis for a coalgebra. As before, we could try using type families to make this look more like the mathematical definition:

type TensorProd k u v =
    (u ~ Vect k a, v ~ Vect k b) => Vect k (Tensor a b)

class Coalgebra2 k c where
    counit2 :: c -> k
    comult2 :: c -> TensorProd k c c

In this definition, c is the coalgebra itself. However, I'm not going to pursue this approach for now.

What do coalgebras look like? Well, they look a bit like algebras would, if you looked at them through the wrong end of the telescope.

More specifically, given any basis b, define a dual basis as follows:

newtype Dual basis = Dual basis deriving (Eq,Ord)

instance Show basis => Show (Dual basis) where
    show (Dual b) = show b ++ "'"

(For those who know what a dual vector space is - this is it. For those who don't, I'll explain in a minute.)

Then, given an Algebra instance on some finite-dimensional basis b, we can define a Coalgebra instance on Dual b as follows:

1. Write out the unit and mult operations in the algebra as matrices.

For example, in the case of the quaternions, we have
unit:

     1 i j k

1 -> 1 0 0 0

mult:

        1  i  j  k

1⊗1 ->  1  0  0  0
1⊗i ->  0  1  0  0
1⊗j ->  0  0  1  0
1⊗k ->  0  0  0  1
i⊗1 ->  0  1  0  0
i⊗i -> -1  0  0  0
i⊗j ->  0  0  0  1
i⊗k ->  0  0 -1  0
j⊗1 ->  0  0  1  0
j⊗i ->  0  0  0 -1
j⊗j -> -1  0  0  0
j⊗k ->  0  1  0  0
k⊗1 ->  0  0  0  1
k⊗i ->  0  0  1  0
k⊗j ->  0 -1  0  0
k⊗k -> -1  0  0  0

2. Then transpose these two matrices, and use them as the definitions for counit and comult, but replacing each basis element by its dual.

So, in the case of the quaternions, we would get
counit:

      1'

1' -> 1
i' -> 0
j' -> 0
k' -> 0

comult:

     1'⊗1' 1'⊗i' 1'⊗j' 1'⊗k' i'⊗1' i'⊗i' i'⊗j' i'⊗k' j'⊗1' j'⊗i' j'⊗j' j'⊗k' k'⊗1' k'⊗i' k'⊗j' k'⊗k'

1' ->  1     0     0     0      0    -1     0     0     0     0    -1     0     0     0     0    -1
i' ->  0     1     0     0      1     0     0     0     0     0     0     1     0     0    -1     0
j' ->  0     0     1     0      0     0     0    -1     1     0     0     0     0     1     0     0
k' ->  0     0     0     1      0     0     1     0     0    -1     0     0     1     0     0     0

In code, we get:

instance Num k => Coalgebra k (Dual HBasis) where
    counit = unwrap . linear counit'
        where counit' (Dual One) = return ()
              counit' _ = zero
    comult = linear comult'
        where comult' (Dual One) = return (Dual One, Dual One) <+>
                  (-1) *> ( return (Dual I, Dual I) <+> return (Dual J, Dual J) <+> return (Dual K, Dual K) )
              comult' (Dual I) = return (Dual One, Dual I) <+> return (Dual I, Dual One) <+>
                  return (Dual J, Dual K) <+> (-1) *> return (Dual K, Dual J)
              comult' (Dual J) = return (Dual One, Dual J) <+> return (Dual J, Dual One) <+>
                  return (Dual K, Dual I) <+> (-1) *> return (Dual I, Dual K)
              comult' (Dual K) = return (Dual One, Dual K) <+> return (Dual K, Dual One) <+>
                  return (Dual I, Dual J) <+> (-1) *> return (Dual J, Dual I)

unwrap :: Num k => Vect k () -> k
unwrap (V []) = 0
unwrap (V [( (),x)]) = x

(Recall that when we want to think of k as a vector space, we have to represent it as Vect k ().)

We should check that this does indeed define a coalgebra. It's clear, by definition, that counit and comult are linear, as required. Here's a quickcheck property for coassociativity and counit:

prop_Coalgebra x =
    ((comult `tf` id) . comult) x == (assocL . (id `tf` comult) . comult) x  && -- coassociativity
    ((counit' `tf` id) . comult) x == unitInL x                              && -- left counit
    ((id `tf` counit') . comult) x == unitInR x                                 -- right counit
        where counit' = wrap . counit

wrap :: Num k => k -> Vect k ()
wrap 0 = zero
wrap x = V [( (),x)]

(It's a bit awkward that we have to keep wrapping and unwrapping between k and Vect k (). I think we could have avoided this if we had defined counit to have signature Vect k c -> Vect k () instead of Vect k c -> k. To be honest, I think that is probably the right thing to do, so perhaps I'll change it in some future release of HaskellForMaths. What does anyone else think?)

Anyway, here's a quickCheck property to test that the dual quaternions are a coalgebra:

instance Arbitrary b => Arbitrary (Dual b) where
    arbitrary = fmap Dual arbitrary

prop_Coalgebra_DualQuaternion x = prop_Coalgebra x
    where types = x :: Vect Q (Dual HBasis)

> quickCheck prop_Coalgebra_DualQuaternion
+++ OK, passed 100 tests.

So, given an algebra structure on some vector space (basis), we can define a coalgebra structure on the dual vector space (dual basis). But why did we use the dual? Why didn't we just use the above construction to define a coalgebra structure on the vector space itself? For example, for the quaternions, that would look like this:

instance Num k => Coalgebra k HBasis where
    counit = unwrap . linear counit'
        where counit' One = return ()
              counit' _ = zero
    comult = linear comult'
        where comult' One = return (One,One) <+> (-1) *> ( return (I,I) <+> return (J,J) <+> return (K,K) )
              comult' I = return (One,I) <+> return (I,One) <+> return (J,K) <+> (-1) *> return (K,J)
              comult' J = return (One,J) <+> return (J,One) <+> return (K,I) <+> (-1) *> return (I,K)
              comult' K = return (One,K) <+> return (K,One) <+> return (I,J) <+> (-1) *> return (J,I)

Well, we could indeed do that. However, it obscures the underlying mathematics. The point is that an algebra structure on a finite-dimensional vector space gives rise "naturally" to a coalgebra structure on the dual space. It is then also true that a finite-dimensional vector space is naturally isomorphic to its dual, which is I guess why the second construction works.

Okay, so why is the coalgebra structure on the dual vector space "natural"?

Well first, what is a dual vector space anyway? Given a vector space V, the dual space V* is Hom(V,k), the space of linear maps from V to k. If V is finite-dimensional with basis b, then as we have seen, we can define a basis Dual b for V* which is in one-to-one correspondence with b. Given a basis element ei, then the dual basis element Dual ei represents the linear map that sends ei to 1, and any other basis element ej, j /= i, to 0.

It is convenient to define a linear map called the evaluation map, ev: V*⊗V -> k, with ev(a⊗x) = a(x):

ev :: (Num k, Ord b) => Vect k (Tensor (Dual b) b) -> k
ev = unwrap . linear (\(Dual bi, bj) -> delta bi bj *> return ())

-- where delta i j = if i == j then 1 else 0

Then given an element a in V* = Vect k (Dual b), and an element x in V = Vect k b, we can evaluate a(x) by calling ev (a `te` x).

For example:

dual = fmap Dual

> ev $ dual e1 `te` (4 *> e1 <+> 5 *> e2)
4
> ev $ dual e2 `te` (4 *> e1 <+> 5 *> e2)
5
> ev $ dual (e1 <+> 2 *> e2) `te` (4 *> e1 <+> 5 *> e2)
14

Provided V is finite-dimensional, then every element of V* can be expressed in the form dual v, for some v in V.

If we want to turn an element of Vect k (Dual b) into a real Haskell function Vect k b -> k, we can use the following code:

reify :: (Num k, Ord b) => Vect k (Dual b) -> (Vect k b -> k)
reify a x = ev (a `te` x)

For example:

> let f = reify (dual e2)
> f (4 *> e1 <+> 5 *> e2)
5

Now, suppose that we have a linear map f: U -> V between vector spaces. This gives rise to a linear map f*: V* -> U*, by defining:
ev (f* a ⊗ x) = ev (a ⊗ f x)
It turns out that the matrix for f* will be the transpose of the matrix for f.

For, suppose f(ui) = sum [mij vj | j <- ...] and f*(vi*) = sum [m*ij uj* | j <- ...]. Then

ev (f* vk* ⊗ ui) = ev (vk* ⊗ f ui)                                           -- definition of f*
=> ev (sum [m*kl ul* | l <- ...] ⊗ ui) = ev (vk* ⊗ sum [mij vj | j <- ...])  -- expanding f and f*
=> sum [m*kl ev(ul* ⊗ ui) | l <- ...] = sum [mij ev (vk* ⊗ vj) | j <- ...]   -- linearity of ev
=> sum [m*kl (delta l i) | l <- ...] = sum [mij (delta k j) | j <- ...]       -- definition of ev
=> m*ki = mik                                                                 -- definition of delta

So m* is the transpose of m.

Hence we have a contravariant functor * from the category of finite-dimensional vector spaces to itself, which takes a space V to the dual space V*, and a linear map f: U -> V to the dual map f*: V* -> U*. ("Contravariant" just means that it reverses the directions of arrows.)

Hopefully this explains why it is natural to think of an algebra structure on V as giving rise to a coalgebra structure on V*: a coalgebra is basically an algebra but with arrows reversed - and in going from vector spaces to their duals, we reverse arrows.

By the way, the converse is also true: A coalgebra structure on V gives rise to an algebra structure on V*, in the same way.

What is an Algebra?

2011-04-16T20:21:00.001+01:00

Over the last few months, we've spent somewhat longer than I originally expected looking at vector spaces, direct sums and tensor products. I hope you haven't forgotten that the reason we were doing this is because we want to look at quantum algebra, and "quantum groups". What are quantum groups? Well, one thing they are is algebras - so the next thing we need to do is define algebras.

Informally, an algebra is just a vector space which is also a ring (with unit) - or to put it another way, a ring (with unit) which is also a vector space. So a straightforward definition would be, A is an algebra if
(i) A is an additive group (this is required by both vector spaces and rings)
(ii) There is a scalar multiplication smult :: k×A -> A (satisfying some laws, as discussed in a previous post)
(iii) There is a multiplication mult :: A×A -> A, satisfying some laws:
- mult is associative: a(bc) = (ab)c
- mult distributes over addition: a(b+c) = (ab)+(ac), (a+b)c = (ac)+(bc)
(iv) There is a unit :: A, which is an identity for mult: 1a = a = a1

Some examples:
- C is an algebra over R (2-dimensional as a vector space)
- 2×2 matrices over a field k form a k-algebra (4-dimensional)
- polynomials over a field k form a k-algebra (infinite-dimensional)

It would be fairly straightforward to translate these definitions into a Haskell type class as they are. However, we're going to do things slightly differently, for two reasons.

Firstly, we would like to use the language of category theory, and define multiplication and unit in terms of linear maps (arrows in the category of vector spaces). Specifically, we define linear maps:
unit :: k -> A
mult :: A⊗A -> A

We then require that the following diagrams commute:

Associativity:

When we say that this diagram commutes, what it means is that it doesn't matter which way you decide to follow the arrows, the result is the same. Specifically:
mult . mult⊗id == mult . id⊗mult

Unit:

In this case there are two commuting triangles, so we are saying:
mult . unit⊗id == unitOutL
mult . id⊗unit == unitOutR
(where unitOutL, unitOutR are the relevant isomorphisms, which we defined last time.)

But hold on, haven't we forgotten one of the requirements? What about distributivity? Well, notice that the signature of mult was A⊗A -> A, not A×A -> A. A linear map from the tensor product is bilinear in each component (by definition of tensor product) - and bilinearity implies distributivity. So in this version, distributivity is built into the definition of mult. Neat, eh?

The second reason our definition will be different is that in HaskellForMaths, all our vector spaces are free vector spaces over some basis type b: V = Vect k b. Consequently, our algebras will be A = Vect k a, where a is a k-basis for the algebra. Because of this, it will turn out to be more natural to express some things in terms of a, rather than A.

With those forewarnings, here's the HaskellForMaths definition of an algebra:

module Math.Algebras.Structures where

import Math.Algebras.VectorSpace
import Math.Algebras.TensorProduct

class Algebra k a where
    unit :: k -> Vect k a
    mult :: Vect k (Tensor a a) -> Vect k a

In this definition, a represents the basis of the algebra, not the algebra itself, which is A = Vect k a. Recall that we defined a type Tensor a b, which is a basis for Vect k a ⊗ Vect k b.

If we wanted to stay a bit closer to the category theory definition, we could try continuing with the type family approach that we looked at last time:

type TensorProd k u v =
    (u ~ Vect k a, v ~ Vect k b) => Vect k (Tensor a b)

class Algebra2 k a where
    unit2 :: k -> a
    mult2 :: TensorProd k a a -> a

In this definition a is the algebra itself. I'm not going to pursue that approach any further here, but if anyone fancies giving it a go, I'd be interested to hear how they get on.

Anyway, as discussed, unit and mult are required to be linear maps. We could check this using QuickCheck, but in practice we will always define unit and mult in such a way that they are clearly linear.

However, we can write a QuickCheck property to check the other requirements:

prop_Algebra (k,x,y,z) =
    (mult . (id `tf` mult)) (x `te` (y `te` z)) ==
         (mult . (mult `tf` id)) ((x `te` y) `te` z)              && -- associativity
    unitOutL (k' `te` x) == (mult . (unit' `tf` id)) (k' `te` x)  && -- left unit
    unitOutR (x `te` k') == (mult . (id `tf` unit')) (x `te` k')     -- right unit
    where k' = k *> return ()

(Recall that when we wish to consider k as a vector space, we represent it as the free vector space Vect k ().)

When we have an algebra, then we have a ring, so we can define a Num instance:

instance (Num k, Eq b, Ord b, Show b, Algebra k b) => Num (Vect k b) where
    x+y = x <+> y
    negate x = neg x
    x*y = mult (x `te` y)
    fromInteger n = unit (fromInteger n)
    abs _ = error "Prelude.Num.abs: inappropriate abstraction"
    signum _ = error "Prelude.Num.signum: inappropriate abstraction"

This means that when we have an algebra, we'll be able to write expressions using the usual arithmetic operators.

Okay, so how about some examples of algebras. Well, we mentioned the complex numbers as an algebra over the reals, but let's go one better and define the quaternion algebra. This is a four dimensional algebra over any field k, which is generated by {1, i, j, k}, satisfying the relations i^2 = j^2 = k^2 = ijk = -1. (It follows, for example that (ijk)k = (-1)k, so ij(k^2) = -k, so ij = k.)

Here's the code. First, we define our basis. The quaternions are traditionally denoted H, after Hamilton, who discovered them:

data HBasis = One | I | J | K deriving (Eq,Ord)

type Quaternion k = Vect k HBasis

i,j,k :: Num k => Quaternion k
i = return I
j = return J
k = return K

instance Show HBasis where
    show One = "1"
    show I = "i"
    show J = "j"
    show K = "k"

instance (Num k) => Algebra k HBasis where
    unit x = x *> return One
    mult = linear m
         where m (One,b) = return b
               m (b,One) = return b
               m (I,I) = unit (-1)
               m (J,J) = unit (-1)
               m (K,K) = unit (-1)
               m (I,J) = return K
               m (J,I) = -1 *> return K
               m (J,K) = return I
               m (K,J) = -1 *> return I
               m (K,I) = return J
               m (I,K) = -1 *> return J

Note that unit and mult are both linear by definition.

Let's just check that the code works as expected:

> :l Math.Algebras.Quaternions
> i^2
-1
> j^2
-1
> i*j
k

Now, are we sure that the quaternions are an algebra? Well, it's clear from the definition that the left and right unit conditions hold - see the lines m (One,b) = m (b,One) = return b. But it's not obvious that the associativity condition holds, so perhaps we should quickCheck:

instance Arbitrary HBasis where
    arbitrary = elements [One,I,J,K]

prop_Algebra_Quaternion (k,x,y,z) = prop_Algebra (k,x,y,z)
    where types = (k,x,y,z) :: (Q, Quaternion Q, Quaternion Q, Quaternion Q)

> quickCheck prop_Algebra_Quaternion
+++ OK, passed 100 tests.

Ok, how about 2*2 matrices. These form a four dimensional algebra generated by the elementary matrices {e11, e12, e21, e22}, where eij is the matrix with a 1 in the (i,j) position, and 0s elsewhere:

data Mat2 = E2 Int Int deriving (Eq,Ord,Show)

instance Num k => Algebra k Mat2 where
   unit x = x *> V [(E2 i i, 1) | i <- [1..2] ]
   mult = linear mult' where
        mult' (E2 i j, E2 k l) = delta j k *> return (E2 i l)

delta i j | i == j    = 1
          | otherwise = 0

Notice the way that we only have to define multiplication on our basis elements, the elementary matrices eij, and the rest follows by bilinearity. Notice also that unit and mult are linear by definition. Let's just sanity check that this works as expected:

> :l Math.Algebras.Matrix
> unit 1 :: Vect Q Mat2
E2 1 1+E2 2 2
> let a = 2 *> return (E2 1 2) + 3 *> return (E2 2 1)
> a^2
6E2 1 1+6E2 2 2

It's straightforward to define an Arbitrary instance for Mat2, and quickCheck that it satisfies the algebra conditions.

Finally, what about polynomials? Let's go one better, and define polynomials in more than one variable. An obvious basis for these polynomials as a vector space is the set of monomials. For example, polynomials in x,y,z are a vector space on the basis { x^i y^j z^k | i <- [0..], j <- [0..], k <- [0..] }.

Recall that our vector space code requires that the basis be an Ord instance, so we need to define an ordering on monomials. There are many ways to do this. We'll use the graded lex or glex ordering, which says that monomials of higher degree sort before those of lower degree, and among those of equal degree, lexicographic (alphabetical) ordering applies.

Given any set X of variables, we can construct the polynomial algebra over X as the vector space with basis the monomials in X. In the example above, we had X = {x,y,z}. For this reason, we'll allow our glex monomials to be polymorphic in the type of the variables. (In practice though, as you will see shortly, we will often just use String as the type of our variables.)

So a glex monomial over variables v is basically just a list of powers of elements of v:

data GlexMonomial v = Glex Int [(v,Int)] deriving (Eq)
-- The initial Int is the degree of the monomial. Storing it speeds up equality tests and comparisons

For example x^3 y^2 would be represented as Glex 5 [("x",3),("y",2)].

instance Ord v => Ord (GlexMonomial v) where
    compare (Glex si xis) (Glex sj yjs) =
        compare (-si, [(x,-i) | (x,i) <- xis]) (-sj, [(y,-j) | (y,j) <- yjs])
-- all the minus signs are to make things sort in the right order

[There's a bug in the HaskellForMaths v0.3.2 version of this Ord instance - the above code, which fixes it, will be in the next release.]

I won't bore you with the Show instance - it's a bit fiddly.

The Algebra instance uses the fact that monomials form a monoid under multiplication, with unit 1. Given any monoid, we can form the free vector space having the monoid elements as basis, and then lift the unit and multiplication in the monoid into the vector space, thus forming an algebra called the monoid algebra. Here's the construction:

instance (Num k, Ord v) => Algebra k (GlexMonomial v) where
    unit x = x *> return munit
        where munit = Glex 0 []
    mult xy = nf $ fmap (\(a,b) -> a `mmult` b) xy
        where mmult (Glex si xis) (Glex sj yjs) = Glex (si+sj) $ addmerge xis yjs

(In the addmerge function, we ensure that provided the variables were listed in ascending order in both inputs, then they are still so in the output.)

Finally, here's a convenience function for injecting a variable into the polynomial algebra:

glexVar v = V [(Glex 1 [(v,1)], 1)]

Then for example, we can do the following:

type GlexPoly k v = Vect k (GlexMonomial v)

> let x = glexVar "x" :: GlexPoly Q String
> let y = glexVar "y" :: GlexPoly Q String
> let z = glexVar "z" :: GlexPoly Q String
> (x+y+z)^3
x^3+3x^2y+3x^2z+3xy^2+6xyz+3xz^2+y^3+3y^2z+3yz^2+z^3

I hope you were impressed at how easy that all was. The foundation of free vector spaces, tensor products and algebras is only around a hundred lines of code. It then takes just another dozen lines or so to define an algebra: just define a basis, and define how the basis elements multiply - the rest follows by linearity.

Tensor Products, part 2: Monoids and Arrows

2011-03-18T10:23:00.000Z

[New release, HaskellForMaths v0.3.2, available on Hackage]

Last time we looked at the tensor product of free vector spaces. Given A = Vect k a, B = Vect k b, then the tensor product A⊗B can be represented as Vect k (a,b). As we saw, the tensor product is the "mother of all bilinear functions".

In the HaskellForMaths library, I have defined a couple of type synonyms for direct sum and tensor product:

type DSum a b = Either a b
type Tensor a b = (a,b)

This means that in type signatures we can write the type of a direct sum as Vect k (DSum a b), and of a tensor product as Vect k (Tensor a b). The idea is that this will remind us what we're dealing with, and make things clearer.

During development, I initially called the tensor type TensorBasis. In maths, tensor product is thought of as an operation on vector spaces - A⊗B - rather than on their bases. It would be nicer if we could define direct sum and tensor product as operators on the vector spaces themselves, rather than their bases.

Well, we can have a go, something like this:

{-# LANGUAGE TypeFamilies, RankNTypes #-}

type DirectSum k u v =
    (u ~ Vect k a, v ~ Vect k b) => Vect k (DSum a b)

type TensorProd k u v =
    (u ~ Vect k a, v ~ Vect k b) => Vect k (Tensor a b)

type En = Vect Q EBasis

This appears to work:

$ ghci -XTypeFamilies
...
> :l Math.Test.TAlgebras.TTensorProduct
...
> e1 `te` e2 :: TensorProd Q En En
(e1,e2)

I'll reserve judgement. (Earlier in the development of the quantum algebra code for HaskellForMaths, I tried something similar to this, and ran into problems later on - but I can't now remember exactly what I did, so perhaps this will work.)

Okay, so what can we do with tensor products. Well first, given vectors u in A = Vect k a and v in B = Vect k b, we can form their tensor product, u⊗v, an element of A⊗B = Vect k (Tensor a b). To calculate u⊗v, we use the bilinearity of tensor product to reduce the tensor product of arbitrary vectors to a linear combination of tensor products of basis elements:
(x1 a1 + x2 a2)⊗(y1 b1 + y2 b2) = x1 y1 a1⊗b1 + x1 y2 a1⊗b2 + x2 y1 a2⊗b1 + x2 y2 a2⊗b2
Here's the code:

te :: Num k => Vect k a -> Vect k b -> Vect k (Tensor a b)
te (V us) (V vs) = V [((a,b), x*y) | (a,x) <- us, (b,y) <- vs]

This is in essence just the "tensor" function from last time, but rewritten to take its two inputs separately rather than in a direct sum. Mnemonic: "te" stands for "tensor product of elements". Note that the definition respects normal form: provided the inputs are in normal form (the as and bs are in order, and the xs and ys are non-zero), then so is the output.

Associativity

We can form tensor products of tensor products, such as A⊗(B⊗C) = Vect k (a,(b,c)), and likewise (A⊗B)⊗C = Vect k ((a,b),c). These two are isomorphic as vector spaces. This is obvious if you think about it in the right way. Recall from last week that we can think of elements of A⊗B = Vect k (a,b) as 2-dimensional matrices with rows indexed by a, columns indexed by b, and entries in k. Well A⊗B⊗C (we can drop the parentheses as it makes no difference) is the space of three-dimensional matrices, with one dimension indexed by a, one by b, and one by c.

We can define isomorphisms either way with the following Haskell code:

assocL :: Vect k (Tensor u (Tensor v w)) -> Vect k (Tensor (Tensor u v) w)
assocL = fmap ( \(a,(b,c)) -> ((a,b),c) )

assocR :: Vect k (Tensor (Tensor u v) w) -> Vect k (Tensor u (Tensor v w))
assocR = fmap ( \((a,b),c) -> (a,(b,c)) )

It's clear that these functions are linear, since they're defined using fmap. It's also clear that they are bijections, since they are mutually inverse. Hence they are the required isomorphisms.

Unit

Last time we saw that the field k is itself a vector space, which can be represented as the free vector space Vect k (). What happens if we take the tensor product k⊗A of the field with some other vector space A = Vect k a? Well, if you think about it in terms of matrices, Vect k () is a one-dimensional vector space, so Vect k ((),a) will be a 1*n matrix (where n is the number of basis elements in a). But a 1*n matrix looks just the same as an n-vector:

    a1 a2 ...      a1 a2 ...
() ( .  .    ) ~= ( .  .    )

So we should expect that k⊗A = Vect k ((),a) is isomorphic to A = Vect k a. And indeed it is - here are the relevant isomorphisms:

unitInL = fmap ( \a -> ((),a) )

unitOutL = fmap ( \((),a) -> a )

unitInR = fmap ( \a -> (a,()) )

unitOutR = fmap ( \(a,()) -> a )

So tensor product is associative, and has a unit. In other words, vector spaces form a monoid under tensor product.

Tensor product of functions

Given linear functions f: A -> A', g: B -> B', we can define a linear function f⊗g: A⊗B -> A'⊗B' by
(f⊗g)(a⊗b) = f(a)⊗g(b)

Exercise: Prove that f⊗g is linear

Here's the Haskell code:

tf :: (Num k, Ord a', Ord b') => (Vect k a -> Vect k a') -> (Vect k b -> Vect k b')
   -> Vect k (Tensor a b) -> Vect k (Tensor a' b')
tf f g (V ts) = sum [x *> te (f $ return a) (g $ return b) | ((a,b), x) <- ts]
    where sum = foldl add zero

(Mnemonic: "tf" stands for "tensor product of functions".)

Let's just check that this is linear:

prop_Linear_tf ((f,g),k,(a1,a2,b1,b2)) = prop_Linear (linfun f `tf` linfun g) (k, a1 `te` b1, a2 `te` b2)
    where types = (f,g,k,a1,a2,b1,b2) :: (LinFun Q ABasis SBasis, LinFun Q BBasis TBasis, Q,
                                          Vect Q ABasis, Vect Q ABasis, Vect Q BBasis, Vect Q BBasis)

> quickCheck prop_Linear_tf
+++ OK, passed 100 tests.

So we now have tensor product operations on objects and on arrows. In each case, tensor product takes a pair of objects/arrows, and returns a new object/arrow.

There is a product category k-Vect×k-Vect, consisting of pairs of objects and pairs of arrows from k-Vect. The identity arrow is defined to be (id,id), and composition is defined by (f,g) . (f',g') = (f . f', g . g'). Given these definitions, it turns out that tensor product is a functor from k-Vect×k-Vect to k-Vect. (Another way to say this is that tensor product is a bifunctor in the category of vector spaces.)

Recall that a functor is just a map that "commutes" with the category operations, id and . (composition).
So the conditions for tensor product to be a functor are:
id⊗id = id
(f' . f)⊗(g' . g) = (f'⊗g') . (f⊗g)

Both of these follow immediately from the definition of f⊗g that was given above. However, just in case you don't believe me, here's a quickCheck property to prove it:

prop_TensorFunctor ((f1,f2,g1,g2),(a,b)) =
    (id `tf` id) (a `te` b) == id (a `te` b) &&
    ((f' . f) `tf` (g' . g)) (a `te` b) == ((f' `tf` g') . (f `tf` g)) (a `te` b)
    where f = linfun f1
          f' = linfun f2
          g = linfun g1
          g' = linfun g2
          types = (f1,f2,g1,g2,a,b) :: (LinFun Q ABasis ABasis, LinFun Q ABasis ABasis,
                                        LinFun Q BBasis BBasis, LinFun Q BBasis BBasis,
                                        Vect Q ABasis, Vect Q BBasis)

> quickCheck prop_TensorFunctor
+++ OK, passed 100 tests.

We can think of composition as doing things in series, and tensor as doing things in parallel.

Then the second bifunctor condition can be paraphrased as "Doing things in parallel, in series, is the same as doing things in series, in parallel", as represented by the following diagram.

You might recall that there are a couple of Haskell type classes for that kind of thing. The Category typeclass from Control.Category is about doing things in series. Here is the definition:

class Category cat where
    id :: cat a a
    (.) :: cat b c -> cat a b -> cat a c

The Arrow typeclass from Control.Arrow is about doing things in parallel:

class Category arr => Arrow arr where
    arr :: (a -> b) -> arr a b
    first :: arr a b -> arr (a,c) (b,c)
    second :: arr a b -> arr (c,a) (c,b)
    (***) :: arr a b -> arr a' b' -> arr (a,a') (b,b')
    (&&&) :: arr a b -> arr a b' -> arr a (b,b')

Intuitively, linear functions (Vect k a -> Vect k b) are arrows, via the definitions:
id = id
(.) = (.)
arr = fmap
first f = f `tf` id
second f = id `tf` f
f *** g = f `tf` g
(f &&& g) a = (f `tf` g) (a `te` a)

However, in order to define an Arrow instance we'll have to wrap the functions in a newtype.

import Prelude as P
import Control.Category as C
import Control.Arrow

newtype Linear k a b = Linear (Vect k a -> Vect k b)

instance Category (Linear k) where
    id = Linear P.id
    (Linear f) . (Linear g) = Linear (f P.. g)

instance Num k => Arrow (Linear k) where
    arr f = Linear (fmap f)
    first (Linear f) = Linear f *** Linear P.id
    second (Linear f) = Linear P.id *** Linear f
    Linear f *** Linear g = Linear (f `tf2` g)
        where tf2 f g (V ts) = V $ concat
                  [let V us = x *> te (f $ return a) (g $ return b) in us | ((a,b), x) <- ts]
    Linear f &&& Linear g = (Linear f *** Linear g) C.. Linear (\a -> a `te` a)

Note that we can't use tf directly, because it requires Ord instances for a and b, and Haskell doesn't give us a way to require these. For this reason we define a tf2 function, which is equivalent except that it doesn't guarantee that results are in normal form.

There is loads of other stuff I could talk about:
Exercise: Show that direct sum is also a monoid, with the zero vector space as its identity. (Write Haskell functions for the necessary isomorphisms.)
Exercise: Show that tensor product distributes over direct sum - A⊗(B⊕C) is isomorphic to (A⊗B)⊕(A⊗C). (Write the isomorphisms.)
Exercise: Show that given f: A->A', g: B->B', it is possible to define a linear function f⊕g: A⊕B->A'⊕B' by (f⊕g)(a⊕b) = f(a)⊕g(b). (Write a dsumf function analogous to tf.)

There is another arrow related typeclass called ArrowChoice. It represents arrows where you have a choice of doing either one thing or another thing:

class Arrow arr => ArrowChoice arr where
    left :: arr a b -> arr (Either a c) (Either b c)
    right :: arr a b -> arr (Either c a) (Either c b)
    (+++) :: arr a a' -> arr b b' -> arr (Either a b) (Either a' b')
    (|||) :: arr a c -> arr b c -> arr (Either a b) c

Exercise: Show that the dsumf function can be used to give an ArrowChoice instance for linear functions, where the left summand goes down one path and the right summand down another.

Tensor products of vector spaces, part 1

2011-02-21T20:13:00.001Z

A little while back on this blog, we defined the free k-vector space over a type b:

newtype Vect k b = V [(b,k)] deriving (Eq,Ord)
Elements of Vect k b are k-linear combinations of elements of b.

Whenever we have a mathematical structure like this, we want to know about building blocks and new-from-old constructions.

We already looked at one new-from-old construction: given free k-vector spaces A = Vect k a and B = Vect k b, we can construct their direct sum A⊕B = Vect k (Either a b).

We saw that the direct sum is both the product and the coproduct in the category of free vector spaces - which means that it is the object which satisfies the universal properties implied by the following two diagrams:

So we have injections i1, i2 : Vect k a, Vect k b -> Vect k (Either a b), to put elements of A and B into the direct sum A⊕B, and projections p1, p2 : Vect k (Either a b) -> Vect k a, Vect k b to take them back out again.

However, there is another obvious new-from-old construction: Vect k (a,b). What does this represent?

In order to answer that question, we need to look at bilinear functions. The basic idea of a bilinear function is that it is a function of two arguments, which is linear in each argument. So we might start by looking at functions f :: Vect k a -> Vect k b -> Vect k t.

However, functions of two arguments don't really sit very well in category theory, where arrows are meant to have a single source. (We can handle functions of two arguments in multicategories, but I don't want to go there just yet.) In order to stay within category theory, we need to combine the two arguments into a single argument, using the direct sum construction. So instead of looking at functions f :: Vect k a -> Vect k b -> Vect k t, we will look at functions f :: Vect k (Either a b) -> Vect k t.

To see that they are equivalent, recall from last time that Vect k (Either a b) is isomorphic to (Vect k a, Vect k b), via the isomorphisms:

to :: (Vect k a, Vect k b) -> Vect k (Either a b)
to = \(u,v) -> i1 u <+> i2 v
from :: Vect k (Either a b) -> (Vect k a, Vect k b)
from = \uv -> (p1 uv, p2 uv)

So in going from f :: Vect k a -> Vect k b -> Vect k t to f :: Vect k (Either a b) -> Vect k t, we're really just uncurrying.

Ok, so suppose we are given f :: Vect k (Either a b) -> Vect k t. It helps to still think of this as a function of two arguments, even though we've wrapped them up together in either side of a direct sum. Then we say that f is bilinear, if it is linear in each side of the direct sum. That is:
- for any fixed a0 in A, the function f_a0 :: Vect k b -> Vect k t, f_a0 = \b -> f (i1 a0 <+> i2 b) is linear
- for any fixed b0 in B, the function f_b0 :: Vect k a -> Vect k t, f_b0 = \a -> f (i1 a <+> i2 b0) is linear

Here's a QuickCheck property to test whether a function is bilinear:

prop_Bilinear :: (Num k, Ord a, Ord b, Ord t) =>
     (Vect k (Either a b) -> Vect k t) -> (k, Vect k a, Vect k a, Vect k b, Vect k b) -> Bool
prop_Bilinear f (k,a1,a2,b1,b2) =
    prop_Linear (\b -> f (i1 a1 <+> i1 b)) (k,b1,b2) &&
    prop_Linear (\a -> f (i1 a <+> i1 b1)) (k,a1,a2)

prop_BilinearQn f (a,u1,u2,v1,v2) = prop_Bilinear f (a,u1,u2,v1,v2)
    where types = (a,u1,u2,v1,v2) :: (Q, Vect Q EBasis, Vect Q EBasis, Vect Q EBasis, Vect Q EBasis)

What are some examples of bilinear functions?

Well, perhaps the most straightforward is the dot product of vectors. If our vector spaces A and B are the same, then we can define the dot product:

dot0 uv = sum [ if a == b then x*y else 0 | (a,x) <- u, (b,y) <- v]
    where V u = p1 uv
          V v = p2 uv

However, as it stands, this won't pass our QuickCheck property - because it has the wrong type! This has the type dot0 :: Vect k (Either a b) -> k, whereas we need something of type Vect k (Either a b) -> Vect k t.

Now, it is of course true that k is a k-vector space. However, as it stands, it's not a free k-vector space over some basis type t. Luckily, this is only a technicality, which is easily fixed. When we want to consider k as itself a (free) vector space, we will take t = (), the unit type, and equate k with Vect k (). Since the type () has only a single inhabitant, the value (), then Vect k () consists of scalar multiples of () - so it is basically just a single copy of k itself. The isomorphism between k and Vect k () is \k -> k *> return ().

Okay, so now that we know how to represent k as a free k-vector space, we can define dot product again:

dot1 uv = nf $ V [( (), if a == b then x*y else 0) | (a,x) <- u, (b,y) <- v]
    where V u = p1 uv
          V v = p2 uv

This now has the type dot1 :: Vect k (Either a b) -> Vect k (). Here's how you use it:

> dot1 ( i1 (e1 <+> 2 *> e2) <+> i2 (3 *> e1 <+> e2) )
5()

(So thinking of our function as a function of two arguments, what we do is use i1 to inject the first argument into the left hand side of the direct sum, and i2 to inject the second argument into the right hand side.)

So we can now use the QuickCheck property:

> quickCheck (prop_BilinearQn dot1)
+++ OK, passed 100 tests.

Another example of a bilinear function is polynomial multiplication. Polynomials of course form a vector space, with basis {x^i | i <- [0..] }. So we could define a type to represent the monomials x^i, and then form the polynomials as the free vector space in the monomials. In a few weeks we will do that, but for the moment, to save time, let's just use our existing EBasis type, and take E i to represent x^i. Then polynomial multiplication is the following function:

polymult1 uv = nf $ V [(E (i+j) , x*y) | (E i,x) <- u, (E j,y) <- v]
    where V u = p1 uv
          V v = p2 uv

Let's just convince ourselves that this is polynomial multiplication:

> polymult1 (i1 (e 0 <+> e 1) <+> i2 (e 0 <+> e 1))
e0+2e1+e2

So this is just our way of saying that (1+x)*(1+x) = 1+2x+x^2.

Again, let's verify that this is bilinear:

> quickCheck (prop_BilinearQn polymult1)
+++ OK, passed 100 tests.

So what's all this got to do with Vect k (a,b)? Well, here's another bilinear function:

tensor :: (Num k, Ord a, Ord b) => Vect k (Either a b) -> Vect k (a, b)
tensor uv = nf $ V [( (a,b), x*y) | (a,x) <- u, (b,y) <- v]
    where V u = p1 uv; V v = p2 uv

> quickCheck (prop_BilinearQn tensor)
+++ OK, passed 100 tests.

So this "tensor" function takes each pair of basis elements a, b in the input to a basis element (a,b) in the output. The thing that is interesting about this bilinear function is that it is in some sense "the mother of all bilinear functions". Specifically, you can specify a bilinear function completely by specifying what happens to each pair (a,b) of basis elements. It follows that any bilinear function f :: Vect k (Either a b) -> Vect k t can be factored as f = f' . tensor, where f' :: Vect k (a,b) -> Vect k t is the linear function having the required action on the basis elements (a,b) of Vect k (a,b).

For example:

bilinear :: (Num k, Ord a, Ord b, Ord c) =>
    ((a, b) -> Vect k c) -> Vect k (Either a b) -> Vect k c
bilinear f = linear f . tensor

dot = bilinear (\(a,b) -> if a == b then return () else zero)

polymult = bilinear (\(E i, E j) -> return (E (i+j)))

We can check that these are indeed the same functions as we were looking at before:

> quickCheck (\x -> dot1 x == dot x)
+++ OK, passed 100 tests.
> quickCheck (\x -> polymult1 x == polymult x)
+++ OK, passed 100 tests.

So Vect k (a,b) has a special role in the theory of bilinear functions. If A = Vect k a, B = Vect k b, then we write A⊗B = Vect k (a,b) (pronounced "A tensor B").

[By the way, it's possible that this diagram might upset category theorists - because the arrows in the diagram are not all arrows in the category of vector spaces. Specifically, note that bilinear maps are not, in general, linear. We'll come back to this in a moment.]

So a bilinear map can be specified by its action on the tensor basis (a,b). This corresponds to writing out matrices. To specify any bilinear map Vect k (Either a b) -> Vect k t, you write out a matrix with rows indexed by a, columns indexed by b, and entries in Vect k t.

      b1  b2 ...
 a1 (t11 t12 ...)
 a2 (t21 t22 ...)
... (...        )

So this says that (ai,bj) is taken to tij. Then given an element of A⊕B = Vect k (Either a b), which we can think of as a vector (x1 a1 + x2 a2 + ...) in A = Vect k a together with a vector (y1 b1 + y2 b2 + ...) in B = Vect k b, then we can calculate its image under the bilinear map by doing matrix multiplication as follows:

 a1 a2 ...        b1  b2 ...
(x1 x2 ...)  a1 (t11 t12 ...)  b1 (y1)
             a2 (t21 t22 ...)  b2 (y2)
            ... (...        ) ... (...)

(Sorry, this diagram might be a bit confusing. The ai, bj are labeling the rows and columns. The xi are the entries in a row vector in A, the yj are the entries in a column vector in B, and the tij are the entries in the matrix.)

So xi ai <+> yj bj goes to xi yj tij.

For example, dot product corresponds to the matrix:

(1 0 0)
(0 1 0)
(0 0 1)

Polynomial multiplication corresponds to the matrix:

    e0 e1 e2 ...
e0 (e0 e1 e2 ...)
e1 (e1 e2 e3 ...)
e2 (e2 e3 e4 ...)
...

A matrix with entries in T = Vect k t is just a convenient way of specifying a linear map from A⊗B = Vect k (a,b) to T.

Indeed, any matrix, provided that all the entries are in the same T, defines a bilinear function. So bilinear functions are ten-a-penny.

Now, I stated above that bilinear functions are not in general linear. For example:

> quickCheck (prop_Linear polymult)
*** Failed! Falsifiable (after 2 tests and 2 shrinks): 
(0,Right e1,Left e1)

What went wrong? Well:

> polymult (Right e1)
0
> polymult (Left e1)
0
> polymult (Left e1 <+> Right e1)
e2

So we fail to have f (a <+> b) = f a <+> f b, which is one of the requirements of a linear function.

Conversely, it's also important to realise that linear functions (on Vect k (Either a b)) are not in general bilinear. For example:

> quickCheck (prop_BilinearQn id)
*** Failed! Falsifiable (after 2 tests): 
(1,0,0,e1,0)

The problem here is:

> id $ i1 (zero <+> zero) <+> i2 e1
Right e1
> id $ (i1 zero <+> i2 e1) <+> (i1 zero <+> i2 e1)
2Right e1

So we fail to have linearity in the left hand side (or the right for that matter).

Indeed we can kind of see that linearity and bilinearity are in conflict.
- Linearity requires that f (a1 <+> a2 <+> b) = f a1 <+> f a2 <+> f b
- Bilinearity requires that f (a1 <+> a2 <+> b) = f (a1 <+> b) <+> f (a2 <+> b)

Exercise: Find a function which is both linear and bilinear.

Products of lists and vector spaces

2011-02-01T21:47:00.000Z

Last time, we looked at coproducts - of sets/types, of lists, and of free vector spaces. I realised afterwards that there were a couple more things I should have said, but forgot.

Recall that the coproduct of A and B is an object A+B, together with injections i1: A -> A+B, i2: B-> A+B, with the property that whenever we have arrows f: A -> T, g: B -> T, they can be factored through A+B to give an arrow f+g, satisfying f+g . i1 = f, f+g . i2 = g.

Firstly then, I forgot to say why we called the coproduct A+B with a plus sign. Well, it's because, via the injections i1 and i2, it contains (a copy of) A and (a copy of) B. So it's a bit like a sum of A and B.

Second, I forgot to say that in the case of vector spaces, the coproduct is called the direct sum, and has its own special symbol A⊕B.

Okay, so this time I want to look at products.

Suppose we have objects A and B in some category. Then their product (if it exists) is an object A×B, together with projections p1: A×B -> A, p2: A×B -> B, with the following universal property: whenever we have arrows f: S -> A and g: S -> B, then they can be factored through A×B to give an arrow f×g: S -> A×B, such that f = p1 . f×g, g = p2 . f×g.

(The definitions of product and coproduct are dual to one another - the diagrams are the same but with the directions of the arrows reversed.)

In the category Set, the product of sets A and B is their Cartesian product A×B. In the category Hask, of course, the product of types a and b is written (a,b), p1 is called fst, and p2 is called snd.

We can then define the required product map as:
(f .*. g) x = (f x, g x)

Then it should be clear that fst . (f .*. g) = f, and snd . (f .*. g) = g, as required.

Okay, so what do products look like in the category of lists (free monoids)? (Recall that in this category, the arrows are required to be monoid homomorphisms, meaning that f [] = [] and f (xs++ys) = f xs ++ f ys. It follows that we can express f = concatMap f', for some f'.)

Well, the obvious thing to try as the product is the Cartesian product ([a],[b]). Is the Cartesian product of two monoids a monoid? Well yes it is actually. We could give it a monoid structure as follows:

(as1, bs1) ++ (as2, bs2) = (as1++as2, bs1++bs2)
[] = ([],[])

This isn't valid Haskell code of course. It's just my shorthand way of expressing the following code from Data.Monoid:

instance Monoid [a] where
        mempty  = []
        mappend = (++)

instance (Monoid a, Monoid b) => Monoid (a,b) where
        mempty = (mempty, mempty)
        (a1,b1) `mappend` (a2,b2) =
                (a1 `mappend` a2, b1 `mappend` b2)

From these two instances, it follows that ([a],[b]) is a monoid, with monoid operations equivalent to those I gave above. (In particular, it's clear that the construction satisfies the monoid laws: associativity of ++, identity of [].)

But it feels like there's something unsatisfactory about this. Wouldn't it be better for the product of list types [a] and [b] to be another list type [x], for some type x?

Our first thought might be to try [(a.b)]. The product map would then need to be something like \ss -> zip (f ss) (g ss). However, we quickly see that this won't work: what if f ss and g ss are not the same length.

What else might work? Well, if you think of ([a],[b]) as some as on the left and some bs on the right, then the answer should spring to mind. Let's try [Either a b]. We can then define:

p1 xs = [x | Left x <- xs] -- this is doing a filter and a map at the same time
p2 xs = [x | Right x <- xs]

with the product map

f×g = \ss -> map Left (f ss) ++ map Right (g ss)

Then it is clear that p1 . f×g = f and p2 . f×g = g, as required.

What is the relationship between ([a],[b]) and [Either a b]? Well, ([a],[b]) looks a bit like a subset of [Either a b], via the injection i (as,bs) = map Left as ++ map Right bs. However, this injection is not a monoid homomorphism, since
i ([],[b1] ++ ([a1],[]) /= i ([],[b1]) ++ i ([a1],[])
So ([a],[b]) is not a submonoid of [Either a b].

On the other hand, there is a projection p :: [Either a b] -> ([a],[b]), p xs = (p1 xs, p2 xs). This is a monoid homomorphism, so ([a],[b]) is a quotient of [Either a b].

So which is the right answer? Which of ([a],[b]) and [Either a b] is really the product of [a] and [b]?

Well, it depends. It depends which category we think we're working in. If we're working in the category of monoids, then it is ([a],[b]). However, if we're working in the category of free monoids (lists), then it is [Either a b].

You see, ([a],[b]) is not a free monoid. What does this mean? Well, it basically means it's not a list. But how do we know that ([a],[b]) isn't equivalent to some list? And anyway, what does "free" mean in free monoid?

"Free" is a concept that can be applied to many algebraic theories, not just monoids. There is more than one way to define it.

An algebraic theory defines various constants and operations. In the case of monoids, there is one constant - which we may variously call [] or mempty or 0 or 1 - and one operation - ++ or mappend or + or *. Now, a given monoid may turn out to be generated by some subset of its elements - meaning that every element of the monoid can be equated with some expression in the generators, constants, and operations.

For example, the monoid of natural numbers is generated by the prime numbers: every natural number is equal to some expression in 1, *, and the prime numbers. The monoid [x] is generated by the singleton lists: every element of [x] is equal to some expression in [], ++, and the singleton lists. By a slight abuse of notation, we can say that [x] is generated by x - by identifying the singleton lists with the image of x under \x -> [x].

Then we say that a monoid is free on its generators if there are no relations among its elements other than those implied by the monoid laws. That is, no two expressions in the generators, constants, and operators are equal to one another, unless it is as a consequence of the monoid laws.

For example, suppose it happens that
(x ++ y) ++ z = x ++ (y ++ z)
That's okay, because it follows from the monoid laws. On the other hand, suppose that
x ++ y = y ++ x
This does not follow from the monoid laws (unless x = [] or y = []), so is a non-trivial relation. (Thus the natural numbers under multiplication are not a free monoid - because they're commutative.)

What about our type ([a],[b]) then? Well consider the following relations:

(as,[]) ++ ([],bs) = ([],bs) ++ (as,[])
(as1++as2,bs1) ++ ([],bs2) = (as1,bs1) ++ (as2,bs2) = (as1,[]) ++ (as2,bs1++bs2)

We have commutativity relations between the [a] and [b] parts of the product. Crucially, these relations are not implied by the monoid structure alone. So intuitively, we can see that ([a],[b]) is not free.

The "no relations" definition of free is the algebraic way to think about it. However, there is also a category theory way to define it. The basic idea is that if a monoid is free on its generators, then given any other monoid with the same generators, we can construct it as a homomorphic image of our free monoid, by "adding" the appropriate relations.

In order to express this properly, we're going to need to use some category theory, and specifically the concept of the forgetful functor. Recall that given any algebraic category, such as Mon (monoids), there is a forgetful functor U: Mon -> Set, which consists in simply forgetting the algebraic structure. U takes objects to their underlying sets, and arrows to the underlying functions. In Haskell, U: Mon -> Hask consists in forgetting that our objects (types) are monoids, and forgetting that our arrows (functions) are monoid homomorphisms. (As a consequence, U is syntactically invisible in Haskell. However, to properly understand the definition of free, we have to remember that it's there.)

Then, given an object x (the generators), a free monoid on x is a monoid y, together with a function i: x -> U y, such that whenever we have an object z in Mon and a function f': x -> U z, then we can lift it to a unique arrow f: y -> z, such that f' = Uf . i.

When we say that lists are free monoids, we mean specifically that (the type) [x] is free on (the type) x, via the function i = \x -> [x] (on values). This is free, because given any other monoid z, and function f' :: x -> z, then we can lift to a monoid homomorphism f :: [x] -> z, with f' = f . i. How? Well, the basic idea is to use concatMap. The type of concatMap is:
concatMap :: (a -> [b]) -> [a] -> [b]
So it's doing the lifting we want. However this isn't quite right, because this assumes that the target monoid z is a list. So we need this slight variant:

mconcatmap :: (Monoid z) => (x -> z) -> [x] -> z
mconcatmap f xs = mconcat (map f xs)

If we set f = mconcatmap f', then we will have

(f . i) x
= f (i x)
= f [x]
= mconcatmap f' [x]
= mconcat (map f' [x])
= mconcat [f' x]
= foldr mappend mempty [f' x]  -- definition of mconcat
= mappend mempty (f' x)  -- definition of foldr
= f' x  -- identity of mempty

Now, what would it mean for ([a],[b]) to be free? Well, first, what is it going to be free on? To be free on a and b is the same as being free on Either a b (the disjoint union of a and b). Then our function i is going to be

i (Left a) = ([a],[])
i (Right b) = ([],[b])

Then for ([a],[b]) to be free would mean that whenever we have a function f' :: Either a b -> z, with z a monoid, then we can lift it to a monoid homomorphism f : ([a],[b]) -> z, such that f' = f . i.

So can we?

Well, what if our target monoid z doesn't satisfy the a-b commutativity relations that we saw. That is, what if:
f' a1 `mappend` f' b1 /= f' b1 `mappend` f' a1 -- (A)
That would be a problem.

We are required to find an f such that f' = f . i.
We know that i a1 = ([a1],[]), i b1 = ([],[b1]). So we know that i a1 `mappend` i b1 = i b1 `mappend` i a1.
f is required to be a monoid homomorphism, so by definition:

f (i a1 `mappend` i b1) = f (i a1) `mappend` f (i b1)
f (i b1 `mappend` i a1) = f (i b1) `mappend` f (i a1)

But then since the two left hand sides are equal, then so are the two right hand sides, giving:
f (i a1) `mappend` f (i b1) = f (i b1) `mappend` f (i a1) -- (B)

But now we have a contradiction between (A) and (B), since f' = f . i.

So for a concrete counterexample, showing that ([a],[b]) is not free, all we need is a monoid z in which the a-b commutativity relations don't hold. Well that's easy: [Either a b]. Just take f' :: Either a b -> [Either a b], f' = \x -> [x]. Now try to find an f :: ([a],[b]) -> [Either a b], with f' = f . i.

The obvious f is f (as,bs) = map Left as ++ map Right bs.
But the problem is that this f isn't a monoid homomorphism:
f ( ([],[b1]) `mappend` ([a1],[]) ) /= f ([],[b1]) `mappend` f ([a1],[])

Notice the connection between the two definitions of free. It was because ([a],[b]) had non-trivial relations that we couldn't lift a function to a monoid homomorphism in some cases. The cases where we couldn't were where the target monoid z didn't satisfy the relations.

Okay, so sorry, that got a bit technical. To summarise, the product of [a], [b] in the category of lists / free monoids is [Either a b].

What about vector spaces? What is the product of Vect k a and Vect k b?

Well, similarly to lists, we can make (Vect k a, Vect k b) into a vector space, by defining
0 = (0,0)
(a1,b1) + (a2,b2) = (a1+a2,b1+b2)
k(a,b) = (ka,kb)

Exercise: Show that with these definitions, fst, snd and f .*. g are vector space morphisms (linear maps).

Alternatively, Vect k (Either a b) is of course a vector space. We can define:

p1 = linear p1' where
    p1' (Left a) = return a
    p1' (Right b) = zero

p2 = linear p2' where
    p2' (Left a) = zero
    p2' (Right b) = return b

prodf f g = linear fg' where
    fg' b = fmap Left (f (return b)) <+> fmap Right (g (return b))

In this case p1, p2, f×g are vector space morphisms by definition, since they were constructed using "linear". How do we know that they satisfy the product property? Well, this looks like a job for QuickCheck. The following code builds on the code we developed last time:

prop_Product (f',g',x) =
    f x == (p1 . fg) x &&
    g x == (p2 . fg) x
    where f = linfun f'
          g = linfun g'
          fg = prodf f g

newtype SBasis = S Int deriving (Eq,Ord,Show,Arbitrary)

prop_ProductQn (f,g,x) = prop_Product (f,g,x)
    where types = (f,g,x) :: (LinFun Q SBasis ABasis, LinFun Q SBasis BBasis, Vect Q SBasis)

> quickCheck prop_ProductQn
+++ OK, passed 100 tests.

As we did with lists, we can ask again, which is the correct definition of product, (Vect k a, Vect k b), or Vect k (Either a b)?

Well, in this case it turns out that they are equivalent to one another, via the mutually inverse isomorphisms

\(va,vb) -> fmap Left va <+> fmap Right vb
\v -> (p1 v, p2 v)

Unlike in the list case, these are both vector space morphisms (linear functions).

Why the difference? Why does it work out for vector spaces whereas it didn't for lists? Well, I think it's basically because vector spaces are commutative.

(It is also the case that vector spaces are always free on a basis. So since we have an obvious bijection between the bases of (Vect k a, Vect k b) and Vect k (Either a b), then we must have an isomorphism between the vector spaces.)

Now, we're left with a little puzzle. We have found that both the product and the coproduct of two vector spaces is Vect k (Either a b). So we still haven't figured out what Vect k (a,b) represents.

Coproducts of lists and free vector spaces

2011-01-21T16:35:00.001Z

Recently we've been looking at vector spaces. We defined a type Vect k b, representing the free k-vector space over a type b - meaning, the vector space consisting of k-linear combinations of the inhabitants of b - so b is the basis. Like any good mathematical structure, vector spaces admit various new-from-old constructions. Last time I posed the puzzle, what do Vect k (Either a b) and Vect k (a,b) represent? As we're aiming for quantum algebra, I'm going to frame the answers in the language of category theory.

Suppose we have objects A and B in some category. Then their coproduct (if it exists) is an object A+B, together with injections i1: A -> A+B, i2: B -> A+B, with the following universal property: whenever we have arrows f: A -> T and g: B -> T, they can be factored through A+B to give an arrow f+g: A+B -> T, such that f = f+g . i1, g = f+g . i2.

Notice that this definition does not give us a construction for the coproduct. In any given category, it doesn't tell us how to construct the coproduct, or even if there is one. Even if we have a construction for the coproduct in one category, there is no guarantee that it, or something similar, will work in another related category.

In the category Set, the coproduct of sets A and B is their disjoint union. In order to see this, we can work in the category Hask of Haskell types. We can regard Hask as a subcategory of Set, by identifying a type with its set of inhabitants. If a and b are Haskell types / sets of inhabitants, then their disjoint union is Either a b. The elements of Either a b can be from either a or b (hence, from their union), and they are kept disjoint in the left and right parts of the union (so that for example Either a a contains two copies of a, not just one). The injections i1 and i2 are then the value constructors Left and Right. Given f :: a -> t, g :: b -> t, we define:

(f .+. g) (Left a) = f a
(f .+. g) (Right b) = g b

Then it should be clear that (f .+. g) . Left = f, and (f .+. g) . Right = g, as required

In a moment, we'll look at coproducts of vector spaces, but first, as a warmup, let's think about the coproducts in a simpler category: lists / free monoids. Recall that a monoid is an algebraic structure having an associative operation ++, and an identity for ++ called []. (That is [] ++ x = x = x ++ [].)

A monoid homomorphism is a function f :: [a] -> [b] such that f [] = [] and f (a1 ++ a2) = f a1 ++ f a2. With a little thought, you should be able to convince yourself that all monoid homomorphisms are of the form concatMap f', where f' :: a -> [b]. (Which is, incidentally, the same as saying that they are of the form (>>= f').) In the category of free monoids, the arrows are constrained to be monoid homomorphisms.

So for our coproduct, we are looking for an object satisfying the universal property shown in the following diagram:

Perhaps the first thing to try is the disjoint union: Either [a] [b]. This is the coproduct of [a] and [b] as sets, but is it also their coproduct as monoids? Well, let's see.

Hmm, firstly, it's not a list (doh!): so you can't apply ++ to it, and it doesn't have a []. However, before we give up on that, let's consider whether we're asking the right question. Perhaps we should only be requiring that Either [a] [b] is (or can be made to be) a Monoid instance. Is the disjoint union of two monoids a monoid? Suppose we try to define ++ for it:
Left a1 ++ Left a2 = Left (a1++a2)
Right b1 ++ Right b2 = Right (b1++b2)
But now we begin to see the problem. What are we going to do for Left as ++ Right bs? There's nothing sensible we can do, because our disjoint union Either [a] [b] does not allow mixed lists of as and bs.

However, this immediately suggests that we would be better off looking at [Either a b] - the free monoid over the disjoint union of a and b. This is a list - and it does allow us to form mixed lists of as and bs.

We can then set i1 = map Left, i2 = map Right, and these are list homomorphisms (they interact with [] and ++ in the required way). Then we can define:

h = concatMap h' where
    h' (Left a) = f' a
    h' (Right b) = g' b

So our suspicion is that [Either a b] is the coproduct, with h the required coproduct map.

Let's just check the coproduct conditions:
h . i1
= concatMap h' . map Left
= concatMap f'
= f
and similarly, h . i2 = g, as required.

Notice that the disjoint union Either [a] [b] is (isomorphic to) a subset of the free monoid [Either a b], via Left as -> map Left as; Right bs -> map Right bs. So we were thinking along the right lines in suggesting Either [a] [b]. The problem was that Either [a] [b] isn't a monoid, it's only a set. We can regard [Either a b] as the closure of Either [a] [b] under the monoid operations. [Either a b] is the smallest free monoid containing the disjoint union Either [a] [b] (modulo isomorphism of Haskell types).

(This is a bit hand-wavy. This idea of closure under algebraic operations makes sense in maths / set theory, but I'm not quite sure how best to express it in Haskell / type theory. If anyone has any suggestions, I'd be pleased to hear them.)

Okay, so what about a coproduct in the category of k-vector spaces. First, recall that the arrows in this category are linear maps f satisfying f (a+b) = f a + f b, f (k*a) = k * f a. Again, it should be obvious that a linear map is fully determined by its action on basis elements - so every linear map f :: Vect k a -> Vect k b can be expressed as linear f' where f' :: a -> Vect k b.

Recall that we defined linear f' last time - it's really just (>>= f'), but followed by reduction to normal form:
linear f v = nf $ v >>= f

Okay, so vector spaces of course have underlying sets, so we will expect the coproduct of Vect k a and Vect k b to contain the disjoint union Either (Vect k a) (Vect k b). As with lists though, we will have the problem that this is not closed under vector addition - we can't add an element of Vect k a to an element of Vect k b within this type.

So as before, let's try Vect k (Either a b). Then we can set i1 = fmap Left, i2 = fmap Right, and they are both linear maps by construction. (We don't need to call nf afterwards, since Left and Right are order-preserving.)

i1 = fmap Left
i2 = fmap Right

Then we can define the coproduct map (f+g) as follows:

coprodf f g = linear fg' where
    fg' (Left a) = f (return a)
    fg' (Right b) = g (return b)

We need to verify that this satisfies the coproduct conditions: f+g . i1 = f and f+g . i2 = g. It would be nice to test this using a QuickCheck property. In order to do that, we need a way to construct arbitrary linear maps. (This is not the same thing as arbitrary functions Vect k a -> Vect k b, so I don't think that I can use QuickCheck's Coarbitrary class - but the experts may know better.) Luckily, that is fairly straightforward: we can construct arbitrary lists [(a, Vect k b)], and then each pair (a,vb) can be interpreted as saying that the basis element a is taken to the vector vb.

type LinFun k a b = [(a, Vect k b)]

linfun :: (Eq a, Ord b, Num k) => LinFun k a b -> Vect k a -> Vect k b
linfun avbs = linear f where
    f a = case lookup a avbs of
          Just vb -> vb
          Nothing -> zero

With that preparation, here is a QuickCheck property that expresses the coproduct condition.

prop_Coproduct (f',g',a,b) =
    f a == (fg . i1) a &&
    g b == (fg . i2) b
    where f = linfun f'
          g = linfun g'
          fg = coprodf f g

That property can be used for any vector spaces. Let's define some particular vector spaces to do the test on.

newtype ABasis = A Int deriving (Eq,Ord,Show,Arbitrary) -- GeneralizedNewtypeDeriving
newtype BBasis = B Int deriving (Eq,Ord,Show,Arbitrary)
newtype TBasis = T Int deriving (Eq,Ord,Show,Arbitrary)

instance (Num k, Ord b, Arbitrary k, Arbitrary b) => Arbitrary (Vect k b) where
    arbitrary = do ts <- arbitrary :: Gen [(b, k)] -- ScopedTypeVariables
                   return $ nf $ V ts

(I should emphasize that not all vector space bases are newtypes around Int - we can have finite bases, or bases with other interesting internal structure, as we will see in later installments. For the purposes of this test however, I think this is sufficient.)

prop_CoproductQn (f,g,a,b) = prop_Coproduct (f,g,a,b)
    where types = (f,g,a,b) :: (LinFun Q ABasis TBasis, LinFun Q BBasis TBasis, Vect Q ABasis, Vect Q BBasis)

> quickCheck prop_CoproductQn
+++ OK, passed 100 tests.

So we do indeed have a coproduct on vector spaces. To summarise: The coproduct of free vector spaces is the free vector space on the coproduct (of the bases).

Next time, we'll look at products - where there might be a small surprise.

The free vector space on a type, part 2

2011-01-10T21:26:00.000Z

Last time, I defined the free k-vector space over a type b:

data Vect k b = V [(b,k)]
Elements of Vect k b represent formal sums of scalar multiples of elements of b, where the scalars are taken from the field k. (For example, V [(E1,5),(E3,2)] represents the formal sum 5E1+2E3.) Thus b is the basis for the vector space.

We saw that there is a functor (in the mathematical sense) from the category Hask (of Haskell types and functions) to the category k-Vect (of k-vector spaces and k-linear maps). In maths, we would usually represent a functor by a capital letter, eg F, and apply it to objects and arrows by prefixing. For example, if a, b are objects in the source category, then the image objects in the target category would be called F a and F b. If f :: a -> b is an arrow in the source category, then the image arrow in the target category would be called F f.

Haskell allows us to declare a type constructor as a Functor instance, and give it an implementation of fmap. This corresponds to describing a functor (in the mathematical sense), but with a different naming convention. In our case, we declared the type constructor (Vect k) as a Functor instance. So the functor's action on objects is called (Vect k) - given any object b in Hask (ie a Haskell type), we can apply (Vect k) to get an object Vect k b in k-Vect (ie a k-vector space). However, the functor's action on arrows is called fmap - given any arrow f :: a -> b in Hask (ie a Haskell function), we can apply fmap to get an arrow fmap f :: Vect k a -> Vect k b in k-Vect (ie a k-linear map).

Haskell allows us to declare only type constructors as functors. In maths, there are many functors which are not of this form. For example, the simplest is the forgetful functor. Given any algebraic category A, we have the forgetful functor A -> Set, which simply forgets the algebraic structure. The forgetful functor takes the objects of A to their underlying sets, and the arrows of A to the underlying functions.

For example, in our case, the forgetful functor k-Vect -> Hask consists in forgetting that the objects Vect k b are vector spaces (with addition, scalar multiplication etc defined), and considering them just as Haskell types; and forgetting that the arrows Vect k a -> Vect k b are linear maps, and considering them just as Haskell functions.

(Notice that when working in Haskell, the category Hask acts as a kind of stand-in for the category Set. If we identify a type with its inhabitants (which form a set), then Hask is something like the computable subcategory of Set.)

So we actually have functors going both ways. We have the free functor F: Hask -> k-Vect - which takes types and functions to free k-vector spaces and linear maps. And we have the forgetful functor G: k-Vect -> Hask - which takes free k-vector spaces and linear maps, and just forgets their algebraic structure, so that they're just types and functions again.

Note that these two functors are not mutual inverses. (G . F) is not the identity on Hask - indeed it takes b to Vect k b, both considered as objects in Hask (and similarly with arrows). The two functors are however adjoint. (I'm not going to explain what this means, but see Wikipedia, or most books on category theory.)

Whenever we have an adjunction, then in fact we also have a monad. Here's the definition for our case:

instance Num k => Monad (Vect k) where
    return a = V [(a,1)]
    V ts >>= f = V $ concat [ [(b,y*x) | let V us = f a, (b,y) <- us] | (a,x) <- ts]

This monad is most easily understood using the monad as container analogy. A free k-vector space over b is just a container of elements of b. Okay, so it's a slightly funny sort of container. It most resembles a multiset or bag - an unordered container in which you can have more than one of each element. However, free k-vector spaces go further. A free Q-vector space is a container in which we're allowed to have fractional or negative amounts of each basis element, such as 1/2 e1 - 3 e2. In a free C-vector space, we're even allowed imaginary amounts of each element, such as i e3.

These oddities aside though, free k-vector spaces are monads in much the same way as any other container. For example, let's compare them to the list monad.
- For a container monad, return means "put into the container". For List, return a is [a]. For k-Vect, return a is 1 a (where 1 is the scalar 1 in our field k).
- For container monads, it's most natural to look next at the join operation, which combines a container of containers into a single container. For List, join = concat. For k-Vect, join combines a linear combination of linear combinations into a single linear combination (in the obvious way).
- Finally there is bind or (>>=). bind can be defined in terms of join and fmap:
x >>= f = join ((fmap f) x)
For lists, bind is basically concatMap. For k-Vect, bind corresponds to extending a function on basis elements to a function on vectors "by linearity". That is, if f :: a -> Vect k b, then (>>= f) :: Vect k a -> Vect k b is defined (in effect, by structural induction) so as to be linear, by saying that (>>= f) 0 = 0, (>>= f) (k a) = k (f a), (>>= f) (a + b) = f a + f b.

So k-Vect is like just a strange sort of list. Think of them as distant cousins. Incidentally, this is not only because they're both containers (if you believed my story about k-Vect being a container). There's another reason: Both the List and k-Vect monads arise from free-forgetful adjunctions. In the case of lists, the list datatype is the free monoid, and the list monad arises from the free-forgetful adjunction for monoids.

That's most of what I wanted to say for now. However, there's one small detail to add. Last time we defined a normal form for elements of Vect k b, in which the basis elements are in order, without repeats, and none of the scalars are zero. In order to calculate this normal form, we require an Ord instance for b. Unfortunately, Haskell doesn't let us specify that when defining the Monad instance for (Vect k). So whenever we use (>>=), we should call nf afterwards, to put the result in normal form.

For this reason, we define a convenience function:

linear :: (Ord b, Num k) => (a -> Vect k b) -> Vect k a -> Vect k b
linear f v = nf $ v >>= f

Given f :: a -> Vect k b, linear f :: Vect k a -> Vect k b is the extension of f from basis elements to vectors "by linearity". Hence, linear f is guaranteed to be linear, by construction. We can confirm this on an example using the QuickCheck properties we defined last time.

> let f (E 1) = e1; f (E 2) = e1 <+> e2; f _ = zero
> (linear f) (e1 <+> e2)
2e1+e2
> quickCheck (prop_LinearQn (linear f))
+++ OK, passed 100 tests.

Acknowledgement: I'm partially retreading in Dan Piponi's steps here. He first described the free vector space monad here. When I come to discuss the connection between quantum algebra and knot theory, I'll be revisiting some more material that Dan sketched out here.

Exercise for next time: What does Vect k (Either a b) represent? What does Vect k (a,b) represent?

The free vector space on a type, part 1

2010-12-13T19:53:00.000Z

As I mentioned last time, I want to spend the next few posts talking about quantum algebra. Well, we've got to start somewhere, so let's start with vector spaces.

You probably know what a vector space is. It's what is sounds like: a space of vectors, that you can add together, or multiply by scalars (that is, real numbers, or more generally, elements of some field k). Here's the official definition:
An additive (or Abelian) group is a set with a binary operation called addition, such that
- addition is associative: x+(y+z) = (x+y)+z
- addition is commutative: x+y = y+x
- there is an additive identity: x+0 = x = 0+x
- there are additive inverses: x+(-x) = 0 = (-x)+x
A vector space over a field k is an additive group V, together with an operation k * V -> V called scalar multiplication, such that:
- scalar multiplication distributes over vector addition: a(x+y) = ax+ay
- scalar multiplication distributes over scalar addition: (a+b)x = ax+bx
- associativity: (ab)x = a(bx)
- unit: 1a = a

There are some obvious examples:
- R^2 is a 2-dimensional vector space over the reals R, R^3 is a 3-dimensional vector space. (I'm not going to define dimension quite yet, but hopefully it's intuitively obvious what it means.)
- R^n is an R-vector space for any n.
- Indeed, k^n is a k-vector space for any field k.

Some slightly more interesting examples:
- C is a 2-dimensional vector space over R. (The reason it's more interesting is that of course C possesses additional algebraic structure, beyond the vector space structure.)
- 2*2 matrices over k form a 4-dimensional k-vector space.
- Polynomials in X with coefficients in k form an (infinite dimensional) k-vector space.

If we wanted to code the above definition into Haskell, probably the first idea that would come to mind would be to use type classes:

class AddGrp a where
    add :: a -> a -> a
    zero :: a
    neg :: a -> a -- additive inverse

class (Field k, AddGrp v) => VecSp k v where
    smult :: k -> v -> v -- scalar multiplication

(Type classes similar to these are defined are in the vector-space package.)

For most vector spaces that one encounters in "real life", there is some set of elements, usually obvious, which form a "basis" for the vector space, meaning that all elements can be expressed as linear combinations of basis elements. For example, in R^3, the obvious basis is {(0,0,1), (0,1,0), (1,0,0)}. Any element (x,y,z) of R^3 can be expressed as the linear combination x(1,0,0)+y(0,1,0)+z(0,0,1).

(Mathematicians would want to stress that there are other bases for R^3 that would serve equally well, and indeed, that a significant part of the theory of vector spaces can go through without even talking about bases. However, for our purposes - we want to write code to calculate in vector spaces - then working with a basis is natural.)

Okay, so we want a way to build a vector space from a basis. (More specifically, a k-vector space, for some given field k.) What sorts of things shall we allow as our basis? Well, why not just allow any type, whatsoever:

module Math.Algebras.VectorSpace where

import qualified Data.List as L

data Vect k b = V [(b,k)] deriving (Eq,Ord)

This says that a k-vector space over basis b consists of a linear combination of elements of b. (So the [(b,k)] is to be thought of as a sum, with each (b,k) pair representing a basis element in b multiplied by a scalar coefficient in k.)

For example, we can define the "boring basis" type, which just consists of numbered basis elements:

newtype EBasis = E Int deriving (Eq,Ord)

instance Show EBasis where show (E i) = "e" ++ show i

e i = return (E i) -- don't worry about what "return" is doing here for the moment
e1 = e 1
e2 = e 2
e3 = e 3

Then a typical element of Vect Double EBasis is:

> :load Math.Algebras.VectorSpace
> V [(E 1, 0.5), (E 3, 0.7)]
0.5e1+0.7e3

So of course, the Show instances for EBasis (see above), and Vect k b (not shown) are coming into play here.

How do we know that this is a vector space? Well actually, it's not yet, because we haven't defined the addition and scalar multiplication operations on it. So, without further ado:

infixr 7 *>
infixl 6 <+>

-- |The zero vector
zero :: Vect k b
zero = V []

-- |Addition of vectors
add :: (Ord b, Num k) => Vect k b -> Vect k b -> Vect k b
add (V ts) (V us) = V $ addmerge ts us

-- |Addition of vectors (same as add)
(<+>) :: (Ord b, Num k) => Vect k b -> Vect k b -> Vect k b
(<+>) = add

addmerge ((a,x):ts) ((b,y):us) =
    case compare a b of
    LT -> (a,x) : addmerge ts ((b,y):us)
    EQ -> if x+y == 0 then addmerge ts us else (a,x+y) : addmerge ts us
    GT -> (b,y) : addmerge ((a,x):ts) us
addmerge ts [] = ts
addmerge [] us = us

-- |Negation of vector
neg :: (Num k) => Vect k b -> Vect k b
neg (V ts) = V $ map (\(b,x) -> (b,-x)) ts

-- |Scalar multiplication (on the left)
smultL :: (Num k) => k -> Vect k b -> Vect k b
smultL 0 _ = zero -- V []
smultL k (V ts) = V [(ei,k*xi) | (ei,xi) <- ts]

-- |Same as smultL. Mnemonic is "multiply through (from the left)"
(*>) :: (Num k) => k -> Vect k b -> Vect k b
(*>) = smultL

A few things to mention:
- First, note that we required a Num instance for k. Strictly speaking, as we stated that k is a field, then we should have required a Fractional instance. However, on occasion we are going to break the rules slightly.
- Second, note that for addition, we required an Ord instance for b. We could have defined addition using (++) to concatenate linear combinations - however, the problem with that is that it wouldn't then easily follow that e1+e3 = e3+e1, or that e1+e1 = 2e1. By requiring an Ord instance, we can guarantee that there is a unique normal form in which to express any vector - namely, list the basis elements in order, combine duplicates, remove zero coefficients.
- Finally, note that I didn't define Vect k b as an instance of a vector space type class. That's just because I didn't yet see a reason to.

It will turn out to be useful to have a function that can take an arbitrary element of Vect k b and return a vector in normal form:

-- |Convert an element of Vect k b into normal form. Normal form consists in having the basis elements in ascending order,
-- with no duplicates, and all coefficients non-zero
nf :: (Ord b, Num k) => Vect k b -> Vect k b
nf (V ts) = V $ nf' $ L.sortBy compareFst ts where
    nf' ((b1,x1):(b2,x2):ts) =
        case compare b1 b2 of
        LT -> if x1 == 0 then nf' ((b2,x2):ts) else (b1,x1) : nf' ((b2,x2):ts)
        EQ -> if x1+x2 == 0 then nf' ts else nf' ((b1,x1+x2):ts)
        GT -> error "nf': not pre-sorted"
    nf' [(b,x)] = if x == 0 then [] else [(b,x)]
    nf' [] = []
    compareFst (b1,x1) (b2,x2) = compare b1 b2

Okay, so we ought to check that the Vect k b is a vector space. Let's write some QuickCheck properties:

prop_AddGrp (x,y,z) =
    x <+> (y <+> z) == (x <+> y) <+> z && -- associativity
    x <+> y == y <+> x                 && -- commutativity
    x <+> zero == x                    && -- identity
    x <+> neg x == zero                   -- inverse

prop_VecSp (a,b,x,y,z) =
    prop_AddGrp (x,y,z) &&
    a *> (x <+> y) == a *> x <+> a *> y && -- distributivity through vectors
    (a+b) *> x == a *> x <+> b *> x     && -- distributivity through scalars
    (a*b) *> x == a *> (b *> x)         && -- associativity
    1 *> x == x                            -- unit

instance Arbitrary EBasis where
    arbitrary = do n <- arbitrary :: Gen Int
                   return (E n)

instance Arbitrary Q where
    arbitrary = do n <- arbitrary :: Gen Integer
                   d <- arbitrary :: Gen Integer
                   return (if d == 0 then fromInteger n else fromInteger n / fromInteger d)

instance Arbitrary (Vect Q EBasis) where
    arbitrary = do ts <- arbitrary :: Gen [(EBasis, Q)]
                   return $ nf $ V ts

prop_VecSpQn (a,b,x,y,z) = prop_VecSp (a,b,x,y,z)
    where types = (a,b,x,y,z) :: (Q, Q, Vect Q EBasis, Vect Q EBasis, Vect Q EBasis)

> quickCheck prop_VecSpQn
+++ OK, passed 100 tests.

(I'm using Q instead of R as my field in order to avoid false negatives caused by the fact that arithmetic in Double is not exact.)

So it looks like Vect k b is indeed a vector space. In category theory, it is called the free k-vector space over b. "Free" here means that there are no relations among the basis elements: it will never turn out, for example, that e1 = e2+e3.

Vector spaces of course form a category, specifically an algebraic category (there are other types, as we'll see in due course). The objects in the category are the vector spaces. The arrows or morphisms in the category are the functions between vector spaces which "commute" with the algebraic structure. Specifically, they are the functions f such that:
- f(x+y) = f(x)+f(y)
- f(0) = 0
- f(-x) = -f(x)
- f(a.x) = a.f(x)

Such a function is called linear, the idea being that it preserves lines. This is because it follows from the the conditions that f(a.x+b.y) = a.f(x)+b.f(y) .

We can write a QuickCheck property to check whether a given function is linear:

prop_Linear f (a,x,y) =
    f (x <+> y) == f x <+> f y &&
    f zero == zero             &&
    f (neg x) == neg (f x)     &&
    f (a *> x) == a *> f x

prop_LinearQn f (a,x,y) = prop_Linear f (a,x,y)
    where types = (a,x,y) :: (Q, Vect Q EBasis, Vect Q EBasis)

For example:

> quickCheck (prop_LinearQn (2 *>))
+++ OK, passed 100 tests.

We won't need to use this quite yet, but it's handy to have around.

Now, in category theory we have the concept of a functor, which is a map from one category to another, which commutes with the category structure. Specifically, a functor F consists of a map from the objects of one category to the objects of the other, and from the arrows of one category to the arrows of the other, satisfying:
- F(id_A) = id_F(A)
- F(f . g) = F(f) . F(g) (where dot denotes function composition)

How does this relate to the Functor type class in Haskell? Well, the Haskell type class enables us to declare that a /type constructor/ is a functor. For example, (Vect k) is a type constructor, which acts on a type b to construct another type Vect k b. (Vect k) is indeed a functor, witness the following declaration:

instance Functor (Vect k) where
    fmap f (V ts) = V [(f b, x) | (b,x) <- ts]

This says that if we have a function f on our basis elements, then we can lift it to a function on linear combinations of basis elements in the obvious way.

In mathematics, we would think of the free vector space construction as a functor from Set (the category of sets) to k-Vect (the category of k-vector spaces). In Haskell, we need to think of the (Vect k) construction slightly differently. It operates on types, rather than sets, so the source category is Hask, the category of Haskell types.

What is the relationship between Hask and Set? Well, if we identify a type with the set of values which inhabit it, then we can regard Hask as a subcategory of Set, consisting of those sets and functions which can be represented in Haskell. (That would imply for example that we are restricted to computable functions.)

So (Vect k) is a functor from Hask to the subcategory of k-Vect consisting of vector spaces over sets/types in Hask.

So just to spell it out, the (Vect k) functor:
- Takes an object b in Hask/Set - ie a type, or its set of inhabitants - to an object Vect k b in k-Vect
- Takes an arrow f in Hask (ie a function f :: a -> b), to an arrow (fmap f) :: Vect k a -> Vect k b in k-Vect.

Now, there's just one small fly in the ointment. In order to get equality of vectors to work out right, we wanted to insist that they were expressed in normal form, which meant we needed an Ord instance for b. However, in the Functor instance for (Vect k), Haskell doesn't let us express this constraint, and our fmap is unable to use the Ord instance for b. What this means is that fmap f might return a vector which is not in normal form - so we need to remember to call nf afterwards. For example:

newtype FBasis = F Int deriving (Eq,Ord)

instance Show FBasis where show (F i) = "f" ++ show i

> let f = \(E i) -> F (10 - div i 2)
> let f' = fmap f :: Vect Q EBasis -> Vect Q FBasis
> f' (e1 <+> 2 *> e2 <+> e3)
f10+2f9+f9
> let f'' = nf . fmap f :: Vect Q EBasis -> Vect Q FBasis
> f'' (e1 <+> 2 *> e2 <+> e3)
3f9+f10

So it might be fairer to say that it is the combination of nf and fmap that forms the functor on arrows.

The definition of a functor requires that the target arrow is an arrow in the target category. In this case, the requirement is that it is a linear function, rather than just any function between vector spaces. So let's just check:

> quickCheck (prop_LinearQn f'')
+++ OK, passed 100 tests.

That's enough for now - more next time.

New modules - Quantum Algebra

2010-11-07T18:16:00.000Z

I've put up a new version of HaskellForMaths on Hackage, v0.3.1. It's quite a significant update, with more than a dozen new modules, plus improved documentation of several existing modules. I wrote the new modules in the course of reading Kassel's Quantum Groups. The modules are about algebras, coalgebras, bialgebras, Hopf algebras, tensor categories and quantum algebra.

The new modules fall into two groups:

Math.Algebras.* - Modules about algebras (and co-, bi- and Hopf algebras) in general
Math.QuantumAlgebra.* - Modules specifically about quantum algebra

In (slightly) more detail, here are the modules:

Math.Algebras.VectorSpace - defines a type for the free k-vector space over a basis set b
Math.Algebras.TensorProduct - defines tensor product of two vector spaces
Math.Algebras.Structures - defines a number of additional algebraic structures that can be given to vector spaces: algebra, coalgebra, bialgebra, Hopf algebra, module, comodule
Math.Algebras.Quaternions - a simple example of an algebra
Math.Algebras.Matrix - the 2*2 matrices - another simple example of an algebra
Math.Algebras.Commutative - commutative polynomials (such as x^2+3yz) - another algebra
Math.Algebras.NonCommutative - non-commutative polynomials (where xy /= yx) - another algebra
Math.Algebras.GroupAlgebra - a key example of a Hopf algebra
Math.Algebras.AffinePlane - the affine plane and its symmetries - more Hopf algebras, preparing for the quantum plane
Math.Algebras.TensorAlgebra
Math.Algebras.LaurentPoly - we use Laurent polynomials in q (that is, polynomials in q and q^-1) as our quantum scalars

Math.QuantumAlgebra.QuantumPlane - the quantum plane and its symmetries, as examples of non-commutative, non-cocommutative Hopf algebras
Math.QuantumAlgebra.TensorCategory
Math.QuantumAlgebra.Tangle - The tangle category (which includes knots, links and braids as subcategories), and some representations (from which we derive knot invariants)
Math.QuantumAlgebra.OrientedTangle

The following diagram is something like a dependency diagram, with "above" meaning "depends on".

Each layer depends on the layers beneath it. There are also a few dependencies within layers, that I've hinted at by proximity.

Some of these modules overlap somewhat in content with other modules that already exist in HaskellForMaths. In particular, there are already modules for commutative algebra, non-commutative algebra, and knot theory. For the moment, those existing modules still offer some features not offered by the new modules (for example, calculation of Groebner bases).

For the next little while in this blog, I want to start going through these new modules, and investigating quantum algebra. I should emphasize that this is still work in progress. For example, I'm intending to add modules for quantum enveloping algebras - but I have some reading to do first.

(Oh, and I know that I did previously promise to look at finite simple groups, and Coxeter groups, in this blog. I'll probably still come back to those at some point.)

Word length in the Symmetric group

2010-10-14T20:03:00.000+01:00

Previously on this blog, we saw how to think about groups abstractly via group presentations, where a group is given as a set of generators satisfying specified relations. Last time, we saw that questions about the length of reduced words in such a presentation can be visualised as questions about the length of paths in the Cayley graph of the group (relative to the generators).

This time, I want to focus on just one family of groups - the symmetric groups Sn, as generated by the adjacent transpositions {si = (i i+1)}. Here's the Haskell code defining this presentation of Sn:

newtype SGen = S Int deriving (Eq,Ord)

instance Show SGen where
    show (S i) = "s" ++ show i

_S n = (gs, r ++ s ++ t) where
    gs = map S [1..n-1]
    r = [([S i, S i],[]) | i <- [1..n-1]]
    s = [(concat $ replicate 3 [S i, S (i+1)],[]) | i <- [1..n-2]]
    t = [([S i, S j, S i, S j],[]) | i <- [1..n-1], j <- [i+2..n-1]]

The three sets of relations say: each generator si squares to the identity; if i, j are not adjacent, then si and sj commute; if i, j are adjacent, then (si*sj)^3 is the identity.

Here is the Cayley graph for S4 under this presentation:

The vertices are labeled with the reduced words in the generators si. How can we find out which permutations these correspond to? Well, that's easy:

fromTranspositions ts = product $ map (\(S i) -> p [[i,i+1]]) ts
For example:

> :load Math.Algebra.Group.CayleyGraph
> fromTranspositions [S 1, S 2, S 1, S 3, S 2, S 1]
[[1,4],[2,3]]

This is the permutation that reverses the list [1..4]

> map (.^ it) [1..4]
[4,3,2,1]

What about the other way round? Suppose we are given a permutation in Sn. How do we find its expression as a product of the transpositions si? Well the answer is (roughly): use bubblesort!

Here's bubblesort in Haskell:

bubblesort [] = []
bubblesort xs = bubblesort' [] xs where
    bubblesort' ls (r1:r2:rs) = if r1 <= r2 then bubblesort' (r1:ls) (r2:rs) else bubblesort' (r2:ls) (r1:rs)
    bubblesort' ls [r] = bubblesort (reverse ls) ++ [r]

So we sweep through the list from front to back, swapping any pairs that we find out of order - and then repeat. At the end of each sweep, we're guaranteed that the last element (in sort order) has reached the end of the list - so for the next sweep, we can leave it at the end and only sweep through the earlier elements. Hence the list we're sweeping through is one element shorter each time, so we're guaranteed to terminate. (We could terminate early, the first time a sweep makes no swaps - but I haven't coded that.)

Just to prove that the code works:

> bubblesort [2,3,1]
[1,2,3]

How does this help for turning a permutation into a sequence of transpositions? Well it's simple - every time we swap two elements, we are performing a transposition - so just record which swaps we perform. So here's a modified version of the above code, which records the swaps:

-- given a permutation of [1..n] (as a list), return the transpositions which led to it
toTrans [] = []
toTrans xs = toTrans' 1 [] [] xs where
    toTrans' i ts ls (r1:r2:rs) =
        if r1 <= r2
        then toTrans' (i+1) ts (r1:ls) (r2:rs)         -- no swap needed
        else toTrans' (i+1) (S i : ts) (r2:ls) (r1:rs) -- swap needed
    toTrans' i ts ls [r] = toTrans (reverse ls) ++ ts

Notice that the ts are returned in reverse to the order that they were used. This is because we are using them to undo the permutation - so we are performing the inverse of the permutation we are trying to find. Since each generator is its own inverse, we can recover the permutation we are after simply by reversing. In the code, we reverse as we go along.

For example:

> toTrans [2,3,1]
[s1,s2]
> toTrans [4,3,2,1]
[s1,s2,s1,s3,s2,s1]

Now, there's only one problem. As you can see, this code takes as input a rearrangement of [1..n]. This is a permutation, yes, but considered passively. Whereas in this blog we have been more accustomed to thinking of permutations actively, as something a bit like a function, which has an action on a graph, or other combinatorial structure, or if you like, just on the set [1..n]. In other words, our Permutation type represents the action itself, not the outcome of the action. (Recall that the implementation uses a Data.Map of (from,to) pairs.)

But of course it's easy to convert from one viewpoint to the other. So here's the code to take a permutation in cycle notation and turn it into a sequence of transpositions:

-- given a permutation action on [1..n], factor it into transpositions
toTranspositions 1 = []
toTranspositions g = toTrans [i .^ (g^-1) | i <- [1..n] ] where
    n = maximum $ supp g

For example:

> toTranspositions $ p [[1,4],[2,3]]
[s1,s2,s1,s3,s2,s1]

Why does the code have [i .^ (g^-1) | i <- [1..n]], rather than [i .^ g | i <- [1..n]]?
Well, suppose i .^ g = j. This says that g moves i to the j position. But we want to know what ends up in the i position. Suppose that j .^ g = i, for some j. Applying g^-1 to both sides, we see that j = i .^ (g^-1).

Okay, so given a permutation, in either form, we can reconstruct it as a reduced word in the generators.

We saw last time that the length of a reduced word is also the length of the shortest path from 1 to the element in the Cayley graph. Distance in the Cayley graph is a metric on the group, so the length of a reduced word tells us "how far" the element is from being the identity.

If it's only this distance that we're interested in, then there is a more direct way to work it out. Given a permutation g of [1..n], then an inversion is a pair (i,j) with i < j but i .^ g > j .^ g. In Haskell:

inversions g = [(i,j) | i <- [1..n], j <- [i+1..n], i < j, i .^ g > j .^ g]
    where n = maximum $ supp g

For example:

> inversions $ fromList [1,4,3,2]
[(2,3),(2,4),(3,4)]

With a little thought, you should be able to convince yourself that the number of inversions is equal to the length of the reduced word for g - because each swap that we perform during bubblesort corrects exactly one inversion.

Okay, so this is all very nice, but what use is it? Well, of course, maths doesn't have to be useful, any more than any other aesthetic pursuit. However, as it happens, in this case it is.

In statistics, the Kendall tau test gives an indicator of the correlation between two measured quantities (for example, the height and weight of the test subjects). Suppose that we are given a list of pairs (eg (height,weight) pairs), and we want to know how strongly correlated the first and second quantities are.

Ok, so what we do is, we rank the first quantities from lowest to highest, and replace each quantity by its rank (a number from 1 to n). We do the same for the second quantities. So we end up with a list of pairs of numbers from 1 to n. Now, we sort the list on the first element, and then count the number of inversions in the second element.

For example, suppose our original list was [(1.55m, 60kg), (1.8m, 80kg), (1.5m, 70kg), (1.6m, 72kg)]. Converting to ranks, we get [(2nd,1st),(4th,4th),(1st,2nd),(3rd,3rd)]. Sorting on fst, we get [(1,2),(2,1),(3,3),(4,4)]. Looking at snd, we see that we have just one inversion. The idea is that the fewer inversions we have, the better correlated the two quantities. (Of course in reality there's a bit more to it than that - to convert the number of inversions into a probability, we need to know the distribution of word lengths for Sn, where n is the number of pairs of test data that we have.)

So you can think of the Kendall tau test as saying: What permutation has been applied in moving from the first quantities (the heights) to the second quantities (the weights)? How far is that permutation from the identity (on the Cayley graph)? What proportion of all permutations (in Sn) lie at that distance or less from the identity? (Even more concretely, we can imagine colouring in successive shells of the Cayley graph, working out from the identity, until we hit the given permutation, and then asking what proportion of the "surface" of the graph has been coloured.)

Cayley graphs of groups

2010-09-20T21:06:00.000+01:00

[New version HaskellForMaths 0.2.2 released here]

Recently, we've been looking at group presentations, where a group is presented as a set of generators together with a set of relations that hold between those generators. Group elements are then represented as words in the generators.

One can then ask questions about these words, such as: What is the longest (reduced) word in the group? How many (reduced) words are there of each length?

This week I want to look at Cayley graphs, which are a way of visualising groups. Questions about word length translate to questions about path distance in the Cayley graph.

So, the Cayley graph of a group, relative to a generating set gs, is the graph
- with a vertex for each element of the group
- with an edge from x to y just whenever x*g = y for some generator g in gs

Notice that as we have defined it, the edges are directed (from x to y), so this is a directed graph, or digraph.

Here's the Haskell code:

module Math.Algebra.Group.CayleyGraph where

import Math.Algebra.Group.PermutationGroup as P
import Math.Algebra.Group.StringRewriting as SR
import Math.Combinatorics.Graph

import qualified Data.List as L
import qualified Data.Set as S


data Digraph a = DG [a] [(a,a)] deriving (Eq,Ord,Show)

-- Cayley digraph given a group presentation of generators and relations
cayleyDigraphS (gs,rs) = DG vs es where
    rs' = knuthBendix rs
    vs = L.sort $ nfs (gs,rs') -- reduced words
    es = [(v,v') | v <- vs, v' <- nbrs v ]
    nbrs v = L.sort [rewrite rs' (v ++ [g]) | g <- gs]

-- Cayley digraph given group generators as permutations
cayleyDigraphP gs = DG vs es where
    vs = P.elts gs
    es = [(v,v') | v <- vs, v' <- nbrs v ]
    nbrs v = L.sort [v * g | g <- gs]

As an example, let's look at the Cayley digraph of the dihedral group D8 (the symmetries of a square), generated by a rotation and a reflection:

> :load Math.Algebra.Group.CayleyGraph
> let a = p [[1,2,3,4]]
> let b = p [[1,2],[3,4]]
> a^3*b == b*a
True
> cayleyDigraphS (['a','b'],[("aaaa",""),("bb",""),("aaab","ba")])
DG ["","a","aa","aaa","aab","ab","b","ba"] [("","a"),("","b"),("a","aa"),("a","ab"),("aa","aaa"),("aa","aab"),("aaa",""),("aaa","ba"),("aab","aa"),("aab","ab"),("ab","a"),("ab","b"),("b",""),("b","ba"),("ba","aaa"),("ba","aab")]
> cayleyDigraphP [a,b]
DG [[],[[1,2],[3,4]],[[1,2,3,4]],[[1,3],[2,4]],[[1,3]],[[1,4,3,2]],[[1,4],[2,3]],[[2,4]]] [([],[[1,2],[3,4]]),([],[[1,2,3,4]]),([[1,2],[3,4]],[]),([[1,2],[3,4]],[[1,3]]),([[1,2,3,4]],[[1,3],[2,4]]),([[1,2,3,4]],[[2,4]]),([[1,3],[2,4]],[[1,4,3,2]]),([[1,3],[2,4]],[[1,4],[2,3]]),([[1,3]],[[1,4,3,2]]),([[1,3]],[[1,4],[2,3]]),([[1,4,3,2]],[]),([[1,4,3,2]],[[1,3]]),([[1,4],[2,3]],[[1,3],[2,4]]),([[1,4],[2,3]],[[2,4]]),([[2,4]],[[1,2],[3,4]]),([[2,4]],[[1,2,3,4]])]

These are of course the same Cayley digraph, just with different vertex labels. Here's a picture:

The picture might remind you of something. You can think of a Cayley digraph as a state transition diagram, where the states are the group elements, and the transitions are multiplication (on the right) by g, for each generator g. It might help to think of each edge as being labelled by the generator that caused it.

A few things to notice.

First, Cayley digraphs are always regular: the out-degree of each vertex, the number of edges leading out of it, will always equal the number of generators; and similarly for the in-degree, the number of edges leading into each vertex. (Exercise: Prove this.) In fact, we can say more - the graph "looks the same" from any vertex - this follows from the group properties. (Exercise: Explain.)

Second, notice how some of the edges come in pairs going in opposite directions. Why is that? In this case, it's because one of our generators is its own inverse (which one?) - so if it can take you from x to y, then it can take you back again. In general, whenever our set of generators contains a g such that g^-1 is also in the set, then the edges corresponding to g, g^-1 will come in opposing pairs.

Given this, we can omit the arrows on the edges if we adopt the convention that whenever we are given a set of generators, their inverses are also implied. In this way, we obtain an undirected or simple graph. Here's the code:

toSet = S.toList . S.fromList

-- The Cayley graph on the generators gs *and their inverses*, given relations rs
cayleyGraphS (gs,rs) = graph (vs,es) where
    rs' = knuthBendix rs
    vs = L.sort $ nfs (gs,rs') -- all reduced words
    es = toSet [ L.sort [v,v'] | v <- vs, v' <- nbrs v ] -- toSet orders and removes duplicates
    nbrs v = [rewrite rs' (v ++ [g]) | g <- gs]

cayleyGraphP gs = graph (vs,es) where
    vs = P.elts gs
    es = toSet [ L.sort [v,v'] | v <- vs, v' <- nbrs v ]
    nbrs v = [v * g | g <- gs]

For example:

> cayleyGraphS (['a','b'],[("aaaa",""),("bb",""),("aaab","ba")])
G ["","a","aa","aaa","aab","ab","b","ba"] [["","a"],["","aaa"],["","b"],["a","aa"],["a","ab"],["aa","aaa"],["aa","aab"],["aaa","ba"],["aab","ab"],["aab","ba"],["ab","b"],["b","ba"]]

One important point to note is that the Cayley graph of a group is relative to the generators. For example, we saw last time that the dihedral groups can also be generated by two reflections. In the case of D8, we can set r = (1 2)(3 4), s = (1 3).

Before scrolling down, see if you can guess what the Cayley graph looks like. I'll give you a clue: Cayley graphs are always regular - what is the valency of each vertex in this case?

> let r = p [[1,2],[3,4]]
> let s = p [[1,3]]
> (r*s)^4
[]
> cayleyGraphS (['r','s'],[("rr",""),("ss",""),("rsrsrsrs","")])
G ["","r","rs","rsr","rsrs","s","sr","srs"] [["","r"],["","s"],["r","rs"],["rs","rsr"],["rsr","rsrs"],["rsrs","srs"],["s","sr"],["sr","srs"]]

So the point to emphasize is that this graph and the one shown previously are Cayley graphs of the same group. The vertices represent the same elements (considered as permutations). However, because we have taken different sets of generators, we get different edges, and hence different graphs.

Ok, so what does the Cayley graph tell us about the group? Well, as an example, consider the Cayley graph of the Rubik cube group, as generated by the face rotations f, b, l, r, u, d (as defined previously). The vertices of the graph can be identified with the possible positions or states of the cube. The group element 1 corresponds to the solved cube. The edges correspond to single moves that can be made on the cube. If someone gives you a scrambled cube to solve, they are asking you to find a path from that vertex of the Cayley graph back to 1.

Given any graph, and vertices x and y, the distance from x to y is defined as the length of the shortest path from x to y. On the Rubik graph (ie, the Cayley graph of the Rubik cube), the distance from x to 1 is the minimum number of moves needed to solve position x. The HaskellForMaths library provides a distance function on graphs. Thus for example:

> let graphD8rs = cayleyGraphS (['r','s'],[("rr",""),("ss",""),("rsrsrsrs","")])
> distance graphD8rs "" "rsr"
3

The distance from 1 to an element g is of course the length of the reduced word for g.

The diameter of a graph is defined as the maximum distance between vertices. So the diameter of the Rubik graph is the maximum number of moves that are required to solve a scrambled position. It has recently been shown that this number is twenty.

> diameter graphD8rs
4

The diameter of a Cayley graph is the length of the longest reduced word.

That's really all I wanted to say for the moment. Next time, we'll take some of these ideas further.

Group presentations

2010-07-18T20:56:00.001+01:00

Last time we looked at string rewriting systems and the Knuth-Bendix completion algorithm. The motivation for doing that was to enable us to think about groups in a more abstract way than before.

The example we looked at last time was the symmetry group of the square. We found that this group could be generated by two elements a, b, satisfying the relations:
a^4 = 1
b^2 = 1
a^3 b = b a

This way of thinking about the group is called a group presentation, and the usual notation would be:
<a,b | a^4=1, b^2=1, a^3 b = b a>

In our Haskell code, we represent this as:

( ['a','b'], [("aaaa",""),("bb",""),("aaab","ba")] )

We saw how to use the Knuth-Bendix algorithm to turn the relations into a confluent rewrite system:

> :load Math.Algebra.Group.StringRewriting
> mapM_ print $ knuthBendix [("aaaa",""),("bb",""),("aaab","ba")]
("bb","")
("bab","aaa")
("baa","aab")
("aba","b")
("aaab","ba")
("aaaa","")

The rewrite system itself isn't particularly informative. Its importance lies in what it enables us to do. Given any word in the generators, we reduce it as follows: wherever we can find the left hand side of one of rules as a subword, we replace it by the right hand side of the rule. If we keep doing this until there are no more matches, then we end up with a normal form for the element - that is, another word in the generators, which represents the same group element, and is the smallest such word relative to the shortlex ordering. Several things follow from this.

First, the ability to find the shortest word is sometimes useful in itself. If we could do this for the Rubik cube group (taking the six face rotations as generators), then we would be able to code "God's algorithm" to find the shortest solution to any given cube position.

Second, any two words that represent the same group element will reduce to the same normal form. Hence, given any two words in the generators, we can tell whether they represent the same element. This is called "solving the word problem" for the group.

Third, this enables us to list (the normal forms of) all the elements of the group - and hence, among other things, to count them.

Fourth, it enables us to do arithmetic in the group:
- To multiply two elements, represented as words w1, w2, just concatenate them to w1++w2, then reduce using the rewrite system
- The identity element of the group is of course the empty word ""
- But what about inverses?

Strings (lists) under concatenation form a monoid, not a group. So what do we do about inverses?

Well, one possibility is to include them as additional symbols. So, suppose that our generators are a,b. Then we should introduce additional symbols a^-1, b^-1, and consider words over the four symbols {a,b,a^-1,b^-1}. (For brevity, it is customary to use the symbols A, B for a^-1, b^-1.)

If we take this approach, then we will need to add some new rules too. We will need the rules a a^-1 = 1, etc. We will probably also need the "inverses" of the relations in our presentation. For example, if we have a^4 = 1, then we should also have a rule (a^-1)^4 = 1.

It's going to be a bit of a pain. (And it's probably going to cause Knuth-Bendix to get indigestion, in some cases at least.) Luckily, for finite groups, we don't really need this. In a finite group, each generator must have finite order: in our example a^4 = 1, b^2 = 1. So the inverse of each generator is itself a power of that generator - a^-1 = a^3, b^-1 = b. So for a finite group - or in fact any group where the generators are all of finite order - then the inverses are already there, expressible as words in the generators.

So for most purposes, we have no need to introduce the inverses as new symbols. For example, if we want to list the elements of a finite group, or tell whether two words in the generators represent the same element, then we are fine. When it will matter is if we are specifically interested in the length of the words. For example, if we want God's algorithm for solving Rubik's cube, we are interested in the length of words in the generators and their inverses - the clockwise and anti-clockwise face rotations.

There is one situation when even this won't matter - and that is if the generators are their own inverses. If we have a generator g such that g^2 = 1, then it follows that g^-1 = g. Such an element is called an involution.

Are there groups which can be generated by involutions alone? Yes, there are. Let's have a look at a couple.

Consider the symmetry group of a regular polygon, say a pentagon. Consider the two reflections shown below. 'a' is the reflection in an axis through a vertex, and 'b' is the reflection in an axis through the midpoint of an adjacent edge. Hence the angle between the axes is pi/5 (or for an n-gon, pi/n).

It should be clear that ab is a 1/5 rotation of the pentagon. It follows that a,b generate the symmetry group of the pentagon, with a^2 = b^2 = (ab)^5 =1.

> elts (['a','b'], [("aa",""), ("bb",""), ("ababababab","")])
["","a","b","ab","ba","aba","bab","abab","baba","ababa"]
> length it
10

Next, consider the symmetric group S4 (or in general, Sn). It can be generated by the transpositions s1 = (1 2), s2 = (2 3), and s3 = (3 4), which correspond to the diagrams below:

Now, multiplication in the group corresponds to concatenation of the diagrams, going down the page. For example:

In fact, each of those diagrams represents the identity element - as you can check by following each point along the lines down the page, and checking that it ends up where it started. Hence each diagram represents a relation for S4. The diagrams show that s1^2 = 1, (s1s2)^3=1, and (s1s3)^2 = 1.

In the general case, it's clear that Sn can be generated by n-1 transpositions s_i of the form (i i+1), and that they satisfy the following relations:
si^2 = 1
(si sj)^3 = 1 if |i-j| = 1
(si sj)^2 = 1 if |i-j| > 1

Here's some Haskell code to construct these presentations of Sn. (Did I mention that all of the string rewriting code works on arbitrary lists, not just strings?)

newtype S = S Int deriving (Eq,Ord)

instance Show S where
    show (S i) = "s" ++ show i

_S n = (gs, r ++ s ++ t) where
    gs = map S [1..n-1]
    r = [([S i, S i],[]) | i <- [1..n-1]]
    s = [(concat $ replicate 3 [S i, S (i+1)],[]) | i <- [1..n-2]]
    t = [([S i, S j, S i, S j],[]) | i <- [1..n-1], j <- [i+2..n-1]]

And just to check:

> _S 4
([s1,s2,s3],[([s1,s1],[]),([s2,s2],[]),([s3,s3],[]),([s1,s2,s1,s2,s1,s2],[]),([s2,s3,s2,s3,s2,s3],[]),([s1,s3,s1,s3],[])])
> elts $ _S 4
[[],[s1],[s2],[s3],[s1,s2],[s1,s3],[s2,s1],[s2,s3],[s3,s2],[s1,s2,s1],[s1,s2,s3],[s1,s3,s2],[s2,s1,s3],[s2,s3,s2],[s3,s2,s1],[s1,s2,s1,s3],[s1,s2,s3,s2],[s1,s3,s2,s1],[s2,s1,s3,s2],[s2,s3,s2,s1],[s1,s2,s1,s3,s2],[s1,s2,s3,s2,s1],[s2,s1,s3,s2,s1],[s1,s2,s1,s3,s2,s1]]
> length it
24

Anyway, that's it for now. Where I'm heading with this stuff is finite reflection groups and Coxeter groups, but I might take a couple of detours along the way.

String rewriting and Knuth-Bendix completion

2010-05-28T22:08:00.001+01:00

Previously in this blog we have been looking at symmetry groups of combinatorial structures. We have represented these symmetries concretely as permutations - for example, symmetries of graphs as permutations of their vertices. However, mathematicians tend to think about groups more abstractly.

Consider the symmetry group of the square (the cyclic graph C4). It can be generated by two permutations:

> :load Math.Algebra.Group.PermutationGroup
> let a = p [[1,2,3,4]]
> let b = p [[1,2],[3,4]]

We can list various relations that are satisfied by these generators:
a^4 = 1
b^2 = 1
a^3 b = b a

Of course, there are other relations that hold between the generators. However, the relations above are in fact sufficient to uniquely identify the group (up to isomorphism).

Since a and b generate the group, any element in the group can be expressed as a product of a's and b's (and also their inverses, but we'll ignore that for now). However, there are of course an infinite number of such expressions, but only a finite number of group elements, so many of these expressions must represent the same element. For example, since b^2=1, then abba represents the same element as aa.

Given two expressions, it would obviously be helpful to have a method for telling whether they represent the same group element. What we need is a string rewriting system. We can think of expressions in the generators as words in the symbols 'a' and 'b'. And we can reinterpret the relations above as rewrite rules:
"aaaa" -> ""
"bb" -> ""
"aaab" -> "ba"

Each of these rules consists of a left hand side and a right hand side. Given any word in the generator symbols, if we find the left hand side anywhere in the word, we can replace it by the right hand side. For example, in the word "abba", we can apply the rule "bb" -> "", giving "aa".

So, the idea is that given any word in the generator symbols, we repeatedly apply the rewrite rules until we can go no further. The hope is that if two words represent the same group element, then we will end up with the same word after rewriting. We'll see later that there's a bit more to do before that will work, but for the moment, let's at least write some Haskell code to do the string rewriting.

So the first thing we are going to need to do is try to find the left hand side of a rule as a subword within a word. Actually, we want to do a bit more than that - if X is our word, and Y the subword, then we want to find the A and B such that X = AYB.

import qualified Data.List as L

splitSubstring xs ys = splitSubstring' [] xs where
    splitSubstring' ls [] = Nothing
    splitSubstring' ls (r:rs) =
        if ys `L.isPrefixOf` (r:rs)
        then Just (reverse ls, drop (length ys) (r:rs))
        else splitSubstring' (r:ls) rs

Then if our rewrite rule is L -> R, then a single application of the rule consists in replacing L by R within the word:

rewrite1 (l,r) xs =
    case xs `splitSubstring` l of
    Nothing -> Nothing
    Just (a,b) -> Just (a++r++b)

Okay, so suppose we have a rewrite system (that is, a collection of rewrite rules), and a word. Then we want to repeatedly apply the rules until we find that no rule applies:

rewrite rules word = rewrite' rules word where
    rewrite' (r:rs) xs =
        case rewrite1 r xs of
        Nothing -> rewrite' rs xs
        Just ys -> rewrite' rules ys
    rewrite' [] xs = xs

For example:

> :load Math.Algebra.Group.StringRewriting
> rewrite [("aaaa",""),("bb",""),("aaab","ba")] "abba"
"aa"

So far, so good. However, there are some problems with the rewrite system that we constructed above. Suppose that the word we wanted to reduce was "aaabb".
If we apply the rule "aaab" -> "ba", then we have "aaabb" -> "bab".
However, if we apply the rule "bb" -> "", then we have "aaabb" -> "aaa".
Neither "bab" nor "aaa" reduces any further. So we have two problems:
- The same starting word can end up at different end words, depending on the order in which we apply the rules
- We can see from the example that the words "bab" and "aaa" actually represent the same element in our group, but our rewrite system can't rewrite either of them

What can we do about this? Well here's an idea. Let's just add "bab" -> "aaa" as a new rewrite rule to our system. We know that they are equal as elements of the group, so this is a valid thing to do.

That's good, but we still have problems. What about the word "aaaab"?
If we apply the rule "aaaa" -> "", then "aaaab" -> "b"
On the other hand, if we apply the rule "aaab" -> "ba", then "aaaab" -> "aba"

So let's do the same again, and add a new rule "aba" -> "b".

What we're doing here is called the Knuth-Bendix algorithm. Let's take a step back. So in each case, I found a word that could be reduced in two different ways. How did I do that? Well, what I was looking for is two rules with overlapping left hand sides. That is, I was looking for rules L1 -> R1, L2 -> R2, with
L1 = AB
L2 = BC
A pair of rules like this is called a critical pair. If we can find a critical pair, then by looking at the word ABC, we see that
ABC = (AB)C = L1 C -> R1 C
ABC = A(BC) = A L2 -> A R2
So we are justified in adding a new rule R1 C -> A R2

So the Knuth-Bendix algorithm basically says, for each critical pair, introduce a new rule, until there are no more critical pairs. There's a little bit more to it than that:
- We want the rewrite system to reduce the word. That means that we want an ordering on words, and given a pair, we want to make them into a rule that takes the greater to the lesser, rather than vice versa. The most obvious ordering to use is called shortlex: take longer words to shorter words, and if the lengths are equal, use alphabetical ordering.
- Whenever we introduce a new rule, it might be that the left hand side of some existing rule becomes reducible. In that case, the existing rule becomes redundant, since any word that it would reduce can now be reduced by the new rule.

Here's the code:

-- given two strings x,y, find if possible a,b,c with x=ab y=bc
findOverlap xs ys = findOverlap' [] xs ys where
    findOverlap' as [] cs = Nothing
    findOverlap' as (b:bs) cs =
        if (b:bs) `L.isPrefixOf` cs
        then Just (reverse as, b:bs, drop (length (b:bs)) cs)
        else findOverlap' (b:as) bs cs

shortlex x y = compare (length x, x) (length y, y)

ordpair x y =
    case shortlex x y of
    LT -> Just (y,x)
    EQ -> Nothing
    GT -> Just (x,y)

knuthBendix1 rules = knuthBendix' rules pairs where
    pairs = [(lri,lrj) | lri <- rules, lrj <- rules, lri /= lrj]
    knuthBendix' rules [] = rules
    knuthBendix' rules ( ((li,ri),(lj,rj)) : ps) =
        case findOverlap li lj of
        Nothing -> knuthBendix' rules ps
        Just (a,b,c) -> case ordpair (rewrite rules (ri++c)) (rewrite rules (a++rj)) of
                        Nothing -> knuthBendix' rules ps -- they both reduce to the same thing
                        Just rule' -> let rules' = reduce rule' rules
                                          ps' = ps ++
                                                [(rule',rule) | rule <- rules'] ++
                                                [(rule,rule') | rule <- rules']
                                      in knuthBendix' (rule':rules') ps'
   reduce rule@(l,r) rules = filter (\(l',r') -> not (L.isInfixOf l l')) rules

For example:

> knuthBendix1 [("aaaa",""), ("bb",""), ("aaab","ba")]
[("baa","aab"),("bab","aaa"),("aba","b"),("aaaa",""),("bb",""),("aaab","ba")]

A few words about the Knuth-Bendix algorithm
- It is not guaranteed to terminate. Every time we introduce a new rule, we have the potential to create new critical pairs, and there are pathological examples where this goes on forever
- The algorithm can be made slightly more efficient, by doing things like choosing to process shorter critical pairs first. In the HaskellForMaths library, a more efficient version is given, called simply "knuthBendix"

Back to the example. So Knuth-Bendix has found three new rules. The full system, with these new rules added, has no more critical pairs. As a consequence, it is a confluent rewrite system - meaning that if you start at some given word, and reduce it using the rules, then it doesn't matter in what order you apply the rules, you will always end up at the same word. This word that you end up with can therefore be used as a normal form.

This allows us to "solve the word problem" for this group. That is, given any two words in the generator symbols, we can find out whether they represent the same group element by rewriting them both, and seeing if they end up at the same normal form. For example:

> let rules = knuthBendix [("aaaa",""), ("bb",""), ("aaab","ba")]
> rewrite rules "aaaba"
"aab"
> rewrite rules "baabb"
"aab"
> rewrite rules "babab"
"b"

So we see that "aaaba" and "baabb" represent the same group element, whereas "babab" represents a different one. (If you want, you could go back and check this using the original permutations.)

We can even list (the normal forms of) all elements of the group. What we do is start with the empty word (which represents the identity element of the group), and then incrementally build longer and longer words. At each stage, we look at all combinations that can be formed by pre-pending a generator symbol to a word from the preceding stage. However, if we ever come across a word which can be reduced, then we know that it - and any word that could be formed from it at a later stage - is not a normal form, and so can be discarded. Here's the code:

nfs (gs,rs) = nfs' [[]] where
    nfs' [] = [] -- we have run out of words
    nfs' ws = let ws' = [g:w | g <- gs, w <- ws, (not . isNewlyReducible) (g:w)]
              in ws ++ nfs' ws'
    isNewlyReducible w = any (`L.isPrefixOf` w) (map fst rs)

elts (gs,rs) = nfs (gs, knuthBendix rs)

For example:

> elts (['a','b'], [("aaaa",""), ("bb",""), ("aaab","ba")])
["","a","b","aa","ab","ba","aaa","aab"]

As expected, we have eight elements.

That's enough for now. Next time (hopefully) I'll look at some more examples.

Block systems and block homomorphism

2010-04-25T20:55:00.000+01:00

Recently in this blog, we looked at the strong generating set (SGS) algorithm for permutation groups, and how we can use it to investigate the structure of groups. Last time, we saw how to partially "factor" intransitive groups, using the transitive constituent homomorphism. (Recall that by "factoring" a group G, we mean finding a proper normal subgroup K, and consequently also a quotient group G/K - which is equivalent to finding a proper homomorphism from G.) This time, I want to do the same for imprimitive groups. So, what is an imprimitive group?

Well, given a permutation group acting on a set X, it can happen that X consists of "blocks" Y1, Y2, ... of points which always "move together". That is, a subset Y of X is a block if for all g in G, Y^g (the image of Y under the action of g) is either equal to Y or disjoint from it. A full set of blocks (that is, blocks Y1, Y2, ... which are disjoint, and whose union is the whole of X) is called a block system.

For example, suppose that X is the vertices of the hexagon. The symmetry group of the hexagon is the dihedral group D12, generated by a rotation and a reflection:

> :load Math.Algebra.Group.Subquotients
> mapM_ print $ _D 12
[[1,2,3,4,5,6]]
[[1,6],[2,5],[3,4]]

A block system for the hexagon is shown below. The blocks are the pairs of opposite vertices. You can verify that they satisfy the definition of blocks: any symmetry must take a pair of opposite points either to itself, or to another pair disjoint from it.

A given group can have more than one block system. Here is another block system for the hexagon. The blocks are the two equilateral triangles formed by the vertices.

There are also the trivial block systems, consisting of either just one block containing all the points, or a block for each point. From now on, we will silently exclude these.

So, I was meant to be telling you what an imprimitive group is. Well, it's just a group which has a non-trivial block system. Conversely, a primitive group is one which has no non-trivial block system.

When we have an imprimitive group, we will be able to form a homomorphism - and hence factor the group - by considering the induced action of the group on the blocks. But I'm jumping ahead. First we need to write some Haskell code - to find block systems.

The idea is to write a function that, given a pair of points Y = {y1,y2} in X (or indeed any subset Y of X), can find the smallest block containing Y. The way it works is as follows. We start by supposing that each point is in a block of its own, except for the points in Y. We initialise a map, with the points in X as keys, and the blocks as values, where we represent a block by its least element.

Now, suppose that we currently think that the minimal block is Y = {y1,y2,...}. What we're going to do is work through the elements of Y, and work through the generators of G, trying to find a problem. So suppose that we have got as far as some element y of Y, and some generator g of G. We know that y is in the same block as y1, and what we have to check is that y^g is in the same block as y1^g. So we look up their representatives in the map, and check that they're the same. If they're not, then we need to merge the two classes. Here's the code (it's a little opaque, but it's basically doing what I just described).

minimalBlock gs ys@(y1:yt) = minimalBlock' p yt gs where
    xs = foldl union [] $ map supp gs
    p = M.fromList $ [(yi,y1) | yi <- ys] ++ [(x,x) | x <- xs \\ ys]
    minimalBlock' p (q:qs) (h:hs) =
        let r = p M.! q         -- representative of class containing q
            k = p M.! (q .^ h)  -- rep of class (q^h)
            l = p M.! (r .^ h)  -- rep of class (r^h)
        in if k /= l -- then we need to merge the classes
           then let p' = M.map (\x -> if x == l then k else x) p
                    qs' = qs ++ [l]
                in minimalBlock' p' (q:qs') hs
           else minimalBlock' p (q:qs) hs
    minimalBlock' p (q:qs) [] = minimalBlock' p qs gs
    minimalBlock' p [] _ =
        let reps = toListSet $ M.elems p
        in L.sort [ filter (\x -> p M.! x == r) xs | r <- reps ]

Once we have this function, then finding the block systems is simple - just take each pair {x1,xi} from X, and find the minimal block containing it.

blockSystems gs
    | isTransitive gs = toListSet $ filter (/= [x:xs]) $ map (minimalBlock gs) [ [x,x'] | x' <- xs ]
    | otherwise = error "blockSystems: not transitive"
    where x:xs = foldl union [] $ map supp gs

If we have an SGS for G, then we can do slightly better. For suppose that within the stabiliser G_x1, there is an element taking xi to xj. Then clearly xi and xj must be in the same minimal block. So in fact, we need only consider pairs {x1,xi}, with xi the minimal element of each orbit in G_x1. (Of course, the point is that if we have an SGS for G, then we can trivially list a set of generators for G_x1.)

blockSystemsSGS gs = toListSet $ filter (/= [x:xs]) $ map (minimalBlock gs) [ [x,x'] | x' <- rs ]
    where x:xs = foldl union [] $ map supp gs
          hs = filter (\g -> x < minsupp g) gs -- sgs for stabiliser Gx
          os = orbits hs
          rs = map head os ++ (xs \\ L.sort (concat os)) -- orbit representatives, including singleton cycles

Let's test it:

> mapM_ print $ blockSystems $ _D 12
[[1,3,5],[2,4,6]]
[[1,4],[2,5],[3,6]]

Okay, so given a group, we can find its non-trivial block systems, if any. What next? Well, as I hinted earlier, this enables us to factor the group. For if there is a non-trivial block system, then the action of the group on the points induces a well-defined action on the blocks. This induced action gives us a homomorphism from our original group G, a subgroup of Sym(X), to another group H, a subgroup of Sym(B), where B is the set of blocks.

So as we did last time, we can find the kernel and image of the homomorphism, and thus factor the group. How do we do that?

Well, it's simple. In the following code, the function lr takes a group element acting on the points, and returns a group element acting on the blocks (in the Left side) and the points (in the Right side) in an Either union. If we do this to all the group generators, and then find an SGS, then as the Left blocks sort before the Right points, then the SGS will split neatly into two parts:
- The initial segment of the SGS will consist of elements which move the Left blocks. If we restrict their action to just the blocks, we will have an SGS for the image of the homomorphism, acting on the blocks.
- The final segment of the SGS will consist of elements which fix all the Left blocks. These elements move points but not blocks, so they form an SGS for the kernel of the homomorphism.

blockHomomorphism' gs bs = (ker,im) where
    gs' = sgs $ map lr gs
    lr g = fromPairs $ [(Left b, Left $ b -^ g) | b <- bs] ++ [(Right x, Right y) | (x,y) <- toPairs g]
    ker = map unRight $ dropWhile (isLeft . minsupp) gs' -- stabiliser of the blocks
    im = map restrictLeft $ takeWhile (isLeft . minsupp) gs' -- restriction to the action on blocks

blockHomomorphism gs bs
    | bs == closure bs [(-^ g) | g <- gs] -- validity check on bs
        = blockHomomorphism' gs bs

Let's try it out on our two block systems for the hexagon:

> blockHomomorphism (_D 12) [[1,4],[2,5],[3,6]]
([[[1,4],[2,5],[3,6]]],
 [[[[1,4],[2,5],[3,6]]],[[[2,5],[3,6]]]])

I've formatted the output for clarity. The first line is (an SGS for) the kernel, consisting of elements of D12 which permute points within the blocks, without permuting the blocks. In this case, the kernel is generated by the 180 degree rotation, which swaps the points within each pair. The second line is (an SGS for) the image, consisting of the induced action of D12 on the blocks. In this case, we have the full permutation group S3 acting on the three pairs of points.

> blockHomomorphism (_D 12) [[1,3,5],[2,4,6]]
([[[1,5,3],[2,6,4]],[[2,6],[3,5]]],
 [[[[1,3,5],[2,4,6]]]])

In this case the kernel is generated by a 120 degree rotation and a reflection, and consists of all group elements which send odd points to odd and even points to even, thus preserving the blocks. The image has only one non-trivial element, which just swaps the two blocks.

Armed with this new tool, let's have another look at Rubik's cube. Recall that we labelled the faces of the cube as follows:

Last time, we split the Rubik cube group into two homomorphic images - a group acting on just the corner faces, and a group acting on just the edge faces. Let's look for block systems in these groups:

> :load Math.Projects.Rubik
> let [cornerBlocks] = blockSystems imCornerFaces
> let [edgeBlocks] = blockSystems imEdgeFaces
> cornerBlocks
[[1,17,23],[3,19,41],[7,29,31],[9,33,47],[11,21,53],[13,43,51],[27,37,59],[39,49,57]]
> edgeBlocks
[[2,18],[4,26],[6,44],[8,32],[12,52],[14,22],[16,42],[24,56],[28,34],[36,48],[38,58],[46,54]]

It's obvious really - in the corner group, we have a block system with blocks consisting of the three corner faces that belong to the same corner piece, and in the edge group, we have a block system with blocks consisting of the two edge faces that belong to the same edge piece. Furthermore, these are the only block systems.

So we can form the kernel and image under the block homomorphism:

> let (kerCornerBlocks,imCornerBlocks) = blockHomomorphism imCornerFaces cornerBlocks
> let (kerEdgeBlocks,imEdgeBlocks) = blockHomomorphism imEdgeFaces edgeBlocks

If we look at the sizes of these groups, the structure will be obvious:

> orderSGS kerCornerBlocks
2187
> orderSGS imCornerBlocks
40320

These are 3^7, and 8! respectively. The kernel is the permutations of the corner faces which leave the corner blocks where they are. It turns out that whenever you twist one corner block, you must untwist another. So when you have decided what to do with seven corners, the eighth is determined - hence 3^7. For the image, we have eight blocks, and 8! permutations of them, so this must be the full symmetry group S8 - meaning that we can perform any rearrangement of the corner blocks that is desired.

> orderSGS kerEdgeBlocks
2048
> orderSGS imEdgeBlocks
479001600

These are 2^11 and 12! respectively. For the kernel, whenever we flip one edge piece we must also flip another. So when we have decided what to do with eleven edges, the twelfth is determined - hence 2^11. For the image, we have twelve pieces, and 12! permutations of them, so we have the full symmetry group S12 on edge blocks.

That's it.

Incidentally, my references for this material are:
- Holt, Handbook of Computational Group Theory
- Seress, Permutation Group Algorithms
both of which are very good - but expensive.

These books, particularly the latter, go on to describe further algorithms that can be used to factor even transitive primitive groups, enabling us to arrive at a full decomposition of a group into simple groups. Unfortunately, the algorithms get a bit more complicated after this, and I haven't yet implemented the rest in HaskellForMaths.

Transitive constituent homomorphism

2010-03-23T20:44:00.000Z

[New version HaskellForMaths 0.2.1 available here]

Last time, we looked at what it means to be a strong generating set (SGS) for a permutation group, and how to find one. We saw that with an SGS we can easily calculate the number of elements in the group, and test whether an arbitrary permutation is a member of the group.

This time, I want to start showing how strong generating sets can be used as the basis for further algorithms to investigate the structure of groups. In a previous post, I explained about normal subgroups and quotient groups. Over the next couple of posts, I want to show how, for groups having the appropriate form, we can find some of the easier normal subgroups and quotient groups. Specifically, in this post I want to look at intransitive groups, and in the next at imprimitive groups.

Many of the groups we have looked at so far - such as the symmetries of the pentagon, cube, and so on - have been "transitive" - meaning that given any two points, there is a group element which takes one to the other. (Or to put it another way, all the points "look the same".)

However, it is easy to find groups that are not transitive on their points. For example, if we take any graph which is not regular (that is, it has vertices of different valencies), then the symmetry group is not transitive (because no symmetry can take a vertex of one valency to a vertex of another valency).

In order to demonstrate the code that I want to look at this week, I'm going to need an example - so here are the generators of a group G:

> :load Math.Algebra.Group.Subquotients
> let gs = map p [ [[1,2,3]], [[4,5,6]], [[1,2],[4,5]] ]

We can think of this group as the symmetry group of some kind of plastic puzzle, as shown below. The puzzle consists of a blue and a red triangle. Valid moves consist of rotating the blue triangle, rotating the red triangle, or performing a double flip, which swaps a pair of blue points and at the same time a pair of red points.

It's a small group, so let's have a look at the elements:

> mapM_ print $ elts gs
[]
[[1,2],[4,5]]
[[1,2],[4,6]]
[[1,2],[5,6]]
[[1,2,3]]
[[1,2,3],[4,5,6]]
[[1,2,3],[4,6,5]]
[[1,3,2]]
[[1,3,2],[4,5,6]]
[[1,3,2],[4,6,5]]
[[1,3],[4,5]]
[[1,3],[4,6]]
[[1,3],[5,6]]
[[2,3],[4,5]]
[[2,3],[4,6]]
[[2,3],[5,6]]
[[4,5,6]]
[[4,6,5]]

Notice that the group is not transitive: there is no move that takes a blue point to a red point, or vice versa.

So we have a group G, acting on a set X. If G is not transitive on X, then we can write X as a disjoint union of subsets A, B, C, ..., such that G is transitive on each of the subsets. In our example, we can take A = [1,2,3], B = [4,5,6], and we see that G is transitive on A, and on B. A and B are called the transitive constituents (or the orbits) of G.

It's fairly straightforward to write Haskell code to find the orbits of a permutation group. This allows us to find:

> orbits gs
[[1,2,3],[4,5,6]]

[The implementation is left as an exercise.]

Okay, so what's so interesting about intransitive groups? Well, an intransitive group always has a non-trivial normal subgroup, so we can always "factor" it into smaller groups. How?

Well, take one of the transitive constituents. In our example, let's take [1,2,3]. Then we can consider the restriction of the action of the group to just this constituent. For our generators, the restriction gives us:

[[1,2,3]] -> [[1,2,3]]
[[4,5,6]] -> []
[[1,2],[4,5]] -> [[1,2]]

It might look as if the first generator got taken to itself - but it is important to realise that these are elements of different groups. This would have been clearer if we had written out the permutations in row notation:

1 2 3 4 5 6 -> 1 2 3
2 3 1 4 5 6    2 3 1

Or if we had included singleton cycles in the cycle notation:

[[1,2,3],[4],[5],[6]] -> [[1,2,3]]

So this restriction of the action is a map from our group - a subgroup of Sym([1..6]) - to another group - a subgroup of Sym([1..3]).

If we call this map f, then it should be obvious that f has the following properties:
- f(1) = 1 (where 1 here means the identity element in the respective groups, not the point 1)
- f(g*h) = f(g)*f(h) (consider their action on a point x)
- f(g^-1) = f(g)^-1

Hence f is a homomorphism - a function between groups that preserves the group structure.

The kernel of a homomorphism is the set {g <- G | f(g) = 1}. The kernel is always a normal subgroup:
It's obviously a subgroup. (Why?)
And f(h^-1 g h) = f(h)^-1 f(g) f(h) (by the homomorphism properties).
So if f(g) = 1, then f(h^-1 g h) = 1.
In other words, if g is in the kernel, then so is h^-1 g h - which is just the condition for the kernel to be a normal subgroup.

The image of a homomorphism is the set {f(g) | g <- G}. This is naturally isomorphic to the quotient group G/K, where K is the kernel: since for k <- K, f(gk) = f(g) f(k) = f(g) (since f(k) = 1) - so elements of the image are in one-to-one correspondence with cosets of the kernel, and it is easy to check that the group operations correspond too.

Hence any (non-trivial) homomorphism induces a "factorisation" of the group into smaller groups K and G/K - the kernel and the image.

In our example, the kernel is < [[4,5,6]] > <= Sym([1..6]), and the image is < [[1,2,3]], [[1,2]] > = Sym([1..3]).

We would like a Haskell function that can find the kernel and image of a transitive constituent homomorphism for us, like so:

> transitiveConstituentHomomorphism gs [1,2,3]
([[[4,5,6]]],[[[1,2,3]],[[1,2]]])

So how do we do that? Well, the basic idea is that we will convert all our permutations on X into permutations on Either A B, where A is the transitive constituent we are restricting to, and B is the rest. In our group, we have X = [1..6], A = [1,2,3], B = [4,5,6], so for example we will convert [[1,2],[4,5]] into [[Left 1, Left 2], [Right 4, Right 5]].

Next, we construct a strong generating set for the new group. Since the SGS algorithm looks for point stabilisers /in order/ - it will first find point stabilisers for all the Left points, and then for all the Right points. So our SGS will split into two parts:
- A first part consisting of elements that move some of the Left points. If we restrict this part to its action on the Left points, we will have an SGS for the image.
- A second part consisting of elements that move only the Right points. This part is an SGS for the kernel.

For example:

> mapM_ print $ sgs $ map p [ [[Left 1, Left 2, Left 3]], [[Right 4, Right 5, Right 6]], [[Left 1, Left 2],[Right 4, Right 5]] ]
[[Left 1,Left 2,Left 3],[Right 4,Right 6,Right 5]]
[[Left 2,Left 3],[Right 4,Right 5]]
[[Right 4,Right 6,Right 5]]

[Note that as our SGS algorithm uses randomness, we might get a different set of generators each time we call it.]

Okay, so we know what to do. Time for some code. We need first of all a few utility functions for packing and unpacking our group elements into the Left / Right parts of the Either type.

isLeft (Left _) = True
isLeft (Right _) = False

isRight (Right _) = True
isRight (Left _) = False


unRight = fromPairs . map (\(Right a, Right b) -> (a,b)) . toPairs

restrictLeft g = fromPairs [(a,b) | (Left a, Left b) <- toPairs g]
-- note that this is doing a filter - taking only the left part of the action - and a map, unLefting

Then the code is simple:


transitiveConstituentHomomorphism gs delta
    | delta == closure delta [(.^ g) | g <- gs] -- delta is a transitive constituent
       = transitiveConstituentHomomorphism' gs delta

transitiveConstituentHomomorphism' gs delta = (ker, im) where
    gs' = sgs $ map (fromPairs . map (\(a,b) -> (lr a, lr b)) . toPairs) gs
    -- as delta is a transitive constituent, we will always have a and b either both Left or both Right
    lr x = if x `elem` delta then Left x else Right x
    ker = map unRight $ dropWhile (isLeft . minsupp) gs' -- pointwise stabiliser of delta
    im = map restrictLeft $ takeWhile (isLeft . minsupp) gs' -- restriction of the action to delta

That's it.

Okay, so the transitive constituent homomorphism gives us a way to "take apart" intransitive groups. Let's use it to take another look at the Rubik cube group.

Remember how we labelled the faces of the Rubik cube:

Then the generators of the Rubik cube group are:

f = p [[ 1, 3, 9, 7],[ 2, 6, 8, 4],[17,41,33,29],[18,44,32,26],[19,47,31,23]]
b = p [[51,53,59,57],[52,56,58,54],[11,27,39,43],[12,24,38,46],[13,21,37,49]]
l = p [[21,23,29,27],[22,26,28,24],[ 1,31,59,11],[ 4,34,56,14],[ 7,37,53,17]]
r = p [[41,43,49,47],[42,46,48,44],[ 3,13,57,33],[ 6,16,54,36],[ 9,19,51,39]]
u = p [[11,13,19,17],[12,16,18,14],[ 1,21,51,41],[ 2,22,52,42],[ 3,23,53,43]]
d = p [[31,33,39,37],[32,36,38,34],[ 7,47,57,27],[ 8,48,58,28],[ 9,49,59,29]]

rubikCube = [f,b,l,r,u,d]

What orbits does it have?

> :load Math.Projects.Rubik
> orbits rubikCube
[[1,3,7,9,11,13,17,19,21,23,27,29,31,33,37,39,41,43,47,49,51,53,57,59],[2,4,6,8,12,14,16,18,22,24,26,28,32,34,36,38,42,44,46,48,52,54,56,58]]

It's obvious really: the odd numbered points are corner faces, and the even numbered points are edge faces - and there is no valid move that can take a corner face to an edge face, or vice versa. So the group is intransitive. So we can define:

[cornerFaces,edgeFaces] = orbits rubikCube

We then have two different ways to split the group, depending which transitive constituent we take:

(kerCornerFaces,imCornerFaces) = transitiveConstituentHomomorphism rubikCube cornerFaces

(kerEdgeFaces,imEdgeFaces) = transitiveConstituentHomomorphism rubikCube edgeFaces

Note that:

> orderSGS imCornerFaces * orderSGS imEdgeFaces
86504006548979712000
> let rubikSGS = sgs rubikCube
> orderSGS rubikSGS
43252003274489856000

Why are they not equal? It is because in fact we can't operate on corners and edges totally independently. The following examples show that we can't swap a pair of corners without also swapping a pair of edges, and vice versa:

> isMemberSGS rubikSGS (p [[1,3],[17,41],[19,23],[2,6],[18,44]])
True
> isMemberSGS rubikSGS (p [[1,3],[17,41],[19,23]])
False
> isMemberSGS rubikSGS (p [[2,6],[18,44]])
False

(This is similar to the way that in our original example, we couldn't swap a pair of blue points without also swapping a pair of red points.)

That's it for now. Next time, I'll show how we can factor the corner and edge groups still further, to arrive at a fairly complete understanding of the structure of the Rubik cube group.

(By the way, I would like to acknowledge a debt to the following example from the GAP system: http://www.gap-system.org/Doc/Examples/rubik.html )

How to find a strong generating set

2010-02-14T20:10:00.001Z

Nearly three months since my last post - sorry! I have an excuse - we've been moving house - but the truth is I've had a bit of writer's block over this post. Anyway...

Previously in this blog we have been looking at permutation groups, and especially those which arise as the symmetry groups of graphs or other combinatorial objects. We developed a "graphAuts" function, which finds generators for the group of symmetries of a graph. In fact, "graphAuts" finds a strong generating set or SGS, which is a generating set of a special form, making various calculations in permutation groups relatively easy. What I want to look at this week, is how we can find a strong generating set for an arbitrary permutation group, given generators for the group.

Let's remind ourselves what a strong generating set looks like. Here's q 3, the graph of the cube.

The graphAuts function returns us a strong generating set for the symmetry group of the cube:

> :load Math.Combinatorics.GraphAuts
> mapM_ print $ graphAuts $ q 3
[[0,1],[2,3],[4,5],[6,7]]
[[0,2,3,1],[4,6,7,5]]
[[0,4,6,7,3,1],[2,5]]
[[1,2],[5,6]]
[[1,4,2],[3,5,6]]
[[2,4],[3,5]]

So what is special about this set of symmetries? Well, first of all, the strong generating set is a generating set for the group - all symmetries of the cube can be expressed as products of these generators (and their inverses). But notice how the SGS is composed of a number of "levels". The first level consists of some elements that move 0. The second level consists of some elements that fix 0 but move 1. The third level consists of an element that fixes 0 and 1 but moves 2. So we can think of there being a 0-level, a 1-level, and a 2-level. The sequence [0,1,2] is called the base for the SGS.

Let's call our group G. Now, if we discard the first level of the SGS (the 0-level), the elements that remain form an SGS for the subgroup of elements that fix 0, which is called the point stabiliser of 0, denoted G₀. If we also discard the second level (the 1-level), the element that remains is an SGS for the subgroup of elements that fix both 0 and 1, which is called the pointwise stabiliser of {0,1}, denoted G_{0,1}. So an SGS defines a sequence of subgroups called a point stabiliser chain for the group:

G_0 = G                    >  G_1 = G_{0}        >  G_2 = G_{0,1}  >  G_{0,1,2} = 1
[[0,1],[2,3],[4,5],[6,7]]     [[1,2],[5,6]]         [[2,4],[3,5]]
[[0,2,3,1],[4,6,7,5]]         [[1,4,2],[3,5,6]]
[[0,4,6,7,3,1],[2,5]]         [[2,4],[3,5]]
[[1,2],[5,6]]
[[1,4,2],[3,5,6]]
[[2,4],[3,5]]

Now, consider a pair of successive subgroups within this chain, for example, G_{0} and G_{0,1}. G_{0} consists of some elements that fix 1, some elements that send 1 to 2, and some elements that send 1 to 4. The elements that fix 1 are just the subgroup G_{0,1}, the point stabiliser of 1 within G_{0}. The elements that send 1 to 2, and those that send 1 to 4, are cosets of G_{0,1} within G_{0}. We look for a representative of each coset - that is, an element that sends 1 to 2 - [[1,2],[5,6]], an element that sends 1 to 4 - [[1,4,2],[3,5,6]], and for an element that sends 1 to 1 (fixes 1), we can take the identity. Once we have done that, every element of G_{0} can be expressed the product of an element of G_{0,1} with one of these representatives.

Just to show you what I mean, here are some representatives for each link in the chain for the cube group:

U_0                              U_1                      U_2
(0,[])                           (1,[])                   (2,[])
(1,[[0,1],[2,3],[4,5],[6,7]])    (2,[[1,2],[5,6]])        (4,[[2,4],[3,5]])
(2,[[0,2,3,1],[4,6,7,5]])        (4,[[1,4,2],[3,5,6]])
(3,[[0,3],[1,2],[4,7],[5,6]])
(4,[[0,4,6,7,3,1],[2,5]])
(5,[[0,5,6,3],[1,4,7,2]])
(6,[[0,6,3],[1,4,7]])
(7,[[0,7],[1,6],[2,5],[3,4]])

At each link in the point stabiliser chain, we have a base element b that we are going to stabilise. For each point within the orbit of b, we find a group element that takes b to that point. This set of group elements, U_b, is a set of coset representatives - a transversal for the cosets.

The Haskell code to find these transversals, given an SGS, is fairly straightforward:

baseTransversalsSGS gs = [let hs = filter ( (b <=) . minsupp ) gs in (b, cosetRepsGx hs b) | b <- bs]
    where bs = toListSet $ map minsupp gs

cosetRepsGx gs x = cosetRepsGx' gs M.empty (M.singleton x 1) where
    cosetRepsGx' gs interior boundary
        | M.null boundary = interior
        | otherwise =
            let interior' = M.union interior boundary
                boundary' = M.fromList [(p .^ g, h*g) | g <- gs, (p,h) <- M.toList boundary] M.\\ interior'
            in cosetRepsGx' gs interior' boundary'

Here it is in action:

> :load Math.Algebra.Group.RandomSchreierSims
> let gs = map p [ [[0,1],[2,3],[4,5],[6,7]], [[0,2,3,1],[4,6,7,5]], [[0,4,6,7,3,1],[2,5]], [[1,2],[5,6]], [[1,4,2],[3,5,6]], [[2,4],[3,5]] ]
> mapM_ print $ baseTransversalsSGS gs
(0, fromList [(0,[]), (1,[[0,1],[2,3],[4,5],[6,7]]), (2,[[0,2,3,1],[4,6,7,5]]), (3,[[0,3],[1,2],[4,7],[5,6]]), (4,[[0,4,6,7,3,1],[2,5]]), (5,[[0,5,6,3],[1,4,7,2]]), (6,[[0,6,3],[1,4,7]]), (7,[[0,7],[1,6],[2,5],[3,4]])])
(1, fromList [(1,[]), (2,[[1,2],[5,6]]), (4,[[1,4,2],[3,5,6]])])
(2, fromList [(2,[]), (4,[[2,4],[3,5]])])

Okay, so what do these transversals give us? Well, hopefully it is obvious (by induction) that every element in our original group G can be expressed as a product u2*u1*u0, taking an element ui from each transversal. Moreover, this expression is unique.

This is the key to the usefulness of strong generating sets. For example, it means we can easily calculate the order of the group (the number of elements). Just multiply the sizes of the transversals. For the cube group, this tells us that we have 8 * 3 * 2 = 48 elements. We can also use the transversals to create a search tree for the group (more on this some other time). Of particular importance, we can test an arbitrary permutation for membership in the group, by "sifting" it through the transversals, as follows.

Given a permutation g, we look at its action on the first base, b1, and see if we can find an element in the b1-transversal with the same action. If we can, say u1, then we replace g by g*u1^-1, and proceed to the second base. We now look at the action of g*u1^-1 on b2, and see if there is an element in the b2-transversal with the same action. If g is in our group G, then g = un*...*u1, so at the end of this process, we will get 1. Conversely, if g is not in our group, then at some stage we will fail to find a transversal element matching the action of our g (including possibly the case where we go through all the transversals, but are still left at the end with a g /= 1).

Here's the code:

isMemberSGS gs h = sift bts h where
    bts = baseTransversalsSGS gs
    sift _ 1 = True
    sift ((b,t):bts) g = case M.lookup (b .^ g) t of
                         Nothing -> False
                         Just h -> sift bts (g * inverse h)
    sift [] _ = False

Okay, so we've seen what a strong generating set looks like, and why they're useful. Now, suppose someone gives us a set of generators for a group. How do we find a strong generating set for the group?

Well, it's actually fairly simple. We start with a set of empty transversals, corresponding to an empty SGS. What we're going to try to do is add elements to the SGS (and hence to the transversals), until we end up with a full SGS for the group. First of all, we need a slightly modified version of the sift function from above:

sift _ 1 = Nothing
sift ((b,t):bts) g = case M.lookup (b .^ g) t of
                     Nothing -> Just g
                     Just h -> sift bts (g * inverse h)
sift [] g = Just g -- g == 1 case already caught above

In this version, instead of returning either True or False, we return either Nothing, indicating success, or Just g, indicating failure, where g is the part-sifted element that we had at the point of failure.

Okay, so to find an SGS, what we do is feed random elements of the group through the transversals, using the sifting procedure. (How do we generate random elements of the group? - well roughly, we just form random products of the generators - but see below.) If our random element sifts through, then we move on to the next random element. If our random element doesn't sift through, then at some transversal, we have a group element whose action on a base is not matched by any element in the transversal. If this happens, what we have to do is add this group element as a new element to the SGS, and recalculate the transversals. In that way, we will gradually grow the SGS. What we're hoping is that if we look at enough random elements, we'll end up with an SGS for the whole group.

This algorithm is called random Schreier-Sims. It's a Monte-Carlo algorithm - meaning it can return the wrong answer: We might stop too soon, and return an SGS which only generates a subgroup, not the whole group. However, we can make this unlikely by choosing enough random elements.

Okay, so first of all, given some generators for a group, we need to be able to generate random elements of the group. One way to do this is called the "product replacement algorithm". What we will do is create an array containing 10 random elements of the group. Then every time we're asked for another random element, we will randomly replace one element of the array, by multiplying it either on left or right, by either another element or its inverse. The replacement element will thus be a new random element of the group. Ok, but how do we get 10 random elements in the array in the first place? Well, what we do is start off by putting the generators into the array, repeated as many times as necessary - eg a,b,c,a,b,c,a,b,c,a. Then we just run the product replacement a few times to "mix up" the array.

Here's the code. (Note, in this version, we maintain an eleventh element (the zeroth) in the array as an accumulator for results we are going to return - this has been found empirically to give better results.)

initProdRepl :: (Ord a, Show a) => [Permutation a] -> IO (Int, IOArray Int (Permutation a))
initProdRepl gs =
    let n = length gs
        r = max 10 n -- if we have more than 10 generators, we have to use a larger array
        xs = (1:) $ take r $ concat $ repeat gs
    in do xs' <- newListArray (0,r) xs
          replicateM_ 60 $ nextProdRepl (r,xs') -- perform initial mixing
          return (r,xs')

nextProdRepl :: (Ord a, Show a) => (Int, IOArray Int (Permutation a)) -> IO (Maybe (Permutation a))
nextProdRepl (r,xs) = -- r will usually be 10
    do s <- randomRIO (1,r)
       t <- randomRIO (1,r)
       u <- randomRIO (0,3 :: Int)
       out <- updateArray xs s t u
       return out

updateArray xs s t u =
    let (swap,invert) = quotRem u 2 in
    if s == t
    then return Nothing
    else do
        x_0 <- readArray xs 0
        x_s <- readArray xs s
        x_t <- readArray xs t
        let x_s' = mult (swap,invert) x_s x_t
            x_0' = mult (swap,0) x_0 x_s'
        writeArray xs 0 x_0'
        writeArray xs s x_s'
        return (Just x_0')
    where mult (swap,invert) a b = case (swap,invert) of
                                   (0,0) -> a * b
                                   (0,1) -> a * b^-1
                                   (1,0) -> b * a
                                   (1,1) -> b^-1 * a

One thing to point out: We have to select a random element to replace - x_s - and another random element to multiply by - x_t. We're not allowed to have s=t. The way I get round this is to wrap Maybe around the return type, and return Nothing if I happen to pick s=t. With a little bit of work I could probably do better.

Okay, so we can generate random elements from the generators. Now, to find a strong generating set, we start out with an empty set of transversals, and keep sifting random elements through the transversals, updating the transversals when we find elements that don't sift:

sgs :: (Ord a, Show a) => [Permutation a] -> [Permutation a]
sgs gs = toListSet $ concatMap snd $ rss gs

rss gs = unsafePerformIO $
    do (r,xs) <- initProdRepl gs
       rss' (r,xs) (initLevels gs) 0

rss' (r,xs) levels i
    | i == 25 = return levels -- stop if we've had 25 successful sifts in a row
    | otherwise = do g <- nextProdRepl (r,xs)
                     let (changed,levels') = updateLevels levels g
                     rss' (r,xs) levels' (if changed then 0 else i+1)
-- if we currently have an sgs for a subgroup of the group, then it must have index >= 2
-- so the chance of a random elt sifting to identity is <= 1/2

initLevels gs = [((b,M.singleton b 1),[]) | b <- bs]
    where bs = toListSet $ concatMap supp gs

updateLevels levels Nothing = (False,levels)
updateLevels levels (Just g) =
    case sift (map fst levels) g of
    Nothing -> (False, levels)
    Just g' -> (True, updateLevels' [] levels g' (minsupp g'))

updateLevels' ls (r@((b,t),s):rs) h b' =
    if b == b'
    then reverse ls ++ ((b, cosetRepsGx (h:s) b), h:s) : rs
    else updateLevels' (r:ls) rs h b'

That's it.

This algorithm makes it possible to work with much larger permutation groups. We already saw it in use in an earlier post, where I found an SGS for the Rubik's cube group, in order to calculate the number of possible positions (approximately 4*10^19). Over the next few weeks I want to look at how we can use strong generating sets to investigate the structure of groups.

Three new modules in HaskellForMaths

2009-11-18T21:40:00.000Z

It's been a little while since I posted here - and the reason is that I've been too busy writing code. The fruits of my labour are three new modules for HaskellForMaths, which I've uploaded in version 0.2.0. In due course, each of these deserves a blog post in its own right, but for the moment, I thought I'd give a quick overview.

Math.Algebra.Group.RandomSchreierSims

In our investigation of graph automorphisms, we've come across the concept of a strong generating set (SGS) for a permutation group. This is a generating set of a special form, which makes certain calculations particularly easy - for example, calculating the order of the group (the number of elements it contains). Our graphAuts algorithm naturally produced a strong generating set for the group. But what if we are just given any old set of generators for a permutation group.

The random Schreier-Sims algorithm is a fast Monte-Carlo algorithm for finding a strong generating set, given a set of generators. Monte-Carlo means that there is a small probability that it will return the wrong answer - in this case, it might return an SGS for a subgroup of the group we are interested in, instead of the whole group. In practice, this isn't an issue. First, we can make the probability of error as small as we like. Second, we usually know the order of the group, so can check whether the SGS is right.

The main function provided in Math.Algebra.Group.RandomSchreierSims is sgs, which takes a set of generators and returns a strong generating set.

(We already have an sgs function in Math.Algebra.Group.SchreierSims. This is guaranteed to give the right answer - but can be significantly slower. So, you have a choice which one to use.)

Math.Projects.MiniquaternionGeometry

A projective plane is an incidence structure of points and lines satisfying the following rules:
- Given any two distinct points, there is exactly one line which passes through both
- Given any two distinct lines, there is exactly one point which lies on both
- There exists a quadrangle: a set of four points, no three of which are collinear

(There is another definition of projective plane, in terms of designs, but we haven't covered those yet.)

For any prime power q = p^n, the projective geometry PG(2,Fq) over the finite field Fq is a projective plane. However, it turns out that there are projective planes which are not of this form, associated with algebraic systems called near-fields, which are nearly fields but not quite.

In Math.Projects.MiniquaternionGeometry I construct PG(2,F9), together with three "non-Desarguesian" planes of order 9, based on the near-field J9 of order 9. The name Miniquaternion Geometry comes from the book of that name by Room and Kirkpatrick, because the multiplicative structure of J9 is the same as that of the unit quaternions.

Math.Combinatorics.LatinSquares

A latin square is an n*n grid, with each cell containing a letter from an alphabet of n letters, such that in each row and in each column, each letter appears exactly once. For example:

Two latin squares are orthogonal if, when you superimpose them, each of the n^2 possible letter-pairs appears exactly once in the grid. For example:

1a 2b 3c 4d
2d 1c 4b 3a
3b 4a 1d 2c
4c 3d 2a 1b

It turns out that the existence of mutually orthogonal n*n latin squares is related to the existence of projective planes of order n. Euler's thirty six officers problem asks for a pair of orthogonal latin squares of order 6. It turns out that there is no solution, and this is related to the fact that there is no projective plane of order 6.

(There are projective planes of order q for any prime power q = p^n, and these give rise to mutually orthogonal q*q latin squares. However 6 is not a prime power. It is not known whether there exist projective planes whose order is not a prime power, but it is known by exhaustive computer search that there is no projective plane of order 6, and no pair of mutually orthogonal latin squares of order 6.)

This module lets you construct mutually orthogonal latin squares from projective planes, and other nice stuff.

Simple groups, the atoms of symmetry

2009-10-26T20:36:00.001Z

[New release, HaskellForMaths-0.1.9, available here. Contains documentation improvements, new version of graphAuts using equitable partitions, and the code used in this blog post.]

Over the last few months on this blog, we have been looking at symmetry. We started out by looking at symmetries of graphs, then more recently we've been looking at symmetries of finite geometries. Last time, I said that there was more to say about finite geometries. There is, but first of all I realised that I need to say a bit more about groups.

We've come across groups already many times. Whenever we have an object, then the collection of all symmetries of that object is called a group. (Recall that a symmetry is a change that leaves the object looking the same.) Suppose that we call this group G, and arbitrary elements in the group g and h. Then G satisfies the following properties:
- There is a binary operation, *, defined on G. (g*h is the symmetry that you get if you do g then h.)
- If g and h are in G, then so is g*h. (Doing one symmetry followed by another is again a symmetry.)
- There is an element 1 in G, the identity, such that 1*g = g = g*1. (1 is the symmetry "do nothing".)
- If g is in G, then there is an element g^-1 in G, called the inverse of g, such that g*g^-1 = g^-1*g = 1. (Doing g^-1 is the same as undoing g.)
- (g*h)*k = g*(h*k)

Mathematicians sometimes (usually) turn things the other way round: They define a group by the properties above, and then observe that the collection of symmetries of an object fits the definition. The problem with this is that it covers over the original intuitions behind the concept.

For us, the key point to take away from the definition is that groups are "closed". We've been looking at symmetries, represented by permutations. But symmetries and permutations just are the sort of thing that you can multiply (do one then another), that have an identity (do nothing), and that have inverses (undo). Those all come for free. So when we say that a collection of symmetries, or of permutations, is a group, the only new thing that we're claiming is that the collection is closed - for g, h in G, g*h, g^-1, and 1 are also in G.

Okay, so now we know what a group is. Now you know what mathematicians are like - as soon as they've discovered some new type of structure, they always ask themselves questions like the following:
- can we classify all structures of this type?
- can smaller structures of this type be combined into larger structures of the same type?
- can larger structures of this type be broken down into smaller structures?
- can we classify all "atomic" structures of this type, that can't be broken down any further?

What about groups? Well, the answer is, for finite groups, mostly yes. Specifically, groups are made out of "atoms" - but there isn't a complete understanding of all the different ways they can be built from the atoms. Anyway, for the moment I'm going to skip the first two questions (about building groups up), and concentrate on the last two (breaking them down).

Okay, so can groups be broken down into smaller groups? Well, groups can have smaller groups contained within them. For example, let's consider the group D10 - the symmetries of the pentagon. (Somewhat confusingly, it's called D10 rather than D5, because it has 10 elements.)

> :load Math.Algebra.Group.PermutationGroup
> mapM_ print $ elts $ _D 10
[]
[[1,2],[3,5]]
[[1,2,3,4,5]]
[[1,3,5,2,4]]
[[1,3],[4,5]]
[[1,4],[2,3]]
[[1,4,2,5,3]]
[[1,5,4,3,2]]
[[1,5],[2,4]]
[[2,5],[3,4]]

Now, any subset of these elements, which is closed, is again a group. For example, the set { [], [[1,2],[3,5]] } is closed - 1, and all products and inverses of elements in the set, are in the set. This is called a subgroup of D10.

Here's Haskell code to find all subgroups of a group. (The group is passed in as a list of generators, and the subgroups are returned as lists of generators.)

> subgps gs = [] : subgps' S.empty [] (map (:[]) hs) where
>     hs = filter isMinimal $ elts gs
>     subgps' found ls (r:rs) =
>         let ks = elts r in
>         if ks `S.member` found
>         then subgps' found ls rs
>         else r : subgps' (S.insert ks found) (r:ls) rs
>     subgps' found [] [] = []
>     subgps' found ls [] = subgps' found [] [l ++ [h] | l <- reverse ls, h <- hs, last l < h]

> -- g is the minimal elt in the cyclic subgp it generates
> isMinimal 1 = False
> isMinimal g = all (g <=) primitives -- g == minimum primitives
>     where powers = takeWhile (/=1) $ tail $ iterate (*g) 1
>           n = orderElt g -- == length powers + 1
>           primitives = filter (\h -> orderElt h == n) powers

By the way, most of the code we'll be looking at this week, including the above, should only be used with small groups.

Okay, so let's try it out:

> mapM_ print $ subgps $ _D 10
[]
[[[1,2],[3,5]]]
[[[1,2,3,4,5]]]
[[[1,3],[4,5]]]
[[[1,4],[2,3]]]
[[[1,5],[2,4]]]
[[[2,5],[3,4]]]
[[[1,2],[3,5]],[[1,2,3,4,5]]]

The last in the list is D10 itself, generated by a reflection and a rotation. Then there are five subgroups generated by five different reflections. (For example, recall that [[2,5],[3,4]] is the reflection that swaps 2 with 5, and 3 with 4 - in other words, reflection in the vertical axis.). Then there is also a subgroup generated by the rotation [[1,2,3,4,5]] - the clockwise rotation that moves 1 to 2, 2 to 3, 3 to 4, 4 to 5, and 5 to 1. Finally, [] is the subgroup with no generators, which by convention is the group consisting of just the single element 1.

So can we break a group down into its subgroups? Well, not quite. When we break something down, we should end up with two (or more) parts. Subgroups are parts. But we may not always be able to find another matching part to make up the group.

I'm probably being a bit cryptic. Let's consider an analogy. We know that 15 = 3*5. We also have 5 = 15/3. So 15 = 3 * (15/3). Given a group G, and a subgroup H, we would like to be able to build G as G = H * (G/H). I'm not going to dwell on the * in this statement. What about the / ?

Well, it turns out that we can form a quotient G/H of a group by a subgroup - but only if the subgroup satisfies some extra conditions.

Given any subsets X and Y of G, we can define a product XY:

> xs -*- ys = toListSet [x*y | x <- xs, y <- ys]

For X or Y consisting of just a single element, we can form the products xY = {x}Y, or Xy = X{y}:

> xs -*  y  = L.sort [x*y | x <- xs] -- == xs -*- [y]
> x   *- ys = L.sort [x*y | y <- ys] -- == [x] -*- ys

Now it turns out that we can use this multiplication to define a quotient G/H. We let the elements of G/H be the (right) "cosets" Hg, for g in G. For example:

> let hs = [1, p [[1,2],[3,5]]]
> let x = p [[1,2,3,4,5]]
> let y = p [[2,5],[3,4]]
> hs -* x
[[[1,2,3,4,5]],[[1,3],[4,5]]]
> hs -* y
[[[1,5,4,3,2]],[[2,5],[3,4]]]

Here we have taken a subgroup H in D10, and then found the cosets Hx, Hy, for a couple of elements x, y in D10.

Now, to get our quotient group working, what we'd like is to have (Hx) (Hy) = H(xy). Unfortunately, this isn't guaranteed to be true. For example:

> (hs -* x) -*- (hs -* y)
[[],[[1,2],[3,5]],[[1,4,2,5,3]],[[1,5],[2,4]]]
> hs -* (x*y)
[[[1,4,2,5,3]],[[1,5],[2,4]]]

Oh dear. But there are some subgroups for which this construction does work. Suppose that we had a subgroup K, such that for all g in G, g^-1 K g = K. Then Kx Ky = Kx (x^-1Kx)y = KKxy = Kxy. (KK=K follows from the fact that K is a subgroup, hence closed.)

A subgroup K of G such that g^-1 K g = K for all g in G is called a normal subgroup. For a normal subgroup, we can form the quotient G/K, and so we have broken G down into two parts, K and G/K.

Here's the code:

> isNormal gs ks = all (== ks') [ (g^-1) *- ks' -* g | g <- gs]
>     where ks' = elts ks

> normalSubgps gs = filter (isNormal gs) (subgps gs)

> quotientGp gs ks
>     | ks `isNormal` gs = gens $ toSn [action cosetsK (-* g) | g <- gs]
>     | otherwise = error "quotientGp: not well defined unless ks normal in gs"
>     where cosetsK = cosets gs ks

> gs // ks = quotientGp gs ks

For example:

> mapM_ print $ normalSubgps $ _D 10
[]
[[[1,2,3,4,5]]]
[[[1,2],[3,5]],[[1,2,3,4,5]]]

Two of these are "trivial". [] is the subgroup {1}, which will always be normal. The last subgroup listed is D10 itself. A group will always be a normal subgroup of itself. So the only "proper" normal subgroup is [ [[1,2,3,4,5]] ] - the subgroup of rotations of the pentagon.

In this case, there are only two cosets:

> mapM_ print $ cosets (_D 10) [p [[1,2,3,4,5]]]
[[],[[1,2,3,4,5]],[[1,3,5,2,4]],[[1,4,2,5,3]],[[1,5,4,3,2]]]
[[[1,2],[3,5]],[[1,3],[4,5]],[[1,4],[2,3]],[[1,5],[2,4]],[[2,5],[3,4]]]

Not surprisingly then, the quotient group is rather simple:

> _D 10 // [p [[1,2,3,4,5]]]
[[[1,2]]]

Just to explain a little: In quotientGp gs ks, we form the cosets of ks, and then we look at the action of the gs on the cosets by right multiplication - g sends Kh to Khg. So the elements of the quotient group are permutations of the cosets. However, if we printed this out, it would be hard to see what was going on. So we call "toSn", which labels the cosets 1,2,..., and rewrites the elements as permutations of the numbers. In the case above, we see that the quotient group can be generated by a single element, which simply swaps the two cosets.

Let's look at another example. S 4 is the group of all permutations of [1..4] - which also happens to be the symmetry group of the tetrahedron. S 4 has lots of subgroups:

> mapM_ print $ subgps $ _S 4
[]
[[[1,2]]]
[[[1,2],[3,4]]]
[[[1,2,3]]]
[[[1,2,3,4]]]
[[[1,2,4,3]]]
[[[1,2,4]]]
[[[1,3],[2,4]]]
[[[1,3,2,4]]]
[[[1,3]]]
[[[1,3,4]]]
[[[1,4],[2,3]]]
[[[1,4]]]
[[[2,3]]]
[[[2,3,4]]]
[[[2,4]]]
[[[3,4]]]
[[[1,2]],[[1,2],[3,4]]]
[[[1,2]],[[1,2,3]]]
[[[1,2]],[[1,2,3,4]]]
[[[1,2]],[[1,2,4]]]
[[[1,2]],[[1,3],[2,4]]]
[[[1,2],[3,4]],[[1,2,3]]]
[[[1,2],[3,4]],[[1,2,3,4]]]
[[[1,2],[3,4]],[[1,2,4,3]]]
[[[1,2],[3,4]],[[1,3],[2,4]]]
[[[1,3],[2,4]],[[1,3]]]
[[[1,3]],[[1,3,4]]]
[[[1,4],[2,3]],[[1,4]]]
[[[2,3]],[[2,3,4]]]

However, only a few of them are normal:

> mapM_ print $ normalSubgps $ _S 4
[]
[[[1,2]],[[1,2,3,4]]]
[[[1,2],[3,4]],[[1,2,3]]]
[[[1,2],[3,4]],[[1,3],[2,4]]]

The first two are 1 and S4 itself. The third is the group of rotations of the tetrahedron. The fourth is just the rotations which move all the points.

Why are these subgroups normal, but the others aren't? Well, let's have a look at some of the ones that aren't.

We have a subgroup [ [[1,2]] ], and another [ [[1,3]] ], and several others that "look the same" - they just swap two points. We also have a subgroup [ [[1,2,3]] ], a subgroup [ [[1,2,4]] ], and several others that "look the same" - they just just rotate three points while leaving the fourth fixed. In fact, in both cases we have a set of conjugate subgroups. What I mean is, for a subgroup to be normal, we required that g^-1 K g = K. If a subgroup is not normal, then we can find a g such that g^-1 H g /= H. Then g^-1 H g will also be a subgroup (exercise: why?), and it will "look the same" as H.

By contrast, a normal subgroup is often the only subgroup that "looks like that".

Okay, so we've seen how groups can be broken down into smaller groups. Now, what about "atoms"? Well, a group always has itself and 1 as normal subgroups. If it has other normal subgroups, then we can break it down as K, G/K. We can keep going, and see if either of these in their turn has non-trivial normal subgroups. Eventually, we will end up with groups which have no non-trivial normal subgroups. Such a group is called a simple group. This is a bit like a prime number - which has no other factors besides itself and 1. Mark Ronan, in Symmetry and the Monster, calls simple groups the "atoms of symmetry". They are the pieces out of which all symmetry groups are made.

(Note that starting from a group G, with normal subgroups K1, K2, ..., we have more than one choice of how to break it down. Luckily, it turns out that if we keep going, then at the end we will have the same collection of "atoms", no matter which order we chose.)

So, what do simple groups look like?

Well, first, here's some code to detect them:

> isSimple gs = length (normalSubgps gs) == 2

Then the finite simple groups are as follows:

- The cyclic groups of order p, p prime (the rotation group of a regular p-gon):

> _C n | n >= 2 = [p [[1..n]]]

For example:

> isSimple $ _C 5
True
> isSimple $ _C 6
False

Exercise: Why is C n simple when n is prime, and not when it is not?

- The alternating groups A n, n >= 5. (A n is the subgroup of S n consisting of those elements of S n which are "even" - see wikipedia.)

> isSimple $ _A 4
False
> isSimple $ _A 5
True

(As it happens, A4 is the group of rotations of the tetrahedron, and A5 is the group of rotations of the dodecahedron. Unfortunately, An, n>5, doesn't have any such intuitive interpretation.)

- The projective special linear groups PSL(n,Fq) (except for PSL(2,F2), PSL(2,F3) and PSL(3,F2)). This is the subgroup of PGL(n,Fq) consisting of projective transformations with determinant 1. (Recall that last time we looked at PΓL(n,Fq), the group of symmetries of the projective geometry PG(n-1,Fq), which includes both projective transformations and field automorphisms.)

- Several more infinite families of subgroups of PΓL(n,Fq), with geometric significance.

- Finally, there are 26 "sporadic" simple groups, which don't belong to any of the infinite families described above.

I should just say a thing or two about the classification of finite simple groups:
- The proof of the classification was a monster effort by a whole generation of group theorists. The proof runs to 10000 pages.
- The part that mathematicians find most interesting is the sporadic simple groups. In many mathematical classifications, one only finds infinite families (for example, the primes), so the sporadic groups feel like little miracles that have dropped out of the sky. Why are they there?

Two of the things I'm hoping to do in future blog posts is describe and construct some of the other "simple groups of Lie type" (subgroups of PΓL(n,Fq), and describe and construct some of the sporadic simple groups.

Symmetries of PG(n,Fq)

2009-10-08T19:55:00.000+01:00

Previously in this blog we looked at the affine geometries AG(n,Fq), and their symmetries. Following that, we looked at the points and lines in the projective geometries PG(n,Fq). This week I want to look at the symmetries of PG(n,Fq).

However, first, something I forgot to mention last week. When we looked at the points in PG(n,Fq), we saw that they can be thought of as the points of AG(n,Fq), plus some additional points "at infinity". What I forgot to do last week was to discuss how the lines in PG(n,Fq) relate to the lines in AG(n,Fq).

Here's PG(2,F3) again:

Recall that the points in PG(2,F3) are (by definition) the lines through the origin in F3^3, each of which can be represented, as here, by one of the non-zero points on it. (Note that for technical reasons, the coordinates are shown as zyx instead of xyz.) We see that these "points" of PG(2,F3) consist of something looking very much like AG(2,F3) - the blue points - together with some additional points, which are said to be "at infinity".

Then, the lines of PG(2,F3) are (by definition) the planes through 0 in F3^3. What is the relationship between lines in AG(2,F3) and "lines" in PG(2,F3)? Well, any line in (the copy of) AG(2,F3) (embedded in PG(2,F3)) is contained in a plane through 0 in F3^3. So each line in AG(2,F3) gives rise to a line in PG(2,F3). In PG(2,F3), however, the line will have an additional point "at infinity".

For example, there is a line in AG(2,F3) consisting of the points {(0,0),(1,0),(2,0)}. These embed into PG(2,F3) as {(1:0:0), (1:1:0), (1:2:0)}. We can ask what is the closure in PG(2,F3) of these points - the least "flat" (point, line, plane, etc) that contains them:

> closurePG [ [1,0,0],[1,1,0],[1,2,0] ] :: [[F3]]

[[0,1,0],[1,0,0],[1,1,0],[1,2,0]]

So in PG(2,F3), the line gains a point at infinity, (0:1:0).

Similarly, there is a line in AG(2,F3) consisting of the points {(0,1),(1,1),(2,1)}. These embed into PG(2,F3) as {(1:0:1), (1:1:1), (1:2:1)}.

> closurePG [ [1,0,1],[1,1,1],[1,2,1] ] :: [[F3]]

[[0,1,0],[1,0,1],[1,1,1],[1,2,1]]

Notice that both lines from AG(2,F3) gained the same point at infinity. This is because the two lines were parallel. In PG(2,F3), "parallel lines meet at infinity".

So each line in AG(2,F3) extends to a line in PG(2,F3), gaining an extra point at infinity. In addition to these lines, there is one more line in PG(2,F3) - a line at infinity consisting of the points {(0:0:1), (0:1:0), (0:1:1), (0:1:2)}, corresponding to the plane z = 0 in F3^3.

So, that was all stuff that I should really have talked about last time, but I forgot.

--*--

Okay, so back to the symmetries of PG(n,Fq). So we consider PG(n,Fq) as an incidence structure between points and lines, and we ask which permutations of the points preserve collinearity. That is, given a permutation g of the points of PG(n,Fq), we say that it is a symmetry (or automorphism, or collineation) of PG(n,Fq) if, whenever a and b are collinear points, then so are a^g and b^g, and vice versa.

We can find the symmetries of PG(n,Fq) by forming its incidence graph. This is the bipartite graph having as vertices: the points of PG(n,Fq) on the left (or in blue), the lines of PG(n,Fq) on the right (or in red), with an edge joining a point vertex to a line vertex if the point is incident with the line.

Here's the code:

incidenceGraphPG n fq = G vs es where

points = ptsPG n fq

lines = linesPG n fq

vs = L.sort $ map Left points ++ map Right lines

es = L.sort [ [Left x, Right b] | b <- lines, x <- closurePG b]

Let's look at an example. Here's PG(2,F2):

The lines in PG(2,F2), corresponding to planes through 0 in F2^3, are:

> mapM_ (print . closurePG) $ linesPG 2 f2

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

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

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

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

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

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

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

PG(2,F2) is called the Fano plane, and is often represented by the picture below:

You can easily check that this is right. Notice the embedded AG(2,F2) in the bottom left, with parallel lines meeting "at infinity" as expected.

The incidence graph of the Fano plane is called the Heawood graph, and it is often represented like this:

The blue vertices are the points of the Fano plane, the red vertices are the lines. For each red vertex, you can check that the three blue vertices it is connected to are indeed a line of the Fano plane.

Then we can find the automorphisms of the Fano plane as follows:

> let heawood = incidenceGraphPG 2 f2

> let auts = incidenceAuts heawood

> mapM_ print auts

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

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

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

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

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

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

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

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

> orderSGS auts

168

Now the symmetries of PG(n,Fq) actually have a very simple description in terms of matrices. Consider the invertible matrices of degree n+1 over Fq - the general linear group GL(n+1,Fq). These matrices act on vectors in Fq^(n+1) by multiplication, and in doing so, they permute the lines and planes through 0 in Fq^(n+1), whilst preserving containment of lines within planes. In other words, they permute the points and lines of PG(n,Fq) and preserve collinearity. However, any scalar matrix within GL(n+1,Fq) - that is, kI, where k is in Fq\{0} and I is the identity - acts as the identity on PG(n,Fq) - since if (x0:x1:...:xn) is a "point" of PG(n,Fq), then (kx0:kx1:...:kxn) represents the same point. So the group of projective transformations of PG(n,Fq) is actually GL(n+1,Fq) factored out by scalar multiplications. This group is called the projective general linear group, PGL(n+1,Fq).

If we are given a symmetry of PG(n,Fq), expressed as a permutation, then we can express it as a matrix by looking at what it does to the basis vectors. For example, consider the symmetry [[[0,0,1],[0,1,0]],[[1,0,1],[1,1,0]]] of PG(2,F2) from above. It sends [1,0,0] to [1,0,0], [0,1,0] to [0,0,1], and [0,0,1] to [0,1,0]. So it can be represented by the matrix

[1 0 0]

[0 0 1]

[0 1 0]

We already worked out the order of GL(n,Fq) in the orderGL function. There are q-1 non-zero scalars to factor out by. So the order of PGL(n,Fq) will be given by the following:

orderPGL n q = orderGL n q `div` (q-1)

Let's test it (remembering that we want PGL(n+1,Fq) for PG(n,Fq)):

> orderPGL 3 2

168

> orderSGS $ incidenceAuts $ incidenceGraphPG 2 f3

5616

> orderPGL 3 3

5616

> orderSGS $ incidenceAuts $ incidenceGraphPG 2 f4

120960

> orderPGL 3 4

60480

With F4, we get twice as many automorphisms as we might have expected. If you've been paying attention, you'll remember that this is because F4 has a field automorphism (the Frobenius automorphism) of order 2. The extended group, PGL(n,Fq) plus field automorphisms, is called PGammaL(n,Fq) (or better - but I'm not sure this will work for everybody - PΓL(n,Fq)).

Now, this is kind of a trivial point, but I think it's significant: PG(n,Fq) has more symmetries than AG(n,Fq). Somehow, by adding the points at infinity, PG(n,Fq) has made AG(n,Fq) "whole", making it more symmetrical.

Anyway, that'll do for now. Next week, just a few more odds and ends about finite geometries.

Finite geometries, part 4: Lines in PG(n,Fq)

2009-09-30T20:15:00.000+01:00

[Apologies about the formatting - blogger seems to have won the struggle this time.]

Last time we saw that the points in the projective geometry PG(n,Fq) are defined to be the lines (through the origin) in the vector space Fq^(n+1). What about lines in PG(n,Fq)? Well, it's simple - they are defined as the planes (through 0) in Fq^(n+1).

Let's look at an example, to try to clarify our intuitions. Here are the points of PG(2,F3) again.

> :load Math.Combinatorics.FiniteGeometry

> ptsPG 2 f3

[[0,0,1],[0,1,0],[0,1,1],[0,1,2],[1,0,0],[1,0,1],[1,0,2],[1,1,0],[1,1,1],[1,1,2],[1,2,0],[1,2,1],[1,2,2]]

(Recall that each point in PG(2,F3) represents a line through 0 in F₃³. We have chosen the representative which has 1 as its first non-zero coordinate. This has the consequence that the coordinates in the diagram are best read back-to-front - they're zyx instead of xyz.)

So the lines in PG(2,F3) correspond to the planes through 0 in F₃³. For example, the plane x = 0 (the left hand face of the cube) gives rise to a line in PG(2,F3) containing the points 010, 100, 110, 120. The plane x = z (a diagonal cut through the cube from front left to back right) gives rise to the line containing the points 010, 101, 111, 121. Just to be clear, we are only looking at the planes in F₃³ that contain 0. So for example, 012, 101, 111, 121 - the plane x+z=2 - is not a line in PG(2,F3). Also, remember that F₃ wraps around - so for example 010, 102, 112, 122 is a line in PG(2,F3).

Okay, onwards. As we did with AG(n,Fq), we can represent a line in PG(n,Fq) by any pair of points on it. We would then like to be able to do the following:

Given a line of PG(n,Fq), list the points of PG(n,Fq) which are on it
Given a line of PG(n,Fq), find a canonical representation for it
Find all lines of PG(n,Fq)

(These are basically just the same things that we did when we looked at AG(n,Fq) a couple of weeks back.)

Okay, so suppose that we are given two distinct points of PG(n,Fq), representing a line, and we are asked to list all the points of PG(n,Fq) that are on the line. Well, the two points are also points in F_qⁿ⁺¹, and as they are distinct in PG(n,Fq), they are linearly independent in F_qⁿ⁺¹ - hence they generate a plane through 0. We are being asked to find all lines through 0 that are contained in the plane. Well, that's easy:

Generate all linear combinations of the two points - all points in the plane in F_qⁿ⁺¹.
Each non-zero point in the plane represents a line through 0, and hence a point of PG(n,Fq)
However, we can filter out all except those that are in the canonical form for pts of PG(n,Fq) - having a 1 as their first non-zero coordinate - in order to avoid double counting the lines.

So we need an auxiliary function to detect whether a point is in canonical form - what I call "projective normal form" or PNF:

ispnf (0:xs) = ispnf xs

ispnf (1:xs) = True

ispnf _ = False

Then, given two distinct points, to list all points on the line they generate:


linePG [p1,p2] = toListSet $ filter ispnf [(a *> p1) <+> (b *> p2) | a <- fq, b <- fq]
    where fq = eltsFq undefined

(Recall that c *> v is scalar multiplication of a vector, and u <+> v is addition of vectors. Note that in HaskellForMaths <= 0.1.9, you must use the more general "closurePG" function instead, which has the same behaviour.)

For example, here are the two lines of PG(2,F3) that we discussed earlier:

> linePG [[0,1,0],[1,0,0]] :: [[F3]]

[[0,1,0],[1,0,0],[1,1,0],[1,2,0]]

> linePG [[0,1,0],[1,0,1]] :: [[F3]]

[[0,1,0],[1,0,1],[1,1,1],[1,2,1]]

A line in PG(n,Fq) can be represented by any pair of distinct points on the line. Is there a canonical pair to pick? Yes, there is. Given any pair of points of PG(n,Fq), we can consider them as the rows of a matrix. Any matrix which can be obtained from this matrix by elementary row operations represents the same plane in F_qⁿ⁺¹. For our canonical form, we can take the reduced row echelon form for the matrix (see wikipedia). However, I'm not going to dwell on this - use the "reducedRowEchelonForm" function if you'd like to experiment.

What I'd like to consider instead is how we can list all the lines in PG(n,Fq). The answer is simply to list all the reduced row echelon forms. We can do this in two steps:

List all the reduced row echelon "shapes"
Fill in all possible combinations of values from Fq, to get all reduced row echelon forms

An example will make this clearer. Suppose we want to find the lines in PG(3,Fq). Thus we are looking for the 2-dimensional subspaces of Fq^4. The reduced row echelon shapes are the following:


[1 0 * *]                   [1 * 0 *]                   [1 * * 0]
[0 1 * *]                   [0 0 1 *]                   [0 0 0 1]

[0 1 0 *]                   [0 1 * 0]                   [0 0 1 0]
[0 0 1 *]                   [0 0 0 1]                   [0 0 0 1]

To get the reduced row echelon forms, we take each shape in turn, and fill in the stars with the values from Fq in all possible ways.

The Haskell code finds row echelon shapes for k-dimensional subspaces of Fq^n. (For lines, set k = 2):


data ZeroOneStar = Zero | One | Star deriving (Eq)

instance Show ZeroOneStar where
   show Zero = "0"
   show One  = "1"
   show Star = "*"

rrefs n k = map (rref 1 1) (combinationsOf k [1..n]) where
   rref r c (x:xs) =
       if c == x
       then zipWith (:) (oneColumn r) (rref (r+1) (c+1) xs)
       else zipWith (:) (starColumn r) (rref r (c+1) (x:xs))
   rref _ c [] = replicate k (replicate (n+1-c) Star)
   oneColumn r = replicate (r-1) Zero ++ One : replicate (k-r) Zero
   starColumn r = replicate (r-1) Star ++ replicate (k+1-r) Zero

Thus we can calculate the shapes we saw above as follows:


> mapM_ print $ rrefs 4 2
[[1,0,*,*],[0,1,*,*]]
[[1,*,0,*],[0,0,1,*]]
[[1,*,*,0],[0,0,0,1]]
[[0,1,0,*],[0,0,1,*]]
[[0,1,*,0],[0,0,0,1]]
[[0,0,1,0],[0,0,0,1]]

Next, to find all lines in PG(n,Fq), we need to substitute values from Fq for the stars. The following code does the trick, although I suspect that there might be a better way to write it:


flatsPG n fq k = concatMap substStars $ rrefs (n+1) (k+1) where
   substStars (r:rs) = [r':rs' | r' <- substStars' r, rs' <- substStars rs]
   substStars [] = [[]]
   substStars' (Star:xs) = [x':xs' | x' <- fq, xs' <- substStars' xs]
   substStars' (Zero:xs) = map (0:) $ substStars' xs
   substStars' (One:xs) = map (1:) $ substStars' xs
   substStars' [] = [[]]

linesPG n fq = flatsPG n fq 1

For example:


> mapM_ print $ linesPG 2 f3
[[1,0,0],[0,1,0]]
[[1,0,0],[0,1,1]]
[[1,0,0],[0,1,2]]
[[1,0,1],[0,1,0]]
[[1,0,1],[0,1,1]]
[[1,0,1],[0,1,2]]
[[1,0,2],[0,1,0]]
[[1,0,2],[0,1,1]]
[[1,0,2],[0,1,2]]
[[1,0,0],[0,0,1]]
[[1,1,0],[0,0,1]]
[[1,2,0],[0,0,1]]
[[0,1,0],[0,0,1]]

It occurs to me that I should explain why PG(n,Fq) is worth looking at. Well, I mentioned before that the symmetry groups of PG(n,Fq) are some of the "atoms of symmetry", from which all symmetry groups are composed. In addition, projective geometry is in fact more fundamental than affine geometry (at least from the point of view of algebraic geometry, although I'm not sure how you would justify that statement in general). Well, maybe you'll just have to take my word for it that PG(n,Fq) is a beautiful thing, for the moment.

Next time, symmetries of PG(n,Fq).

Finite geometries, part 3: Points in PG(n,Fq)

2009-09-25T09:00:00.000+01:00

[New release HaskellForMaths-0.1.8 - contains bugfix to Fractional instance for ExtensionField, plus documentation improvements.]

Over the last two weeks, we have looked at the finite affine geometries AG(n,Fq). Hopefully, this has all been fairly straightforward - affine geometry is just the familiar Euclidean geometry of points and lines, but without angles or distances, and all we have done is replace the reals R by the finite fields Fq.

This week we're going to look at finite projective geometries. Now, I know from my own experience that it can take a little while to "get" projective geometry. There are probably several reasons for this, but one of them is that mathematicians define projective geometry one way, but then most of the time think about it in a different way. So, I'll do my best to explain it, but don't worry if you don't get it at first, one day it will all click into place.

Okay, so let's start with the points. The points of PG(n,Fq) are defined to be the lines through the origin (ie the one-dimensional subspaces) in Fq^(n+1). That should sound a bit strange at first reading - how can the points (in PG(n,Fq)) be lines (in Fq^(n+1))? Bear with me - we'll see why it makes sense to think of them as points.

For example, the points of PG(2,F3) are the lines through 0 in F3^3. So the line {(0,0,0), (0,1,0), (0,2,0)} is a point of PG(2,F3), as is the line {(0,0,0), (1,0,2), (2,0,1)}. (Recall that arithmetic in F3 is modulo 3.) As a line in F3^3 is a one-dimensional subspace, it is generated by any non-zero point on the line, and consists of all scalar multiples of such a point. So we could represent the two lines we just mentioned as <(0,1,0)> and <(1,0,2)>, with angle brackets meaning "the line generated by". Any non-zero point on the line will do to represent the line, but it would be good to have a way of choosing a canonical representative. For the moment, let's say that we will choose the point whose last non-zero coordinate is 1. Here, then, are the points of PG(2,F3):

So the blue points are the points of the form (x,y,1). Every point of the form (x,y,1) represents a line through 0 in F3^3. However, this is not all the lines through 0, for it misses out those lines which are in the plane z = 0. The green points are the points of the form (x,1,0). Each point of the form (x,1,0) represents a line through 0 in the plane z = 0. However, this is still not all, for we are missing the line y = z = 0. This is represented by the red point (1,0,0). And that's it - these are all the "points" of PG(2,F3). (Confirm for yourself that every line through zero in F3^3 is represented by one of these "points".)

Since any scalar multiple of a point in F3^3 represents the same line through 0, and hence the same point of PG(2,F3), mathematicians often write the representatives as (x:y:z), with the colons indicating that it is only the ratios we are interested in. Thus (0:1:0) and (0:2:0) represent the same line in F3^3, and hence the same point of PG(2,F3). (Note that (0:0:0) is not a point of PG(2,F3).)

Okay, so I mentioned that mathematicians define PG(n,Fq) one way, but then think about it another way. So how do mathematicians think about it?

Well, the blue points in PG(2,F3) look just like a copy of AG(2,F3), don't they. Then the green points are called "the line at infinity". Finally, the red point is called "the point at infinity". So mathematicians think of PG(2,F3) as being like AG(2,F3), but with some additional points "at infinity". In particular, although PG(2,F3) was defined in terms of lines through 0 in F3^3, mathematicians actually think of it as consisting of points (not lines).

The idea that the additional points are "at infinity" comes from perspective drawing. Imagine an artist sitting in front of a canvas, looking along the z-axis. Their eye is the origin. The canvas is the plane {(x,y,1)}. To draw what they see, the artist needs to project a line from their eye, through the canvas, until it hits something. Given a large enough canvas, the artist can draw anything that is in front. However, things that are above or to the sides would have to be drawn "at infinity".

One thing I should emphasize is that, contrary to appearance, the line and point at infinity are no different from the other points. Indeed, I could have had my artist looking along the x-axis instead of the z-axis. In that case, the embedded affine plane would have been the points {(1,y,z)}, the line at infinity would have been {(0,1,z)}, and the point at infinity would have been (0,0,1). If we revert to thinking about lines in F3^3 for a moment, it should be clear that the points at infinity are not distinguished in any way from other lines in F3^3. Their appearance of being distinguished has arisen solely from our choice of coordinate system.

Okay, time for some code. In the above, we said that the canonical representative for a line through 0 would be the point with 1 as its last non-zero coordinate - corresponding to looking along the z-axis. That made things easier to explain. However, it turns out that it's usually more convenient to choose the point with 1 as its first non-zero coordinate (ie, looking along the x-axis). Here is the code to list the points of PG(n,Fq):

ptsPG 0 _ = [[1]]
ptsPG n fq = map (0:) (ptsPG (n-1) fq) ++ map (1:) (ptsAG n fq)

For example:


> ptsPG 2 f3
[[0,0,1],[0,1,0],[0,1,1],[0,1,2],[1,0,0],[1,0,1],[1,0,2],[1,1,0],[1,1,1],[1,1,2],[1,2,0],[1,2,1],[1,2,2]]
> map reverse it
[[1,0,0],[0,1,0],[1,1,0],[2,1,0],[0,0,1],[1,0,1],[2,0,1],[0,1,1],[1,1,1],[2,1,1],[0,2,1],[1,2,1],[2,2,1]]

If we reverse the coordinates, you can see that we have the point at infinity, the line at infinity, and the embedded copy of AG(2,F3), as before.

That's probably enough to take in at one sitting. Next time, we'll look at lines in PG(n,Fq).

Finite geometries, part 2: Symmetries of AG(n,Fq)

2009-09-18T13:00:00.000+01:00

Last time, we met the finite affine geometries AG(n,Fq), which are analogous to n-dimensional Euclidean geometry, but defined over the finite fields Fq instead of the reals R. This time I want to talk about the symmetries of AG(n,Fq).

Recall that when we were looking at graphs, we defined a symmetry of a graph as a permutation of the vertices which left the edges (collectively) in the same places. For the moment, we will think of a finite geometry as a configuration of points and lines, and define a symmetry as a permutation of the points which leaves the lines (collectively) in the same places.

Here's AG(2,F3) again:

There's an obvious four-fold rotational symmetry, which moves all the points except the middle point, moves all the lines, but leaves the lines collectively in the same places. Less obvious perhaps is that reflection in the middle line, or translating everything one point to the right (with wraparound), are also symmetries. Remember that the straightness or otherwise of the lines doesn't matter - the only thing we're interested in is which points they are incident with.

So how can we find all the symmetries? Well, when we were looking at graphs, we used depth-first search, trying to create a pairing of source and target vertices, and backtracking whenever the adjacency between the source vertices didn't match the adjacency between the target vertices. Unfortunately, this method doesn't easily generalize to finite geometries. Exercise: Why not?

Instead, we will change the problem into a problem about graphs. Given a collection of points and lines, we can create an "incidence graph". This is the bipartite graph constructed as follows. On the left, we have a vertex for every point. On the right, we have a vertex for every line. We have an edge joining a point on the left to a line on the right just in the case that the point is incident with the line in our geometry.

Here's the code to create the incidence graph of AG(n,Fq):

incidenceGraphAG n fq = G vs es where
    points = ptsAG n fq
    lines = linesAG n fq
    vs = L.sort $ map Left points ++ map Right lines
    es = L.sort [ [Left x, Right b] | b <- lines, x <- closureAG b]

(Recall that a line is represented by two distinct points on it. "closureAG" calculates all the points on the line.)

Now, the trick is, that we can find the symmetries of the finite geometry by finding the symmetries of its incidence graph. For clearly, any symmetry of the finite geometry permutes the points among themselves, permutes the lines among themselves, and preserves incidence - hence it gives rise to a symmetry of the incidence graph. So in principle, we can find the symmetries of the finite geometry by doing the following:

> mapM_ print $ graphAuts $ incidenceGraphAG 2 f3
[[Left [0,0],Left [0,1]],[Left [1,0],Left [2,2],Left [1,1],Left [2,1],Left [1,2],Left [2,0]],[Right [[0,0],[1,0]],Right [[0,1],[1,0]],Right [[0,0],[1,1]],Right [[0,1],[1,1]],Right [[0,0],[1,2]],Right [[0,1],[1,2]]],[Right [[0,2],[1,0]],Right [[0,2],[1,2]],Right [[0,2],[1,1]]],[Right [[1,0],[1,1]],Right [[2,0],[2,1]]]]
[[Left [0,0],Left [0,2],Left [0,1]],[Left [2,0],Left [2,1],Left [2,2]],[Right [[0,0],[1,0]],Right [[0,2],[1,0]],Right [[0,1],[1,0]]],[Right [[0,0],[1,1]],Right [[0,2],[1,1]],Right [[0,1],[1,1]]],[Right [[0,0],[1,2]],Right [[0,2],[1,2]],Right [[0,1],[1,2]]]]
...

However, there are a couple of minor problems with this:

We only really want to know what happens to the points. What happens to the lines follows from that. So we only want to know about the Left parts, not the Right parts. (Furthermore, if we were to show only the Left part, we then wouldn't need all those Lefts and Rights which are cluttering up the output.)
It is just possible that the incidence graph might have some symmetries which interchange points and lines. This is a very interesting situation, which we'll discuss in due course. But it is a symmetry of the graph which does not correspond to a symmetry of the geometry.

For these reasons, HaskellForMaths provides an "incidenceAuts" function, which works just like graphAuts, except that it knows that the input is an incidence graph, so it avoids wasting time looking for point-line crossover symmetries, and outputs only the permutations of the points.

> mapM_ print $ incidenceAuts $ incidenceGraphAG 2 f3
[[[0,0],[0,1]],[[1,0],[2,2],[1,1],[2,1],[1,2],[2,0]]]
[[[0,0],[0,2],[0,1]],[[2,0],[2,1],[2,2]]]
[[[0,0],[1,0],[0,1]],[[0,2],[2,0],[2,2]],[[1,1],[2,1],[1,2]]]
[[[0,1],[0,2]],[[1,1],[1,2]],[[2,1],[2,2]]]
[[[0,1],[1,0]],[[0,2],[2,0]],[[1,2],[2,1]]]
[[[0,1],[1,1],[1,2],[2,0],[0,2],[2,2],[2,1],[1,0]]]
[[[1,0],[1,1],[1,2]],[[2,0],[2,2],[2,1]]]
[[[1,0],[2,0]],[[1,1],[2,1]],[[1,2],[2,2]]]

This is a strong generating set for the symmetries of AG(2,F3). We could now go on to investigate further. We could ask what different types of symmetries there are, using the conjClassReps function, or how many symmetries there are, using the orderSGS function.

For example, among the elements of the SGS, some have obvious interpretations:

[[[1,0],[2,0]],[[1,1],[2,1]],[[1,2],[2,2]] - remember that this means that it swaps [1,0] and [2,0], swaps [1,1] and [2,1], and swaps [1,2] and [2,2]. So it can be thought of as a reflection in the line x = 0, or as a 2* stretch of the x-coordinate. It is given by the matrix [[2,0],[0,1]]. (In F3, 2 = -1, which is why this can be thought of as either a reflection or a stretch.)

[[[1,0],[1,1],[1,2]],[[2,0],[2,2],[2,1]]] is the shear given by the matrix [[1,1],[0,1]].

[Later: You might question whether a shear is a symmetry. We have defined a symmetry informally as a change that leaves things looking the same. But, for example, a shear in R^2 doesn't leave the unit square looking the same. Is it really a symmetry? Well, when we're studying geometries within combinatorics, we're only interested in the incidence between points and lines, and not, for example, distances. From this point of view, a shear is indeed a symmetry, because it preserves incidence. We will later look at what happens when we require symmetries to preserve more than just incidence.]

[[[0,0],[1,0],[0,1]],[[0,2],[2,0],[2,2]],[[1,1],[2,1],[1,2]]] cannot be represented by a 2*2 matrix, since it moves the origin. However, it can be represented by a 3*3 matrix, as follows:

[x']   [2 2 1] [x]
[y'] = [1 0 0] [y]
[ 1]   [0 0 1] [1]

This is the same as saying:

[x'] = [2 2] [x] + [1]

[y'] [1 0] [y] [0]

So we see that it is a composite of some sort of shear, followed by a translation.

Indeed, we could conjecture that all symmetries of AG(n,Fq) are of this form - a linear transformation followed by a translation. Let's see whether we can confirm this conjecture by counting symmetries.

We have seen that some symmetries can be represented by a 2*2 matrix. In fact, the matrix must be non-singular, otherwise we won't have a permutation. The group of 2*2 non-singular matrices over Fq is the general linear group GL(2,Fq). How many elements does it have? Well, for the first row, we can choose any non-zero vector in Fq^2, of which there are q^2-1. For the second row, we can choose any vector in Fq^2 which is not linearly dependent on the one we already chose, of which there are q^2-q. Generalising to GL(n,Fq), the following code calculates the number of elements of GL(n,Fq):

orderGL n q = product [q^n - q^i | i <- [0..n-1] ]

Then, in addition to these linear transformations, all of which leave the origin fixed, we can do a translation. The number of translations of Fq^n is of course q^n.

A transformation consisting of a translation and a linear transformation is called an affine transformation. The group of all such is called the affine group, Aff(n,Fq). It is straightforward to calculate the number of elements of this group:

orderAff n q = q^n * orderGL n q

(Note that by multiplying the number of translations by the number of linear transformations, I have made an assumption that they are "semi-independent" of one another. This assumption is valid, but at this stage I don't want to spell out exactly what it means.)

Okay, so our conjecture is that all symmetries of AG(n,Fq) are affine transformations. Let's test our conjecture:

> orderSGS $ incidenceAuts $ incidenceGraphAG 2 f3
432
> orderAff 2 3
432
> orderSGS $ incidenceAuts $ incidenceGraphAG 2 f4
5760
> orderAff 2 4
2880

Uh-oh - what's going on here? We seem to have twice as many symmetries as we expected.

It took me quite a while to figure this one out actually. What's going on?

Well, remember that when we looked at extension fields, I mentioned that there can be field automorphisms. (For example, in the complex numbers, we have complex conjugation.) When we looked at the finite fields Fq, for q = p^n a prime power, I mentioned the Frobenius automorphism x -> x^p. Okay, well what's going on is that the field automorphisms of Fq give rise to further symmetries of AG(n,Fq). If q = p^n, then the number of such automorphisms is n, so the total number of symmetries of AG(n,Fq) will be n * orderAff n q.

And now that really does account for all symmetries of AG(n,Fq). Apart from the following exceptions:

> orderSGS $ incidenceAuts $ incidenceGraphAG 3 f2
40320
> orderAff 3 2
1344

What's going on here? Well the problem is that AG(n,F2) is degenerate. A line in AG(n,F2) has just two points - and every pair of points forms a line. So in fact, AG(n,F2) is just the complete graph K (2^n), and hence has (2^n)! symmetries.

That's it for now. Next time, finite projective geometries.

Finite geometries, part 1: AG(n,Fq)

2009-09-13T14:00:00.000+01:00

Time to recap a little. Through June and July, we looked at graphs and their symmetries. Then through August, we looked at finite fields. Now I want to look at finite geometries. Why?

Well, it's a continuation of my earlier investigation of symmetry. I chose graphs as the starting point because they're easy to understand. However, graphs are not the only combinatorial structure having symmetry - far from it. Finite geometries are important for several reasons:
1. We'll see that all symmetry groups turn out to be composed out of "atoms of symmetry", called "simple groups". There are several infinite families of finite simple groups - and a few "sporadic" ones that don't fit into any family. One of the infinite families are symmetries of finite projective geometries, which we will look at in due course. (I got the term "atoms of symmetry" from Ronan, Symmetry and the Monster, which is a nice easy read on this stuff.)
2. The symmetry groups of infinite geometries are very important in physics (especially particle physics). So if you like, you can pretend that what we're looking at here will help you understand physics.

This time, I want to look at affine geometries, but what we're aiming towards is projective geometries.

Within mathematics, there are probably many different definitions of geometry. Within combinatorics, we're just interested in the combinatorial structure, so roughly speaking, a geometry is just a collection of points, lines, planes and so on, together with a relation, called "incidence", saying which points are on which lines, and so on. Most of the time it's natural to think of a line as just being a set of points, in which case we can think of the incidence relation as just being the membership relation. However, this isn't always the best way to represent things in the code (or in the maths), as we shall see.

The most familiar geometry is the Euclidean plane, R^2. In this geometry, every pair of lines meets in at most one point. Then of course, there are Euclidean spaces of three or more dimensions. There are also non-Euclidean geometries. For example, spherical geometry is geometry on the surface of a sphere, with the lines being great circles. In this case, every pair of lines meets in exactly two points.

If you know about vector spaces, then you might think that Euclidean geometry is basically about the vector space R^n, with points, lines, planes being the zero-, one-, two-dimensional subspaces, and so on. However, this isn't quite right, because a vector space has a distinguished point, the origin, and every subspace of a vector space contains the origin. In Euclidean geometry there is no such distinguished point - we can have points, lines, planes which don't contain the origin. So we need to take the vector space R^n, and then "forget" the origin. If we do this, we get what is called an affine space. The subspaces of an affine space are the subspaces of the underlying vector space and their translates. The collection of points, lines, planes, and so on in an affine space, together with the incidence relation, is called an affine geometry.

We can construct finite affine geometries simply by replacing R by a finite field Fq in the above. In this way, we obtain the affine geometries AG(n,Fq). So in AG(n,Fq), we will have points, lines, planes and so on, but there will be only a finite number of each. Let's start with the points. Well, that's easy. The points of AG(n,Fq) are just the points Fq^n. The HaskellForMaths code goes like this:


ptsAG 0 fq = [[]]
ptsAG n fq = [x:xs | x <- fq, xs <- ptsAG (n-1) fq]

Actually, I think the stylish way to do this in Haskell would be:


ptsAG n fq = sequence $ replicate n fq

In any case, here's an example:


> :load Math.Combinatorics.FiniteGeometry
> ptsAG 2 f3
[[0,0],[0,1],[0,2],[1,0],[1,1],[1,2],[2,0],[2,1],[2,2]]

In other words, the points of AG(2,F3) form a 3*3 grid.

What about the lines? Well, the points of AG(2,F3) are the vector space F3^2. The lines are just the one-dimensional subspaces of this vector space, and their translates. Here's a picture:

Don't worry about the fact that the lines aren't straight. All we're concerned about is which points are incident with which lines. You might just want to convince yourself that every one-dimensional subspace of F3^2 is shown, as well as every translate.

Hopefully this picture makes the analogy between finite geometries and graphs pretty obvious. Points are analogous to vertices; lines are analogous to edges. The main difference is that unlike edges, lines can have more than two points on them. Indeed, if we consider just the points and lines (and not planes, hyperplanes, etc), then a finite geometry gives rise to a "hypergraph" - a graph, but where the "hyperedges" are allowed to contain more than two points. However, finite geometries have some additional structure, since they can also contain planes, hyperplanes etc.

Anyway, we want to calculate the symmetries of finite geometries, and in order to do that, we're going to need to work out what the lines are. First, let's think about how to represent the lines in code. There are several ways we could do it, including at least the following:
- We could represent a line as a list/set of the points on it. We will occasionally want to list all the points on a line, but as a way to represent a line, it's not very efficient.
- A line in the plane has an equation ax+by=c. We could represent a line by the triple (a,b,c). However, that would only work in two dimensions.
- We could represent a line by any pair of distinct points on the line. This will work in any number of dimensions, and will generalise to planes (represented by three points), and so on.

This last approach is the one we will take. So we would like two functions:
- given a pair of points in AG(n,Fq), return all points on the line that they generate
- given n, Fq, return all lines in AG(n,Fq) (represented as a pair of points)

The first is straightforward:


lineAG [p1,p2] = L.sort [ p1 <+> (c *> (p2 <-> p1)) | c <- fq ] where    fq = eltsFq undefined

There are a couple of unfamiliar things here. <+>, <->, *> are the HaskellForMaths functions for adding two vectors, subtracting, and multiplying by a scalar. The second line is just a little bit of phantom type trickery, using type inference to magic up the elements of Fq. Note that the points are returned sorted.

Now to find all the lines. The naive way to do this is as follows:
- list all pairs of points in AG(n,Fq)
- but any pair of points on a line generates the same line, so each line will be listed many times
- so keep only those pairs of points which are the first two points on the line

Here's the code:


linesAG1 n fq = [ [x,y] | [x,y] <- combinationsOf 2 (ptsAG n fq),
                          [x,y] == take 2 (lineAG [x,y]) ]

(This is not a very efficient way to find the lines, but we don't yet have the background to do it more efficiently.)

For example:


> mapM_ print $ linesAG1 2 f3
[[0,0],[0,1]]
[[0,0],[1,0]]
[[0,0],[1,1]]
[[0,0],[1,2]]
[[0,1],[1,0]]
[[0,1],[1,1]]
[[0,1],[1,2]]
[[0,2],[1,0]]
[[0,2],[1,1]]
[[0,2],[1,2]]
[[1,0],[1,1]]
[[2,0],[2,1]]

Or we can list all the points on all the lines:


> mapM_ (print . lineAG) $ linesAG1 2 f3
[[0,0],[0,1],[0,2]]
[[0,0],[1,0],[2,0]]
[[0,0],[1,1],[2,2]]
[[0,0],[1,2],[2,1]]
[[0,1],[1,0],[2,2]]
[[0,1],[1,1],[2,1]]
[[0,1],[1,2],[2,0]]
[[0,2],[1,0],[2,1]]
[[0,2],[1,1],[2,0]]
[[0,2],[1,2],[2,2]]
[[1,0],[1,1],[1,2]]
[[2,0],[2,1],[2,2]]

You can verify for yourself that these are the lines shown in the picture above.

That'll do for now. Next time we'll look at how to work out the symmetries of AG(n,Fq).

Finite fields, part 2

2009-09-02T20:00:00.000+01:00

Okay, so last time, we saw how to extend Q, the field of rational numbers, by adjoining a new element which is the zero of a polynomial in Q[X]. For example, we can adjoin sqrt 2, a zero of X^2-2, to obtain the field Q(sqrt2).

The time before that, we saw that for each prime p, there is a finite field Fp, consisting of the set {0,1,...,p-1}, with the arithmetic operations carried out modulo p.

This week we're going to bring these two ideas together, and look at algebraic extensions of Fp.

Algebraic extensions of Fp work just the same way as algebraic extensions of Q. First, find an irreducible polynomial over the base field K - that means, a polynomial f(x) in K[x], which can't be expressed as f(x) = g(x) * h(x), for g(x), h(x) in K[x] of lower degree. For example, whenever d is an element of K without a square root in K, then x^2-d is irreducible. Then, form the quotient ring K[x]/ - that means, work in K[x], but use the division with remainder algorithm to set f=0, which has the effect of making x act like a zero of f. For f(x)=x^2-d, this means x ends up acting like sqrt d.

For example, in F3, and also in F5, there is no square root of 2. So we can form the fields F3(sqrt2) or F5(sqrt2). We need to be a bit careful about terminology. 2 in F3, 2 in F5, and 2 in Q, are all different things - consequently sqrt2 in F3, in F5, and in Q are all different things. We mustn't make the mistake of thinking they're the same thing.

We saw last time that Q(sqrt2) consists of { a + b sqrt2 | a, b in Q }. Similarly, F5(sqrt2) consists of { a + b sqrt2 | a, b in F5 }. In particular, it follows that F5(sqrt2) has 25 elements. In general, a simple algebraic extension of Fp will have p^n elements, where n is the degree of the irreducible polynomial.

Over Q, we were able to form the fields Q(sqrt2), Q(sqrt3), Q(i), and so on. Somewhat surprisingly, over Fp, it turns out that any algebraic extension of degree n is isomorphic to any other. For example, F5(sqrt2) is isomorphic to F5(sqrt3). For this reason, we typically just talk about F25 (meaning, the field with 25 elements) - or in general, Fq, where q = p^n is a prime power.

However, in order to do arithmetic in F25 or in Fq, we clearly need to have some particular irreducible polynomial in mind. Otherwise, we won't know whether to reduce x^2 to 2, or to 3, for example. So for each p, n, we need to agree on an irreducible polynomial of degree n over Fp. There doesn't appear to be any natural way to choose, so an artificial way to choose has been devised, called the Conway polynomials.

So, using a little phantom type trickery as before, here are the first few finite fields of prime power order, as defined in the HaskellForMaths library:


data ConwayF4
instance PolynomialAsType F2 ConwayF4 where pvalue _ = convert $ x^2+x+1
type F4 = ExtensionField F2 ConwayF4
f4 = map Ext (polys 2 f2) :: [F4]
a4 = embed x :: F4

data ConwayF8
instance PolynomialAsType F2 ConwayF8 where pvalue _ = convert $ x^3+x+1
type F8 = ExtensionField F2 ConwayF8
f8 = map Ext (polys 3 f2) :: [F8]
a8 = embed x :: F8

data ConwayF9
instance PolynomialAsType F3 ConwayF9 where pvalue _ = convert $ x^2+2*x+2
type F9 = ExtensionField F3 ConwayF9
f9 = map Ext (polys 2 f3) :: [F9]
a9 = embed x :: F9

data ConwayF16
instance PolynomialAsType F2 ConwayF16 where pvalue _ = convert $ x^4+x+1
type F16 = ExtensionField F2 ConwayF16
f16 = map Ext (polys 4 f2) :: [F16]
a16 = embed x :: F16

data ConwayF25
instance PolynomialAsType F5 ConwayF25 where pvalue _ = convert $ x^2+4*x+2
type F25 = ExtensionField F5 ConwayF25
f25 = map Ext (polys 2 f5) :: [F25]
a25 = embed x :: F25

As we did with the fields Fp, we provide the functions f4, f8, f9, etc, which just return a list of the elements of the corresponding fields. a4, a8, a9 etc return the element that has been adjoined to the underlying field to make the extension field.

For example:


> f8
[0,a^2,a,a+a^2,1,1+a^2,1+a,1+a+a^2] :: [F8]
> a8 ^ 3
1+a :: F8

We can partially verify the claim that F25 = F5(sqrt2) = F5(sqrt3), by exhibiting square roots of 2 and 3 in F25:


> (2+a25)^2
2 :: F25
> (1+3*a25)^2
3 :: F25

Last time, I mentioned that extension fields can have automorphisms, without really spelling out what that means. Well, if L is an extension of K, then f is an automorphism of the extension if for a,b in L, c in K, f(a+b) = f(a)+f(b), f(a*b) = f(a)*f(b), etc, and f(c*a)= c*f(a). In other words, f(L) looks the same as L, viewed from K.

The finite fields Fq (q = p^n) have an automorphism called the Frobenius automorphism, defined by f(x)=x^p. For example:


> f9
[0,a,2a,1,1+a,1+2a,2,2+a,2+2a]
> map frobeniusAut f9
[0,1+2a,2+a,1,2+2a,a,2,2a,1+a]

Finally, a word about efficiency. The HaskellForMaths implementation of fields Fq is very natural from a mathematical point of view, but it isn't very efficient.

Firstly, we are working with polynomials with coefficients in Fp. To multiply two degree n polynomials together, we have to do O(n^2) additions and multiplications of coefficients. If the coefficients are in Fp, each one of these involves a (`mod` p) operation. It would be better to do the polynomial arithmetic over the integers, and only do the `mod` p at the end. We would then only have to do O(n) (`mod` p) operations.

Secondly, we could consider implementing the fields Fq using Zech logarithms instead.

Anyway, it doesn't matter too much for the moment, as I only want to use the finite fields to construct combinatorial structures. But sometime I might try to write a more efficient implementation.

Next time: finite geometries.

Extension fields

2009-08-27T20:00:00.000+01:00

[New version HaskellForMaths 0.1.7 uploaded here.]

Last week, we looked at the finite fields of prime order, Fp. A field, remember, is a set in which you can do addition, subtraction, multiplication and division. For a given prime p, the field Fp consists of the set {0, 1, ..., p-1}, with arithmetic operations done modulo p. Why does p have to be a prime? We can do addition, subtraction, and multiplication modulo n for any n, prime or not. However, for division, we require p to be prime. (Why?)

You might think, therefore, that the fields Fp of prime order are the only finite fields. However, you would be wrong. We shall see that there are in fact fields Fq of order q (ie, with q elements) for every prime power q = p^n. However, before we can do that we need to understand about algebraic extensions of a field.

It can happen that a field is contained within another field. For example, we have Q

At step 0, start with Q
At step 1, add a
At step n, add all elements that can be obtained by combining two elements from step n-1 using one of the arithmetic operators.

So for example, at step 2, we will have 1+a, 2*a, a*a, 1/a, etc. Anything which appears in the list at step n for some n is in Q(a).

The field Q(a), obtained by adjoining a single element, is called a simple extension of Q. (We could then go on to adjoin another element b, to get Q(a,b), and so on.) Among simple extensions, there are two ways it can go.

Suppose that a is a zero of some polynomial with coefficients in Q - for example x^2-2. In that case, we will construct the whole of Q(a) after only a finite number of steps. For example, in the case where a = sqrt 2, every element of Q(sqrt 2) can be expressed as c + d sqrt 2, for some c, d in Q. For example:
(sqrt 2)^2 = 2
1/(sqrt 2) = (sqrt 2) / 2
1/(1 + sqrt 2) = (1 - sqrt 2) / (1 - 2) = -1 + sqrt 2
This is called an algebraic extension.

Alternatively, if a is not the zero of a polynomial over Q (for example: e, pi), then the construction of Q(a) in steps will never finish. This is called a transcendental extension.

Algebraic extensions of Q are fascinating things. Initially, the main motivation for studying them came from number theory. Some other time, I'd like to take a closer look at them. However, for the moment, I just want to show how to construct them in Haskell, as what we're really aiming for is algebraic extensions of the finite fields Fp.

Let's suppose that we're trying to construct Q(sqrt2). The polynomial we're interested in is x^2-2. Then the basic idea is:

Form the polynomial ring Q[x], consisting of polynomials in x with coefficients in Q
Represent elements of Q(a) as polynomials in Q[x]
Do addition, subtraction, multiplication in Q(a) using the underlying operations in Q[x]
After any arithmetic operation, replace the result by its remainder on division by x^2-2.

Step 4 is the key step. Suppose that we end up with a polynomial f. We use division with remainder to write f = q(x^2-2)+r, and we replace f by r. In effect, this means that we're setting x^2-2 = 0, which is the same as setting x = sqrt 2. If we do this consistently, then the x ends up acting as sqrt2, and we end up with Q(sqrt2).

To construct Q(sqrt3), Q(i), etc, just replace x^2-2 by the appropriate polynomial.

Okay, rather belatedly, time for some code. First, we need to define a type for (univariate) polynomials:


newtype UPoly a = UP [a] deriving (Eq,Ord)

UP [c0,c1,...cn] is to be interpreted as the polynomial c0 + c1 x + ... + cn x^n.

Exercise: Define a Num instance for UPoly a.

If we then define


x = UP [0,1] :: UPoly Integer

together with a suitable Show instance, then we can do things like the following:


> (1+x)^3
1+3x+3x^2+x^3

Exercise: Write quotRemUP, to perform division with remainder in UPoly k, on the assumption that k is a field (that is, a Fractional instance).

So we now have a type, UPoly Q, representing Q[x]. Next, we want to wrap this in another type to represent extension fields Q(a). Rather than have to do this over again each time for Q(sqrt2), Q(sqrt3), Q(i), and so on, we use a little bit of phantom type trickery again. Last time, we used phantom types to represent integers; this time, we're going to use phantom types to represent the polynomials that define the fields (x^2-2, x^2-3, x^2+1, etc).


class PolynomialAsType k poly where
    pvalue :: (k,poly) -> UPoly k

data ExtensionField k poly = Ext (UPoly k) deriving (Eq,Ord)

Here, k represents the field we are extending - Q, to begin with - and if a is the element that we want to adjoin, then poly represents the polynomial over Q of which a is a zero.

From here, it's a short step to define some extension fields:


data Sqrt a = Sqrt a

-- n should be square-free
instance IntegerAsType n => PolynomialAsType Q (Sqrt n) where
    pvalue _ = convert $ x^2 - fromInteger (value (undefined :: n))

type QSqrt2 = ExtensionField Q (Sqrt T2)
sqrt2 = embed x :: QSqrt2

type QSqrt3 = ExtensionField Q (Sqrt T3)
sqrt3 = embed x :: QSqrt3

And now we can do arithmetic in Q(sqrt2), Q(sqrt3), etc. For example:


> :set +t
> (1+sqrt2)^2
3+2a :: QSqrt2
> (1+sqrt3)^2
4+2a :: QSqrt3

As you see, the show function of ExtensionField k poly shows the adjoined element as "a", regardless of which field we're working in. With a little more type hackery, we could have made that a parameter to the type too, so that it would show "sqrt2", "sqrt3", "i", depending which field we were in. However, that would have obscured the code even more. In practice, when we use this code we're only going to be working over one field at a time, so it's fine as it is.

Another limitation of this code is that it's going to be a bit unwieldy to construct Q(a,b) as the type ExtensionField (ExtensionField k poly1) poly2. Luckily, this doesn't matter, as any algebraic extension Q(a,b) is equal to Q(c) for some c.

Exercise: Show that Q(sqrt2, sqrt3) = Q(sqrt2 + sqrt3), and find the polynomial of which sqrt2 + sqrt3 is a zero.

I apologise that I've only really explained the bare bones of how this works. Hopefully you can fill in the gaps. (They're all in the HaskellForMaths source - see link at top of page.) This post has already taken me far too long to write, but there is just one other thing I ought to mention.

Field extensions can have automorphisms (symmetries). Recall that a symmetry is a change that leaves something looking the same. In the case of Q(sqrt2), there is a non-trivial symmetry that sends c + d sqrt 2 to c - d sqrt2 (conjugation). Field automorphisms have a very important place in the history of mathematics: Galois invented group theory in order to study these symmetries, during his investigations into the unsolvability (in general) of the quintic.

Finite fields, part 1

2009-08-17T21:24:00.000+01:00

Okay, so you could think of this as the beginning of chapter two of the imaginary book that I'm writing in this blog. Chapter one was about graphs, groups and symmetry. In chapter two I want to talk about some other highly symmetric combinatorial objects: finite geometries and (probably) designs. But first, I need to talk about finite fields.

A field is roughly a set where you can add, subtract, multiply and divide. You're probably familiar with the following fields:

Q, the rational numbers
R, the real numbers
C, the complex numbers

And also with the following non-fields:

N, the natural numbers. Not a field, because it doesn't contain additive inverses (the negative numbers), so you can't always subtract.
Z, the integers. Not a field, because it doesn't contain multiplicative inverses, so you can't always divide.

I can't resist mentioning a little curiosity. To a number theorist / discrete mathematician like myself, the most interesting of the above are Q and Z. I like to think that this is why they were given the most interesting letters. Indeed, the interestingness of the mathematical objects appears to be correlated with the Scrabble value of the letters: in Scrabble, Q and Z are worth 10, indicating highly interesting; R and N are worth 1, indicating a bit boring; C is worth 3 - moderately interesting. Going further, H, the quaternions, is worth 4 - slightly more interesting than C. However, O, the octonions, is only worth 1, which looks like an anomaly.

Anyway, back to fields. Before we go any further, here's the HaskellForMaths version of Q, the rationals:

import Data.Ratio

newtype Q = Q Rational deriving (Eq,Ord,Num,Fractional)

instance Show Q where
   show (Q x) | b == 1    = show a
              | otherwise = show a ++ "/" ++ show b
              where a = numerator x
                    b = denominator x

Silly, I know, but I just got a bit annoyed with those percentage signs all over the place when using Haskell's built-in rationals.

Anyway, in addition to the fields Q, R, and C, there are many other types of field. The ones I want to look at here are the finite fields of prime order. These are the fields you get if you do arithmetic modulo p, where p is a prime. For each different prime p, we have a field Fp, consisting of the set {0,1, ... p-1}, with arithmetic performed modulo p. For example, in the field F5, we have 3*3 = 4; while in the field F7, we have 3*3 = 2.

Okay, so how to represent these fields in Haskell?

Well, in an earlier version of the HaskellForMaths code, I tried this:


data Fp = F Integer Integer

The idea is that the first integer is p, the prime we are working over, and the second is the value. For example, F 5 3 would represent the value 3 in the field F5.

Next we either derive or define instances of Eq, Ord, Show, Num, Fractional. The Num instance looks something like this:


instance Num Fp where
    F p x + F q y | p == q  = F p $ (x+y) `mod` p
    F p x * F q y | p == q  = F p $ (x*y) `mod` p
    fromInteger _ = error "Fp.fromInteger: not well defined"

The drawbacks of this approach are probably pretty obvious:

Whenever I have two Fp values, F p x and F q y, I first have to check that p = q. This is a run-time check of something that would be better handled at compile time by the type system
When defining numerical types in Haskell, the fromInteger function gives us implicit type conversion from integers to our type, so that we can just write code like 3 * 3, and have the type system convert them into the correct type. However, in this case, we can't give a sensible definition of fromInteger, because we don't know which prime p the user wants to work over. This means that everything is going to have to be written longhand as F 5 3 * F 5 3, and so on.
It feels like a waste of memory to be carrying that first Integer p around all the time, when if I'm using the code correctly, all the ps in a given calculation will be the same.

For these reasons, I decided to try another approach. The aim is to make F2, F3, F5, etc be different types, so that adding 3 :: F5 to 3 :: F7 is a type error, which can be spotted at compile time, rather than a runtime error.

It's fairly simple really.

First, define a whole bunch of phantom types representing the primes:


data T2
data T3
data T5
...

Now define a mechanism for retrieving the values from the types:

class IntegerAsType a where
   value :: a -> Integer

instance IntegerAsType T2 where value _ = 2
instance IntegerAsType T3 where value _ = 3
instance IntegerAsType T5 where value _ = 5
...

This enables us to do the following:

> value (undefined :: T2)
2
> value (undefined :: T3)
3

Now, we use the phantom types to parameterise a finite field data type:


newtype Fp n = Fp Integer

Here, Fp (the one on the left of the equals sign) is a type constructor which takes a type parameter n. This type parameter will be one of the phantom types we just defined. For example:


type F2 = Fp T2
type F3 = Fp T3
type F5 = Fp T5

As before, we need to derive or define Eq, Ord, Show, Num, and Fractional instances. The Num instance looks like this:


instance IntegerAsType n => Num (Fp n) where
    Fp x + Fp y = Fp $ (x+y) `mod` p where p = value (undefined :: n)
    Fp x * Fp y = Fp $ (x*y) `mod` p where p = value (undefined :: n)
    fromInteger m = Fp $ m `mod` p where p = value (undefined :: n)
    ...

This is all fairly standard phantom type trickery, but it now enables us to do things like the following:


> :load Math.Algebra.Field.Base
> 3*3 :: F5
4
> 3*3 :: F7
2
> (3 :: F5) * (3 :: F7)
#type error#

We still have to tell Haskell which prime p we're working over, but now we do it via a type annotation, rather than an argument to the constructor. This means that Haskell can catch errors where we mix values from different fields at compile time, rather than run time.

The module Math.Algebra.Field.Base provides types for the finite fields F2, F3, F5, ..., F97. The module also provides functions to list the elements of each field:


f2 = map fromInteger [0..1] :: [F2]
f3 = map fromInteger [0..2] :: [F3]
f5 = map fromInteger [0..4] :: [F5]
...

And that's about it.

However, we must recognise, in all honesty, that there are also drawbacks to this implementation of the fields Fp.

First, whenever we want to work over a new prime p, we have to define a new type. For example, if we want to work over F101, we are going to have to define a new type for it.

(For the avoidance of doubt: No, type level Peano arithmetic is not what I'm after.)

Second, I can't decide at run time what field to work over. For certain applications it would be very useful to choose my p at run time. I could do that with the first implementation, but I can't do it with the second.

What I really want then, is for Fp and Fq to be different types when p /= q, but to be able to decide at run-time which p to work over. For example, I might have calculated the p I want to work over as the result of a previous calculation. This is impossible in Haskell.

However, there is a technology called "dependent types", in which what I want is possible. Dependent types means roughly that types can depend upon values. So in a Haskell-like language that supported dependent types, I might be able to write code like the following (not valid Haskell):


> 3 * 3 :: F 5
4

Here, the value 5 (not the type T5) is being passed as a type parameter to the type constructor F.

For a while, I thought this was the holy grail. I even considered defecting from Haskell. I did look into languages that support dependent types, such as Coq. However, I concluded that, sadly, this is not the answer. The problem is that - I think - if you have types dependent on values, then you lose the distinction between compile time and run time, which means - I think - that you can't compile the language any more.

Or does someone else know better?

Where we've been, and where we're going

2009-08-05T21:42:00.000+01:00

It feels like about time to take a step back and survey where we've been and where we're going in this blog.

What I'm trying to do in this blog is a kind of guided tour of my HaskellForMaths library. In turn, the library itself, although it contains efficient implementations of fundamental algorithms in computer algebra, is really an educational tool. I wrote it to help me understand maths I wanted to learn about, and especially, to help develop my intuitions about the objects of study. My hope is that by retracing my steps in this blog, I can help others to do the same.

Up to now, we've been looking at graphs and their symmetries, considered as permutation groups. Group theory is sometimes described as the study of symmetry, and I hope that by introducing groups through symmetries of graphs, I've made it clear why.

We've really only looked at quite a small amount of code. The main functions we've seen (forgetting a few stepping stones along the way) are:

graph (vs,es) - for constructing a graph with given vertices and edges
graphAuts g - returns a strong generating set for the group of symmetries / automorphisms of the graph
orderSGS sgs - given a strong generating set, return the order of the group - for example, the total number of symmetries of the graph
conjClassReps - given a generating set for a group, return representatives and sizes for the conjugacy classes - the classes of elements which are "the same, just viewed from a different angle"
p [[1,2,3],[4,5]] - given a list of cycles, return the corresponding permutation
g*h, 1, g^-1 - multiplication, identity and inverses for group elements
v .^ g - the action of a permutation on a vertex
e -^ g - the action of a permutation on an edge
sgs gs - return a strong generating set, given a set of generators

This is the core of the HaskellForMaths code for graphs, graph automorphisms, and permutation groups. Incidentally, I think this code really shows off Haskell to advantage: code to do the same things would be much more complicated in C++, say. I'm especially pleased with the ability to define graphs and groups over arbitrary types, which enables very natural constructions of various classes of graphs and groups. (In this respect, though not in others, I think the library probably beats even dedicated packages such as GAP and MAGMA.)

However, the graph and permutation group code represents, I would say, only about one sixth of what's in the HaskellForMaths library. (Not to mention that I have quite a lot more code in various stages of development back in the lab.) So what else is there still to look at?

Well, what I want to do next is stick with permutation groups a little longer, but widen the set of combinatorial objects that we consider, to include finite geometries, designs, and perhaps other incidence structures.

One of my aims is to take a look at a few of the sporadic finite simple groups. Finite simple groups are the "atoms of symmetry". Most of them fall in one of several infinite families (the cyclic groups of prime order, the alternating groups, and the finite groups of Lie type, aka Chevalley groups). However, it turns out that there are also 26 "sporadic" finite simple groups, which don't belong to any of these families. The HaskellForMaths library contains code for investigating both the infinite families, and some of the smaller sporadic groups.

Then there are root systems and Coxeter groups, including string rewriting and the Knuth-Bendix algorithm. Root systems have some connection to Lie algebras, and are therefore important in physics, though I don't claim to fully understand all that stuff.

Another major part of the library, which I'll get round to eventually, is commutative algebra and Groebner bases. Groebner bases provide a way of calculating with ideals in polynomial rings, which is equivalent to doing algebraic geometry.

Finally, there's some code for working in non-commutative algebra. As an example, there's some code which uses Temperley-Lieb and Hecke algebras to calculate polynomial invariants of knots (the Jones polynomial and the HOMFLY polynomial). In due course, I'm hoping to expand this part of the library to do group algebras, representation theory, and other stuff.

In case I've managed to interest anyone, I should probably mention a few of the books that I used for inspiration:

Ronan, Symmetry and the Monster
Godsil and Royle, Algebraic Graph Theory
Cameron, Combinatorics
Seress, Permutation Group Algorithms
Holt, Handbook of Computational Group Theory

Next time - back to the maths.

How to count the number of positions of Rubik's cube

2009-08-01T19:00:00.000+01:00

Previously on this blog, we have been looking at the symmetries of graphs, which we defined as permutations of the vertices which leave the edges (collectively) in the same places. We think of permutations as active rearrangements of the vertices, rather than as static sequences. For example, we think of the static permutation [2,1,4,3] as the active rearrangement [[1,2],[3,4]], which swaps the 1 and 2 positions, and swaps the 3 and 4 positions. Given two permutations g and h, the active viewpoint enables us to ask what happens when you do g then h, which we think of as multiplication g*h; and how to undo g, which we think of as an inverse g^-1.

The symmetries of a graph form a group. This means, that if g and h are symmetries, then so are g*h, and g^-1. So a group is a set of permutations which is closed under multiplication and inverses.

We can use permutation groups to study all sorts of things besides graph symmetries. As an example, this week I'm going to look at Rubik's cube.

First, we need to work out a way to represent cube moves as permutations. That's easy - we just number the little squares as follows:

Then a clockwise rotation of the front face (the blue face, with an F in the centre), can be represented as the following permutation:

f = [[1,3,9,7],[2,6,8,4],[17,41,33,29],[18,44,32,26],[19,47,31,23]]

This says, that f sends the 1 square to the 3 position, 3 to 9, 9 to 7, 7 to 1 (the corners of the blue face); 2 to 6, 6 to 8, 8 to 4, 4 to 2 (the edges of the blue face); together with three other cycles (the corners and edges adjacent to the blue face).

It's then an easy matter to write down permutations corresponding to the clockwise rotations of the other faces:


b = p [[51,53,59,57],[52,56,58,54],[11,27,39,43],[12,24,38,46],[13,21,37,49]]
l = p [[21,23,29,27],[22,26,28,24],[ 1,31,59,11],[ 4,34,56,14],[ 7,37,53,17]]
r = p [[41,43,49,47],[42,46,48,44],[ 3,13,57,33],[ 6,16,54,36],[ 9,19,51,39]]
u = p [[11,13,19,17],[12,16,18,14],[ 1,21,51,41],[ 2,22,52,42],[ 3,23,53,43]]
d = p [[31,33,39,37],[32,36,38,34],[ 7,47,57,27],[ 8,48,58,28],[ 9,49,59,29]]

Note that I didn't bother to number the centre square in each face, or write down any permutations that move a middle slice. For example, we could have written down:


x = [[4,44,54,24],[5,45,55,25],[6,46,56,26]]

However, doing x is the same as doing u^-1 * d^-1 and then turning the whole cube round. What we're interested in is the possible positions of the 54 squares relative to one another, so we don't really want to count rotations of the whole cube. For this reason, we leave out rotations of the middle slices.

Okay, so we have written down six permutations, corresponding to clockwise rotations of the six faces. Hopefully it's clear that these six elements generate the Rubik's cube group. So we are now in a position to perform calculations within the Rubik's cube group.

Within a group, the order of an element g is defined as the least n>0 such that g^n = 1. In other words, the order of g is the number of times you have to do g to get back where you started. (Recall that the order of a group is the number of elements it contains. The order of an element is equal to the order of the group generated by the element, since this is {g, g^2, ... , g^n-1, g^n = 1}, having n elements.)

For permutations written in cycle notation, it is easy to see the order by inspection. For example:

The order of g = [[1,2]] is 2 - you have to do g twice to get back where you started.
The order of h = [[1,2,3]] is 3 - you have to do h three times to get back where you started.
The order of k = [[1,2,3],[4,5]] is 6. At k^2, we have [[1,3,2]] - the 4 and 5 are in the correct positions, but the 1, 2, 3 are not. At k^3 we have [[4,5]] - now the 1, 2, 3 are in the correct positions, but the 4 and 5 are not. Only at k^6 does everything get back to the starting position.

Hopefully it's clear that the order of an element will always be the least common multiple of the lengths of the cycles:


orderElt g = foldl lcm 1 $ map length $ toCycles g

For example, in the Rubik's cube group we have:


> f*l
[[1,3,9,37,53,17,41,33,27,21,23,19,47,59,11],[2,6,8,34,56,14,4],[7,31,29],[18,44,32,28,24,22,26]]
> orderElt it
105

That is, if you do a front turn followed by a left turn, you will have to repeat 105 times before you get back to the starting position.

When we were looking at graph symmetries, we found that it was particularly useful to have our generators in the form of a strong generating set, since it was then straightforward, for example, to calculate the order of the group.

Rubik's cube isn't a graph. Luckily, the HaskellForMaths library has a function to calculate a strong generating set from any set of generators:


> mapM_ print $ sgs [f,b,l,r,u,d]
[[1,3,9,7],[2,6,8,4],[17,41,33,29],[18,44,32,26],[19,47,31,23]]
[[1,21,51,41],[2,22,52,42],[3,23,53,43],[11,13,19,17],[12,16,18,14]]
[[1,31,59,11],[4,34,56,14],[7,37,53,17],[21,23,29,27],[22,26,28,24]]
[[2,34,22,26,32,44,16,12,46,38,24],[3,31,47,13,49,37,21],[4,8,6,42,52,54,58,56,18,28,14],[7,9,43,39,27,11,41],[19,29,33,51,57,59,53]]
[[3,13,57,33],[6,16,54,36],[9,19,51,39],[41,43,49,47],[42,46,48,44]]
[[4,38,32,46,16,44,34,24,36,14,12],[6,28,56,48,22,52,26,58,8,54,42],[7,31,29],[9,59,57,51,53,33,37,49,13,21,47,27,39,43,11]]
[[6,52,56,54],[9,53,57,51],[11,39,43,47],[12,24,46,44],[13,33,21,49]]
[[7,47,57,27],[8,48,58,28],[9,49,59,29],[31,33,39,37],[32,36,38,34]]
[[8,56,16,52,54],[9,33,47],[12,46,32,24,42],[13,51,43]]
[[9,59],[12,38],[13,49],[24,46],[27,47],[33,37],[39,43],[51,57],[52,58],[54,56]]
[[11,27,39,43],[12,24,38,46],[13,21,37,49],[51,53,59,57],[52,56,58,54]]
[[11,37,39,43],[13,21,59,49],[16,56,54,58],[24,46,38,42],[27,57,51,53]]
[[12,38,46,42,24],[16,56,52,58,54],[27,37,59],[39,57,49]]
[[12,56,58,36],[16,54],[24,38,48,52],[27,37,59],[39,57,49],[42,46]]
[[13,39,43,57,51,49],[16,24,46,38,42,56,54,58],[27,37,59],[36,48]]
[[13,43,51],[27,37,59]]
[[14,28,38,48,54,16,56],[22,34,58,36,46,42,24],[27,59,37],[39,49,57]]
[[16,56,58],[24,38,42],[27,37,59],[39,57,49]]
[[24,36,56,48],[27,37,59],[38,54,58,46],[39,57,49]]
[[24,46,38],[27,37,59],[39,57,49],[54,58,56]]
[[27,39],[36,48],[37,49],[38,46,58,54],[57,59]]
[[27,59,37],[39,49,57]]
[[28,54,58],[34,46,38]]
[[36,38,54],[46,48,58]]
[[36,58,54],[38,46,48]]
[[38,58],[46,54]]

Okay, so strong generating sets aren't all that interesting in themselves - although the last few entries in this one are a little bit interesting, as they show that there are moves which just rotate a pair of corners, or a pair of edges.

Remember though, that SGS allow us to easily calculate the order of a group. In this case:


> orderSGS $ sgs [f,b,l,r,u,d]
43252003274489856000

That's about 4 * 10^19 different positions.

The sgs function, to calculate a strong generating set, uses the Schreier-Sims algorithm. I'll probably explain how it works some time, but not just yet.

Strong generating sets for graph symmetries

2009-07-25T20:00:00.000+01:00

[New release: HaskellForMaths 0.1.6 here]

Previously in this blog, we've been using two functions for finding generators for the group of symmetries of a graph. Both graphAuts2 and graphAuts3 use depth first search to find a transversal generating sequence - the only difference is that graphAuts3 does more pruning of the search tree, and so is faster.

A transversal generating sequence, remember, means this: If our vertices are labelled 1 to n, then we first try to find a symmetry that takes 1 to 2, then another that takes 1 to 3, and so on up to n; then looking only at those that leave 1 fixed, another which takes 2 to 3, then another which takes 2 to 4, and so on; then looking only at those that fix 1 and 2, another which takes 3 to 4, and so on; and so on. So the answer we get always takes the form of a series of levels - first, a set of symmetries taking 1 to some of [2..n], then a set of symmetries fix 1, and take 2 to some of [3..n], and so on.

Okay, so now suppose that we're doing our depth first search. Suppose that we have already found [[1,2],[3,4]], and [[1,3,2]]. Then in fact we don't need to search for a symmetry taking 1 to 4, because the two symmetries we know already generate one. Specifically:


> p [[1,3,2]] * p [[1,2],[3,4]]
[[1,4,3]]

So the point is that, as it happens, it was already possible to take 1 to 4 by repeated application of the symmetries we had already found to take 1 to 2 and 3.

Given a set of group elements, the orbit of a vertex is defined as those vertices that we can send it to by repeated application of the elements:


orbitV gs x = closure [x] [ (.^ g) | g <- gs ]

(This is using the closure algorithm that we defined previously. Recall that v .^ g means the image of v after action by g.)

This idea enables us to write an even faster graphAuts function. As before, we start by looking for symmetries that send 1 to 2, 1 to 3, and so on. However, this time, we don't bother to look for symmetries sending 1 to a vertex which is already in the orbit of the symmetries we have found.

An example will probably make this clearer. First of all, here is the full transversal generating set of symmetries of the cube:


> mapM_ print $ graphAuts3 $ q 3
[[0,1],[2,3],[4,5],[6,7]]
[[0,2,3,1],[4,6,7,5]]
[[0,3],[4,7]]
[[0,4,6,7,3,1],[2,5]]
[[0,5,3],[2,4,7]]
[[0,6,5,3],[1,2,4,7]]
[[0,7],[1,3],[2,5],[4,6]]
[[1,2],[5,6]]
[[1,4,2],[3,5,6]]
[[2,4],[3,5]]

The list consists of three levels: symmetries moving 0 to 1, 2, 3, 4, 5, 6, 7; symmetries fixing 0 and moving 1 to 2, 4; and a symmetry fixing 0 and 1 and moving 2 to 4. The sequence [0,1,2] is called the base for the TGS.

Now, here's what happens if we skip the search for vertices that are already in the orbit, using a new graphAuts function:


> mapM_ print $ graphAuts $ q 3
[[0,1],[2,3],[4,5],[6,7]]
[[0,2,3,1],[4,6,7,5]]
[[0,4,6,7,3,1],[2,5]]
[[1,2],[5,6]]
[[1,4,2],[3,5,6]]
[[2,4],[3,5]]

We still have three levels, with the same base, [0,1,2]. But in the first level, we haven't needed to find as many symmetries. After we found the 0 to 1 and 0 to 2 symmetries, we didn't need to find a 0 to 3 symmetry, because 3 was already in the orbit of the symmetries we had found. Then, after we had found the 0 to 4 symmetry, we didn't need to find the 0 to 5, 6, 7 symmetries.


> orbitV [ p [[0,1],[2,3],[4,5],[6,7]], p [[0,2,3,1],[4,6,7,5]], p [[0,4,6,7,3,1],[2,5]] ] 0
[0,1,2,3,4,5,6,7]

This new graphAuts function (I won't give the code, as it's really only a minor variation on the graphAuts3 function) is faster still, because we're now pruning the search tree even more.

However, it has lost the nice property of the graphAuts2 and graphAuts3 functions, that the returned list was a transversal generating sequence. Recall that a TGS made it particularly easy to work out the order of the group, or list its elements.

However all is not lost. We can easily reconstruct a TGS from the output of the graphAuts function. Consider that first level in the cube symmetries again. In the TGS, we had symmetries taking 0 to 1, 2, 3, 4, 5, 6, 7. The graphAuts function only returned symmetries taking 0 to 1, 2, 4. That was because it turned out that 3, 5, 6, 7 were already in the orbit of 0 under these symmetries. To reconstruct the TGS, what we need to do is, calculate the orbit of 0 under the 0 to 1, 2, 4 symmetries, but this time, keep track of the group elements as we go. For example, the reason 3 is in the orbit is that:


> p [[0,2,3,1],[4,6,7,5]] * p [[0,1],[2,3],[4,5],[6,7]]
[[0,3],[4,7]]

Here's the code. We find the base "bs" by looking at the minimum supports (the least vertex that is moved) of the inputs. We then sort the inputs into levels, using this base. Finally, for each level, we use a modified version of the closure algorithm, that tracks not only where we've got to in the orbit, but also how we got there.


tgsFromSgs sgs = concatMap transversal bs where
    bs = toListSet $ map minsupp sgs
    transversal b = closure b $ filter ( (b <=) . minsupp ) sgs
    closure b gs = closure' M.empty (M.fromList [(b, 1)]) where
        closure' interior boundary
            | M.null boundary = filter (/=1) $ M.elems interior
            | otherwise =
                 let interior' = M.union interior boundary
                     boundary' = M.fromList [(x .^ g, h*g) | (x,h) <- M.toList boundary, g <- gs] M.\\ interior'
                 in closure' interior' boundary'

A set of generators from which we can reconstruct a TGS in this way is called a strong generating set, or SGS. (Strictly speaking, an SGS is relative to a base - we've been using the base that is implied by the Ord instance and the minimum supports of the elements.)

In practice, we prefer to work with SGS than TGS, because they're shorter. Since you can so easily reconstruct a TGS from them, they're just as useful.

I should admit that I'm sidestepping a few subtleties here. The take-home message is:

A strong generating set is a set of generators for a group of a particularly useful form
The graphAuts function gives us a strong generating set by construction

That means that for symmetries of graphs, we are doing pretty well. We have an efficient algorithm (graphAuts) for finding a strong generating set for the symmetry group.

However, what we would like is, given just any set of generators for a group, to be able to construct a strong generating set. For that, we will need the Schreier-Sims algorithm. Next time, I'll show why this would be so useful, by looking at Rubik's cube.

Faster graph symmetries using distance partitions

2009-07-20T21:25:00.000+01:00

Up to now, we've been using the graphAuts2 function to find us a generating set (in fact a transversal generating sequence) for the symmetries of a graph. Unfortunately, graphAuts2 simply isn't fast enough for some of the larger graphs we would like to investigate. (I'll give an example at the end.) This week therefore, we're going to look at a more efficient way to find graph automorphisms, based on distance partitions.

Given a graph, the distance between vertices x and y is defined as the length of the shortest path between them - that is, the number of edges on the path. For example, in the cube shown below, the distance between the 0 and 7 vertices is 3. There are several routes, but each route involves passing along at least 3 edges. In HaskellForMaths, we can use the "distance" function to find this out:


> :load Math.Combinatorics.GraphAuts
> distance (q 3) 0 7
3

In the picture, we see that the vertices fall into four levels, depending on whether their distance from 0 is 0, 1, 2, or 3. Thus, distance from a given vertex can be used to partition the vertices in a graph. This is called the distance partition:


> distancePartition (q 3) 0
[[0],[1,2,4],[3,5,6],[7]]

(The implementation of distancePartition is a variant on the closure algorithm that we have seen a couple of times before, and is left as an exercise.)

Okay, so how do distance partitions help us with finding graph symmetries?

Well, recall that the way graphAuts2 works is, using depth-first search, try to find a symmetry that sends 0 to 1, another that sends 0 to 2, another that sends 0 to 3, and so on, then looking only at those that fix 0, another that sends 1 to 2, another that sends 1 to 3, and so on, then looking only at those that fix 0 and 1, another that sends 2 to 3, another that sends 2 to 4, and so on, and so on.

Well, the first step is to realise that graph symmetries must preserve distance. That is, if g is a symmetry, then we must have: distance x y = distance (x .^ g) (y .^ g). So, suppose that we are looking for a graph symmetry that sends 0 to 1, and we're wondering where to send x. Well, we need only consider those y such that distance 0 x = distance 1 y.

In terms of distance partitions, this means that, once we have decided to map 0 to 1, then we must map each cell in the distance partition of 0 to the corresponding cell in the distance partition of 1.


> distancePartition (q 3) 0
[[0],[1,2,4],[3,5,6],[7]]
> distancePartition (q 3) 1
[[1],[0,3,5],[2,4,7],[6]]

Specifically, we must map {1,2,4} to {0,3,5} (though not necessarily in that order), and {3,5,6} to {2,4,7}, and 7 to 6.

This is already going to cut down our search space considerably, perhaps especially on graphs with many vertices. However, we can go further.

Suppose now that we are looking for symmetries that send 0 to 0 and 1 to 1. Because they send 0 to 0, they must preserve the cells in the distance partition from 0. So they must send {1,2,4} to {1,2,4} (but not necessarily in that order), {3,5,6} to {3,5,6}, and 7 to 7. Also, because they send 1 to 1, they must preserve the cells in the distance partition from 1. So they must send {0,3,5} to {0,3,5}, {2,4,7} to {2,4,7}, and 6 to 6.

I think of this a being a bit like triangulation. Suppose we're wondering where to send 2 to. Well, looking at the distance partition from 0, we see that it must go to one of {1,2,4}. On the other hand, looking at the distance partition from 1, we see that it must go to one of {2,4,7}. So that actually narrows it down to {2,4}. (Yes, I know, in this case we could already see that, but you get the idea.) If we take another fix from a third point, that will narrow it down even further, and so on.

So the idea of our new graphAuts3 function will be as follows.

We will still start by trying to send 0 to 1, 2, 3, 4, 5, 6, 7, and finally 0.
If we're trying to send 0 to x, then when we come to wonder where to send 1, we'll consider only those y with distance x y = distance 0 1. In terms of the distance partition, this means that if 1 falls into cell d of the distance partition from 0, then y must fall into cell d of the distance partition from x.
As we successively decide where to send 1, 2, etc, we will "refine" the cells of the partition by triangulation with the new point.

The code for refining partitions is very simple:


refine p1 p2 = concat [ [c1 `intersect` c2 | c2 <- p2] | c1 <- p1]

For example:


> distancePartition (q 3) 0 `refine` distancePartition (q 3) 1
[[],[0],[],[],[1],[],[2,4],[],[],[3,5],[],[6],[],[],[7],[]]

We get quite a few empty lists in the refinement - they're a necessary evil, but we will remove them as we go along.

Okay, so here's our graphAuts3 function:


graphAuts3 g@(G vs es) = graphAuts' [] [vs] where
   graphAuts' us ((x:ys):pt) =
       let px = refine (ys : pt) (dps M.! x)
           p y = refine ((x : L.delete y ys) : pt) (dps M.! y)
           uus = zip us us
           p' = L.sort $ filter (not . null) $ px
       in concat [take 1 $ dfs ((x,y):uus) px (p y) | y <- ys]
        ++ graphAuts' (x:us) p'
   graphAuts' us ([]:pt) = graphAuts' us pt
   graphAuts' _ [] = []
   dfs xys p1 p2
       | map length p1 /= map length p2 = []
       | otherwise =
            let p1' = filter (not . null) p1
                p2' = filter (not . null) p2
            in if all isSingleton p1'
               then let xys' = xys ++ zip (concat p1') (concat p2')
                    in if isCompatible xys' then [fromPairs' xys'] else []
               else let (x:xs):p1'' = p1'
                        ys:p2'' = p2'
                    in concat [dfs ((x,y):xys)
                                   (refine (xs : p1'') (dps M.! x))
                                   (refine ((L.delete y ys):p2'') (dps M.! y))
                              | y <- ys]
   isCompatible xys = and [([x,x'] `S.member` es') == (L.sort [y,y'] `S.member` es') | (x,y) <- xys, (x',y') <- xys, x < x']
   dps = M.fromList [(v, distancePartition g v) | v <- vs]
   es' = S.fromList es

It looks complicated (and perhaps I could tidy it up a little), but what it's doing is really fairly straightforward. Basically, we're still doing depth first search, in levels, as before. However, we're now maintaining two partitions as we go, the source partition p1 and the target partition p2, and we're constrained to map cells in the source partition to the corresponding cells in the target partition. If we ever find that the "shape" (map length) of the source and target partitions are different, then we know we have gone wrong and need to backtrack. Having satisfied ourselves that the shapes are the same, we can remove those pesky empty lists. Finally, as soon as we find that every cell is a singleton, then we can can shortcut any further search - although we still need to check that the implied mapping is a valid symmetry, to avoid false positives.

Okay, so how about a brief demonstration of its power. Time for a confession - on most of the smallish graphs that we've considered so far, graphAuts3 is actually slightly slower than graphAuts2. However, it's not too difficult to find larger graphs where it wins out.

The Kneser graphs are defined as follows:


kneser n k | 2*k <= n = graph (vs,es) where
   vs = combinationsOf k [1..n]
   es = [ [v1,v2] | [v1,v2] <- combinationsOf 2 vs, disjoint v1 v2]

So the Kneser graph has as vertices the k-subsets of [1..n], with edges joining subsets which are disjoint.

We've already met one of them - the Petersen graph is kneser 5 2. Kneser 7 3 is shown below - it has 35 vertices and 70 edges:

So here's an experiment you can try at home. Compare the running times of graphAuts2 and graphAuts3 on kneser 7 3. On my laptop, graphAuts3 manages to find a transversal generating set of 49 symmetries in less than a second. On the other hand, graphAuts2 had only managed to find 1 symmetry after 10 minutes, at which point I gave up.

Incidentally:


> orderTGS $ graphAuts3 $ kneser 7 3
5040

... which is 7 factorial. That's because the action of S 7 on [1..7] induces an action on kneser 7 3, and it turns out that every symmetry of kneser 7 3 arises from an underlying permutation of [1..7].

Anyway, with graphAuts3, and the orderTGS function from last time (for calculating the number of symmetries, given a transversal generating sequence), we are beginning to have the tools to investigate some larger graphs. However, there is another improvement we can make, which will lead us to strong generating sets and the Schreier-Sims algorithm - next time.

Counting symmetries using transversals

2009-07-15T08:01:00.000+01:00

Previously, we've been using the HaskellForMaths library's graphAuts2 function to find a generating set for the symmetries of a graph. The generating set returned by graphAuts2 has a particularly useful form, which I want to explore this week. Let's just remind ourselves how it works.

> :load Math.Combinatorics.GraphAuts
> mapM_ print $ graphAuts2 $ q 3
[[0,1],[2,3],[4,5],[6,7]]
[[0,2,3,1],[4,6,7,5]]
[[0,3],[4,7]]
[[0,4,6,7,3,1],[2,5]]
[[0,5,3],[2,4,7]]
[[0,6,5,3],[1,2,4,7]]
[[0,7],[1,3],[2,5],[4,6]]
[[1,2],[5,6]]
[[1,4,2],[3,5,6]]
[[2,4],[3,5]]

What graphAuts2 does is, it tries to find a symmetry which sends 0 to 1, another which sends 0 to 2, another which sends 0 to 3, and so on, then looking only at symmetries which leave 0 where it is, another which sends 1 to 2, another which sends 1 to 3, and so on, then leaving 0 and 1 where they are, another which sends 2 to 3, another which sends 2 to 4, and so on, and so on. Most of the time, it won't find all the symmetries it's looking for - as in this case.

Now, a few weeks back I set a puzzle: Given just the list returned by graphAuts2, how can I instantly tell how many symmetries there are in total?

Well, think of it like this.

First of all, I'm trying to find all the different places that I can move the 0 vertex to. Hopefully it's obvious that I can move the 0 vertex to any of the eight vertices. The first seven elements returned by graphAuts2 are symmetries which move 0 to 1,2,3,4,5,6,7 respectively - plus of course I can always just leave the cube as it is and leave 0 at 0.

Now, suppose that I have decided to leave 0 where it is. Next, I try to find all the different places that I can move 1 to, but leaving 0 where it is. If you look at the picture, I hope you can see that having fixed 0, our only choices for 1 are 1, 2, and 4. The next two elements returned by graphAuts2 are symmetries which move 1 to 2 and 4 respectively - plus of course, the identity element leaves the cube as it is, so leaves 0 at 0 and 1 at 1.

Now, suppose that I have decided to leave 0 and 1 where they are. Next I try to find all the different places that I can move 2 to. If you look at the picture, you'll see that if we fix 0 and 1, our only choices for 2 are 2 and 4. The identity leaves 2 where it is, and the last element returned by graphAuts2 moves 2 to 4.

Now, suppose that we decide to leave 0, 1 and 2 where they are. Next I try to find all the different places I can move 3 to. graphAuts2 didn't return any more elements. What it is telling us is that once I have decided to fix 0, 1, and 2, I have no choice but to fix 3 and all the other vertices too.

Okay, so how many symmetries are there in total? Well, I had 8 choices for the first vertex, then 3 choices for the second vertex, then 2 choices for the third vertex, then no more choices. So that gives us 8*3*2 = 48. Let's just check:


> length $ elts $ graphAuts2 $ q 3
48

We can write code to do this for us.


orderTGS tgs =
    let transversals = map (1:) $ L.groupBy (\g h -> minsupp g == minsupp h) tgs
    in product $ map L.genericLength transversals

This code needs a bit of explaining:

The order of a group means the number of elements. That's what we're trying to calculate.
TGS stands for transversal generating set. This is my name for a generating set of the type returned by graphAuts2. (I'll explain the name in a little bit.)
Given a TGS, the first thing we need to do is group the elements into levels based on what I call the minimum support - the least vertex that they move. In the case of the cube, we have three levels - seven elements that move 0, two elements that fix 0 and move 1, and one element that fixes 0 and 1 but moves 2
Next we add the identity element (called 1) to each level, to show that we can also just leave the vertex where it is
Finally, to calculate the order, we multiply the number of elements in each level.


> orderTGS $ graphAuts2 $ q 3
48

Let's think a little more about how this works. What are transversals, and what's so special about "transversal generating sets"? Well, consider the following picture:

On the left, the 48 dots represent the 48 symmetries of the cube. The eight columns divide the elements according to whether they send 0 to 0, 1, 2, 3, 4, 5, 6, or 7. The red dot in the top left is the identity, which sends 0 to 0. Then graphAuts2 finds us one dot (the blue one) in each column. So we have a representative from each column. This set of representatives is called a transversal.

Next, on the right, we confine our attention to only the first column, the elements that send 0 to 0. The three blocks divide the elements by whether they send 1 to 1, 2, or 4. The red dot in the top block is again the identity. Then graphAuts2 finds us a blue dot in each of the other two blocks. So we have another transversal.

Finally (not shown), we divide the top left block by whether the elements send 2 to 2 or 4. In this case we have just two divisions, consisting of a single element each - the identity, and the last element found by graphAuts2. Taken together, these two elements constitute our third transversal.

So a transversal generating set - it should really be called a transversal generating sequence, because the order matters - is a sequence of transversals (but omitting the identity), starting at the outermost layer and going successively inwards, that generates the group.

And now it should be totally clear why this method of calculating the order of the group works. The only thing which perhaps isn't clear is, how do we know that there are the same number of elements in each column, and in each block, at each stage? (This is clearly required, if the method is to work.) Well, the answer is simply that multiplication by the blue dot in any column gives a one-to-one correspondence between the elements in the left hand column and the elements in that column, so the number of elements must be the same.

As this last remark perhaps suggests, we can also use these transversals to list the elements of a group:


eltsTGS tgs =
    let transversals = map (1:) $ L.groupBy (\g h -> minsupp g == minsupp h) tgs
    in map product $ sequence transversals

So once again, we construct some transversals by grouping the TGS into levels, and then adding the identity to each level. Then what the last line is saying is, to get an element of the group, just pick an element from each of the transversals (levels), and multiply them together. To get all the elements of the group, consider all possible such choices.

The reason that transversal generating sets and the orderTGS function are important, is that they are the first step towards being able to work with large or very large groups. Previously, if we wanted to know how many elements there are in a group, we would have had to use the "elts" function to generate them all, and then count them. This is okay for small groups, but for graphs with thousands or millions of symmetries, it is not going to be practical. With the orderTGS function, we will be able to find out how many symmetries a graph has, even if it is in the thousands or millions, so long as we can find a TGS.

Our next stumbling block is that the graphAuts2 function itself is not efficient enough to handle these larger graphs. So next time, we'll look at a more efficient way to find a TGS for a graph. After that, we'll look at something called a strong generating set, and an algorithm called the Schreier-Sims algorithm, that can always find us a strong generating set, given only some generators for the group.

Conjugacy classes, part 2

2009-07-08T08:03:00.000+01:00

Last time we saw that the symmetries of a graph can be divided into conjugacy classes, where the elements of each class are somehow the same symmetry, just viewed from a different angle. This time we're going to look at the conjugacy classes of symmetries of a few simple graphs, to get more of a feel for how this works.

Let's start with the graph of the cube, otherwise known (in graph theory) as q 3. q n has as vertices all points in {0,1}ⁿ, with edges between pairs of vertices which differ in only one position:


q' k = graph (vs,es) where
    vs = sequence $ replicate k [0,1]
    es = [ [u,v] | [u,v] <- combinationsOf 2 vs, hammingDistance u v == 1 ]
    hammingDistance as bs = length $ filter id $ zipWith (/=) as bs

For example:


> :load Math.Combinatorics.GraphAuts
> q' 3
G [[0,0,0],[0,0,1],[0,1,0],[0,1,1],[1,0,0],[1,0,1],[1,1,0],[1,1,1]]
[[[0,0,0],[0,0,1]],[[0,0,0],[0,1,0]],[[0,0,0],[1,0,0]],[[0,0,1],[0,1,1]],
[[0,0,1],[1,0,1]],[[0,1,0],[0,1,1]],[[0,1,0],[1,1,0]],[[0,1,1],[1,1,1]],
[[1,0,0],[1,0,1]],[[1,0,0],[1,1,0]],[[1,0,1],[1,1,1]],[[1,1,0],[1,1,1]]]

However, this is a bit hard to read, a jumble of 0s, 1s and square brackets, so we normally work over integers instead, by considering the sequences of 0s and 1s as the bits in the binary representation of a number:


> q 3
G [0,1,2,3,4,5,6,7]
[[0,1],[0,2],[0,4],[1,3],[1,5],[2,3],[2,6],[3,7],[4,5],[4,6],[5,7],[6,7]]

Okay, so we would like to investigate the symmetries of the cube, and as we saw last time, what we're really interested in is just the different classes of symmetry, rather than listing every single symmetry.


> mapM_ print $ conjClassReps $ graphAuts2 $ q 3
([],1)
([[0,1],[2,3],[4,5],[6,7]],3)
([[0,1],[2,5],[3,4],[6,7]],6)
([[0,1,3,2],[4,5,7,6]],6)
([[0,1,3,7,6,4],[2,5]],8)
([[0,3],[1,2],[4,7],[5,6]],3)
([[0,3,6,5],[1,2,7,4]],6)
([[0,3,6],[1,7,4]],8)
([[0,3],[4,7]],6)
([[0,7],[1,6],[2,5],[3,4]],1)

What we need to do is look at the representatives of each class, and try to understand what they are doing, so as to be able to give a textual description:

[] is the identity permutation, which leaves everything where it is
[[0,1],[2,3],[4,5],[6,7]] is a reflection in a plane midway between opposing faces. There are three such elements, because there are three pairs of opposing faces.
[[0,1],[2,5],[3,4],[6,7]] is a 180 degree rotation about an axis joining the midpoints of opposite edges. There are six elements, because there are six pairs of opposite edges.
[[0,1,3,2],[4,5,7,6]] is a 90 degree rotation about an axis through the centres of opposing faces. There are six faces (or three pairs of faces, with clockwise and anti-clockwise turns).
[[0,1,3,7,6,4],[2,5]] - this is an interesting one. There is a 6-cycle, so there's some sort of 60 degree rotation going on. But there is also a 2-cycle, so some sort of reflection or 180 degree rotation. If you think about it, you'll see that what is happening is that you're holding the cube by opposite corners (2 and 5 here), and spinning it by 60 degrees - but at the same time, reflecting it in the plane midway between.
Exercise: Finish off the list

Previously we saw that the automorphism group of the complete graph k n is the symmetric group S n.

The conjugacy classes of S n are particularly easy to understand. Let's look at an example:


> mapM_ print $ conjClassReps $ _S 5
([],1)
([[1,2]],10)
([[1,2],[3,4]],15)
([[1,2],[3,4,5]],20)
([[1,2,3]],20)
([[1,2,3,4]],30)
([[1,2,3,4,5]],24)

k 5 doesn't have an obvious spatial interpretation (unlike the pentagon and the cube that we have considered previously), so we're going to need a more abstract description of the classes.

Let's have a look at one of the classes:

> conjClass (_S 5) (p [[1,2]])
[[[1,2]],[[1,3]],[[1,4]],[[1,5]],[[2,3]],[[2,4]],[[2,5]],[[3,4]],[[3,5]],[[4,5]]]

There might be a slight surprise here. Surely [[1,2]] and [[1,3]] don't belong in the same class, you might think, because 1 and 2 are next to each other, and 1 and 3 aren't. Or to take another example, we have [[1,2,3,4,5]], and [[1,3,5,2,4]] in the same class - surely that can't be right - we saw when looking at c 5 that one is a 1/5 rotation, and the other is a 2/5 rotation.

Well, the first thing to point out is that we can be misled by the spatial representation of a graph. In the picture of k 5, it might look like the relationship between 1 and 2 is not the same as the relationship between 1 and 3. However, from the point of the view of the graph, it is. An alien who saw in graph space rather than in 2-d or 3-d space wouldn't be able to tell a difference.

Another thing to emphasize is that conjugacy classes are relative to the group you're working in. For example, the elements [[1,2,3,4,5]] and [[1,3,5,2,4]] occur in both S 5 and D 10 (the symmetry group of c 5). In S 5, they are in the same conjugacy class, but in D 10, they are not. Specifically, we have:


> p [[1,2,3,4,5]] ~^ p [[2,3,5,4]]
[[1,3,5,2,4]]

[[2,3,5,4]] is an element of S 5, so [[1,2,3,4,5]] and [[1,3,5,2,4]] are conjugate in S 5. However, [[2,3,5,4]] is not an element of D 10, so they are not conjugate in D 10.

We can translate this back to our intuition about conjugate elements "looking the same but from a different angle". We have our alien who sees in graph space. When our alien is sitting in the graph space k 5, then [[2,3,5,4]], being a symmetry of the graph, is a transformation that moves the graph to a different angle. On the other hand, in the graph space c 5, [[2,3,5,4]] is not a symmetry - instead of moving the graph to a different angle, it just crumples the graph up.

Let's look a little more closely at that conjugation operation again:


> p [[1,2,3,4,5]] ~^ p [[2,3,5,4]]
[[1,3,5,2,4]]

Let's compare where we started - [[1,2,3,4,5]] - with where we ended - [[1,3,5,2,4]] - what has changed? Well, the 2 turned into a 3, the 3 turned into a 5, the 5 turned into a 4, and the 4 turned into a 2. So conjugating by [[2,3,5,4]] has the same effect as actually applying [[2,3,5,4]] to each of the numbers in the cycle notation for [[1,2,3,4,5]].

Why is this? Well, remember that:


g ~^ h = h^-1 * g * h

So if h is [[2,3,5,4]], what this says is:

First, undo [[2,3,5,4]]. That is, put 3 into the 2 position, 5 into the 3 position, and so on.
Next, do g. But if g says to do something to 2, we'll actually be doing it to 3, which is in the 2 position, and so on.
Finally, do [[2,3,5,4]]. So put whatever's in the 2 position back into the 3 position, put what's in the 3 position back in the 5 position, and so on.

For example, if g = [[1,2,3,4,5]], then the overall effect on 3 is: put into the 2 position (h^-1), then put into the 3 position (g), then put into the 5 position (h). Notice the way that 3 gets off the merry-go-round at one point, but gets on again at a different point, because someone rearranged the chairs in the middle.

Well anyway, I'm sure if you think about it hard enough, you'll get it.

Getting back to the conjugacy classes of S n, it turns out that a conjugacy class in S n just consists of all elements having the same "cycle shape". For example, the class of [[1,2]] consists of all elements having a single 2-cycle. The class of [[1,2],[3,4]] consists of all elements having two 2-cycles. The class of [[1,2],[3,4,5]] consists of all elements having a 2-cycle and a 3-cycle.

The reason for this should now be obvious. To get from [[1,2],[3,4,5]] to [[a,b],[c,d,e]], you just need to find the permutation that sends 1 to a, 2 to b, 3 to c, 4 to d, and 5 to e. Since S 5 contains all permutations of [1..5], this permutation must be among them.

Hope that all made sense.