Iteration through finite Coxeter groups

[stumpc5]

The algorithm in chevie is very different from the algorithm we currently have:

ForEachElement:=function(W,f)local l,g;
  if not IsFinite(W) then Error("only for finite Coxeter groups");
  elif W.nbGeneratingReflections=0 then f(W.identity);return;
  fi;
  l:=List([0..W.nbGeneratingReflections],i->
        ReflectionSubgroup(W,W.reflectionsLabels{[1..i]}));
  l:=List([1..Length(l)-1],i->ReducedRightCosetRepresentatives(l[i+1],l[i]));
  g:=function(x,v)local y;
    if Length(v)=0 then f(x);
    else for y in v[1] do g(x*y,v{[2..Length(v)]});od;
    fi;
  end;
  g(W.identity,l);
end;

We should also implement that algorithm (maybe even with hard-coded coset representatives in the finite case) so that we can see if this is indeed faster.

The two needed methods are

#############################################################################
##
#F  ReducedInRightCoset( <W> , <w> )  . . . . . reduced element in coset W.w
##
##  w is an automorphism of a parent Coxeter group of W.
##  'ReducedInRightCoset' returns the unique element in the right
##  coset W.w which is W-reduced.
##

AbsCoxOps.ReducedInRightCoset:=function(W,w)local i;
  while true  do
    i := FirstLeftDescending(W, w);
    if i=false then return w; fi;
    w := W.reflections[i] * w;
  od;
end;

#############################################################################
##
#F  ReducedRightCosetRepresentatives( <W>, <H> )  . . . . . . . . . . . . .
#F  . . . . . . . . . . . . . . .  distinguished right coset representatives
##
##  'ReducedRightCosetRepresentatives'  returns a list of reduced words in the
##  Coxeter group <W> that are distinguished right  coset representatives for
##  the right cosets H\W where  H is a reflection subgroup  of W.
##
ReducedRightCosetRepresentatives:=function(W, H)local res, totest, new;
  totest:=Set([W.identity]);
  res:=Set([W.identity]);
  repeat
    new:=Concatenation(List(totest,w->List(
      W.reflections{W.generatingReflections},s->ReducedInRightCoset(H, w*s))));
    UniteSet(res,totest);
    totest:=Difference(new,res);
  until Length(totest)=0;
  InfoChevie2("#I nb. of cosets: ",Length(res),"\n");
  SortBy(res,x->CoxeterLength(W,x));
  return res;
end;

[stumpc5]

###########################################################################
#
# N.B. a difference with older versions of Chevie is that IsLeftDescending and
# FirstLeftDescending  return  an  index  in  the  list  of  generators, not a
# reflectionName,  for  efficiency.  LeftDescentSet  still  returns  names  of
# reflections.
#
###########################################################################

##  generic reduction of FirstLeftDescending to IsLeftDescending
##  in practical implementations of Coxeter groups this routine can often
##  be overriden by a more efficient one.

AbsCoxOps.FirstLeftDescending:=function(W, x)local i,ILD;
  ILD:=W.operations.IsLeftDescending; # avoid dispatching overhead
  for i in W.generatingReflections do if ILD(W,x,i) then return i; fi; od;
  return false;
end;

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​a2ef2e4	added a few sentences to the doc
​dc01066	started to add 'optional - gap3'
​fd0c9cc	added optional gap3
​1537a79	added optional gap3
​60acebc	Merge branch 'u/stumpc5/11187' into u/stumpc5/20445
​3d5bd19	first version of the algorithm in chevie

[stumpc5]

Last 10 new commits:

​8810c41	fixed one doctest + work on the invariant form
​2ba3740	fixed another doctest + a typo
​e28bbec	edited the length function to always use the reduced word
​8f79ffa	just playing with first methods
​a2ef2e4	added a few sentences to the doc
​dc01066	started to add 'optional - gap3'
​fd0c9cc	added optional gap3
​1537a79	added optional gap3
​60acebc	Merge branch 'u/stumpc5/11187' into u/stumpc5/20445
​3d5bd19	first version of the algorithm in chevie

[stumpc5]

Travis, I implemented a first version of the algoritm, getting all needed cosets in E7 takes 2/3 of our iteration time, but constructing the elements in the end takes about 30secs. If you find the time, I am sure you see how to speed the algorithm. The function is parabolic_iteration.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​a60cf31	added some missing optional gap3
​aaf59b1	Merge branch 'u/stumpc5/11187' into u/stumpc5/20445
​eb75f71	turned it into an iterator to make every product only once. somhow results in infinite iteration

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​feebf97	finalized the invariant form
​4955d68	added the missing optional gap3 in the cython file
​c5c43bc	the next round of optional gap3 insertions
​6d788c9	Merge branch 'u/stumpc5/11187' into u/stumpc5/20445
​26d89e2	- made the computation of reduced_coset_repesentatives really fast ~0.1sec in E8

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​e7eb93c

pushed a new version of the parabolic algorithm that is 30% faster than ours

[stumpc5]

sage: timeit("for w in W.iteration('depth',False): pass",number=5)
5 loops, best of 3: 6.33 s per loop
sage: from sage.combinat.root_system.reflection_group_c import parabolic_iteration
sage: timeit("for w in parabolic_iteration(W): pass",number=5)    
5 loops, best of 3: 4.28 s per loop

Okay, I can now get down quite a bit from our last weeks algorithm. Drawback is that it needs quite some memory (for E8, we have to keep E7 in memory). I provide an alternative in which order the parabolic is computed, but that doesn't seem to speed the computation.

• If @tscrim looks at the code and improves it further, we might beat the 5 minutes for E8!
• The code can also perfectly be paralellized!

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​7eebd4b

minor tweak

[stumpc5]

I locally reimplemented the multiplication of PermutationGroupElement. This resulted in

sage: W = ReflectionGroup(['E',7])                                                
sage: from sage.combinat.root_system.reflection_group_c import parabolic_iteration
sage: timeit("for w in parabolic_iteration(W): pass",number=5)                    
5 loops, best of 3: 2.4 s per loop

If we provide a very restricted implementation of permutations ourselves, we might get down further quite a bit. Remark: the algorithm without the multiplication and creation of new permutation group elements only takes .5sec, so there is a lot to improve still.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​88ec6a0

implemented a local perm grp elt multiplication

[stumpc5]

I just redid the computations and it seems to be about twice as fast as in gap3:

sage: W = ReflectionGroup(['E',7])
sage: from sage.combinat.root_system.reflection_group_c import parabolic_iteration
sage: timeit("for w in parabolic_iteration(W): pass",number=5)
5 loops, best of 3: 1.56 s per loop
sage: timeit("for w in parabolic_iteration(W): pass",number=5)
5 loops, best of 3: 1.43 s per loop
sage: W = ReflectionGroup(['E',8])
sage: %time for w in parabolic_iteration(W): pass

CPU times: user 9min 12s, sys: 3.81 s, total: 9min 16s
Wall time: 9min 16s

So we are slowly approaching the 5min...

[tscrim]

Took 3.5Gb of RAM with this ticket, but on my laptop, this is what I got:

sage: W = ReflectionGroup(['E',8])
sage: from sage.combinat.root_system.reflection_group_c import parabolic_iteration
sage: %time for w in parabolic_iteration(W): pass
CPU times: user 4min 23s, sys: 604 ms, total: 4min 23s
Wall time: 4min 23s

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​244e6ff	Merge branch 'public/11187' of trac.sagemath.org:sage into public/reflection_groups-11187
​f541349	Some last little beautification.
​bc3c926	Merge branch 'u/stumpc5/11187' into u/stumpc5/20445

[stumpc5]

Replying to tscrim:

CPU times: user 4min 23s, sys: 604 ms, total: 4min 23s Wall time: 4min 23s

Great, I believe you can claim to be the first person ever who iterated through E8 in less than 5 minutes! (I was expecting that since 1. your computer was about twice as fast as mine last week, and 2. you took ~8 minutes to iter through E8 in chevie, and my code is about twice as fast as the chevie code.)

1. I would really like to see how to implement our own permutation group elements with only what we need. I would hope that to again result in some speedup. One question there: we have w(-\beta) = -w(\beta), so we would not need to record the complete permutation on all roots, but on the positive roots would be enough. That would speed several computations such as creating new elements and testing for equality), but it would have the drawback that we constantly need to work mod N (N=nr of positive roots), e.g., we would have such checks and mod's when multiplying two permutations.

2. We should get your parallelization to work with this so that people can then use many cores to actually do stuff with the elements in type E8.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​0db28d6	fixed doctests that have not been run last week!
​7afc02c	Fixing doctests in combinat/root_system/coxeter_group.py.
​0269425	Blanket fix bare except block.
​9b0ef92	Making it an AttributeError.
​0985005	Merge branch 'public/11187' into u/stumpc5/20445
​e740893	uses the local multiplication now also in the other iterator

[stumpc5]

Travis, could you also now run

sage: W = ReflectionGroup(['E',8])
sage: %time for w in W.iteration("depth",False): pass

It now also uses the local multiplication, so it should also speed quite a bit. On my machine, it's only 1.5 times as slow.

[tscrim]

With doing other things on my laptop:

sage: W = ReflectionGroup(['E',8])
sage: %time for w in W.iteration("depth",False): pass
CPU times: user 6min 39s, sys: 22.5 ms, total: 6min 39s
Wall time: 6min 39s

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​3e25053	Merge branch 'public/11187' into u/stumpc5/20445
​abfff5f	Merge branch 'develop' into u/stumpc5/20445

[stumpc5]

@Nicolas: I just now see that you (I think) accidentally edited my comment (comment 5 above from here) instead of replying three weeks ago. Could you quickly recheck?

[nthiery]

Replying to stumpc5:

1. I would really like to see how to implement our own permutation group elements with only what we need. I would hope that to again result in some speedup. One question there: we have w(-\beta) = -w(\beta), so we would not need to record the complete permutation on all roots, but on the positive roots would be enough. That would speed several computations such as creating new elements and testing for equality), but it would have the drawback that we constantly need to work mod N (N=nr of positive roots), e.g., we would have such checks and mod's when multiplying two permutations.

This business sounds of the same nature as what we have for affine permutations (in window notation). Would there be a way to use the same implementation behind the scene?

Cheers,

Nicolas

[nthiery]

Replying to stumpc5:

@Nicolas: I just now see that you (I think) accidentally edited my comment (comment 5 above from here) instead of replying three weeks ago. Could you quickly recheck?

Oops, reverted!

[stumpc5]

Replying to nthiery:

Replying to stumpc5: This business sounds of the same nature as what we have for affine permutations (in window notation). Would there be a way to use the same implementation behind the scene?

And also with colored permutations, isn't it? The main difference is that here (and also in signed permutations) one works mod N, for colored permutations one works "+N mod kN", and for affine permutation one does not work mod anything.

I would propose to first work out the implementation here and then see if we can use it also in the other places. I only don't see how to actually do the implementation in an optimal way, so some support of yours and/or Travis is appreciated.

Some concrete questions:

1. It seems that we should use the same data structure as for PermutationGroupElement:

       self.perm = <int *>sig_malloc(sizeof(int) * self.N)
   

   Do you agree? Can we even get anything significantly better than sticking to PermutationGroupElement if we do it ourselves? This also asks whether we can do better when multiplying elements, I do not see what

    cdef PermutationGroupElement prod = PermutationGroupElement.__new__(PermutationGroupElement)
   

   does or how long it takes, see the method _new_mul_ in reflection_group_c.pyx.

3. It seems that we are using PermutationGroupElement in a few places (when talking to GAP3}), but this might just be that we need the cycle string representation for that.

[stumpc5]

@Travis: Could you have a brief look at the _new_mul_ in reflection_group_c.pyx to check whether I am missing something, or whether I could use that as a first improvement in this algorithm.

I indeed plan to have this ticket only contain the improvements for now (then having plenty of options for iterating through finite Coxeter groups), postponing the work and testing for a new data structure to a new ticket.

[tscrim]

Replying to stumpc5:

I would propose to first work out the implementation here and then see if we can use it also in the other places. I only don't see how to actually do the implementation in an optimal way, so some support of yours and/or Travis is appreciated.

I'm really starting to consider that what we should do is create our own separate project where we write all of this independently (in, say, C/C++). At this stage, I'm somewhat concerned with the additional overhead that Cython could impose and the lack of complete memory control. Although I cannot commit serious time to working on this for then next two weeks (I will be in grading and math mode). Over the summer starting in June, I should be able to do so.

Some concrete questions:

1. It seems that we should use the same data structure as for PermutationGroupElement:

       self.perm = <int *>sig_malloc(sizeof(int) * self.N)
   

   Do you agree? Can we even get anything significantly better than sticking to PermutationGroupElement if we do it ourselves?

I think we can avoid a bit of overhead of maintaining the GAP and category information. Although it is hard to tell how much of an impact this will have on things.

This also asks whether we can do better when multiplying elements, I do not see what

 cdef PermutationGroupElement prod = PermutationGroupElement.__new__(PermutationGroupElement)

does or how long it takes, see the method _new_mul_ in reflection_group_c.pyx.

This creates a new element in memory, but it does not call the __init__. It is essential and done in Python kernel, so we won't get any better.

3. It seems that we are using PermutationGroupElement in a few places (when talking to GAP3}), but this might just be that we need the cycle string representation for that.

GAP4 doesn't store things as cycle strings AFAIK, and so I doubt GAP3 does either.Replying to stumpc5:

I would propose to first work out the implementation here and then see if we can use it also in the other places. I only don't see how to actually do the implementation in an optimal way, so some support of yours and/or Travis is appreciated.

I'm really starting to consider that what we should do is create our own separate project where we write all of this independently (in, say, C/C++). At this stage, I'm somewhat concerned with the additional overhead that Cython could impose and the lack of complete memory control. Although I cannot commit serious time to working on this for then next two weeks (I will be in grading and math mode). Over the summer starting in June, I should be able to do so.

Some concrete questions:

1. It seems that we should use the same data structure as for PermutationGroupElement:

       self.perm = <int *>sig_malloc(sizeof(int) * self.N)
   

   Do you agree? Can we even get anything significantly better than sticking to PermutationGroupElement if we do it ourselves?

I think we can avoid a bit of overhead of maintaining the GAP and category information. Although it is hard to tell how much of an impact this will have on things.

This also asks whether we can do better when multiplying elements, I do not see what

 cdef PermutationGroupElement prod = PermutationGroupElement.__new__(PermutationGroupElement)

does or how long it takes, see the method _new_mul_ in reflection_group_c.pyx.

This creates a new element in memory, but it does not call the __init__. It is essential and done in Python kernel, so we won't get any better.

3. It seems that we are using PermutationGroupElement in a few places (when talking to GAP3}), but this might just be that we need the cycle string representation for that.

GAP4 doesn't store things as cycle strings AFAIK, and so I doubt GAP3 does either. We could very likely replace these calls with something better (if we are calling GAP3).

[stumpc5]

Replying to tscrim:

Replying to stumpc5: I'm really starting to consider that what we should do is create our own separate project where we write all of this independently (in, say, C/C++). At this stage, I'm somewhat concerned with the additional overhead that Cython could impose and the lack of complete memory control. Although I cannot commit serious time to working on this for then next two weeks (I will be in grading and math mode). Over the summer starting in June, I should be able to do so.

Let me quickly ask: by "create our own separate project where we write all of this independently (in, say, C/C++)" you mean implementing this permutation group element there and then use it from our sage implementation, right?

I'd be happy to follow and help with than whenever you find the time, but at first I will likely only be watching you doing it...

[nthiery]

Replying to tscrim:

I'm really starting to consider that what we should do is create our own separate project where we write all of this independently (in, say, C/C++). At this stage, I'm somewhat concerned with the additional overhead that Cython could impose and the lack of complete memory control. Although I cannot commit serious time to working on this for then next two weeks (I will be in grading and math mode). Over the summer starting in June, I should be able to do so.

We can discuss this in Meudon.

I think we can avoid a bit of overhead of maintaining the GAP and category information. Although it is hard to tell how much of an impact this will have on things.

For the record: I just checked, and the data structure for permutation group elements is just a reference to the parent and a C-list of ints; plus a few slots which are unused by default (e.g. a reference to a gap element). So the only overhead is copying over the reference to the parent, and a bit of extra memory allocation

Of course, the parent itself has stuff about GAP and categories, but that's initialized once for all, and does not influence arithmetic on elements.

Cheers,

Nicolas

[tscrim]

Replying to nthiery:

Replying to tscrim:

I'm really starting to consider that what we should do is create our own separate project where we write all of this independently (in, say, C/C++). At this stage, I'm somewhat concerned with the additional overhead that Cython could impose and the lack of complete memory control. Although I cannot commit serious time to working on this for then next two weeks (I will be in grading and math mode). Over the summer starting in June, I should be able to do so.

We can discuss this in Meudon.

Sounds good.

I think we can avoid a bit of overhead of maintaining the GAP and category information. Although it is hard to tell how much of an impact this will have on things.

For the record: I just checked, and the data structure for permutation group elements is just a reference to the parent and a C-list of ints; plus a few slots which are unused by default (e.g. a reference to a gap element). So the only overhead is copying over the reference to the parent, and a bit of extra memory allocation

Of course, the parent itself has stuff about GAP and categories, but that's initialized once for all, and does not influence arithmetic on elements.

I am partially worried about how much these extra copy operations cost on the E₈ iteration scale, as well as the inherent overhead of the generated code from Cython. Plus I like having code where we completely control the memory allocations. It might end up that we really don't see much of an improvement, but I think it is worth trying.

[tscrim]

New commits:

​7b8b94b	Merge branch 'develop' into u/stumpc5/20445
​2cef8bf	Merge branch 'develop' into u/stumpc5/20445
​111d92a	Merge branch 'u/stumpc5/20445' of trac.sagemath.org:sage into public/combinat/finite_coxeter_iteration-20445

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​f061066

Doing some cleanup and now ready for review.

[chapoton]

some failing doctests, missing optional gap3

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​06e5a46

Fixing doctests.

[tscrim]

Fixed.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​d622df2	started the discriminant of a reflection group
​3d4f036	implemented the discriminant up- and downstairs
​1c82a89	Merge branch 'public/combinat/finite_coxeter_iteration-20445' of git://trac.sagemath.org/sage into 20445
​348ece0	removed a few comments

[git]

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

​d3ef5a8

removed a few comments

[stumpc5]

looked through Travis' changes...