Sage: Ticket #20127: zetaderiv: numerically unstable
https://trac.sagemath.org/ticket/20127
<p>
The implementation of the derivative of the zeta function <code>zetaderiv</code> is numerically unstable and very slow (for large negative values):
</p>
<pre class="wiki">sage: zetaderiv(1, CIF(-600))
Traceback (most recent call last):
...
NoConvergence: zeta: too much cancellation
</pre><pre class="wiki">sage: timeit('zetaderiv(1, CIF(-500))')
5 loops, best of 3: 3.08 s per loop
</pre><p>
Could either the current implementation be improved or an alternative numerical implementation be used?
</p>
<p>
---
</p>
<p>
This also causes errors when testing relations like
</p>
<pre class="wiki">sage: bool(gamma(k*pi) * zetaderiv(1, k*pi*I)/log(2)^2 == 0)
Traceback (most recent call last):
...
NoConvergence: zeta: too much cancellation
</pre>en-usSagehttps://trac.sagemath.org/chrome/site/logo_sagemath_trac.png
https://trac.sagemath.org/ticket/20127
Trac 1.1.6behacklMon, 29 Feb 2016 06:20:58 GMTcommit, branch set
https://trac.sagemath.org/ticket/20127#comment:1
https://trac.sagemath.org/ticket/20127#comment:1
<ul>
<li><strong>commit</strong>
set to <em>4ddf10c1c65acabc21bae40d80b51735cf26be24</em>
</li>
<li><strong>branch</strong>
set to <em>u/behackl/symbolic/test_relation/noconvergence</em>
</li>
</ul>
<p>
New commits:
</p>
<table class="wiki">
<tr><td><a class="ext-link" href="http://git.sagemath.org/sage.git/commit/?id=4ddf10c1c65acabc21bae40d80b51735cf26be24"><span class="icon"></span>4ddf10c</a></td><td><code>catch NoConvergence error in test_relation</code>
</td></tr></table>
TicketrwsSat, 05 Mar 2016 07:13:30 GMT
https://trac.sagemath.org/ticket/20127#comment:2
https://trac.sagemath.org/ticket/20127#comment:2
<p>
Looks fine, can you please add a doctest?
</p>
TicketbehacklTue, 08 Mar 2016 11:43:59 GMT
https://trac.sagemath.org/ticket/20127#comment:3
https://trac.sagemath.org/ticket/20127#comment:3
<p>
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.
</p>
TicketrwsTue, 08 Mar 2016 14:29:31 GMT
https://trac.sagemath.org/ticket/20127#comment:4
https://trac.sagemath.org/ticket/20127#comment:4
<p>
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?
</p>
TicketbehacklTue, 08 Mar 2016 14:43:39 GMT
https://trac.sagemath.org/ticket/20127#comment:5
https://trac.sagemath.org/ticket/20127#comment:5
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/20127#comment:4" title="Comment 4">rws</a>:
</p>
<blockquote class="citation">
<p>
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?
</p>
</blockquote>
<p>
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.
</p>
<pre class="wiki">sage: bool(gamma(2 + k*pi*I/log(2)) * zetaderiv(1, 1+k*pi*I/log(2))/log(2)^10 == 0)
</pre>
TicketrwsTue, 08 Mar 2016 14:48:52 GMT
https://trac.sagemath.org/ticket/20127#comment:6
https://trac.sagemath.org/ticket/20127#comment:6
<p>
I get 10 times immediately False with develop. I can even do
</p>
<pre class="wiki">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
</pre>
TicketrwsTue, 08 Mar 2016 15:07:39 GMT
https://trac.sagemath.org/ticket/20127#comment:7
https://trac.sagemath.org/ticket/20127#comment:7
<p>
I think the difference to your result is that I had a command <code>k = var('k', domain='complex')</code> instead of <code>k = var('k')</code>. With the latter I can confirm the issue. I thought I needed complex because I used <code>k</code> as integer before that. Now isn't the default complex anyway? It gets more mysterious.
</p>
TicketdkrennTue, 08 Mar 2016 15:24:47 GMT
https://trac.sagemath.org/ticket/20127#comment:8
https://trac.sagemath.org/ticket/20127#comment:8
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/20127#comment:5" title="Comment 5">behackl</a>:
</p>
<blockquote class="citation">
<p>
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.
</p>
<pre class="wiki">sage: bool(gamma(2 + k*pi*I/log(2)) * zetaderiv(1, 1+k*pi*I/log(2))/log(2)^10 == 0)
</pre></blockquote>
<pre class="wiki">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
</pre><p>
returns
</p>
<pre class="wiki">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
</pre><p>
Seems like I am lucky to get that much exceptions ;)
</p>
TicketrwsTue, 08 Mar 2016 15:28:43 GMT
https://trac.sagemath.org/ticket/20127#comment:9
https://trac.sagemath.org/ticket/20127#comment:9
<p>
See <a class="closed ticket" href="https://trac.sagemath.org/ticket/14305" title="defect: Doctest: Immediate simplifications of symbolic powers (closed: fixed)">#14305</a>. Variables may by default be complex in Pynac, but they are not in Maxima, which is (still) used by <code>__nonzero__</code> (see <a class="needs_info ticket" href="https://trac.sagemath.org/ticket/19040" title="enhancement: rewrite Expression.__nonzero__() (needs_info)">#19040</a>).
</p>
TicketrwsTue, 08 Mar 2016 16:06:50 GMT
https://trac.sagemath.org/ticket/20127#comment:10
https://trac.sagemath.org/ticket/20127#comment:10
<p>
The exception happens when executing <code>test_relation()</code>. With <code>domain='complex'</code> <code>test_relation()</code> is not called at all because then the assumption <code>(k, 'complex')</code> has to be taken into account so Maxima is called for proof which, remarkably, is faster here and correct.
</p>
<p>
However what happens with the original case is that in <code>test_relation()</code> random substitutions into the relation are performed to see if the resulting values differ. One substitution that causes the exception is for example <code>{k: -312.363505734545? + 152.010201758413?*I}</code> so you can get the exception with the code <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()</code> reliably.
</p>
TicketbehacklTue, 08 Mar 2016 16:30:30 GMT
https://trac.sagemath.org/ticket/20127#comment:11
https://trac.sagemath.org/ticket/20127#comment:11
<p>
Coincidentally, this solves my original problem. However, I still think that the fact that this
</p>
<pre class="wiki">sage: zetaderiv(1, CIF(-600)).n()
</pre><p>
is both very slow and raises the <code>NoConvergence</code>-error is a problem.
</p>
<p>
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.
</p>
TicketrwsWed, 09 Mar 2016 09:06:52 GMT
https://trac.sagemath.org/ticket/20127#comment:12
https://trac.sagemath.org/ticket/20127#comment:12
<p>
Change of title and component suggested so interested people find it better.
</p>
TicketbehacklWed, 09 Mar 2016 09:14:35 GMTcomponent, description, summary changed
https://trac.sagemath.org/ticket/20127#comment:13
https://trac.sagemath.org/ticket/20127#comment:13
<ul>
<li><strong>component</strong>
changed from <em>symbolics</em> to <em>numerical</em>
</li>
<li><strong>description</strong>
modified (<a href="/ticket/20127?action=diff&version=13">diff</a>)
</li>
<li><strong>summary</strong>
changed from <em>test_relation: uncaught NoConvergence</em> to <em>zetaderiv: numerically unstable</em>
</li>
</ul>
Ticketfredrik.johanssonWed, 09 Mar 2016 12:02:11 GMT
https://trac.sagemath.org/ticket/20127#comment:14
https://trac.sagemath.org/ticket/20127#comment:14
<p>
Arb can compute derivatives of the zeta function without difficulty. E.g. with my own python-flint interface, I can do
</p>
<pre class="wiki">>>> 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)
</pre><p>
which takes 0.1 milliseconds.
</p>
<p>
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 <code>acb_poly_zeta_series</code> and <code>arb_poly_zeta_series</code>.
</p>
<p>
In the left half plane, <code>mpmath.diff(mpmath.zeta, s, n)</code> could also be used instead of <code>mpmath.zeta(s, 1, n)</code>.
</p>
<p>
It's a bit worrying that <code>zetaderiv</code> 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*:
</p>
<pre class="wiki">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
</pre><p>
Are there other Sage functions that treat intervals as carelessly?
</p>
Ticket