[SciPy-user] Orthogonal polynomial evaluation

[Lou Pecora]

I recall that stable routines to calculate Bessel functions use descending recurrence relations (descending in indices), rather than ascending recurrence which is unstable.  Only trick is that the values have to be renormalized at the end (or as you go to keep the numbers in floating point range, I think) since descending recurrence does not set the scale initially.  I wonder if this is the case for other recurrence relations.  That is, if one recurrence is unstable (e.g. ascending), then the opposite (descending) will be stable.  Perhaps the scipy routines can be easily reworked, if so.  Just a thought.

-- Lou Pecora,   my views are my own.

--- On Mon, 4/13/09, Anne Archibald <peridot.faceted at gmail.com> wrote:

From: Anne Archibald <peridot.faceted at gmail.com>

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
_______________________________________________
SciPy-user mailing list
SciPy-user at scipy.org
http://mail.scipy.org/mailman/listinfo/scipy-user



      
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20090413/81d35b60/attachment.html>