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

[Charles R Harris]

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

> 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.
>
>
It's not a special syntax, it's just using a class instead of a number as
the argument, the same (python) code works for both. Well, it doesn't work
for the poly1d class because poly1d exports an array interface and behaves
badly. But it works for a properly constructed class.


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

You can always call the roots method, which works for Chebyshev series as
well as for isolated terms. It doesn't provide the weights, though.


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


Interesting thought.


> 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.
>
>
I would like to make the change, but I would like some sort of consensus
first. Of course, I will need to rewrite all the code/documentation if it
happens ;)

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