speedup central orthogonal idempotents for SymmetricGroupAlgebra

[MeeSeongIm]

The current implementation of central_orthogonal_idempotents is inherited from FiniteDimensionalSemisimpleAlgebrasWithBasis. The combinatorial formula in terms of standard tableaux is much faster (implemented as the method cpis).

This ticket does two/four things:

• replaces the default implementation of central_orthogonal_idempotent/s with one based on the code in cpi/s;

• fixes a bug in cpi/s, which assumes base ring contains QQ;

• implements primitive central orthogonal idempotents in any characteristic, using Nakayama's Conjecture;

• adds a deprecated_function_alias for the nondescriptive cpi/s methods.

CC: @tscrim @nthiery @alauve @canozanoguz @sagetrac-meeseongim @srdoty

Component: combinatorics

Keywords: sagedays@icerm, central idempotents

Author: Mee Seong Im, Aaron Lauve

Branch/Commit: 77c38e2

Reviewer: Travis Scrimshaw

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

[alauve]

New commits:

78de7b5

overwrote default implementation of orthogonal idempotents with the faster cpis method.

[alauve]

Branch: public/combinat/speedup_idempotents_for_SymmetricGroupAlgebras_25942

[alauve]

Commit: 78de7b5

[alauve]

comment:2

The current implementation (develop branch) comes from FiniteDimensionalSemisimpleAlgebrasWithBasis and is quite slow. E.g.,

sage: %time z = SymmetricGroupAlgebra(QQ, 5).central_orthogonal_idempotents()
CPU times: user 3.04 s, sys: 65.8 ms, total: 3.11 s
Wall time: 3.11 s

On this branch's ticket, the method central_orthogonal_idempotents just points to cpis, which is much faster:

sage: %time z = SymmetricGroupAlgebra(QQ, 5).central_orthogonal_idempotents()
CPU times: user 177 ms, sys: 38.4 ms, total: 215 ms
Wall time: 608 ms

[56b43a0c-d98d-40ca-8d50-ad4f2594114b]

comment:3

Minor change. Needs review.

[tscrim]

comment:4

There is a test failure:

sage -t src/sage/categories/finite_dimensional_semisimple_algebras_with_basis.py
**********************************************************************
File "src/sage/categories/finite_dimensional_semisimple_algebras_with_basis.py", line 89, in sage.categories.finite_dimensional_semisimple_algebras_with_basis.FiniteDimensionalSemisimpleAlgebrasWithBasis.ParentMethods.central_orthogonal_idempotents
Failed example:
    idempotents
Expected:
    (1/6*() + 1/6*(2,3) + 1/6*(1,2) + 1/6*(1,2,3) + 1/6*(1,3,2) + 1/6*(1,3),
     2/3*() - 1/3*(1,2,3) - 1/3*(1,3,2),
     1/6*() - 1/6*(2,3) - 1/6*(1,2) + 1/6*(1,2,3) + 1/6*(1,3,2) - 1/6*(1,3))
Got:
    [1/6*() + 1/6*(2,3) + 1/6*(1,2) + 1/6*(1,2,3) + 1/6*(1,3,2) + 1/6*(1,3),
     2/3*() - 1/3*(1,2,3) - 1/3*(1,3,2),
     1/6*() - 1/6*(2,3) - 1/6*(1,2) + 1/6*(1,2,3) + 1/6*(1,3,2) - 1/6*(1,3)]
**********************************************************************

The test in the file is no longer actually testing the method there. So you need a different semisimple finite-dimensional algebra with a (distinguished) basis.

Also, for consistency (and speed), you should return a tuple for cpis and make it an @cached_method.

[tscrim]

Reviewer: Travis Scrimshaw

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

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

d646607	added caching; fixed an error in coefficient computation within cpi method
7e6c7ef	work on docstrings

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

Changed commit from 78de7b5 to 7e6c7ef

[alauve]

comment:6

What is the correct way to check for semisimplicity?

(To maintain the old behavior, the method central_orthogonal_idempotents shouldn't appear for the user in the bad cases.)

[tscrim]

comment:7

Replying to @alauve:

What is the correct way to check for semisimplicity?

You will need some form of a mathematical test, but the easiest way is:

sage: S3 = SymmetricGroupAlgebra(QQ, 3)
sage: S3 in Algebras.Semisimple
True
sage: S3p = SymmetricGroupAlgebra(GF(3), 3)
sage: S3p in Algebras.Semisimple
False

(To maintain the old behavior, the method central_orthogonal_idempotents shouldn't appear for the user in the bad cases.)

That is way too difficult and complicated. Just raise a ValueError when it is not semisimple in cpis.

[MeeSeongIm]

comment:8

User gh-canozanoguz added.

[56b43a0c-d98d-40ca-8d50-ad4f2594114b]

comment:9

Replying to @alauve:

What is the correct way to check for semisimplicity?

An algebra is semisimple if the Jacobson radical of the algebra is trivial, and an element x in R is in the Jacobson radical of R if 1 + RxR is a unit, or rx is right quasi-regular for all r in R (rx is right quasi-regular if there exists an element z in R such that rx + z + rxz = 0).

If R has an identity element, then rx is right quasi-regular if and only if 1 + rx has a right inverse in R.

Here is an alternative way to check for semisimplicity:

sage: S3 = SymmetricGroupAlgebra(QQ,3)
sage: S3.radical()
sage: S3.radical_basis()
()
sage: S3p = SymmetricGroupAlgebra(GF(3), 3)
sage: S3p.radical()
sage: S3p.radical_basis()
([1, 2, 3] + 2*[3, 1, 2],
 [1, 3, 2] + 2*[3, 2, 1],
 [2, 1, 3] + 2*[3, 2, 1],
 [2, 3, 1] + 2*[3, 1, 2])

[tscrim]

comment:10

I think Aaron meant more in terms of code. Mathematically, the radical is the way to do it, but when it is known by other (computationally easier)- methods (since this is a group algebra), it is better to do those. Computing the radical can be computationally intensive, so is something generally to be avoided (at least by default).

[alauve]

comment:11

That's right. I was speaking less about "how best" rather than "how".

I don't know how the semisimplicity check gets done...

I don't declare the axiom "Semisimple" at any point when I run

sage: S3 = SymmetricGroupAlgebra(QQ, 3)

versus

sage: S3mod = SymmetricGroupAlgebra(GF(3), 3)

yet somehow Sage knows it.

When I browse through the combinat.symmetric_group_algebra.py file, I don't see any mention of semisimple. (Nor in the algebras.group_algebra.py file, for that matter.) Yet---in the development branch, at least---the method .central_orthogonal_idempotents is not made available to the user for S3mod but it is made available for S3.

[alauve]

comment:12

I've made two more changes that I'd like an experienced contributor to chime-in on:

Change No.1/

The method .cpis was not cached. Okay, but probably it should be, since .central_orthogonal_idempotents was and we're replacing it. However, .cpis simply calls .cpi for each partition of n. It seems more natural to cache these.

That way, when the user wants to study specht modules of hook-shape, or some such (i.e., not ever calling .cpis), then s/he would not have to repeat expensive computations.

Problem: .cpi is happy to receive a list as argument (instead of a proper Partition). But you cannot add @cached_method in this case. So, do I:

• cache individual idempotents to be friendly to Professor Hook-Shape, or
• allow lists as inputs to be friendly to Professor Hook-Shape?

Change No.2/

The .cpi code was incorrect, creating expressions with coefficients in QQ even if the group algebra was defined over GF(3). I have added a check that denominators are not zero, and separately coerced the coefficients into self.base_ring(). (Note that some idempotents are well-defined even when the ring is not semisimple, so this seems like a better solution than just disallowing use of .cpi in these cases.)

[nthiery]

comment:13

That's right. I was speaking less about "how best" rather than "how".

In addition, there is the distinction between:

• Is -- mathematically speaking -- the algebra semisimple?
• Does Sage know that the algebra is semi-simple?

Travis' suggestion is about the latter, and this is was we want in such situations: based on the current knowledge which algorithm is it best to use.

I don't know how the semisimplicity check gets done...

The relevant statement is that the algebra of a finite group is semisimple if the characteristic does not divide the order of the group, right? This is thus implemented in Groups.Finite.Algebra. And since one needs to have information about the order of the group, this goes in Groups.Finite.Algebra.ParentMethods.__init_extra__ which is called on any group algebra of a finite group upon it's initialization.

[alauve]

comment:14

From what you say, Nicolas, I gather this means that one way I could successsfully not provide central_orthogonal_idempotents in the bad cases is by moving my definition (see the one-line change in commit:​78de7b5) to the end of __init__, wrapped in an if statement?:

if self in Algebras.Semisimple:
    self.central_orthogonal_idempotents = self.cpis

[nthiery]

comment:15

Hmm, interesting use case. As long as this is after the Parent.__init__ (at the end of which the __init_extra__ methods are called), this would indeed work.

This is polluting a bit the __init__ method though. Also we would want to deprecate cpis as well.

One day, we really should have this @conditionally_defined attribute.

For now, I would:

• mv cpis to central_orthogonal_idempotents
• make cpi a deprecated_function_alias
• raise a ValueError at the beginning of central_orthogonal_idempotents if the group algebra is not semisimple

[alauve]

comment:16

Do you mean to make cpis a deprecated_function_alias or cpi?

I'm inclined to leave cpi available (with new name or no), as one can sometimes build these idempotents even when the algebra is not semisimple. (Of course the "p"rimitive part is dubious in that case.)

[nthiery]

comment:17

Replying to @alauve:

Do you mean to make cpis a deprecated_function_alias or cpi? I'm inclined to leave cpi available (with new name or no), as one can sometimes build these idempotents even when the algebra is not semisimple. (Of course the "p"rimitive part is dubious in that case.)

Agreed: deprecating both cpis and cpi, while renaming the later to something relevant; maybe a name xxx that suggests how it is constructed, which makes it meaningful even when the algebra is not semisimple.

central_orthogonal_idempotents could then return a LazyFamily(Partitions(n), self.xxx) in the semi-simple case.

[nthiery]

comment:18

(with presumably a cache on xxx; maybe on central_orthogonal_idempotents)

[tscrim]

comment:19

Replying to @alauve:

That's right. I was speaking less about "how best" rather than "how".

I don't know how the semisimplicity check gets done...

I don't declare the axiom "Semisimple" at any point when I run

sage: S3 = SymmetricGroupAlgebra(QQ, 3)

versus

sage: S3mod = SymmetricGroupAlgebra(GF(3), 3)

yet somehow Sage knows it.

This is because a finite group algebra can quickly check Maschke's theorem, which is done in the finite_groups.py (look for the Algebras category's extra_super_categories). (Group algebras is too generic as it doesn't special case when it comes from a finite group.)

[tscrim]

comment:20

You can also cache things that take an arbitrary input:

@cached_method(key=lambda self,x: _Partitions(x))
def cpi(self, foo):
    # Need to redo the cast as the first was only for the key, not (local) usage
    foo = _Partitions(foo)

[tscrim]

comment:21

I don't see why we need to deprecate cpi/s and not just keep them as an alias.

Also, +1 to returning a Family indexed by Partitions(n).

[nthiery]

comment:22

I don't see why we need to deprecate cpi/s and not just keep them as an alias.

Oh, just because they are incomprehensible names that pollute (a bit) the namespace without helping the discovery the functionalities. If it was a standard shorthand notation in the community that would be different.

Cheers

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

Changed commit from 7e6c7ef to 6e2d52d

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

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

44f7d8d	Merge branch 'develop' into public/combinat/speedup_idempotents_for_SymmetricGroupAlgebras_25942
6e2d52d	First draft of implementing Nakayama's Conjecture (now a theorem) to correct

[alauve]

comment:24

Code needs a bit more work, but I think we're on the right track now with the changes in this ticket. I still need to read carefully Murphy's proof of Nakayama's Conjecture to be sure I'm defining idempotents correctly (and completely).

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

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

a441a59

allow p-cores as input for central_orthogonal_idempotent

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

Changed commit from 6e2d52d to a441a59

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

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

efa01f7

setup deprecated_function_alias

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

Changed commit from a441a59 to efa01f7

[alauve]

comment:28

I expect this is ready for "needs_review" but I'll way to see if anybody has suggestions before toggling it.

The main result to absorb to understand the changes I've made is

Murphy (J. Alg., 1983), Theorem 2.8, which gives a constructive proof of Nakayama's Conjecture.

---

New commits:

a441a59	allow p-cores as input for central_orthogonal_idempotent
efa01f7	setup deprecated_function_alias

[alauve]

Description changed:

--- 
+++ 
@@ -1 +1,11 @@
-the current implementation of `central_orthogonal_idempotents` is inherited from `FiniteDimensionalSemisimpleAlgebrasWithBasis`. The combinatorial formula in terms of standard tableaux is much faster (implemented as the method `cpis`).
+The current implementation of `central_orthogonal_idempotents` is inherited from `FiniteDimensionalSemisimpleAlgebrasWithBasis`. The combinatorial formula in terms of standard tableaux is much faster (implemented as the method `cpis`).
+
+This ticket does two/four things:
+
+* replaces the default implementation of `central_orthogonal_idempotent/s` with one based on the code in `cpi/s`;
+
+* fixes a bug in `cpi/s`, which assumes base ring contains `QQ`;
+
+* implements primitive central orthogonal idempotents in any characteristic, using Nakayama's Conjecture; 
+
+* adds a `deprecated_function_alias` for the nondescriptive `cpi/s` methods.

[tscrim]

comment:30

Some comments:

• Why are you not displaying all 3 idempotents in the test you changed in categories/finite_dimensional_semisimple_algebras_with_basis.py?
• You don't need to say central_orthogonal_idempotents overrides a generic implementation (if you insist on linking to that, then you should link to the explicit method).
• In central_orthogonal_idempotent, core will all be in math mode. I think you should just use c(\lambda).

• To better follow conventions:

      INPUT:
  
      - ``la`` -- a partition of ``self.n`` or a
        ``self.base_ring().characteristic()``-core of such
        a partition
  
      - ``block`` -- boolean (default: ``True``); when ``False``,
        this returns the classical idempotent associated to ``la``
        (defined over `\QQ`)
  
      OUTPUT:
  
      If ``block=False`` and the corresponding coefficients are
      not defined over ``self.base_ring()``, then return ``None``.
      Otherwise return an element of ``self``.

• You can remove the (taking ``block=False``) part as it is clear from the example.

• I think your error message is better written as

  raise ValueError("the {1}-core of {0} is not a {1}-core of a partition of {2}".format(la, p, self.n))

  I don't see the benefit of showing the user mu.

• The cpi = sum() should not be in the try-except block. Similarly for calling self._blocks_dictionary(). That way you are less likely to catch errors that should be propagating up.
• Better to use _Partitions instead of Partitions() as the former is already an instance.
• In the same vein, you probably can just use _Partitions as the cache key for cpi_over_QQ.
• Avoid calling keys() when you can just iterate over the dict.

• Why are you using gap and not libgap here:

  eval(gap.eval("Display(Irr(SymmetricGroup(%d)));"%n))

  Also, why are you sending it through eval and not using the .sage() of the GAP objects? Actually, why not just use SymmetricGroup(n).character_table()?

• cpi_over_QQ is a function, not a method, so you should use @cached_function.
• The error message in cpi_over_QQ says it should be a partition of 4, which doesn't seem correct.
• You don't need the list of all partitions to get the index. Just get the cardinality C and run {la: C-1-i for i,la in enumerate(Partitions(n))}.

[alauve]

comment:31

Hi Travis,

Most of your suggestions/questions seem reasonable. I'll get to them soon.

Let me focus here on the gap stuff:

• Why are you using gap and not libgap here:

  eval(gap.eval("Display(Irr(SymmetricGroup(%d)));"%n))

  Also, why are you sending it through eval and not using the .sage() of the GAP objects? Actually, why not just use SymmetricGroup(n).character_table()?

I didn't write these lines of code. They are hold-overs from the original .cpi method. I'm guessing you expect a speedup from doing things differently? I'm a novice when it comes to interacting with gap. (Stay tuned, there might be a call for help on the horizon.)

[tscrim]

comment:32

Replying to @alauve:

• Why are you using gap and not libgap here:

  eval(gap.eval("Display(Irr(SymmetricGroup(%d)));"%n))

  Also, why are you sending it through eval and not using the .sage() of the GAP objects? Actually, why not just use SymmetricGroup(n).character_table()?

I didn't write these lines of code. They are hold-overs from the original .cpi method. I'm guessing you expect a speedup from doing things differently? I'm a novice when it comes to interacting with gap. (Stay tuned, there might be a call for help on the horizon.)

Using libgap will definitely be a speedup as it doesn't have to go through a pexpect interface (i.e., convert to/from strings), but it also is cleaner code. Perhaps a better way would be to directly call libgap to avoid more intermediate objects:

sage: [c.sage() for c in libgap.Irr(libgap.SymmetricGroup(4))]
[[1, -1, 1, 1, -1],
 [3, -1, -1, 0, 1],
 [2, 0, 2, -1, 0],
 [3, 1, -1, 0, -1],
 [1, 1, 1, 1, 1]]

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

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

15c3139

addressing reviewer comments

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

Changed commit from efa01f7 to 15c3139

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

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

90a2093

moving from gap to libgap

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

Changed commit from 15c3139 to 90a2093

[alauve]

comment:35

Okay, libgap seems to be a bit faster:

sage: for n in [3,4, 8,9, 12,13]:
....:     print n
....:     print "gap.eval:"
....:     timeit('character_table = eval(gap.eval("Display(Irr(SymmetricGroup(%d)));"%n))')
....:     print "libgap:"
....:     timeit('_character_table = [c.sage() for c in libgap.Irr(libgap.SymmetricGroup(n))]')
....:     print
....:     
3
gap.eval:
5 loops, best of 3: 6.84 ms per loop
libgap:
125 loops, best of 3: 5.52 ms per loop

4
gap.eval:
25 loops, best of 3: 7.08 ms per loop
libgap:
25 loops, best of 3: 8.29 ms per loop

8
gap.eval:
5 loops, best of 3: 48.3 ms per loop
libgap:
5 loops, best of 3: 49.5 ms per loop

9
gap.eval:
5 loops, best of 3: 97.8 ms per loop
libgap:
5 loops, best of 3: 72.6 ms per loop

12
gap.eval:
5 loops, best of 3: 971 ms per loop
libgap:
5 loops, best of 3: 273 ms per loop

13
gap.eval:
5 loops, best of 3: 3.72 s per loop
libgap:
5 loops, best of 3: 573 ms per loop

---

New commits:

15c3139	addressing reviewer comments
90a2093	moving from gap to libgap

---

New commits:

15c3139	addressing reviewer comments
90a2093	moving from gap to libgap

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

Changed commit from 90a2093 to 5b752c4

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

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

d77425d	Merge branch 'public/combinat/speedup_idempotents_for_SymmetricGroupAlgebras_25942' of git://trac.sagemath.org/sage into public/combinat/speedup_idempotents_for_SymmetricGroupAlgebras_25942
86bfd32	Changing the caching of SGA central_orthogonal_idempotent(s).
5b752c4	Some tweaks to the changed test in fin-dim semiseimple algebras w/ basis.

[tscrim]

comment:38

The big change I made was that I altered how the caching was done for the idempotents. Specifically, I moved the caching into the SGA itself with a custom cache. Although I also write one that uses @cached_method with a slightly intricate key to extract the p-core when needed.

Advantages:

• Things are not permanently nailed in a (non-publically accessible) function cache. So when you stop working with a particular SGA, then it will be garbage collected.
• You do not have to redo the computation across partitions that have the same p-core.
• The character table and indices does not have to constantly be recomputed across a block.
• We can work with lower-level objects (e.g., a dict) safely.

Disadvantages:

• The cache of cpi's over Q is no longer shared. So if you create a bunch of SGAs with different fields (say, with different characteristics) but with the same n, the Q computation will have to be redone for each of the SGAs.

I think the advantages, in particular the better locality of the cache, outweighs the disadvantages. Now we can get around the looks-like-a-memory-leak issue by making the separate function cache a @weak_cache_function, but then it is possible (maybe even likely) the results will be GC'ed and then have to be redone.

Some timings:

sage: S3 = SymmetricGroup(10).algebra(GF(3))
sage: _ = S3._blocks_dictionary  # So this is created and will not affect the timings
sage: %time idem = S3.central_orthogonal_idempotent([1])
CPU times: user 1.2 s, sys: 12 ms, total: 1.21 s
Wall time: 1.21 s
sage: %time idem = S3.central_orthogonal_idempotent([2,1,1])
CPU times: user 332 ms, sys: 44 ms, total: 376 ms
Wall time: 333 ms
sage: %time idem = S3.central_orthogonal_idempotent([3,1])
CPU times: user 416 ms, sys: 36 ms, total: 452 ms
Wall time: 405 ms

vs caching _blocks_dictionary and other optimizations:

sage: %time idem = S3.central_orthogonal_idempotent([1])
CPU times: user 2.51 s, sys: 248 ms, total: 2.76 s
Wall time: 2.77 s
sage: %time idem = S3.central_orthogonal_idempotent([2,1,1])
CPU times: user 1.26 s, sys: 236 ms, total: 1.49 s
Wall time: 1.49 s
sage: %time idem = S3.central_orthogonal_idempotent([3,1])
CPU times: user 1.2 s, sys: 192 ms, total: 1.4 s
Wall time: 1.39 s

vs your branch

sage: %time idem = S3.central_orthogonal_idempotent([1])
CPU times: user 2.7 s, sys: 260 ms, total: 2.96 s
Wall time: 3 s
sage: %time idem = S3.central_orthogonal_idempotent([2,1,1])
CPU times: user 1.84 s, sys: 312 ms, total: 2.15 s
Wall time: 2.15 s
sage: %time idem = S3.central_orthogonal_idempotent([3,1])
CPU times: user 1.92 s, sys: 312 ms, total: 2.23 s
Wall time: 2.23 s

---

I also did a bunch of other small fixes and tweaks (in particular, I constructed cubed roots by using a cyclotomic field). If my changes are good, then positive review.

[tscrim]

Changed author from Mee Seong Im to Mee Seong Im, Aaron Lauve

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

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

3223375

Further optimizations of idempotent computation in positive char.

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

Changed commit from 5b752c4 to 3223375

[tscrim]

comment:40

I ended up seeing another way to get further speed out in positive characteristic: do not compute the cycle types over and over again. This leads to the following timings:

sage: %time idem = S3.central_orthogonal_idempotent([1])
CPU times: user 360 ms, sys: 36 ms, total: 396 ms
Wall time: 393 ms
sage: %time idem = S3.central_orthogonal_idempotent([2,1,1])
CPU times: user 120 ms, sys: 8 ms, total: 128 ms
Wall time: 120 ms
sage: %time idem = S3.central_orthogonal_idempotent([3,1])
CPU times: user 192 ms, sys: 24 ms, total: 216 ms
Wall time: 195 ms

Perhaps the best thing to do would be to iterate over the partitions in reverse order and then iterate over the conjugacy classes. However, this is sufficient for now I think as it is already a 10x speedup.