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 G0. 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.
