some thoughts about the Odlyzko–Schönhage algorithm
David Bernier
evaluations of the Riemann zeta function on the critical
line was first presented around 1988:
Odlyzko, A. M.; Schönhage, A. (1988), "Fast algorithms for multiple evaluations
of the Riemann zeta function", Trans. Amer. Math. Soc. 309 (2): 797–809.
If I recall correctly, this article is now available without
a subscription in electronic format from the Web-site
of the journal Trans. Amer. Math. Soc., in the "Past Issues"
section.
It's apparent to me that the authors wanted to give a
detailed treatment of storage requirements and
upper bounds on time complexity. Also, a good discussion
of error analysis. I'm fairly confident that one of
their primary objectives was in part to verify the
Riemann Hypothesis over a wide range, and also to
accumulate a lot of data about normalized spacings between
consecutive non-trivial zeros going up the critical line
Re(s) = 1/2 in the upper half-plane. Some of
A. Odlyzko's computational work has been in "testing"
the Montgomery pair-correlation conjecture, which in its
simplest form has to do with normalized spacings of
consecutive non-trivial zeros of Riemann zeta,
arranged by increasing values of the imaginary part.
[ Riemann zeta is real-valued on the line Im(s)=0,
with a simple pole at s=1; so the
Schwarz reflection principle applies, and
Riemann zeta in the lower half-plane in
some precise sense mirrors zeta in the upper half-plane.]
I mention all this because I think one productive way
to read the paper is to concentrate more on the two
or three methods/transformations that allow this fast
computation.
(1) The grid. zeta is computed
at equally spaced points from 1/2 + N*i to
1/2 + N*i + c_N * i , where c_N is roughly
of the order sqrt(N). The spacing h*i
was made "largish", maybe about the length
of a Gram interval at height N (Im(s) ~= N),
so as to minimize the number of evaluations
of the Riemann-Siegel formula (up to perhaps
two or three small corrective terms, which
get more complicated for 4th, 5th and higher
corrective terms). The main part of
the Riemann-Siegel formula is the evaluation
of what's close to a complex trigonometric sum,
keeping in mind that the n'th trigonometric term
gets divided by sqrt(n);
Cf.:
< http://en.wikipedia.org/wiki/Z_function#The_Riemann-Siegel_formula >
BTW: from the value of Z(t), it's easy to obtain zeta(1/2 +i*t), t real.
At this point, I'll mention what I might want to do eventually.
With relatively large spacings between the grid points,
the interpolation formulas for points in between grid points
amount to rather sophisticated signal processing methods
( maybe sinc(.), Nyquist formula, and all that).
If grid points were 1/10 th or 1/20 th of a Gram interval,
it seems to me that simpler interpolation formulas might
exist. There have been computations of large values
of |zeta(1/2 +i*t)| (compared to t), notably by Tadej
Kotnik who published an article in Math. Comp. a few years
ago. He may have used Mathematica, I'm not sure.
Also, I don't know how Mathematica goes about computing
zeta on the critical line. Optimistically, one might hope
to outperform Mathematica by using an Odlyzko–Schönhage
type algorithm with grid points separated by about 1/10 th or a
a bit less of a Gram interval.
The Odlyzko–Schönhage algorithm doesn't improve on
Riemann-Siegel for single-point, wherever, sparse
evaluations of zeta. But in searching for extrema
of |zeta(1/2 + i*t)|, it (or a simpler but slower
variant as I sketched above, might be useful).
There are also conjectures on moments of zeta.
This involves integrating |zeta(1/2+i*t)|^k
over quite long intervals. The problem had
was that with k = 8 or so, one has to compose
with "little mountains" in the graph of Z(t).
See e.g.
< http://berniermath.net/ZlocalMax.jpg > .
(2) The Fast Fourier Transform. This is an essential
ingredient in the Odlyzko–Schönhage algorithm,
as it permits fast evaluation of zeta at
the grid points through the intermediary
of methods very similar to
the Fast Multipole Method (or methods)
used in potential theory in computational
physics.
Cf.:
< http://en.wikipedia.org/wiki/Fast_multipole_method >.
According to Wikipedia,
<< The FMM, introduced by Rokhlin and Greengard, has been acclaimed as one of
the top ten algorithms of the 20th century[2]. >>
One application in mathematical physics is the study of
n-body problems in celestial mechanics, with 'n' large.
With 1000 massive bodies, one might think that computing
the gravitational field at a random point in the neighborhood
of the 1000 bodies could involve summing up to
1000 space-dependent vectors, representing the acceleration
due to body #1, #2, ... #1000. Apparently, one can do better
and I believe this type of question is addressed by
the method(s) that were first given by Rokhlin and Greengard.
So as next item, I'll put fast evaluation of rational functions
( the FFT allows one to appeal to such methods).
(3) Super-fast evaluation of some sums of rational functions.
[ the Odlyzko–Schönhage method for this.]
[ this may be close in spirit to Fast Multipole Methods
from scientific computations in physics and
chemistry.]
-----
I really recommend the paper by Xavier Gourdon on
their "variation" on the theme of the
Odlyzko–Schönhage algorithm.
As regards super-fast evaluation of some sums of
rational functions, see for example Section
2.3.2 of Gourdon's paper (also section 2.3.1,
same topic).
Also, the first half-page of Section 2.3 of Gourdon's paper
shows how, using the discrete Fourier Transform, one
can easily go back and forth between computing
F(t) on the grid, and the rational function
f(z) defined in equation (9). It's worth keeping in
mind, I think, that k_1 - k_0 in (9) can easily
be over 10^6 at high heights ( ~ 10^12 to 10^13 or so).
Finally, the F(t) of Gourdon is not precisely the Riemann-Siegel
long sum, but rather a slightly truncated form of it.
Nevertheless, computation of F(t) without
re-casting it as a sum of rational functions as in
eqn. (9) would be far more time-consuming (perhaps
not outperforming Riemann-Siegel by much).
So my impression is that the vast increase in speed of
computation with the Odlyzko–Schönhage algorithm
is due in large measure to using very efficient methods
for computing some kinds of (long) sums of rational functions;
still, I'd be curious to know how much "band-interpolation techniques"
contribute to the ~ O(t log(t)) time-expense at height ~ t.
Here is a link to Gourdon's paper on the Web:
<
http://numbers.computation.free.fr/Constants/Miscellaneous/zetazeros1e13-1e24.pdf >
.
David Bernier
--
$ gpg --fingerprint davi...@videotron.ca
pub 2048D/653721FF 2010-09-16
Key fingerprint = D85C 4B36 AF9D 6838 CC64 20DF CF37 7BEF 6537 21FF
uid David Bernier (Biggy) <davi...@videotron.ca>
David Bernier
[...]
[...]
In the same area, there's the article:
V. Y. Pan, M. A. Tabanjeh, Z. Q. Chen, E. I. Landowne and A. Sadikou;
"New transformations of Cauchy matrices and Trummer's problem",
Computers & Mathematics with Applications, Volume 35, Issue 12,
June 1998, Pages 1-5 .
It's (c) Elsevier Science and not available at ScienceDirect, where
I read the abstract. The formulation of "Trummer's problem" by
someone whose name I forget was followed by the discovery
of the Greengard-Rokhlin method.
->
< http://comet.lehman.cuny.edu/vpan/pdf/pan150.pdf > (V. Y. Pan).
Interestingly, they cite the 1988 paper by Odlyzko and Schönhage (ref. no. 9).
David Bernier
This document based on a talk of Odlyzko is also worthwhile; he says
what the main ideas of the OS algorithm are in the abstract.
< http://www.dtc.umn.edu/~odlyzko/doc/arch/zeta.fn.supercomp.pdf >
The name of the file is zeta.fn.supercomp.pdf
Chip Eastham
> pub 2048D/653721FF 2010-09-16
> Key fingerprint = D85C 4B36 AF9D 6838 CC64 20DF CF37 7BEF 6537 21FF
Interesting topic, David. Odlyzko's talk precis you linked
above mentions that one possible improvement would use the
Greengard-Roklin algorithm for fast multipole method. Are
you aware if they explored this further?
regards, chip
David Bernier
>>> June 1998, Pages 1-5 .
>>
>>> It's (c) Elsevier Science and not available at ScienceDirect, where
>>> I read the abstract. The formulation of "Trummer's problem" by
>>> someone whose name I forget was followed by the discovery
>>> of the Greengard-Rokhlin method.
>>
>>> ->
>>> <http://comet.lehman.cuny.edu/vpan/pdf/pan150.pdf> (V. Y. Pan).
>>
>>> Interestingly, they cite the 1988 paper by Odlyzko and Schönhage (ref.
>>> no. 9).
>>
>> This document based on a talk of Odlyzko is also worthwhile; he says
>> what the main ideas of the OS algorithm are in the abstract.
>>
>> <http://www.dtc.umn.edu/~odlyzko/doc/arch/zeta.fn.supercomp.pdf>
>> The name of the file is zeta.fn.supercomp.pdf
>>
>> David Bernier
>> --
>> $ gpg --fingerprint david...@videotron.ca
>> pub 2048D/653721FF 2010-09-16
>> Key fingerprint = D85C 4B36 AF9D 6838 CC64 20DF CF37 7BEF 6537 21FF
>> uid David Bernier (Biggy)<david...@videotron.ca>
>
>
> Interesting topic, David. Odlyzko's talk precis you linked
> above mentions that one possible improvement would use the
> Greengard-Roklin algorithm for fast multipole method. Are
> you aware if they explored this further?
Hello Chip. I'm not aware either way (that they did/(did not) )
use something like the Greengard-Rokhlin method.
There's a pre-print at arXiv about numerical integrations related
to conjectures on the asymptotics of moments (as in probability theory)
of |zeta(1/2 +i*t|^(2k), k = 1, 2, 3 ... about 8.
The pre-print by Ghaith Hiary and Odlyzko is available from:
< http://arxiv.org/abs/1008.2173 > .
They write that Odlyzko took pre-computations stored on tape from
AT & T Labs - Research to Minnesota when he moved to a new place and
job.
These were OS pre-computations. These could be values of zeta
at grid points for various intervals of Im(s) at very high
heights. Integration is tricky because of very steep mole-hill (or valley)
regions in the graph of Z(t), which differs by a factor of modulus 1
from zeta(1/2 + i*t).
They also write that Hiary and Odlyzko have a book manuscript
in preparation, their reference [HO].
---
I was wondering about band-limited interpolation.
In the Odlyzko–Schönhage algorithm, they don't use the
sine cardinal function because it decays too slowly.
There was a mention somewhere of kernels for
band-limited interpolation. I guess if
the sampling rate is high enough, kernels
that decay faster than sinc(.) can be used to
get good interpolated values.
Would you happen to know of good survey articles
for sine cardinal or other DSP-related
interpolation methods?
regards,
David
David Bernier
[...]
> I was wondering about band-limited interpolation.
> In the Odlyzko–Schönhage algorithm, they don't use the
> sine cardinal function because it decays too slowly.
From reading Xavier Gourdon's 37-page paper about verifying RH up to
10^13 zeros, I got a reference that relates the interpolation
method used by Odlyzko; it's in an unpublished document or paper:
"The 10^20-th zero of the Riemann zeta function and 175 million of its
neighbors", 1992 version.
In section 4.4 of that document, the topic is
band-limited interpolation. The interpolation
formula, choices of parameters and the kernel h(u)
are given in equations 4.4.15 through 4.4.20,
more or less. That's on page number 91 as
printed on the page, or page 94 according
to the Adobe PDF file. Odlyzko writes that
the kernel h(u) was suggested to him by
B.F. Logan .
He might be the same person as this B. F. Logan:
< http://en.wikipedia.org/wiki/Benjamin_F._Logan > .
David Bernier
David Bernier
I understand the Greengard-Rokhlin method a little bit. For a
collection of sources in a unit interval or unit square,
there's a "far-field" approximation [to the 3D potential],
say as a Laurent series in r, r = the distance to the center
of the interval or square. "Far-field" refers to
the fact that the approximation is very good at
distance 10 or more from the center of the
unit interval of unit square.
[ P.S. I don't think the Odlyzko and Schönhage approach
for getting very good approximations to
some sums of rational functions is much easier to
understand/implement . ]
The far-field expansion is very efficient in computational
complexity because even for a million sources in the
unit interval or square, the Laurent series only needs to
be computed once (complexity ~ 1 million), produces
a Laurent series with, say, 40 terms and this Laurent
series can be used for any point outside a ball of radius
10 about the center of the unit interval or square.
If we're at (x,y) in R^2, and the center of the interval/square
is (a,b):
- Let r := sqrt( (a-x)^2 + (b-y)^2 );
- Substitute r in the 40-term Laurent series.
That's the easier part. For a point P where we want the electric
potential, some sources are far enough that they can all be
taken together to produce a far-field expansion,
which will be very good if we're not too far from P.
The remainder of the electric potential comes from sources quite near to P.
If we move P around sufficiently, the sources that can be
put together into a far-field Laurent series expansion change.
With all that, I'm still far from grasping in detail
the Greengard-Rokhlin method.
Above, I was assuming an inverse-square force law which leads
to a potential for a point source of A/r, r being the
distance to the source.
In the complex analysis case, we need to
approximate sum_{ 2<=k < 10^7} a_k/(b_k - z) , where z varies
over the unit circle, I think. The a_k and b_k are fixed constants.
That would be for s = 1/2 + it, t~= 4.0E+15, and where we want
zeta(1/2 + iu), u real and near t.
David Bernier
David Bernier
>>> June 1998, Pages 1-5 .
>>
>>> It's (c) Elsevier Science and not available at ScienceDirect, where
>>> I read the abstract. The formulation of "Trummer's problem" by
>>> someone whose name I forget was followed by the discovery
>>> of the Greengard-Rokhlin method.
>>
>>> ->
>>> <http://comet.lehman.cuny.edu/vpan/pdf/pan150.pdf> (V. Y. Pan).
>>
>>> Interestingly, they cite the 1988 paper by Odlyzko and Schönhage (ref.
>>> no. 9).
>>
>> This document based on a talk of Odlyzko is also worthwhile; he says
>> what the main ideas of the OS algorithm are in the abstract.
>>
>> <http://www.dtc.umn.edu/~odlyzko/doc/arch/zeta.fn.supercomp.pdf>
>> The name of the file is zeta.fn.supercomp.pdf
>>
>> David Bernier
>> --
>> $ gpg --fingerprint david...@videotron.ca
>> pub 2048D/653721FF 2010-09-16
>> Key fingerprint = D85C 4B36 AF9D 6838 CC64 20DF CF37 7BEF 6537 21FF
>> uid David Bernier (Biggy)<david...@videotron.ca>
>
>
> Interesting topic, David. Odlyzko's talk precis you linked
> above mentions that one possible improvement would use the
> Greengard-Roklin algorithm for fast multipole method. Are
> you aware if they explored this further?
For two vectors u and v in R^n or C^n, there's the notion of
the Cauchy matrix C(u, v) which is nxn :
C(u,v)_{i, j} = 1/(u_i - v_j) if i =/= j,
C(u,v)_{i, j} = 0 if i=j.
Then the paper by Pan, M. A. Tabanjeh, Z. Q. Chen, E. I. Landowne and A. Sadikou
mentioned above ,
< http://comet.lehman.cuny.edu/vpan/pdf/pan150.pdf >
has an excellent introduction in the sense that it explains how to
reformulate in linear algebra "Trummer's problem", and one sees the
link between fast Cauchy matrix product by a vector X question and
(a) the fast multi-evaluation of a sum of 'n' rational functions
in the Odlyzko-Schönhage algorithm
[on the one part]
and
(b) the Greengard-Rokhlin method
[ which is more general than the "1/r" potential problem, I think].
I believe part of the difficulties for (a) is when z_i1 and
z_i2 on the unit circle are very close, say within
(2pi/n)/K for K>1000, say.
Perhaps Cauchy matrices is a more concrete way to think of
the problem solved in (a) of the OS algorithm ...
Regards,
David Bernier
consumer
> Chip Eastham wrote:
> > On Sep 18, 6:34 am, David Bernier<david...@videotron.ca> wrote:
> >> David Bernier wrote:
> >>> David Bernier wrote:
> >>> [...]
>
> >>>> (2) The Fast Fourier Transform. This is an essential
> >>>> as it permits fast evaluation of zeta at
> >>>> the grid points through the intermediary
> >>>> of methods very similar to
> >>>> the Fast Multipole Method (or methods)
> >>>> used in potential theory in computational
> >>>> physics.
>
> >>>> Cf.:
> >>>> <http://en.wikipedia.org/wiki/Fast_multipole_method>.
>
> >>>> According to Wikipedia,
> >>>> << The FMM, introduced by Rokhlin and Greengard, has been acclaimed as
> >>>> one of the top ten algorithms of the 20th century[2].>>
>
> >>>> One application in mathematical physics is the study of
> >>>> n-body problems in celestial mechanics, with 'n' large.
> >>>> With 1000 massive bodies, one might think that computing
> >>>> the gravitational field at a random point in the neighborhood
> >>>> of the 1000 bodies could involve summing up to
> >>>> 1000 space-dependent vectors, representing the acceleration
> >>>> due to body #1, #2, ... #1000. Apparently, one can do better
> >>>> and I believe this type of question is addressed by
> >>>> the method(s) that were first given by Rokhlin and Greengard.
>
> >>>> So as next item, I'll put fast evaluation of rational functions
> >>>> ( the FFT allows one to appeal to such methods).
>
> >>>> (3) Super-fast evaluation of some sums of rational functions.
> >>> [...]
>
> >>> In the same area, there's the article:
>
> >>> V. Y. Pan, M. A. Tabanjeh, Z. Q. Chen, E. I. Landowne and A. Sadikou;
> >>> "New transformations of Cauchy matrices and Trummer's problem",
> >>> Computers& Mathematics with Applications, Volume 35, Issue 12,
> >>> June 1998, Pages 1-5 .
>
> >>> It's (c) Elsevier Science and not available at ScienceDirect, where
> >>> I read the abstract. The formulation of "Trummer's problem" by
> >>> someone whose name I forget was followed by the discovery
> >>> of the Greengard-Rokhlin method.
>
> >>> ->
> >>> <http://comet.lehman.cuny.edu/vpan/pdf/pan150.pdf> (V. Y. Pan).
>
> [on the one part]
>
> and
> (b) the Greengard-Rokhlin method
> [ which is more general than the "1/r" potential problem, I think].
>
> I believe part of the difficulties for (a) is when z_i1 and
> z_i2 on the unit circle are very close, say within
> (2pi/n)/K for K>1000, say.
>
> Perhaps Cauchy matrices is a more concrete way to think of
> the problem solved in (a) of the OS algorithm ...
I've been thinking that I might try to understand the
rational function approximate evaluations using a
small test-case.
There's an unpublished paper by Odlyzko,
"The 10^20-th zero of the Riemann zeta function and 175 million of its
neighbors",
revised in 1992. The relevant section for
me is Section 4.3 . The rational function
to be evaluated, f(z), is given
explicitly in formula 4.3.8 , page 81.
Then in each term a_k / (z - b_k) ,
for the desired values of z, I think
both z = omega^h and the b_k (defined in 4.3.6)
lie on the unit circle.
David Bernier