

Hi, I noticed some frustrating inconsistencies in the various ways to evaluate polynomials using numpy. Numpy has three ways of evaluating polynomials (that I know of) and each of them has a different syntax:
numpy.polynomial.polynomial.Polynomial: You define a polynomial by a list of coefficients in order of increasing degree, and then use the class’s call() function.
np.polyval: Evaluates a polynomial at a point. First argument is the polynomial, or list of coefficients in order of decreasing degree, and the second argument is the point to evaluate at.
np.polynomial.polynomial.polyval: Also evaluates a polynomial at a point, but has more support for vectorization. First argument is the point to evaluate at, and second argument the list of coefficients in order of increasing degree.
Not only the order of arguments is changed between different methods, but the order of the coefficients is reversed as well, leading to puzzling bugs (in my experience). What could be the reason for this madness? As polyval is a shameless ripoff of Matlab’s function of the same name anyway, why not just use matlab’s syntax (polyval([c0, c1, c2...], x) ) across the board?
_______________________________________________
NumPyDiscussion mailing list
[hidden email]
https://mail.python.org/mailman/listinfo/numpydiscussion


Here's my take on this, but it may not be an accurate summary of the history.
`np.poly<func>` is part of the original matlabstyle API, built around `poly1d` objects. This isn't a great design, because they represent: p(x) = c[0] * x^2 + c[1] * x^1 + c[2] * x^0
For this reason, among others, the `np.polynomial` module was created, starting with a clean slate. The core of this is `np.polynomial.Polynomial`. There, everything uses the convention p(x) = c[0] * x^0 + c[1] * x^1 + c[2] * x^2
It sounds like we might need clearer docs explaining the difference, and pointing users to the more sensible `np.polynomial.Polynomial`
_______________________________________________
NumPyDiscussion mailing list
[hidden email]
https://mail.python.org/mailman/listinfo/numpydiscussion


Thanks, that explains a lot! I didn't realize the reverse ordering actually originated with matlab's polyval, but that makes sense given the onebased indexing. I see why it is the way it is, but I still think it would make more sense for np.polyval() to use conventional indexing ( c[0] * x^0 + c[1] * x^1 + c[2] * x^2). np.polyval() can be convenient when a polynomial object is just not needed, but if a single program uses both np.polyval() and np.polynomail.Polynomial, it seems bound to cause unnecessary confusion.
Max
_______________________________________________
NumPyDiscussion mailing list
[hidden email]
https://mail.python.org/mailman/listinfo/numpydiscussion


> if a single program uses both np.polyval() and np.polynomail.Polynomial, it seems bound to cause unnecessary confusion.
Yes, I would recommend definitely not doing that!
> I still think it would make more sense for np.polyval() to use conventional indexing
Unfortunately, it's too late for "making sense" to factor into the design. `polyval` is being used in the wild, so we're stuck with it behaving the way it does. At best, we can deprecate it and start telling people to move from `np.polyval` over to `np.polynomial.polynomial.polyval`. Perhaps we need to make this namespace less cumbersome in order for that to be a reasonable option.
I also wonder if we want a more lightweight polynomial object without the extra domain and range information, which seem like they make `Polynomial` a more questionable dropin replacement for `poly1d`.
Eric Thanks, that explains a lot! I didn't realize the reverse ordering actually originated with matlab's polyval, but that makes sense given the onebased indexing. I see why it is the way it is, but I still think it would make more sense for np.polyval() to use conventional indexing ( c[0] * x^0 + c[1] * x^1 + c[2] * x^2). np.polyval() can be convenient when a polynomial object is just not needed, but if a single program uses both np.polyval() and np.polynomail.Polynomial, it seems bound to cause unnecessary confusion.
Max
_______________________________________________
NumPyDiscussion mailing list
[hidden email]
https://mail.python.org/mailman/listinfo/numpydiscussion
_______________________________________________
NumPyDiscussion mailing list
[hidden email]
https://mail.python.org/mailman/listinfo/numpydiscussion


I think restricting polynomials to time series is not a generic way and quite specific.
Apart from the series and certain filter design actual usage of polynomials are always presented with decreasing order (control and signal processing included because they use powers of s and inverse powers of z if needed). So if that is the use case then probably it should go under a namespace of `TimeSeries` or at least require an option to present it in reverse. In my opinion polynomials are way more general than that domain and to everyone else it seems to me that "the intuitive way" is the decreasing powers.
For the design
> This isn't a great design, because they represent: > p(x) = c[0] * x^2 + c[1] * x^1 + c[2] * x^0
I don't see the problem actually. If I ask someone to write down the coefficients of a polynomial I don't think anyone would start from c[2].
_______________________________________________
NumPyDiscussion mailing list
[hidden email]
https://mail.python.org/mailman/listinfo/numpydiscussion


“the intuitive way” is the decreasing powers.
An argument against this is that accessing the ith power of x is spelt:
x.coeffs[i] for increasing powers
x.coeffs[i1] for decreasing powers
The former is far more natural than the latter, and avoids a potential offbyone error
If I ask someone to write down the coefficients of a polynomial I don’t think anyone would start from c[2]
You wouldn’t? I’d expect to see
rather than
Sure, I’d write it starting with the highest power, but I’d still number my coefficients to match the powers.
Eric
_______________________________________________
NumPyDiscussion mailing list
[hidden email]
https://mail.python.org/mailman/listinfo/numpydiscussion


Interesting, I wasn't aware that both conventions were widely used.
Speaking of series with inverse powers (i.e. Laurent series), I wonder how useful it would be to create a class to represent expressions with integral powers from m to n. These come up in my work sometimes, and I usually represent them with coefficient arrays ordered like this:
c[0]*x^0 + ... + c[n]*x^n + c[n+1]x^m + ... + c[n+m+1]*x^1
Because then with negative indexing you have:
c[m]*x^m + ... + c[n]*x^n
Still, these objects can't be manipulated as nicely as polynomials because they aren't closed under integration and differentiation (you get log terms).
Max
_______________________________________________
NumPyDiscussion mailing list
[hidden email]
https://mail.python.org/mailman/listinfo/numpydiscussion


Since the one of the arguments for the decreasing order seems to just be textual representation  do we want to tweak the repr to something like
Polynomial(lambda x: 2*x**3 + 3*x**2 + x + 0)
(And add a constructor that calls the lambda with Polynomial(1) )
Eric
“the intuitive way” is the decreasing powers.
An argument against this is that accessing the ith power of x is spelt:
x.coeffs[i] for increasing powers
x.coeffs[i1] for decreasing powers
The former is far more natural than the latter, and avoids a potential offbyone error
If I ask someone to write down the coefficients of a polynomial I don’t think anyone would start from c[2]
You wouldn’t? I’d expect to see
rather than
Sure, I’d write it starting with the highest power, but I’d still number my coefficients to match the powers.
Eric
_______________________________________________
NumPyDiscussion mailing list
[hidden email]
https://mail.python.org/mailman/listinfo/numpydiscussion


Say we add a constructor to the polynomial base class that looks something like this:
 @classmethod def literal(cls, f): def basis_function_getter(self, deg): coefs = [0]*deg + [1] return lambda _: cls(coefs) basis = type('',(object,),{'__getitem__': basis_function_getter})() return f(basis, None) 
Then the repr for, say, a Chebyshev polynomial could look like this:
>>> Chebyshev.literal(lambda T,x: 1*T[0](x) + 2*T[1](x) + 3*T[2](x))
Does this sound like a good idea to anyone?
Max
_______________________________________________
NumPyDiscussion mailing list
[hidden email]
https://mail.python.org/mailman/listinfo/numpydiscussion


I think the `x` is just noise there, especially if it's ignored (that is, `T[0](x*2)` doesn't do anything reasonable). Chebyshev.literal(lambda T: 1*T[0] + 2*T[1] + 3*T[2])
Would work, but honestly I don't think that provides much clarity. I think the value here is mainly for "simple" polynomials.
Say we add a constructor to the polynomial base class that looks something like this:
 @classmethod def literal(cls, f): def basis_function_getter(self, deg): coefs = [0]*deg + [1] return lambda _: cls(coefs) basis = type('',(object,),{'__getitem__': basis_function_getter})() return f(basis, None) 
Then the repr for, say, a Chebyshev polynomial could look like this:
>>> Chebyshev.literal(lambda T,x: 1*T[0](x) + 2*T[1](x) + 3*T[2](x))
Does this sound like a good idea to anyone?
Max
_______________________________________________
NumPyDiscussion mailing list
[hidden email]
https://mail.python.org/mailman/listinfo/numpydiscussion
_______________________________________________
NumPyDiscussion mailing list
[hidden email]
https://mail.python.org/mailman/listinfo/numpydiscussion


Ok I see what you mean. If people really want mathlike symbolic representations for everything it’s probably better to use sympy or something I think the `x` is just noise there, especially if it's ignored (that is, `T[0](x*2)` doesn't do anything reasonable). Chebyshev.literal(lambda T: 1*T[0] + 2*T[1] + 3*T[2])
Would work, but honestly I don't think that provides much clarity. I think the value here is mainly for "simple" polynomials.
Say we add a constructor to the polynomial base class that looks something like this:
 @classmethod def literal(cls, f): def basis_function_getter(self, deg): coefs = [0]*deg + [1] return lambda _: cls(coefs) basis = type('',(object,),{'__getitem__': basis_function_getter})() return f(basis, None) 
Then the repr for, say, a Chebyshev polynomial could look like this:
>>> Chebyshev.literal(lambda T,x: 1*T[0](x) + 2*T[1](x) + 3*T[2](x))
Does this sound like a good idea to anyone?
Max
_______________________________________________
NumPyDiscussion mailing list
[hidden email]
https://mail.python.org/mailman/listinfo/numpydiscussion
_______________________________________________
NumPyDiscussion mailing list
[hidden email]
https://mail.python.org/mailman/listinfo/numpydiscussion
_______________________________________________
NumPyDiscussion mailing list
[hidden email]
https://mail.python.org/mailman/listinfo/numpydiscussion

