The representation of polynomials with Polynomial class using "domain" and "window" is not clearly explained in the class documentation

[alexeymuranov]

>>> from numpy.polynomial import Polynomial as P
>>> p = P([1,1], domain=[0,1])
>>> p.linspace(2)
(array([ 0.,  1.]), array([ 0.,  2.]))
>>> p.linspace(2, domain=[0,1])
(array([ 0.,  1.]), array([ 0.,  2.]))

Expected result: (array([ 0., 1.]), array([ 1., 2.])).

This is with Python 3.6.0, NumPy 1.12.0.

---

Update. A simpler example of unexpected behaviour:

>>> from numpy.polynomial import Polynomial as P
>>> p = P([0,1], domain=[0,1])
>>> p(0)
-1.0

The documentation of domain and window parameters in Polynomial:

domain : (2,) array_like, optional
Domain to use. The interval [domain[0], domain[1]] is mapped to the interval [window[0], window[1]] by shifting and scaling. The default value is [-1, 1].

window : (2,) array_like, optional
Window, see domain for its use. The default value is [-1, 1].

[charris]

Unexpected, but correct ;) The polynomial coefficients give 1 + x, which evaluated on the default window [-1, 1] ranges from 0 to 2. The domain [0, 1] is mapped to the window, where the polynomial is evaluated, so p(0) == 0 and p(1) == 2. Things are done that way for numerical reasons. For what you expect, define the window equal to the domain so the map is the identity.

In [1]: from numpy.polynomial import Polynomial as P

In [2]: P([1, 1], domain=[0, 1], window=[0, 1]).linspace(2)
Out[2]: (array([ 0.,  1.]), array([ 1.,  2.]))

Or don't bother with the domain at all

In [3]: P([1, 1]).linspace(2, domain=[0, 1])
Out[3]: (array([ 0.,  1.]), array([ 1.,  2.]))

[alexeymuranov]

If 1 + x is evaluated on [-1, 1] (despite domain=[0, 1]), then the first array should be array([ -1., 1.]). Currently the values in the two arrays do not match.

[alexeymuranov]

Ok, i see my mistake: currently

>>> from numpy.polynomial import Polynomial as P
>>> p = P([1,1], domain=[0,1])
>>> p(0)
0.0

Then i consider this a documentation issue, the current documentation is quite unclear:

domain : (2,) array_like, optional
Domain to use. The interval [domain[0], domain[1]] is mapped to the interval [window[0], window[1]] by shifting and scaling. The default value is [-1, 1].

window : (2,) array_like, optional
Window, see domain for its use. The default value is [-1, 1].

I'll change the title of this issue to make it about documentation and will think about a documentation patch.

Nevertheless, i would expect the window to match the domain if only domain is given explicitly.

[mhvk]

@alexeymuranov - if you can clarify the documentation, that would be great; I've been confused about this a lot. Frankly, easiest might simply be to describe the scaling that actually happens in the form of an equation, and then give an example or two. And/or add a link to polyutils.mapdomain, which implements exactly this (and has the equations, though I think it may not actually be included in the documentation by default). Or, perhaps, just a link to where the explanations actually are: https://docs.scipy.org/doc/numpy/reference/routines.polynomials.classes.html#fitting

p.s. Note that the same documentation would be helpful in all the polynomical classes. Possibly it is best added to _polybase.ABCPolyBase with the other classes adjusted so they just add some class-specific stuff rather than repeat the whole docstring.

[charris]

Nevertheless, i would expect the window to match the domain if only domain is given explicitly.

That would ruin the whole purpose of using [-1, 1] as the window, which is that the power series is better numerically conditioned over that interval. That is clearer if you look at the Chebyshev or Legendre polynomials where the weight has support [-1, 1], so that is the interval in which they possess orthogonality and other good properties.

[alexeymuranov]

@charris, the special role of [-1, 1] for Chebyshev or some other polynomial classes IMO should not affect the interface of entire Polynomial class that intends to implement arbitrary polynomials. I think

>>> from numpy.polynomial import Polynomial as P
>>> p = P([0,1], domain=[0,1])
>>> p(0)
-1.0

is quite confusing.

A different but related issue is how polynomials are printed by default. In the above example one gets:

>>> print(p)
poly([ 0.  1.])
>>> print("p(%g) = %g" % (0, p(0)))
p(0) = -1

which looks like an utter nonsense :).

I suppose i was confused by the suggestive domain parameter name. I didn't realize that it only makes sense in combination with window parameter. I used domain because i hoped that this would liberate me from having to supply domain parameter to p.linspace.

[alexeymuranov]

@mhvk, the current behaviour is in fact quite easy to explain: the internal representation of a polynomial p in Polynomial class keeps it in the form of a composition of an affine polynomial m and some other polynomial q: p(X) = q(m(X)); the affine polynomial m in turn is stored as a pair of intervals [a, b], [c, d] such that m(a) = c and m(b) = d, and the polynomial q is stored in the "usual" way as the list of its coefficients. I would need to think how to write it down better.

In other words,

p = Polynomial([2, 0, 1], domain=[0,1], window=[-1,1]))

represents the actual polynomial p(X) = 2 + (-1 + 2X)^2 = 3 + 4X + 4X^2, and print(p) should output something reminiscent of 3 + 4X + 4X^2, or otherwise not skip the domain and window values in the output.

[charris]

which looks like an utter nonsense :).

Looks exactly like I would expect ;) You are arguing apriori assumptions, I would argue that such assumptions should be adjusted pending experience with the actual mechanism, for which there is, IMHO, good reason.

[alexeymuranov]

Looks exactly like I would expect ;)

The value at 0 of the polynomial with coefficient 0, 1 (the polynomial X) is -1??

[charris]

On Thu, Aug 10, 2017 at 8:55 AM, Alexey Muranov ***@***.***> wrote: Looks exactly like I would expect ;) The value at 0 of the polynomial with coefficient 0, 1 (the polynomial X) is -1??

Assumes facts not in evidence. Chuck

[charris]

For instance, here is the condition number of a Vandermonde matrix of degree 20 sampled at 100 equally spaced points over [-x, x]. It demonstrates why using the interval [-1, 1] makes sense.

[alexeymuranov]

I am complaining about the interfaces (the API and the print output) and about documentation, not about algorithms.

[mhvk]

@alexeymuranov - changing the API is really not an option at this point, it would just break code. But adding to the documentation something that would help someone like you, for whom the current logic is not obvious, as well as for, who has been continuously confused by what the "domain" and "window" really refer to (it is mostly which is which that I find confusing), is very helpful.

[eric-wieser]

Note that print(repr(p)) already shows the domain and window information

[charris]

If you want plain old polynomials, do not specify the domain. That is the default and was intentional. If you want to use more sophisticated things, then domain and window are available, probably in that order. For things like fitting, where the domain is crucial for accurate results, the domain is adusted by the algorithm, but can be overridden.

[eric-wieser]

I think I can see an argument for requiring (read: warning otherwise) that the domain and window are either not specified or both specified? That seems to be where the confusion comes from.

[charris]

@eric-wieser The window argument was a late add on that came in with the Hermite and Laguerre polynomials in an attempt to deal consistently with with the already existing domain. I would rate it as mostly a failure, and if I had it to do over, and may do in future, I would use scale and offset for those functions. In any case, the window argument is almost never used for power, Chebyshev, and Lagrange series and it is rather less important on the user side, as opposed to the internal side, than domain. I would rather it not be required.

[StanczakDominik]

While I agree with the need for numerical stability, as someone used to np.polyfit and wanting to use Polynomial.fit to try out the new code as recommended, best practices and all that, I don't think #9536 fixed this moment of confusion (or rather, one hour of confusion, because I had no idea the issue lay squarely with numpy and not with some weird interaction between numpy and xarray):

>>> np.polyfit(x.time, x, 1)                                                                                                                          
array([-1.00082196e+00,  5.79168388e-10])
# coefficients as I had intended them                                                                                                         
>>> np.polynomial.Polynomial.fit(x.time, x, 1)                                                                                                        
Polynomial([-2.49626322e-07, -2.50205490e-07], domain=[0.e+00, 5.e-07], window=[-1.,  1.])
# wut?

It is not at all obvious that the first object in the repr is not the coefficients in the data domain, but the coefficients in the scaled domain. I think this deserves something along the lines of

Please note that Polynomial.coef are coefficients in the scaled domain, and that to get coefficients from the data domain, you should use Polynomial.convert().coef.

right in the middle of docs for the Polynomial class AND in the fit method.

This should probably be a short PR instead of a comment rant, shouldn't it...?