Errors in complex hyp2f1 values (Trac #1034)

[scipy-gitbot]

Original ticket http://projects.scipy.org/scipy/ticket/1034 on 2009-10-25 by @pv, assigned to @pv.

Complex-valued stuff from gh-1286:

>>> mpmath.hyp2f1(1,1,4,3+4j) 
mpc(real='0.49234384000963544', imag='0.60513406166123973')
>>> special.hyp2f1(1,1,4,3+4j)
(0.65925460474847919+1.3277887264473263j) 

[scipy-gitbot]

Milestone changed to Unscheduled by @rgommers on 2011-06-13

[rgommers]

There's still some issue for real values as well:

In [4]: special.hyp2f1(0.5, 1 - 270.5, 1.5, 0.999**2)
Out[4]: nan

Here's a plot showing the region of inputs for which nans are returned:

import numpy as np
from scipy import special
import matplotlib.pyplot as plt


x = np.r_[np.linspace(-1, -0.99, num=1e2),
          np.linspace(-0.99, 0.99, num=1e2),
          np.linspace(0.99, 1, num=1e2)]
c = np.linspace(0, 2000, num=600)
xx, cc = np.meshgrid(x, c)
res = np.ones((x.shape[0], c.shape[0]))
for ii, xi in enumerate(x):
    res[ii, :] = special.hyp2f1(0.5, 1 - c / 2.0, 1.5, xi**2)

ix = np.isnan(res)
res[ix] = 1
res[~ix] = 0
plt.imshow(res, cmap=plt.cm.bone_r)

plt.show()

This was found by investigating gh-1285 in more detail.

[gertingold]

In the original complex example, HYGFZ is making use of http://dlmf.nist.gov/15.8#E8 together with http://dlmf.nist.gov/5.5#E4. If c-a equals an integer, divergencies arise in http://dlmf.nist.gov/15.8#E8. Instead of handling this divergence as described after http://dlmf.nist.gov/15.8#E9, c is shifted by 1e-8 (cf. lines 6126 and 6323 in specfun.f). This trick does not really solve the problem as can be seen from the following plot:

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import hyp2f1

a, b, c = 1, 1, 4
z = 3+4j

delta_c = np.logspace(-15, -5, 500)
y = hyp2f1(a, b, c-delta_c, z)
plt.plot(delta_c, y)
plt.xscale("log")
plt.ylim([0, 2])
plt.show()

[pv]

The HYGFZ implementation in SPECFUN may be lacking also in other respects, I'm a bit worried about its behavior for large a, b, or c. The systematic mpmath comparison test scipy/special/tests/test_mpmath.py:TestSystematic.test_hyp2f1_complex displays a number of problems, but I haven't yet got around to looking into them in detail.

It may be faster to rewrite the routine from scratch (in C or C++) following some systematic approach to 2F1 evaluation, than to try to address the problems in the F77 code. I believe this problem has been extensively studied in the literature.

I'm not aware of license-compatible 2F1 implementations that we could use.

[ev-br]

It's been studied indeed: Somebody even wrote an MSc thesis on exactly this. http://people.maths.ox.ac.uk/porterm/research/pearson_final.pdf

[argriffing]

Somebody even wrote an MSc thesis on exactly this.

I've implemented a tiny part of this, and it seems to work well enough to solve #3479. A newer publication http://arxiv.org/pdf/1407.7786v1.pdf has more details, but it has some typos in the formulas that could be annoying if you are just directly transcribing it -- there are fencepost errors in the indexing and at least one of the yN should be y0. The arxiv links to some code http://datashare.is.ed.ac.uk/handle/10283/607 which is in a repository marked as http://opendatacommons.org/licenses/by/. I'm not sure whether that code is actually licensed under that license, and if so, I'm not sure whether that license is compatible with scipy.

from __future__ import print_function, division

import numpy as np
import scipy.linalg
from scipy.special import hyp2f1

def recurrence_00p(a, b, c, z, n):
    denom = (c-a+n)*(c-b+n)*z
    an = (c+n)*(c+n-1)*(z-1) / denom
    bn = (c+n)*(c+n-1-(2*(c+n)-a-b-1)*z) / denom
    return an, bn

def solve_recurrence(a, b, c, z, N, recurrence):
    # This could be improved in a couple of ways.
    # First, this uses an upper triangular matrix with
    # small bandwidth that doesn't need to be constructed
    # and solved naively. Second, it could be solved using
    # Olver's Algorithm instead of Miller's Algorithm
    # to avoid using arbitary N.
    M = np.identity(N+1)
    for i in range(N-1):
        an, bn = recurrence(a, b, c, z, i+1)
        M[i, i+0] = an
        M[i, i+1] = bn
        M[i, i+2] = 1
    v = np.zeros(N+1)
    v[-2] = 1
    return scipy.linalg.solve(M, v)

def main():
    a = 10
    b = 5
    c = -300.5
    k = 300 # We want a recursion reducing |c| by this amount.
    z = 0.5
    N = 400 # An arbitrary number a bit bigger than k.
    y = solve_recurrence(a, b, c, z, N, recurrence_00p)
    fk = hyp2f1(a, b, c+k, z) # Base case.
    f0_approx = fk * (y[0] / y[k])
    print(f0_approx)

main()

-3.85202708152e+32

[argriffing]

There is also http://arxiv.org/abs/0708.0116. The code is very safe in http://www.cpc.cs.qub.ac.uk/ where Elsevier has top men working on it right now.

[matthew-brett]

CPC - amazing website! Code apparently has license http://www.cpc.cs.qub.ac.uk/licence/licence.html. C++ and Fortran 90 code. I guess we could ask the authors to relicense.

[pv]

The real-valued cephes hyp2f1 actually knows it's returning garbage in this case, because it also calculates error estimates:

>>> from scipy.special import hyp2f1, errprint
>>> errprint(1)
0
>>> hyp2f1(10, 5, -300.5, 0.5)
__console__:1: SpecialFunctionWarning: scipy.special/hyp2f1: loss of precision
-6.5184944973518031e+156

One probably should add something like if (err > 1e-7) return NPY_NAN; in order to avoid returning garbage.

[argriffing]

This is a bit off-topic, but here's a recent 2F1 blog post by the mpmath author who is now developing a C arbitrary precision library http://fredrikj.net/blog/2015/10/the-2f1-bites-the-dust/.

[h-vetinari]

Still relevant as of 1.10, unfortunately not fixed by @steppi's work in the 1.8 cycle.

[steppi]

Still relevant as of 1.10, unfortunately not fixed by @steppi's work in the 1.8 cycle.

I’m going to start working on this again shortly. There’s still a couple of other cases where the Fortran needs to be replaced and then I’ll implement Miller’s algorithm or Olver’s algorithm to handle large a, b, c.