more accurate representation of polynomials

[pv]

numpy.poly1d represents polynomials by storing the polynomial coefficients. I believe that for high-order polynomials, this is not always the optimal way to represent them.

Consider one of the orthogonal polynomial functions in scipy.special:

>>> import scipy, scipy.special
>>> for n in xrange(0,70,10): print n, scipy.special.chebyt(n)(scipy.cos(z)) - scipy.cos(n*z)
0 0.0
10 3.80806497446e-14
20 1.94978477808e-10
30 2.65364360286e-06
40 -0.0237953834481
50 -114.714051465
60 272370.465462

It can be seen that as the order of the polynomial increases, the accuracy of evaluation decreases very rapidly, and for n > 40 the result is basically numerical noise.

The scipy.special.chebyt function generates a poly1d object by passing the roots of the polynomial to poly1d. In poly1d.init, polynomial coefficients are calculated from the roots. For high-order Chebyshev polynomials, the coefficients are large and of varying sign which leads to loss of precision when evaluating the polynomial.

It might be useful if there was a poly1d-compatible object that would retain a more accurate representation of the polynomial. (Using the roots, or a Horner scheme?) This would come useful when more information than the coefficients is initially available.

See also http://scipy.org/scipy/scipy/ticket/581

[charris]

Polynomials represented as power series are inherently ill conditioned. If you want accuracy for high degree polynomials use Chebyshev polynomial expansions. Don't convert the Chebyshev representation to a polynomial, you will lose all the numerical advantages if you do so, instead evaluate the expansion using backward recursion. Backward recursion, of which Horner's method is a special case, is useful for evaluting polynomial series whose terms are defined recursively.

What you need is a complete Chebyshev package. I posted my own Chebyshev package on the Numpy developer's list sometime back. If you search the list you may find it. There were actually two postings, you will want the last one. Perhaps that package could be adapted to use parts of the SciPy? Chebyshev package, at the moment it is written in pure Python/Numpy. In any case, such a package shouldn't be part of Numpy. Polynomials really don't belong in Numpy either.

[pv]

This is probably good advice, I'm no expert on this subject. This ticket was filed more as a response to scipy ticket #581 than a personal need. Anyway, I now also filed an enhancement ticket #610 to scipy.special which may be a better place for improvement if modifying poly1d is a wontfix.

[charris]

I'm closing this because the suggested enhancement should be part of a package in Scipy, not part of Numpy.