Maruri-Aguilar, H

Trandafir, PC

English

Queen Mary Research Online

Sequential barycentric interpolationHugo Maruri-Aguilar and Paula Camelia TrandafirAbstract Polynomial interpolators may exhibit oscillating behavior which oftenmakes them inadequate for modeling functions. A well known correction to thisproblem is to use Chebyshev design points. However, in a sequential strategy it isnot very clear how to add points, while still improving polynomial interpolation.We present a sequential design alternative by allocating anextra observation wherethe difference between consecutive interpolators is largest. Our proposal is inde-pendent of the response and does not require distributionalassumptions. In simu-lated examples, we show the good interpolation performanceof our proposal and itsasymptotical convergence to the Chebyshev distribution.1 IntroductionIn classical optimal design theory, a connection can be establi hed between a certainoptimality criterion, a linear polynomial model and a design constructed with zerosof Chebyshev polynomials, i.e. Chebyshev points. This connection was first notedby Studden (1968). Subsequent research led to various articles and books which alsoexhibited designs whose points are zeros of Jacobi, Laguerre, or Hermite polyno-mials, among others. Among some of the better known examplesof the connectionbetween polynomial models and Chebyshev points, we may cite, classified by op-timality criteria, Pukelsheim and Torsney (1991) forA-optimal designs, Fedorov(1972, pp.85), Pázman (1986, pp.178), Pukelsheim (1993, pp.214–216) and Kar-Hugo Maruri-AguilarSchool of Mathematical Sciences, Queen Mary, University ofLondon, Mile End Road E1 4NSUK, e-mail: H.Maruri-Aguilar@qmul.ac.ukPaula Camelia TrandafirDepartamento de Estadı́stica e Investigación Operativa,Universidad de Valladolid, Prado de laMagdalena s/n, 47005, Valladolid, Spain, e-mail: camelia@eio.uva.es12 Hugo Maruri-Aguilar and Paula Camelia Trandafirlin and Studden (1966) forD-optimality, Dette (1993) forDs-optimality; and Dette(1993) and Heiligers (1998) forE-optimality.In the analysis of computer experiments, usually there is nora dom error asso-ciated with the response and models interpolate the observed response values. Avariety of models are available for the analysis of such experiments. Spline modelscan be used, but also models based on radial bases and krigingare have becomewidely used, see Müller (2001), O’Hagan (2006) and Fedorovand Müller (2007).We are concerned with sequential polynomial interpolation. Polynomials aresimple and potentially effective models, often with a straightforward interpretation.However, they have a tendency to oscillate between design points. Those oscilla-tions are called theRunge phenomenon, which in some cases can only get worseas the number of data points increases. It is a well-known classic l result that theoscillation caused by Runge’s phenomenon can be minimized by interpolating atthe Chebyshev nodes. In Epperson (1987), additional conditi s for mitigating thephenomenon are studied. In the literature of approximationtheory there are severalproposals which use Chebyshev points and may be used to interpolate. In Boyd andOng (2009) and Boyd and Xu (2009) the authors use a subsample of the uniform dis-tribution to generate “mock Chebyshev” (i.e. approximate Ch byshev) points, whilein Platte and Driscoll (2005), interpolation points are selected using more generalpolynomial approximation techniques. The methods producegood interpolation re-sults, although the strategies are non-sequential.The aim of the present paper is to introduce a sequential design strategy for uni-variate polynomial interpolation. Our proposal is based onan equivalent form ofthe Lagrange interpolator called barycentric interpolator. In Section 2 we reviewthe Lagrange and barycentric interpolators. In Section 3 wepresent our sequentialproblem, introduce a design algorithm and prove that it doesn t depend on responsevalues. Section 4 presents simulated examples to evaluate the performance of our al-gorithm. Conclusions, future work and a conjecture are presented in Section 5.2 Barycentric Lagrange interpolationLet Dn = {d1, . . . ,dn} be a design ofn distinct univariate points. The valuesf1, . . . , fn are observations, one for every design point. Those values are assumedto be evaluations of a deterministic (but unknown) functionwhich is to be interpo-lated. A classic solution is the Lagrange interpolatorgn(x) = ∑nj=1 f j ∏ni=1,i6= jx−did j−di.Lagrange interpolation exhibits numerical and computational drawbacks and an al-ternative form of it is available, known as the barycentric interpolator:gn(x) = mn(x)n∑j=1wn, j f jx−d j, (1)Sequential barycentric interpolation 3with barycentric weights defined bywn, j =(∏ni=1,i6= j(d j −di))−1and mn(x) =∏ni=1(x− di). The first subindex inwn, j denotes number of design points used forcomputing it, while the second subindex relates the weight to a design point.3 Sequential interpolationBarycentric formulæ allow sequential interpolation of data, that is, adding an extraobservation to an existing data set and updating the barycentri interpolator. Sequen-tial updating can be made part of an adaptive procedure, in which using informationfrom the interpolation process helps in selecting a new design point.3.1 Response based updateConsider two interpolators, one of which is considered to bemore accurate thanthe other, as it is built with one extra observation. We postulate that the differencebetween them can be used as a guide to future experimentation. In ther words, themore accurate interpolator may be used to validate the less accur te fit and to insertanother design point where this difference is largest over an arbitrary design region.We set the design region to[0,1] but it can be adapted to other design region[a,b].The starting point isgn(x), the barycentric interpolator as defined in Equation(1). Consider an extra design pointd∗n+1 and its corresponding observation. We termd∗n+1 a dummy pointand only require it to be different from existing design points.Denote byGn+1(x) the interpolator constructed with the temporary design consist-ing of the original design plus the dummy point,Dn∪d∗n+1. A new design point isselected according todn+1 = arg maxx∈[0,1]|Gn+1(x)−gn(x)|. (2)After the search, the dummy pointd∗n+1 is discarded and the original designDn isaugmented toDn+1 = Dn∪dn+1. The search problem is well posed, i.e. maximisa-tion of a bounded function over a closed compact set.Example 1.Consider the functionf (x) = 1/(1+25(2x−1)2), which is to be inter-polated using the designD10 = {0, 19, . . . ,1}. Let g10(x) be the interpolator functionconstructed with observations off (x) atD10. The dummy pointd∗11=12 is added tobuild an updated interpolatorG11(x). The next point is selected where the absolutedifference between the interpolatorsG11(x) andg10(x) is largest over[0,1]; this oc-curs at the points 0. 325 and 0.9675. Any of these two points can be added toD10andd∗11 is discarded.4 Hugo Maruri-Aguilar and Paula Camelia Trandafir3.2 Sequential design algorithmThe sequential design procedure described above simplifiesto a response-independentalternate maximization and update ofmn(x). We now describe the algorithm.Input An initial designDn of n distinct pointsd1, . . . ,dn ⊂ [0,1]; a numberk ofextra design points required.Output A set of additional runsdn+1, . . . ,dn+k ⊂ [0,1].Initialization Setmn(x) := ∏ni=1(x−di); set j := 0.Step 1 Maximize|mn+ j(x)|with respect tox, in the interval[0,1], i.e. letdn+ j+1 :=argmaxx∈[0,1] |mn+ j(x)|.Step 2 Updatemn+ j+1(x) := mn+ j(x)(x−dn+ j+1) and j := j +1. If j < k, repeatfrom Step 1.The algorithm does not depend on actual response values observed, but only onthe design points. Additionally, it does not depend on the actual location of thedummy point. These two characteristics are implied by Theorem 1, which is provenin the Appendix.Theorem 1. For n> 0, let gn(x) and Gn+1(x) for n> 0, be two barycentric interpo-lators, where gn(x) is defined as in Equation (1); and Gn+1(x) is constructed withan additional dummy design point d∗n+1. ThenGn+1(x)−gn(x) = mn(x)n+1∑j=1wn+1, j f j . (3)The right hand side of Equation (3) is the product ofmn(x), which depends onx and onDn (but not on the dummyd∗n+1), and a second quantity∑n+1j=1 wn+1, j f jthat depends on design points and responses (including dummy data), but not onx and thus it can be ignored when searching for the new design point. Theorem1 makes the search for a new design point independent of the response, indeed itmakes Equation (2) equivalent todn+1 = argmaxx∈[0,1] |mn(x)|.Example 2.Consider again the design of Example 1. The sequential algorithm isapplied fork= 10 extra runs and pointsd11 to d20 are sequentially obtained: 0.0325,0.9684, 0.9335, 0.0662, 0.8306, 0.1677, 0.4999, 0.0099, 0.9902 and 0.2813.A special condition arises from Equation (3), when∑n+1j=1 wn+1, j f j = 0 holds. Thisimplies thatGn+1(x)−gn(x) ≡ 0, in other words, that the dummy pointd∗n+1 doesnot update the interpolator and consequently, this step does n t yield information forthe next design pointdn+1. This condition appears, for instance, when allf j valuesare equal. This could occur when sampling a constant functioor a periodic functionat the same point in every period. A different instance appears when response datatruly comes from a polynomial of degree at mostn−1. In any of the above situa-tions, any point in the interval[0,1] could be selected as new design point. However,in all the examples we tried, none of them occurred and we suggest that they shouldnot be cause of concern.Sequential barycentric interpolation 54 Performance and large sample propertiesIn this section we first evaluate the performance of our sequential design strategy tointerpolate. We then study the large sample properties of our sequential designs.4.1 Interpolating performanceThe accuracy of polynomial interpolators with our sequential design algorithmwas assessed in a simulation study. The following four functio s with domain[0,1] were used as test functions:1(x) is the function of Example 1;s2(x) =120 exp(u1/3)sin(u/2)φ(u) with φ(u) the Heaviside unit step function andu =60x−30;s3(x) = sin(10x) ands4(x) = 21−cos(v)v2with v= 35x−15. The functionswere selected to exhibit features which are not easy for modeling with polynomials,such as flat regions followed by regions with sharp change, orperiodic behavior.Fig. 1 Maximum distancebetween simulators andbarycentric interpolators,plotted against number ofextra pointsk added.1100.110−210−310−410−5s2s1s4s3n+k11 15 20 25From a uniform design of size ten, fifteen points were sequentially added, to-talling 25 points. At every step, barycentric interpolators were fitted independentlyfor each function. The maximum distance between the true functio s1(x), . . . ,s4(x)and its barycentric interpolator, over the design region, was recorded. Decreasingvalues of this distance show good approximation, while or inc easing values pointto the presence of Runge phenomenon.The results are plotted in Figure 1, where a decreasing trendis evident, thusshowing good approximation to simulators for all cases. Convergence to the truefunctions3 was faster than the other cases, while convergence was slower f r s2.4.2 Large sample properties of pointsThe points generated with our algorithm cluster in the borders of the design region.We studied whether the points converge asymptotically to a known distribution.6 Hugo Maruri-Aguilar and Paula Camelia TrandafirTo study large sample behavior, points were sequentially added to each of the fol-lowing eight initial designs of sizen: uniform designs 0, 1n−1, . . . ,1 (termed UI) and1n+1, . . . ,nn+1 (termed UII); firstn points of Sobol’s space filling sequence (termedS), see Bratley and Fox (1988); Chebyshev type I and II points(CI and CII, re-spectively), see Berrut and Trefethen (2004); the designs labelled TI and TII weregenerated by transforming UI and UII to the symmetric triangular distribution withmode in 12; and a design with random points (termed R). We used initial designsizesn = 5,10,35,50,100,150, in each case sequentially adding points with ouralgorithm up to one thousand. Two statistics were computed:a) Quantile-Quantile(QQ) plot and b) goodness of fit Kolmogorov-Smirnoff (KS) statistic. The Beta dis-tributionβ (12,12) (also known as the Chebyshev or arcsine distribution) was used incomputations. This choice was suggested by literature on polynomial interpolationconvergence Berrut and Trefethen (2004).Fig. 2 QQ plot for design UI.We show results forn= 50, which are representative of the results for other initialsizes. Figure 2 shows QQ plot for UI, which converges to Chebyshev distributionwith 20 extra points. QQ plots for other initial designs exhibit similar pattern.0.010.10.30.550 55 60 65 70 75 80TIIRSCIIBBB0.0010.010.10.30.550 80 100 200 500 800 1000TIIRSCIIFig. 3 KS statistic vs. design size.The left-hand panel is a close up view of the right-hand panel.The KS statistic “grouped” designs according to their initial performance, withbest (i.e. low) values obtained by Chebyshev points (CI, CII). In second place wereuniform designs (UI, UII, S), followed by random design R. The worst values wereSequential barycentric interpolation 7observed for designs with points clustered in the centre of the design region (TI, TII).Figure 3 shows the evolution of KS statistic for one design for each of the observed“groups”: CII, R, S and TII. After adding about two times as many points as theinitial design size, designs behave similarly, showing non-monotonous decreasinglinear trend (in the log-log scale) for the KS statistic. Simulation results suggest avalue for the slope of the linear trend between−0.8 and−1, see also Figure 3.5 Discussion and future workWe introduced a univariate sequential adaptive design algorithm. In the exampleswe tried, the algorithm produced good points for polynomialinterpolation, whichconverged rapidly to the Chebyshev distribution and lead tothe following claim.Conjecture 1.For any initial design in[0,1], as the number of extra pointsk tends toinfinity, the algorithm of Section 3.2 produces samples fromthe distributionβ (12,12);and the KS statistic is of orderO(kα ), with α a suitable constant.Acknowledgements The first author was supported by Research Council UK grant EP/D048893/1(MUCM), while the second author by Spain’sMinisterio de Ciencias e Innovación(ECO2008-02358/ECON) andConserjerı́a de Educación de la Junta de Castilla y León(VA004B08).Appendix A: Proof of Theorem 1The barycentric interpolatorGn+1(x) isGn+1(x) = mn(x)(x−d∗n+1)(n∑j=1wn+1, j f jx−d j+wn+1,n+1 fn+1x−d∗n+1),where barycentric weightswn+1, j are computed using design and dummy point. Wehave thatGn+1(x)−gn(x) = mn(x)A, whereA= (x−d∗n+1)(n∑j=1wn+1, j f jx−d j+wn+1,n+1 fn+1x−d∗n+1)−n∑j=1wn, j f jx−d j.We now show thatA does not depend onx. After simplifying and using the updatingformula of barycentric weightswn+1, j = wn, j(d j −d∗n+1)−1, we haveA =n∑j=1f jx−d j((x−d∗n+1)wn, j(d j −d∗n+1)−wn, j)+wn+1,n+1 fn+1=n∑j=1wn, j f j(d j −d∗n+1)+wn+1,n+1 fn+1 =n+1∑j=1wn+1, j f j .8 Hugo Maruri-Aguilar and Paula Camelia TrandafirReferencesJ.-P. Berrut and L. N. Trefethen. Barycentric Lagrange interpolation. SIAM Rev.,46(3):501–517 (electronic), 2004.J. P. Boyd and J. R. Ong. Exponentially-convergent strategies for defeating theRunge phenomenon for the approximation of non-periodic functio s. I. Single-interval schemes.Commun. Comput. Phys., 5(2-4):484–497, 2009.J. P. Boyd and F. Xu. Divergence (Runge phenomenon) for least-squares polynomialapproximation on an equispaced grid and Mock-Chebyshev subset interpolation.Appl. Math. Comput., 210(1):158–168, 2009.P. Bratley and B. L. Fox. ALGORITHM 659 Implementing Sobol’squasirandomsequence generator.ACM Trans. Math. Soft., 14(1):88–100, 1988.H. Dette. A note onE-optimal designs for weighted polynomial regression.Ann.Statist., 21(2):767–771, 1993.H. Dette. On a mixture of theD- andD1-optimality criterion in polynomial regres-sion. J. Statist. Plann. Inference, 35(2):233–249, 1993.J. F. Epperson. On the Runge example.Amer. Math. Monthly, 94(4):329–341, 1987.V. V. Fedorov. Theory of optimal experiments. Academic Press, New York, 1972.Translated from the Russian and edited by W. J. Studden and E.M. Klimko, Proba-bility and Mathematical Statistics, No. 12.V. V. Fedorov and W. G. Müller. Optimum design for correlated fields via covariancekernel expansions. InmODa 8—Advances in model-oriented design and analysis,Contrib. Statist., pages 57–66. Physica-Verlag/Springer, H idelberg, 2007.B. Heiligers.E-optimal designs for polynomial spline regression.J. Statist. Plann.Inference, 75(1):159–172, 1998.S. Karlin and W. J. Studden. Optimal experimental designs.Ann. Math. Statist.,37:783–815, 1966.W. G. Müller. Collecting spatial data. Contributions to Statistics. Physica-Verlag,Heidelberg, revised edition, 2001. Optimum design of experim nts for randomfields.A. O’Hagan. Bayesian analysis of computer code outputs: a tutorial. ReliabilityEngineering & System Safety91 1290–1300, 2006.A. Pázman. Foundations of optimum experimental design, volume 14 ofMath-ematics and its Applications (East European Series). D. Reidel Publishing Co.,Dordrecht, 1986. Translated from the Czech.R. B. Platte and T. A. Driscoll. Polynomials and potential theory for Gaussian radialbasis function interpolation.SIAM J. Numer. Anal., 43(2):750–766 (electronic),2005.F. Pukelsheim.Optimal design of experiments. John Wiley & Sons Inc., New York,1993.F. Pukelsheim and B. Torsney. Optimal weights for experimental designs on linearlyindependent support points.Ann. Statist., 19(3):1614–1625, 1991.W. J. Studden. Optimal designs on Tchebycheff points.Ann. Math. Statist, 39:1435–1447, 1968.

Sequential Barycentric Interpolation

https://qmro.qmul.ac.uk/xmlui/bitstream/123456789/85887/2/Aguilar%20Sequential%20Barycentric%20Interpolation%202010%20Accepted.pdf

Sequential Barycentric Interpolation

Abstract

Similar works

Full text

Available Versions

Queen Mary Research Online