Numerical differentiation on scattered data through multivariate polynomial interpolation

We discuss a pointwise numerical differentiation formula on multivariate scattered data, based on the coefficients of local polynomial interpolation at Discrete Leja Points, written in Taylor's formula monomial basis. Error bounds for the approximation of partial derivatives of any order compatible with the function regularity are provided, as well as sensitivity estimates to functional perturbations, in terms of the inverse Vandermonde coefficients that are active in the differentiation process. Several numerical tests are presented showing the accuracy of the approximation.

Let us consider a set (1.1) σ = {x 1 , . . . , x m } of m = d+s s pairwise distinct points in R s and let us assume that they are unisolvent for Lagrange interpolation in Π d (R s ), that is for any choice of y 1 , . . . , y m ∈ R, there exists and it is unique p ∈ Π d (R s ) satisfying x any point in R s and by fixing the basis the Vandermonde matrix centered at x and by setting, using matrix notation, where c = [c α ] T |α|≤d ∈ R m is the solution of the system (1.6). This approach, for s = 2 and x the barycenter of the node set σ, has been recently proposed in [11] in connection with the use of the P A = LU factorization of the matrix V x (σ).
The main goal of the paper is to provide a pointwise numerical differentiation method of a target function f sampled at scattered points, by locally using the interpolation formula (1.7).
The key tools are the connection to Taylor's formula via the shifted monomial basis (1.3), suitably preconditioned by local scaling to reduce the conditioning of the Vandermonde matrix, together with the extraction of Leja-like local interpolation subsets from the scattered sampling set via basic numerical linear algebra. Our approach is complementary to other existing techniques, based on least-square approximation or on different function spaces, see for example [6,7,1,14] with the references therein. In Section 2 we provide error bounds in approximating function and derivative values at a given point x, as well as sensitivity estimates to perturbations of the function values, and in Section 3 we conduct some numerical experiments to show the accuracy of the proposed method.

Error bounds and sensitivity estimates
In the following we assume that Ω ⊂ R s is a convex body containing σ and that the sampled function f : Ω → R is of class C d, 1 (Ω), that is f ∈ C d (Ω) and all its partial derivatives of order d are Lipschitz continuous in Ω. Let K ⊆ Ω compact convex: we equip the space C d,1 (K) with the semi-norm [12] (2.1) and by R T [f, x] (x) the corresponding remainder term in integral form [17] with the multi-indices β following the order specified in Section 1. Let us denote by ℓ i (x) the i th bivariate fundamental Lagrange polynomial. Since in order to control the conditioning, it is useful to consider the scaled canonical polynomial basis centered at x in such a way that (x − x) /h belongs to the unit disk (cf. [11]), and to rewrite the expression (2.5) in the scaled basis (2.7) as where a i α,h = a i α h |α| . The interpolation polynomial p [y, σ] (1.7) can also be expressed in the basis (2.7) as In the following we denote by V x,h (σ) the Vandermonde matrix in the scaled basis (2.7) and by B h (x) the ball of radius h centered at x. Proposition 2.2. Let x ∈ Ω and f ∈ C d,1 (Ω). Then for any x ∈ K = B h (x) ∩ Ω and for any ν ∈ N s 0 such that |ν| ≤ d, we have Proof. Since By representing f (x) and f (x i ) = y i in truncated Taylor series of order d centered at x (2.2) with integral remainder (2.3), we obtain On the other hand since interpolation of degree d at the nodes in σ reproduces exactly polynomials of total degree less than or equal to d. Therefore by (2.4) we have Consequently, by substituting (2.14) in (2.13), we get By applying the differentiation operator D ν to the expression (2.15) and by using the triangular inequality, we obtain where [12, Lemma 2.1] The estimates in Theorem 2.2 can be written in terms of the Lebesgue constant of interpolation at the node set σ, defined by 1 (Ω) and B h (x) ⊂ Ω. Then for any x ∈ K = B h (x) and for any ν ∈ N s 0 such that |ν| ≤ d, we have where k j = s j (j−1)! for j > 0, k 0 = 1, and In particular, for x = x, the following inequality holds is a polynomial of total degree less than or equal to d, by repeatedly applying the Markov inequality [18] for a ball with radius h in the form and recalling that each partial derivative lowers the degree by one, we easily obtain Based on the inequality (2.11) in Theorem 2.2, we have and since x ∈ B h (x), it follows that The inequalities between multi-indices are interpreted componentwise, that is α ≤ β if and only if α i ≤ β i , i = 1, . . . , s. In particular, for x = x, the following inequality holds Proof. By using the expression of the fundamental Lagrange polynomial ℓ i (x) in the scaled basis (2.8) and by applying the differentiation operator D ν to the expression (2.15), we obtain (2.25) By taking the modulus of both sides of (2.25) and by using the triangular inequality, we get Therefore, (2.16), (2.26) and (2.20) imply Evaluating (2.27) at x yields to (2.24).
In line with [11] we can write the bounds in Proposition 2.4 in terms of the 1-norm condition number, say cond h (σ), of the Vandermonde matrix V x,h (σ).
Corollary 2.5. Let x ∈ Ω and f ∈ C d,1 (Ω). Then for any x ∈ K = B h (x) ∩ Ω and for any ν ∈ N s 0 such that |ν| ≤ d, we have The inequalities between multi-indices are interpreted componentwise, that is α ≤ β if and only if α i ≤ β i , i = 1, . . . , s. In particular, for x = x, the following inequality holds Proof. By applying the triangular inequality to (2.27), we get .
The results of Table 1 in Section 3 show that the bounds (2.19) and (2.29) are much larger than (2.24), which is only based on the "active coefficients" in the differentiation process. Therefore, in the analysis of the sensitivity to the perturbation of the function values, we use only the "active coefficients".
where k j = s j (j−1)! for j > 0, k 0 = 1. The inequalities between multi-indices are interpreted componentwise, that is α ≤ β if and only if α i ≤ β i , i = 1, . . . , s. In particular, for x = x, the following inequality holds Proof. Denoting by y = y + ∆y, where ∆y = [∆y i ] i=1,...,m corresponds to the perturbation on the function values y = [y i ] i=1,...,m , by (2.12) we get Then for any ∆y such that |∆y| ≤ ε, Since, using (2.26), we have Remark 2.7. It is worth observing that the quantity is the "stability constant" of pointwise differentiation via local polynomial interpolation, namely the value at the center of the "stability function" for the ball, that is Notice also that in view of (2.24) and (2.33) the overall numerical differentiation error, in the presence of perturbations on the the function values of size not exceeding ε, can be estimated as For the purpose of illustration, in Table 2 and Figure 14 of Section 3 we show the magnitude of the stability constant (2.33) relative to some numerical tests.
Remark 2.8. The previous results are useful to estimate the error of approximation in several processes of scattered data interpolation which use polynomials as local interpolants, like for example triangular Shepard [10,3], hexagonal Shepard [9] and tetrahedral Shepard methods [4]. They are also crucial to realize extensions of those methods to higher dimensions [8].

Numerical experiments
In this section we provide some numerical tests to support the above theoretical results in approximating function, gradient and second order derivative values. We fix s = 2, Ω = [0, 1] 2 and we take different positions for the point x in Ω: at the center, and near/on a side and a corner of the square. We use different distributions of scattered points in Ω, namely Halton points [13] and uniform random points. We focus on the scattered points in the ball B r (x) centered at x for different radii, from which we extract an interpolation subset σ of m = d+2 2 Discrete Leja Points computed through the algorithm proposed in [2] (see Figure 1).
The reason for adopting Discrete Leja Points is twofold. We recall that they are extracted from a finite set of points (in this case the scattered points in the ball) by LU factorization with row pivoting of the corresponding rectangular Vandermonde matrix. Indeed, Gaussian elimination with row pivoting performs a sort of greedy optimization of the Vandermonde determinant, by searching iteratively the new row (that is selecting the new interpolation point) in such a way that the modulus of the augmented determinant is maximized. In addition, if the polynomial basis is lexicographically ordered, the Discrete Leja Points form a sequence, that is the first k+s s are the Discrete Leja Points for interpolation of degree k in s variables, 1 ≤ k ≤ d; see [2] for a comprehensive discussion.
Then, on one hand Discrete Leja Points provide, with a low computational cost, a unisolvent interpolation set, since a nonzero Vandermonde determinant is automatically seeked. On the other hand, since they are computed by a greedy maximization, one can expect, as a qualitative guideline, that the elements of the corresponding inverse Vandermonde matrix (that are cofactors divided by the Vandermonde determinant), and thus also the relevant sum in the error bound (2.24) as well as the condition number, are not allowed to increase rapidly. These results are in line with those shown in [11, Table 1]. In addition, using Discrete Leja Points has also the effect of trying to minimize the sup-norm of the fundamental Lagrange polynomials ℓ i (which, as it is well-known, can be written as ratio of determinants, cf. [2]) and thus the Lebesgue constant, which is relevant to estimate (2.19). Nevertheless, it is clear from Table 1 that the bounds involving the Lebesgue constant and the condition number are much larger than (2.24) which rests only on the "active coefficients" in the differentiation process. Further numerical experiments show that, while decreasing r, for each value of |ν|, the first and third rows in Table 1  and, to measure the error of approximation, we compute the relative errors using the following bivariate test functions where f 1 is the well known Franke's function and f 3 is an oscillating function (see Figure 6) both in Renka's test set [15], whereas f 2 is obtained by a superposition of the univariate exponential with an inner product and then is constant on the parallel hyperplanes x + y = q, q ∈ R (ridge function).
For each test function we approximate D ν f (x) by  where c ν,h , with h ≤ r defined in (2.6), are the coefficients of the interpolating polynomial (2.9) at the point x. Interpolation is made at Discrete Leja Points in B r (x) for r = 1 2 , 3 8 , 1 4 , 1 8 at a sequence of degrees d. We stress that for a fixed radius r, unisolvence of interpolation is possible only for a finite number of degrees, that is until there are enough scattered points in the ball.
In the first experiment we start from 1000, 2000 and 4000 Halton and uniform random points and we set x = (0.5, 0.5) (see Figure 1). For the test function f 1 , the numerical results are displayed in Figures 2-3.
In the second experiment, we start from 2000 Halton points for the test function f 2 , again with x = (0.5, 0.5). The numerical results are displayed in Figure 4.
In the third experiment, for the test function f 3 , we start from 4000 Halton points choosing x at the center, then close to the right side and finally close to the north-east corner (see Figures  5-6). The numerical results are displayed in Figures 7, 8 and 10. We repeat the same experiments choosing x on the right side and at the north-east corner and we report the results in Figures 9 and  11.  Figure 3. As in Figure 2 starting from uniform random points. In the last experiment, for each test function f k , k = 1, 2, 3, we include a random noise in the m function values, namely (3.5)ỹ = y + U (−ε, ε),   where U (−ε, ε) denotes the multivariate uniform distribution in [−ε, ε] m . In Figure 12 we display the relative error   for the gradient at x = (0.5, 0.5), computed using 1000 Halton points with exact function values (3.2) and perturbed function values (3.6) for ε = 10 −6 . In Figure 13 we display the relative sensitivity in computing the gradient of the interpolating polynomial p under the perturbation of  the function values (ε = 10 −6 ) together with its estimate involving the stability constant of the gradient (2.33) Notice that the relative errors gep in Figure 12 and sensitivity gs in Figure 13 are of the same order of magnitude when the errors ge become negligible with respect to gs. Moreover, gse turns out to be a slight overestimate of the relative sensitivity gs. This fact, as we can see in Table 2 where x = (0.5, 0.5), is due to the relative small values of the stability constant that varies slowly while decreasing the radius r or increasing the degree of the interpolating polynomial. Finally, it is worth stressing that the stability constant (2.33) is a function of x and then it depends on the position of x in the square. More precisely, for a fixed interpolation degree d, it slowly varies except for a neighborhood of the boundary where it increases rapidly, expecially near the vertices, as can be observed in Figure 14, where for clarity we restrict the stability constant to lines. On the other hand, it can be noticed that by increasing the degree, the stability constant tends to increase, much more rapidly at the boundary. Such a behavior, that explains the worsening of the accuracy at the boundary (see  and at the vertex (see Figures 10-11) of the square, can be ascribed to the fact that the stability functions (2.34) increase rapidly near the boundary of the local interpolation domains B h (x) ∩ Ω. This phenomenon, that resembles the behavior at the boundary of Lebesgue functions of univariate interpolation at equispaced points [16], is worth of further investigation.