Trac #19821: Increase speed for Coxeter groups, Weyl groups, and quan…

…tum Bruhat graph

The primary goal of this ticket is to improve the creation speed for the
quantum Bruhat graph. We do this in a number of ways:

- Better management of data associated to
`lattice.nonparabolic_positive_roots`.
- Implement a (temporary) cache of the lengths of elements.

In addition, we also provide some general speedups to all matrix groups
and Coxeter groups that came from looking into the above improvements.
The net result is over 12x speedup of the creation of the quantum Bruhat
graph:
{{{
sage: W = WeylGroup(['D',5], prefix='s')
sage: %time G = W.quantum_bruhat_graph()
CPU times: user 14 s, sys: 60.6 ms, total: 14 s
Wall time: 14 s
}}}
whereas previously this took over 3 minutes to compute. The downside is
this has a larger memory footprint because of the temporary cache, but
repeatedly computing the lengths of the elements was far too expensive.

This also includes a speedup of iterating over the entire Coxeter/Weyl
group.

URL: http://trac.sagemath.org/19821
Reported by: tscrim
Ticket author(s): Travis Scrimshaw
Reviewer(s): Frédéric Chapoton

src/sage/categories/weyl_groups.py

@@ -122,15 +122,18 @@ def pieri_factors(self, *args, **keywords):
             raise NotImplementedError("Pieri factors for type {}".format(ct))
 
         @cached_method
-        def quantum_bruhat_graph(self, index_set = ()):
+        def quantum_bruhat_graph(self, index_set=()):
             r"""
-            Returns the quantum Bruhat graph of the quotient of the Weyl group by a parabolic subgroup `W_J`.
+            Return the quantum Bruhat graph of the quotient of the Weyl
+            group by a parabolic subgroup `W_J`.
 
             INPUT:
 
-            - ``index_set`` -- a tuple `J` of nodes of the Dynkin diagram (default: ())
+            - ``index_set`` -- (default: ()) a tuple `J` of nodes of
+              the Dynkin diagram
 
-            By default, the value for ``index_set`` indicates that the subgroup is trivial and the quotient is the full Weyl group.
+            By default, the value for ``index_set`` indicates that the
+            subgroup is trivial and the quotient is the full Weyl group.
 
             EXAMPLES::
 
@@ -141,22 +144,69 @@ def quantum_bruhat_graph(self, index_set = ()):
                 sage: g.vertices()
                 [s2*s3*s1*s2, s3*s1*s2, s1*s2, s3*s2, s2, 1]
                 sage: g.edges()
-                [(s2*s3*s1*s2, s2, alpha[2]), (s3*s1*s2, s2*s3*s1*s2, alpha[1] + alpha[2] + alpha[3]),
-                (s3*s1*s2, 1, alpha[2]), (s1*s2, s3*s1*s2, alpha[2] + alpha[3]),
-                (s3*s2, s3*s1*s2, alpha[1] + alpha[2]), (s2, s1*s2, alpha[1] + alpha[2]),
-                (s2, s3*s2, alpha[2] + alpha[3]), (1, s2, alpha[2])]
+                [(s2*s3*s1*s2, s2, alpha[2]),
+                 (s3*s1*s2, s2*s3*s1*s2, alpha[1] + alpha[2] + alpha[3]),
+                 (s3*s1*s2, 1, alpha[2]),
+                 (s1*s2, s3*s1*s2, alpha[2] + alpha[3]),
+                 (s3*s2, s3*s1*s2, alpha[1] + alpha[2]),
+                 (s2, s1*s2, alpha[1] + alpha[2]),
+                 (s2, s3*s2, alpha[2] + alpha[3]),
+                 (1, s2, alpha[2])]
                 sage: W = WeylGroup(['A',3,1], prefix="s")
                 sage: g = W.quantum_bruhat_graph()
                 Traceback (most recent call last):
                 ...
-                ValueError: The Cartan type ['A', 3, 1] is not finite
+                ValueError: the Cartan type ['A', 3, 1] is not finite
             """
             if not self.cartan_type().is_finite():
-                raise ValueError("The Cartan type {} is not finite".format(self.cartan_type()))
+                raise ValueError("the Cartan type {} is not finite".format(self.cartan_type()))
+
+            # Find all the minimal length coset representatives
+            WP = [x for x in self if all(not x.has_descent(i) for i in index_set)]
+
+            # This is a modified form of quantum_bruhat_successors.
+            # It does not do any error checking and also is more efficient
+            #   with how it handles memory and checks by using data stored
+            #   at this function level rather than recomputing everything.
+            lattice = self.cartan_type().root_system().root_lattice()
+            NPR = lattice.nonparabolic_positive_roots(index_set)
+            NPR_sum = sum(NPR)
+            NPR_data = {}
+            full_NPR_sum = lattice.nonparabolic_positive_root_sum(())
+            for alpha in NPR:
+                ref = alpha.associated_reflection()
+                alphacheck = alpha.associated_coroot()
+                NPR_data[alpha] = [self.from_reduced_word(ref), # the element
+                                   len(ref) == full_NPR_sum.scalar(alphacheck) - 1, # is_quantum
+                                   NPR_sum.scalar(alphacheck)] # the scalar
+            # We also create a temporary cache of lengths as they are
+            #   relatively expensive to compute and needed frequently
+            len_cache = {}
+            def length(x):
+                # It is sufficient and faster to use the matrices as the keys
+                m = x.matrix()
+                if m in len_cache:
+                    return len_cache[m]
+                len_cache[m] = x.length()
+                return len_cache[m]
+            def succ(x):
+                w_length_plus_one = length(x) + 1
+                successors = []
+                for alpha in NPR:
+                    elt, is_quantum, scalar = NPR_data[alpha]
+                    wr = x * elt
+                    wrc = wr.coset_representative(index_set)
+                    # coset_representative returns wr if nothing gets changed
+                    if wrc is wr and length(wr) == w_length_plus_one:
+                        successors.append((wr, alpha))
+                    elif is_quantum and length(wrc) == w_length_plus_one - scalar:
+                        successors.append((wrc, alpha))
+                return successors
+
             from sage.graphs.digraph import DiGraph
-            WP = [x for x in self if x==x.coset_representative(index_set)]
-            return DiGraph([[x,i[0],i[1]] for x in WP for i in x.quantum_bruhat_successors(index_set, roots = True)],
-                           name="Parabolic Quantum Bruhat Graph of %s for nodes %s"%(self, index_set))
+            return DiGraph([[x,i[0],i[1]] for x in WP for i in succ(x)],
+                           name="Parabolic Quantum Bruhat Graph of %s for nodes %s"%(self, index_set),
+                           format="list_of_edges")
 
     class ElementMethods:
 
@@ -560,24 +610,26 @@ def inversion_arrangement(self, side='right'):
 
         def bruhat_lower_covers_coroots(self):
             r"""
-            Returns all 2-tuples (``v``, `\alpha`) where ``v`` is covered by ``self`` and `\alpha`
-            is the positive coroot such that ``self`` = ``v`` `s_\alpha` where `s_\alpha` is
+            Return all 2-tuples (``v``, `\alpha`) where ``v`` is covered
+            by ``self`` and `\alpha` is the positive coroot such that
+            ``self`` = ``v`` `s_\alpha` where `s_\alpha` is
             the reflection orthogonal to `\alpha`.
 
             ALGORITHM:
 
-            See :meth:`.bruhat_lower_covers` and :meth:`.bruhat_lower_covers_reflections` for Coxeter groups.
+            See :meth:`.bruhat_lower_covers` and
+            :meth:`.bruhat_lower_covers_reflections` for Coxeter groups.
 
             EXAMPLES::
 
                 sage: W = WeylGroup(['A',3], prefix="s")
                 sage: w = W.from_reduced_word([3,1,2,1])
                 sage: w.bruhat_lower_covers_coroots()
-                [(s1*s2*s1, alphacheck[1] + alphacheck[2] + alphacheck[3]), (s3*s2*s1, alphacheck[2]), (s3*s1*s2, alphacheck[1])]
-
+                [(s1*s2*s1, alphacheck[1] + alphacheck[2] + alphacheck[3]),
+                 (s3*s2*s1, alphacheck[2]), (s3*s1*s2, alphacheck[1])]
             """
-
-            return [(x[0],x[1].reflection_to_coroot()) for x in self.bruhat_lower_covers_reflections()]
+            return [(x[0],x[1].reflection_to_coroot())
+                    for x in self.bruhat_lower_covers_reflections()]
 
         def bruhat_upper_covers_coroots(self):
             r"""
@@ -594,41 +646,52 @@ def bruhat_upper_covers_coroots(self):
                 sage: W = WeylGroup(['A',4], prefix="s")
                 sage: w = W.from_reduced_word([3,1,2,1])
                 sage: w.bruhat_upper_covers_coroots()
-                [(s1*s2*s3*s2*s1, alphacheck[3]), (s2*s3*s1*s2*s1, alphacheck[2] + alphacheck[3]), (s3*s4*s1*s2*s1, alphacheck[4]), (s4*s3*s1*s2*s1, alphacheck[1] + alphacheck[2] + alphacheck[3] + alphacheck[4])]
-
+                [(s1*s2*s3*s2*s1, alphacheck[3]),
+                 (s2*s3*s1*s2*s1, alphacheck[2] + alphacheck[3]),
+                 (s3*s4*s1*s2*s1, alphacheck[4]),
+                 (s4*s3*s1*s2*s1, alphacheck[1] + alphacheck[2] + alphacheck[3] + alphacheck[4])]
             """
+            return [(x[0],x[1].reflection_to_coroot())
+                    for x in self.bruhat_upper_covers_reflections()]
 
-            return [(x[0],x[1].reflection_to_coroot()) for x in self.bruhat_upper_covers_reflections()]
-
-        def quantum_bruhat_successors(self, index_set = None, roots = False, quantum_only = False):
+        def quantum_bruhat_successors(self, index_set=None, roots=False, quantum_only=False):
             r"""
-            Returns the successors of ``self`` in the parabolic quantum Bruhat graph.
+            Return the successors of ``self`` in the quantum Bruhat graph
+            on the parabolic quotient of the Weyl group determined by the
+            subset of Dynkin nodes ``index_set``.
 
             INPUT:
 
-            - ``self`` -- a Weyl group element, which is assumed to be of minimum length in its coset with respect to the parabolic subgroup
+            - ``self`` -- a Weyl group element, which is assumed to
+              be of minimum length in its coset with respect to the
+              parabolic subgroup
 
-            - ``index_set`` -- (default: None) indicates the set of simple reflections used to generate the parabolic subgroup;
-               the default value indicates that the subgroup is the identity
+            - ``index_set`` -- (default: ``None``) indicates the set of
+              simple reflections used to generate the parabolic subgroup;
+              the default value indicates that the subgroup is the identity
 
-            - ``roots`` -- (default: False) if True, returns the list of 2-tuples (``w``, `\alpha`) where ``w`` is a
-               successor and `\alpha` is the positive root associated with the successor relation.
+            - ``roots`` -- (default: ``False``) if ``True``, returns the
+              list of 2-tuples (``w``, `\alpha`) where ``w`` is a successor
+              and `\alpha` is the positive root associated with the
+              successor relation
 
-            - ``quantum_only`` -- (default: False) if True, returns only the quantum successors
-
-            Returns the successors of ``self`` in the quantum Bruhat graph on the parabolic
-            quotient of the Weyl group determined by the subset of Dynkin nodes ``index_set``.
+            - ``quantum_only`` -- (default: ``False``) if ``True``, returns
+              only the quantum successors
 
             EXAMPLES::
 
                 sage: W = WeylGroup(['A',3], prefix="s")
                 sage: w = W.from_reduced_word([3,1,2])
                 sage: w.quantum_bruhat_successors([1], roots = True)
-                [(s3, alpha[2]), (s1*s2*s3*s2, alpha[3]), (s2*s3*s1*s2, alpha[1] + alpha[2] + alpha[3])]
+                [(s3, alpha[2]), (s1*s2*s3*s2, alpha[3]),
+                 (s2*s3*s1*s2, alpha[1] + alpha[2] + alpha[3])]
                 sage: w.quantum_bruhat_successors([1,3])
                 [1, s2*s3*s1*s2]
                 sage: w.quantum_bruhat_successors(roots = True)
-                [(s3*s1*s2*s1, alpha[1]), (s3*s1, alpha[2]), (s1*s2*s3*s2, alpha[3]), (s2*s3*s1*s2, alpha[1] + alpha[2] + alpha[3])]
+                [(s3*s1*s2*s1, alpha[1]),
+                 (s3*s1, alpha[2]),
+                 (s1*s2*s3*s2, alpha[3]),
+                 (s2*s3*s1*s2, alpha[1] + alpha[2] + alpha[3])]
                 sage: w.quantum_bruhat_successors()
                 [s3*s1*s2*s1, s3*s1, s1*s2*s3*s2, s2*s3*s1*s2]
                 sage: w.quantum_bruhat_successors(quantum_only = True)
@@ -638,11 +701,10 @@ def quantum_bruhat_successors(self, index_set = None, roots = False, quantum_onl
                 Traceback (most recent call last):
                 ...
                 ValueError: s2*s3 is not of minimum length in its coset of the parabolic subgroup generated by the reflections (1, 3)
-
             """
             W = self.parent()
             if not W.cartan_type().is_finite():
-                raise ValueError("The Cartan type {} is not finite".format(W.cartan_type()))
+                raise ValueError("the Cartan type {} is not finite".format(W.cartan_type()))
             if index_set is None:
                 index_set = []
             else:

src/sage/combinat/root_system/weyl_group.py

@@ -711,11 +711,10 @@ def __init__(self, parent, g, check=False):
             sage: TestSuite(s1).run()
         """
         MatrixGroupElement_gap.__init__(self, parent, g, check=check)
-        self.__matrix = self.matrix()
         self._parent = parent
 
     def __hash__(self):
-        return hash(self.__matrix)
+        return hash(self.matrix())
 
     def domain(self):
         """
@@ -803,7 +802,7 @@ def __eq__(self, other):
         """
         return self.__class__ == other.__class__ and \
                self._parent   == other._parent   and \
-               self.__matrix  == other.__matrix
+               self.matrix()  == other.matrix()
 
     def _cmp_(self, other):
         """
@@ -850,7 +849,7 @@ def action(self, v):
         """
         if v not in self.domain():
             raise ValueError("{} is not in the domain".format(v))
-        return self.domain().from_vector(self.__matrix*v.to_vector())
+        return self.domain().from_vector(self.matrix()*v.to_vector())
 
 
     ##########################################################################
@@ -859,12 +858,14 @@ def action(self, v):
 
     def has_descent(self, i, positive=False, side = "right"):
         """
-        Tests if self has a descent at position `i`, that is if self is
+        Test if ``self`` has a descent at position ``i``.
+
+        An element `w` has a descent in position `i` if `w` is
         on the strict negative side of the `i^{th}` simple reflection
         hyperplane.
 
-        If positive is True, tests if it is on the strict positive
-        side instead.
+        If ``positive`` is ``True``, tests if it is on the strict
+        positive side instead.
 
         EXAMPLES::
 
@@ -921,33 +922,54 @@ def has_descent(self, i, positive=False, side = "right"):
             self = ~self
 
         if use_rho:
-            s = self.action(L.rho()   ).scalar(L.alphacheck()[i]) >= 0
+            s = self.action(L.rho()).scalar(L.alphacheck()[i]) >= 0
         else:
             s = self.action(L.alpha()[i]).is_positive_root()
 
         return s is positive
 
-    def has_left_descent(self,i):
+    def has_left_descent(self, i):
         """
-        Tests if self has a left descent at position `i`.
+        Test if ``self`` has a left descent at position ``i``.
 
         EXAMPLES::
 
             sage: W = WeylGroup(['A',3])
             sage: s = W.simple_reflections()
-            sage: [W.one().has_descent(i) for i in W.domain().index_set()]
+            sage: [W.one().has_left_descent(i) for i in W.domain().index_set()]
             [False, False, False]
-            sage: [s[1].has_descent(i) for i in W.domain().index_set()]
+            sage: [s[1].has_left_descent(i) for i in W.domain().index_set()]
             [True, False, False]
-            sage: [s[2].has_descent(i) for i in W.domain().index_set()]
+            sage: [s[2].has_left_descent(i) for i in W.domain().index_set()]
             [False, True, False]
-            sage: [s[3].has_descent(i) for i in W.domain().index_set()]
+            sage: [s[3].has_left_descent(i) for i in W.domain().index_set()]
+            [False, False, True]
+            sage: [(s[3]*s[2]).has_left_descent(i) for i in W.domain().index_set()]
             [False, False, True]
-            sage: [s[3].has_descent(i, True) for i in W.domain().index_set()]
-            [True, True, False]
         """
         return self.has_descent(i, side = "left")
 
+    def has_right_descent(self, i):
+        """
+        Test if ``self`` has a right descent at position ``i``.
+
+        EXAMPLES::
+
+            sage: W = WeylGroup(['A',3])
+            sage: s = W.simple_reflections()
+            sage: [W.one().has_right_descent(i) for i in W.domain().index_set()]
+            [False, False, False]
+            sage: [s[1].has_right_descent(i) for i in W.domain().index_set()]
+            [True, False, False]
+            sage: [s[2].has_right_descent(i) for i in W.domain().index_set()]
+            [False, True, False]
+            sage: [s[3].has_right_descent(i) for i in W.domain().index_set()]
+            [False, False, True]
+            sage: [(s[3]*s[2]).has_right_descent(i) for i in W.domain().index_set()]
+            [False, True, False]
+        """
+        return self.has_descent(i, side="right")
+
     def apply_simple_reflection(self, i, side = "right"):
         s = self.parent().simple_reflections()
         if side == "right":