zetaderiv: numerically unstable

[behackl]

The implementation of the derivative of the zeta function zetaderiv is numerically unstable and very slow (for large negative values):

sage: zetaderiv(1, CIF(-600))
Traceback (most recent call last):
...
NoConvergence: zeta: too much cancellation

sage: timeit('zetaderiv(1, CIF(-500))')
5 loops, best of 3: 3.08 s per loop

Could either the current implementation be improved or an alternative numerical implementation be used?

---

This also causes errors when testing relations like

sage: bool(gamma(k*pi) * zetaderiv(1, k*pi*I)/log(2)^2 == 0)
Traceback (most recent call last):
...
NoConvergence: zeta: too much cancellation

CC: @cheuberg @dkrenn @rwst

Component: numerical

Author: Benjamin Hackl

Branch/Commit: u/behackl/symbolic/test_relation/noconvergence @ 4ddf10c

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

[behackl]

Branch: u/behackl/symbolic/test_relation/noconvergence

[behackl]

Commit: 4ddf10c

[behackl]

New commits:

4ddf10c

catch NoConvergence error in test_relation

[rwst]

comment:2

Looks fine, can you please add a doctest?

[behackl]

comment:3

In principle I could, yes. However, the simplest (reliable) doctest I can construct takes about 25 sec., and that's not really ideal. I'd really like to understand why this takes so much time.

[rwst]

comment:4

Which code would that be? Have you tried callgrind? I tried to confirm the random behaviour of the above but the first 16 tries all gave False. Do you have a reliable example?

[behackl]

comment:5

Replying to @rwst:

Which code would that be? Have you tried callgrind? I tried to confirm the random behaviour of the above but the first 16 tries all gave False. Do you have a reliable example?

Could you try it with the following expression? This is the original one which brought me to this error. If it works (which is rarely; maybe one out of 20 times on my laptop here), it takes several seconds.

sage: bool(gamma(2 + k*pi*I/log(2)) * zetaderiv(1, 1+k*pi*I/log(2))/log(2)^10 == 0)

[rwst]

comment:6

I get 10 times immediately False with develop. I can even do

sage: %timeit assert not bool(gamma(2 + k*pi*I/log(2)) * zetaderiv(1,1+k*pi*I/log(2))/log(2)^10 == 0)
10 loops, best of 3: 45.3 ms per loop

[rwst]

comment:7

I think the difference to your result is that I had a command k = var('k', domain='complex') instead of k = var('k'). With the latter I can confirm the issue. I thought I needed complex because I used k as integer before that. Now isn't the default complex anyway? It gets more mysterious.

[dkrenn]

comment:8

Replying to @behackl:

Could you try it with the following expression? This is the original one which brought me to this error. If it works (which is rarely; maybe one out of 20 times on my laptop here), it takes several seconds.

sage: bool(gamma(2 + k*pi*I/log(2)) * zetaderiv(1, 1+k*pi*I/log(2))/log(2)^10 == 0)

from datetime import datetime
var('k')

for _ in range(10):
    tic = datetime.now()
    try:
        result = bool(gamma(2 + k*pi*I/log(2)) * zetaderiv(1, 1+k*pi*I/log(2))/log(2)^10 == 0)
    except Exception as e:
        print datetime.now()-tic, 'exception:', e
    else:
        print datetime.now()-tic, 'successful:', result

returns

0:00:11.445669 exception: zeta: too much cancellation
0:00:08.250215 exception: zeta: too much cancellation
0:00:03.706793 exception: zeta: too much cancellation
0:00:02.428823 successful: False
0:00:00.757367 successful: False
0:00:16.853168 exception: zeta: too much cancellation
0:00:00.228982 successful: False
0:00:03.824248 successful: False
0:00:37.852729 exception: zeta: too much cancellation
0:00:15.264686 exception: zeta: too much cancellation

Seems like I am lucky to get that much exceptions ;)

[rwst]

comment:9

See #14305. Variables may by default be complex in Pynac, but they are not in Maxima, which is (still) used by __nonzero__ (see #19040).

[rwst]

comment:10

The exception happens when executing test_relation(). With domain='complex' test_relation() is not called at all because then the assumption (k, 'complex') has to be taken into account so Maxima is called for proof which, remarkably, is faster here and correct.

However what happens with the original case is that in test_relation() random substitutions into the relation are performed to see if the resulting values differ. One substitution that causes the exception is for example {k: -312.363505734545? + 152.010201758413?*I} so you can get the exception with the code (gamma(2 + k*pi*I/log(2)) * zetaderiv(1,1+k*pi*I/log(2))/log(2)^10 == 0).subs({k: CIF(-312.363505734545 + 152.010201758413*I)}).n() reliably.

[behackl]

comment:11

Coincidentally, this solves my original problem. However, I still think that the fact that this

sage: zetaderiv(1, CIF(-600)).n()

is both very slow and raises the NoConvergence-error is a problem.

Is it possible to use a numerically more stable version of the derivative of the zeta function? Maybe someone knowing more about this numerical stuff could help with this.

[rwst]

comment:12

Change of title and component suggested so interested people find it better.

[fredrik-johansson]

comment:14

Arb can compute derivatives of the zeta function without difficulty. E.g. with my own python-flint interface, I can do

>>> ctx.cap = 10
>>> acb_series([-600,1]).zeta()
([7.82232679749e+928 +/- 8.22e+916])*x + ([-3.56689160315e+929 +/- 4.12e+917])*x^2 + ([7.8112800125e+929 +/- 5.34e+918])*x^3 + ([-1.08969439943e+930 +/- 6.02e+918])*x^4 + ([1.07928682824e+930 +/- 9.55e+918])*x^5 + ([-7.957487571e+929 +/- 4.23e+919])*x^6 + ([4.390792240e+929 +/- 6.15e+919])*x^7 + ([-1.700332217e+929 +/- 6.11e+919])*x^8 + ([3.04336993e+928 +/- 6.00e+919])*x^9 + O(x^10)

which takes 0.1 milliseconds.

This would be easier to wrap with a Sage wrapper for Arb power series in place, but it should not be too hard to do directly either: see acb_poly_zeta_series and arb_poly_zeta_series.

In the left half plane, mpmath.diff(mpmath.zeta, s, n) could also be used instead of mpmath.zeta(s, 1, n).

It's a bit worrying that zetaderiv currently accepts CIF input and outputs a nonrigorous CIF without warning, by going through a plain numerical computation. It is easy to produce examples where the output is plain wrong:

sage: q = CIF("2.46316186945432128587439505331", "23.2983204927628579020109616266")
sage: zetaderiv(1,q)
-3.8826886735960628?e-17 - 7.4180200774526877?e-17*I
sage: q = ComplexIntervalField(128)("2.46316186945432128587439505331", "23.2983204927628579020109616266")
sage: zetaderiv(1,q)
2.809208149461043895562049836274827424167?e-31 + 4.678424144202694674839595616043132108038?e-32*I

Are there other Sage functions that treat intervals as carelessly?