Laguerre’s Method for Polynomial Roots

This documentation was generated with the assistance of AI. Please report any inaccuracies.

Beyond degree three, no closed-form solution exists in general — a polynomial’s roots must be found iteratively. Laguerre’s method is the most straightforward general-purpose technique that reliably finds all the roots of a polynomial of any degree — real, complex, single, or multiple — converging to a root from essentially any starting point. It is described in Press et al., 2007, Section 9.5 (page 463), and implemented as LaguerrePolynomialRootsEstimator in the com.irurueta.numerical.roots package.

At normal scale, linear, quadratic, and cubic behavior near a root all look like one clean crossing; only high magnification reveals that the quadratic actually has two nearby roots and the cubic is a flat, ill-conditioned crossing

This is why root-finding for polynomials of degree three and up deserves its own method: at the scale a root is first located, a close pair of roots or a flattened multiple root is often indistinguishable from a single, well-conditioned one. Only after zooming in — exactly what iterative refinement and polishing do — does the true structure become visible.

Motivating the formula

Writing a degree- polynomial in terms of its roots, , the log-derivative and its derivative give two useful sums over all the roots:

Laguerre’s method makes a deliberately drastic simplifying assumption: that the one root being sought, , sits a distance from the current guess , while every other root sits at the same distance . Substituting and (for all other roots) into and gives two equations in and , which solve to:

where the sign is chosen to make the denominator’s magnitude as large as possible (for numerical stability), and the quantity inside the square root may be negative — so , and hence the root estimate, is computed using complex arithmetic even when searching for a real root. The method simply iterates until becomes negligibly small.

Despite this simplifying assumption looking crude, Laguerre’s method behaves remarkably well: for a polynomial with entirely real roots, it is guaranteed to converge to some root from any starting point. For complex roots, convergence is not rigorously proven, but empirically failures are rare — and, when they do happen, almost always take the form of a nonconverging limit cycle rather than outright divergence, which this package’s implementation specifically detects and escapes.

Breaking limit cycles

Every 10th iteration (MT), instead of the normal step, the implementation takes a step scaled by one of a fixed sequence of fractions (0, 0.5, 0.25, 0.75, 0.13, 0.38, 0.62, 0.88, 1.0) instead of the full Laguerre correction. This nudges the iteration out of the rare pathological cycles that would otherwise repeat forever — for example, a high-degree polynomial (degree ) whose roots are all just outside the unit circle, roughly equally spaced around it.

Finding all the roots: deflation and polishing

To find every root of a degree- polynomial, Laguerre’s method is applied repeatedly, each time to a deflated polynomial: once a root is found, the polynomial is divided by (synthetic division), producing a polynomial of one lower degree whose remaining roots are exactly the ones still being sought. Each Laguerre search starts from , which favors converging to whichever remaining root has the smallest magnitude.

Because each deflation step introduces a small amount of numerical error that compounds over successive deflations, this package — following Press et al., 2007 — treats the roots found this way as only tentative, then optionally polishes every one of them by running Laguerre’s method again, starting from each tentative root, but against the original, non-deflated polynomial. This removes most of the drift that deflation alone would otherwise leave behind.

flowchart TD A["Polynomial of degree n"] --> B["Run Laguerre's method from x=0\non the current (possibly deflated) polynomial"] B --> C{"every 10th iteration:\nuse a fractional step instead,\nto escape limit cycles"} C --> D["Converged: record this root"] D --> E["Deflate: divide the polynomial\nby (x - root)"] E --> F{"any roots remaining?"} F -->|"yes"| B F -->|"no, all n roots found"| G{"polishing enabled?"} G -->|"yes"| H["Re-run Laguerre's method on each\ntentative root against the ORIGINAL polynomial"] G -->|"no"| I["Sort roots by real part; done"] H --> I

Library implementation

LaguerrePolynomialRootsEstimator extends PolynomialRootsEstimator, and — unlike Closed-Form Polynomial Roots — is the only polynomial root finder in this package that accepts complex coefficients (a Complex[] array, constant term first), matching the fact that Laguerre’s method itself always uses complex arithmetic internally, even when every root turns out to be real.

  • setPolishRootsEnabled(boolean) / areRootsPolished() — whether to re-run Laguerre’s method on the original polynomial after deflation; enabled (DEFAULT_POLISH_ROOTS = true) by default.

Up to MAXIT = MT * MR = 8000 iterations are allowed per root (MT = 100 steps between fractional-step attempts, MR = 80 distinct fractional values available). After estimate() completes, the inherited getRoots() returns a Complex[] of length equal to the polynomial’s degree, sorted by real part; any root whose imaginary part is negligible relative to its real part is snapped to exactly zero.

Press et al., 2007 warns that some polynomials are exceedingly ill-conditioned — tiny changes in the coefficients can send the roots sprawling across the complex plane (Wilkinson’s classic example is one such case). This is a property of the polynomial itself, not a limitation of Laguerre’s method specifically; no root-finding algorithm can be expected to do much better in that situation.

Example

// x^3 - 6x^2 + 11x - 6 = (x-1)(x-2)(x-3), coefficients [a0, a1, a2, a3]
final var coefficients = new Complex[]{
        new Complex(-6.0), new Complex(11.0), new Complex(-6.0), new Complex(1.0)
};

final var estimator = new LaguerrePolynomialRootsEstimator(coefficients);
estimator.estimate();

final var roots = estimator.getRoots(); // Complex[]{1+0i, 2+0i, 3+0i}