Writing up the barrier computation

[teorth]

Looking through the writeup, the main things still left to do are to describe the numerical verification of the barrier (Claim (ii) in Section 9 of the writeup) and of the finite zero free region to the right of the barrier (Claim (iii)). We have notes of sorts for claim (iii), but right now claim (ii) only has the bound on the derivatives.

In the most recent blog post (https://terrytao.wordpress.com/2018/09/06/polymath15-tenth-thread-numerics-update/) there is a brief description of the method - selecting a non-adaptive mesh for the rectangle at each time t, bounding the derivative by a conservative bound, and using Taylor expansions to compute the mesh - but it isn't at the level of detail needed to describe what we do precisely. Would one of you be willing to write a more detailed description of the verification to put in the writeup? Alternatively, if you could provide a link to the code used, and maybe just give an informal summary of the verification procedure, I can try to read the code directly and write up a description.

[km-git-acc]

@teorth @rudolph-git-acc

I will try to describe in more detail the code by tomorrow.
There are both parigp and arb versions of the code, and the former should be easier to follow.

https://github.com/km-git-acc/dbn_upper_bound/blob/master/dbn_upper_bound/pari/barrier_multieval_t_agnostic.txt

While going through the code, it's helpful to start from the last block (between the two getwalltime functions).
We run through increments of t, updating ddzbound and ddtbound at every step using the ddzboundint and ddtboundint functions (which replace the summation in the paper with integrals). Mesh points are decided using the ddzbound value. Post calculating derivatives, we supply the storedsums matrix and required parameters to the abbeff_multieval_symmetric_rectangle function, which evaluates |a+b| at the mesh points, and computes the winding number. Finally, key metrics at each t step are printed, and the next time step is chosen.
(Number of taylorterms used for generating the storedsums matrix is currently set to 50 in the code, but it can be determined through trial and error through the ttermmagnitude3 function earlier in the code. The arb code has a automatic way of doing this).
(Links to the Arb codes
https://github.com/km-git-acc/dbn_upper_bound/blob/master/dbn_upper_bound/arb/StoredSumSinglematv1.c

https://github.com/km-git-acc/dbn_upper_bound/blob/master/dbn_upper_bound/arb/TloopSinglematv2.c)

[km-git-acc]

@teorth @rudolph-git-acc

Attempting a more detailed explanation

1. At time t, df/dz is calculated according to Lemma 9.3, and the number of mesh points is set as ceiling(df/dz). Since the barrier location has been chosen in a way such that the minimum of f_t on any mesh point is expected to be above 1 (normally much above 1), the distance between consecutive mesh points is 1/ceiling(df/dz), hence f_t should not drop to 0 between the mesh points. (When f_t on actual calculation on the mesh points turns out to be above 1, this expectation is confirmed).

df/dz is calculated using numerical integrals instead of summation : summand_(n=1) + integral_(n=1 to N), which gives a conservative estimate of the sum (since summation becomes a significant bottleneck at large N values, while the integral calculation is instantaneous). The numerical integration is stable since df/dz shows non oscillatory behavior with x,y and t.

2. Again at time t, once f_t is calculated at all the mesh points, the winding number is computed and the actual minimum f_t value observed on the mesh points is stored. Using Lemma 9.3 and numerical integration as above, df/dt is calculated. The next t value is found using
   t_current + (f_t_mesh_min - 0.5)/(df/dt)

3. A storedsum matrix of Taylor coefficients (of a z,t polynomial) has been generated earlier outside loop 1-2 to be supplied for the f_t evaluations at each time step.
   For the f_t evaluations, firstly the z coordinates of all the mesh points are generated as a vector.
   The storedsums matrix is semi evaluated by passing it the t parameter, thus converting it to a z polynomial (generated as a vector of coefficients).
   Post this, using vector and matrix operations, the f_t values are generated as a vector.

Point 3) is present as the function abbeff_multieval_symmetric_rectangle in the parigp code.

[teorth]

Thanks for this! I've ended up being overwhelmed with other urgent duties at present, but will get around to incorporating this description into the writeup soon.

[teorth]

Beginning to read through the pari code, starting with barrier_multieval.txt . My first question is for the barrier_strip_search routine. It seems that the default number of primes used here is 5 (so I guess one uses the primes 2,3,5,7,11). In the writeup for some reason (on page 38) it says the primes up to 29 are used, but this could perhaps be a typo. Do either of you remember how many primes were used to locate the strip 6 x 10^10 + 83952 \pm 0.5? I guess basically one just needs to review the discussion on page 38 for accuracy.

[rudolph-git-acc]

We have definitely used 30 primes for the search. Note that the default settings in: barrier_strip_search(bfile,X,numx=10^5,nprimes=5) always get overruled by the parameters set when the function is called.

Attached is the ARB-output of the latest Barrier search that we did to find the Barrier locations for the 'Na = Nb' runs.

BarrierlocationselectionNoTriangle.txt

[teorth]

Just to clarify ... did you use the first 30 primes, or the primes up to 30 (of which there are 10)?

I have now read through the abbapproxarr routine. There may be a tiny issue with it: this routine approximates a sum from 1 to N by subdividing into intervals of length H; as implemented it subdivides into intervals [ n_0(v) - H/2+1, n_0(v) + H/2 ] where n_0(v) = H/2 + (v-1) H and v runs from 1 to (N-H/2) / H. (I guess H must then be an even integer, but you set it to 10, which is OK.) But the problem is that this doesn't quite match the range from 1 to N if N is not a multiple of H (and in our case N = 69098 and H=10 so I think the sum undercounts by 8 terms in this case if I did the math correctly). Presumably the contribution of these last few terms is very small so it should not seriously affect the calculations, but the code may need to be fixed (perhaps by truncating N to a multiple of H and evaluating the final few terms separately, presumably this can be done naively in a reasonable amount of time).

(ADDED LATER) I am now rewriting Section 8 of the writeup to match the actual implementation in the pari code. I wanted to point out that direct numerical calculation of all the expression f_t(x+iy) without the stored sums method would be computationally prohibitive. For this I would need to know the rough sizes of N (the length of the sum), T (the total number of times t one evaluates at), and S (the number of values of s = x+iy one evaluates at). I know that N = 69098, but what are T and S? I guess perhaps S is t-dependent, so maybe the more useful statistic is the total number M of pairs (t,s) one needs to evaluate f_t(x+iy) at.

[teorth]

Actually, it may be easier to just estimate the final few terms as part of the error term instead of computing them. If the size of the final term is much less than H=10 times the margin of error then this would be safe. Do we know roughly how large these quantities are?

(ADDED LATER) If we overcount the 69098 terms rather than undercount then the error consists of just 2 terms rather than 8 so that should lead to smaller errors. Alternatively one could tweak the abbapproxarr routine by reducing the range of h in the n0sumsa and n0sumsb calculation for the final value of v.

(ADDED YET LATER) In any case, I would need to know what the margin of error is (how far away the output of abbapproxarr is) in order to rigorously show that the effect of the tail of the Taylor expansion (the values of expo beyond expterms, which seems to be set to 4 in this code) is negligible. As I recall these errors were extremely small numerically, but we would need to show this rigorously. (There may be an issue here from the very small values of n or n0, where Taylor expansion is not that precise.)

[rudolph-git-acc]

It is definitely the first 30 primes.

Your analysis of the subdivision in the abbapproxarr loop looks correct and this should be easy to fix. I expect the error of the missing 8 terms to be small. Re-computing the 0.22 case is not an issue since a full rerun could be done in less than 2 minutes and better to get it fully correct (will do it tomorrow).

On your N, T, S question; agree computations without the stored sums will become prohibitive. I am afraid I would need some help from KM to produce T & S (he is sorting out some personal issues at the moment but will return soon). Will try to take a shot at it tomorrow.

[rudolph-git-acc]

On your (ADDED YET LATER) see also the Nov. 19 comment above from KM:

(Number of taylorterms used for generating the storedsums matrix is currently set to 50 in the code, but it can be determined through trial and error through the ttermmagnitude3 function earlier in the code. The arb code has a automatic way of doing this).

In ARB we indeed automatically establish the required Taylor/Exp terms based on the number of digits accuracy required (we always used 20 digits for our computations, except at 10^20 and 10^21 runs where we reduced to 10 digits). To check this all works correctly, we did some manual checks on specific ABB_eff values that are 'recreated' from the stored sums matrix during the Barrier evaluation process.

[rudolph-git-acc]

After rereading some of the old discussions with KM, I now see that I have informed you incorrectly about the number of primes used for the Barrier_strip_search calculations. My apologies.

I now found this summary from KM that clearly indicates the first 10 primes have been used (and not 30). Let's therefore stick to the number 29 (=10 primes) that is also fully aligned with KM's comments in his Sep 24 post in the 'Corrections to figures' thread.

(P.S. KM incorrectly mentions 83951.5 instead of 83952 /pm 0.5 in the copied below):

[rudolph-git-acc]

Below are some data that could help show that a direct evaluation of f_t(x+iy) (=ABB_eff) at each mesh point in the Barrier evaluation process is dramatically slower than using the stored sums approach. At this level (6xE10) direct evaluation is maybe not yet entirely 'prohibitive' for computations, but will quickly become so for each factor 10 increase of x.

Took a deeper dive into the abbapproxarr loop and although it should be easy to fix the error by adding a single line (the undershoot was 13):

abbapproxarr(X,y,t,H)={
    N = floor(sqrt((X+Pi*t/4)/(4*Pi)));
    numn0 = (N - H/2)\/H;
    n0 = vector(numn0,v,H/2+(v-1)*H);
>  n0[numn0] = N;
(...)

I then suddenly realised that barrier_multieval.txt is an older script that has become obsolete after we moved to barrier_multieval_t_agnostic.txt (the one KM refers to at the beginning of this thread in his Nov 19 comment and can be found here: https://github.com/km-git-acc/dbn_upper_bound/blob/master/dbn_upper_bound/pari/barrier_multieval_t_agnostic.txt)

In this later version H was simply set equal to N and n0 to H/2, so the subdivisions disappeared (as did the error) and the script got simpler. This is also the logic used in ARB and should therefore be the code to use for your walk through.

[teorth]

Thanks for this! This is indeed simpler and I will adjust the writeup accordingly.

Technically, one has to rigorously bound the error term in the taylor expansion to 50 terms, but this would just add some uninteresting and lengthy computations to the paper and it may be better just to give convincing numerical evidence that the error term is negligible. Would it be possible to give some table or graphics that show, for a given value of t (probably t=0 is the worst case), the rectangular mesh used, the evaluation of both the Taylor expanded abeff and the exact value of abeff on this contour (hopefully in a way that makes the visually obvious that the winding number is zero), and some numerical data or plot on the error between the two? Given that the full plot for all ts takes 78.5 hours, I'm assuming that if one just wants to plot the t=0 data that one could actually compute that in a reasonable amount of time. If the relative error is many orders of magnitude better than what we need then that should be convincing enough, I think.

[teorth]

OK, I have now written up a description of the barrier computation strategy in Sections 8 and 9 of the writeup. It needs to be checked for accuracy, and actual numerical data added (e.g. how many mesh points are needed, what the derivative bounds look like, etc.) but I think we are getting there. Basically we should provide enough numerical data to be convincing that even with very conservative error estimates we can verify the barrier in a short amount of time, at least at the X ~ 6 * 10^10 level.

After this is done, I'll take a look at the verification of Dirichlet series stuff that covers larger values of x. You've mentioned that for some ranges of parameters we can actually skip this step - I take it though that for the original unconditional bound of 0.22 this step is still needed?

[rudolph-git-acc]

Yes, such a graph at t=0 is feasible. There are 11076 mesh points to be processed at t=0 and I could compute the exact ABB_eff-value for each of them and then plot the difference (i.e. the error) that should be better than 20 digits. Will start working on it right away.

You are indeed correct that for x=6x10^10, y=t=0.2 _> DBN = 0.22, we unfortunately can't skip the 'lower Lemma bound > error upper bound' verification step. The reason is that the x has been fixed by the height of Platt's RH-verification work and therefore can't be optimised by playing with y_0 and t_0 to obtain a positive Triangle bound at the Barrier location.

Also happy to provide numerics on the derivative bounds, if you could specify what you like to see in that area.

[rudolph-git-acc]

Here is the graph of the differences between Taylor expanded and exact evaluation of ABB_eff for all 11076 all mesh points at t=0. Accuracy looks pretty good and is based 54 Taylor terms that where automatically established to yield 20 digits accuracy.

The sides of the xy-rectangle are 'walked' in the following sequence:

(83951.5, 0.2) -> (83951.5, 1) -> (83952.5, 1) -> (83952.5, 0.2) -> (83951.5, 0.2)

so it seems the largest 'swings' occur when y gets closer to 0.2.