Hecke series for overconvergent modular forms

[sagetrac-lauder]

Computes to any required p-adic precision the characteristic series of the
Atkin U_p operator on the space of p-adic overconvergent modular forms of
given weight k and level N.

Apply:

1. attachment: trac_12043-all.patch
2. attachment: 12043-fix.patch
3. attachment: 12043_long_time.patch

Depends on #12724
Depends on #12740

CC: @loefflerd

Component: modular forms

Author: Alan Lauder

Reviewer: David Loeffler, Jan Vonk

Merged: sage-5.3.beta0

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

[sagetrac-janv]

comment:2

Reviewed this patch, everything works perfectly and looks nice.

[sagetrac-janv]

comment:3

Ignore the previous comment, I have not reviewed this. Sorry!

[loefflerd]

comment:5

(For the record: I am in the process of reviewing this.)

[loefflerd]

comment:6

Hi Alan,

This is really nice code: I'm very happy that you chose to make this Sage implementation available.

I am preparing a second reviewer patch, to be applied on top of your patch, which mainly cleans up the documentation and removes some functions which already existed in Sage under different names (such as a lot of the routine linear algebra stuff). However, I've not completed my review, because I think there's a serious lurking bug.

The problem is this. The code in sage/modular/modform/find_generators.py is only guaranteed to return generators for modular forms as a Q-algebra. You've carefully trapped and dealt with the case when the generators aren't p-adically integral; but even if the generators happen to be integral, there's no guarantee that they generate the ring of p-adically integral modular forms as a Zp-algebra. I experimented by making ModularFormsRing.generators() multiply its output by 25 before returning it, which is completely consistent with that function's specifications; and then doing 5-adic computations with modformsring=True just looped forever. So I don't think we can use this code without also adjusting ModularFormsRing.generators() to add the option to do its linear algebra integrally.

There are a few other (relatively minor) things that puzzle me about this implementation, as well, and I'd be grateful if you could explain them to me.

• Is the algorithm, and the implementation, still OK if k is a negative integer? It seems to be, based on all of the tests I ran -- the output for a negative k looks entirely consistent with the output for positive integers p-adically close to k. If so, then what does it mean to say that the run time is polynomial in the logarithm of the weight? Is it polynomial in the absolute value of the weight?

• Should the level N be prime to p? Jan Vonk made this observation at the Sage Days a few weeks back. If you run the code with N divisible by p, it doesn't complain, and returns the same output as for the prime-to-p part of N, but much more slowly! So we should trap this case. I added code to do this to my reviewer patch.

• For N > 1, a great deal of time is wasted in the routine "row_reduced_form". For instance, I ran "hecke_series(5, 3, 4, 15, modformsring=True)" with the Sage profiler. The command took 29 seconds, of which 24 seconds was spent in 213 calls to the method "row_reduced_form". But only 0.052 seconds of that was spent in the matrix echelon form routine (Sage's mod p linear algebra is pretty quick). The rest is all spent, as far as I can see, in translating back and forth between power series and vectors.

• In the modformsring = False version of the algorithm, it seems to be randomized twice over: firstly, the low weight bases are randomized (using random_low_weight_bases); then a randomized algorithm using these as input (random_new_basis_modp) is used to compute the complementary spaces. What is the benefit of randomizing twice in this way?

[loefflerd]

comment:7

Here's my reviewer patch. Comments and questions in the code are marked with XXX. I emphasize that it's incomplete, because I stopped once I realized the problem with Z versus Q generators. Still it does a large part of the work of hooking your code into Sage (adding the new module to the global namespace, to the reference manual, etc) and thus I'd be grateful if you worked on top of this patch for the next round of corrections.

[sagetrac-lauder]

comment:8

Dear David

Many thanks for your message, and for looking so carefully at the code. This is greatly appreciated!

I am a novice at the SAGE review process, so will wait to get advice from Sebastian (Pancratz) on exactly what to do next with
the reviewer's patch. In the meantime though perhaps I could address your comments.

On the Z versus Q basis, I was aware of that problem. In fact, when modformsring=false the algorithm will not terminate
either if the space of modular forms in low weight does not "generate" a certain Zp-space of modular forms. (In the magma implementation preamble I mention that the algorithm is not guaranteed to terminate: I should have made this more explicit in the sage one.) However, when modformsring=false in practice if the algorithm does fail to terminate (which I have never
encountered with the default setting of weightbound = 6) then one can increase "weightbound". (Also, the output is in any case I believe provably correct, so it never terminates with a false output.)

When modformsring=true, as you noticed the algorithm defaults to modformsring = false when the ModularFormsRings.generators() gives non p-adically integral generators. I did ask William (Stein) about
the possibility of having an integral_generators() function in this case: looking at the source code I could not see an
obvious obstruction to this. The key point for me though was the following: The generators() function only gives a set
of modular forms which generate over Q all forms up to maxweight (default 20). If one increased maxweight to a value which
guaranteed that (over Q at least) one had all the modular forms you need (weight k0 + n(p-1) in the notation in my paper)
then the generators() function would take much much longer. So even if one had an integral_generators() function, to
absolutely guarantee that the algorithm terminates one would need "maxweight = k0 + n(p-1)" and this would seriously
impact on the time taken to find the generators. So I thought it is better to live with the possibility of the algorithm not terminating, since it seems unobservably small (at least, I have not observed it when modformsring = false and weightbound = 6, although I have not experimented so much when modformsring = true).

That said though I think your point is a very good one. In particular, if one had an integral_generators() function then
one could just set modformsring = true "permanently" (and perhaps replace the "weightbound" optional parameter
with "maxweight" in the main function). In any case, the fact that the algorithm may in principle fail to terminate needs to be clear in the documentation.

On your other points:

• Yes, negative k is fine and the algorithm has running time polynomial log(|k| + 1) . (I do clarify this point at the start of
  Section 3.2.2 in my paper, but have got into the habit of just writing polynomial in log(k).)

• Yes, for the theoretical analysis of the algorithm I need N prime to p, although as you observed it does seem to work in practice (and perhaps in theory) regardless. Again, the magma implementation picks this up: thanks for sorting it out here.

• This is a very interesting observation: I had not realised so much time was taken up doing this! I would have to think about this more carefully, but I think the main point of the call to row_reduced_form is simply to detect when one has increased the rank of the space of mod p modular forms. Since "TotalBasisModp" was already in "echelon form" (as a list of q-series)
  before "TotalBasisi" was appended, the fact that the matrix echelon form is so quick is not surprising: only the final row needs to be done. It sounds to me that it might be much better to avoid matrices altogether, and just write a little function which
  given a list of q-series in "echelon form" modulo p, appends a new q-series and "echelonises". Many thanks for spotting this!
  (I must learn how to use the sage profiler.)

• The problem I encountered when not randomising the low weight bases was as follows: the low weight bases are given by default in echelon form, at least a union of the cuspidal and eisenstein spaces in that form. When you multiply random
  elements from such a basis together, they have a tendency to vanish (to higher order) at q = 0. I noticed in an early
  implementation that I could generate the spaces more quickly by randomising the bases to start with. So I have always done this: it is probably worth checking again though whether it is really helping.

Best wishes

Alan.

Replying to @loefflerd:

Hi Alan,

This is really nice code: I'm very happy that you chose to make this Sage implementation available.

I am preparing a second reviewer patch, to be applied on top of your patch, which mainly cleans up the documentation and removes some functions which already existed in Sage under different names (such as a lot of the routine linear algebra stuff). However, I've not completed my review, because I think there's a serious lurking bug.

The problem is this. The code in sage/modular/modform/find_generators.py is only guaranteed to return generators for modular forms as a Q-algebra. You've carefully trapped and dealt with the case when the generators aren't p-adically integral; but even if the generators happen to be integral, there's no guarantee that they generate the ring of p-adically integral modular forms as a Zp-algebra. I experimented by making ModularFormsRing.generators() multiply its output by 25 before returning it, which is completely consistent with that function's specifications; and then doing 5-adic computations with modformsring=True just looped forever. So I don't think we can use this code without also adjusting ModularFormsRing.generators() to add the option to do its linear algebra integrally.

There are a few other (relatively minor) things that puzzle me about this implementation, as well, and I'd be grateful if you could explain them to me.

• Is the algorithm, and the implementation, still OK if k is a negative integer? It seems to be, based on all of the tests I ran -- the output for a negative k looks entirely consistent with the output for positive integers p-adically close to k. If so, then what does it mean to say that the run time is polynomial in the logarithm of the weight? Is it polynomial in the absolute value of the weight?

• Should the level N be prime to p? Jan Vonk made this observation at the Sage Days a few weeks back. If you run the code with N divisible by p, it doesn't complain, and returns the same output as for the prime-to-p part of N, but much more slowly! So we should trap this case. I added code to do this to my reviewer patch.

• For N > 1, a great deal of time is wasted in the routine "row_reduced_form". For instance, I ran "hecke_series(5, 3, 4, 15, modformsring=True)" with the Sage profiler. The command took 29 seconds, of which 24 seconds was spent in 213 calls to the method "row_reduced_form". But only 0.052 seconds of that was spent in the matrix echelon form routine (Sage's mod p linear algebra is pretty quick). The rest is all spent, as far as I can see, in translating back and forth between power series and vectors.

• In the modformsring = False version of the algorithm, it seems to be randomized twice over: firstly, the low weight bases are randomized (using random_low_weight_bases); then a randomized algorithm using these as input (random_new_basis_modp) is used to compute the complementary spaces. What is the benefit of randomizing twice in this way?

[loefflerd]

comment:9

Replying to @sagetrac-lauder:

When modformsring=true, as you noticed the algorithm defaults to modformsring = false when the ModularFormsRings.generators() gives non p-adically integral generators. I did ask William (Stein) about
the possibility of having an integral_generators() function in this case: looking at the source code I could not see an
obvious obstruction to this.

Are you asking for integral_generators to return (a) modular forms which have integral q-expansions and are Q-algebra generators for the ring of all modular forms, or (b) Z-algebra generators for the ring of integral modular forms? The former would be completely trivial to do. The latter is certainly possible, but would be harder; in fact it would be a great thing to have anyway, but it would be much slower than the current implementation of course.

The key point for me though was the following: The generators() function only gives a set
of modular forms which generate over Q all forms up to maxweight (default 20). If one increased maxweight to a value which
guaranteed that (over Q at least) one had all the modular forms you need (weight k0 + n(p-1) in the notation in my paper)
then the generators() function would take much much longer. So even if one had an integral_generators() function, to
absolutely guarantee that the algorithm terminates one would need "maxweight = k0 + n(p-1)" and this would seriously
impact on the time taken to find the generators. So I thought it is better to live with the possibility of the algorithm not terminating, since it seems unobservably small (at least, I have not observed it when modformsring = false and weightbound = 6, although I have not experimented so much when modformsring = true).

I'm pretty sure it's true that for all N, the ring of modular forms of level N is generated over Q by forms of weight at most 6. Working out the details is a fiddly exercise which I should probably just sit down and do one of these days. But of course to get integral generators you might need larger weights.

Perhaps using some more algebraic geometry one can bound above the max weight needed for generators over Zp (at least for p not dividing the level). There is some yoga with ample line bundles that I've never quite got my head around.

That said though I think your point is a very good one. In particular, if one had an integral_generators() function then
one could just set modformsring = true "permanently" (and perhaps replace the "weightbound" optional parameter
with "maxweight" in the main function). In any case, the fact that the algorithm may in principle fail to terminate needs to be clear in the documentation.

I'm not sure what the policy is on including algorithms which aren't known to terminate, I'd have to check. That would simplify the code a lot, which would be nice. We can certainly make the generators code return generators over Q which were integral (but not necessarily integral generators), which would mean we could remove the fallback code.

On your other points:

• Yes, negative k is fine and the algorithm has running time polynomial log(|k| + 1) . (I do clarify this point at the start of
  Section 3.2.2 in my paper, but have got into the habit of just writing polynomial in log(k).)

OK, thanks; I missed that in the paper.

• Yes, for the theoretical analysis of the algorithm I need N prime to p, although as you observed it does seem to work in practice (and perhaps in theory) regardless. Again, the magma implementation picks this up: thanks for sorting it out here.

Well, thanks here are due to Jan, not me :-)

• This is a very interesting observation: I had not realised so much time was taken up doing this! I would have to think about this more carefully, but I think the main point of the call to row_reduced_form is simply to detect when one has increased the rank of the space of mod p modular forms. Since "TotalBasisModp" was already in "echelon form" (as a list of q-series)
  before "TotalBasisi" was appended, the fact that the matrix echelon form is so quick is not surprising: only the final row needs to be done. It sounds to me that it might be much better to avoid matrices altogether, and just write a little function which
  given a list of q-series in "echelon form" modulo p, appends a new q-series and "echelonises". Many thanks for spotting this!

There is a choice to be made here: one can either work consistently with power series objects, or consistently with vectors. You seem to be suggesting the first, but I'd strongly recommend the second: this allows you to use Sage's very fast echelon form routines, which are written in C, rather than using hand-rolled Python versions which would inevitably be much slower. There is some code in modforms/find_generators.py which does something quite similar to this, and that works with vectors where possible.

(I must learn how to use the sage profiler.)

It's really simple: all I did was to type

sage: %prun hecke_series(5, 3, 4, 15, modformsring=True)

and that does the job.

• The problem I encountered when not randomising the low weight bases was as follows: the low weight bases are given by default in echelon form, at least a union of the cuspidal and eisenstein spaces in that form. When you multiply random
  elements from such a basis together, they have a tendency to vanish (to higher order) at q = 0. I noticed in an early
  implementation that I could generate the spaces more quickly by randomising the bases to start with. So I have always done this: it is probably worth checking again though whether it is really helping.

I see. Yes, that makes a lot of sense -- if your random products of the basis in random_new_basis_modp are landing in a specific subspace with high probability then it's clearly not good.

You'll notice that in the reviewer patch I got rid of a call to GL(n, p).random_element(), which was overwhelmingly dominating the run time in the modformsring=False case; on my machine GL(4, 3).random_element() takes nearly two whole seconds! This is a bug, and should probably be reported as such.

[sagetrac-lauder]

comment:10

Dear David

Thanks for your message.

On integral_generators, ideally I would like (b).

Regarding the algorithm terminating, I thought very briefly about bounding the number of random q-series
generated during the attempt to construct the complementary spaces modulo p, with the algorithm then extending the
low weight bases if this was exceeded. This would of course guarantee eventual termination! The funny thing is though that if you get the code to print the rank "rk" when generating the bases modulo p, you will see it looks rather uneven, with the routine getting "stuck" in certain subspaces and then whizzing off again. This didn't fit with my heuristic idea of how the rank
should be increasing. So in short, since I couldn't think of a mathematically elegant bound on how many series to generate, and coding it up this way looked a bit fiddly, I didn't bother. Probably it is better to wait and see what we do with integral_generators.

I followed your suggestion, and cut row_reduced_form altogether, replacing the list of q-series modulo p with a matrix. This
is much better, but in larger examples doesn't give such a dramatic speed-up. (The same idea applied to the magma implementation doubled the speed though!) I have not applied this to your patch yet though, just my old sage copy.

Best wishes

Alan.

Replying to @loefflerd:

[sagetrac-lauder]

comment:11

Dear David

I have applied your patch (with Sebastian's help) and then made my changes. I created a new patch myself: hopefully I did the right thing. In any case, I will email you my new source code directly (just in case I messed up).

Note that your new version failed on a few tests on my machine: there was a stray "print" statement you must have put in sometime; also, it complained about some alteration you made to the little function for testing if a weight list was valid. I did not understand what was happening with the latter, so I have not touched this: when I tested my new version it passed everything OK except that weight list one.

I tried to address all of your smaller concerns. The only substantial change was that I altered the random_new_basis_modp function so that it does not waste time changing between matrices and lists: this also involved a small change to the complementary_spaces_modp function.

I meant to check why the new version was slower for large p, but forgot. Perhaps you could have a look at whether using eisenstein series over Q is slowing down the computation of E_(p-1) for, say, p = 79, N = 1, k = 2 and m = 3.

Anyway, many thanks again for all the improvements you made to the code!

Best wishes

Alan.

[loefflerd]

Reviewer: David Loeffler

[loefflerd]

Changed author from lauder to Alan Lauder

[loefflerd]

comment:12

Replying to @sagetrac-lauder:

I have applied your patch (with Sebastian's help) and then made my changes. I created a new patch myself: hopefully I did the right thing. In any case, I will email you my new source code directly (just in case I messed up).

No problems there, the patch applied fine.

Note that your new version failed on a few tests on my machine: there was a stray "print" statement you must have put in sometime; also, it complained about some alteration you made to the little function for testing if a weight list was valid. I did not understand what was happening with the latter, so I have not touched this: when I tested my new version it passed everything OK except that weight list one.

I see what's happened there -- silly mistake on my part, I'll fix it.

I tried to address all of your smaller concerns. The only substantial change was that I altered the random_new_basis_modp function so that it does not waste time changing between matrices and lists: this also involved a small change to the complementary_spaces_modp function.

I meant to check why the new version was slower for large p, but forgot. Perhaps you could have a look at whether using eisenstein series over Q is slowing down the computation of E_(p-1) for, say, p = 79, N = 1, k = 2 and m = 3.

The new version is actually about 50 times faster. For instance, computing 10000 terms of E_{78} mod 79 takes 10.23 seconds with your old code, and 0.49 seconds with the new code.

Anyway, many thanks again for all the improvements you made to the code!

Best wishes

Alan.

No problem. I am a bit busy at the moment, but when I get a chance I will look at your changes, sort out the last few minor issues and then I think it'll be good to go.

Best regards, David

[sagetrac-lauder]

comment:13

Dear David

Many thanks for your message. I'm glad it worked!

On the level 1 code, on my machine setting p = 79, N = 1, k = 2, m = 3, the original version took 40s and the new version 43.5s. The algorithm is not randomised here, so something genuinely appears to be taking longer. Anyway, perhaps when you come to try out the new code you could have a quick look at this.

Best wishes

Alan.

[loefflerd]

comment:14

Whatever is causing the slowdown, it's not my changes to the Eisenstein series code. With the parameters you suggest, your original code takes 0.987 seconds per call to "eisen_series", my first round of changes in attachment: trac_12043-reviewer.DRAFT.patch gets this down to 0.080, and with some minor tinkering elsewhere in the Sage library I can speed this up another whole order of magnitude, to 0.009 seconds.

However, there is indeed a small and persistent slowdown, which cancels out the speedup in eisen_series. Looking at the profiler output, this slowdown is in the function level1_UpGj (rather than any of the functions it calls); it is this routine that overwhelmingly dominates the calculation in this example, with 95% of the running time spent here. I'm puzzled because I haven't actually changed anything nontrivial in this function! I'll keep looking.

[loefflerd]

comment:15

OK, so I got to the bottom of this.

Something like 90% of the running time for 'hecke_series(p = 79, N = 1, klist = 2, m = 3)' is taken up by the following 3 lines of code, the first part of Step 4 of the algorithm:

    Ep1p = Ep1(q**p)
    Ep1pinv = Ep1p**(-1)
    G = Ep1*Ep1pinv

These are slowed down slightly by my changes, since my first reviewer patch added an O(q^x) error placeholder onto the end of Ep1, which slightly increases the time spent on this computation. But one can vastly speed up these lines by:

• truncating Ep1, reducing its precision by a factor of p, before substituting q^p in

• using the method "V" (Verschiebung?) of Sage power series objects, which does the map f(q) |--> f(q^p) much more quickly than the generic power series composition code.

With the code above changed to

    Ep1p = (Ep1.add_bigoh(ceil(Ep1.prec() / ZZ(p)))).V(p)
    Ep1pinv = Ep1p**(-1)
    G = Ep1*Ep1pinv

the execution time drops from 43.24 s to 2.68 s.

It's a bit of a mystery to me why the first version is quite so slow; I guess it's calling some general-purpose routine for composition of power series, not realizing that what we're substituting is so simple.

How does the new timing compare with Magma?

[sagetrac-lauder]

comment:16

Dear David

Many thanks for looking into this!

By the way: I am in Japan for the next five weeks, with very slow internet access and a keyboard that keeps switching to Japanese (then I have to start all over again). So it has taken about an hour to get this far ...

That timing is about the same as magma now. It is great you spotted that! I did try a larger example: p = 59, N = 1, k = 2 and m = 5. Here magma was about 7 seconds and sage 28. So I think there is probably still a big difference, even in level 1. (Magma code is on my website, in case you have access to magma and want to try it.)

Incidentally: you are of course quite right that the end of the Miller basis gives the complementary spaces, but the start has nothing particularly to do with the image ... Please correct the blurb appropriately. (See my silly comment in the patch.)

I am going to finish now, before the keyboard goes funny and I have to start again.

Best wishes

Alan.

[loefflerd]

comment:17

Here's what I get with all my optimizations:

sage: time hecke_series(p = 59, N = 1, klist = 2, m = 5)
690689577*x^10 + 227559932*x^9 + 426975979*x^8 + 473462905*x^7 + 665599061*x^6 + 356193732*x^5 + 408559793*x^4 + 585093835*x^3 + 432796757*x^2 + 22614222*x + 1
Time: CPU 1.34 s, Wall: 1.52 s

So the Sage implementation has pulled ahead of the Magma one now, at least for these level 1 cases.

I think there is only one more thing we need to do, if we're willing to live with the unproven conjecture about weights of generators. That is, to modify the ModularFormsRing code so the generators it returns are always integral, and so that the ZZ-submodule they generate is saturated in weights up to the given bound. This is not too hard to do -- I already have code that does it, I just need to merge it into the Sage library. Then I think we can get this in (modulo yourself or someone else reviewing my contributions). I'm not saying that there isn't room for further optimizations, but those can come in subsequent tickets.

[loefflerd]

All patches folded into one

[loefflerd]

Dependencies: #12724, #12740

[loefflerd]

comment:18

Attachment: trac_12043-all.patch.gz

I've uploaded a new revised version, in two forms: attachment: trac_12043-part4.patch is the difference from the previous version, and attachment: trac_12043-all.patch is a single patch combining the effects of the four patches so far.

This version is dependent on two new tickets #12724 and #12740. The first of these (already positively reviewed, thanks Alex!) adds options for different normalizations to the top-level Eisenstein series function, incorporating the functionality from the eisen_series function in Alan's original patch. The second is a complete overhaul of the ModularFormsRing code, which means (among other things) that it can now be required to return generators for the mod p modular forms; this removes a possible cause of the algorithm failing to terminate (if one of the generators returned by ModularFormsRing happened to be divisible by p, which could easily have happened with the old code). I've also applied a number of optimizations, particularly for level 1 modular forms, as described in my last few comments.

I think this is ready to go now. All that's required is for someone to have a look at my part4 patch and verify that the changes look reasonable.

[loefflerd]

comment:19

Apply trac_12043-all.patch

(for the patchbot)

[loefflerd]

minor fix -- apply over previous folded patch

[loefflerd]

comment:20

Attachment: 12043-fix.patch.gz

Apply trac_12043-all.patch, 12043-fix.patch

(for patchbot). This patch just changes the output type when a list of spaces is called for, all of which are zero (returning a list of ones).

[sagetrac-janv]

comment:21

Thanks for fixing this David, everything works absolutely impeccably.
I ran all the tests, and believe that this is now ready to go.

[jdemeyer]

Merged: sage-5.3.beta0

[jdemeyer]

Additional patch

[jdemeyer]

comment:23

Attachment: 12043_long_time.patch.gz

Additional patch attachment: 12043_long_time.patch needs review. This is a blocker patch.

[jdemeyer]

Changed reviewer from David Loeffler to David Loeffler, Jan Vonk

[loefflerd]

comment:25

Looks harmless enough to me. I'd be happy to reinstate the positive review, but I don't have sufficient trac privileges to reopen the ticket in order to do so.

[jdemeyer]

comment:26

No need to reopen the ticket, I'll just merge the additional patch.