Bairstow's method initial quadratic approximations

Question

Bairstow's root finding method needs very good initial approximations for the quadratic factors in order to converge.

I tried various constants, random numbers, fractions out of the trailing coefficient (-a1/a2, -a0/a2; by Lin?) to no avail.

Please, does anyone know of a good method for choosing the factors?

For example:

1*x^8 + 118*x^7 + 1*x^6 + 2*x^5 - 2*x^4 - 3*x^3 + 3*x^2 + 2*x + 1

It take 3x as much time to find the root with the initial approximations 0.1, 0.2 than it does with 0.2, 2.0.

Or:

1*x^8 - 36*x^7 + 546*x^6 - 4536*x^5 + 22449*x^4 - 67284*x^3 + 118124*x^2 - 109584*x + 40320

takes slightly longer (~50%) with 0.1, 1.2 than with 0.1, 0.1

---

Trying to use Cauchy's bound for the initial quadratic approximation:

R=0
for i in range(1,n+1):
    R=max(abs(a[i]/a[0]),R)
R=1+R
phi=2*pi*random()
x1=complex(R*cos(phi),R*sin(phi))
x2=complex(x1.real,-x1.imag)
r=-x1.real-x2.real
s=(x1*x2).real

Unfortunately, this does not really speed-up the converge.

Answer 1

As I've worked on the article and supplied the pictures I can tell that you really do not need that good approximations.

The most important initial step is to reduce the polynomial to even degree, as told in the article. After that, you can not do wrong, there should be almost global convergence. To be sure, the same as in Newton's method applies: If there is no noticeable sign of convergence after 10 steps, restart with a different initial point.

It is of course sensible to compute some outer root radius and choose the initial quadratic factor to have roots inside this radius.

---

See the javascript implementation in the source code of http://catc.ac.ir/mazlumi/jscodes/bairstow.php for a "naive" or "vanilla" but seemingly robust implementation. No reduction to even degree, no care for coefficient/root magnitudes,...

---

The example is efficiently an odd degree polynomial within the unit disk with one root -117.9917 at virtual infinity. For the initialization in every step one should compute the outer root radius, which is R=119 in the "1+max of abs of coefficients" (leading coefficient 1) version. Then initialize with x^2-R^2 or phi=2*pi*random(); and x^2+R^2*cos(phi)*x+R^2*sin(phi) or something similar.