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$
  • 2
    $\begingroup$ Can you say which specific parameter values cause underflows? $\endgroup$ – Kirill Aug 8 '14 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 '14 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$ – becko Aug 15 '14 at 15:30
  • 1
    $\begingroup$ If it doesn't have to be exact, you can try a continued fraction expansion. $\endgroup$ – Christian Clason Aug 15 '14 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 '14 at 15:55
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$
  • $\begingroup$ I edited the question to consider large positive values of the parameters. $\endgroup$ – becko Aug 18 '14 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$ – becko Aug 18 '14 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 '14 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$ – becko Aug 19 '14 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 '14 at 15:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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