Multivariate polynomial interpolation on lower sets

Available online at www.sciencedirect.com

ScienceDirect Journal of Approximation Theory 177 (2014) 34–42 www.elsevier.com/locate/jat

Full length article

Multivariate polynomial interpolation on lower sets Nira Dyn a , Michael S. Floater b,∗ a School of Mathematical Sciences, Tel Aviv University, Ramat Aviv, Tel Aviv 69978, Israel b Department of Mathematics, University of Oslo, PO Box 1053, Blindern, 0316 Oslo, Norway

Received 12 April 2013; received in revised form 8 August 2013; accepted 22 September 2013 Available online 11 October 2013 Communicated by Paul Nevai

Abstract In this paper we study multivariate polynomial interpolation on lower sets of points. A lower set can be expressed as the union of blocks of points. We show that a natural interpolant on a lower set can be expressed as a linear combination of tensor-product interpolants over various intersections of the blocks that define it. c 2013 Elsevier Inc. All rights reserved. ⃝ Keywords: Multivariate polynomial interpolation; Lower sets; Tensor-product interpolants

1. Introduction Lower sets are subsets of Cartesian grids of points in Rd that admit unique polynomial interpolation from a natural, associated space of polynomials, and include rectangular and triangular grids as special cases. Such interpolation has been studied, for example, in [1,3–6, 8–13]. This natural polynomial space is the span of a collection of monomials with powers taken from a similar lower set of indices. The resulting interpolant is the “least interpolant”, introduced in [3], and studied further in [2]. A lower set of points can be expressed as the union of blocks (rectangular grids) of points. The purpose of this paper is (1) to make the observation that the interpolant can be expressed as a linear combination of tensor-product interpolants, each corresponding to an intersection of some ∗ Corresponding author.

E-mail addresses: [email protected] (N. Dyn), [email protected] (M.S. Floater). c 2013 Elsevier Inc. All rights reserved. 0021-9045/$ - see front matter ⃝ http://dx.doi.org/10.1016/j.jat.2013.09.008

N. Dyn, M.S. Floater / Journal of Approximation Theory 177 (2014) 34–42

35

of these blocks, and (2) to give an explicit formula for the linear combination, and to show how the formula simplifies in various important cases, including the 2D case, and the total degree case. It was brought to our attention, after submitting the paper, that similar relations were obtained for piecewise polynomial interpolants on sparse grids, for the numerical solution of PDE’s. In particular, our formula for the total degree polynomial interpolant is similar to the ‘combination technique’ for sparse grids (see e.g. [7]). Here is an outline of the paper. In Section 2 we introduce some notation, and discuss tensorproduct polynomial interpolation. In Section 3 we present the interpolation problem on lower sets, and give an explicit expression for the interpolants in terms of a Newton-type formula. An incremental double sum formula for these interpolants, based on blocks, is derived in Section 4, while its simplifications are obtained in Section 5 for the 2D case. For the arbitrary dimension case, a general formula in terms of tensor-product interpolants is derived in Section 6, and again simplified. The simplified formula is applied in Section 7 to interpolation by polynomials of a fixed total degree. 2. Tensor-product interpolation We consider a Cartesian grid of points X ⊂ Rd . For each j ∈ {1, . . . , d}, let x j,k , k ∈ N0 , be a sequence of distinct real points. Then the points of X are xα = (x1,α1 , x2,α2 , . . . , xd,αd ),

α ∈ Nd0 ,

where we use the multi-index notation α = (α1 , α2 , . . . , αd ) ∈ Nd0 , with |α| := α1 + · · · + αd , and the convention that α ≤ β means that α j ≤ β j for j = 1, . . . , d. Given any β ∈ Nd0 , we call the set of indices Bβ = {α ∈ Nd0 : 0 ≤ α ≤ β}, a block, and with it we can associate the set of grid points X β = {xα : α ∈ Bβ }. Note that unlike Bβ , the points of X β are not necessarily consecutive in the grid X , as illustrated in Fig. 1, showing an example of the set X (1,1) corresponding to the block B(1,1) . If we denote a monomial in Rd by x α = x1α1 · · · xdαd , for α ∈ Nd0 , then Pβ = span{x α : α ∈ Bβ } is an interpolation space for X β . In other words, for every function f : Rd → R there is a unique tensor-product polynomial p ∈ Pβ such that p(xα ) = f (xα ) for all α ∈ Bβ . One way of expressing p is the Newton form; see [8], Chap. 5. We define, for each j = 1, . . . , d, the univariate polynomials ω j,0 (y) = 1 and ω j,k (y) =

k−1  i=0

(y − x j,i ),

k ≥ 1, y ∈ R,

36

N. Dyn, M.S. Floater / Journal of Approximation Theory 177 (2014) 34–42

Fig. 1. An example of the set X (1,1) .

and the d-variate polynomial ωα (x) = ω1,α1 (x1 ) · · · ωd,αd (xd ),

α ∈ Nd0 .

We further define the divided difference ∆α,β f , for 0 ≤ α ≤ β, by the recursion ∆α,β f :=

∆α+e j ,β f − ∆α,β−e j f x j,β j − x j,α j

if α j < β j , and otherwise ∆α,β f := f (xα ). j

Here, e j ∈ Nd0 is the multi-index with ei = δi j . This kind of divided difference is the usual tensor-product one, i.e., it is applied to the function coordinate-wise. Then we can express p as  p(x) = ωα (x)∆0,α f, x ∈ Rd . (1) α∈Bβ

In later sections it will be convenient to denote the interpolant p by p(Bβ ). 3. Interpolation on lower sets We call a finite set L ⊂ Nd0 a lower set if whenever µ ∈ L and 0 ≤ α ≤ µ then α ∈ L. We call a point β ∈ L a maximal point if there is no µ ∈ L such that µ ̸= β and β ≤ µ. There is a natural space of polynomials associated with L, PL = span{x α : α ∈ L}. The monomials x α , α ∈ L, form a basis for PL , and so the dimension of PL is |L|, the cardinality of L. Corresponding to a lower set L, let Y L denote the set of points Y L = {xα : α ∈ L}. Due to the arbitrary ordering of the coordinates x j,k for each j = 1, . . . , d, the set Y L can take on different configurations. This is illustrated in Fig. 2 for the lower set L = B(1,3) ∪ B(3,1) where three different configurations of Y L are shown. The following theorem has been established in various ways in various papers in various cases. The earliest reference appears to be Kuntzmann [9] for d = 2. Here, for the sake of completeness, we prove it in the general case.

N. Dyn, M.S. Floater / Journal of Approximation Theory 177 (2014) 34–42

37

Fig. 2. Three configurations of the set Y L .

Theorem 1. There is a unique polynomial p ∈ PL that interpolates f on Y L . This polynomial can be expressed as  p(x) = ωα (x)∆0,α f, x ∈ Rd . (2) α∈L

Proof. To see that p interpolates f , suppose x = xµ in (2) for some µ ∈ L. Then, since ωα (xµ ) = 0 if α ̸≤ µ,  p(xµ ) = ωα (xµ )∆0,α f. 0≤α≤µ

But the polynomial  q(x) = ωα (x)∆0,α f, 0≤α≤µ

in view of (1), is the tensor-product interpolant to f in Pµ at the points X µ , and so p(xµ ) = q(xµ ) = f (xµ ). Having shown that p interpolates f , and since f is arbitrary, it follows from elementary linear algebra that p is also unique.  In later sections we will sometimes refer to the interpolant p as p(L). 4. Block interpolation A lower set L can, by definition, be expressed as the union of the blocks associated with its maximal points:  L= Bα , (3) α∈V

where V ⊂ L is the set of maximal points of L. Conversely, every finite union of blocks is a lower set. The intersection of blocks, on the other hand, is itself a block. In fact, the intersection of the two blocks Bα and Bβ is the block Bγ where γ = min(α, β), namely γ j = min(α j , β j ) for j = 1, . . . , d. In what follows we will show that interpolation on L in the space PL can be expressed as linear combinations of tensor-product interpolants over the intersections of the blocks Bα , α ∈ V . The main ingredient in this is the following identity. Lemma 1. If L and M are lower sets, then p(L ∪ M) = p(L) + p(M) − p(L ∩ M).

(4)

38

N. Dyn, M.S. Floater / Journal of Approximation Theory 177 (2014) 34–42

Proof. From the Newton form in Theorem 1,  p(L ∪ M)(x) = ωα (x)∆0,α f, x ∈ Rd , α∈L∪M

and so the result follows from the fact that     = + − .  α∈L∪M

α∈L

α∈M

α∈L∩M

In order to interpolate on a lower set, we can express it as the union of several blocks, and apply recursion, based on the lemma, with M a block, to build up the interpolant incrementally. From now on, when referring to a block we will sometimes write Bi when it is the ith block in a sequence of blocks. Given blocks Bi , i = 1, 2, . . . , let L r = ∪ri=1 Bi . Then Lemma 1 gives p(L n ) = p(L n−1 ) + p(Bn ) − p(L n−1 ∩ Bn ). This formula leads to a double sum formula p(L n ) =

n 

p(Bi ) −

i=1

n 

p(L i−1 ∩ Bi ).

(5)

i=2

5. Two dimensions Blocks in R2 are rectangles. Suppose B1 , . . . , Bn are n rectangles, such that Bi = Bβ i with ∈ N20 the multi-index β i = (β1i , β2i ). Let us further assume that β i ̸≤ β k for any i ̸= k. Then by reordering the rectangles, if necessary, we can assume that

βi

0 ≤ β11 < β12 < · · · < β1n ,

β21 > β22 > · · · > β2n ≥ 0.

(6)

Thus the configuration of the rectangles can be thought of as a kind of staircase shape; see Fig. 3, in which n = 5 and β 1 = (2, 10), β 2 = (3, 8), β 3 = (5, 5), β 4 = (8, 4), β 5 = (9, 2). Theorem 2. In the staircase configuration the interpolant reduces to p(L n ) =

n 

p(Bi ) −

i=1

n 

p(Bi−1 ∩ Bi ).

(7)

i=2

Proof. This follows from the double sum formula (5), since by (6), L i−1 ∩ Bi = {α : 0 ≤ α ≤ (β1i−1 , β2i )} = Bi−1 ∩ Bi .  6. Arbitrary dimension Consider now general dimension d. Repeated use of the double sum formula (5) leads to the following general formula. Theorem 3. Let L n be the union of the n blocks B1 , B2 , . . . , Bn . Then p(L n ) =

n  (−1)k−1 k=1

 1≤i 1
p(Bi1 ∩ Bi2 ∩ · · · ∩ Bik ).

(8)

N. Dyn, M.S. Floater / Journal of Approximation Theory 177 (2014) 34–42

39

Fig. 3. A lower set in N20 . The points β i , i = 1, . . . , 5, are shown as black circles, the points (β1i−1 , β2i ), i = 2, . . . , 5, as white circles.

Proof. The proof is by induction on n. The formula trivially holds when n = 1. For n > 1, since L i−1 ∩ Bi = (B1 ∩ Bi ) ∪ · · · ∪ (Bi−1 ∩ Bi ),

i = 2, . . . , n,

we conclude that L i−1 ∩ Bi is the union of i − 1 blocks, and i − 1 ≤ n − 1. Thus we can apply the induction hypothesis to the double sum formula (5) to deduce that p(L n ) =

n 

p(Bi ) −

n  i−1  (−1)k−1

i=1

i=2 k=1



p(Bi1 ∩ · · · ∩ Bik ∩ Bi ).

1≤i 1 <···
The second sum above can be rewritten as −

n  (−1)k−1 k=2



p(Bi1 ∩ · · · ∩ Bik ),

1≤i 1 <···
which when subtracted from the first sum gives (8).



With the shorthand p = p(L n ) and pi1 ,...,ik := p(Bi1 ∩ · · · ∩ Bik ), the formulas in the cases n = 2, 3, 4 are p = ( p1 + p2 ) − p12 , p = ( p1 + p2 + p3 ) − ( p12 + p13 + p23 ) + p123 , p = ( p1 + p2 + p3 + p4 ) − ( p12 + p13 + p14 + p23 + p24 + p34 ) + ( p123 + p124 + p134 + p234 ) − p1234 . If n > d some of the block intersections in the sum (8) are over more than d blocks, in which case they reduce to intersections of at most d of the blocks. Thus, in general, many of the terms in the sum (8) are repeated. It follows from Theorem 3 that there must be integer coefficients cα ,

40

N. Dyn, M.S. Floater / Journal of Approximation Theory 177 (2014) 34–42

α ∈ L, such that  p(L) = cα p(Bα ).

(9)

α∈L

We now use Lemma 1 to derive a general formula for the coefficients. For any set L ⊂ Nd0 , we denote by χ (L) : Nd0 → {0, 1} its characteristic function defined by  1 if α ∈ L; χ (L)(α) = 0 otherwise. Theorem 4. The coefficients in (9) are given by  cα = (−1)|ϵ| χ (L)(α + ϵ), α ∈ L .

(10)

ϵ∈{0,1}d

We note that the sum in (9) can be extended to all α ∈ Nd0 if necessary because if α ̸∈ L then cα = 0. This follows from the fact that if α ̸∈ L then α + ϵ ̸∈ L for all ϵ ∈ {0, 1}d in the sum in (10). This fact will facilitate the proof of the theorem by induction. Proof. The proof of the formula is by induction on n ≥ 1, the number of blocks defining L. Suppose first that n = 1 and that L = Bβ for some β ∈ Nd0 . For any α ∈ Bβ define γ ∈ {0, 1}d by γ j = 1 if α j < β j and γ j = 0 if α j = β j , j = 1, . . . , d. Then for ϵ ∈ {0, 1}d , α + ϵ ∈ L if and only if 0 ≤ ϵ ≤ γ , and so  cα = (−1)|ϵ| . 0≤ϵ≤γ

Therefore, if α = β, γ = 0, and cα = (−1)0 = 1. On the other hand, if α ̸= β, γ ̸= 0, and, letting k = |γ |,   k  i k cα = (−1) = (1 − 1)k = 0, i i=0 which proves (10) for n = 1. Otherwise, suppose that L = L n ∪ B, where L n is the union of n blocks and B is a block. We can assume by the induction hypothesis that the formula holds for the union of at most n blocks. By Lemma 1, p(L) = p(L n ) + p(B) − p(L n ∩ B), and so, by applying (9)–(10) to the three terms on the right hand side, it follows that (9) holds for L with  cα = (−1)|ϵ| (χ (L n ) + χ (B) − χ (L n ∩ B))(α + ϵ). ϵ∈{0,1}d

This reduces to (10) because χ (L n ) + χ (B) − χ (L n ∩ B) = χ (L).



N. Dyn, M.S. Floater / Journal of Approximation Theory 177 (2014) 34–42

41

Let 1d := (1, 1, . . . , 1) ∈ Nd0 . We call a point α ∈ L a boundary point if α + 1d ̸∈ L, and an interior point otherwise. Corollary 1. For α ∈ L an interior point, cα = 0. Proof. Because L is a lower set, if α + 1d belongs to L, so does α + ϵ for any ϵ ∈ {0, 1}d and so (10) implies d     d |ϵ| (−1) = cα = (−1)i = 0.  i d i=0 ϵ∈{0,1}

It follows that a point α ∈ L with a non-zero coefficient cα must be both an intersection of some of the blocks defining L and a boundary point of L. This means that in the staircase configuration (6), where d = 2, α must belong to either {β i : i = 1, . . . , n} or {(β1i−1 , β2i ) : i = 2, . . . , n}. Since Theorem 4 gives cα = 1 in the first case, and cα = −1 in the second, it thus recovers Theorem 2. 7. Interpolation of total degree As an example of interpolation on a lower set, consider polynomial interpolation of total degree m ≥ 0, so that L = {α ∈ Nd0 : |α| ≤ m}. Then the set of maximal points is V = {α ∈ Nd0 : |α| = m}, and L is the union of the blocks Bα with |α| = m. In the special case d = 2, the interpolant is a staircase interpolant (7),   p(L) = p(Bα ) − p(Bα ). |α|=m

|α|=m−1

For general d, we apply Theorem 4 and Corollary 1. Suppose α ∈ L with |α| ≤ m − d. Then |α + 1d | ≤ m, and so α + 1d ∈ L and α is an interior point, and by Corollary 1, cα = 0. Suppose otherwise that |α| = m − k for some k ∈ {0, 1, . . . , d − 1}, in which case it is a boundary point. Then Theorem 4 gives   k   d cα = (−1)|ϵ| = (−1)i , i d i=0 ϵ∈{0,1} :|ϵ|≤k

and so cα = 1 +

      k  d −1 d −1 d −1 (−1)i + = (−1)k . i −1 i k i=1

42

N. Dyn, M.S. Floater / Journal of Approximation Theory 177 (2014) 34–42

So, for example, the interpolant for d = 3 is   p(Bα ) + p(L) = p(Bα ) − 2 |α|=m

|α|=m−1



p(Bα ),

|α|=m−2

where, as usual, an empty sum has value 0. References [1] V. Barthelmann, E. Novak, K. Ritter, High dimensional polynomial interpolation on sparse grids, Adv. Comput. Math. 12 (2000) 273–288. [2] C. de Boor, On the error in multivariate polynomial interpolation, Appl. Numer. Math. 10 (1992) 297–305. [3] C. de Boor, A. Ron, On multivariate interpolation, Const. Approx. 6 (1990) 287–302. [4] A. Chkifa, A. Cohen, C. Schwab, High-dimensional adaptive sparse polynomial interpolation and applications to parametric PDEs, Found. Comput. Math., in press. [5] M. Gasca, T. Sauer, On the history of multivariate polynomial interpolation, J. Comput. Appl. Math. 122 (2000) 23–35. [6] M. Gasca, T. Sauer, Polynomial interpolation in several variables, Adv. Comput. Math. 12 (2000) 377–410. [7] M. Hegland, The combination technique and some generalisations, Linear Algebra Appl. 420 (2007) 249–275. [8] E. Isaacson, H.B. Keller, Analysis of Numerical Methods, Dover, 1994. [9] J. Kuntzmann, M´ethodes num´eriques: interpolation, d´eriv´ees, Dunod, Paris, 1959. [10] G.G. Lorentz, R.A. Lorentz, Solvability problems of bivariate interpolation I, Const. Approx. 2 (1986) 153–169. [11] G. M¨uhlbach, On multivariate interpolation by generalized polynomials on subsets of grids, Computing 40 (1988) 201–215. [12] T. Sauer, Lagrange interpolation on subgrids of tensor product grids, Math. Comp. 73 (2004) 181–190. [13] H. Werner, Remarks on Newton type multivariate interpolation for subsets of grids, Computing 25 (1980) 181–191.