Lagrange polynomials and numpy Polynomials

[Jebby993]

My issue is about the documentation of lagrange interpolation here: https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.lagrange.html
lagrange returns a numpy.poly1d that contains "The polynomial’s coefficients, in decreasing powers"

Then, the documentation shows to pass this result to numpy Polynomials. But this class accepts as parameters "Polynomial coefficients in order of increasing degree"
https://numpy.org/doc/stable/reference/generated/numpy.polynomial.polynomial.Polynomial.html

Does it seem like a conceptual error?
Coefficients are the same, but the order is reversed so when you want to interpolate things do not work well

Reproducing code example:

from numpy.polynomial.polynomial import Polynomial
from scipy.interpolate import lagrange
x = np.array([0, 1, 2])
y = x**3
poly = lagrange(x, y)

poly1d([ 3., -2.,  0.]) #3x^2 -2x

Polynomial(poly).coef
array([ 3., -2.,  0.]) #-2x +3

Scipy/Numpy/Python version information:

import sys, scipy, numpy; print(scipy.__version__, numpy.__version__, sys.version_info)
1.7.0 1.20.3 sys.version_info(major=3, minor=7, micro=3, releaselevel='final', serial=0)

[ev-br]

Yup, it needs to be Polynomial(poly.coef[::-1]), cf "Transition guide" at https://numpy.org/doc/stable/reference/routines.polynomials.html (I'd expect Polynomial constructor to recognize its argument to be a poly1d instance, but apparently it does not, probably for good reasons.)

Want to send a quick PR with a fix @Jebby993 ?
Bonus point can be an example with a plot of Polynomial(poly.coef[::-1])(x) vs 3*x**2 - 2*x

[Jebby993]

Thank you @ev-br ,I could try (I'm quite new here on github).
However, I think it would be better to skip the use of Polynomial that could be misleading and exploit numpy.poly1d attributes that are sufficient to show the polynomial coefficient, evaluate the polynomial and eventually compute the roots, without the need to import and use the Polynomial class.

[ev-br]

poly1d is considered legacy and its use in new code is discouraged. This function predates numpy Polynomial, so we cannot just change its return type. All in all, best is to show the recommended usage --- which this example tried to do :-).

[tupui]

@ev-br can we mark these kind of things to deprecate them slowly?

[ev-br]

Problem is, I don't see a clean transition path. Do you?

Just changing the return type is bad; doc notes do not change that; deprecating the whole function and inventing a new one with a different return type feels like an overkill. So my gut feeling is just keep the status quo and improve the documentation for users to show best practices.

[tupui]

Problem is, I don't see a clean transition path. Do you?

Just changing the return type is bad; doc notes do not change that; deprecating the whole function and inventing a new one with a different return type feels like an overkill. So my gut feeling is just keep the status quo and improve the documentation for users to show best practices.

What I mean is that in case of functions that we say are legacy and should not be used in new code, then I would decorate them to raise a deprecation warning for a few release cycle and then delete them... I know we try to not do this, but we could have long cycles to mitigate issues with users.

[charris]

For polynomials whose weight functions have compact support, returning a Chebyshev series would be a good idea. There is a chebinterpolate function that takes a function and returns an Chebyshev series, but only good over the interval [-1, 1]. Something like that could be added to the Chebyshev class as a class method to do the scaling. Chebyshev series tend to be more numerically stable than power series.

[ev-br]

What I mean is that in case of functions that we say are legacy and should not be used in new code, then I would decorate them to raise a deprecation warning for a few release cycle and then delete them... I know we try to not do this, but we could have long cycles to mitigate issues with users.

Yes, but here our function is not deprecated. Instead it's the return type being declared legacy by numpy.

For polynomials whose weight functions have compact support, returning a Chebyshev series would be a good idea.

Agreed, it'd be great. Not in scipy.interpolate.lagrange function however I suspect :-).

[tupui]

What I mean is that in case of functions that we say are legacy and should not be used in new code, then I would decorate them to raise a deprecation warning for a few release cycle and then delete them... I know we try to not do this, but we could have long cycles to mitigate issues with users.

Yes, but here our function is not deprecated. Instead it's the return type being declared legacy by numpy.

To be sure, I am talking about your comment on poly1d here.

[ev-br]

Sure. Chuck @charris , are there any plans w. r. t. np.poly1d? I guess we'd need to coordinate with scipy. special. orthopoly1d if there's action.

[tupui]

Sure. Chuck @charris , are there any plans w. r. t. np.poly1d? I guess we'd need to coordinate with scipy. special. orthopoly1d if there's action.

At some point we have to plan for SciPy 2.0 (just in 2 releases). Such things could be good candidates.

[charris]

are there any plans w. r. t. np.poly1d

Nothing planned that I am aware of. The Chebyshev polynomials have the advantage that the coefficients are easy to compute from values sampled at the Chebyshev points, so would work well with the scipy functions that just compute the polynomial values.