[SciPy-user] hyp2f1 - bad performance for some values of the argument
Anne Archibald
peridot.faceted@gmail....
Fri Jun 1 15:57:15 CDT 2007
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.
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...)
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).
Thanks,
Anne M. Archibald
More information about the SciPy-user
mailing list