[SciPy-user] Orthogonal polynomial evaluation
Anne Archibald
peridot.faceted@gmail....
Mon Apr 13 03:59:39 CDT 2009
2009/4/13 Andrew York <Andrew.G.York+scipy@gmail.com>:
> I am running a simulation with Scipy that needs to numerically
> evaluate Jacobi polynomials. I have been using the function
> scipy.special.orthogonal.jacobi() heavily. Evaluations of this
> function are >90% of my simulation running time. Much more troubling,
> I begin to see unphysical results when I need to evaluate high-order
> Jacobi polynomials, and I suspect I am running into the problem warned
> against here:
>
> http://docs.scipy.org/doc/scipy/reference/special.html
>
> "Warning Evaluating large-order polynomials using these functions can
> be numerically unstable."
Yes, unfortunately all scipy's orthogonal polynomials are computed
using generic recurrence relations. This allows a quick flexible
implementation, but as you and I both discovered it suffers from
numerical instability.
> Another Scipy user seems to have run into a similar problem:
>
> http://thread.gmane.org/gmane.comp.python.scientific.user/11206
>
> So, it seems the correct course of action is clear: I should identify
> a stable algorithm to numerically evaluate Jacobi polynomials, I
> should implement this algorithm in something fast like C, and then
> call it from within python. I found some promising information on a
> possible algorithm:
>
> http://portal.acm.org/citation.cfm?id=363393
> http://portal.acm.org/citation.cfm?id=362686.362700
> https://people.sc.fsu.edu/~burkardt/f77_src/toms332/toms332.html
>
> But I do not speak FORTRAN. I am also new to python and not yet
> familiar with how to integrate a fast low-level language with python,
> nor am I especially proficient in any low-level language. I suspect if
> numerical evaluation of these functions was faster or more robust,
> other scipy users would also benefit. Should I attempt this myself (in
> C with weave?)? If I do, would others be interested in the results?
>
> Is this an easy task for anyone else? I would greatly appreciate any
> help or guidance.
If you are comfortable with python but not FORTRAN or
python-to-other-language linking, then I recommend using cython to
implement compiled extensions. It can be very fast, and the path to a
working implementation is all modest steps:
* write a python point-at-a-time implementation based on the articles
* write a python unit testing module (makes the rest much easier!)
* compile the point-at-a-time implementation with cython
* add enough type annotations that you can see that the generated C
code contains no python calls
* write a wrapper function in cython that calls the point-at-a-time
function on each element of an array
Ideally, I'd like to see the special functions module modified so that
the orthogonal polynomials could use better implementations where
available, falling back to the recurrence relation. But it's kind of a
maze.
Good luck,
Anne
More information about the SciPy-user
mailing list