Gonnet, Pedro

Trefethen, Lloyd N.

Webb, Marcus

English

The barycentric interpolation formula defines a stable algorithm for evaluation atpoints in [−1, 1] of polynomial interpolants through data on Chebyshev grids. Here it is shown that for evaluation at points in the complex plane outside [−1, 1], the algorithm becomes unstable and should be replaced by the alternative modified Lagrange or “first barycentric” formula dating to Jacobi in 1825. This difference in stability confirms the theory published by N. J. Higham in 2004 [IMA J. Numer. Anal., 24 (2004), pp. 547–556] and has practical consequences for computation with rational functions

The University of Manchester - Institutional Repository

The University of Manchester ResearchStability of Barycentric interpolation formulas forextrapolationDOI:10.1137/110848797Document VersionAccepted author manuscriptLink to publication record in Manchester Research ExplorerCitation for published version (APA):Webb, M., Trefethen, L. N., & Gonnet, P. (2012). Stability of Barycentric interpolation formulas for extrapolation.SIAM Journal on Scientific Computing, 34(6), A3009-A3015. https://doi.org/10.1137/110848797Published in:SIAM Journal on Scientific ComputingCiting this paperPlease note that where the full-text provided on Manchester Research Explorer is the Author Accepted Manuscriptor Proof version this may differ from the final Published version. If citing, it is advised that you check and use thepublisher's definitive version.General rightsCopyright and moral rights for the publications made accessible in the Research Explorer are retained by theauthors and/or other copyright owners and it is a condition of accessing publications that users recognise andabide by the legal requirements associated with these rights.Takedown policyIf you believe that this document breaches copyright please refer to the University of Manchester’s TakedownProcedures [http://man.ac.uk/04Y6Bo] or contact uml.scholarlycommunications@manchester.ac.uk providingrelevant details, so we can investigate your claim.Download date:09. Jun. 2022STABILITY OF BARYCENTRIC INTERPOLATION FORMULASMARCUS WEBB∗, LLOYD N. TREFETHEN† , AND PEDRO GONNET‡Abstract. The barycentric interpolation formula defines a stable algorithm for evaluation atpoints in [−1, 1] of polynomial interpolants through data on Chebyshev grids. Here it is shownthat for evaluation at points in the complex plane outside [−1, 1], the algorithm becomes unstableand should be replaced by the alternative modified Lagrange or “first barycentric” formula datingto Jacobi in 1825. This difference in stability confirms the theory published by N. J. Higham in2004 (IMA J. Numer. Anal., v. 24) and has practical consequences for computation with rationalfunctions.Key words. Barycentric interpolation, Chebfun, rational approximation, Bernstein ellipse,Chebfun ellipseAMS subject classifications. 65D05, 41A201. The effect. In experiments involving numerical analytic continuation in thecomplex plane, we encountered an instability that surprised us [18]. On investigationit proved to have an elegant explanation, and the effect is worth knowing about forthose who compute with polynomials or rational functions. It has led to a significantcorrection in Chebfun [17].  −2 −1 0 1 2−2−1012−30−20−100102030  −2 −1 0 1 2−2−1012−30−20−100102030first barycentric formula (2.1) second barycentric formula (2.6)Fig. 1.1. Numerically computed values Rep(x) in the complex x-plane, where p is thepolynomial interpolant in 42 Chebyshev points in [−1, 1] to f(x) = tanh(πx/2). Since the aimis to show both sign and magnitude of Rep(x), the precise quantity plotted is sign(Rep(x)) ·log(1 + |Rep(x)|). On the left, the computation uses the formula (2.1) and the picture iscorrect, with oscillatory behavior outside a certain ellipse enclosing [−1, 1]. On the right,the computation is based on the more usual formula (2.6) and the picture is entirely wrongoutside approximately the same ellipse, showing seemingly random values that fail to grow inmagnitude as |x| → ∞. The errors are significant inside the ellipse, too, though this is notapparent in the plot.Figure 1.1 illustrates the effect. The function f(x) = tanh(πx/2) has been ap-proximated to 16 digit accuracy on [−1, 1] by interpolation by a polynomial p(x) ofdegree n in n+1 Chebyshev points xj = cos(πj/n), 0 ≤ j ≤ n; the value that achieves∗Worcester College, Oxford OX1 2HB, UK (webb@maths.ox.ac.uk)†Oxford University Mathematical Institute, 24–29 St Giles’, Oxford OX1 3LB, UK(trefethen@maths.ox.ac.uk, http://people.maths.ox.ac.uk/trefethen)‡Oxford University Mathematical Institute, 24–29 St Giles’, Oxford OX1 3LB, UK(gonnet@maths.ox.ac.uk, http://people.comlab.ox.ac.uk/gonnet)12 M. WEBB, L. N. TREFETHEN AND P. GONNETthe prescribed accuracy is n = 41. The figure shows contour plots of the numericallycomputed function Rep(x) in the complex x-plane. The result labeled “first barycen-tric formula” is correct, but the result labeled “second barycentric formula” is entirelywrong outside an ellipse enclosing [−1, 1], and also quite inaccurate for values of xinside the ellipse. The ellipse in question is a Bernstein ellipse, that is, with foci ±1,the largest such ellipse inside of which f is analytic.2. Two variants of the barycentric formula. By the first barycentric for-mula, we mean the modified Lagrange interpolation formulap(x) = `(x)n∑j=0λjfjx − xj,(2.1)where `(x) is the node polynomial`(x) =n∏k=0(x − xk)(2.2)and the numbers λj are the barycentric weightsλj =1`′(xj).(2.3)For the Chebyshev points, the weights turn out to beλj =2n−1n(−1)j, 1 ≤ j ≤ n − 1,(2.4)and half these values for j = 0 and n. For derivations, see [1] or [16]. Formulas (2.1)–(2.3) originate with Jacobi in 1825, and (2.4) is due to Marcel Riesz in 1916 [6, 11].We are grateful to Folkmar Bornemann of the Technical University of Munich fordrawing our attention to these early papers.The second barycentric formula is obtained as follows. If we represent the constantfunction f(x) = 1 by (2.1), we get1 = `(x)n∑j=0λjx − xj,(2.5)and dividing (2.1) by (2.5) givesp(x) =n∑j=0λjfjx − xj/n∑j=0λjx − xj.(2.6)Equation (2.6) originates with Taylor and Dupuy in the 1940s, and the special caseof Chebyshev points was first treated by Salzer in 1972 [2, 13, 15]. The terms “first”and “second” come from Rutishauser [12]. Formula (2.6) has a special elegance aboutit, since the factor `(x) of (2.1) has dropped out, and this feature has practical con-sequences. That factor has size approximately 2−n for x ∈ [−1, 1], making (2.2)susceptible to underflow in floating point arithmetic for large n; similarly (2.4) showsthat the weights λj have size approximately 2n, leading to a risk of overflow when λjis defined by (2.4). Moreover, all these numbers change with nth powers if [−1, 1] isBARYCENTRIC INTERPOLATION FORMULAS 3rescaled to a general interval [a, b]; the number 2 arises for [−1, 1] because the log-arithmic capacity of this interval is 1/2. With the formula (2.6), however, we cancancel the common factor 2n−1/n from (2.4), taking the weights in both numeratorand denominator to be simply ±1 in the interior and ± 12at the endpoints. This makes(2.6) scale-invariant and circumvents all problems of underflow and overflow.3. Explanation of the instability. Both formulas (2.1) and (2.6) work beau-tifully for interpolating a function on a Chebyshev grid in [−1, 1] and evaluating theinterpolant at points in [−1, 1] (assuming the underflow/overflow problem of (2.1) isaddressed for large n, which can be done by reformulating it via logarithms). Thisnumerical stability of (2.6) has been emphasized over the years by Salzer, Rutishauser,and other authors, including Henrici [4], and it is relied upon by Chebfun for poly-nomial interpolants even in millions of points. In terms of theoretical support, twoimportant contributions are a paper by Rack and Reimer in 1982 [10] and a definitivework by N. J. Higham in 2004 [5].Yet Higham’s analysis identifies an Achilles heel in (2.6). He shows that (2.1)has the gold standard property of backward stability: under standard assumptions offloating point arithmetic, the computed result p̂(x) delivered by (2.1) is the exactlycorrect value for a set of data {f̂j} that differ from {fj} by relative perturbations nogreater than about 5nu, where u is the unit roundoff (typically u ≈ 10−16). Formula(2.6), on the other hand, is not backward stable, satisfying only a more restrictiveforward stability bound. For x ∈ [−1, 1] and interpolation in Chebyshev points, whichis the most familiar case from applications such as ordinary Chebfun computations,there is little difference in the two bounds, and (2.6) is stable. As x moves away fromfrom [−1, 1], however, the forward stability bound grows rapidly, and (2.6) becomesunstable.Like many numerical instabilities, this one is a consequence of cancellation. Theaccuracy of (2.6) depends on the accuracy of (2.5), which we can rewrite asn∑j=0λjx − xj=1`(x).(3.1)Higham calls this “a mathematical identity that does not necessarily hold in floatingpoint arithmetic.” The problem is that as x moves farther from [−1, 1], 1/`(x) shrinksrapidly, so the validity of (3.1) relies increasingly on cancellation. In floating pointarithmetic, the result is rapid loss of accuracy, until soon there are no accurate digitsat all.Making the argument more quantitative explains why ellipses arise. The weightsλj on the left of (3.1) have size approximately 2n, and for x ∈ [−1, 1], the right-handside is of the same order since `(x) is approximately of size 2−n. Thus cancellationis not an issue. As x moves away from [−1, 1], however, `(x) grows to order approx-imately 2−nρn, where ρ is the parameter of the Bernstein ellipse passing through x.(That is, x = (w + w−1)/2 for some w with |w| = ρ [16].) This means that (3.1)relies on cancellation of magnitude ρn to occur, so we must expect loss of accuracyby approximately this factor. If ρn is of size u−1 or larger, all the accuracy will belost. If |x| is increased still further, so that ρn grows beyond u−1, then (3.1) will failto decrease any further and the computed p(x) from (2.6) will fail to increase anyfurther. As is apparent in the second image of Figure 1.1, we find ourselves with apolynomial p̂(x) that is bounded as |x| → ∞!Suppose a function f of size approximately 1 is interpolated in Chebyshev pointsto machine precision on [−1, 1], the basic operation of Chebfun. It is customary to4 M. WEBB, L. N. TREFETHEN AND P. GONNET  −2 −1 0 1 2−2−1012−30−20−100102030first barycentric formula second barycentric formulaFig. 3.1. Repetition of Figure 1.1 for f(x) = exp(x) sin(15x). Again the erroneous resultsappear approximately outside an ellipse, the Chebfun ellipse for this function.  −1 −0.5 0 0.5 1−1−0.8−0.6−0.4−0.200.20.40.60.81−16−14−12−10−8−6−4−2024  −1 −0.5 0 0.5 1−1−0.8−0.6−0.4−0.200.20.40.60.81−16−14−12−10−8−6−4−2024first barycentric formula second barycentric formulaFig. 3.2. Contour plot of numerically computed values log10|f(x) − p(x)| for formulas(2.1) and (2.6) with f(x) = exp(x) sin(15x). The first plot is unaffected by rounding errors,so it is a correct picture of interpolation errors only, showing accuracy precisely within theChebfun ellipse. The errors in the second image are significantly larger because of roundingerrors amplified by instability.define the Chebfun ellipse for f as that Bernstein ellipse whose parameter ρ satisfiesρn = u−1, where n is the degree of p, and it is inside this ellipse that p can beexpected to have some accuracy as an approximation to f . Because the Chebfunellipse is determined by the factor ρn, we can expect the numerical accuracy of thebarycentric formula to be lost entirely, approximately, if x lies outside it. In otherwords, for just the values of x for which p(x) has no accuracy as an approximation off , the evaluation of p(x) by (2.6) is also completely inaccurate because of numericalinstability.For many functions, the Chebfun ellipse is approximately the Bernstein ellipsethat passes through the closest singularity of f to [−1, 1], since that singularity con-trols the value of n needed to represent f to machine precision. However, the Cheb-fun ellipse may be smaller than this, for example if f is an entire function (analyticthroughout the complex plane). Figure 3.1 shows such a case, repeating Figure 1.1for the entire function f(x) = exp(x) sin(15x).Figure 3.2 compares the two barycentric formulas from another angle, looking atthe computed differences |f(x) − p(x)| between the function and its interpolant.BARYCENTRIC INTERPOLATION FORMULAS 5−4 −3 −2 −1 0 1 2 3 4−4−3−2−101234Re(z)Im(z)  12345quotient of polynomials rational barycentric formulaFig. 4.1. Contour plots of |r(z)| for type (9, 6) rational linearized least-squares fits in 42Chebyshev points in [−1, 1] to tanh(πx/2), based on calling Chebfun’s ratinterp commandwith (m, n) = (10, 10). The unstable computation on the right uses a rational barycentricformula, whereas on the left r(x) = p(x)/q(x) is computed by applying (2.1) separately top(x) and q(x).4. Practical implications. Confirming Higham’s theory of [5], we have shownthat the standard barycentric interpolation formula (2.6) should not be used for eval-uating a polynomial outside the interval of interpolation. One should use (2.1) in-stead, after reformulating it via logarithms to avoid over- and underflow. We close bymentioning two contexts in which this and a related observation may have practicalimportance.First, suppose one wants to evaluate a rational function r in the complex plane.This task arises frequently, because Padé approximations and their relatives offer thebest known techniques for many problems of extrapolation and analytic continuation.One of the attractions of the barycentric formula (2.6) is that it generalizes to rationalfunctions, as pointed out first by Schneider and Werner [1, 14], and this is the methodused by Chebfun up till now to evaluate the rational functions that result from in-terpolation and least-squares, minimax, Chebyshev–Padé, and Carathéodory–Fejérapproximation (Chebfun commands ratinterp, remez, chebpade, cf). We find thatthe instability described in this paper applies to the rational analogues of (2.6) too,and so it is not a good idea to evaluate r(x) in this way. Since it is not clear thata rational analogue of (2.1) exists, we recommend evaluating r(x) = p(x)/q(x) bytreating the numerator and denominator separately with (2.1) and then just takingthe quotient. This adjustment was crucial to the success of the explorations reportedin [18]. An example of stable and unstable evaluation of a rational function is shownin Figure 4.1.We will not give details for rational approximants in this paper, because this topicis considerably more involved than the polynomial case and we are currently inves-tigating how best to treat it. The rational interpolation scheme used for Figure 4.1,implemented in Chebfun’s ratinterp command, is an analogue for Chebyshev pointsin [−1, 1], as in [8], of the scheme presented in [3] for robust rational interpolationand least-squares fitting in roots of unity on the unit disk.Second, questions arise about evaluation of both polynomial rational functionseven within [−1, 1] when the interpolation points are far from the Chebyshev distri-bution. For polynomials, this was the example focused on by Higham to illustratethe possible instability of (2.6). As it happens, barycentric formulas have been used6 M. WEBB, L. N. TREFETHEN AND P. GONNETapparently to great advantage by Pachón for working with the irregular reference setsthat arise in the Remez algorithm for computing best polynomial and rational ap-proximations [7, 9]. It would appear that it may be worth revisiting these algorithmsin the light of our new practical experience of the instability of (2.6).Acknowledgments. This research originated with a EPSRC Vacation Bursaryawarded to MW to support a summer undergraduate research project.REFERENCES[1] J.-P. Berrut and L. N. Trefethen, Barycentric Lagrange interpolation, SIAM Rev. 46 (2004),501–517.[2] M. Dupuy, Le calcul numérique des fonctions par l’interpolation, C. R. Acad. Sci. 226 (1948),158–159.[3] P. Gonnet, R. Pachón and L. N. Trefethen, Robust rational interpolation and least-squares, Elect.Trans. Numer. Anal. 38 (2011), 146–167.[4] P. Henrici, Essentials of Numerical Analysis, Wiley, 1982.[5] N. J. Higham, The numerical stability of barycentric Lagrange interpolation, IMA J. Numer.Anal. 24 (2004), 547–556.[6] C. G. J. Jacobi, Disquisitiones Analyticae de Fractionibus Simplicibus, thesis, Berlin, 1825.[7] R. Pachón, Algorithms for Polynomial and Rational Approximation in the Complex Domain,DPhil thesis, Oxford University Mathematical Institute, 2010.[8] R. Rachón, P. Gonnet and J. van Deun, Fast and stable rational interpolation in roots of unityand Chebyshev points, SIAM J. Numer. Anal., submitted.[9] R. Pachón and L. N. Trefethen, Barycentric-Remez algorithms for best polynomial approximationin Chebfun, BIT Numer. Math. 49 (2009), 721–741.[10] H.-J. Rack and M. Reimer, The numerical stability of evaluation schemes for polynomials basedon the Lagrange interpolation form, BIT 22 (1982), 101–107.[11] M. Riesz, Eine trigonometrische Interpolationsformel und einige Ungleichungen für Polynome,Jahresber. Deutsch. Math.-Ver. 23 (1914), 354–368.[12] H. Rutishauser, Lectures on Numerical Mathematics, M. H. Gutknecht et al, eds., Birkhäuser,1990.[13] H. E. Salzer, Lagrangian interpolation at the Chebyshev points xn,ν = cos(νπ/n), ν = 0(1)n;some unnoted advantages, Computer J. 15 (1972), 156–159.[14] C. Schneider and W. Werner, Some new aspects of rational interpolation, Math. Comp. 47(1986), 285–299.[15] W. J. Taylor, Method of Lagrangian curvilinear interpolation, J. Res. Nat. Bur. Stand. 35(1945), 151–155.[16] L. N. Trefethen, Approximation Theory and Approximation Practice, book in preparation.[17] L. N. Trefethen and others, Chebfun Version 4.1, The Chebfun Development Team, 2011,http://www.maths.ox.ac.uk/chebfun/.[18] M. Webb, article in preparation to be submitted to SIAM Undergraduate Research Online.

Stability of Barycentric interpolation formulas for extrapolation

https://pure.manchester.ac.uk/ws/files/149436252/webbbarycentric.pdf

Stability of Barycentric interpolation formulas for extrapolation

Abstract

Similar works

Full text

Available Versions

The University of Manchester - Institutional Repository