Wednesday, January 27, 2010

I've written a very basic mpmath context for computing "natively" with Sage's real and complex numbers. It can be tried out by upgrading mpmath in Sage to the svn trunk, applying this patch this patch to fix the helper code in Sage, and then importing this file.

Unfortunately, there are some issues (besides the fact that some methods are missing, so not everything works yet).

This context doesn't provide variable precision, so the user has to manually allocate sufficient extra precision to compensate for rounding errors.

I first tried to do it, but variable precision is very inconvenient to implement using Sage's way of managing precision. There is no direct way to perform operations with a given target precision (independent of the inputs), and the second best option (which is to perform coercions to the target precision everywhere) is very slow, besides losing accuracy. The only way to implement variable precision in a reasonable way is to bypass Sage's RealNumber and ComplexNumber (at least their public interfaces) and wrap MPFR directly using a precision model similar to what MPFR and mpmath uses, where the precision of the result is always specified and independent of the inputs.

Secondly, there is not that much of a speedup (in the examples above, the speedup is at most about 2x). This is mainly due to the fact that the context uses wrapper classes for RealNumber and ComplexNumber, with all interface and conversion code written in Python. So the overhead is about the same as for the corresponding code in vanilla mpmath (where it accounts for about 50% of the total overhead). The reason RealNumber and ComplexNumber can't be used directly, even in a fixed-precision setting, is that mpmath in many places multiplies numbers by floats (usually exact float literals like 0.5), and Sage always coerces to the number with less precision. This could presumably be fixed by replacing all float and complex constants in mpmath, but I'm not in the mood to do that right now.

There should be a significant performance benefit if direct-from-Cython classes and conversion methods were used. It's also very important to optimize certain core functions; for example, quadrature would benefit greatly from a fast fdot implementation.

All things considered, I'm probably not going to continue much more down this particular path. It's better to write a fully compatible Cython context with classes designed directly for mpmath. Trying to wrap Sage's numbers did however help identify a few problems with the existing interfaces, so I might extend the work on this context a little more just to find more such issues.

Tuesday, January 19, 2010

I'm very grateful to Juan Arias de Reyna who has contributed a module implementing zeta function evaluation using the Riemann-Siegel formula in mpmath. This finally permits rapid evaluation of the Riemann zeta function for arguments with large imaginary part.

To follow tradition on this blog, pictorial proof shall be given. Here is a plot of a segment of the critical line, ζ(1/2+ti) with t between 1010+50 and 1010+55:A complex plot of showing the critical strip, t ranging between 108+40 and 108+45 (note: the y scale is reversed):

Juan Arias de Reyna, who is a professor of mathematics at the University of Seville, has done a thorough job with this code. He has even proved rigorous error bounds for his algorithm (subject to assumptions that the underlying functions in mpmath being well-implemented). The news is that the bounds -- documented in an as-yet unpublished paper -- are valid off the critical line. The code also computes derivatives (up to 4th derivatives), although not as rigorously but still very robustly.

I integrated the code (and added a few optimizations) during the last few days. The zeta function in mpmath now automatically switches between Borwein's algorithm (close to the real line), Euler-Maclaurin summation, and the Riemann-Siegel expansion.

Some example values and timings on my laptop, computing ζ(1/2+10ni) for n up to 12:

The implementation also supports use of the fp context. Double precision unavoidably becomes insufficient as the imaginary part approaches 1015, but it has the advantage of speed in the range where it works:

We have done some comparison with Mathematica, and the mpmath version appears to be about as fast (a bit faster or a bit slower, sometimes substantially faster, depending on circumstance). The most expensive part of the computation occurs in a simple internal function that adds lots of n−s terms. I think for Sage, it will be very easy to switch to a Cython version of this function which should improve speed by a large factor.

But most importantly, Mathematica's Zeta[] is notoriously buggy for large imaginary arguments. As a first example, here Mathematica 7.0 computes two entirely different values at different precisions:

That's all for now. I'm back in school again, so maybe I won't have as much time for programming in the near future. But my schedule is flexible, so we'll see. A new release of mpmath shouldn't be far away.

Update: I should point out that the bugs in Mathematica's Zeta[] only seem to occur in recent versions. Mathematica 5 and older versions do not seem to have these problems.

The code also works with the double precision fp context, which obviously is much faster for large matrices. When computing square roots and logarithms, most of the time is spent on matrix inversion, which can be accelerated by substituting in numpy:

It could be much faster still by doing everything in numpy. Probably in the future fp should be improved to seamlessly use numpy internally for all available matrix operations. Patches are welcome!

Borel summation of divergent hypergeometric series

The hypergeometric series of degree p > q+1 are divergent for |z| > 0. Nevertheless, they define asymptotic expansions for analytic functions which exist in the sense of Borel summation. Previously, mpmath knew only how to evaluate 2F0 (whose Borel sum has a closed form), but now it can evaluate the functions of higher degree as well.

To illustrate, the cosine integral Ci(z) has an asymptotic series representation in terms of 3F0:

This representation is efficient for very large values of |z|, where the argument of the hypergeometric series is small. But with z = 0.5, say, the series does not yield a single digit. With a value of z = 10 or so, the series only gives about 3 digits with an optimal truncation:

Of course, there are ways to compute the cosine integral using convergent series. But if we pretend that we only know about the asymptotic series, then the Borel regularization means that we can evaluate the function to high precision even for small z:

It's instructive to visualize the optimal truncations of an asymptotic series compared to the exact solution. Notice how, for an alternating series like this, the truncations alternate between over- and undershooting:

The catch with this feature? Computing the Borel regularization involves evaluating (nested) numerical integrals with a hypergeometric function in the integrand. Except for special parameters (where degree reductions happen internally -- the Ci series above happens to be such a case), it's not very fast.

Optional consistency with Mathematica

I often try to follow Mathematica's conventions regarding limits, special values and branch cuts of special functions. This simplifies testing (I can directly compare values with Mathematica), and usually Mathematica's conventions seem well-reasoned.

There are exceptions, however. One such exception concerns Mathematica's interpretation of 2F1(a,b,c,z) for b = c a negative integer. Mathematica interprets this as a polynomial, which doesn't make much sense since the denominator of the zero term is zero. It's not even consistent with how Mathematica evaluates this function for a symbolic variable:

In[3]:= Hypergeometric2F1[3,-2,-2,2]

Out[3]= 31

In[4]:= Hypergeometric2F1[3,b,b,2]

Out[4]= -1

Maybe there is a good reason for doing what Mathematica does, but it's at the very least not documented anywhere. I've now changed mpmath to instead interpret the two parameters as eliminating each other and giving a 1F0 function (which is what Mathematica does for a generic value). The Mathematica-compatible value can be recovered by explicitly disabling parameter elimination:

On a related note, I've fixed the Meijer G-function to switch between its z and 1/z forms automatically to follow Mathematica's definition of the function. This introduces discontinuities in the function for certain orders. The user can explicitly choose which form to use (so as to obtain a continuous function) with series=1 or series=2.

Internal reorganization

In a long overdue change, I've moved many of the modules in mpmath into subdirectories. The multiprecision library code is now located in mpmath/libmp; special functions are in mpmath/functions, calculus routines are in mpmath/calculus, and routines for matrices and linear algebra are in mpmath/matrices.

This shouldn't affect any documented interfaces, but it will break external code that directly uses mpmath internal functions. The mpmath interfaces in SymPy and Sage will need some editing for the next version update. The breakage should be straightforward to fix (mostly just a matter of changing imports).

Since the next version of mpmath is definitely going to break some code, I might use the opportunity to do some other cosmetic interface changes as well. Ideally, after the next version, the interface will be stable until and beyond mpmath 1.0 (whenever that happens). But the version number is still at 0.x for a reason -- no compatibility guarantees.

Thursday, January 7, 2010

Holiday season means more free time for coding (yay!). In this commit, I've fixed the remaining major accuracy issues with hypergeometric functions in mpmath. The effect, generally speaking, is that mpmath will now deal much more robustly with large parameters of hypergeometric-type functions.

Previously, you would obtain an accurate value with parameters of magnitude ≈ 10−1 - 101 (say), and often for large parameters as well, but for some functions the accuracy would quickly deteriorate if parameters were increased (or moved closer to singularities). You could still obtain any accuracy by simply increasing the working precision, but you had to manually figure out the amount of extra precision required; the news is that mpmath now automatically gives a value with full accuracy even for large parameters (or screams out loud if it fails to compute an accurate value).

This doesn't mean that you can't trick mpmath into computing a wrong value by choosing sufficiently evil parameters, but it's much harder now than simply plugging in large values. Abnormally large values will now mainly accomplish abnormal slowness (while mpmath implements asymptotic expansions with respect to the argument of hypergeometric functions, evaluation for asymptotically large parameters is a much harder problem as far as I'm aware -- so mpmath in effect accomplishes it by brute force.)

The most trivial change was to change the series summation to aim for control of the relative error instead of the absolute error. This affects alternating series, where large parameters lead to very small function values. Previously, something like the following would happen:

This isn't disastrously bad since the absolute error is controlled (summing this series naively in floating-point arithmetic might give something like 1e+234 due to catastrophic cancellation). But now, you get the relatively accurate answer right away, which is much nicer:

Of course, if the value is too small, the evaluation loop must abort eventually. The default behavior now is to raise an exception if relative accuracy cannot be obtained. The user can force a return value by either permitting a higher internal precision before aborting, or specifying a size 2−zeroprec below which the value is small enough to be considered zero:

Since exceptions are inconvenient, it will be necessary to add some more symbolic tests for exact zeros for certain functions (particularly orthogonal polynomials) where exact zeros arise as important special cases.

Also, cancellation detection in hypercomb is now the default for all functions. This fixes, among other things, a bug in the Meijer G-function reported by a user via email. Before:

Another important change is more correct handling parameters very close to negative integers, particularly those appearing in denominators of series. Previously, unless the integer n was small enough or the precision was set high enough, this situation would yield a bogus value (that of a prematurely truncated series):

Finally, and probably most importantly, the limit evaluation code has been changed to adaptively decrease the size of perturbation until convergence. Under fairly general assumptions, the maximum accurate perturbation at precision p can easily be shown to be 2−(p+C); the problem is that the parameter-dependent constant C isn't generally known. Previously C was just set to a small value (10), and naturally this would break down for some functions when parameters were increased beyond roughly that magnitude.

I briefly considered the idea of estimating C analytically in terms of the parameters, and while I think this can be done rigorously, it seems difficult -- especially to do it tightly (grossly overestimating C would murder performance). The algorithm implemented now is quasi-rigorous, and although there is some slowdown (sometimes by a fair amount), the improved reliability is definitely worth it. A user can also manually supply a perturbation size of their own taste, thereby overriding adaptivity. If the user supplies an intelligent value, this gives both the speed of the old code and full rigor. Probably some intelligent choices for particular functions could be made automatically by mpmath too, to recover the old speed in common cases.

An example of a situation where this becomes necessary is for Legendre functions with certain large parameter combinations. With the old code:

In related news, I've given all hypergeometric-type functions the ability to accept kwargs and pass them on to hypercomb and hypsum. This means that tuning and error control options will be available for a large number of functions (Bessel functions, etc), which could be useful for some users who need to do advanced calculations. I have yet to document these features thoroughly (the interface will also perhaps be tweaked before the next release).