Using Cauchy bound along with Graeffe's method for largest polynomial root bound

107 Views Asked by At

I'm applying Cauchy bound to the example from this http://www.numdam.org/item/M2AN_1990__24_6_693_0.pdf paper, and I'm having the results different from the written one. I may have some misunderstanding of either Graeffe's method or Cauchy bound https://en.wikipedia.org/wiki/Geometrical_properties_of_polynomial_roots

Cauchy bound:

Cauchy bound

The example from the paper: Example of applying Cauchy bound

So, applying Cauchy bound:
For P: max(a_k) = 6, Cauchy bound = 7
For P1: max(a_k) = 52, Cauchy bound = (52+1)^(1/2) = 7.28011
For P2: max(a_k) = 2951, Cauchy bound = (2951+1)^(1/4)= 7.37105
For P3: max(a_k) = 31 116 689, Cauchy bound = 31 116 689 ^ (1/8) = 8.6422 
For P4: max(a_k) = 36 706 467 938 304, Cauchy bound = 36 706 467 938 304 ^ (1/16) = 7.04363

I'm ignoring +1 in bound starting from P3. The power of root I'm using is due to Greaffe's method powering.

The result I have is different from what's written in the paper, so what I'm doing wrong?

1

There are 1 best solutions below

1
On BEST ANSWER

The root bounds $$ R=1+\max_{k<n}\left|\frac{a_k}{a_n}\right|~~\text{ and }~~ R=\max\left(1,\frac{|a_{n-1}|+...+|a_1|+|a_0|}{|a_n|}\right) $$ and similar are upper bounds for the first and only positive root of $$ q(R)=|a_n|R^n-|a_{n-1}|R^{n-1}-...-|a_1|R-|a_0|. $$ The bound called "Cauchy" and used in the article seems to be exactly this root.

Computing this root via Newton is relatively simple as there is only one sign variation in the coefficient sequence and $\tilde q(S)=-S^nq(1/S)$ is convex and increasing from the value $\tilde q(0)=-1$. The sequence of Newton iterates is thus falling from the second element on. One of the simple bounds gives a good starting value.

def graeffe(p):
    alt = np.ones_like(p); alt[1::2]=-1;
    p1 = np.convolve(p, alt*p)[::2]
    return p1

def root_bound(p):
    q = -abs(p); q[0]*=-1;
    R = np.roots(q)[0]
    return R.real

p = np.array((1,1,6,-5,3,0,2))

for k in range(1,6):
    p = graeffe(p)
    print(f"p[{k}] = ",p)
    R = root_bound(p)**(2**-k);
    print("R = ",R)

which prints

p[1] =  [ 1 11 52 15 33 12  4]
R =  3.825462972418143

p[2] =  [   1  -17 2440 2951 1145  120   16]
R =  2.7729087018852683

p[3] =  [       1     4591  6056224 -3116689   680865    22240      256]
R =  2.9451738688536895

p[4] =  [             1       -8964833 36706467938304 -1467012622369
   605308261633     -146014720          65536]
R =  2.7700750552861755

p[5] =  [                   1       -6955294841281  2131130497528204288
  6693824706615912383 -8910505650865396223    58018666012082176
           4294967296]
R =  2.5195464036912854

Rounding up the bound for the 4th iteration gives indeed $\rho<2.7701<2.771$.