## Mathematica and hypergeometric functions

Recently I did some numerical experiments in Mathematica involving the hypergeometric function $_2F_1$. The results were clearly wrong (a positive-definite matrix having negative eigenvalues, for example), so I spent a couple of hours checking the code. Finally, I discovered a bug, but it was not in my program!

This post is not about the hypergeometric function, it is enough to know that it accepts four arguments, usually typed as $_2F_1(a,b;c;x)$. For some parameters $a$, $b$, $c$ it reduces to elementary functions, for example $_2F_1(\frac{1}{2},\frac{1}{2};\frac{3}{2};-x) =$ $\mathrm{arsinh}(\sqrt{x}) / \sqrt{x}$. Of course, Wolfram Research products are aware of that (link). A sligth modification of the parameters of the hypergeometric function does not change the value of the function significantly, it is smooth in all four arguments. But according to Wolfram Alpha, $_2F_1(\frac{1}{2},\frac{1}{2}-\frac{1}{1000};\frac{3}{2}-\frac{1}{1000};-x)$ is not even continuous! Click on an image to see the plot at the Wolfram Alpha site.

There is clearly something wrong with the way Wolfram Alpha (and its kernel, Mathematica) computes $_2F_1$. This issue is related to finite precision calculations: forcing Mathematica to work with higher precision helps a lot, but it does not resolve the problem completely. Numerical instability is nothing spectacular, it is a common issue in approximate computations. Here, however, the error cannot be explained simply by accumulation of numerical errors, the result is just completely wrong: zero instead of somethign close to 0.8. This kind of bug is quite unexpected in such an advanced application.

Let us play around this bug for a while. Wolfram Alpha asked about $_2F_1(\frac{1}{2},\frac{1}{2}-\frac{1}{1000};\frac{3}{2}-\frac{1}{1000};-2)$, answers correctly (link). The numerical approximation (N[…]) is “result: 0.81072238…”, which is fine; the next line, however, reads “number name: zero” (link). Things get even worse if we request five digits only: “decimal approximation: 0.810722, result: 0.0” (link).

I was in fact quite lucky to find this issue: there are not too many parameters for which $_2F_1$ is not computed accurately. In the contour plot of $_2F_1(\frac{1}{2},\frac{1}{2}+p;\frac{3}{2}+q;-2)$ they are all concentrated near the line $p = q$, when $p$ is close to $1$. And this is exactly the case I was interested in.

Nearly a week ago, I submitted this bug to Wolfram Support Center. Unfortunately, apart from an automated confirmation letter (claiming that they typically answer within three business days), I have received no reply.

The above bug might be caused by Mathematica using a wrong hypergeometric identity (perhaps this one?) in numerical evaluation. If this is the case, it will be rather easy to fix. Until then, one can simply use another hypergeometric relation (like this one) explicitly. For an example, see the last two (theoretically identical) plots on the right.

Let me end this post with a few words about what I was trying to do when I discovered this bug. Unfortunately, this will be a little bit technical for those not used to non-local operators and jump processes.

Let $D$ be an open set in $\mathbf{R}^d$. Consider a Brownian motion $X_t$ (with variance $\mathrm{Var} \, X_t = 2t$) starting at a fixed point $x \in D$, and let $\tau_D$ denote the (random) first time when $X_t$ hits the complement of $D$. The mean amount of time that $X_t$ spent in a Borel set $A \subseteq D$ before $\tau_D$ can be expressed as $\int_A G_D(x, y) dy$ for some function $G_D(x, y)$. It can be proved that $G_D(x, y)$ is the Green function for $-\Delta_D$, the Laplace operator in $D$ with zero Dirichlet boundary conditions. In fact, this is probably the easiest way to construct the Green function for a general open set.

Replace in the above paragraph the Brownian motion $X_t$ with the isotropic $\alpha$-stable Lévy process (here $\alpha \in (0, 2)$), and $G_D(x, y)$ will become the Green function for the fractional power of the Laplace operator $(-\Delta)^{\alpha/2}$, with zero exterior condition (this is no longer a local operator, so instead of boundary conditions, one must give exterior condition). This operator is usually denoted as $((-\Delta)^{\alpha/2})_D$, and it is important to note that this is something else than $(-\Delta_D)^{\alpha/2}$. Fractional powers of $-\Delta$ are inverse operators to Riesz potentials, and they are important examples of pseudo-differential operators. It is not my point to give a detailed introduction to probabilistic potential theory and Riesz potential theory herem and there is much literature on these topics, inlcuding classical textbooks by Landkoff, Bliedtner and Hansen, and Blumenthal and Getoor.

If $D$ is a ball, then the formula for the Green function for $-\Delta_D$ is relatively easy to find using Kelvin transform. Suprisingly, similar method works for $((-\Delta)^{\alpha/2})_D$, due to the calculation of M. Riesz in 1938, in this case the formula for the Green function for a ball is:

\displaystyle \begin{aligned} G_D(x, y) & = \frac{\Gamma(\frac{d}{2}) |x - y|^{\alpha - d}}{2^\alpha \pi^{d/2} (\Gamma(\frac{\alpha}{2}))^2} \int_0^{\frac{(1 - x^2)(1 - y^2)}{|x - y|^2}} \frac{s^{\alpha/2 - 1}}{(1 + s)^{d/2}} \, ds \\ & = \frac{\Gamma(\frac{d}{2}) (1 - x^2)^{\alpha / 2} (1 - y^2)^{\alpha / 2}}{2^\alpha \pi^{d/2} \Gamma(\frac{\alpha}{2}) \Gamma(1 + \frac{\alpha}{2}) |x - y|^{d}} \, {_2 F_1}\Bigl(\frac{\alpha}{2}, \, \frac{d}{2}; \, 1 + \frac{\alpha}{2}; \, -\frac{(1 - x^2)(1 - y^2)}{|x - y|^2}\Bigr) . \end{aligned}

While I was working on numerical bounds for the eigenvalues of $((-\Delta)^{\alpha/2})_D$ for $D = (-1, 1)$ (that is, a one-dimensional ball), I needed to compute an array of values of $G_D(x, y)$. Numerical approximation to the integral worked fine, but it was very slow. For that reason, I tried the formula with the hypergeometric function. It worked much faster, but the results were different (and plainly wrong). While I was checking the code, I plotted the graph of $G_D(x, y)$ for $\alpha = 0.99$. The correct picture should be like the one below (prepared using a better hypergeometric identity).

Instead, Mathematica produced the following image.

The erroneous plot can be reproduced in lower resolution using Wolfram Alpha (link).

I will update this post as soon as Wolfram contacts me.

Update, Mar 11, 2011:

I have just received a copy of Mathematica 8 and istalled it on my computer. The first thing I checked was, clearly, whether Wolfram fixed the hypergeometric bug. And they did! A huge surprise to me, as they did not respond to any of my two bug submissions. And the bug is still present on Wolfram Alpha! Does Wolfram Alpha run an older version of Mathematica then?

Update, Mar 23, 2011:

I have just received a kind email from Wolfram Research confirming that the issue has been resolved in Mathematica 8.