Algorithm for Approximating Complex Polynomial Zeros

Algorithm for Approximating Complex Polynomial Zeros

Victor Pan

Lehman College, CUNY

Algorithms Seminar

June 8, 1998

[summary by Bruno Salvy]

A properly typeset version of this document is available in postscript and in pdf.

If some fonts do not look right on your screen, this might be fixed by configuring your browser (see the documentation here).

Abstract

An algorithm for approximating complex polynomial zeros is presented. Its complexity is optimal up to polylogarithmic factors and holds the current record.

Finding roots of a complex polynomial numerically in a guaranteed way with a fixed prescribed accuracy is difficult when no approximation is known in advance. This task cannot be performed in a fixed precision environment and implementations in computer algebra systems (where arbitrary precision is available) are seldom able to treat polynomials of degree a few hundreds. However, polynomials of very high degree arise frequently when solving a polynomial system by elimination. The work summarized here provides an algorithm supporting the following theorem.

Theorem 1 Let p(x) be a monic polynomial of degree n and z₁,...,z_n its zeros, with |z_i|£1, i=1,...,n. For a fixed positive b, approximations z_i^* satisfying

|z_i-z

|<2^-b, i=1,...,n (1)

can be computed at a cost bounded by O^~(n) arithmetic operations and O^~(n²(b+n)) boolean operations. The notation O^~ means that factors log n, log b or smaller are neglected.

Much more precise statements, proofs and parallel complexity estimates can be found in [5] and a pedagogical introduction to this area is [6].

The statement of the theorem can be modified to accommodate polynomials which are not monic (by first scaling the coefficients) or with roots of modulus larger than 1 by computing a bound on the moduli (see below) and then scaling the polynomial.

1 Lower Bounds

It is clear that the arithmetical complexity O^~(n) is optimal, since n coefficients of the input polynomial have to be treated. The boolean complexity O^~(n²(b+n)) is optimal in the very frequent case n=O(b).

Actually, O(n²b) is even a lower bound for the computation of one root of polynomials of degree n. This bound follows from the high susceptibility of the roots of a polynomial with respect to the coefficients. For instance, the polynomial xⁿ-a with a small a>0 has for root a^1/n. If this root is of order 2^-b, changing a to 0 is a change of the nb-th bit of a coefficient that changes the b-th bit of the root. This reasoning extends to other coefficients: let p=O(n) and consider xⁿ-a x^p. Then again a change of a bit at position O(nb) modifies the b-th bit of the solution. Thus b bits of the solution depend on O(nb) bits of each of O(n) coefficients, whence the O(n²b) lower bound. This example also illustrates why clusters of zeros defeat many numerical algorithms.

2 Outline of the Algorithm

The algorithm is based on a splitting technique where the polynomial p is split into factors of degree k and n-k with k=a n, for some aÎ(1/2,r), r being fixed. Applying this process recursively, any polynomial can be completely factored in O(log n) steps.

The splitting itself is computed in 3 steps:

Find a ``splitting'' circle not ``too close'' to roots of p and containing a n of them;
Compute the polynomial vanishing at these a n roots;
Divide p by this polynomial to obtain the other factor.

Each of these steps has to be performed in O(n²b+n³) boolean operations to yield the theorem.

The factors p_k and p_n-k of p are computed numerically. The following two lemmas show how the precision with which they are required can be bounded by ensuring that e^* is sufficiently small in the following inequality:

||p(x)-p_k(x)p_n-k(x)||£e

||p(x)||, (2)

where ||q(x)|| denotes the sum of the moduli of the coefficients of a polynomial q.

Lemma 1 [8] If

½½
½½
½½
½½

p(x)-

i=1

(x-z

)

½½
½½
½½
½½

<e||p(x)||,

with -log₂e³ bn+n+2, the inequalities (1) are satisfied.

Lemma 2 [8] Let p(x), f₁(x),...,f_k(x) and f(x),g(x) be polynomials such that

||p(x)-f₁(x)··· f_k(x)||

£e

||p(x)||

(3)

||f₁(x)-f(x)g(x)||

£e_k||f₁(x)||,

(4)

then

||p(x)-f(x)g(x)f₂(x)··· f_k(x)||£e

k+1

||p(x)||

holds, provided

e_k£e

||p(x)||

i=1

||f_i(x)||

From these lemmas follows that it is sufficient to compute the splitting with e^*£e/(n2ⁿ) in (2), where e comes from Lemma 1.

The splitting circle method was introduced by Schönhage [8, 9]. We now review the algorithms used in steps 1 and 2, together with the recent progress due to Victor Pan.

3 Numerical Factorization

To simplify the notation, assume the unit circle is a splitting circle for the polynomial p(x). Let p_k(x) be the monic polynomial whose k roots are those roots of p lying inside the circle. The computation of p_k(x) relies on the following integral representation of the power sums s_j of its zeros:

s_j=

2ip

ó
õ

|z|=1

p'(z)

p(z)

z^j dz

This idea originates in [2] and was refined by [8] to produce error bounds, i.e., to bound Q such that the s_j's can be computed by the discretization

Q-1

q=0

w^(i+1)q

p'(w^q)

p(w^q)

The value of Q depends on a lower bound for |p(z)| on the unit circle, which in turns is related to a bound on the distance from this circle to the closest root of p, hence the need for a circle ``not too close'' to the roots in Step 1 of the algorithm.

Efficiency is attained at the price of quite technical developments [8]. If the closest root to the circle is at distance O(1/n), a value of Q of order O(n²) is used¹ and the corresponding p'(w^q) and p(w^q) are computed by a discrete Fourier transform. From there, the sums s_j^* for j=0,...,K are computed by DFT, K being the smallest power of 2 larger than k=s₀. An approximation of the factor p_k(x) can then be recovered efficiently by a variant of Newton-Hensel's lifting (see [1, p. 34]). Then the other factor is obtained by division. In order to reach the right level of complexity, it is necessary to compute only O(n) bits for these steps and then refine the factorization by another Newton like algorithm as follows. Starting from the approximate factorization

||p(x)-p_k⁽⁰⁾(x)p_n-k⁽⁰⁾(x)||£e,

where p_k⁽⁰⁾ has degree k, the aim is to find a refinement p_k⁽¹⁾=p_k⁽⁰⁾+q_k, p_n-k⁽¹⁾=p_n-k⁽⁰⁾+q_n-k with deg q_i<i, improving the error. Since

p-p_k⁽¹⁾p_n-k⁽¹⁾=(p-p_k⁽⁰⁾p_n-k⁽⁰⁾)-p_k⁽¹⁾p_n-k⁽⁰⁾ -p_k⁽⁰⁾p_n-k⁽¹⁾-p_k⁽¹⁾p_n-k⁽¹⁾,

the Newton iteration is obtained by satisfying

(p-p_k⁽⁰⁾p_n-k⁽⁰⁾)=p_k⁽¹⁾p_n-k⁽⁰⁾+p_k⁽⁰⁾p_n-k⁽¹⁾, (5)

which determines p_k⁽¹⁾ and p_n-k⁽¹⁾ uniquely. These polynomials could be found by Euclid's algorithm, but this is too expensive. Instead, one also computes an inverse q⁽ⁱ⁾ of p_n-k⁽ⁱ⁾ modulo p_k⁽ⁱ⁾ by a second, parallel, Newton iteration and then p_k⁽ⁱ⁾ is given by q⁽ⁱ⁾p=q⁽ⁱ⁾(p-p_k⁽ⁱ⁾p_n-k⁽ⁱ⁾)mod p_k⁽ⁱ⁾. A similar formula gives p_n-k⁽ⁱ⁾. Then the required precision is obtained after a few iteration at a cost bounded by O(nloge^*).

4 Finding Splitting Circles

The basic technique to find discs containing a known number of roots of a polynomial is the iteration of Graeffe's method (see [3]). Starting from p(x) of degree n, one performs the following iteration:

p_i+1(x²)=(-1)ⁿp_i(x)p_i(-x),

which transforms the polynomial p_i(x) into a polynomial p_i+1(x) whose roots are the squares of the roots of p_i(x). This process emphasizes the differences of moduli between the roots. The coefficients of these iterates are Newton sums from which precise information about the different moduli of the roots of the original polynomial can be recovered at a low cost. More precisely, one gets the following lemma.

Lemma 3 Let z₁,...,z_n be the roots of p(x), satisfying |z₁|£...£|z_n|£1. Given c>0 and d³0, it is possible to compute r₁,r₁,...,r_n,r_n such that r_k£|z_k|£r_k=(1+c/n^d)r_k, k=1,...,n with O^~(n) arithmetic operations.

This iteration is applied after having first shifted the origin to the center of gravity of the roots, which is given by the first two coefficients of the polynomial. When it follows from this computation that there is a k=a n, a in a fixed interval (r,1-r), with some r<1 such that |z_k+1|/|z_k|³1+c/n for some c fixed in advance, then this yields a splitting circle and the factoring algorithm of the previous section can be applied.

It is when no such circle can be found that progress has been made by Victor Pan recently. In this case, there is an annulus centered at 0 which contains most of the roots of the polynomial. Now the idea is to shift the origin to each of r'=2r_11n/12 and ir', and apply the same method. Then either a good splitting circle is found, or there is a small circle which is easily computed and contains the intersection of these three annuli, itself containing an important cluster of zeros (at least half of the zeros of p if c=1/100). In this case, the idea is that one of the zeros of a derivative of p of high order (for instance, one can take p^{(ë n/2û +1)}) is either the center of a good splitting circle or makes it possible to isolate a massive cluster of zeros, where more than half of the zeros of p are at distance less than the desired accuracy 2^-b. In both cases, the polynomial can then be factored numerically and the computation proceeds on those factors that do not correspond to a massive cluster. Many refinements are given in [5], in particular it is shown that it is not necessary to compute all the zeros of p^{(ë n/2û +1)}.

Conclusion

This summary is a very rough sketch of a very detailed study given in [5]. For practical polynomial solving, other algorithms are known to perform extremely well, but their complexity analysis has yet to be done.

The talk also mentioned extensions to the multivariate case, this is described in [4].

References

[1]: Bini (Dario) and Pan (Victor Y.). -- Polynomial and matrix computations. Vol. 1. -- Birkhäuser Boston Inc., Boston, MA, 1994, Progress in Theoretical Computer Science, xvi+415p. Fundamental algorithms.
[2]: Delves (L. M.) and Lyness (J. N.). -- A numerical method for locating the zeros of an analytic function. Mathematics of Computation, vol. 21, 1967, pp. 543--560.
[3]: Henrici (Peter). -- Applied and computational complex analysis. -- Wiley-Interscience, New York, 1974, Pure and Applied Mathematics, vol. Volume 1: Power series, integration, conformal mapping, location of zeros, xv+682p.
[4]: Mourrain (Bernard) and Pan (Victor Y.). -- Asymptotic acceleration of solving multivariate polynomial systems of equations. In Proceedings of the Thirtieth Annual ACM Symposium on Theory of Computing. pp. 488--496. -- ACM Press, 1998.
[5]: Pan (V. Y.). -- Optimal and nearly optimal algorithms for approximating polynomial zeros. Computers & Mathematics with Applications, vol. 31, n°12, 1996, pp. 97--138.
[6]: Pan (Victor). -- Solving polynomials with computers. American Scientist, vol. 86, 1998, pp. 62--69.
[7]: Pan (Victor Y.). -- Optimal (up to polylog factors) sequential and parallel algorithms for approximating complex polynomial zeros. In SIAM Journal on Computing. pp. 741--750. -- ACM Press, 1995.
[8]: Schönhage (Arnold). -- The fundamental theorem of algebra in terms of computational complexity. -- Technical report, Mathematisches Institut der Universität Tübingen, 1982. Preliminary report.
[9]: Schönhage (Arnold). -- Equation solving in terms of computational complexity. In Proceedings of the International Congress of Mathematicians, pp. 131--153. -- 1987. Berkeley, California, 1986.

1: More precise values are given in [8, p. 35].

This document was translated from L^AT_EX by H^EV^EA.