[SciPy-user] hyp2f1 - bad performance for some values of the argument
David M. Cooke
cookedm@physics.mcmaster...
Sun Jun 3 14:38:48 CDT 2007
On Fri, Jun 01, 2007 at 04:57:15PM -0400, Anne Archibald wrote:
> On 30/05/07, Anne Archibald <peridot.faceted@gmail.com> wrote:
>
> > So this is not generally useful, but I have cured my problem through
> > an application of one of the "quadratic transformations" (equation
> > 15.3.32 in Abramowitz & Stegun - which is online!). It gives me
> > accuracies on the order of one part in 10^13, not as good as I was
> > hoping but better than the 1 in 10^8 I was getting from the averaging
> > shortcut. Good enough to get positive definite matrices out of it
> > anyway.
>
> Perhaps I spoke too soon.
>
> It appears that scipy evaluates hyp2f1 very slowly for certain values
> of the argument:
>
> In [92]: x=0.1; n=1000; m=100;
> timeit.Timer("hyp2f1(1/2.,1.,3.2,x)",setup="from scipy.special import
> hyp2f1; from numpy import ones;
> x=%f*ones(%d)"%(x,m)).timeit(n)/float(n*m)
> Out[92]: 2.2619605064392088e-06
>
> In [93]: x=0.999; n=1000; m=100;
> timeit.Timer("hyp2f1(1/2.,1.,3.2,x)",setup="from scipy.special import
> hyp2f1; from numpy import ones;
> x=%f*ones(%d)"%(x,m)).timeit(n)/float(n*m)
> Out[93]: 0.00085240376949310298
>
> (uh, no I don't write real code like this; those numbers are seconds
> for a single function evaluation. I use the vectorized version because
> for small x the function-call overhead swamps the evaluation time.)
>
> Presumably it works by transforming the function until |x|<1 and then
> using the series. But the series converges very slowly for |x| close
> to 1. It is possible to keep transforming the function until |x|<1/2,
> where the series converges at a respectable speed.
For c < a+b, |x| < 1, it first tries the power series. If the
accumulated error is too large, it uses the recurrence from 15.2.27 to
move the value of c. (I don't know how effective that is).
For c > a+b, |x| < 1, it again tries the power series first, then uses
15.3.6 to transform x to 1-x, and tries the power series again.
> Of course these transformations suffer from the same difficulties as
> the ones used to get |x|<1, namely, they are singular for certain
> values of the arguments; while there are special-case formulas in A&S
> for these situations, they do not handle points *near* the
> singularities (although conceivably a derivative formula could be used
> there...)
Except that you need derivatives with respect to the parameters, for
which I don't know of any nice expressions except for doing it on the
power series
> Any suggestions on how to efficiently evaluate hyp2f1(0.5,1,a+2,x) for
> x near 1? It turns out to be the limiting factor for performance of my
> program (dwarfing the linear algebra on 300 by 300 matrices).
AMS 15.3.6? That turns hyp2f1(0.5, 1, a+2, x) (for x < 1) into
GAMMA(a+2)*GAMMA(a+0.5)/(GAMMA(a+1.5)*GAMMA(a+1)) * hyp2f1(0.5, 1,
0.5-a, 1-x)
+ (1-x)^(a+0.5)*GAMMA(a+2)*GAMMA(-0.5-a)/sqrt(pi) * hyp1f0(a+1,1-x)
hyp1f0(a+1,1-x) = x^(-a-1), so the last term is
(1/x-1)^a * sqrt(1-x) * GAMMA(a+2)*GAMMA(-0.5-a)/sqrt(pi)
Further Maple manipulation suggests that total is
(2*a+2)/(2*a+1) * hyp2f1(1/2,1,1/2-a, 1-x)
- sqrt(Pi)*GAMMA(a+2)/GAMMA(a+3/2) * sec(Pi*a) * x^(-a-1)*(1-x)^(a+1/2)
which is probably as good as it's going to get. There is a (removable)
singularity at a ~ (2*n+1)/2; I'm not sure how to avoid that except by
clever rewriting of the power series for hyp2f1(1/2,1,1/2-a,1-x).
--
|>|\/|<
/--------------------------------------------------------------------------\
|David M. Cooke http://arbutus.physics.mcmaster.ca/dmc/
|cookedm@physics.mcmaster.ca
More information about the SciPy-user
mailing list