Google Summer of Code 2012: Week 10

Posted by: Александър Макелов | July 30, 2012

Google Summer of Code 2012: Week 10

Hi all,

here’s a brief summary of what I’ve been doing during the 10th week of my GSoC.

Though I fixed a bug in the SUBGROUPSEARCH function during the week, I ran some more comprehensive tests as I had planned to, and some of them broke the function. If you’re particularly interested, something like that will work:

In [87]: S = SymmetricGroup(5)
In [88]: prop_fix_3 = lambda x: x(3) == 3
In [89]: %autoreload
In [90]: S.subgroup_search(prop_fix_3)
---------------------------------------------------------------------------
StopIteration                             Traceback (most recent call last)
<ipython-input-90-6b85aa1285b8> in <module>()
----> 1 S.subgroup_search(prop_fix_3)

/home/alexander/workspace/sympy/sympy/combinatorics/perm_groups.py in subgroup_search(self, prop, base, strong_gens, tests, init_subgroup)
2660
2661                 # this function maintains a partial BSGS structure up to position l
-> 2662                 _insert_point_in_base(res, res_base, res_strong_gens, l, new_point, distr_gens=res_distr_gens, basic_orbits=res_basic_orbits, transversals=res_transversals)
2663                 # find the l+1-th basic stabilizer
2664                 new_stab = PermutationGroup(res_distr_gens[l + 1])

/home/alexander/workspace/sympy/sympy/combinatorics/util.py in _insert_point_in_base(group, base, strong_gens, pos, point, distr_gens, basic_orbits, transversals)
423     # baseswap with the partial BSGS structures. Notice that we need only
424     # the orbit and transversal of the new point under the last stabilizer
--> 425     new_base, new_strong_gens = group.baseswap(partial_base, strong_gens, pos, randomized=False, transversals=partial_transversals, basic_orbits=partial_basic_orbits, distr_gens=partial_distr_gens)
426     # amend the basic orbits and transversals
427     stab_pos = PermutationGroup(distr_gens[pos])

/home/alexander/workspace/sympy/sympy/combinatorics/perm_groups.py in baseswap(self, base, strong_gens, pos, randomized, transversals, basic_orbits, distr_gens)
2472             # ruling out member of the basic orbit of base[pos] along the way
2473             while len(current_group.orbit(base[pos])) != size:
-> 2474                 gamma = iter(Gamma).next()
2475                 x = transversals[pos][gamma]
2476                 x_inverse = ~x

StopIteration:

The reason is certainly the change of base performed on line 11 in the pseudocode (this is also indicated in my code on my local week6 branch here ). The use of the function BASESWAP there is what gets us into trouble. It is meant to be applied to base and a strong generating set relative to it, switch two consecutive base points and change the generating set accordinly. However, in subgroup_search the goal is to change a base $(b_1, b_2, \ldots, b_l, \ldots, b_k)$ to $(b_1, b_2, \ldots, b_l', \ldots, b_k)$ where $b_l'$ is a new point. The book ([1]) mentions that this is done by using BASESWAP but doesn’t provide any details. My strategy is the following: I cut the base so that it becomes $(b_1, b_2,\ldots, b_l)$ and cut the correponding data structures – the strong generators strong_gens, the basic_orbits, the transversals, and the strong generators distributed according to membership in basic stabilizers distr_gens (I know, I still have to rename this to strong_gens_distr). Then I append the point $b_l'$ so that the base is $(b_1, b_2, \ldots, b_l, b_l')$ and calculate an orbit and transversal for $b_l’$ under the stabilzier of $b_1, b_2, \ldots, b_l$ . Finally I apply BASESWAP to this new base in order to switch the two rightmost points. Then I go back to $(b_1, b_2, \ldots, b_l', \ldots, b_k)$ by appending what I had cut in the start and calculating a transversal/orbit for $b_{l+1}$ under the stabilizer just found, that of $b_1, \ldots, b_l'$ . Obviously, the resulting BSGS structures are valid only up to position $l$ , and that’s all the information we can acquire without another application of baseswap or finding another stabilizer ( and in general, finding a stabilizer is a computationally hard task relative to calculating orbits/transversals). The entire purpose of this use of BASESWAP in SUBGROUPSEARCH is to obtain generators for the stabilizer of $b_1, b_2, \ldots, b_l'$ and maintain a base/strong generating set that are valid up to a certain position. There are many such base changes performed on the same base throughout the course of the function and something goes wrong along the way. I still have to figure out why and where.

The good news: There is a straightforward alternative to using BASESWAP: maintain a list of generators for each of the basic stabilizers in $(b_1, b_2, \ldots, b_k)$ and change it accordingly as the base is changed, using the function stabilizer() in sympy/combinatorics/perm_groups.py. For each base change we have to calculate one more stabilizer, so that’s not terrible. It is also sort of suggested in “Notes on Computational Group Theory” by Alexander Hulpke (page 34). The problem with this approach is that stabilizer() tends to return a group with many generators, and repeated applications keep increasing this number. However, using this removed the bug from SUBGROUPSEARCH. As before, more comprehensive tests are on the way : )
Yet another alternative : we can use the incremental Schreier-Sims algorithm with the new base $(b_1, \ldots, b_l', \ldots, b_k)$ and the strong generating set for $(b_1, \ldots, b_l, \ldots, b_k)$ . There will likely be redundant generators after that, and it will probably involve more computation than finding a single stabilizer. However, in the long run (since there are many base changes performed) this might perform faster (due to the increasing number of generators that stabilizer() tends to create). I have not tried that approach yet.
Other than that, I had my latest major pull request merged! Thanks a lot to Stefan and my mentor David for the review! That was the largest one so far…
I started reading about some of the applications of subgroup search; subgroup intersection seems to be the easiest to implement, so I’ll probably go for it first.