[SciPy-dev] Generic polynomials class (was Re: Volunteer for Scipy Project)

[Charles R Harris]

On Tue, Oct 6, 2009 at 9:02 AM, Anne Archibald <peridot.faceted at gmail.com>wrote:

> 2009/10/6 Charles R Harris <charlesr.harris at gmail.com>:
> >
> >
> > On Mon, Oct 5, 2009 at 10:14 PM, Anne Archibald <
> peridot.faceted at gmail.com>
> > wrote:
> >>
> >> 2009/10/5 Charles R Harris <charlesr.harris at gmail.com>:
> >> >
> >> >
> >> > On Mon, Oct 5, 2009 at 7:29 PM, Anne Archibald
> >> > <peridot.faceted at gmail.com>
> >> > wrote:
> >> >>
> >> >> 2009/10/5 Charles R Harris <charlesr.harris at gmail.com>:
> >> >> >
> >> >> >
> >> >> > On Mon, Oct 5, 2009 at 3:20 PM, Anne Archibald
> >> >> > <peridot.faceted at gmail.com>
> >> >> > wrote:
> >> >> >> For examples where I think a bit of lit review plus implementation
> >> >> >> work might help, I'd say that the orthogonal polynomials could use
> >> >> >> some work - the generic implementation in scipy.special falls
> apart
> >> >> >> rapidly as you go to higher orders. I always implement my own
> >> >> >> Chebyshev polynomials using the cos(n*arccos(x)) expression, for
> >> >> >> example, and special implementations for the others might be very
> >> >> >> useful.
> >> >> >>
> >> >> >
> >> >> > At what order does the scipy implementation of the Chebyshev
> >> >> > polynomials
> >> >> > fall apart? I looked briefly at that package a long time ago, but
> >> >> > never
> >> >> > used
> >> >> > it. I ask so I can check the chebyshev module that is going into
> >> >> > numpy.
> >> >>
> >> >> By n=30 they are off by as much as 0.0018 on [-1,1]; n=31 they are
> off
> >> >> by 0.1, and by n=35 they are off by four - not great for values that
> >> >> should be in the interval [-1,1]. This may seem like an outrageously
> >> >> high degree for a polynomial, but there's no reason they need be this
> >> >> bad, and one could quite reasonably want to use an order this high,
> >> >> say for function approximation.
> >> >>
> >> >
> >> > Hmm, I get an maximum error of about 1e-13 for n=100 using my routine,
> >> > which
> >> > isn't great and can be improved a bit with a few tricks, but is
> probably
> >> > good enough.  That's using simple Clenshaw recursion for the
> evaluation.
> >> > There must be something seriously wrong with the scipy version. I
> >> > routinely
> >> > use fits with n > 50 because truncating such a series gives a good
> >> > approximation to the minmax fit and it's also nice to see how the
> >> > coefficients fall off to estimate the size of the truncation error.
> >>
> >> Upon closer inspection, it seems that scipy starts from the recurrence
> >> relation but then uses it only to extract roots, which it then uses to
> >> construct a poly1d instance. (Which presumably then goes on to convert
> >> it into the 1, x, x**2, ... basis.) This is numerically disastrous.
> >>
> >> The reason it does this is because scipy doesn't (didn't) have a
> >> general-purpose polynomial class that can keep the polynomial in other
> >> representations. What do you think of the idea of a polynomial class
> >> that supports several different bases for the set of polynomials? (At
> >> least, the power basis and Chebyshev polynomials, but possibly also
> >> the other orthogonal polynomials. Translated and scaled bases might
> >> come in handy too.) Too general? Not general enough (barycentric
> >> representations, representation in terms of roots)?
> >
> > Barycentric representation using the type II Chebyshev points would be
> the
> > way to go for *all* the functions without singularities at the ends,
> except
> > I couldn't figure out how to get the remainder during division, and the
> > remainder is needed for converting to the coefficients of the various
> types
> > of series -- it's just like converting a decimal representation to binary
> by
> > repeated division by two and keeping the remainders. That doesn't mean it
> > can't be done, however...
>
> We even already have compiled code for evaluation. But there are
> several polynomial operations I can't see an easy way to handle using
> this. It seems to me that with a polynomial class one should be able
> to easily:
>
>

> * evaluate
>

check


> * add/subtract
>

check


> * multiply by another polynomial
>

of the same family, check

* divide with remainder
>

check

* extract derivatives and antiderivatives, both at a point and as
> polynomials
>

I'd split this. into derivatives/antiderivatives of the series, then use
evaluation at a point (or array of points.


> * express in another polynomial basis
>

see below, but it may be difficult to guarantee stability.

* measure degree
>

len() returns the length of the coefficient array. We could throw a deg
method in there.


> (Have I missed any?)
>
> With barycentric interpolation, it's not clear how to divide or take
> antiderivatives (and hence converting bases becomes awkward).


Yeah, the barycentric form would be a natural for rational functions.
Polynomials are more like integers. Hmm, it strikes me that the terms y_j/(x
- x_j) in the barycentric form are just remainders. Maybe something *can* be
done here...


> On the
> other hand, multiplication is O(n) rather than O(n^2). And I think
> antidifferentiation could be worked out with some cleverness.
>
>
The conversion between polynomial types can also be done by function
composition. For instance, conversion of a polynomial to a Chebyshev series
would be mypoly.val(mycheb), where mycheb is a Chebyshev object
representing  'x'. It goes the other way also. This wouldn't work well for
the Barycentric form because there are still divisions. Or so it seems at
first thought...


> Nevertheless I think one would want to be able to work with orthogonal
> polynomial bases, for example when dealing with an infinite (or
> unspecified) interval.
>
Plus one needs to record all the various extra
> information about each family of orthogonal polynomials (roots,
> weights for integration, etc.)
>
>
Ah, now that becomes more specialized. The weights/zeros for integration
should probably come from the Gaussian integration package. I feel that the
packages in numpy should cover the common uses and scipy should have the
parts for which specialized routines are available.

I'd still like to settle the coefficient access order for the Chebyshev
class. As is, it is like the poly1d class. The constructor and coef
attribute contain the coefficients in high to low order, while indexing like
mypoly[i] goes the other way.

Chuck
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/scipy-dev/attachments/20091006/44e9b93d/attachment.html>