hyp1f1 explodes for large arguments. Recurrence relation unstable.

[cossio]

hyp1f1 explodes for modestly large arguments. For example:

scipy.special.hyp1f1(30, 70, 20j)

returns

(-1748.2243027213672+51.441200608645886j)

which is completely wrong (whenever 1 < a < b and z is purely imaginary, abs(1F1(a,b,z)) <= 1).

This doesn't seem to be insoluble, since the package mpmath (also open sourced here on GitHub) returns what seems to be the correct answer:

mpmath.hyp1f1(30, 70, 20j)
mpc(real='-0.32066129006007188', imag='0.38180277645752586')

(I assume this is correct because Mathematica returns the same result. It would be very unlikely if they were both wrong and returned the same answer).

Moreover, I did a simple test. Simply using the Taylor series expansion gives the correct result. The problem with scipy's implementation seems to be in using a recurrence relation to lower the value of 'a' until it is < 2. The recurrence relation relation used is Abramowitz 13.4.1, but I think this is unstable.

[ev-br]

Mpmath's ways are not always directly applicable since we only have double precision. That being said, hyperbolichypergeometric functions have multiple issues and can definitely be improved. The main problem so far is that we're short on person-power, so patches are welcome!

[ev-br]

Possibly related: gh-3593, gh-3492, gh-2957, gh-2282. Also see the discussion and links in #1561.

[cossio]

It seems that people solve this recurrence unstability using Miller’s
algorithm. I'll let you know if I find the time to implement it. Thanks
anyway.

On Mon, Oct 12, 2015 at 12:13 PM, Evgeni Burovski notifications@github.com
wrote:

Mpmath's ways are not always directly applicable since we only have double
precision. That being said, a hyperbolic functions have multiple issues and
can definitely be improved. The main problem so far is that we're short on
person-power, so patches are welcome!

—
Reply to this email directly or view it on GitHub
#5349 (comment).

[argriffing]

It seems that people solve this recurrence unstability using Miller’s
algorithm.

If it helps, in one of @ev-br's links I've used Miller's algorithm for an edge case of hyp2f1.
#1561 (comment)

[matthew-brett]

It would be great if you had some time to think about this - we have run into this problem too - doing diffusion MRI imaging analysis.

[zhou13]

You may be interested in using cython with gsl as a workaround. The following code can be used as a starter:

# distutils: libraries = gsl

cdef extern from "gsl/gsl_sf_hyperg.h":
    double gsl_sf_hyperg_1F1(double, double, double)
cdef extern from "gsl/gsl_errno.h":
    void *gsl_set_error_handler_off()

def hyp1f1(double a, double b, double c):
    return gsl_sf_hyperg_1F1(a, b, c)

gsl_set_error_handler_off()

From my test, this function is as accurate as mpmath.hyp1f1 with positive parameter combination. It is 10 times faster than the scipy.special.hyp1f1 and about 100 times faster than mpmath.hyp1f1.

[cossio]

GSL's hyperg_1F1 doesn't support complex arguments.

[pv]

GSL is also license-incompatible, so it is not very helpful for the present discussion.

[argriffing]

There is a nice public domain C implementation at https://github.com/fredrik-johansson/arbcmath. But I think it depends on GPL libraries in ways that are not compatible with distributing scipy under a more permissive license.

[argriffing]

Linking a recent mailing list thread: https://mail.scipy.org/pipermail/scipy-dev/2016-January/021164.html.

[samuelstjean]

Well, seems like plain arb itself does it (along with other similarly related functions) [1] and it's lgpl [2], so it could be worth having a look if it's fast enough to replace the old fortran versions from scipy. There is also a set of rudimentary wrappers to python for it [3].

But yeah for anyone wanting to use it for a personal project, the fastest version still remains wrapping the gsl lib.

[1] http://fredrikj.net/arb/acb_hypgeom.html?highlight=hyp#c.acb_hypgeom_m
[2] https://github.com/fredrik-johansson/arb
[3] https://github.com/fredrik-johansson/python-flint

[person142]

Porting over the necessary arbitrary precision machinery is not my preferred way of doing this. (And the algorithms used in arb do not work in double precision for reasonable ranges of the parameters.) A better approach is improving the parabolic cylinder functions to the point that asymptotic expansions like

http://dlmf.nist.gov/13.8

can be used, but that effort is ongoing.

That being said, evaluating the spheroidal wave functions accurately in any reasonable range probably requires arbitrary precision, so it might be necessary to port over an arbitrary precision library at some point anyway. Boost is a good choice, it provides its own implementation but also for linking against gsl.

[samuelstjean]

Such a shame then, I was more thinking of evaluating at arbitrary precision, then recasting to double as it probably uses its own internal representation, but that might end up way too slow for practical use then again. I tried doing just that with mpmath for some time, but that was really too slow for any practical use when one wants a large number of evaluations.

Well anyway, hopefully some numerical workaround for those imprecise regime will turn up someday.

[samuelstjean]

Oh wow, I forgot I had actually looked at it semi seriously, but anyway, well, while the subject came back up elsewhere and I realized maybe people didn't know, the newer version from 0.17 works for more large arguments, but instead of being plain wrong it is now having some precision issues unfortunately, so I'll add this one here

In [13]: hyp1f1(-0.5, 61, -247207.56154023242) - mp.hyp1f1(-0.5, 61, -247207.56154023242)
Out[13]: mpf('-0.10859700032764903')

Assuming mpmath knows what it is doing of course (the gsl version is the same up to 1e-12 in this case, so I'd pin the problem on the scipy fortran version not being precise enough still).