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

[Anne Archibald]

2009/10/6 Charles R Harris <charlesr.harris at gmail.com>:
>
>
> 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.

I put them together because I thought that on the one hand one would
want convenience functions, and on the other hand in some
representations it might be more efficient to evaluate the
derivative/antiderivative at a point than it is ton construct the
polynomial representation. Though come to think of it, that's probably
only true if the derivative/antiderivative evaluation code is written
in C. And in any case I doubt it will be a big win.

>>
>> * express in another polynomial basis
>
>
> see below, but it may be difficult to guarantee stability.

Sure, some polynomial bases are going to give trouble. But I think
that it shouldn't be too hard to convert without introducing much
additional instability beyond what the representations themselves
imply.

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

Sure. I was imagining an implementation (for example) that uses
barycentric interpolation at a power-of-two number of points, so that
raising the number of interpolating points is a rare operation and you
can always reuse the values you have. But then getting the degree
becomes a numerically-unstable operation, which isn't very nice.

>>
>> (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...

I don't know that it's such a common operation that we need special
syntax for it. But if we had classes for polynomial representations
(ChebyshevBasis, PowerBasis, Barycentric) it would indeed be natural
to use call syntax on these classes to transform polynomials.

>> 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.

Well, it's reasonable, if one has Chebyshev polynomials as basis
functions, to want to know where their zeros and maxima/minima are.
It's also not a difficult thing to do (i.e. no special functions
needed), since they can be expressed in terms of cosines.
Programmatically it's a little awkward, since they are associated to
Chebyshev polynomials in general, not to any particular Chebyshev
polynomial.

> 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.

Personally I think that the high-to-low order in poly1d is a mistake,
but at this point I think it may be one we're stuck with. But since
chebyshev polynomials are all new, they do offer the chance to start
fresh. Low-to-high order also has the advantage that you could
implement access to higher coefficients returning zero, if you wanted,
so that they could be treated something like their mathematical notion
of series that happen to terminate.  That said, I haven't used poly1d
much at all, so I'm not necessarily the best person to ask.

Anne