Not really blogging

A. Makelov

Google Summer of Code 2012: Week 9

5 Comments Posted by on July 23, 2012

Hi all, here’s a brief summary of what I’ve been doing for the 9th week of my GSoC.

This week saw (and still has to see) some exciting new additions:

I. The incremental Schreier-Sims algorithm.

This is a version of the Schreier-Sims algorithm that takes a sequence of points B and a generating set S for a group G as input, and extends B to a base and S to a strong generating set relative to it. It is described in [1], pp.87-93. The default value of B is [], and that of S is \text{G.generators}. Here’s an example:


In [41]: S = SymmetricGroup(5)
In [42]: base = [3, 4]
In [43]: gens = S.generators
In [44]: x = S.schreier_sims_incremental(base, gens)
In [45]: x
Out[45]:
([3, 4, 0, 1],
[Permutation([1, 2, 3, 4, 0]),
Permutation([1, 0, 2, 3, 4]),
Permutation([4, 0, 1, 3, 2]),
Permutation([0, 2, 1, 3, 4])])
In [46]: from sympy.combinatorics.util import _verify_bsgs
In [47]: _verify_bsgs(S, x[0], x[1])
Out[47]: True

The current implementation stores the transversals for the basic orbits explicitly (the alternative is to use Schreier vectors to describe the orbits – this saves a lot of space, but requires more time in order to compute transversal elements whenever they are needed. This feature is still to be implemented, and this probably won’t happen in this GSoC). The current implementation of the Schreier-Sims algorithm on the master branch uses Jerrum’s filter (for more details and comparisons of the incremental version and the one using Jerrum’s filter, go here) as an optimization, and also stores the transversals explicitly. The incremental version seems to be asymptotically faster though. Here’s several comparisons of the current version on the master branch and the incremental one which can be found on a local branch of mine which is somewhat inadequately called week6):

For symmetric groups:


In [50]: groups = []
In [51]: for i in range(20, 30):
....:     groups.append(SymmetricGroup(i))
....:
In [52]: for group in groups:
....:     %timeit -r1 -n1 group.schreier_sims()
....:
1 loops, best of 1: 590 ms per loop
1 loops, best of 1: 719 ms per loop
1 loops, best of 1: 981 ms per loop
1 loops, best of 1: 1.35 s per loop
1 loops, best of 1: 1.66 s per loop
1 loops, best of 1: 2.19 s per loop
1 loops, best of 1: 2.74 s per loop
1 loops, best of 1: 3.37 s per loop
1 loops, best of 1: 4.28 s per loop
1 loops, best of 1: 5.37 s per loop
In [53]: for group in groups:
....:     %timeit -r1 -n1 group.schreier_sims_incremental()
....:
1 loops, best of 1: 612 ms per loop
1 loops, best of 1: 737 ms per loop
1 loops, best of 1: 927 ms per loop
1 loops, best of 1: 1.15 s per loop
1 loops, best of 1: 1.41 s per loop
1 loops, best of 1: 1.72 s per loop
1 loops, best of 1: 2.1 s per loop
1 loops, best of 1: 2.52 s per loop
1 loops, best of 1: 3.02 s per loop
1 loops, best of 1: 3.58 s per loop

For alternating groups:


In [54]: groups = []
In [55]: for i in range(20, 40, 2):
....:     groups.append(AlternatingGroup(i))
....:
In [56]: for group in groups:
%timeit -r1 -n1 group.schreier_sims()
....:
1 loops, best of 1: 613 ms per loop
1 loops, best of 1: 1.03 s per loop
1 loops, best of 1: 1.77 s per loop
1 loops, best of 1: 2.65 s per loop
1 loops, best of 1: 3.51 s per loop
1 loops, best of 1: 5.31 s per loop
1 loops, best of 1: 7.71 s per loop
1 loops, best of 1: 11.1 s per loop
1 loops, best of 1: 15.3 s per loop
1 loops, best of 1: 19.1 s per loop
In [57]: for group in groups:
%timeit -r1 -n1 group.schreier_sims_incremental()
....:
1 loops, best of 1: 504 ms per loop
1 loops, best of 1: 787 ms per loop
1 loops, best of 1: 1.23 s per loop
1 loops, best of 1: 1.9 s per loop
1 loops, best of 1: 2.8 s per loop
1 loops, best of 1: 3.99 s per loop
1 loops, best of 1: 5.48 s per loop
1 loops, best of 1: 7.45 s per loop
1 loops, best of 1: 10 s per loop
1 loops, best of 1: 13.2 s per loop

And for some dihedral groups of large degree (to illustrate the case of small-base groups of large degrees):


In [58]: groups = []
In [59]: for i in range(100, 2000, 200):
....:     groups.append(DihedralGroup(i))
....:
In [60]: for group in groups:
%timeit -r1 -n1 group.schreier_sims()
....:
1 loops, best of 1: 29.6 ms per loop
1 loops, best of 1: 108 ms per loop
1 loops, best of 1: 278 ms per loop
1 loops, best of 1: 527 ms per loop
1 loops, best of 1: 861 ms per loop
1 loops, best of 1: 1.29 s per loop
1 loops, best of 1: 1.83 s per loop
1 loops, best of 1: 2.39 s per loop
1 loops, best of 1: 3.06 s per loop
1 loops, best of 1: 3.83 s per loop
In [61]: for group in groups:
%timeit -r1 -n1 group.schreier_sims_incremental()
....:
1 loops, best of 1: 20.8 ms per loop
1 loops, best of 1: 52.8 ms per loop
1 loops, best of 1: 121 ms per loop
1 loops, best of 1: 223 ms per loop
1 loops, best of 1: 365 ms per loop
1 loops, best of 1: 548 ms per loop
1 loops, best of 1: 766 ms per loop
1 loops, best of 1: 1 s per loop
1 loops, best of 1: 1.25 s per loop
1 loops, best of 1: 1.51 s per loop

In addition to this algorithm I implemented a related function _remove_gens in sympy.combinatorics.util which removes redundant generators from a strong generating set (since there tend to be some redundant ones after schreier_sims_incremental() is run):


In [68]: from sympy.combinatorics.util import _remove_gens
In [69]: S = SymmetricGroup(6)
In [70]: base, strong_gens = S.schreier_sims_incremental()
In [71]: strong_gens
Out[71]:
[Permutation([1, 2, 3, 4, 5, 0]),
Permutation([1, 0, 2, 3, 4, 5]),
Permutation([0, 5, 1, 2, 3, 4]),
Permutation([0, 1, 2, 3, 5, 4]),
Permutation([0, 1, 2, 4, 3, 5]),
Permutation([0, 1, 3, 2, 4, 5]),
Permutation([0, 1, 2, 5, 4, 3]),
Permutation([0, 1, 5, 3, 4, 2])]
In [72]: new_gens = _remove_gens(base, strong_gens)
In [73]: new_gens
Out[73]:
[Permutation([1, 0, 2, 3, 4, 5]),
Permutation([0, 5, 1, 2, 3, 4]),
Permutation([0, 1, 2, 4, 3, 5]),
Permutation([0, 1, 2, 5, 4, 3]),
Permutation([0, 1, 5, 3, 4, 2])]
In [74]: _verify_bsgs(S, base, new_gens)
Out[74]: True

II. Subgroup search.
This is an algorithm used to find the subgroup K of a given group G of all elements of G satisfying a given property P. It is described in [1], pp.114-118 and is quite sophisticated (the book is right when it says “The function SUBGROUPSEARCH is rather complicated and will require careful study by the reader.”). On the other hand, it is one of the most interesting additions to the groups module to date since it can do so much. The idea is to do a depth-first search over all group elements and prune large parts of the search tree based on several different criteria. It’s currently about 150 lines of code and works in many cases but still needs debugging. It can currently do some wonderful stuff like this:


In [77]: S = SymmetricGroup(6)
In [78]: prop = lambda g: g.is_even
In [79]: G = S.subgroup_search(prop)
In [80]: G == AlternatingGroup(6)
Out[80]: True

to find the alternating group as a subgroup of the full symmetric group by the defining property that all its elements are the even permutations, or this:


In [81]: D = DihedralGroup(10)
In [82]: prop_true = lambda g: True
In [83]: G = D.subgroup_search(prop_true)
In [84]: G == D
Out[84]: True

to find the dihedral group D_{10} as a subgroup of itself using the trivial property that always returns \text{True}; or this:


In [106]: A = AlternatingGroup(4)
In [107]: G = A.subgroup_search(prop_fix_23)
In [108]: G == A.stabilizer(2).stabilizer(3)
Out[108]: True

to find the pointwise stabilizer of \{2,3\}. And so on and so on. What is more wonderful is that you can specify the base used for G in advance, and the generating set returned for K will be a strong generating set with respect to that base!


In [119]: A = AlternatingGroup(5)
In [120]: base, strong_gens = A.schreier_sims_incremental()
In [121]: G = A.subgroup_search(prop_fix_1, base=base, strong_gens=strong_gens)
In [122]: G == A.stabilizer(1)
Out[122]: True
In [123]: _verify_bsgs(G, base, G.generators)
Out[123]: True

The bad news is that the function breaks somewhere. For example:


In [125]: S = SymmetricGroup(7)
In [126]: prop_true = lambda g: True
In [127]: G = S.subgroup_search(prop_true)
In [128]: G == S
Out[128]: False

This needs some really careful debugging, but overall it looks promising since it works in so many cases – so the bug is hopefully small : ).

So, that’s it for now!

[1] Derek F. Holt, Bettina Eick, Bettina, Eamonn A. O’Brien, “Handbook of computational group theory”, Discrete Mathematics and its Applications (Boca Raton). Chapman & Hall/CRC, Boca Raton, FL, 2005. ISBN 1-58488-372-3

All, Google Summer of Code 2012

5 responses to “Google Summer of Code 2012: Week 9”

  1. Mario July 23, 2012 at 9:48 am

    I find that in the symmetric group case schreier_sims() is faster
    than schreier_sims_incremental();
    I did the benchmarking using your week6 branch;
    in the case use_named_groups=False an efficient SGS is used;
    in the case use_named_groups=True it is used SymmetricGroup;
    notice that the sgs is much longer in this case

    use_named_groups=False
    n schreier_sims schreier_sims_incremental len(sgs)
    20 0.025 0.12 19
    40 0.22 1.4 39
    60 0.85 6.3 59
    80 2.3 19 79
    100 5.0 44 99
    120 9.7 86 119

    use_named_groups=True
    n schreier_sims schreier_sims_incremental len(sgs)
    20 0.32 0.62 108
    40 20 31 569
    60 240 303 1416

    here is the benchmark script:

    from time import time
from sympy.combinatorics.permutations import Permutation, cyclic
from sympy.combinatorics.perm_groups import PermutationGroup
from sympy.combinatorics.named_groups import SymmetricGroup
def get_symmetric_group_sgs(n):
    gens = [Permutation(Permutation(cyclic([(i,i+1)], n)).array_form) for i in range(1, n)]
    return gens


def bench_sims(n, use_incremental=False, use_named_groups=False):
    gens = get_symmetric_group_sgs(n)
    G1 = PermutationGroup(get_symmetric_group_sgs(n))
    t0 = time()
    if use_incremental:
        G1.schreier_sims_incremental()
    else:
        G1.schreier_sims()
    # check that the SGS is the same
    sgs1 = G1.strong_gens
    t1 = time()
    assert gens == sgs1
    print 'n=%d len(sgs1)=%d  %.3f' %(n, len(sgs1), t1 - t0),
    if not use_named_groups:
        print ''
        return
    t0 = time()
    G2 = SymmetricGroup(n)
    if use_incremental:
        G2.schreier_sims_incremental()
    else:
        G2.schreier_sims()
    sgs2 = G2.strong_gens
    t1 = time()
    assert G1 == G2
    print 'len(sgs2)=%d  %.3f' %(len(sgs2), t1 - t0)

bench_sims(20, 1, 1)

    Like

    Reply
  2. amakelov July 23, 2012 at 11:17 am

    Notice that the current implementation of schreier_sims_incremental() returns the base and strong generating set as a tuple and doesn’t set any of the attributes of the group:

    In [2]: SymmetricGroup(5).schreier_sims_incremental()
Out[2]: 
([0, 1, 3, 2],
 [Permutation([1, 2, 3, 4, 0]),
  Permutation([1, 0, 2, 3, 4]),
  Permutation([0, 4, 1, 2, 3]),
  Permutation([0, 1, 2, 4, 3]),
  Permutation([0, 1, 3, 2, 4]),
  Permutation([0, 1, 4, 3, 2])])

    Instead, the current implementation of schreier_sims() (the one on the week6 branch) sets attributes like self.base, self.strong_gens and so on. When we call for one of these and it is not already set, for example strong_gens as in your script, we actually run schreier_sims() on the group because we have this in the sympy.combinatorics.perm_groups.py file:

        @property
    def strong_gens(self):
        """
        ...
        """
        if self._strong_gens == []:
            self.schreier_sims()
        return self._strong_gens

    So in your script, when you call for G1.strong_gens, you actually call schreier_sims() if use_incremental was set to True. Thus the strong generating sets you are referring to are always obtained by schreier_sims(). The ones obtained by schreier_sims_incremental() are of reasonable length, and they can be decreased further by using _remove_gens:

    In [3]: from sympy.combinatorics.util import _remove_gens
In [4]: base, strong_gens = SymmetricGroup(40).schreier_sims_incremental()
In [5]: len(strong_gens)
Out[5]: 76
In [6]: new_gens = _remove_gens(base, strong_gens)
In [7]: len(new_gens)
Out[7]: 39

    The incremental version of Schreier-Sims is definitely faster than schreier_sims() for symmetric groups (using SymmetricGroup) on my machine. I’ll look further into this issue of speed…

    Like

    Reply
  3. mario July 23, 2012 at 1:43 pm

    Thanks for the explanation.

    new_gens = _remove_gens(base, strong_gens)

    Is new_gens a SGS ? _remove_gens lacks the docstring.

    Like

    Reply
    • amakelov July 23, 2012 at 1:58 pm

      Yep, it’s a strong generating set – the idea of _remove_gens is to remove redundant generators from a strong generating set while keeping it a strong generating set. I realize that it must be quite painful to work with my current week6 branch since it is a work in progress, so if there’s something unclear don’t hesitate to ask 🙂 I’ll try to write docstrings for Schreier Sims and some of the new functions in sympy.combinatorics.util and get them up on my branch.

      Like

      Reply
  4. amakelov July 24, 2012 at 4:49 pm

    Edit: I believe I just fixed the bug in subgroup_search. It was a silly one, in this fragment:

            # line 6: initializations
        c = [1]*base_len
        u = [identity]*base_len
        sorted_orbits = [None]*base_len
        for i in range(base_len):
            sorted_orbits[i] = basic_orbits[i][:]
            sorted_orbits[i].sort(key = lambda point: base_ordering[point])

    The line 

            c = [1]*base_len

    should be replaced by

            c = [0]*base_len

    For yet another time, I have forgotten that counting starts from zero, as opposed to the book where it starts from 1 : )
    Anyway, I’ll fix this and write some comprehensive tests for the function in the near future.

    Like

    Reply

Leave a comment