Opened 5 years ago
Last modified 5 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 )
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 5 years ago by
- Branch set to u/behackl/symbolic/test_relation/noconvergence
- Commit set to 4ddf10c1c65acabc21bae40d80b51735cf26be24
comment:2 Changed 5 years ago by
Looks fine, can you please add a doctest?
comment:3 Changed 5 years ago by
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: ↓ 5 Changed 5 years ago by
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: ↓ 8 Changed 5 years ago by
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 5 years ago by
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 5 years ago by
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 5 years ago by
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 5 years ago by
comment:10 Changed 5 years ago by
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 5 years ago by
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 5 years ago by
Change of title and component suggested so interested people find it better.
comment:13 Changed 5 years ago by
- Component changed from symbolics to numerical
- Description modified (diff)
- Summary changed from test_relation: uncaught NoConvergence to zetaderiv: numerically unstable
comment:14 Changed 5 years ago by
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?
New commits:
catch NoConvergence error in test_relation