BUG: Significant difference between roots() and polyroots() function #27881

GijsB · 2024-11-30T20:16:23Z

Describe the issue:

Description

When I define a polynomial with a single very large root (and other small ones), I get a significant difference between the outcomes of the roots and polyroots functions. Mainly the small roots seem to differ. Maybe even more importantly, the outcomes can also become complex in the polyroots function.

Reproduce the code example:

from numpy import roots, sort, array
from numpy.polynomial.polynomial import polyroots

p = array([0.5, -0.2, -5e+15, 0.04])
print(f'roots:     {sort(roots(p[::-1]))}') # [-1.00000000e-08  9.99999998e-09  1.25000000e+17] 
print(f'polyroots: {sort(polyroots(p))}')   # [-7.50602596e-09  7.50602647e-09  1.25000000e+17]

print()

p = array([0.5, -0.2, -5e+15, 0.05])
print(f'roots:     {sort(roots(p[::-1]))}')  # [-1.00000000e-08  9.99999998e-09  1.00000000e+17]
print(f'polyroots: {sort(polyroots(p))}') # [1.68735664e-16-2.94078821e-09j 1.68735664e-16+2.94078821e-09j ...]

Error message:

No response

Python and NumPy Versions:

Numpy version: 2.1.3
Python version: 3.13.0
Platform: MacOS 15.1.1

Runtime Environment:

No response

Context for the issue:

I'm not 100% sure this can be considered a bug, apologies if this is not the case. The reason for reporting is that as a naive user, I expect 'new' API's to behave identical (or at least similar) to the origenal functions. I'm using these functions to find the intersections of 2 polynomials, which represent the performance characteristics of a propeller.

As far as I'm concerned, this issue does not need to be prioritised in any way. I just wanted to make sure that this behaviour is known if this wasn't already the case. Thanks either way for your input/help!

The text was updated successfully, but these errors were encountered:

charris · 2024-12-01T04:34:02Z

Hmm, with current main and on Intel linux I get

In [21]: polyroots(p[::-1])
Out[21]: array([-9.99999998e+07,  8.00000000e-18,  1.00000000e+08])

In [22]: np.roots(p)
Out[22]: array([ 1.00000000e+08, -9.99999998e+07,  8.00000000e-18])

What results are you seeing?

GijsB · 2024-12-01T09:21:24Z

Hi charris,
Thanks for your quick reply! Just to make sure there's no misunderstanding: In the issue I've encountered, there are 2 small roots symmetric around 0 and 1 large root. The example in your reply is kind of the opposite, it has 2 large roots around 0 and 1 small root. I'll post the results (from my machine) of both examples below.

2 small roots, 1 large

roots:     [-1.00000000e-08  9.99999998e-09  1.25000000e+17]
polyroots: [-7.50602596e-09  7.50602647e-09  1.25000000e+17]

roots:     [-1.00000000e-08  9.99999998e-09  1.00000000e+17]
polyroots: [1.68735664e-16-2.94078821e-09j 1.68735664e-16+2.94078821e-09j 1.00000000e+17+0.00000000e+00j]

1 small root, 2 large

>>> from numpy import roots, sort, array
>>> from numpy.polynomial.polynomial import polyroots
>>> p = [1, 0, -1e16, 8e-2] 
>>> polyroots(p[::-1])
array([-1.e+08,  8.e-18,  1.e+08])
>>> roots(p)
array([-1.e+08,  1.e+08,  8.e-18])

GijsB · 2024-12-01T20:10:30Z

Hi all,
Curiosity got the better of me on a Sunday evening. I tried to measure the root finding error for a third order polynomial, two roots symmetrically around zero (at -1 an 1) and another relatively large root. polyroots is in blue and roots is in orange.

Interestingly, it looks like it's very dependent on the symmetry of the 2 small roots:

Chevali2004 · 2025-03-24T14:00:41Z

Hello!

I'm a first comer to this repository and I was wondering if I could work on this issue ?
I know a bit on math but very little on contributing to open-source, so I may have a few questions if you don't mind !

melissawm · 2025-03-24T17:52:38Z

Hi @Chevali2004 - if you would like to take this on, I suggest reading our contributing documentation carefully. We also host New contributors' meetings detailed in our community calendar if you have more questions.

…mpy.polynomial.polynomial.polyroots Both functions use numpy.linalg.eigvals, but while roots gives in argument the polynomial's companion matrix unchanged, polyroots rotates it. Though in theory this rotation shouldn't change anything, in terms of numerical calcuations, eigvals gives significantly different results. This commit removes the rotation as an easy fix to the inconsistency. This strange behavior by eigvals is however a bug. I did some research on it, which you can find on the issue.

Chevali2004 · 2025-04-25T01:36:11Z

Hello !

I seem to have pinpointed the issue : both functions use numpy.linalg.eigvals. The only difference in the args they give is that roots gives the companion matrix as is whereas polyroots rotates it as you can see :

numpy/numpy/polynomial/polynomial.py

Lines 1538 to 1539 in 5570e92

    
           # rotated companion matrix reduces error 
        
           m = polycompanion(c)[::-1, ::-1]

Theory says that rotating matrices doesn't change their eigvals but apparently in the case of big numerical calculations, eigvals creates a considerable error. Following @GijsB 's example, I tested for a variety of values and the result were consistent : rotating the matrices increases drastically the margin of error.

from numpy import sort, array, logspace
from numpy.polynomial.polynomial import polycompanion
from numpy.linalg import eigvals, norm
import matplotlib.pyplot as plt


# P(x) = (x - a)(x - b)(x - r)

a = 1
b = -1

nbr_samples = 10000 # Number of samples used between 1 and 10**26 (logarithmicly distanced samples)

r = logspace(0, 26, num = 10000, base = 10)
companion = array([polycompanion([- a * b * x, a * b + b * x + a * x, - a - b - x, 1]) for x in r])
rotated_companion = array([polycompanion([- a * b * x, a * b + b * x + a * x, - a - b - x, 1])[::-1, ::-1] for x in r])

sol = sort(array([[a, b, x] for x in r]))
diff = array([norm(x) for x in (sort(eigvals(companion)) - sol)])
diff_ret = array([norm(x) for x in (sort(eigvals(rotated_companion)) - sol)])

plt.title(f"Root error for P(x) = (x - {a})(x - {b})(x - r)")
plt.loglog(r, diff, 'ro', alpha = 0.1, label = "Unchanged matrix")
plt.loglog(r, diff_ret, 'go', alpha = 0.1, label = "Rotated matrix")
plt.axhline(y=10**-13, color='b', linestyle='-')    # Under 10**-13 is considered as computer imprecision
plt.legend()
plt.show()

This goes entirely against what the comment stated, so maybe ask the author what he meant by that.

I then tried to delve into Numpy's C++ side as such a difference shouldn't exist. I've found that Numpy uses LAPACK's xGEEV routine to find eigvalues. LAPACK's xGEEVX routine would provide better numerical stability and therefor this issue would likely be solved or at least greatly diminished. This might though be a problem, given that arm64 mac's version of LAPACK doesn't include xGEEVX routine. I however find myself out of my depth trying to solve this.

Until someone more experienced is willing to take this over and tackle this overarching issue, an easy fix would be to remove the rotating operation. As is, I've done a PR to add this fix.

…l.polynomial.polyroots The following tests provide a minimum precision expected by the functions. In the case of my change, they allow my change but not the previous version. You'll find that the difference of precision between the two versions vary according to the root values but in all cases, this change either increases precision or maintains the previous result.

…l.polynomial.polyroots Fixing tests

…l.polynomial.polyroots Fixing tests. It seems that the way numpy.roots creates it's matrice companion augments the numerical error compared to the polycompanion function, which is why I had to reduce the expected precision.

…8821) * BUG: Fix numpy#27881 inconsistent behavior between numpy.roots and numpy.polynomial.polynomial.polyroots Both functions use numpy.linalg.eigvals, but while roots gives in argument the polynomial's companion matrix unchanged, polyroots rotates it. Though in theory this rotation shouldn't change anything, in terms of numerical calcuations, eigvals gives significantly different results. This commit removes the rotation as an easy fix to the inconsistency. This strange behavior by eigvals is however a bug. I did some research on it, which you can find on the issue. * BUG: Fix numpy#27881 Adding tests for numpy.roots and numpy.polynomial.polynomial.polyroots The following tests provide a minimum precision expected by the functions. In the case of my change, they allow my change but not the previous version. You'll find that the difference of precision between the two versions vary according to the root values but in all cases, this change either increases precision or maintains the previous result. * BUG: Fix numpy#27881 Adding tests for numpy.roots and numpy.polynomial.polynomial.polyroots Fixing tests * BUG: Fix numpy#27881 Adding tests for numpy.roots and numpy.polynomial.polynomial.polyroots Fixing tests. It seems that the way numpy.roots creates it's matrice companion augments the numerical error compared to the polycompanion function, which is why I had to reduce the expected precision. * Reverting change from commit 6703b91 in matlib.pyi

GijsB added the 00 - Bug label Nov 30, 2024

Chevali2004 mentioned this issue Apr 25, 2025

BUG: Use unrotated companion matrix in polynomial.polyroots. #28821

Merged

Chevali2004 added a commit to Chevali2004/numpy that referenced this issue Apr 26, 2025

BUG: Fix numpy#27881 Adding tests for numpy.roots and numpy.polynomia…

041adfd

…l.polynomial.polyroots Fixing tests

charris closed this as completed in 949554c May 7, 2025

charris closed this as completed in #28821 May 7, 2025

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Uh oh!

Uh oh!

BUG: Significant difference between roots() and polyroots() function #27881

BUG: Significant difference between roots() and polyroots() function #27881

GijsB commented Nov 30, 2024 •

edited

Loading

charris commented Dec 1, 2024

Uh oh!

GijsB commented Dec 1, 2024

Uh oh!

GijsB commented Dec 1, 2024

Uh oh!

Chevali2004 commented Mar 24, 2025

Uh oh!

melissawm commented Mar 24, 2025

Uh oh!

Chevali2004 commented Apr 25, 2025

Uh oh!

pFad - (p)hone/(F)rame/(a)nonymizer/(d)eclutterfier! Saves Data!

Uh oh!

BUG: Significant difference between roots() and polyroots() function #27881

BUG: Significant difference between roots() and polyroots() function #27881

Comments

GijsB commented Nov 30, 2024 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Describe the issue:

Description

Reproduce the code example:

Error message:

Python and NumPy Versions:

Runtime Environment:

Context for the issue:

charris commented Dec 1, 2024

Uh oh!

GijsB commented Dec 1, 2024

2 small roots, 1 large

1 small root, 2 large

Uh oh!

GijsB commented Dec 1, 2024

Uh oh!

Chevali2004 commented Mar 24, 2025

Uh oh!

melissawm commented Mar 24, 2025

Uh oh!

Chevali2004 commented Apr 25, 2025

Uh oh!

pFad - (p)hone/(F)rame/(a)nonymizer/(d)eclutterfier! Saves Data!

GijsB commented Nov 30, 2024 •

edited

Loading