[SciPy-user] Orthogonal polynomial evaluation

[Anne Archibald]

2009/4/13 Pauli Virtanen <pav at iki.fi>:
> Mon, 13 Apr 2009 06:52:16 -0700, Lou Pecora wrote:
>>
>> 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.
>
> Yes, the errors in polynomial evaluation come from a source somewhat
> similar to the reason why Bessel recurrence relations have to be
> evaluated in a certaion direction: loss of precision.
>
> The problem here is that how the routines in Scipy currently work is that
> they compute the polynomial coefficients and produce numpy.poly1d
> objects. All poly1d objects evaluate the polynomials using the *same*
> Horner recurrence scheme, but it is not numerically stable for high-order
> polynomials. The problem is that each of the orthogonal polynomial would
> need to be evaluated using a specialized recurrence relation.
>
> We don't currently have implementations of stable evaluation algorithms.
> So, a complete rewrite of the orthogonal polynomial routines is required.

It's not just a question of polynomial evaluation schemes; the
polynomial coefficients themselves present severe problems. I'd say
that 1, x, ..., x**n is often a very poor basis for expressing
polynomials. It would actually be nice to have a polynomials class
that could use other bases - in fact families of orthogonal
polynomials are natural bases for polynomials.

But let's leave aside poly1d entirely. There are algorithms, and I
thought there was code in scipy, to evaluate orthogonal polynomials,
their zeros, and their derivatives and integrals, in terms of
recurrence relations. These are quite general, fairly efficient, and
reasonably accurate - for small orders. But many orthogonal
polynomials have special-purpose evaluation algorithms: for example
the cos(n*arccos(x)) expression for Chebyshev polynomials. Such
algorithms are usually fast and accurate, but require separate
implementation for each family that has one. I think making it
possible to simply override the evaluation of a particular family of
orthogonal polynomials would be a better solution to the OP's problem.

Of course, once one had efficient accurate evaluation, then one could
think about implementing 1D polynomials represented in other bases...

Anne
> Ticket: http://projects.scipy.org/scipy/ticket/921