Stabiliser Chains

In the previous chapter we used GAP’s built-in Orbit, Stabilizer, and in on the ball group. GAP’s internal representation of a permutation group is a stabiliser chain, and all three of those operations lean on it. This chapter builds one from scratch.

We will keep using the ball group as our running example:

L := (1,2,5,4);;
R := (2,3,6,5);;
G := Group(L, R);;

Orbit calculation

We start by computing the orbit of a point. The orbit of a point \(i\) under a group \(G\) is the set \(i^G = \{i^g : g \in G\}\) — every value that some element of \(G\) maps \(i\) to. The orbit always contains \(i\) itself, because the identity is in \(G\) and fixes \(i\).

In the ball group, the orbit of 1 is all six positions. Starting from 1, we can reach 2, 5, and 4 by repeatedly applying \(L\) (recall \(L = (1,2,5,4)\)). We reach 3 by first applying \(L\) to get to 2, then \(R\), which sends 2 to 3. And we reach 6 by \(L\) followed by \(R\) twice.

That exploration is essentially the whole algorithm. Every element of a finite group is a product of generators, so rather than enumerate all 120 elements of \(G\) and see where each sends 1, we start from 1, apply every generator to every point we have found, and keep going until no new points appear.

I find this is easier to see in code, so let’s write some GAP. The function below is a little inefficient, but will do as a first attempt.

CalculateOrbitSlow := function(G, point)
    local knownOrbit, p, g, gens;
    gens := GeneratorsOfGroup(G);

1    knownOrbit := [point];
    for p in knownOrbit do
        for g in gens do
            if not(p^g in knownOrbit) then
                Add(knownOrbit, p^g);
            fi;
        od;
    od;
    return knownOrbit;
end;

1: Start knownOrbit with the only point we know is in the orbit — point.

function( G, point ) ... end

As we extend knownOrbit, the for p in knownOrbit loop keeps running over the new points too. This is fine — it’s the standard breadth-first exploration.

CalculateOrbitSlow(G, 1);

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

CalculateOrbitSlow(G, 3);

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

This function has two limitations. First, it is inefficient: the line p^g in knownOrbit searches through the whole orbit we have seen so far. Second, we often want to know which permutation takes the base point to each point in the orbit — but we throw that information away.

Let’s do an improved version, called a Schreier vector. Pick a point \(b\) — the base point — and we will compute the orbit of \(b\), recording, at every other orbit point \(x\), a generator that moves \(x\) one step closer to \(b\). Chain those generators together and we get a permutation in \(G\) which maps \(x\) all the way back to \(b\).

This is easier to follow once the code is in front of you, so here it is:

CalculateSchreier := function(G, point)
    local knownOrbit, vec, img, p, g, gens;
    gens := GeneratorsOfGroup(G);

    vec := [];
    knownOrbit := [point];
    vec[point] := ();

    for p in knownOrbit do
        for g in gens do
            img := p^(g^-1);
            if not(IsBound(vec[img])) then
                vec[img] := g;
                Add(knownOrbit, img);
            fi;
        od;
    od;

    return rec(generators := gens,
               orbit := knownOrbit,
               transversal := vec);
end;

function( G, point ) ... end

A quick note on img := p^(g^-1). The caret ^ is doing two different jobs here. On the right, g^-1 inverts the permutation g. On the left, x^h — for a point x and a permutation h — returns the image of x under h. So p^(g^-1) is the point whose image under g is p: equivalently, the point one step “further from the base” along the generator g.

Why compute p^(g^-1) rather than the more obvious p^g? So that walking back to the base is cheap. On the next line we store g at vec[img]. When we later walk back from img, we apply g directly and land at p, which is one step closer to the base. Had we stored things the other way round (img := p^g with vec[img] := g), the walk back would need g^-1 at every step, and we’d be inverting permutations in the inner loop. Computing p^(g^-1) up front moves that cost to the BFS, where each inversion happens once.

So what is a Schreier vector, formally? The record CalculateSchreier returns holds two pieces: the orbit itself (in the order we discovered the points), and a transversal. A transversal of a subgroup \(H \le G\) is a set of coset representatives — one element of \(G\) per coset of \(H\). The subgroup we care about here is the stabiliser \(G_b\). Each orbit point corresponds to exactly one right coset of \(G_b\): two elements \(g_1, g_2 \in G\) send \(b\) to the same place precisely when \(g_1 g_2^{-1} \in G_b\). Our Schreier vector doesn’t store the representatives directly. Instead, at each orbit point it records one generator that moves the point one step closer to the base, and chaining those generators together produces a representative on demand.

That is exactly what RepresentativePerm does. Given a Schreier vector sch with base \(b\), and a point \(x\) in the orbit, it returns some \(g \in G\) with \(x^g = b\). If \(x\) is not in the orbit, it returns fail.

RepresentativePerm := function(schvec, val)
    local ret, gen;
    if not(IsBound(schvec.transversal[val])) then
        return fail;
    fi;

    ret := ();

    while val <> schvec.orbit[1] do
        gen := schvec.transversal[val];
        val := val^gen;
        ret := ret*gen;
    od;

    return ret;
end;

function( schvec, val ) ... end

Let’s see it in action:

sch := CalculateSchreier(G, 1);;
RepresentativePerm(sch, 2);

(1,4,5,2)

RepresentativePerm(sch, 3);

(1,5,3)(2,6,4)

RepresentativePerm(sch, 6);

(1,5,4,2,3,6)

Each of these should send its input point back to the base (1). Let’s check directly:

2^RepresentativePerm(sch, 2);

3^RepresentativePerm(sch, 3);

6^RepresentativePerm(sch, 6);

Reducing to the stabiliser

Now we have RepresentativePerm, let’s use it to tackle an important group problem: checking whether a given permutation is in \(G\).

The key is that groups are closed. If \(p\) and \(q\) are in the group, so is \(p \cdot q\). This also works backwards: if \(q\) and \(p \cdot q\) are both in the group, then so is \(p\), because \(p = (p \cdot q) \cdot q^{-1}\), and inverses are in the group too.

Let’s see how this helps. Consider the permutation \((1,4)(2,5)(3,6)\), which maps 1 to 4. First ask whether 1 can be mapped to 4 by any element of \(G\) — that is, whether 4 is in the orbit of 1. If it isn’t, well, this permutation certainly isn’t in \(G\), and we are done.

For the ball group, 4 is in the orbit of 1, so we continue. RepresentativePerm gives us a specific element \(r \in G\) that maps 4 back to 1. Now look at \((1,4)(2,5)(3,6) \cdot r\). This product:

is in \(G\) if and only if \((1,4)(2,5)(3,6)\) is (by the closure argument above, since \(r\) is in \(G\)), and
maps 1 to 1, because \((1,4)(2,5)(3,6)\) sends 1 to 4 and \(r\) sends 4 back.

So if \((1,4)(2,5)(3,6) \cdot r\) is in \(G\), it lies in \(G_1\) — the stabiliser of 1 in \(G\), the subgroup of elements that fix 1. And that’s a smaller group (20 elements, not 120), so we have reduced to an easier version of the same question.

This trick is exactly what MapToBase does:

MapToBase := function(schvec, perm)
    local map;
    map := RepresentativePerm(schvec, schvec.orbit[1]^perm);
    if map = fail then return fail; fi;
    return perm * map;
end;

function( schvec, perm ) ... end

It multiplies the input by a permutation in \(G\) that cancels its action on the base point. The result fixes the base — or fail, when the input sends the base outside the orbit (in which case the input cannot possibly be in \(G\)).

MapToBase(sch, (1,4)(2,5)(3,6));

(2,4)(3,6)

MapToBase(sch, (3,6));

(3,6)

Both results should fix 1. Let’s check:

1^MapToBase(sch, (1,4)(2,5)(3,6));

1^MapToBase(sch, (3,6));

The fail branch doesn’t trigger for the ball group because \(G\) is transitive on \(\{1, \ldots, 6\}\) — every permutation of these six points sends 1 somewhere in the orbit. In a non-transitive setting, fail would fire whenever the input sent the base to a point \(G\) cannot reach.

Building the chain

So how do we solve the problem of checking whether a permutation is in \(G_1\)? We do the same thing again — build a Schreier vector for some point not already fixed by \(G_1\), and reduce to the stabiliser of that point. Eventually we will find ourselves in a subgroup that fixes everything: the trivial group, whose only element is the identity. At that point membership is obvious.

This nested sequence of Schreier vectors is a stabiliser chain. Each level of the chain holds:

A base point — for now, we just pick the smallest point the current group moves.
A Schreier vector for that base point.
Optionally, if there is still some group to go, another level below — describing the stabiliser of the base point.

The sequence \(b_1, b_2, \ldots, b_k\) of base points we pick is called a base for \(G\). Formally, a base is a list \([b_1, \ldots, b_k]\) whose pointwise stabiliser \(G_{b_1, \ldots, b_k}\) is trivial.

Every finite permutation group has a base — in the worst case, fixing every point the group acts on always leaves the trivial group. But there is no unique base. For the ball group, \([1, 2, 3]\) is a base (the one our chain is about to build), but so are \([3, 1, 5]\) and \([6, 2, 4]\), among many others. Different algorithms prefer different bases: some want the shortest possible base so the chain has few levels, some want a longer base with small orbits at each step. We are going to pick the smallest point the current group moves at every step, via SmallestMovedPoint — a simple default rather than an optimal choice.

Let’s put that in GAP. We have to call our function MyStabilizerChain, because GAP already has a StabilizerChain built in:

MyStabilizerChain := function(G)
    local root, pnt, Gstab;
    pnt := SmallestMovedPoint(G);
    root := CalculateSchreier(G, pnt);
    Gstab := Stabilizer(G, pnt);
    if not IsTrivial(Gstab) then
        root.stabilizer := MyStabilizerChain(Gstab);
    fi;
    return root;
end;

function( G ) ... end

A quick honest note here. We are leaning on GAP’s built-in Stabilizer to build our own stabiliser chain. That is circular — Stabilizer is itself implemented using stabiliser chains inside GAP.

Building a chain from scratch means finding not just a base but also a strong generating set: a set of generators for \(G\) which also generates each stabiliser in the chain. A base together with a strong generating set is called a base and strong generating set, or BSGS, and is the standard way to present a permutation group for computation. Finding a BSGS efficiently, and verifying that it is correct, is a substantial algorithm of its own (Schreier–Sims and its descendants), which we won’t cover here. We just need to know what a stabiliser chain is and what it lets us compute.

Let’s build one for the ball group and look at it. The base point at each level is the first entry of that level’s orbit:

chain := MyStabilizerChain(G);;
[chain.orbit[1],
 chain.stabilizer.orbit[1],
 chain.stabilizer.stabilizer.orbit[1]];

[ 1, 2, 3 ]

chain.orbit;

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

chain.stabilizer.orbit;

[ 2, 5, 6, 3, 4 ]

chain.stabilizer.stabilizer.orbit;

[ 3, 6, 4, 5 ]

Three levels, with base \([1, 2, 3]\) and orbit sizes 6, 5, 4. The chain bottoms out at the trivial group after that.

Using the chain

Now we are in a position to check whether a permutation lies in \(G\). The function follows a basic recursive design:

Look where the permutation \(p\) maps the base point. Is that outside the orbit of the base point? Then \(p\) is not in \(G\).
Otherwise, calculate a permutation \(q\) in \(G\) such that \(p \cdot q\) maps the base point to itself, then recurse in the stabiliser.

The terminating case is when we reach the bottom of the chain. At that point the permutation we are tracking has been reduced, level by level, to fix each base point in turn — so it fixes every \(b_i\) in our base, and therefore lies in the pointwise stabiliser \(G_{b_1, \ldots, b_k}\). By the defining property of a base, that subgroup is trivial. So our surviving permutation is either \(()\) — in which case the original was in \(G\) — or something non-identity, in which case the original was not in \(G\), because each MapToBase step preserves membership.

PermInGroup := function(chain, perm)
    local basemap;
    basemap := MapToBase(chain, perm);
    if basemap = fail then return false; fi;
    if not IsBound(chain.stabilizer) then
        return basemap = ();
    fi;
    return PermInGroup(chain.stabilizer, basemap);
end;

function( chain, perm ) ... end

Let’s test it against permutations that are and aren’t in \(G\):

PermInGroup(chain, (1,4)(2,5)(3,6));

true

PermInGroup(chain, (3,6));

false

PermInGroup(chain, (1,2,3,4,5,6));

true

The first is the row-swap from the introduction. The second is the “swap 3 and 6” permutation we already noted is unreachable. The third is a 6-cycle — also in \(G\), which you probably wouldn’t have guessed from looking at \(L\) and \(R\) alone. As a product of the generators it is \(R^{-1} L^{-1} R^{-1} L^{-1} R\) (five moves, found by GAP). Membership is genuinely non-obvious — which is why we built the chain.

Group size

What else can we do with a stabiliser chain? Lots of things — almost anything which explores a group. Here is one easy example: finding the size of the group.

The orbit–stabiliser theorem says \(|G| = |b^G| \cdot |G_b|\), because (as we saw above) the orbit points biject with the right cosets of the stabiliser. Applied recursively down the chain, this says \(|G|\) is the product of orbit sizes at each level.

GroupSize := function(stabchain)
    if not IsBound(stabchain.stabilizer) then
        return Length(stabchain.orbit);
    fi;
    return Length(stabchain.orbit) *
           GroupSize(stabchain.stabilizer);
end;

GroupSize(chain);

function( stabchain ) ... end

For the ball group this is \(6 \times 5 \times 4 = 120\), agreeing with Size(G).

The chain also describes \(G\) compactly. Its three Schreier vectors hold \(6 + 5 + 4 = 15\) entries between them, compared to 120 elements if we listed the group. For \(S_{10}\) the gap is dramatic: a chain of \(10 + 9 + \ldots + 2 = 54\) entries versus \(10! = 3{,}628{,}800\) group elements.

Many more interesting calculations build on stabiliser chains — including group intersection, which is where we are heading. Except, efficient intersection needs more ideas than just stabiliser chains. In the next few chapters we will build intuition for one of the key ideas, refinement-based search, by looking at McKay’s algorithm for graph isomorphism. We will then return to intersection.