Tell me more ×
Computational Science Stack Exchange is a question and answer site for scientists using computers to solve scientific problems. It's 100% free, no registration required.
up vote 4 down vote favorite

I am looking for a very fast polynomial root finder, hopefully with a matlab code. I don't need very accurate results, 2-3 decimal places would be fine. Also the method should be able to optionally take as input an initial guess of the roots. If the initial guess is close, this should speed the process up.

What do I need this for? I'm working on images and the method I'm implementing requires the roots of a polynomial to be found, at every pixel in the image. Given typical images are from 512x512 and up, the method has to be fast at the expense of accuracy if need be.

share|improve this question
What degree is your polynomial? Do you have to find all roots or just a particular close one? Are the roots real? – Jed Brown Jan 16 at 4:22
Performance-wise, how do the builtin functions fzero and roots fare? – Nat Wilson Jan 16 at 5:22
Thanks @JedBrown and @NatWilson for the comments. The polynomial is a real valued trigonometric polynomial in $e^{i\theta}$ so the roots will be complex. Degree of the polynomial for typical usage will vary between 8 and 16. On my computer 512 x 512 runs of roots takes approx. 22 seconds which is a little too long. All roots need to be found however. – geometrikal Jan 16 at 5:43
With being able to input a starting guess, I figure that the location of the roots is only going to change a little from pixel to pixel if the underlying signal is relatively smooth. Thus I'm hoping that be using the result of a nearby pixel as the initial guess should speed things up. – geometrikal Jan 16 at 5:48

3 Answers

up vote 4 down vote

If you have a reasonable starting guess and the function you're trying to find a root of is not overly complex, then you can't beat Newton's method. It's also easy to implement -- a 10 line exercise in most languages.

share|improve this answer
Thanks for the answer, will this handle cases of multiple roots, or if one root disappears or a new one appears? – geometrikal Jan 16 at 5:50
These are questions you can find answered in most introductory books on numerical methods. – Wolfgang Bangerth Jan 16 at 12:09
You can also modify the standard newton method so that it converges globally (i.e. from any initial guess) if you include a line search with backtracking. I don't think newton's method can be used to find all roots though. – Paul♦ Jan 16 at 18:08
If he knows that the degree at the new point is the same degree as the polynomial from the last, it is likely that running newton's method at all the old roots would capture all the new roots (since they are expected to be shifted slightly). There is the issue of finding them at the first point, and what happens when one root disappears and another appears elsewhere. – Godric Seer Jan 16 at 19:59
up vote 3 down vote

This may seem like a somewhat exotic suggestion, but if your polynomials are of fixed, low degree, why not use the eigenvalues of the companion/colleague-matrix as an initial value for the roots?

This may be expensive, but if you have to repeat the operation for slightly modified polynomials, you could use one or two steps of inverse iteration to update the roots. Although each iteration involves a matrix inverse, recall that the colleague/companion matrix, plus the inverse iteration shift, is only a bi-diagonal matrix with an extra column attached, and the system can thus be solved in $\mathcal O(n)$ operations.

share|improve this answer
1  
Not exotic at all; it's the basis for Jenkins-Traub which is a standard algorithm in polynomial root finding. – Jed Brown Jan 16 at 14:01
@JedBrown: Thanks for pointing that out, I had never heard of it before! – Pedro Jan 16 at 17:32
Thanks @Pedro, I believe this is how matlab's roots function does it. Inverse iteration looks promising. – geometrikal Jan 17 at 5:20
this is exactly what "roots" does. the problem is that it might not be optimal for small degree polynomials (overhead for small matrices). – GertVdE Jan 17 at 8:47
@geometrikal: Yes, that is true for the first part. The iterative refinement part is what I had considered "exotic". – Pedro Jan 17 at 9:35
show 1 more comment
up vote 2 down vote

If you don't have initial guesses ready, I would try Bairstow's method. And Newton refinement if needed (but probably not if you're happy with 2-3 digits).

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.