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.

[rudolph-git-acc]

Below are two graphs on the derivative bounds vs actuals, used for the Barrier evaluation at X=83852 \pm 0.5, y=0.2 at each t-step (xy-rectangle). We have apparently used bounds that are so conservative, that I had to put logs on both to adequately visualise them... :)

Let me know if you like to make any changes. The d/dt actuals are pretty smooth, however the d/dz-actuals show a dip at small t.

I'll now go through Sections 8 and 9 of the write-up and provide the missing numerical data.

UPDATED GRAPHS

[teorth]

These pictures are great! My only suggestion is to make more explicit in the second set of pictures that the y-axis is log-scale (e.g., by labeling it by log|\partial f / \partial z|, etc. By the way I prefer \partial to \delta for partial derivatives.) If you could add them to the writeup while also going through the sections that would be great.

Another picture that would be worthwhile is a picture of a mesh (basically a dotted rectangle in x-y space, together with its image under f_t in the complex plane, say for t=0, with the polygon connecting the points drawn, hopefully in such a way that it becomes obvious that the winding number vanishes... something like Kalpesh's plots in https://terrytao.wordpress.com/2018/03/18/polymath15-sixth-thread-the-test-problem-and-beyond/#comment-494202 .

Finally it would be nice to have some graphic or table displaying the values of t used and the numebr of mesh points used per value of t. In particular it could illustrate that it is the t=0 case that requires the most precision in both the mesh size and the t increment, and one can use much larger steps and mesh sizes for larger values of t (maybe a scatter plot of the choices of t versus the mesh size?)

In parallel, I can start looking at the code used to establish zero free regions right of the barrier using Euler mollifiers. What is the latest version of the code? I found the PARI/GP code to be quite readable, haven't tried looking at ARB code yet.

[rudolph-git-acc]

Have updated the graphs in the previous post and will include them in the write-up later. Will first start working on the other two graphs. I recall that KM had even made an animated version of these polygons that was highly addictive to watch and even induced palpitations when it occasionally came closer to rounding the origin :)

The latest pari/gp version to explore that Na...Nb is zero-free is this one: https://github.com/km-git-acc/dbn_upper_bound/blob/master/dbn_upper_bound/pari/abbeff_largex_bounds.txt

Best to start from the while loop at the bottom.

Pari/gp is indeed a quite readable symbolic language and we used it mostly for testing the various formulae and scripts. All 'production' runs were done in ARB that is much faster, but also more cumbersome to write & read. Some have even compared it to assembler, however once you get used it, it's not bad at all. So far I have been quite impressed by the richness and speed of this relatively new mathematical software (e.g. the speed of evaluating zeta(s) at 'higher heights' is truly amazing).

[rudolph-git-acc]

Here is a plot of the number of mesh points required for each t. As expected most mesh points 11076 need to be evaluated at t=0 dropping to 56 at t=0.2.

The winding number plot has been really 'winding' me up, since it doesn't seem to rotate as much as KM's plots that were done for a larger x-ranges in the y=t=0.4 days. Then suddenly realised that the x-range at the Barrier is so small and we have also selected a barrier location that has a peak value in ABB-eff, that a single rotation (not around the origin) could actually all there is in this domain. Is my reasoning correct or did I do something wrong in making the plot? (the plot currently connects all sequential x, y ABB_eff values of the mesh points on the xy-rectangle at t=0). Grateful for your steer on this.

[teorth]

Thanks for the plots! They would be great additions to the writeup.

The lack of oscillation is a little strange. Possibly the choice of barrier location is indeed reducing some of the oscillation, but it is surprising that so much of it is gone. One could of course try drawing a similar graph at a randomly selected nearby location (e.g. X = 6 x 10^{10}) to see if there is a difference. One can use fewer mesh points for a faster plot here as this is just an exploratory plot (though if there is a significant difference it might be worth putting both plots in the paper to demonstrate the utility of choosing a good barrier location). One could also try reproducing KM's plots (i.e. setting t=0.4 etc.) using the same code to see if the problem is code related.

[rudolph-git-acc]

Have done some checks on the (lack of) oscillation in the polygon plots and it indeed all seems to boil down to the small size of delta x at the Barrier. When mesh points are calculated across a range of e.g. x...x+30, the oscillations become clearly visible. Also checked what would happen when we'd move the Barrier to a less optimal location, but do get a similar (i.e. one rotation) graph without any oscillation.

Will incorporate all graphs in the write-up shortly after Christmas!

[teorth]

OK, great! Looking back at Kalpesh's plots it looks like he was covering the range 0 < x < 300 which explains all the oscillation. Each term in the bn^t / n^{s*} expansion has a wavelength of something like 4 \pi / log n in the x variable, so it makes sense that there is relatively little oscillation on the scale of a single barrier. In retrospect it is also consistent with the observation you made earlier that our bounds for the derivatives are way too conservative.

[km-git-acc]

Prof. Tao and Rudolph, firstly sorry for not being active for some time. Rudolph has indeed made all the clarifications. Just to reiterate, H had to be chosen as a divisor of N in the barrier_multieval script to avoid tail end behavior, which was not explicitly mentioned in the code. However, we had then shifted to making H=N in the later script barrier_multieval_t_agnostic which made things simpler. We have been doing a lot of prototyping work in Parigp, but Rudolph had become quite comfortable in Arb which is often way faster, so we have been doing production runs through those. We have been synchronizing the logic in both and comparing their respective outputs to ensure we remain on track.

Some of the extreme production runs in the unverified RH range e.g. X near 10^20 were done through grid computing where volunteer computers (a few hundred in number) participated. What is generally the recommended way to acknowledge them in the writeup?

Also, Rudolph has mentioned he will be incorporating the newer graphs in the writeup, so will wait for that before making any changes.

[teorth]

We should definitely acknowledge the grid computing volunteers in Section 11 of the writeup; also there is a web page to list participants at http://michaelnielsen.org/polymath1/index.php?title=Polymath15_grant_acknowledgments, one can add a section there for the grid computing contributors. (The two of you should add your own contact information there, and any funding support if any.)

[rudolph-git-acc]

Have now incorporated the Barrier graphs together with some texts in the write-up. Pictures still jump a bit across the pages, but this is probably best fixed in the final review step. Let us know which other data is needed.

[km-git-acc]

Great. I will make the changes regarding the grid computing contribution as well. EDIT: Have added the contribution from grid computing to section 11. Also added a section on the Grant acknowledgements wiki page.