optimize modular symbols decomposition algorithm

[williamstein]

In short, the decomposition function on spaces of modular symbols is mysteriously way slower than it should be. Why?

Consider this:

sage: M = ModularSymbols(1000,2,sign=1).new_subspace().cuspidal_subspace()
sage: time d = M.decomposition(3)
CPU times: user 3.21 s, sys: 0.11 s, total: 3.33 s
Wall time: 3.37 s
sage: t3 = M.hecke_matrix(3)
sage: time d = t3.decomposition()
CPU times: user 0.11 s, sys: 0.00 s, total: 0.11 s
Wall time: 0.11 s
sage: time d = t3.decomposition(algorithm='multimodular', height_guess=1)
CPU times: user 0.06 s, sys: 0.00 s, total: 0.06 s
Wall time: 0.06 s

This huge timing discrepancy isn't due to caching:

^bsd:matrix was$ sage
----------------------------------------------------------------------
| Sage Version 3.2, Release Date: 2008-11-20                         |
| Type notebook() for the GUI, and license() for information.        |
----------------------------------------------------------------------

sage: M = ModularSymbols(1000,2,sign=1).new_subspace().cuspidal_subspace()
sage: t3 = M.hecke_matrix(3)
sage: time d = t3.decomposition(algorithm='multimodular', height_guess=1)
CPU times: user 0.07 s, sys: 0.01 s, total: 0.08 s
Wall time: 0.08 s

For comparison:

sage: magma.eval("M := ModularSymbols(1000,2,1);")
''
sage: magma.eval("S := NewSubspace(CuspidalSubspace(M)); time D := Decomposition(S, 3);")
'Time: 0.050'

So Sage is nearly the same as Magma at the decomposition part of the computation, but is getting totally killed by using the wrong algorithm or doing something really dumb that it shouldn't even bother doing. I.e., above 3.2 seconds is spent doing something probably unnecessary, and only 0.08 is spent doing what should be the dominant step.

There are of course numerous other similar examples. For concreteness, I think to close this ticket one should just worry about making it so that the above example completes in < 0.2 seconds instead of 3.3 seconds.

Depends on:

1. 10987

Apply:

1. trac-4578-optimize_modular_symbols_decomposition-v2.patch
2. trac-4578-optimize_modular_symbols_decomposition-doctest.patch

CC: @sagetrac-GeorgSWeber @aghitza

Component: modular forms

Author: Martin Raum

Reviewer: David Loeffler

Merged: sage-4.7.2.alpha2

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

[18d65770-dc1e-4813-89a9-4828e4aae4a9]

comment:3

There are two reasons why this is much slower than we expect it to be: First, the restriction of subspaces, when decomposing them, is checked. This is not necessary and already this increases the performance. Second, sorting the resulting Hecke modules depends on computing several Hecke matrices. This we cannot change without introducing major incompatibilities. But I added an option, that sorts the modules by their basis structure. This makes the function again much fast.

sage: M = ModularSymbols(1000,2,sign=1).new_subspace().cuspidal_subspace()
sage: %time d  = M.decomposition(3)
CPU times: user 1.62 s, sys: 0.00 s, total: 1.63 s
Wall time: 1.92 s

and

sage: M = ModularSymbols(1000,2,sign=1).new_subspace().cuspidal_subspace()
sage: %time d  = M.decomposition(3, sort_by_basis = True)
CPU times: user 0.94 s, sys: 0.00 s, total: 0.94 s
Wall time: 1.59 s

Note, that the bare decomposition given in the description is not sufficient: We need at least compute two further decompositions, since two of the resulting modules are not irreducible with respect to T(2).

[18d65770-dc1e-4813-89a9-4828e4aae4a9]

Author: Martin Raum

[18d65770-dc1e-4813-89a9-4828e4aae4a9]

Description changed:

--- 
+++ 
@@ -43,4 +43,7 @@

 So Sage is nearly the same as Magma at the decomposition part of the computation, but is getting totally killed by using the wrong algorithm or doing something really dumb that it shouldn't even bother doing.  I.e., above 3.2 seconds is spent doing something probably unnecessary, and only 0.08 is spent doing what should be the dominant step. 

-There are of course numerous other similar examples.   For concreteness, I think to close this ticket one should just worry about making it so that the above example completes in < 0.2 seconds instead of 3.3 seconds. 
+There are of course numerous other similar examples.   For concreteness, I think to close this ticket one should just worry about making it so that the above example completes in < 0.2 seconds instead of 3.3 seconds.
+
+**Apply:**
+1. trac-4578-optimize_modular_symbols_decomposition.patch

[18d65770-dc1e-4813-89a9-4828e4aae4a9]

Description changed:

--- 
+++ 
@@ -45,5 +45,8 @@

 There are of course numerous other similar examples.   For concreteness, I think to close this ticket one should just worry about making it so that the above example completes in < 0.2 seconds instead of 3.3 seconds.

+**Depends on:**
+1. #10987
+
 **Apply:**
 1. trac-4578-optimize_modular_symbols_decomposition.patch

[11d1fc49-71a1-44e1-869f-76be013245a0]

comment:5

Hang on. Surely this patch doesn't remotely match the ticket description? If called with the default arguments this patch is a complete no-op!

Did you perhaps want to call decomposition_of_subspace with the argument "check_restrict=False"? Because Sage's Hecke algebras are always commutative, this check is probably redundant. That would, I imagine, result in a substantial speedup.

[18d65770-dc1e-4813-89a9-4828e4aae4a9]

Description changed:

--- 
+++ 
@@ -50,3 +50,4 @@

 **Apply:**
 1. trac-4578-optimize_modular_symbols_decomposition.patch
+2. trac-4578-optimize_modular_symbols_decomposition-doctest.patch

[18d65770-dc1e-4813-89a9-4828e4aae4a9]

comment:6

Attachment: trac-4578-optimize_modular_symbols_decomposition-doctest.patch.gz

You're completely right. Thank you for catching this! I replaced the old patch with the right one. Doing this it turned out, that I needed to change one doctest. The change is justified as the element alpha - 1 and 1/2 alpha + 1/2 have the same minpolys. They are, thought, elements of number fields defined with completely different polynomial (They are isomorphic!). This happens because the number field is now generated by a different piece of code in the linalg modules of Sage.

The old and new parent are for me

OLD:
Number Field in alpha with defining polynomial x^2 + 4*x + 1

NEW:
Number Field in alpha with defining polynomial x^2 - 2*x - 11

[11d1fc49-71a1-44e1-869f-76be013245a0]

comment:7

I'm happy with that change. As Martin has pointed out, the timing comparison with Magma is misleading: Sage is computing much more information than Magma (decomposing a bunch more Hecke operators). So I advocate merging Martin's patch and closing this ticket.

[11d1fc49-71a1-44e1-869f-76be013245a0]

Reviewer: David Loeffler

[jdemeyer]

comment:9

There is an issue with the documentation:

dochtml.log:/mnt/usb1/scratch/jdemeyer/merger/sage-4.7.2.alpha1/local/lib/python2.6/site-packages/sage/modular/hecke/module.py:docstring of sage.modular.hecke.module.HeckeModule_free_module.T:7: (WARNING/2) Bullet list ends without a blank line; unexpected unindent.

[18d65770-dc1e-4813-89a9-4828e4aae4a9]

Description changed:

--- 
+++ 
@@ -49,5 +49,5 @@
 1. #10987

 **Apply:**
-1. trac-4578-optimize_modular_symbols_decomposition.patch
+1. trac-4578-optimize_modular_symbols_decomposition-v2.patch
 2. trac-4578-optimize_modular_symbols_decomposition-doctest.patch

[18d65770-dc1e-4813-89a9-4828e4aae4a9]

comment:10

Attachment: trac-4578-optimize_modular_symbols_decomposition-v2.patch.gz

This was a very stupid typo. I have checked that it occurs with the old patch and does not occur with the new one. The only thing that I changes is two spaces, one added, one removed. So I set this back to positive review; Sorry for the inconvenience, Jeroen!

[jdemeyer]

Merged: sage-4.7.2.alpha2