ENH: Add polyvalfromroots to numpy.polynomial
#7542
Conversation
|
Couple of questions.
Apart from that, might look at adapting the way the polynomial from roots is constructed, it tries to group and multiply the terms in a numerically better way. Whether that matters or not for numerical evaluation, I don't know. Note that high degree polynomials grow really fast outside of the set of zeros. |
|
Degrees of a few dozen isn't terribly uncommon in use cases I've come across. Many signal processing filters that are "ideal" in some way or another are specified as as the complex poles and zeros of a rational function on the complex plane (e.g. Chebyshev / Elliptic filters), in so called "zeros, poles, gain" format. Thus, evaluating the frequency response of the filter comes down to evaluating the polynomials in the numerator and denominator. (Currently I'm not sure that the manner of conversion from roots to coefficients is the issue; rather, at higher orders, I think that any conversion is poorly conditioned, such that errors in the coefficients at machine precision levels can lead to a substantial error in the location of the roots, and thereby from the intended numerical values. |
|
I was suggesting doing the same conversion as for creating the polynomial but using the actual roots rather than EDIT: actual |
|
BTW, with centering and scaling you might be able to get a useful polynomial up to about degree 50. I experimented a bit a while back and got up about degree sixty with equally spaced roots on a line. Probably sensitive to root location though. |
|
Should also add something to the 1.12 release notes and |
|
Thanks for the feedback! I've added the requisite edits. I think part of what makes this tough for things like sharp bandpass filters is that the root locations can span multiple orders of magnitude. |
| over the columns of `r` for the evaluation. This keyword is useful | ||
| when `r` is multidimensional. The default value is True. | ||
|
|
||
| .. versionadded:: 1.12 |
There was a problem hiding this comment.
This should be after the summary, it is here for polyval because it was added in a later release.
|
Comments largely addressed, thanks for taking the time to look at this |
polyrootval to numpy.polynomialpolyvalfromroots to numpy.polynomial
|
|
||
| If `r` is a 1-D array, then `p(x)` will have the same shape as `x`. If | ||
| `r` is multidimensional, then the shape of the result depends on the | ||
| value of `tensor`. If `tensor` is true the shape will be r.shape[1:] + |
|
Thanks for the feedback! When making these changes, I also found that it failed for scalar Incidentally, this brought to my attention that |
| return np.ones_like(x) | ||
| if np.isscalar(x) or isinstance(x, (tuple, list)): | ||
| x = np.array(x, ndmin=1, copy=0) + 0.0 | ||
| if x.size == 0: |
There was a problem hiding this comment.
If x is not an array it doesn't have the attribute size. Same for x.ndim below. Note how polyval handles this, just assume tensor=True. I'll guess that you don't need to check for scalar above if you do this. You should also test the scalar case for x. Case where x is neither a scalar nor tuple/list would be x an instance of Polynomial.
|
Hmm, and I must apologize. Now that I've understood what I did with |
|
☔ The latest upstream changes (presumably #6872) made this pull request unmergeable. Please resolve the merge conflicts. |
As one can easily encounter when working with high-order signal processing filters, converting a high-order polynomial from its roots to its polynomial coefficients can be quite lossy, leading to inaccuracies in the filter's properties. This PR adds a new function, `polyrootval` - based on `polyval` - that evaluates a polynomial given a list of its roots. The benefit of calculating it this way can be seen at scipy/scipy:6059. Some tests are included, as well.
|
Thanks for the clarification about |
| >>> polyvalfromroots(b, r, tensor=False) | ||
| array([-0., 0.]) | ||
| """ | ||
| r = np.array(r, ndmin=1, copy=0) + 0.0 |
There was a problem hiding this comment.
I think this should be like
r = np.array(r, ndmin=1, copy=0)
if r.dtype.char in '?bBhHiIlLqQpP':
# astype fails with NA
r = r + 0.0
if isinstance(x, (tuple, list)):
x = np.asarray(x)
if isinstance(x, np.ndarray) and tensor:
r = r.reshape(r.shape + (1,)*x.ndim)
if isinstance(x, np.ndarray) and x.size == 0 and not tensor:
raise ValueError("Empty x cannot broadcast when tensor == False")
return np.prod(x - r, axis=0)
The error situation is actually a bit tricky, as sometimes broadcasting empty x does work, but it is a sometimes thing, for instance r = [1] works, but r = [1,2] does not. The first only works because the number of elements in the row is 1, which is a special "broadcasting" value. So probably better to just error out in that case. The prod function already returns 1 for empty first axis, so no need to special case that. For testing you can reshape empty arrays to get various other empty arrays. For instance
In [47]: array([]).reshape(0,2,3).shape
Out[47]: (0, 2, 3)
There was a problem hiding this comment.
Maybe make that
if isinstance(x, np.ndarray) and x.ndim >= r.ndim:
raise ValueError("x.ndim must be < r.ndim when tensor == False")
As one can easily encounter when working with high-order signal processing
filters, converting a high-order polynomial from its roots to its polynomial
coefficients can be quite lossy, leading to inaccuracies in the filter's
properties.
I'm proposing a new function,
polyrootval- based onpolyval- thatevaluates a polynomial given a list of its roots. The benefit of calculating it
this way can be seen at scipy Issue #6059. Some tests are included, as well.
There are a few things where I'm not 100% sure what the best choice is, though
zeros_like(x); maybe there's another value that makes more sense, or even an error?-x**n?roots=[0]*ncan only give youx**n.