Accelerated Approximation of the Complex Roots of a Univariate Polynomial (Extended Abstract)
aa r X i v : . [ c s . S C ] A p r Accelerated Approximation of the Complex Roots of aUnivariate Polynomial
Extended Abstract
Victor Y. Pan
Depts. of Mathematics and Computer ScienceLehman College and Graduate Centerof the City University of New YorkBronx, NY 10468 USA [email protected] http://comet.lehman.cuny.edu/vpan/
Elias P. Tsigaridas
INRIA, Paris-Rocquencourt Center,PolSys ProjectUPMC, Univ Paris 06, LIP6CNRS, UMR 7606, LIP6Paris, France [email protected]
Highly efficient and even nearly optimal algorithms havebeen developed for the classical problem of univariate poly-nomial root-finding (see, e.g., [6], [7], [4], and the bibliogra-phy therein), but this is still an area of active research. Bycombining some powerful techniques developed in this areawe devise new nearly optimal algorithms, whose substantialmerit is their simplicity, important for the implementation.We first recall the basic concept of the isolation ratio , cen-tral also for [6], [7]. Assume a real or complex polynomial p = p ( x ) = P di =0 p i x i = p n Q dj =1 ( x − z j ) , p d = 0 , of de-gree d , an annulus A ( X, R, r ) = { x : r ≤ | x − X | ≤ R } onthe complex plane with a center X and the radii r and R ofthe boundary circles. Then the internal disc D ( X, r ) = { x : | x − X | ≤ r } is R/r - isolated and R/r is its isolation ratio if the polynomial p has no roots in the annulus. Next wereproduce [13, Corollary 4.5]. It shows that Newton’s iter-ation converges quadratically to a single simple root of p ifis initiated at the center of a 5 d -isolated disc that containsjust this root. Theorem 1.
Suppose both discs D ( c, r ) and D ( c, r/s ) for s ≥ d contain a single simple root α of a polynomial p = p ( x ) of degree d . Then Newton’s iteration x k +1 = x k − p ( x k ) /p ′ ( x k ) , k = 0 , , . . . (1) converges quadratically to the root α right from the startprovided x = c . Now suppose that we are given a disc with a single simplezero of p having an isolation ratio 1 + η for a fixed constant η >
0. Can we increase the ratio to 5 d ? Yes, we just needto apply a technique already used in [15] for the computationof the power sums of the roots lying inside such a disc. Inour case this is a single root, the power sum is the root itself,and we just need its approximation c within an error at most Permission to make digital or hard copies of all or part of this work forpersonal or classroom use is granted without fee provided that copies arenot made or distributed for profit or commercial advantage and that copiesbear this notice and the full citation on the first page. To copy otherwise, torepublish, to post on servers or to redistribute to lists, requires prior specificpermission and/or a fee.Copyright 20XX ACM X-XXXXX-XX-X/XX/XX ...$15.00. ∆ such that rη/ ∆ ≥ d . Indeed in this case ∆ ≤ . rη/d ,and so the disc D ( c, ∆) is 5 d -isolated.We can shift and scale the variable x , and so wlog we as-sume dealing with a (1 + t ) -isolated disc D (0 , r ) for r =1 / (1 + t ) for a fixed t >
0, and with polynomial p having asingle simple root z in this disc. Recall the Laurent expan-sion, p ′ ( x ) p ( x ) = d X j =1 x − z j = − ∞ X k =1 S k x k − + ∞ X k =0 s k x − k − = ∞ X h = −∞ c h x h . (2)Here | x | = 1, s = 1 , s k = z k , S k = P di =2 z − ki , k = 1 , , . . . Consequently s k = z k , whereas S k is the k th power sum ofthe zeros of the reverse polynomial p rev ( x ) that lie in the disc D (0 , r ). The leftmost equation of (2) is verified by the dif-ferentiation of p ( x ) = p n Q dj =1 ( x − z j ). The middle equationis implied by the decompositions x − z = x P ∞ h =0 (cid:0) z x (cid:1) h and x − z i = − z i P ∞ h =0 (cid:16) xz i (cid:17) h for i >
1, provided | x | = 1 for all i . We cover the case of any positive integer k , although weonly need the case where k = 1. For a fixed positive integer q we compute the approximations s ∗ k ≈ s k as follows, s ∗ k = 1 q q − X j =0 ω j ( k +1) p ( ω j ) /p ′ ( ω j ) , k = 1 , , . . . , q − . (3)Here ω = ω q = exp(2 π √− /q ) is a primitive q th root ofunity. Then the evaluation of the polynomial p ( x ) at the q throots of unity amounts to the same task for a polynomial p q ( x ) of degree at most q − p q,i = P lj =0 p i + jq for l = ⌊ d/q ⌋ obtained by means of less than d additions of the coefficients of p .Having computed the polynomial p q ( x ) we reduce the eval-uation of all the desired approximations s ∗ k for k = 1 , . . . , q − q points,that is to a total of O ( q log( q )) ops. Namely, we apply twoDFTs to compute p ( ω i ) and p ′ ( ω i ) for i = 0 , , . . . , q − ω hi ] q − h,i =0 bythe vector v = [ p ( ω i ) /p ′ ( ω i ] q − i =0 .et us estimate the approximation errors. Equations (2)and (3) imply that s ∗ k = + ∞ X l = −∞ c − k − lq . Moreover, (2) for h = − k − , k ≥ s k = c − k − ,whereas (2) for h = k − , k ≥ S k = − c k − .Consequently s ∗ k − s k = ∞ X l =1 ( c lq − k − + c − lq − k − ) . We assumed in (3) that 0 < k < q −
1. It follows that c − lq − k − = s lq + k and c lq − k − = − S lq − k for l = 1 , , . . . , andwe obtain s ∗ k − s k = ∞ X l =1 ( s lq + k − S lq − k ) . (4)On the other hand | s h | ≤ z h , | S h | ≤ ( d − z h , h = 1 , , . . . where z = max ≤ j ≤ d min( | z j | , / | z j | ), and so z ≤ t in ourcase. Substitute these bounds into (4) and obtain | s ∗ k − s k | ≤ ( z q + k + ( d − z q − k ) / (1 − z q ). Therefore it is sufficient tochoose q of order log( d ) to decrease the error of the approx-imation to the root z by a factor of gd h for any pair ofconstants g and h , and so we can ensure the desired errorbound ∆. To support this computation we only need lessthan d additions, followed by O (log( d )) evaluations of thepolynomial p q ( x ) of degree q − l th roots of unityfor l = O (log( d )). This involves O (log( d ) log(log( d ))) opsoverall. (Here and hereafter “ops” stand for “arithmetic op-erations”.) Summarizing we obtain the following estimates. Theorem 2.
Suppose the unit disc D (0 , r ) = { x : | x | ≤ } is (1 + η ) -isolated for (1 + η ) r = 1 and a fixed η > andcontains a single simple root z of a polynomial p = p ( x ) of a degree d . Then it is sufficient to apply less than d additions and O (log( d ) log(log( d ))) other ops to compute a d -isolated subdisc of D (0 , r ) containing this root. Combine Thm. 1 and 2 and obtain the following result.
Corollary 3.
Under the assumptions of Theorem 2 we canapproximate the root z of the polynomial p ( x ) within a fixedpositive error bound ǫ < by using O (log( d ) log(log( d )) + d log(log(1 /ǫ ))) ops. Corollary 4.
Suppose that we are given d discs, each con-taining a single simple root of a polynomial p = p ( x ) ofdegree d and each being (1 + η ) -isolated for a fixed η > . Then we can approximate all d roots of this polyno-mial within a fixed positive error bound ǫ < by using O ( d log ( d )(1 + log(log(1 /ǫ )))) ops. Proof:
Apply the same algorithm that supports Corollary 3concurrently in all d given discs, but instead of the q th rootsof unity use q equally spaced points at the boundary circleof each input disc (that is dq = O ( d log d ) points overall)and instead of FFT apply the Moenck–Borodin algorithmfor multipoint polynomial evaluation. Also use it at thestage of performing concurrent Newton’s iteration initial-ized at the centers of the 5 d -isolated subdiscs of the d inputdiscs, each subdisc computed by the algorithm that supportsTheorem 2. Here we work with the d th degree polynomial p rather than with the q th degree polynomials p q becauseto support transition to polynomials p q of the degree q for d discs we would need to perform d shifts and scalings ofthe variable x . Instead we employ the Moenck–Borodin al-gorithm, which still enables us to obtain a nearly optimalroot-refiner. Technically, in a relatively minor change of ouralgorithm, we replace the matrix Ω = [ ω j ( k +1) ] j,k in (3) bythe matrix [ c + ω j ( k +1) ] j,k = c [1] j,k + Ω where c is invariantin j and k . The multiplication of the new matrix by a vec-tor v is still reduced to multiplication of the matrix Ω by avector v with the additional 3 d ops for computing the vector c [1] j,k v and adding it to the vector Ω v . (cid:3) The Moenck–Borodin algorithm uses nearly linear arith-metic time, and [2] proved that this algorithm supports mul-tipoint polynomial evaluation at a low Boolean cost as well(see also [14], [10], [3], [11], [8], [9]). Consequently our al-gorithm supporting Corollary 4 can be extended to supporta nearly optimal Boolean cost bound for refining all simpleisolated roots of a polynomial .We can immediately relax the assumption that the rootsare simple because our proof of Theorem 2 applies to a mul-tiple root as well. Furthermore deduce from the Lucas the-orem that the isolation ratio of the basic discs in our algo-rithms does not decrease when we shift from a polynomialto its derivative and higher order derivatives. Therefore wecan just apply Newton’s iteration to the derivative or to ahigher order derivative to approximate a double or multipleroot, respectively.
Boolean cost bounds
Hereafter e O B denotes the bit or Boolean complexity ignor-ing logarithmic factors. To estimate it we apply some resultsfrom [10]–[12], which hold in the general case where the co-efficients of the polynomials are known up to an arbitraryprecision. In our case the input polynomial is known ex-actly; the parameter λ to be specified in the sequel could beconsidered as the working precision.Let p be given as a λ -approximation, ie lg k p − e p k ∞ ≤− λ . We compute p q by using d additions. This produces apolynomial such that lg k p q k ∞ ≤ τ +lg d , and lg k p q − e p q k ∞ ≤− λ + τ lg d + 1 / d + 1 / d = O ( − λ + τ lg d + lg d ).Similar bounds hold for p ′ q , ie lg k p ′ q k ∞ ≤ τ + 2 lg d , andlg k p ′ q − e p ′ q k ∞ ≤ − λ + τ lg d + 3 / d + 1 / d = O ( − λ + τ lg d + lg d ).Recall that | p ′ q ( ω i ) | ≤ τ + 2 lg d + lg lg d + 2 and | p ′ q ( ω i ) − ^ p ′ q ( ω i ) | ≤ − λ + τ lg(2 d ) + 3 / d + 5 / d + lg lg d + 5 forall i , [11, Lemma 16], and similar bounds hold for p q ( ω i ).The divisions p q ( ω i ) /p ′ q ( ω i ) output complex numbers suchthat | p q ( ω ) /p ′ q ( ω ) | ≤ τ +2 lg d +lg lg d +2 with the logarithmof the error ≤ − λ + τ lg(4 d )+3 / d +9 / d +2 lg lg d +11.The final DFT produces numbers such that the logarithmsof their magnitudes are not greater than τ +2 lg d +2 lg lg d +4 and the logarithms of their approximation errors are atmost − λ + τ lg(8 d ) + 3 / d + 13 / d + 4 lg lg d + 18, [11,Lemma 16].To achieve an error within 2 − ℓ in the final result, we per-form all the computations with accuracy λ = ℓ + τ lg(8 d ) +3 / d + 13 / d + 4 lg lg d + 18, that is ℓ = O ( ℓ + τ lg d +lg d ) = e O ( ℓ + τ ).We perform d additions at the cost O B ( dλ ) and performthe rest of computations, that is the 3 DFTs, at the cost O B (lg d lg lg d µ ( λ )) or e O B ( d ( ℓ + τ )) [11, Lemma 16].f the root that we want to refine is not in the unit disc,then we replace τ in our bounds with dτ .We apply a similar analysis from [10, Section 2.3] to theNewton iteration (see also [11, Section 2.3]) and arrive at thesame asymptotic bounds on the Boolean complexity. Onlythe overhead constants change because now we perform com-putations with complex numbers.The overall complexity is e O B ( d τ + dℓ ) and the workingprecision is O ( dτ + ℓ ).Here we assume the exact input, that is assume the co-efficients of the input polynomials known up to arbitraryprecision. For the refinement of the root up to precision of L bits, we arrive at an algorithm with the complexity in e O B ( d τ + dL ).If we are interested in refining all complex roots, we cannotwork anymore with the polynomial p q of degree q = O (lg d )unless we add the cost of d shifts of the initial approxi-mations to the origin. Instead we rely on fast algorithmsfor multipoint evaluation. Initially we evaluate the polyno-mial p of degree d at O ( d lg d ) points, and we assume thatlg k p k ∞ ≤ τ . These d points approximate the roots of p , andso their magnitude is at most ≤ τ .We use the following result of [12, Lemma 21]. Similarbounds appear in [2, 3, 14]. Lemma 5 (Modular representation).
Assume m +1 poly-nomials, F ∈ ( C[ x ] of degree mn and P j ∈ ( C[ x ] of degree n , for j = 1 , . . . , m such that k F k ∞ ≤ τ and all roots ofthe polynomials P j for all j have magnitude of at most ρ .Furthermore assume λ -approximations of F by e F and of P j by e P j such that k F − e F k ∞ ≤ − λ and k P j − e P j k ∞ ≤ − λ .Let ℓ = λ − O ( τ lg m + m n ρ ) . Then we can compute an ℓ -approximations e F j of F j = F mod P j for j = 1 , . . . , m such that k F j − e F j k ∞ ≤ − ℓ in e O B ( m n ( ℓ + τ + m n ρ )) . Using this lemma we bound the overall complexity of mul-tipoint evaluation by e O B ( d ( L + dτ )). The same bounds holdsat the stage where we perform Newton’s iteration. We needto apply Newton’s operator e O (1) for each root. Each ap-plication of the operators consists of two polynomial evalua-tions. We perform the evaluations simultaneously and applyLemma 5 to bound the complexity. On similar estimates forthe refinement of the real roots see [11]. Extensions
The algorithm of [5] computes at nearly optimal cost 64 d -isolated initial discs for all d roots of a polynomial p ( x ). Bycombining this algorithm with ours we obtain a distinct al-ternative algorithm, which like the one of [6], [7], supportsthe record nearly optimal bounds on the Boolean complexityof the approximation of all complex polynomial roots , but hasthe advantage of allowing substantially simpler implementa-tion.Finally the same algorithm of [15] approximate the powersums of any number m of roots (forming, e.g., a single clusteror a number of clusters) in an isolated disc. The algorithmruns at about the same cost, already stated and dependingjust on the isolation ratio. Having the power sums availablewe can readily compute the coefficients of the factor f of p ofdegree d f , whose roots are exactly the roots of p in this disc:this is a numerically stable algorithm using O ( m log( m )) ops(cf. [1, pages 34–35]). Acknowledgments.
VP is supported by NSF Grant CCF–1116736.ET is partially supported by GeoLMI (ANR 2011 BS03 011 06), HPAC (ANR ANR-11-BS02-013) and an FP7 Marie Curie Career IntegrationGrant.
1. REFERENCES [1] D. Bini and V. Pan.
Polynomial and Matrix Computations ,volume 1: Fundamental Algorithms. Birkh¨auser, Boston, 1994.[2] P. Kirrinnis, Polynomial Factorization and Partial FractionDecomposition by Simultaneous Newton’s Iteration,
J. ofComplexity , 378–444 (1998).[3] A. Kobel and M. Sagraloff, Fast Approximate PolynomialMultipoint Evaluation and Applications, arXiv:1304.8069v1[cs.NA] 30 April 2013.[4] J. M. McNamee. and V. Y. Pan, Numerical Methods for Rootsof Polynomials , Part 2 (XXII + 718 pages), Elsevier (2013).[5] K. Mehlhorn, M., Sagraloff, P. Wang, From ApproximateFactorization to Root Isolation with Application to CylindricalAlgebraic Decomposition, in
Proc. Intl. Symp. on Symbolicand Algebraic Computations (ISSAC) , Boston, (M. Kauers,editor), 283–290, ACM Press, New York (2013).[6] V. Y. Pan, Optimal (up to Polylog Factors) Sequential andParallel Algorithms for Approximating Complex PolynomialZeros,
Proc. 27th Ann. ACM Symp. on Theory of Computing(STOC ’95) , ACM Press, New York, 741–750 (1995).[7] V. Y. Pan, Univariate Polynomials: Nearly Optimal Algorithmsfor Factorization and Rootfinding,
Journal of SymbolicComputations , 33, 5, 701–733, 2002.[8] V. Y. Pan, Transformations of Matrix Structures Work Again,accepted by
Linear Algebra and Its Applications , available atarXiv:1311.3729v1 [math.NA] 15 Nov 2013.[9] V. Y. Pan, Fast Approximation Algorithms for Computationswith Cauchy Matrices and Extensions, Tech. ReportTR-2014005,
PhD Program in Comp. Sci. , Graduate Center,CUNY , 2014Available at http://tr.cs.gc.cuny.edu/tr/techreport.php?id=469Proc. version in Proceedings of CSR 2014 (E.A. Hirsch et al.(Eds.)), LNCS 8476, pp. 287-300, 2014 (Springer InternationalPublishing, Switzerland 2014).[10] V. Y. Pan and E. P. Tsigaridas, in
Proc. International Symp.on Symbolic and Algebraic Computations (ISSAC 2013) ,Boston, Massachusetts, June 2013, (M. Kauers, editor),299–306, ACM Press, New York (2013).[11] V. Y. Pan and E. P. Tsigaridas, Nearly Optimal Refinement ofReal Roots of a Univariate Polynomial, Tech. Report, INRIA(2013). url: http://hal.inria.fr/hal-00960896,[12] V. Y. Pan and E. P. Tsigaridas, Nearly Optimal Computationswith Structured Matrices, Tech. Report, INRIA (2014)(Submitted).[13] J. Renegar, On the worst-case arithmetic complexity ofapproximating zeros of polynomials,
J. of Complexity
3, 23, 2