Evaluation of Chebyshev polynomials on intervals and application to root finding
EEvaluation of Chebyshev polynomials on intervalsand application to root finding
Viviane Ledoux , and Guillaume Moroz INRIA, France
[email protected] École Normale Supérieure
[email protected] . Abstract.
In approximation theory, it is standard to approximate func-tions by polynomials expressed in the Chebyshev basis. Evaluating apolynomial f of degree n given in the Chebyshev basis can be donein O ( n ) arithmetic operations using the Clenshaw algorithm. Unfortu-nately, the evaluation of f on an interval I using the Clenshaw algorithmwith interval arithmetic returns an interval of width exponential in n . Wedescribe a variant of the Clenshaw algorithm based on ball arithmeticthat returns an interval of width quadratic in n for an interval of smallenough width. As an application, our variant of the Clenshaw algorithmcan be used to design an efficient root finding algorithm. Keywords:
Clenshaw algorithm · Chebyshev polynomials · Root finding · Ball arithmetic · Interval arithmetic.
Clenshaw showed in 1955 that any polynomial given in the form p ( x ) = n (cid:88) i =0 a i T i ( x ) (1)can be evaluated on a value x with a single loop using the following functionsdefined by recurrence: u k ( x ) = if k = n + 1 a n if k = n xu k +1 ( x ) − u k +2 ( x ) + a k if ≤ k < nxu ( x ) − u ( x ) + a if k = 0 (2)such that p ( x ) = u ( x ) .Unfortunately, if we use Equation (2) with interval arithmetic directly, theresult can be an interval of size exponentially larger than the input, as illustratedin Example 1. Example 1.
Let ε > be a positive real number, and let x be the interval [ − ε, + ε ] of width ε . Assuming that a n = 1 , we can see that u n − is an a r X i v : . [ c s . S C ] D ec nterval of width ε . Then by recurrence, we observe that u n − k is an interval ofwidth at least εF k where ( F n ) n ∈ N denotes the Fibonacci sequence, even if all a i = 0 for i < n .Note that the constant below the exponent is even higher when x is closer to . These numerical instabilities also appear with floating point arithmetic near and − as analyzed in [4].To work around the numerical instabilities near and − , Reinsch suggesteda variant of the Clenshaw algorithm [4,7]. Let d n ( x ) = a n and u n ( x ) = a n , andfor k between and n − , define d k ( x ) and u k ( x ) by recurrence as follows: (cid:40) d k ( x ) = 2( x − u k +1 ( x ) + d k +1 ( x ) + a k u k ( x ) = d k ( x ) + u k +1 (3)Computing p ( x ) with this recurrence is numerically more stable near . However,this algorithm does not solve the problem of exponential growth illustrated inExample 1.Our first main result is a generalization of Equation 3 for any value in theinterval [ − , . This leads to Algorithm 1 that returns intervals with tighterradii, as analyzed in Lemma 2. Our second main result is the use of classicalbackward error analysis to derive Algorithm 2 which gives an even better radii.Then in Section 3 we use the new evaluation algorithm to design a root solverfor Chebyshev series, detailed in Algorithm 3. In this section we assume that we want to evaluate a Chebyshev polynomial onthe interval I . Let a be the center of I and r be its radius. Furthermore, let γ and γ be the conjugate complex roots of the equation: X − aX + 1 = 0 . (4)In particular, using Vieta’s formulas that relate the coefficients to the roots of apolynomial, γ satisfies γ + γ = 2 a and γγ = 1 .Let z n ( x ) = a n and u n ( x ) = a n , and for k between and n − , define z k ( x ) and u k ( x ) by recurrence as follows: (cid:40) z k ( x ) = 2( x − a ) u k +1 ( x ) + γz k +1 ( x ) + a k u k ( x ) = z k ( x ) + γu k +1 ( x ) (5)Using Equation (4), we can check that the u k satisfies the recurrence relation u k ( x ) = 2 xu k +1 ( x ) − u k +2 ( x ) + a k , such that p ( x ) = xu ( x ) − u ( x ) + a .Let ( e k ) and ( f k ) be two sequences of positive real numbers. Let B R ( a, r ) and B R ( u k ( a ) , e k ) represent the intervals [ a − r, a + r ] and [ u k ( a ) − e k , u k ( a ) + e k ] .Let B C ( z k ( a ) , f k ) be the complex ball of center z k ( a ) and radius f k .ur goal is to compute recurrence formulas on the e k and the f k such that: (cid:40) z k ( B R ( a, r )) ⊂ B C ( z k ( a ) , f k ) u k ( B R ( a, r )) ⊂ B R ( u k ( a ) , e k ) . (6) Lemma 1.
Let e n = 0 and f n = 0 and for n > k ≥ : (cid:40) f k = 2 r | u k +1 ( a ) | + 2 re k +1 + f k +1 e k = min( e k +1 + f k , f k √ − a ) if | a | < else e k +1 + f k (7) Then, ( e k ) and ( f k ) satisfy Equation (6) .Proof (sketch). For the inclusion z k ( B R ( a, r )) ⊂ B C ( z k ( a ) , f k ) , note that γ hasmodulus , such that the radius of γz k +1 is the same as the radius of z k +1 whenusing ball arithmetics. The remaining terms bounding the radius of z k followfrom the standard rules of interval arithmetics.For the inclusion u k ( B R ( a, r )) ⊂ B C ( u k ( a ) , e k ) , note that the error segmenton u k is included in the Minkowski sum of a disk of radius f k and a segment ofradius e k +1 , denoted by M . If θ is the angle of the segment with the horizontal,we have cos θ = a . We conclude that the intersection of M with a horizontal lineis a segment of radius at most min( e k +1 + f k , f k √ − a ) . Corollary 1.
Let B R ( u, e ) = BallClenshawForward (( a , . . . , a n ) , a, r ) be theresult of Algorithm 1, then p ( B R ( a, r )) ⊂ B R ( u, e ) Moreover, the following lemma bounds the radius of the ball returned byAlgorithm 1.
Lemma 2.
Let B R ( u, e ) = BallClenshawForward (( a , . . . , a n ) , a, r ) be the re-sult of Algorithm 1, and let M be an upper bound on | u k ( a ) | for ≤ k ≤ n .Assume that ε k < M r for ≤ k ≤ n , then e < M n r if n < √ − a e < M n r √ − a if √ − a ≤ n < √ − a r e < M (cid:104) (1 + r √ − a ) n − (cid:105) if √ − a r < n Proof (sketch).
We distinguish cases. First if n < √ − a , we focus on therelation e k ≤ e k +1 + f k + M r , and we prove by descending recurrence that e k ≤ M ( n − k ) r and f k ≤ M r (2( n − k −
1) + 1) .For the case √ − a ≤ n , we use the relation e k ≤ f k √ − a + M r , that we substi-tute in the recurrence relation defining f k to get f k ≤ rM + r √ − a f k +1 + f k +1 + M r √ − a . We can check by recurrence that f k ≤ M √ − a (cid:104) (1 + r √ − a ) n − (cid:105) ,which allows us to conclude for the case √ − a r ≤ n . Finally, when √ − a ≤ n < √ − a r , we observe that (1 + r √ − a ) n − ≤ n exp(1) r √ − a which leads tothe bound for the last case. lgorithm 1 Clenshaw evaluation algorithm, forward error function
BallClenshawForward ( ( a , . . . , a n ) , a , r ) (cid:46) Computation of the centers u k u n +1 ← u n ← a n for k in n − , n − , . . . , do u k ← au k +1 − u k +2 + a k ε k ← bound on the rounding error for u k u ← au − u + a ε ← bound on the rounding error for u (cid:46) Computation of the radii e k f n ← e n ← for k in n − , n − , . . . , do f k ← r | u k +1 | + 2 re k +1 + f k +1 e k ← min ( e k +1 + f k , f k √ − a ) + ε k f ← r | u | + 2 re + f e ← min ( e + f , f √ − a ) + ε return B R ( u , e ) In the literature, we can find an error analysis of the Clenshaw algorithm [3].The main idea is to add the errors appearing at each step of the Clenshawalgorithm to the input coefficients. Thus the approximate result correspond tothe exact result of an approximate input. Finally, the error bound is obtainedas the evaluation of a Chebyshev polynomial. This error analysis can be useddirectly to derive an algorithm to evaluate a polynomial in the Chebyshev basison an interval in Algorithm 2.
Lemma 3.
Let e n = 0 and for n > k ≥ : e k = 2 r | u k +1 ( a ) | + e k +1 (8) and e = r | u ( a ) | + e . Then ( e k ) satisfies u k ( B R ( a, r )) ⊂ B R ( u k ( a ) , e k ) .Proof (sketch). In the case where the computations are performed without er-rors, D. Elliott [3, Equation (4.9)] showed that for γ = ˜ x − x we have: p (˜ x ) − p ( x ) = 2 γ n (cid:88) i =0 u i (˜ x ) T i ( x ) − γu (˜ x ) In the case where ˜ x = a and x ∈ B R ( a, r ) we have γ ≤ r and | T ( x ) | ≤ which implies e k ≤ r | u ( a ) | + (cid:80) ni =2 r | u i ( a ) | . orollary 2. Let B R ( u, e ) = BallClenshawBackward (( a , . . . , a n ) , a, r ) be theresult of Algorithm 2, and let M be an upper bound on | u k ( a ) | for ≤ k ≤ n .Assume that ε k < M r for ≤ k ≤ n , then e < M nr . Algorithm 2
Clenshaw evaluation algorithm, backward error function
BallClenshawBackward ( ( a , . . . , a n ) , a , r ) (cid:46) Computation of the centers u k u n +1 ← u n ← a n for k in n − , n − , . . . , do u k ← au k +1 − u k +2 + a k ε k ← bound on the rounding error for u k u ← au − u + a ε ← bound on the rounding error for u (cid:46) Computation of the radii e k e n ← for k in n − , n − , . . . , do e k ← e k +1 + 2 r | u k +1 | + ε k e ← e + r | u | + ε return B R ( u , e ) For classical polynomials, numerous solvers exist in the literature, such as thosedescribed in [5] for example. For polynomials in the Chebyshev basis, severalapproaches exist that reduce the problem to polynomial complex root finding[1], or complex eigenvalue computations [2] among other.In this section, we experiment a direct subdivision algorithm based on intervalevaluation, detailed in Algorithm 3. This algorithm is implemented and publiclyavailable in the software clenshaw [6].We applied this approach to Chebyshev polynomials whose coefficients are in-dependently and identically distributed with the normal distribution with mean and variance .As illustrated in Figure 1 our code performs significantly better than theclassical companion matrix approach. In particular, we could solve polynomialsof degree in the Chebyshev basis in less than seconds and polynomials ofdegree in . seconds on a quad-core Intel(R) i7-8650U cpu at 1.9GHz. Forcomparison, the standard numpy function chebroots took more than secondsfor polynomials of degree . Moreover, using least square fitting on the tenlast values, we observe that our approach has an experimental complexity closerto Θ ( n . ) , whereas the companion matrix approach has a complexity closer to Θ ( n . ) . lgorithm 3 Subdivision algorithm for root finding
Require: ( a , . . . , a n ) represents the Chebyshev polynomial approximating f ( x )( b , . . . , b n ) represents the Chebyshev polynomial approximating dfdx ( x ) Ensure:
Res is a list of isolating intervals for the roots of f in [ − , function SubdivideClenshaw ( ( a , . . . , a n ) , ( b , . . . , b n ) ) (cid:46) Partition [ − , in intervals where F either has constant sign or is monotonous L ← [ B R (0 , P artition ← [] while L is not empty do B R ( a, r ) ← pop the first element of L B R ( f, s ) ← BallClenshaw (( a , . . . , a n ) , a, r ) B R ( df, t ) ← BallClenshaw (( b , . . . , b n ) , a, r ) if f − s > then append the pair ( B R ( a, r ) , ” plus ”) to P artition else if f + s < then append the pair ( B R ( a, r ) , ” minus ”) to P artition else if g − s > or g + s < then append the pair ( B R ( a, r ) , ” monotonous ”) to P artition else B , B ← subdivide B R ( a, r ) append B , B to L(cid:46)
Compute the sign of F at the boundaries B R ( f, s ) ← BallClenshaw (( a , . . . , a n ) , − , append the pair ( B R ( − , , sign ( B R ( f, s ))) to P artition B R ( f, s ) ← BallClenshaw (( a , . . . , a n ) , , append the pair ( B R (1 , , sign ( B R ( f, s ))) to P artition(cid:46)
Recover the root isolating intervals
P artition ← sort P artitionRes ← the "monotonous" intervals of P artition such that the adjacent intervals have opposite signs return
Res . . . . . . . . . Degree of the input Chebyshev polynomial (in log deg) − − − − S o l v i n g t i m e i n ( i n l og s ec o nd s ) clenshaw (bisection and clenshaw evaluation)numpy (eigenvalues of the companion matrix) Fig. 1.
Time for isolating the roots of a random Chebyshev polynomial, on a quad-coreIntel(R) i7-8650U cpu at 1.9GHz, with 16G of ram
References
1. Boyd, J.: Computing zeros on a real interval through chebyshev expansion and poly-nomial rootfinding. SIAM Journal on Numerical Analysis (5), 1666–1682 (2002).https://doi.org/10.1137/S00361429013983252. Boyd, J.: Finding the zeros of a univariate equation: Proxy rootfinders, chebyshevinterpolation, and the companion matrix. SIAM Review (2), 375–396 (2013).https://doi.org/10.1137/1108382973. Elliott, D.: Error analysis of an algorithm for summing certain finite se-ries. Journal of the Australian Mathematical Society (2), 213–221 (1968).https://doi.org/10.1017/S14467887000052674. Gentleman, W.M.: An error analysis of Goertzel’s (Watt’s) method for com-puting Fourier coefficients. The Computer Journal (2), 160–164 (01 1969).https://doi.org/10.1093/comjnl/12.2.1605. Kobel, A., Rouillier, F., Sagraloff, M.: Computing real roots of real polynomials... and now for real! In: Proceedings of the ACM on International Symposium onSymbolic and Algebraic Computation. pp. 303–310. ISSAC ’16, ACM, New York,NY, USA (2016). https://doi.org/10.1145/2930889.29309376. Moroz, G.: Clenshaw 0.1 (Dec 2019). https://doi.org/10.5281/zenodo.3571248, https://gitlab.inria.fr/gmoro/clenshaw
7. Oliver, J.: An Error Analysis of the Modified Clenshaw Method for EvaluatingChebyshev and Fourier Series. IMA Journal of Applied Mathematics20