Opened 3 years ago

Last modified 3 years ago

#20127 new defect

zetaderiv: numerically unstable

Reported by: behackl Owned by:
Priority: major Milestone: sage-7.1
Component: numerical Keywords:
Cc: cheuberg, dkrenn, rws Merged in:
Authors: Benjamin Hackl Reviewers:
Report Upstream: N/A Work issues:
Branch: u/behackl/symbolic/test_relation/noconvergence (Commits) Commit: 4ddf10c1c65acabc21bae40d80b51735cf26be24
Dependencies: Stopgaps:

Description (last modified by 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

Change History (14)

comment:1 Changed 3 years ago by behackl

  • Branch set to u/behackl/symbolic/test_relation/noconvergence
  • Commit set to 4ddf10c1c65acabc21bae40d80b51735cf26be24

New commits:

4ddf10ccatch NoConvergence error in test_relation

comment:2 Changed 3 years ago by rws

Looks fine, can you please add a doctest?

comment:3 Changed 3 years ago by behackl

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.

comment:4 follow-up: Changed 3 years ago by rws

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?

comment:5 in reply to: ↑ 4 ; follow-up: Changed 3 years ago by behackl

Replying to rws:

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)

comment:6 Changed 3 years ago by rws

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

comment:7 Changed 3 years ago by rws

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.

comment:8 in reply to: ↑ 5 Changed 3 years ago by dkrenn

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 ;)

comment:9 Changed 3 years ago by rws

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

comment:10 Changed 3 years ago by rws

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.

comment:11 Changed 3 years ago by behackl

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.

comment:12 Changed 3 years ago by rws

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

comment:13 Changed 3 years ago by behackl

  • Component changed from symbolics to numerical
  • Description modified (diff)
  • Summary changed from test_relation: uncaught NoConvergence to zetaderiv: numerically unstable

comment:14 Changed 3 years ago by fredrik.johansson

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?

Note: See TracTickets for help on using tickets.