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

Anne Archibald peridot.faceted@gmail....
Wed Oct 14 23:53:55 CDT 2009

2009/10/14 Charles R Harris <charlesr.harris@gmail.com>:
> On Wed, Oct 14, 2009 at 7:09 PM, Anne Archibald <peridot.faceted@gmail.com>
> wrote:
>> I think passing it to the constructor is kind of a side issue (as I
>> said, users are never supposed to call the constructor); if I
>> understand correctly, what you're objecting to is the fact that every
>> polynomial has a .basis attribute that points to the basis its
>> coefficients represent.
> Wait a minute, let's get this straight. Users aren't supposed to use the
> class object, but only instances generated from some primeval instance
> created out of nothingness in the module? That sounds like a factory object,
> but not of the traditional sort. Mathematicians may theoretically start with
> the empty set and a couple of axioms, but most don't work that way ;) How
> does repr work, since it is supposed to generate a string that can be run to
> recreate the object? Doesn't it need to call the constructor? Same with
> unpickling.

Here's a typical use:

b = ChebyshevBasis((0,1))

X = b.polynomial([0,5,0.5])

f = polyfit(basis=b, degree=1, x=[0,0.5,1], y=[1,2,1])

Y = X**2 + f

print Y([0,0.5,1])

The constructor does indeed get called, but by the basis object
(inside its polynomial method). It's still not possible to construct a
polynomial without specifying its basis, of course; after all a
coefficient array is meaningless without a basis in which to interpret
it. So this may not resolve your objection to passing a basis to the

>> Certainly every polynomial needs to know what basis it is in, and many
>> users will also need to know. But it would be possible to create a
>> class for every different basis and store the basis object in the
>> class, so that a polynomial contains only its coefficients and a
>> pointer to its class. Then you could extract a polynomial's basis by
>> doing, say, mypoly.__class__.basis. Apart from the ugliness of doing
>> this, I find it rather unpleasant to have potentially thousands of
>> class objects floating around
> Well, identifying different classes by attaching a name attribute, in this
> case a basis, is how some folks identified classes back in the paleolithic
> period of C++, Borland comes to mind, IIRC, they used to attach an
> identifier to plain old structures too. It works, sorta, but it is just a
> lightweight version of having lots of different classes.  But it does seem
> to me that your contruction is oriented towards Lagrange bases and
> barycentric interpolation and the desire to have it also incorporate graded
> polynomial basis and such has complicated it. It would be simpler to just
> build a different classes for the different sorts of polynomials.

Certainly my design is intended to accommodate the Lagrange
representation. But it is also shaped by a desire to be able to deal
sensibly with families of orthogonal polynomials, providing auxiliary
information about them. There is also a certain amount of code that
can be shared between the Lagrange basis polynomials and other
polynomial implementations - all the type checking plumbing, for
example. I'd be reluctant to pull them apart and duplicate that code.
(I also rather like being able to use isinstance(x,Polynomial) but I
realize that's not as pythonic as one might wish.)

>> - remember PowerBasis(0) is different
>> from PowerBasis(0.32). Also, you have to do some careful work to check
>> polynomial compatibility: people can, and often will, create different
>> basis objects that represent the same object. Either you need to
>> implement a flyweight pattern for your basis objects (or classes -
>> ugh), or you need to have a way of detecting when two polynomials with
>> different classes are actually compatible (because the Basis objects
>> used to manufacture the classes are equal but not the same object).
> How do different classes end up with the same basis. And where do all the
> different
> basis come from if the user doesn't generate them and call the constructor?

Well, in an ideal world, users would simply create one basis and use
it for all the polynomials they produce in that basis. But it's
already clear from writing the test code that we're going to have
things like (in my syntax):

p1 = PowerBasis(center=3).polynomial([1,1])
... many lines of code ...
p2 = PowerBasis(center=3).polynomial([3,5,6])

In this situation p1 and p2 have received different (i.e.
id(A)!=id(B)) basis objects that are actually equal.

As for different classes with the same basis, I think I misunderstood
your proposal. Since your code requires bases to be specified
completely by an implementation and an interval, and you store the
interval in the polynomial object rather than the class, you never do
end up with many classes. If I've understood you correctly, you do
require the user to specify an interval any time they create a
polynomial (inlcuding zero(), one(), and X()), and this is then
checked for compatibility.

I'm not sure I see how your code can be adapted to sensibly handle the
power basis (specifying a center and maybe a scale) or some of the
orthogonal polynomial bases (some of which need an endpoint and scale,
or possibly also one or two parameters specifying the family of

>> It just seems simpler to me to have the basis as a full-fledged
>> attribute of the polynomial class. Given that polynomial classes are
>> supposed to be immutable, that means it has to be passed in the
>> constructor.
> I think is would be simpler to seek less generality. You have an abstract
> Basis class, then a derived GradedBasis class, then a factory function to
> generate a GeneratedBasis class derived from GradedBasis,  and that class
> has to be instantiated somewhere. Then there is a Polynomial class that
> accepts a basis (or basis class?). How is that simpler then a factory
> function that just builds specific classes for specific instances of graded
> polynomials out of a few defined functions and inherits only from object?
> The setup you have looks a bit over designed.

One reason it looks over-designed is that it only has three polynomial
types; it should at least have one more for families of orthogonal
polynomials, and perhaps another for Bernstein polynomials. (The
orthogonal polynomials will inherit from GradedBasis but probably
can't be produced using GeneratedBasis; the Bernstein basis is not
graded.) The GeneratedBasis I would be inclined to do without, except
that if people want to easily make polynomials given collections of
coefficient-operating functions, it will let them do that. I should
point out that, for example, there's no need to implement chebadd and
chebsub, since those implementations are automatically provided by

Leaving aside the GeneratedBasis, I think it's fairly simple.

>From a user's point of view, they pick a polynomial representation
they want to work with by creating an appropriate Basis object, then
use that to create the polynomial, perhaps directly from coefficients
(with basis.polynomial()) or roots (basis.from_roots()), perhaps using
existing ones (basis.X(), basis.convert(other_polynomial)), or by
using some polynomial-returning function (e.g. poly.polyfit(basis,
...)). They can then use that polynomial as an immutable type.

>From an implementer's point of view, adding a new polynomial type is
simple: they define a new basis class inheriting from the most similar
existing basis class, and implement a few operations on coefficient
arrays. Addition and power they can leave out; if they don't care
about division they can leave it out; if they implement division or
companion matrices they get root-finding for free. They'll probably
want to write a routine to convert arbitrary polynomials to their
representation; the other direction is already implemented.

If they've got low-level functions like chebmul to work with, they can
simply write MyBasis = GeneratedBasis(...functions...) and lo and
behold they've got a reasonably functional polynomial class. Or if we
ditch that as too messy, they write a class with some one- and
two-line wrappers.

If you want to compare our two designs, mine uses a single class for
all polynomials and uses the different classes of Basis to provide
different implementations, while yours puts the code directly in the
polynomial classes. Mine uses inheritance to share code between
representations, while yours shares code by writing it into the class
generator function.

I think that to the extent mine is more complicated, it's because it
wants to allow more different implementations of polynomial classes.
Yours is very specialized to orthogonal polynomials defined on a
finite interval. You'll need to generalize it somewhat to handle
power-basis polynomials at least. I find the class-generating code
fairly opaque, and I'm not sure how well it will interact with
debugging and introspection. If these aren't a problem, and we really
don't need more polynomial types, then maybe it does make sense to go
with yours as the simpler.

I'd argue, though, that a general, efficient, numerically-stable
polynomial class could replace KroghInterpolator and
BarycentricInterpolator and serve to make PiecewisePolynomial more
powerful, for example. The orthogonal polynomial code in scipy.special
could return polynomial objects, ideally in their own basis, that
allowed fast evaluation, integration, differentiation, and whose basis
objects carried information about the roots and weights needed for
Gaussian quadrature. (In fact, those basis objects themselves could
live in scipy.special and serve to produce the orthogonal polynomials
as needed.)


More information about the Scipy-dev mailing list