Multivariate Lagrange Interpolation at Sinc Points Error Estimation and Lebesgue Constant

This paper gives an explicit construction of multivariate Lagrange interpolation at Sinc points. A nested operator formula for Lagrange interpolation over an m-dimensional region is introduced. For the nested Lagrange interpolation, a proof of the upper bound of the error is given showing that the error has an exponentially decaying behavior. For the uniform convergence the growth of the associated norms of the interpolation operator, i.e., the Lebesgue constant has to be taken into consideration. It turns out that this growth is of logarithmic nature O((logn)m). We compare the obtained Lebesgue constant bound with other well known bounds for Lebesgue constants using different set of points.


Introduction
One of the fundamental problems in approximation theory is the Lagrange interpolation problem: given a continuous function f , a set of interpolation points X and a polynomial interpolation space Π, and P n an element of Π which is equal to f at the interpolation points X. Although the statement of this problem is direct, many open questions remain on this topic, especially in the case of multivariate interpolation. Two factors affected the quality of approximation for high degree polynomials: the smoothness of the function to be interpolated, and the locations of the interpolation points that completely determine the interpolation operator. The norm of the interpolation operator gives a bound for the interpolation error by Lebesgue constant. Norms of the interpolation operator are usually difficult to compute and their values are very sensitive to the location of the interpolation points.
There is a well developed theory that quantifies the convergence or divergence of polynomial interpolants. A key notion is that of the Lebesgue constant, Λ n , for interpolation in a given set of points. The Lebesgue constant is the 1−norm of the linear mapping from data to the interpolant: where ∥.∥ denotes the ∞-norm on [-1, 1]. In words, if the data values on an (n + 1)-point grid is given, and the data are not greater than 1 in absolute value, what is the largest possible value of the interpolant P somewhere in [-1, 1]?
The interpretation of the Lebesgue constant as a convergence factor can be stated as: Given a set of points X and let P n be the interpolant in these points to a function f , we know that ∥P n ∥ might be as great as Λ n . Thus ∥ f − P n ∥ might be as great as Λ n − 1, showing that a large Lebesgue constant rigorously implies the possibility of a large interpolation error. For example, interpolation using equally spaced points has been shown to yield undesirable divergent behaviors even for analytical functions as the polynomial degree of interpolation increases. This behavior is called Runge phenomenon, see (Runge, 1901). It means that the Lebesgue constants must be large for equispaced interpolation.
This means that a small Lebesgue constant implies that interpolation will be close to the best.
For the 1D Lagrange interpolation problem, the explicit construction of the interpolating polynomial for certain interpolation points and the derivation of remainder formulas, we refer the reader to the works (Smith, 2006), (Youssef et al., 2016), and (Rivilin, 1974) and references therein. This interpretation of Lebesgue constant is not restricted to polynomial approximation but can be used to judge the accuracy and stability of any approximation. For example, a discussion for Lebesgue constant of Sinc approximation can be found in (Stenger et al. 2014).
In case of multivariate interpolation, much less is known about the optimal point sets and the growth of the optimal Lebesgue constant. Theoretical results on the Lebesgue constant are known in only a few special cases, and explicit constructions for nearly optimal point sets for total degree interpolation are discussed only in the case of bivariate interpolation. But few results have been introduced for higher dimensions. One of the earliest result has been obtained by Subbotin in (Subbotin, 1982). He showed that in the m-dimensional case for Lagrange approximation we have: where C 1 and C 2 are constants and independent of n.
Recently some other works discuss the Lebesgue constants for bivariate Lagrange approximation using different sets of points. For example, using Padua and Xu points as interpolation points can yield a rate of O(log(n)) 2 , see , (Bos and Xu et al., 2006), and (Vecchia et al., 2009). Also using Lissajous points can follow the same asymptotic behavior, see (Erb et al., 2015). One of the recent papers that discuss the bivariate and multivariate case is (Gunzburger and Teckentrup, 2014). The paper discusses numerical results for multivariate cases, for moderate dimensions, using Gauss Chebyshev Lobatto points.
In the current paper we discuss the multivariate Lagrange approximation. We give an explicit form of the approximation, prove the upper bound of the error and, the upper bound of the Lebesgue constant.
The paper is organized as follows: Section 2 introduces a short recall on the univariate Lagrange approximation at Sinc points discussed in (Stenger et al., 2013) and (Youssef et al., 2016). Section 3 gives the explicit definition of the multivariate Lagrange approximation and the main theorems of error bound and Lebesgue constant. Section 4 includes the numerical discussion and finally the conclusion in section 5.

Univariate Approximation
In this section we introduce the Lagrange approximation at Sinc points as interpolation points. To define these interpolation points let Z denote the set of all integers. Let R be the real line, and C denote the complex plane. Let h denote a positive parameter and let k ∈ Z, x ∈ C. Let d denote a positive number and let D ⊂ C be a simply connected region defined as: be a conformal map of D onto the strip Let Γ = [a, b] = ϕ −1 (R) be an arc, where a = ϕ −1 (−∞) and b = ϕ −1 (∞) denote the end points of Γ. Then we define the set of Sinc points by Finally, let α ∈ (0, 1] and β ∈(0,1] denote fixed positive numbers and set ρ = e ϕ(x) . Without loss of generality, let us restrict d introduced above to the interval (0, π/2). Let L α,β (D) denote the family of all functions f that are analytic in D, such that for all z ∈ D, we have The space of functions M α,β (D) denotes the set of all functions q defined on D that have finite limits q(a) = lim z→a q(z) and where the limits are taken from within D, and such that f ∈ L α,β (D), where, Now we define a family of polynomial-like approximation that interpolate given Sinc data of the form where the x k are Sinc points. This novel family of Lagrange polynomials was recently derived in (Stenger et al., 2013).
, there exists a unique polynomial P n (x) of degree at most n − 1 satisfying the interpolation condition, In this case P n (x) can be expressed as: with, where, This approximation, like regular Sinc approximation, yields an exceptional accuracy in approximating the function that is known at Sinc points, see (Stenger, 2010). Unlike Sinc approximation, it gives an exponential convergence rate when differentiating the interpolation formula given in (6), (Stenger et al., 2013).
Theorem 1. Error Bound: Let h = π √ N , and let {x k } N k=−N denote the Sinc points as defined in (5). Let f be in M α,β (D), and let P n (x) be defined as in (6). Then there exist two constants A > 0 and B > 0, independent of N, such that Proof. For the proof of (8), see (Stenger et al., 2013).
To define the Lebesgue constant and Lebesgue function, let P n be defined as in (6) then the Lebesgue constant is the infinity-norm of the linear mapping from data to an interpolant and can be written down as: Let λ n (x) be the associated Lebesgue function defined by, then the Lebesgue constant is defined as, In (Youssef et al., 2016) an estimate for the upper bound of the Lebesgue constant for Lagrange approximation at Sinc points (5) has been derived as Λ n ≈ 1 π log(n + 1) + 1.07618.
Next we extend these results of error and Lebesgue constant for the one dimensional case to the multidimensional case. For the rest of the paper, we introduce the notation Λ n,m to denote the Lebesgue constant using n interpolation points in m dimensional spaces.

Multivariate Approximation
In this section we define a multivariate form of Lagrange approximation introduced in the previous section. (3) and (4) with a = −1, b = 1 and, d = π/2. Define D to be the closure of D.
Define Ω to be Define the following space of functions S α,d to be the space of all functions f such that, where, Hol(D) is the family of all functions f that are analytic in a domain D. And a function f is said to be in a class for all points x 1 and x 2 on the interval [a, b].
Define the following one dimensional interpolation operator, where b j are the basis defined as in (7).
Then Lagrange approximation of a function f (X) can be defined by a nested operator as, In (Stenger, 2010), the notation of one dimensional interpolation operators has been used to define 2D Sinc approximation and derive the upper bound of the error for such approximation. Here we intend to do this but for m-dimensional case.
The following two corollaries follow from the univariate case discussed in the previous section.
Corollary 2. If the operator L i is defined as in (11) then, where C i is a constant independent of N.
Proof. The statement follows by inspection of (9).
Corollary 3. If the operator L i is defined as in (11) then, where A i and B i are constants independent of N.
Proof. The claim is a direct consequence of (8) and (11).
Theorem 4. Upper Bound of Error Let h = π √ N and let f (X) be in S α,d (D), and let (P n f ) be defined as in (12). Then there exist two sets of constants C i > 0, γ i > 0, i = 1, ..., m, are constants independent of N, such that Proof. The difference of f (X) and the interpolant P n ( f ) defined in (12) can be written as: .

Now using Corollaries 2 and 3 finally delivers,
From Theorem 4 the next corollary directly follows.
Corollary 5. The multivariate error has the following behavior, Theorem 6. Lebesgue constant: If P n (X) is defined as in (12), then:

Numerical Results
The following numerical experiments examine the calculations of the error and Lebesgue constant for multivariate Lagrange approximation at Sinc points. For the error we give a function example to verify the exponential decaying error in (15). For the Lebesgue constant, we will use data published in literature to compare our computations with these results in (Brutman, 1997), (Vecchia et al., 2009), (Erb et al., 2015) and, (Gunzburger and Teckentrup, 2014). The comparison is based on the availability of data. Most previous studies are dedicated to the bivariate case. However some few sets of data are related to higher dimensions, see (Gunzburger and Teckentrup, 2014). To this end, we first introduce the bivariate case. In this case, we compare with bivariate Lagrange approximations using different sets of points, like Chebyshev, Padua, Xu, Lissajous points discussed in (Rivilin, 1974), (Brutman, 1997), (Bos and Xu et al. 2006), (Caliari et al., 2005) and, (Erb et al., 2015). We then compare the Lagrange approximations at Sinc data with Sinc approximation. The last set of experiments are related to higher dimensions where we compare available Lagrange data with our Lagrange-Sinc interpolation. In total we performed 6 experiments. In all of the calculations we used the barycentric formula for Lagrange approximation at Sinc points introduced in (Youssef et al., 2016).
• Experiment 1: The first experiment illustrates the exponential convergence rate of the Poly-Sinc algorithm (Theorem 4), we examine the function f (x, y) = 1/(1 + r 2 ) defined on Q = [−1, 1] 2 , where r 2 = x 2 + y 2 . For this function we find the Poly-Sinc approximation defined in (12) with m = 2 for different numbers of Sinc points n = 2N + 1. For each n, we compute the norm error. As result a table for each n and the corresponding error has been created. We then use this table in a least square estimation to find the coefficients of the error function as estimated in (15). Specifically, we used the form where ξ, τ, η and, µ are constants. A least square fit to the collected data delivers the constants ξ = 5.3, τ = 40.6, η = 0.4 and, µ = 3.8 that optimally represent the error list. In Fig. 1 the solid line represents (15) using the constants ξ, τ, η and, µ from the least square fit. Dots in Fig.  1 are the discrete norm errors. The graph in Fig. 1 demonstrates that the error of the Poly-Sinc approximation for f (x, y) = 1/(1 + r 2 ) follows an exponentially decay relation.
• Experiment 2: The second experiment is concerned with bivariate Lagrange approximation of a function defined on a square [−1, 1] 2 . For the x and y directions we use the Sinc points (e kh − 1)/(1 + e kh ), (Stenger, 2010). A sample of these Sinc points on the square is given in Fig. 2, top. The evaluation of the Lebesgue function, λ n,2 (x, y), using N = 12 of Sinc points as interpolation points is shown in Fig. 2, bottom left. This graph shows that the Lebesgue function λ n,2 (x, y) oscillates with a certain amplitude over the square of study. The amplitude of the oscillations return to the same peak level and has a minimal value of 1. These symmetric repeated behavior is shown in Fig.  2, bottom right. Finally, the Lebesgue constant Λ n,2 is obtained as the maximum value of the Lebesgue function λ n,2 (x, y) by varying the number of interpolation points. The results for the Lebesgue constants Λ n,2 are shown in Fig. 3. Fig. 3 also includes the predicted behavior from (16) for m = 2, dashed line. The asymptotic growth of Lebesgue constant corresponds to the order ( log(n) ) 2 .
• Experiment 3: The third experiment compares the Lebesgue constant Λ n,2 for 2D Lagrange approximation at Sinc points with the least squares fitting of the Lebesgue constant for Chebyshev points discussed in (Rivilin, 1974) and (Brutman, 1997). Applying the results of theorem 6 to the two limits derived by Rivilin (Rivilin, 1974) and Brutman (Brutman, 1997) for polynomial interpolation as Λ n,2 = ( 2 π log(n + 1) + 1 ) 2 and Λ n,2 = ( 2 π ( log(n + 1) + log where γ = 0.557.. is the Euler-Mascheroni constant. We observe in Fig. 4 (left panel) that a Lagrange-Sinc interpolation is located below the two reference limits of Rivilin (solid) and Brutman (dots). In the right panel of Fig. 4 we also compare with the least squares fitting of the Lebesgue constant for the Padua and the Xu points. As proposed in (Erb et al., 2015), for the Padua points the relation Λ 2n = (2/π log(2n + 1) + 1.1) 2 , in (Caliari et al., 2005) and for the Xu points Λ 2n+1 = (2/π log(2n + 2)) 2 , in (Bos and Xu et al., 2006) should hold. The comparison with these estimated best fits for Padua and Xu points clearly indicates that the Padua points (solid) and the Xu points (dots) are both above the Lagrange-Sinc interpolation.
• Experiment 4: In this experiment we compare the Lebesgue constant Λ n,2 for bivariate Lagrange approximation at Sinc points with the least squares fitting of the Lebesgue constant for Fekete and Quasi-Lebesgue points introduced in (Briani et al., 2012). In (Briani et al., 2012) the data have been given in a  (Gunzburger and Teckentrup, 2014). Also we compare our results with the discrete values for Gauss Lobatto points introduced in (Gunzburger & Teckentrup, 2014). Rhombohedrals again indicate the calculated Lebesgue constants for Sinc points, Fig. 5 (left panel). The dashed line of the left panel represents (16). The right panel in Fig. 5 compares a fit for Gauss Lobatto points using Λ n,2 = 0.58 ( log 2 (n + 1) + 1.64 ) (solid line), (Gunzburger and Teckentrup, 2014), with the Sinc points (dashed line and rhombohedrals).
• Experiment 5: This experiment related to bivariate Lagrange-Sinc approximation and Sinc approximations in (Stenger, 2010).  (Rivilin, 1974) and Brutman's (Brutman, 1997) bounds generalized to the bivariate case. We note without showing the result that for a double exponential Sinc approximation the same asymptotic behavior is present.
• Experiment 6: As we mentioned in the introduction much less is known about the behavior of Lebesgue constants for higher dimensions. In this experiment we introduce the calculations of Lebesgue constant for Lagrange approximation for higher dimensional cases. In each case we compare between the numerically calculated results and the analytic result given by (16). We calculated Lebesgue constant for the case m > 2, namely m = 3, 4, 5, and 6. For the comparison, we used the best fitting derived from Rivilin and Brutman formulas, derived by Rivilin (Rivilin, 1974) and Brutman (Brutman, 1997). We also include a comparison with some numerical results for multivariate Lagrange approximation using Gauss Chebyshev Lobatto points, in (Gunzburger & Teckentrup, 2014). These values of Lebesgue constants in (Gunzburger and Teckentrup, 2014) are given in a discrete form for a small number of interpolation points. The results are given by Fig. 7. In this figure, the dashed line is for the analytic formula given in (16) (Rivilin, 1974), solid line, and Brutman (Brutman, 1997), dots, generalized to the bivariate case (left panel). Lebesgue constants for Sinc points are indicated by rhombohedrals following relation (16). The right panel compares the Lebesgue constants based on Sinc points with Padua points in (Caliari et al., 2005) (solid line) and Xu points in (Bos and Xu et al., 2006) (dots).   (Stenger, 2010). The right panel compares calculated Lebesgue constant for Lagrange (rhombohedrals) and the Lebesgue constant for Sinc given by stars ( * ) and the analytic upper limit in (16). The solid line represents Rivilin (Rivilin, 1974) and the dots for Brutman (Brutman, 1997) bounds generalized to the bivariate case.  (16) for different dimensions (m = 3, 4, 5, and 6). For a comparison we used the Rivilin (Rivilin, 1974) and Brutman (Brutman, 1997) generalization to dimension m, solid lines. The stars represents the numeric Lebesgue constant using Gauss-Lobatto points calculated in (Gunzburger and Teckentrup, 2014).

Conclusion
This paper discussed the multivariate Lagrange approximation at Sinc points. An upper bound of the error and a generalized Lebesgue constant are introduced. A numerical verification of the exponential decaying behavior of the error is tested by a function example. Finally, we compared the convergence/stability factor, Lebesgue constant/function with different types of points. We can conclude that Lagrange approximation at Sinc points delivers approximation results closer to the conjectured optimal approximation than Chebyshev approximations using barycentric formulas.