Precision issues in random families

[rharron]

Now, there's what appears to be a precision problem in the random families function. I wrote a script called RankOne.sage that computes q-expansions of Hida families whose space of ordinary modular symbols of sign -1 is one-dimensional (so I take primes (> 2) that exactly divide the level). It makes it through levels 11 and 14, but crashes on level 15 with p = 3 and p = 5. Namely, after adjusting t, the 0th moment of t is not quite zero (i.e. its valuation is quite high, but not high enough). I was able to fix that problem for p = 5 by first artificially increasing M_in by 2 and then removing that extra precision right before solving the difference equation. But then, the symbol it creates fails the consistency check. For p = 3, increasing M_in by 2 was insufficient, as was 6, though 10 worked. But again the symbol eventually created failed the consistency check.

Here's some simple code to reproduce what's happening:

sage: p = 5 sage: DD = FamiliesOfOverconvergentDistributions(0, base_coeffs=ZpCA(p, 20), prec_cap=[8,8]) sage: MM = FamiliesOfOMS(15, 0, coefficients=DD) sage: Phis = MM.random_element() sage: Phis._consistency_check()

I'll upload some code that does the makeshift precision adjusting (it has a parameter ADD in it).

[rpollack9974]

In this example, if you look how the consistency check is failing, on the screen is outputted the two distributions that are supposed to be equal. They seen exactly the same except that they are off by a single factor of 5!

Looks like the problem is somehow in how the valuations are being stored?

[rharron]

Oh my goodness, I see! Okay, I should be able to fix that pretty quickly. For some reason, I thought what was being printed by consistency_check was the difference between the two distributions...

[rpollack9974]

Ah. There must be some error in all of these shifting of valuations.

[rharron]

No, the problem appears to come from the line ret = self(D) The ordps seem to be lost at that point.

[rharron]

In the init of coeffmod_OMS_Family_element there was an if that should've been an elif. That made it so that if you coerced a family of distributions into a space of such the valuations would be lost... That's all it was. There still remains what I think is an actual precision issue, namely the problem with t not ending up with total measure 0. I think we need to somehow extend precision in one or many places (automorphy factor? CM?). In the old code, Rob, I remember you had just thrown a bunch of extra precision in computing the automorphy factor.

[rpollack9974]

can you post a specific example where this is failing? the one above now works.

[rharron]

The one above works because the ad hoc parameter ADD (defined on line 108 of modsym_OMS_families_space.py) is set to a high enough value. If you set ADD to 0 (which is like it not being there at all), the above example still fails I think (or at least fails most of the time). If the above example still works after that, then just switch to p = 3 and that should definitely fail.

[rpollack9974]

OK...I'm looking back now at the random function for a single OMS (i.e. not a family). In that function, I'm trying to understand the corresponding M_in.

It looks like M_in is something like M + log_p(M) + 3 + ord_p(k) + ord_p(c)

The k term is there because we divide by k. The c term is there because we divide by c. The log_p(M) + 1 term is there because of the difference equation.

Why is the extra +2 there? Is that just for safety sake?

[rharron]

Weird... The docstring I wrote says that given input precision M_in, the output precision is

M = M_in - ceil(log_p(M)) - 3 - ord_p(k),

so for some reason I thought there was a minus 3 there and that presumably is why there is a +3 in the formula for M_in. I don't know why I thought there was a -3, it seems like I was recording in the docstring something I knew to be true. hmmm.

[rharron]

But anyway that extra +2 there doesn't affect the bug that remains in this issue: the bug occurs before solving the difference equation; it's a bug in cancelling out the total measure of t.

[rpollack9974]

But the function _prec_for_solve_diff_eqn isn't really only about the difference equation. We are including the term ord_p(k) because we divide by k in trying to kill off the total measure of t.

[rharron]

Sure, but what I'm saying is that suppose you have pick some precision M_in (whatever number you want) and you have chosen random values at all the generators except D_infty and you sum up along the boundary and you get t, THEN what happens in the code between lines 140 and 189 should adjust t so that its total measure is 0. None of that code has anything to do with the provenance of M_in (i.e. it doesn't care that actually there was some initial M and M_in was defined to be some bigger number than M). Know what I mean? The value of M only starts mattering again after we solve the difference equation because then we need to know that M_in was big enough that the solution to the difference equation is accurate enough. But the process of adjusting the total measure of t doesn't and shouldn't care about M, only M_in.

[rharron]

Or wait maybe I see what you mean.

[rharron]

Yeah, I think I see what you're saying. One component of M_in is how much precision we think we will lose in the process of adjusting t, so that amount of precision should be removed from t before solving the difference equation. So, are you suggesting that the problem is that there might be some systematic loss of precision coming from dividing by K[1] that we haven't accounted for?

[rpollack9974]

Exactly. I think the precision should work like this:

When k==0, the loss of precision comes from solving the difference equation and from dividing by K[1] and so

M_in = M + log_p(M) + 1 + ord_p(c)

Here the ord_p(c) comes from the proof of Theorem 7.1 in our paper. (I think the w-adic precision also drops by 1 here.)

When k!=0, the loss of precision comes only from solving the difference equation and so

M_in = M + log_p(M) + 1

This is because K[0] - 1 is invertible (again see proof of Thm 7.1). However, to know that this is true we should be choosing our gam0 so that if a is the upper-left corner of gam0 then a^k-1 should not be divisible by p. Do we do this???

[rpollack9974]

Also, I tried dropping that "extra" +2 term in prec.... in the OMS case and everything worked okay in a couple of examples. I haven't uploaded this change though.

[rharron]

But that is how our precision currently works: in the line

M_in = ZZ(1 + M[0] + ceil(ZZ(M[0]).log(p))) + gam_shift + ADD

the variable gam_shift is the ord_p(c). In a recent update to the code I added functionality to the ManinRelations class that for each k and p finds a nice gamma such that, as the case may be (depending on k is 0 or not), either ord_p(c) is minimal or omega(a)^k != 1. (I checked a bunch of different p, k, and N, and ord_p(c) was always 1 btw). So, there seems to be some unaccounted for precision loss, at least when k=0.

[rharron]

Okay, so I've tried to isolate the problem in the adjustment of the total measure in the random families code. This is how far I've gotten so far. Suppose you have some nice gam0 and K[1] is the coefficient of z in its automorphy factor (here the weight is r=0). And suppose you have a (family of) distribution(s) t whose total measure is -K[1]. Then, err = -t.moment(0) / K[1] gives 1. Letting err = (0, err, 0, ...), the zeroth moment of err * gam should be K[1]. The code below shows that, despite the initial p-adic precision being 11 in a p-adic ring with precision 20, the zeroth moment of err * gam is equal to K[1] only up to O(5^9) representing a loss of precision of 5^2. The gam_shift in the example should be 1, so the two need to be equal up to O(5^10). Here's the code to reproduce this:

sage: p = 5 sage: a, c, k, M_in, R = -4, 15, 0, 11, PowerSeriesRing(ZpCA(p, 20), 'w', default_prec=11) sage: from sage.modular.pollack_stevens.families_util import automorphy_factor_vector sage: K = automorphy_factor_vector(p, a, c, k, None, M_in, M_in, R) sage: t0 = - K[1] sage: err = -t0 / K[1] sage: err 1 + O(5^20) + O(5^21)_w + O(5^22)_w^2 + O(5^23)_w^3 + O(5^23)_w^4 + O(5^24)_w^5 + O(5^25)_w^6 + O(5^26)_w^7 + O(5^26)_w^8 + O(5^27)_w^9 + O(5^28)_w^10 + O(w^11) sage: CM = FamiliesOfOverconvergentDistributions(0, base_coeffs=ZpCA(p, 20), prec_cap=[M_in,M_in]) sage: R = CM.base_ring() sage: v = [R(0), err] + [R(0)] * (CM.length_of_moments(M_in) - 2) sage: err = CM(v) sage: S0 = CM.action().actor() sage: gam0 = S0([-4,-3,15,11]) sage: err_gam0 = err * gam0 sage: K[1] - err_gam0._moments[0] O(5^11) + O(5^11)_w + (2_5^10 + O(5^11))_w^2 + (5^9 + 2_5^10 + O(5^11))_w^3 + O(5^11)_w^4 + (5^10 + O(5^11))_w^5 + (4_5^10 + O(5^11))_w^6 + (3_5^9 + 5^10 + O(5^11))_w^7 + (4_5^10 + O(5^11))_w^8 + (4_5^10 + O(5^11))_w^9 + (4_5^10 + O(5^11))*w^10 + O(w^11)

[rpollack9974]

Yikes. So the problem appears to simply be in the action on families.

The distribution err satisfies err(z^j) = 0 unless j=1 in which case err(z) = 1.

So (err * gamma)(1) by definition equals K[1] -- the above code appears to be illustrating that this is not happening in the 11-th filtration module. Let's hope this is just a coding error and not a theoretical problem.

[rpollack9974]

Just looked at the the code, I can no longer make sense of how the action is computed, but I do think the problem must be there.

[rpollack9974]

I guess the root of the problem is in the computation of the automorphy factors. This must be a bug (I expect the difference to be divisible by p^11, yes?)

sage: K11 = automorphy_factor_vector(p,a,c,k,None,11,11,R) sage: K12 = automorphy_factor_vector(p,a,c,k,None,12,12,R) sage: K11[1] - K12[1] O(5^21) + (3_5^11 + 3_5^12 + 4_5^13 + 3_5^14 + 4_5^15 + 4_5^16 + 4_5^17 + 4_5^18 + 4_5^19 + 4_5^20 + O(5^21))_w + (2_5^10 + 3_5^12 + 4_5^13 + 4_5^14 + 4_5^15 + 4_5^16 + 4_5^17 + 4_5^18 + 4_5^19 + 4_5^20 + 4_5^21 + O(5^22))_w^2 + (5^9 + 5^10 + 4_5^12 + 4_5^14 + 4_5^15 + 4_5^16 + 4_5^17 + 4_5^18 + 4_5^19 + 4_5^20 + 4_5^21 + 4_5^22 + O(5^23))_w^3 + (5^10 + 2_5^11 + 4_5^12 + 4_5^13 + 3_5^14 + 4_5^15 + 4_5^16 + 4_5^17 + 4_5^18 + 4_5^19 + 4_5^20 + 4_5^21 + 4_5^22 + 4_5^23 + O(5^24))_w^4 + (5^10 + 5^11 + 2_5^12 + 2_5^13 + O(5^24))_w^5 + (4_5^10 + 3_5^11 + 3_5^12 + 4_5^14 + 4_5^15 + 4_5^16 + 4_5^17 + 4_5^18 + 4_5^19 + 4_5^20 + 4_5^21 + 4_5^22 + 4_5^23 + 4_5^24 + O(5^25))_w^6 + (3_5^9 + 3_5^10 + 2_5^13 + 4_5^14 + 4_5^15 + 4_5^16 + 4_5^17 + 4_5^18 + 4_5^19 + 4_5^20 + 4_5^21 + 4_5^22 + 4_5^23 + 4_5^24 + 4_5^25 + O(5^26))_w^7 + (2_5^10 + 3_5^12 + 3_5^13 + 3_5^14 + 4_5^15 + 4_5^16 + 4_5^17 + 4_5^18 + 4_5^19 + 4_5^20 + 4_5^21 + 4_5^22 + 4_5^23 + 4_5^24 + 4_5^25 + 4_5^26 + O(5^27))_w^8 + (4_5^10 + 4_5^11 + 4_5^12 + 3_5^13 + 4_5^14 + 4_5^15 + 4_5^16 + 4_5^17 + 4_5^18 + 4_5^19 + 4_5^20 + 4_5^21 + 4_5^22 + 4_5^23 + 4_5^24 + 4_5^25 + 4_5^26 + 4_5^27 + O(5^28))_w^9 + (4_5^10 + 3_5^11 + 4_5^12 + 5^13 + 4_5^14 + 4_5^15 + 4_5^16 + 4_5^17 + 4_5^18 + 4_5^19 + 4_5^20 + 4_5^21 + 4_5^22 + 4_5^23 + 4_5^24 + 4_5^25 + 4_5^26 + 4_5^27 + O(5^28))_w^10 + (4_5^10 + 3_5^12 + 4_5^13 + O(5^29))_w^11

[rharron]

I found it, give me a sec. (Well, the adjusting the total measure works, but then the consistency check fails).

[rharron]

So, the answer to your question is: "No, the difference should not be divisible by p^11". Instead:

sage: p = 5 sage: a = -4 sage: c = 15 sage: k =0 sage: DD = FamiliesOfOverconvergentDistributions(0, base_coeffs=ZpCA(p, 20), prec_cap=[8,8]) sage: R = DD.base_ring() sage: K11 = automorphy_factor_vector(p,a,c,k,None,DD.length_of_moments(11),11,R) sage: K12 = automorphy_factor_vector(p,a,c,k,None,DD.length_of_moments(12),11,R) sage: from sage.modular.pollack_stevens.coeffmod_OMS_families_element import _padic_val_of_pow_series sage: _padic_val_of_pow_series(K11[1]-K12[1]) 13

And that was the bug basically. In the definition of K in random_element, automorphy_factor_vector was being passed M_in instead of CM.length_of_moments(M_in).

So now, in the initial example with both p = 5 and p = 3, the adjustment of the total measure of t works. Unfortunately, the _consistency_check still fails!

[rharron]

Hmm, this may be a problem with equality... -t - (f[id]MR.gammas[id] - f[id]) is 5^8 \ () + O(w^8).

[rharron]

Yeah, it's a bug in how sage compares power series. One of the terms has some terms with O(5^19) * w^7, while the other doesn't have those terms. And sage then thinks these two power series are not equal!

[rpollack9974]

OK...I think I see the problem. The function logpp_binom returns the power series p^n (log_gamma(1+pz) choose n) which has integral coefficients.

When computing the automorphy factors one takes this power series and plugs in a+bz in for z for some a and b. Unfortunately, a and b are not necessarily divisibly by p. Since we only computing up to finitely many powers of z, we are just dropping terms of the form (a+bz)^N for large N. This is bad because these terms contribute non-trivially to say the constant term of the power series.

Fortunately, the coefficients in front of (a+bz)^N become more and more divisible by p as N grows. If these coefficients were actually divisible by p^N then all would be fine. However, they are only divisible by p^(N c_p). And so we need to add more accuracy in z to get the lower coefficients correct.

[rharron]

I mean, maybe you found some other problem, but I was just about to fix this issue....(!)

[rpollack9974]

Just pushed a posible fix to this by increasing the precision in z by a factor of (p-1)/(p-2) when computing logpp_binom.

Does this fix the problem?

One thing to mention, in my fix I didn't think about whether the scaling goes above the max precision of the ring we are using.

[rpollack9974]

Yikes! I just saw all of your comments.

On Tue, Jun 18, 2013 at 10:44 PM, rharron notifications@github.com wrote:

I mean, maybe you found some other problem, but I was just about to fix this issue....(!)

— Reply to this email directly or view it on GitHubhttps://github.com/lalitkumarj/OMSCategory/issues/50#issuecomment-19659480 .

[rpollack9974]

Now. I'm all confused. Maybe your fix is what is making things work. I'll going to take out my fix (in my code) and see if anything breaks.

[rpollack9974]

When I removed the increased precision for z your code above from your comment starting with "Okay, so I've tried to isolate the problem" still fails. So probably we found two bugs.

Is it working now??

[rharron]

That code above is flawed in the same way that I explained in this comment (https://github.com/lalitkumarj/OMSCategory/issues/50#issuecomment-19657204) above: In the definition of K, it should be K = automorphy_factor_vector(p, a, c, k, None, CM.length_of_moments(M_in), M_in, R) (i.e. there's a CM.length_of_moments(M_in) instead of just M_in). This increases the precision the right amount.

Your change also increases the precision which explains why it would also make the above example work. What you found may indeed be a bug, but it's not this bug!

I'll undo your change for now, but let's keep it in mind in case something comes up.

[rharron]

Or wait, it looks like you're making more related changes.

[rpollack9974]

Got it! I think it's the same fix essentially -- going from moments to length of moments is like scaling by (p-1)/(p-2). In any case, please do revert my recent changes.

So is there still a problem with random symbols?

[rharron]

Well, I thought there were no problems left, but the following now crashes:

sage: p = 3 sage: N = 15 sage: M = 4 sage: var_prec = M sage: DD = FamiliesOfOverconvergentDistributions(0, base_coeffs=ZpCA(p, M), prec_cap=[M,var_prec]) sage: MM = FamiliesOfOMS(N, r, coefficients=DD, sign=-1)

sage: Phis = MM.random_element()

ValueError Traceback (most recent call last)

in () ----> 1 Phis = MM.random_element() /Users/rharron/sages/sage-5.9/local/lib/python2.7/site-packages/sage/modular/pollack_stevens/modsym_OMS_families_space.pyc in random_element(self, M) 203 #print "M_in =", M_in 204 #print "t[0] =", t.moment(0) --> 205 mu = t.solve_diff_eqn() 206 mu_val = mu.valuation() 207 #print "Shift:", shift /Users/rharron/sages/sage-5.9/local/lib/python2.7/site-packages/sage/modular/pollack_stevens/coeffmod_OMS_families_element.pyc in solve_diff_eqn(self) 754 verbose("Abs_prec of self: %s"%(abs_prec)) 755 verbose("mus_abs_prec before reduction: %s"%(mus.precision_absolute())) --> 756 mus = mus.reduce_precision_absolute([abs_prec - abs_prec.exact_log(p) - 1, var_prec]) 757 verbose("mus_abs_prec after reduction: %s"%(mus.precision_absolute())) 758 return mus.normalize() #Is it necessary to normalize? /Users/rharron/sages/sage-5.9/local/lib/python2.7/site-packages/sage/modular/pollack_stevens/coeffmod_OMS_families_element.pyc in reduce_precision_absolute(self, new_prec) 672 var_prec = aprec[1] 673 if p_prec > aprec[0] or var_prec > aprec[1]: --> 674 raise ValueError("Precisions specified must be less than current precisions.") 675 if p_prec == aprec[0] and var_prec == aprec[1]: 676 return self ValueError: Precisions specified must be less than current precisions.

[rharron]

The error is from there being precision lost in solve_diff_eqn. (Note, I'm not using your edits here, I haven't tried it with your edits yet).

[rharron]

I think I've found this bug now. The M we take in solve_diff_eqn (in families) needs to be bigger. Basically, it needs to account for the natural loss of precision in solving the equation in non-families. I'm working on it.

[rharron]

Yet another problem with adjusting the total measure (and this is with your changes, too):

sage: p = 7; M = 4; var_prec = M sage: N = 21 sage: DD = FamiliesOfOverconvergentDistributions(0, base_coeffs=ZpCA(p, M), prec_cap=[M,var_prec]) sage: MM = FamiliesOfOMS(N, r, coefficients=DD, sign=-1)

sage: Phis = MM.random_element()

ValueError Traceback (most recent call last)

in () ----> 1 Phis = MM.random_element() /Users/rharron/sages/sage-5.9/local/lib/python2.7/site-packages/sage/modular/pollack_stevens/modsym_OMS_families_space.pyc in random_element(self, M) 206 verbose("abs prec of t before solving: %s"%(t.precision_absolute()[0])) 207 verbose("output will have precision (should have %s): %s"%(M[0], t.precision_absolute()[0] - t.precision_absolute()[0].exact_log(p) - 1)) --> 208 mu = t.solve_diff_eqn() 209 verbose("mu's abs prec: %s"%(mu.precision_absolute())) 210 mu_val = mu.valuation() /Users/rharron/sages/sage-5.9/local/lib/python2.7/site-packages/sage/modular/pollack_stevens/coeffmod_OMS_families_element.pyc in solve_diff_eqn(self) 719 return CoeffMod_OMS_Families_element(V([]), self.parent(), ordp=(abs_prec - ZZ(abs_prec).exact_log(p) - 1), check=False, var_prec=var_prec) 720 if not self._unscaled_moment(0).is_zero(): --> 721 raise ValueError("Family of distribution must have total measure 0 to be in image of difference operator; total measure is %s"%self.moment(0)) 722 M = ZZ(len(self._moments)) 723 if M == 2: ValueError: Family of distribution must have total measure 0 to be in image of difference operator; total measure is O(7^5) + (2_7 + 7^2 + 3_7^3 + 7^4 + O(7^5))_w + (5_7 + 2_7^3 + 5_7^4 + O(7^5))_w^2 + (3_7 + 6_7^2 + 6_7^3 + 4_7^4 + O(7^5))_w^3 + O(w^4) I'm going to bed.

[rharron]

Well, the problem might be that gam0 is a three-torsion generator. When you worked out how to fix the total measure for families you didn't address the case of torsion generators, so I just guessed at what it might be given what we do for non-families. Guess I was wrong! (Now, I'm going to bed).

[rpollack9974]

Yup. Even in the OMS case the three torsion adjustment is different.

The OMS code looks like:

        if not g0 in manin.reps_with_three_torsion():
            err = -t.moment(0)/(chara * k * a**(k-1) * c)
        else:
            err = -t.moment(0)/(chara * k * (a**(k-1) * c + aa**(k-1) * cc))

I'm trying to fix this now.

[rpollack9974]

OK. The code above now works. But I didn't deal with changing gam_shift etc in this new case.