Graeffe Iteration Research Articles

A near-optimal subdivision algorithm for complex root isolation based on the Pellet test and Newton iteration

We describe a subdivision algorithm for isolating the complex roots of a polynomial F∈C[x]. Given an oracle that provides approximations of each of the coefficients of F to any absolute error bound and given an arbitrary square B in the complex plane containing only simple roots of F, our algorithm returns disjoint isolating disks for the roots of F in B.Our complexity analysis bounds the absolute error to which the coefficients of F have to be provided, the total number of iterations, and the overall bit complexity. It further shows that the complexity of our algorithm is controlled by the geometry of the roots in a near neighborhood of the input square B, namely, the number of roots, their absolute values and pairwise distances. The number of subdivision steps is near-optimal. For the benchmark problem, namely, to isolate all the roots of a polynomial of degree n with integer coefficients of bit size less than τ, our algorithm needs O˜(n3+n2τ) bit operations, which is comparable to the record bound of Pan (2002). It is the first time that such a bound has been achieved using subdivision methods, and independent of divide-and-conquer techniques such as Schönhage's splitting circle technique.Our algorithm uses the quadtree construction of Weyl (1924) with two key ingredients: using Pellet's Theorem (1881) combined with Graeffe iteration, we derive a “soft-test” to count the number of roots in a disk. Using Schröder's modified Newton operator combined with bisection, in a form inspired by the quadratic interval method from Abbot (2006), we achieve quadratic convergence towards root clusters. Relative to the divide-conquer algorithms, our algorithm is quite simple with the potential of being practical. This paper is self-contained: we provide pseudo-code for all subroutines used by our algorithm.

Numerical analysis of a bisection-exclusion method to find zeros of univariate analytic functions

We state precise results on the complexity of a classical bisection-exclusion method to locate zeros of univariate analytic functions contained in a square. The output of this algorithm is a list of squares containing all the zeros. It is also a robust method to locate clusters of zeros. We show that the global complexity depends on the following quantities: the size of the square, the desired precision, the number of clusters of zeros in the square, the distance between the clusters and the global behavior of the analytic function and its derivatives. We also prove that, closed to a cluster of zeros, the complexity depends only on the number of zeros inside the cluster. In particular, for a polynomial which has d simple roots separated by a distance greater than sep, we will prove the bisection-exclusion algorithm needs O ( d 3 log ( d / sep ) ) tests to isolate the d roots and the number of squares suspected to contain a zero is bounded by 4 d . Moreover, always in the polynomial case, we will see the arithmetic complexity can be reduced to O ( d 2 ( log d ) 2 log ( d / sep ) ) using ⌈ log d ⌉ steps of the Graeffe iteration.

Effective Fast Algorithms for Polynomial Spectral Factorization

Let p(z) be a polynomial of degree n having zeros |ξ1|≤⋅⋅⋅≤|ξ m |<1<|ξ m+1|≤⋅⋅⋅≤|ξ n |. This paper is concerned with the problem of efficiently computing the coefficients of the factors u(z)=∏ i=1 m (z−ξ i ) and l(z)=∏ i=m+1 n (z−ξ i ) of p(z) such that a(z)=z −m p(z)=(z −m u(z))l(z) is the spectral factorization of a(z). To perform this task the following two-stage approach is considered: first we approximate the central coefficients x −n+1,. . .x n−1 of the Laurent series x(z)=∑ i=−∞ +∞ x i z i satisfying x(z)a(z)=1; then we determine the entries in the first column and in the first row of the inverse of the Toeplitz matrix T=(x i−j ) i,j=−n+1,n−1 which provide the sought coefficients of u(z) and l(z). Two different algorithms are analyzed for the reciprocation of Laurent polynomials. One algorithm makes use of Graeffe's iteration which is quadratically convergent. Differently, the second algorithm directly employs evaluation/interpolation techniques at the roots of 1 and it is linearly convergent only. Algorithmic issues and numerical experiments are discussed.

Computations with infinite Toeplitz matrices and polynomials

We relate polynomial computations with operations involving infinite band Toeplitz matrices and show applications to the numerical solution of Markov chains, of nonlinear matrix equations, to spectral factorizations and to the solution of finite Toeplitz systems. In particular two matrix versions of Graeffe's iteration are introduced and their convergence properties are analyzed. Correlations between Graeffe's iteration for matrix polynomials and cyclic reduction for block Toeplitz matrices are pointed out. The paper contains a systematic treatment of known topics and presentation of new results, improvements and extensions.

Solving matrix polynomial equations arising in queueing problems

The matrix equation ∑i=0nAiXi=0, where the Ai's are m×m matrices, is encountered in the numerical solution of Markov chains which model queueing problems. We provide here a unifying framework in terms of Möbius' mapping to relate different resolution algorithms having a quadratic convergence. This allows us to compare algorithms like logarithmic reduction (LR) and cyclic reduction (CR), which extend Graeffe's iteration to matrix polynomials, and the invariant subspace (IS) approach, which extends Cardinal's algorithm. We devise new iterative techniques having quadratic convergence and present numerical experiments.

On the Geometry of Graeffe Iteration

A new version of Graeffe's algorithm for finding all the roots of univariate complex polynomials is proposed. It is obtained from the classical algorithm by a process analogous to renormalization of dynamical systems. This iteration is called the renormalized Graeffe iteration. It is globally convergent, with probability 1. All quantities involved in the computation are bounded once the initial polynomial is given (with probability 1). This implies remarkable stability properties for the new algorithm, thus overcoming known limitations of the classical Graeffe algorithm. If we start with a degree-d polynomial, each renormalized Graeffe iteration costs O(d2) arithmetic operations, with memory O(d). A probabilistic global complexity bound is given. The case of univariate real polynomials is briefly discussed. A numerical implementation of the algorithm presented herein allows us to solve random polynomials of degree up to 1000.

Condition Number Bounds for Problems with Integer Coefficients

An explicit a priori bound for the condition number associated to each of the following problems is given: general linear equation solving, least squares, nonsymmetric eigenvalue problems, solving univariate polynomials, and solving systems of multivariate polynomials. It is assumed that the input has integer coefficients and is not on the degeneracy locus of the respective problem (i.e., the condition number is finite). Our bounds are stated in terms of the dimension and of the bit-size of the input. In the same setting, bounds are given for the speed of convergence of the following iterative algorithms: QR iteration without shift for the symmetric eigenvalue problem and Graeffe iteration for univariate polynomials.

Tangent Graeffe iteration

Graeffe iteration was the choice algorithm for solving univariate polynomials in the XIX-th and early XX-th century. In this paper, a new variation of Graeffe iteration is given, suitable to IEEE floating-point arithmetics of modern digital computers. We prove that under a certain generic assumption the proposed algorithm converges. We also estimate the error after N iterations and the running cost. The main ideas from which this algorithm is built are: classical Graeffe iteration and Newton Diagrams, changes of scale (renormalization), and replacement of a difference technique by a differentiation one. The algorithm was implemented successfully and a number of numerical experiments are displayed.