5
$\begingroup$

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.

Improve this question
$\endgroup$
4
  • $\begingroup$ What degree is your polynomial? Do you have to find all roots or just a particular close one? Are the roots real? $\endgroup$
    – Jed Brown
    Jan 16 '13 at 4:22
  • $\begingroup$ Performance-wise, how do the builtin functions fzero and roots fare? $\endgroup$
    – Nat Wilson
    Jan 16 '13 at 5:22
  • $\begingroup$ 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. $\endgroup$
    – geometrikal
    Jan 16 '13 at 5:43
  • $\begingroup$ 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. $\endgroup$
    – geometrikal
    Jan 16 '13 at 5:48
6
$\begingroup$

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.

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

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.

Improve this answer
$\endgroup$
4
  • $\begingroup$ Thanks for the answer, will this handle cases of multiple roots, or if one root disappears or a new one appears? $\endgroup$
    – geometrikal
    Jan 16 '13 at 5:50
  • $\begingroup$ These are questions you can find answered in most introductory books on numerical methods. $\endgroup$
    – Wolfgang Bangerth
    Jan 16 '13 at 12:09
  • $\begingroup$ 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. $\endgroup$
    – Paul ♦
    Jan 16 '13 at 18:08
  • $\begingroup$ 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. $\endgroup$
    – Godric Seer
    Jan 16 '13 at 19:59
2
$\begingroup$

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).

Improve this answer
$\endgroup$
2
$\begingroup$

To complete the enumeration of numerical methods, if the known roots really are good approximations for the roots of the next polynomial, then you can use Weierstraß' method, commonly known as Durand-Kerner method. It updates all roots at once, and is as fast as Newton's method, in fact is is Newton's method for a multidimensional problem.

Aberth-Ehrlich method is related, but it is not certain that the higher computational effort gives an overall faster convergence.

Additionally, since you are doing some kind of ray tracing, the univariate polynomials resulting from the intersection of a polynomial surface with a ray, you can divide the image into 6x6 or 8x8 blocks, compute the roots at a central pixel and then compute the gradients of the surface at the intersection points. You may get inspiration on how to do that by playing with the derivative code generator tapenade.

Using the point $x^*$ and the gradient $\nabla f(x^*)$ for the tangent plane equation $\langle \nabla f(x^*),\,x\rangle=\langle \nabla f(x^*),\,x^*\rangle$, you get better initial approximations by intersecting the rays $x_v(t)=\vec c+t\vec v$ of the other pixels first with the tangent plane.

Check for $\langle \nabla f(x^*),\,v\rangle\approx 0$, this occurs close to the visual contours. These are the critical cases in any approach since they have multiple or very close roots.

Improve this answer
$\endgroup$
1
  • $\begingroup$ Aberth 1973 says that Durand-Kerner is quadratically convergent and Aberth's method is cubically convergent. $\endgroup$
    – endolith
    Jul 14 '16 at 4:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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