Stanley symmetric functions for type A Weyl Group

[trevorkarn]

The Stanley symmetric function was implemented in several types in #8810.

A theorem of Reiner-Shimozono (1995) allows us to compute these more efficiently in type-A.

This ticket implements that algorithm.

Depends on #34260 Depends on #34343 Depends on #34339

CC: @tscrim @anneschilling @nthiery

Component: combinatorics

Keywords: stanley symmetric-function weyl-group gsoc2022

Author: Trevor K. Karn

Branch/Commit: u/tkarn/stanley-sym-fcn-34335 @ a2431a2

Issue created by migration from https://trac.sagemath.org/ticket/34335

[trevorkarn]

comment:1

sage: W = WeylGroup(['A',4])
sage: w = W.from_reduced_word([1, 2, 3, 4, 1, 2, 3, 1, 2, 1])
sage: %prun w.stanley_symmetric_function()

    11225971 function calls (10707894 primitive calls) in 21.946 seconds

Compared to:

sage: from sage.combinat.diagram import RotheDiagram
sage: sage: w = Permutations(5).from_reduced_word([1, 2, 3, 4, 1, 2, 3, 1, 2, 1])
....: sage: Sym = SymmetricFunctions(QQ)
....: sage: s = Sym.s()
....: sage: m = Sym.m()
....: sage: %prun sum( s[T.shape()] for T in RotheDiagram(w).peelable_tableaux())

         9571 function calls (9082 primitive calls) in 0.014 seconds

   Ordered by: internal time

[trevorkarn]

comment:2

This approach is not universally faster

sage: W = WeylGroup(['A',5])
sage: w = W.from_reduced_word([1, 2, 3, 4, 2, 3, 2, 1])
sage: %timeit w.stanley_symmetric_function()
261 µs ± 81.5 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)

sage: from sage.combinat.diagram import RotheDiagram
sage: Sym = SymmetricFunctions(QQ)
sage: p = Permutations(6).from_reduced_word([1, 2, 3, 4, 2, 3, 2, 1])
sage: s = Sym.s()
sage: %timeit s.sum( s[T.shape()]  for T in RotheDiagram(p).peelable_tableaux())
2.21 ms ± 113 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

[tscrim]

comment:3

Here is the result of

sage: %lprun -f RotheDiagram(p).peelable_tableaux s.sum(s[T.shape()] for T in RotheDiagram(p).peelable_tableaux())

for the example in comment:2:

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
   801         4          2.0      0.5      0.1          if self._n_nonempty_cols == 1:
   802         1         19.0     19.0      0.8              return {Tableau([[i+1] for i, j in self._cells])}
   803                                           
   804         3          7.0      2.3      0.3          first_col = min(j for i, j in self._cells)
   805                                           
   806                                                   # TODO: The next two lines of code could be optimized by only
   807                                                   # looping over self._cells once (rather than two separate times)
   808                                                   # get the diagram without the first column
   809         3         59.0     19.7      2.4          Dhat = NorthwestDiagram([c for c in self._cells if c[1] != first_col])
   810                                           
   811                                                   # get the values from the first column
   812         3          7.0      2.3      0.3          new_vals = sorted(i + 1 for i, j in self._cells if j == first_col)
   813                                           
   814         3          4.0      1.3      0.2          k = self.n_cells() - Dhat.n_cells()
   815                                           
   816         3          2.0      0.7      0.1          peelables = set()
   817                                           
   818         6          4.0      0.7      0.2          for Q in Dhat.peelable_tableaux():
   819                                                       # get the vertical strips
   820         3         50.0     16.7      2.0              mu = Q.shape()
   821         3        918.0    306.0     37.4              vertical_strips = mu.add_vertical_border_strip(k)
   822        12          8.0      0.7      0.3              for s in vertical_strips:
   823         9         72.0      8.0      2.9                  sQ = SkewTableaux()(Q)  # sQ is skew - get it?
   824         9        283.0     31.4     11.5                  new_cells = (s/mu).cells()
   825                                                           # perform the jeu de taquin slides
   826        32         20.0      0.6      0.8                  for c in new_cells:
   827        23        813.0     35.3     33.1                      sQ = sQ.backward_slide(c)
   828                                                           # create the new tableau by filling the columns
   829         9         15.0      1.7      0.6                  sQ_new = sQ.to_list()
   830        32         21.0      0.7      0.9                  for n in range(k):
   831        23         16.0      0.7      0.7                      sQ_new[n][0] = new_vals[n]
   832                                           
   833         9         68.0      7.6      2.8                  T = Tableau(sQ_new)
   834         9         57.0      6.3      2.3                  if T.is_column_strict():
   835         3          5.0      1.7      0.2                      peelables.add(T)
   836                                           
   837         3          3.0      1.0      0.1          return peelables

This suggests that you should optimize the Partition.add_vertical_border_strip() and/or backward_slide(). For the add_vertical_border_strip(), the implementation calls the horizontal version (which also calls the conjugate()!) and does a lot of conjugation. Actually, it seems like the implementations should be swapped: The add_horizontal_border_strip() is working with the conjugate shape. This should be a separate ticket though.

[tscrim]

comment:4

This is now #34339.

[trevorkarn]

comment:5

sage: set_random_seed(12)
....: 
....: P = Partitions(100);
....: mu = P.random_element();
....: len(mu)
18
sage: T = SemistandardTableaux(mu).random_element();  T.pp()
  1  7   8 12 14 15 16 22 28 28 30 32 40 42 55 58 64 72 89 91
  2 10  11 16 20 22 27 27 30 31 36 55 57 57 59 59 93
  4 11  22 25 27 36 44 47 52 52 52 56 61 98
  5 13  23 37 39 50 57 59 62 79 93 93 99
 10 23  40 48 52 63 72 77 83
 23 27  43 78 88
 25 58  69 89
 34 59  80
 35 79  98
 40 80 100
 43 81
 47
 61
 67
 71
 75
 77
 88
sage: S = SkewTableau(T); S.pp()
  1  7  8 12 14 15 16 22 28 28 30 32 40 42 55 58 64 72 89 91
  2 10 11 16 20 22 27 27 30 31 36 55 57 57 59 59 93
  4 11 22 25 27 36 44 47 52 52 52 56 61 98
  5 13 23 37 39 50 57 59 62 79 93 93 99
 10 23 40 48 52 63 72 77 83
 23 27 43 78 88
 25 58 69 89
 34 59 80
 35 79 98
 40 80100
 43 81
 47
 61
 67
 71
 75
 77
 88
sage: %lprun -f  S.backward_slide S.backward_slide((5,5))

Timer unit: 1e-06 s

Total time: 0.001199 s
File: /Users/trevorkarn/Applications/sage/src/sage/combinat/skew_tableau.py
Function: backward_slide at line 911

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
   911                                               def backward_slide(self, corner=None):                                        
   ...
   988         1         35.0     35.0      2.9          new_st = self.to_list()
   989         1        733.0    733.0     61.1          inner_outside_corners = self.inner_shape().outside_corners()
   990         1        112.0    112.0      9.3          outer_outisde_corners = self.outer_shape().outside_corners()
   991         1          1.0      1.0      0.1          if corner is not None:
   992         1          4.0      4.0      0.3              if tuple(corner) not in outer_outisde_corners:
   993                                                           raise ValueError("corner must be an outside corner \
   994                                                                             of the outer shape")
   995                                                   else:
   996                                                       if not outer_outisde_corners:
   997                                                           return self
   998                                                       else:
   999                                                           corner = outer_outisde_corners[0]
  1000                                           
  1001         1          0.0      0.0      0.0          i, j = corner
  1002                                           
  1003                                                   # add the empty cell
  1004                                                   # the column only matters if it is zeroth column, in which
  1005                                                   # case we need to add a new row.
  1006         1          0.0      0.0      0.0          if not j:
  1007                                                       new_st.append(list())
  1008         1          2.0      2.0      0.2          new_st[i].append(None)
  1009                                           
  1010                                           
  1011                                           
  1012        11          8.0      0.7      0.7          while (i, j) not in inner_outside_corners:
  1013                                                       # get the value of the cell above the temporarily empty cell (if
  1014                                                       # it exists)
  1015        10          8.0      0.8      0.7              if i > 0:
  1016         9         10.0      1.1      0.8                  P_up = new_st[i-1][j]
  1017                                                       else:
  1018         1          1.0      1.0      0.1                  P_up = -1  # a dummy value less than all positive numbers
  1019                                           
  1020                                                       # get the value of the cell to the left of the temp. empty cell
  1021                                                       # (if it exists)
  1022        10         10.0      1.0      0.8              if j > 0:
  1023        10         10.0      1.0      0.8                  P_left = new_st[i][j-1]
  1024                                                       else:
  1025                                                           P_left = -1  # a dummy value less than all positive numbers
  1026                                           
  1027                                                       # get the next cell
  1028                                                       # p_left will always be positive, but if P_left
  1029                                                       # and P_up are both None, then it will return
  1030                                                       # -1, which doesn't trigger the conditional
  1031        10          9.0      0.9      0.8              if P_left > P_up:
  1032         5          4.0      0.8      0.3                  new_st[i][j] = P_left
  1033         5          2.0      0.4      0.2                  i, j = (i, j - 1)
  1034                                                       else:  # if they are equal, we slide up
  1035         5          2.0      0.4      0.2                  new_st[i][j] = P_up
  1036         5          4.0      0.8      0.3                  i, j = (i - 1, j)
  1037                                                   # We don't need to reset the intermediate cells inside the loop
  1038                                                   # because the conditional above will continue to overwrite it until
  1039                                                   # the while loop terminates. We do need to reset it at the end.
  1040         1          1.0      1.0      0.1          new_st[i][j] = None
  1041                                           
  1042         1        243.0    243.0     20.3          return SkewTableau(new_st)

[trevorkarn]

Changed dependencies from #34260 to #34260 #34343

[trevorkarn]

Changed dependencies from #34260 #34343 to #34260 #34343 #34339

[trevorkarn]

Branch: u/tkarn/stanley-sym-fcn-34335

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. Last 10 new commits:

6b7906a	Fix __init__
6d1395e	Fix kwargs, repr, and remove double .__init__ call
630017d	Fix accidental check=False
b7970d7	Add references and more examples. All doctests pass
c821a3e	Fix a typo in Rothe diagrams
0473f20	Fix bug in identity partition/empty diagram
cfee8e8	Speed up sliding by skipping classcall
62c7b5f	Remove unneeded creation of skew tableaux
45c372c	Change to list
8d1500d	Add type-A algorithm

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Commit: 8d1500d

[trevorkarn]

comment:10

One doctest fails, giving an incorrect result for the SymmetricGroupElement example. In further investigating, I see that:

sage: G = SymmetricGroup(4)
sage: g = G.from_reduced_word([3,2,3,1])
sage: W = WeylGroup(['A', 3]);
sage: w = W.from_reduced_word([3,2,3,1])
sage: g
(1,2,4)
sage: w
[0 1 0 0]
[0 0 0 1]
[0 0 1 0]
[1 0 0 0]
sage: g.domain()
[2, 4, 3, 1]
sage: w.to_permutation()
(4, 1, 3, 2)

which to me is surprising that they are different? Is this correct? When I do the computation by hand I agree with the Weyl group result for w.

[tscrim]

comment:11

One needs to be careful of multiplication conventions of permutations. That would be my guess why they are different.

[trevorkarn]

comment:12

Replying to @tscrim:

One needs to be careful of multiplication conventions of permutations. That would be my guess why they are different.

Indeed:

sage: SymmetricGroup(4).from_reduced_word([1,3,2,3]).domain()
[4, 1, 3, 2]
sage: W = WeylGroup(['A', 3])
sage: W.from_reduced_word([3,2,3,1]).to_permutation()
(4, 1, 3, 2)

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

8b6ef8a	First round of reviewer suggestions
e973f7d	Assertion -> ValueError
a57b2d0	Add tests and remove Stanley s.f.
28286da	Add methods to interact with polyominos/compositions
75ce24e	Change to list
8c3cb4e	Add type-A algorithm
4a6b2f9	Rewrite to pass to Permutations

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 8d1500d to 4a6b2f9

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

e1f113c	Fix bug in identity partition/empty diagram
798df7a	Speed up sliding by skipping classcall
9227a3c	Remove unneeded creation of skew tableaux
fcb9368	First round of reviewer suggestions
87bbc48	Assertion -> ValueError
766584b	Add tests and remove Stanley s.f.
170c657	Add methods to interact with polyominos/compositions
11bb9df	Add rothe_diagram, n_reduced_words, and Stanley symmetric function
ac13e8e	Add type-A algorithm
30b2471	Rewrite to pass to Permutations

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 4a6b2f9 to 30b2471

[trevorkarn]

comment:16

Rebase off of #34260

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

1eea8bf	Remove `_hash_`, `__contains__`, `_repr_` overrides
1a7eea6	Fix peelable tableaux documentation
c55f72e	Redo loop to avoid looping over self._cells twice
252eafd	Move zero/one inside check
fd97067	Fix pyflakes issue
8873753	Fix pyflakes issue
81a859f	Add .from_* methods to __element_constructor__
2a4a374	Add rothe_diagram, n_reduced_words, and Stanley symmetric function
80eda06	Add type-A algorithm
0c65d4e	Rewrite to pass to Permutations

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 30b2471 to 0c65d4e

[trevorkarn]

comment:18

Replying to @sagetrac-git:

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

1eea8bf	Remove `_hash_`, `__contains__`, `_repr_` overrides
1a7eea6	Fix peelable tableaux documentation
c55f72e	Redo loop to avoid looping over self._cells twice
252eafd	Move zero/one inside check
fd97067	Fix pyflakes issue
8873753	Fix pyflakes issue
81a859f	Add .from_* methods to __element_constructor__
2a4a374	Add rothe_diagram, n_reduced_words, and Stanley symmetric function
80eda06	Add type-A algorithm
0c65d4e	Rewrite to pass to Permutations

Rebased off of rc0 and #34510

[trevorkarn]

comment:19

Fix algorithm block and "type A" not "type-A".

[trevorkarn]

comment:20

Add catch for reduced words.

[trevorkarn]

comment:21

Replying to Trevor Karn:

Add catch for reduced words.

Actually - fix category structure. Try to include Coxeter group elements.

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 0c65d4e to 41d5ed2

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

f079114	Add diagram iterator
5f40788	Add tests and add a catch for subclasses
0f6dfd9	PositiveIntegers -> NonNegativeIntegers
2d1c3cd	Add rothe_diagram, n_reduced_words, and Stanley symmetric function
0a10a5d	Add type-A algorithm
1febc72	Rewrite to pass to Permutations
7a8e64f	Move from category to element class
41d5ed2	Fix basis and add to coxeter group

[trevorkarn]

comment:25

Should the WeylGroups.ElementMethods.stanley_symmetric_function and *_as_polynomial methods be moved to the category of CoxeterGroups? My read of SymmetricFunctions.com is "yes" but that involves more. Would it be OK to leave Coxeter groups off of this ticket and make it a follow up one?

[tscrim]

comment:26

If there is a difference between types B and C, then they cannot move up to Coxeter groups. I am not sure without reading more closely.

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

8c678bb

Remove from Coxeter Groups

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 41d5ed2 to 8c678bb

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 8c678bb to fdaede5

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

fdaede5

Fix author typos

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from fdaede5 to 97d3afb

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

97d3afb

Add examples and fix bug

[trevorkarn]

comment:31

Replying to Travis Scrimshaw:

If there is a difference between types B and C, then they cannot move up to Coxeter groups. I am not sure without reading more closely.

Ok got it. I think this is ready.

[tscrim]

comment:32

Thanks. Two things:

1. Why do you convert the element to the monomial symmetric functions? Is this for consistency with the other types? This seems a bit wasteful (in particular, if you wanted it in the Schur basis). Perhaps we should have this function take a symmetric function basis as an argument that yields the output (with the default being SymmetricFunctions(QQ).m() as it is now).

2. Move the _tableau_contribution into the n_reduced_words method and have it be

   from sage.combinat.tableau import StandardTableaux
   def _tableau_contribution(T):
      r"""
      Get the number of SYT of shape(``T``).
      """    
      return StandardTableaux(T.shape()).cardinality()

   Note that the import is outside of the call (so it isn't done repeatedly every time).

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 97d3afb to feb7d4c

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

c6c9461	Add type-A algorithm
d7be408	Rewrite to pass to Permutations
da0489a	Move from category to element class
37ee809	Fix basis and add to coxeter group
19f2375	Remove from Coxeter Groups
929fef3	Fix author typos
4ccc000	Add examples and fix bug
1692d9e	Rename n_reduced_words -> number_of_reduced_words
2d6b8cb	Move _tableau_contribution
feb7d4c	Add basis option for stanley symmetric function

[trevorkarn]

comment:34

Replying to Travis Scrimshaw:

Thanks. Two things:

1. Why do you convert the element to the monomial symmetric functions? Is this for consistency with the other types? This seems a bit wasteful (in particular, if you wanted it in the Schur basis). Perhaps we should have this function take a symmetric function basis as an argument that yields the output (with the default being SymmetricFunctions(QQ).m() as it is now).

2. Move the _tableau_contribution into the n_reduced_words method and have it be

   from sage.combinat.tableau import StandardTableaux
   def _tableau_contribution(T):
      r"""
      Get the number of SYT of shape(``T``).
      """    
      return StandardTableaux(T.shape()).cardinality()

   Note that the import is outside of the call (so it isn't done repeatedly every time).

Done!

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from feb7d4c to 0316612

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

0316612

Use sum_of_monomials for some speedup

[anneschilling]

comment:36

Thank you for the implementation, Trevor! It was good to see you in Poland this summer and talk to you about this ticket!

One suggestion I have is to allow the user the choice to use the old implementation of the Stanley symmetric function in type A (non affine). You can have your new implementation as the default, but it would be good to still have the ability to access the old method. That way you could also put in some tests that check that the old and new implementation give the same answer!

[tscrim]

comment:37

+1 on Anne's suggestion. Perhaps algorithm=None for the general case, which is typically ignored. Then in type A,

if algorithm is None:
    algorithm = "PT"
if algorithm = "PT":
    return peelable_tableaux
if algorithm = "Pieri":  # if not in the category, otherwise just fall through
    return super().stanley_symmetric_function(...)

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 0316612 to a5e0bd0

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

7e4e524	Add rothe_diagram, n_reduced_words, and Stanley symmetric function
d74232f	Add type-A algorithm
427b392	Rewrite to pass to Permutations
cab87ca	Move from category to element class
d3a4577	Fix basis and add to coxeter group
8f0cc0a	Remove from Coxeter Groups
e37448e	Fix author typos
6056efd	Add examples and fix bug
a66c597	Rename n_reduced_words -> number_of_reduced_words
a5e0bd0	Move _tableau_contribution

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

37692d8	Add type-A algorithm
60e31ab	Rewrite to pass to Permutations
c075aef	Move from category to element class
1cc9b50	Fix basis and add to coxeter group
9bfa258	Remove from Coxeter Groups
b3144bc	Fix author typos
d3f8a3f	Add examples and fix bug
2f47b79	Rename n_reduced_words -> number_of_reduced_words
cdda33c	Move _tableau_contribution
4e0b3dd	Add basis option for stanley symmetric function

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from a5e0bd0 to 4e0b3dd

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

49022ad	Add rothe_diagram, n_reduced_words, and Stanley symmetric function
822906d	Add type-A algorithm
da18880	Rewrite to pass to Permutations
6ed20d4	Move from category to element class
2dcfe9c	Fix basis and add to coxeter group
a1f8e8a	Remove from Coxeter Groups
6f729f4	Add examples and fix bug
04bc80a	Rename n_reduced_words -> number_of_reduced_words
a620662	Move _tableau_contribution
95450f1	Add basis option for stanley symmetric function

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 4e0b3dd to 95450f1

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

82d3587	Add type-A algorithm
58dd99c	Rewrite to pass to Permutations
c7d83d3	Move from category to element class
bf055fb	Fix basis and add to coxeter group
25af462	Remove from Coxeter Groups
ebd122b	Fix author typos
902f784	Add examples and fix bug
d7f0d8e	Rename n_reduced_words -> number_of_reduced_words
8d0192f	Move _tableau_contribution
a2431a2	Add algorithm options