4
$\begingroup$

I need a numerically stable way to compute the following ratio:

$$\frac{_{2}F_{1}(a+1,b;c;x)}{_{2}F_{1}(a,b;c;x)}$$

All the parameters are real numbers, with $a< 0$,$\ $ $b,c > 0$ and $0<x<1$.

Right now I am using GSL's implementation of the hypergeometric function, but I keep getting under/over-flows.

Is there a simplification that I can use? Or an aproximation scheme that remains accurate for a wide range of the parameters (say $-1000 < a,b,c < 1000$)?

As a particular case, I am also interested in positive values of the parameters, in the range $0<a,b,c<1000$. What approximation can I use in that case?

$\endgroup$
7
  • 2
    $\begingroup$ Can you say which specific parameter values cause underflows? $\endgroup$
    – Kirill
    Aug 8, 2014 at 16:23
  • 1
    $\begingroup$ By the way, what error code does gsl give you? Is it definitely GSL_EUNDRFLW? Please post your code, otherwise it is hard to guess whether this is a problem with gsl. $\endgroup$
    – Kirill
    Aug 8, 2014 at 16:54
  • $\begingroup$ @Kirill I get under/over-flows for large values of $a,b,c$ ($\approx 100$). I don't think the problem is with GSL's implementation. The problem is that to compute the ratio, I am computing the numerator and the denominator, which are very small or very large, and thus their ratio is not very accurate. What I need is a way to directly compute the ratio. $\endgroup$
    – a06e
    Aug 15, 2014 at 15:30
  • 1
    $\begingroup$ If it doesn't have to be exact, you can try a continued fraction expansion. $\endgroup$ Aug 15, 2014 at 15:40
  • $\begingroup$ @becko For $(a,b,c,x)=(-100.1, 100.2, 100.3, 0.5)$ gsl tells me "the requested feature is not (yet) implemented". I understand what you want, but I want to point out that gsl might not be giving you an underflow. $\endgroup$
    – Kirill
    Aug 15, 2014 at 15:55

1 Answer 1

4
$\begingroup$

I think the main issue that makes this tricky is that trying to evaluate the power series for the hypergeometric function is numerically unstable due to catastrophic cancellation for one negative parameter and function value $\ll1$. Otherwise, having real parameters, and a real argument in $[0,1]$ should be fine.

If you look at the DLMF entry for hypergeometric functions, you will find a relation between contiguous functions, $F(a,b;c;z)$ and $F(a\pm1,b;c;z)$: $$ f_a = \frac{F(a+1,b;c;z)}{F(a,b;c;z)} $$ $$ (c-a) \frac1{f_{a-1}} + (2a-c+(b-a)z) + a(z-1)f_a = 0. $$

So for large negative $a$, you can find $\bar a=a-\lfloor a\rfloor\geq0$ and then relate $f_a$ to $f_{\bar a}$, which are separated in the above recurrence relation by integer distance $-\lfloor a\rfloor$. This algorithm takes time linear in $|a|$, so it would work well for $|a|<1000$, but poorly for $|a|$ much larger than $10^6$, say.

Since $\bar a$ is positive, evaluating $f_{\bar a}$ is a lot simpler since it involves hypergeometric function of real variable with all positive parameters (so no cancellation). In particular, you can do with power series (or GSL for that matter).

Here is some python code that I wrote to test this idea. (It doesn't work for some values of $a, c$ that cause division by zero, but that's straightforward to fix; you can use limits for example.)

from mpmath import *

def hsum(a, b, c, x):
  ans, ansAbs, t, i = mpf(1), mpf(1), mpf(1), 0
  while abs(t) / abs(ans) > mp.eps():
    t *= (a + i) * (b + i) / ((c + i) * (i+1)) * x
    i += 1
    ans += t
    ansAbs += abs(t)
  print("hsum(a=%f): %f (%f)" % (a, ans, ansAbs / abs(ans)))
  return ans

def h3(a, b, c, x):
  a, b, c, x = mpf(a), mpf(b), mpf(c), mpf(x)

  ans0 = hyp2f1(a+1, b, c, x) / hyp2f1(a, b, c, x)
  if a >= 0:
    return ans0

  def f(a):
    if a >= 0: return hsum(a+1, b, c, x) / hsum(a, b, c, x)
    return - (c-a-1) / ((2*(a+1)-c+(b-a-1)*x) + (a+1)*(x-1)*f(a+1))

  def f2(a):
    n = int(floor(a))
    a0 = a - n
    ans = hsum(a0+1, b, c, x) / hsum(a0, b, c, x)
    def alpha(a): return -(c-a-1)/((a+1) * (x-1))
    def beta(a): return (2*(a+1)-c+(b-a-1)*x) / ((a+1)*(x-1))
    aa = a0
    for i in range(0, -n):
      aa -= 1
      ans = alpha(aa) / (beta(aa) + ans)
    return ans

  ans1 = f2(a)

  if (abs(a) <= 200):
    assert(abs(f2(a)/f(a)-1) < eps() * 1000)

  print("a %s: %s" % (a, float(abs(ans1/ans0-1))))
  return ans1

P.S. You said in a comment "I am computing the numerator and the denominator, which are very small or very large, and thus their ratio is not very accurate". This isn't really right. The numerical stability of evaluating the expression $a/b$ can be found by computing the condition number of the function $(a,b)\mapsto a/b$. It is $$ \left|\frac{d\log(a/b)}{d\log a}\right| = 1, \qquad \left|\frac{d\log(a/b)}{d\log b}\right| = 1. $$ Thus the computation of $a/b$ should be well-conditioned, regardless of their values (assuming no overflow/underflow, which is not what's going on for $|a|,|b|,|c|\approx 100$ anyway). You might be thinking of $a-b$, which does have a singular condition number at $a\approx b$: $$ \left|\frac{d\log(a-b)}{d\log a}\right| = \frac{|a|}{|a-b|}. $$ This doesn't happen to $a/b$. The numerical problems here are of a different sort.

$\endgroup$
5
  • $\begingroup$ I edited the question to consider large positive values of the parameters. $\endgroup$
    – a06e
    Aug 18, 2014 at 19:48
  • $\begingroup$ In this case, I fear this recursion will be too slow. I need to evaluate the hypergeometric many times, with large values of $a$ (say ~600). $\endgroup$
    – a06e
    Aug 18, 2014 at 20:08
  • $\begingroup$ @becko I checked and I think the same recurrence should work fine, unless you have some other constraints. There is no recursion here, it's just a loop. $\endgroup$
    – Kirill
    Aug 18, 2014 at 21:48
  • $\begingroup$ When I said that computing the numerator and the denominator wasn't a good approach, I meant that sometimes, the numerator and the denominator are not representable in floating-point (resulting in under/over-flow), even though their ratio is within floating-point range. This happens even though division is well-conditioned. $\endgroup$
    – a06e
    Aug 19, 2014 at 14:49
  • $\begingroup$ @becko $F(-100.1,100.2;100.3;0.5)=9.55\times 10^{-31}$ is not an underflow, so that's why I thought that was wrong. $\endgroup$
    – Kirill
    Aug 19, 2014 at 15:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.