# Polynomial root finding

I have an univariate polynomial of some degree - how do I numerically find all of its real roots?

I never thought I would ask this question - everyone knows how to find polynomial roots, right..? Please, can someone explain to me how this works? This is what I think I understand so far:

Newton method - it can find only one root at a time (if it converges) so in order to list all the roots I need to somehow guess where to start from. How to guess?

Bisection - also finds only one root and it needs two endpoints with opposite signs. How to guess those endpoints? What if the function simply does not have opposite signs, for example $f(x) = x*x$?

Do I need to know at least the interval I want to search in? Say, that I'm only interested in roots from -1000 to 1000? Or is it possible to list all the roots without any initial guesses?

Then I found about http://en.wikipedia.org/wiki/Sturm%27s_theorem but I don't really understand it yet - is it somehow helpful in this context?

Is the approach basically combination of the above? Count the signs with Descartes (what for?), find the intervals with Sturm, check if there are opposite signs, if no - find the extreme, if yes - bisect?

And I also read somewhere that it is possible to covert the polynomial into a system of linear equations (matrix?) and solve that? Is this possible?

(Note: I search this forum and found this but I don't really understand any of it - "Homotopy Continuation"?)

• I'm not an expert in numerical mathematics, but here's an outline of a possible procedure: First, calculate the square-free part of a polynomial $p$ via $\operatorname{gcd}(p,p')$. Then, use Descartes' Rule and Sturm's Theorem to find small-enough intervals for each zero of the polynomial (since it's square-free, the intervals can be made small enough to contain a single zero each); this can be done by bisection. Also, make the intervals so small that any point inside them gives quadratic convergence for Newton iteration, and do Newton iteration from the center point. – Johannes Kloos Feb 20 '13 at 15:45
• @JohannesKloos: why to calculate square-free part? To eliminate multiple roots? And why to use Descartes' Rule - it only gives the number of sign changes, no? "the intervals can be made small enough" - how? by Sturm's Theorem? – Ecir Hana Feb 20 '13 at 17:33
• The square-free part eliminates multiple roots, as you said. This is advantageous for two reasons: a. It allows you to get intervals containing a single zero using Sturm's Theorem. b. I read somewhere that some numerical methods have problems with multiple roots. – Johannes Kloos Feb 20 '13 at 19:20
• The "small enough" thing you can do as follows: As soon as you have seperated the zeroes into intervals (use Sturm's Theorem), it is guaranteed that there is a sign change inside every interval (because there are no repeated zeroes). You can then use, e.g., bisection to find a sub-interval in which Newton's method converges quadratically. – Johannes Kloos Feb 20 '13 at 19:23
• Sturm's method and the algorithms that are sometimes collectively described as Descartes methods (VCA, VAS, Continued Fraction, ...) are orthogonal - no need to use both. – akobel Jul 29 '15 at 8:58

You can use Sturm's Theorem to find initial guesses for Newton's method, along the lines of Johannes's comment.

The matrix formulation you're referring to is the companion matrix of the polynomial. Roots of the polynomial become eigenvalues of the companion matrix, so you can then use any eigenvalue algorithm to find the roots.

However, in my experience, if you need to find all of the roots of an arbitrary polynomial, the Jenkins-Traub algorithm is easily the best, both in terms of efficiency and robustness, and it's what commercial software like Mathematica uses under the hood when you ask it to numerical solve for polynomial roots. The one downside is that it is far from trivial to implement yourself; fortunately it's not hard to find codes on the Internet, in both Fortran and C.

• I'm a bit lost on this - so should I use the method based on Sturm's Theorem or Jenkins-Traub or companion matrix? Could you please somehow compare these methods, if it is possible at all? Also, Mathematica uses "VAS", according to en.wikipedia.org/wiki/… – Ecir Hana Feb 20 '13 at 17:43
• Knowing nothing about exactly what you are trying to do, as a one-size-fits-all method I recommend Jenkins-Traub. I'm not sure what kind of comparisons you want -- for a comprehensive overview of the different method you can read, e.g. amazon.com/Numerical-Methods-Roots-Polynomials-Computational/dp/… – user7530 Feb 20 '13 at 18:01
• For a polynomial, does each of the method gives the same result as other methods? I mean, for edge cases like multiple roots or no opposite signs, do the methods return the same values? Using floating point calculations, which method is the most stable? For small polynomials (of degree 1 - 10), which method is the fastest? – Ecir Hana Feb 20 '13 at 18:12
• No, each method will give different answers. Which is most stable or fastest will vary widely depending on how the polynomials are sampled, how the methods are implemented, and how the results are measured, (suppose someone asked you, "Which car should I buy that's fastest and safest?") and any kind of definitive partial ordering is way beyond the scope of a comment. I recommend posting a new question along the lines of "what are some recent papers comparing the robustness and performance of different polynomial root-finding algorithms" with tag "reference-request". – user7530 Feb 20 '13 at 23:06

W.r.t. "Do I need to know at least the interval I want to search in?": Yes, typically you do. At least for subdivision-based solvers such as Sturm's method or pretty much anything based on Descartes' rule of sign or its variations. For other algorithms like Durand-Kerner, Aberth-Ehrlich or homotopy continuation methods, you generally don't need bounds.

It is, however, fairly easy to at least find some bound for the initial interval: compute any upper bound $M$ for the absolute value of the roots and start subdivision on $[-M,M]$ (for real root isolation) or $[-M,M]+i\cdot [-M,M]$ (for complex root isolation). Famous formulas are Cauchy's root bound (simple) and Fujiwara's root bound (fairly simple and often very accurate). You can also find specialized root bounds for only the real roots, only positive roots, and certainly also for other specific demands.
Obviously, there is a trade-off between complexity of the root bound and the quality of the estimate: the most accurate root bound computation possible is to isolate all roots and take the largest magnitude of any root. ;-)

I'm totally unqualified to answer this, but one simple algorithm that should work is to find starting points for the bisection method by first (recursively) finding the roots of the polynomial's derivative. The derivative is of a lower order (power rule), so the recursion will eventually reach an order you can solve algebraically. All minima and maxima of the original function will be roots of the derivative. So you can just check the signs of these points, and perform bisections where needed.

The exception to this is potential roots where the function is continually increasing or decreasing to either side of the roots of the derivative. This may be the fatal flaw with this algorithm. In practice, however, moving out exponentially until you find a changed sign and can bisect will work for most inputs.

This is probably a bad algorithm. I don't know a lot about numerical computing. But I'm adding it here because it's fairly easy to understand.