Reproduction of Polynomials by Radial Basis Functions

Robert Schaback

Abstract For radial basis function interpolation of scattered data in IR^{d}, the approximative reproduction of high-degree polynomials is studied. Results include uniform error bounds and convergence orders on compact sets.

§1. Introduction

We consider interpolation of real-valued functions f defined on a set Ω ⊆ IR^{d}, d ≥ 1. These functions are interpolated on a set X := {x_{1}, …, x_{N_{X}} } of N_{X} ≥ 1 pairwise distinct points x_{1}, …, x_{N_{X}} in Ω. Interpolation is done by linear combinations of translates Φ(x − x_{j}) of a single continuous real-valued function Φ defined on IR^{d}. For various reasons it is sometimes necessary to add the space IP_{m}^{d} of d–variate polynomials of order not exceeding m to the interpolating functions. Interpolation is uniquely possible under the requirement

If p ∈ P_{m}^{d} satisfies p(x_{i}) = 0 for all x_{i} ∈ X then p = 0,  (1)

and if Φ is conditionally positive definite of order m (see e.g. [8]):

Definition 1. A function Φ : IR^{d} → IR with Φ(−x) = Φ(x) is conditionally positive definite of order m on IR^{d}, if for all sets X = {x_{1}, …, x_{N_{X}} } ⊂ IR^{d} with N_{X} distinct points and all vectors α := (α_{1}, …, α_{N_{X}}) ∈ IR^{N_{X}} with

∑_{j=1}^{N_{X}} α_{j}p(x_{j}) = 0 for all p ∈ IP_{m}^{d}  (2)

the quadratic form ∑_{j,k=1}^{N_{X}} α_{j}α_{k}Φ(x_{j} − x_{k}) attains nonnegative values and vanishes only if α = 0.

We list a few examples where Φ(x) := ϕ(∥x∥_{2}) is truly radial:

Multiquadrics ϕ(r) = (c^{2} + r^{2})^{β/2} for β ∈ IR_{>−d} \ 2ZZ and 2m > β,
Thin-plate splines ϕ(r) = r^{β} for β ∈ IR_{>0} \ 2ZZ and 2m > β,
Thin-plate splines ϕ(r) = (−1)^{β/2+1}r^{β} log r for β ∈ 2IN, 2m > β,

Sobolew splines ϕ(r) = \frac{2π^{d/2}}{Γ(k)} K_{k−d/2}(r)(r/2)^{k−d/2} for k > d/2 and m ≥ 0, using the Macdonalds (or modified spherical Bessel) function K_{ν}, Gaussians ϕ(r) = e^{−cr^{2}} for c > 0 and m ≥ 0.

There has been some discussion concerning the relative merits of these radial basis functions. Gaussians and other integrable radial basis functions were said to have “a severe disadvantage”, because they are not exact on constant functions [1, p. 93]. The other choices, especially multiquadrics, were said to be “far superior”.

On the other hand, Gaussians have an even better convergence behavior than multiquadrics, when used on scattered data (Madych/Nelson [6]) for c fixed. Thus it is questionable whether the advantage of a method must be necessarily tied to its polynomial reproduction properties. The link between convergence orders and polynomial reproduction properties is clearly visible in theories on grids that make use of Strang and Fix conditions (see the review by W. Light [3]), but it is not at all evident for interpolation of scattered data on irregular domains. Exact reproduction may well be replaced by high-order approximation without any loss in theory and numerical practice.

Numerical results indicate that radial basis function interpolation behaves very well indeed on data from polynomials, provided that the centers are dense enough in and around the domain where the interpolant is evaluated. To add further support to this statement, this paper investigates the polynomial reproduction properties of radial basis function interpolation in the scattered data case.

§2. Spaces for Radial Basis Functions

Each conditionally positive definite function Φ allows two constructions of an inner-product space of functions. Both constructions are based on ideas of Madych and Nelson [4, 5], but we adopt the terminology of [8] here and omit details.

The algebraic approach introduces a space F_{Φ} by direct reference to the conditional positive definiteness of Φ of order m. Functions of the form

f_{α}(x) = ∑_{j=1}^{N} α_{j}Φ(x − x_{j}),  (3)

with (2) are in F_{Φ} and have a norm ∥f_{α}∥_{Φ}^{2} = ∑_{j,k=1}^{N} α_{j}α_{k}Φ(x_{j} − x_{k}). The Fourier transform approach makes use of the fact that all of the radial basis functions given in the introduction have a generalized Fourier transform in the sense of tempered distributions, and the latter coincides with a positive and continuous function φ on IR^{d} \ {0} in the sense of Jones [2]. In this paper we assume the function φ to be continuous and positive on IR^{d} \ {0}, and to satisfy the conditions

\begin{aligned}|φ(ω)| & ≥ C · ∥ω∥^{−d−s_{0}} , ω ∈ U, \\ |φ(ω)| & ≥ c · e^{−γ∥ω∥^{2}} , ω ∉ U, \\ |φ(ω)| & ≤ c · ∥ω∥^{−d−s_{∞}} , ω ∉ U. \end{aligned}  (4)

on a neighborhood U = K_{R}(0) = {x | ∥x∥_{2} < R, x ≠ 0} of the origin in IR^{d} \ {0}. Examples are compiled conveniently in [1] and [8], and we add φ(ω) = (1 + ∥ω∥_{2}^{2})^{−k} for Sobolew radial basis functions of order k. This function, for instance, satisfies (4) for s_{0} = −d, s_{∞} = 2k − d, and an arbitrary γ > 0 related to the radius R of the ball U.

The Fourier transform approach introduces the scalar product

(f, g)_{φ} := (2π)^{−d} ∫_{IR_{d}} \frac{f̂(ω)\overline{ĝ(ω)}}{φ(ω)} dω

and the space

F_{φ} := {f : IR^{d} → IR | f̂ ∈ L_{1}(IR^{d}), ∥f∥_{φ} < ∞},  (5)

that occurs in [9] and dates back to a similar but different space defined by Madych and Nelson [5], using a nonstandard theory of generalized functions. Functions of the form (3) with (2) for some m > s_{0}/2 are also in F_{φ} and have a norm

∥f_{α}∥_{φ}^{2} = (2π)^{−d} ∫_{IR_{d}} \frac{|f̂_{α}(ω)|^{2}}{φ(ω)} dω = ∥f_{α}∥_{Φ}^{2}.

We finally note that the Fourier transform approach adds IP_{q}^{d} to F_{φ} to perform interpolation, where q > s_{0}/2 holds for the smallest possible s_{0} in (4). Reasons for this will become apparent later.

§3. Approximation of Polynomials

We shall now approximate polynomials by functions in the space F_{φ} of (5). First we represent the delta functional δ by the sequence of good

functions γ_{k}(ω) = (\frac{k}{π} )^{d/2} e^{−k∥ω∥^{2}} in the sense of Jones [2]. Then the generalized derivative D^{α}δ for a multi-index α is represented by the sequence D^{α}γ_{k}, and the generalized Fourier transform of this sequence represents the polynomial D̂^{α}δ(x) = (+ix)^{α}. We now fix z and consider the test function f_{k,α,z}(x) := (i(x + z))^{α}γ̂_{k}(x + z) = D̂^{α}γ_{k}(x + z) that satisfies δ(f_{k,α,z}) = f_{k,α,z}(0) = (iz)^{α}γ̂_{k}(z) = (iz)^{α}e^{−∥z∥^{2}/4k} and converge to (iz)^{α} for k → ∞, proving that the δ functional is bounded on these functions when k → ∞. We shall make use of this fact later. The functions f_{k,α,z}(x) converge uniformly to the polynomial (i(x + z))^{α} on compact sets Ω ⊂ IR^{d}, if k tends to infinity. This is due to

γ̂_{k}(u) = e^{−∥u∥^{2}/4k} = 1 − \frac{∥u∥^{2}}{4k} + \frac{∥u∥^{4}}{16k^{2}} ± ⋯ ≥ 1 − \frac{∥u∥^{2}}{4k}

for all u ∈ IR^{d}. More precisely, for each α and each compact set Ω ⊂ IR^{d} there is a constant c_{0}, independent of k ∈ IN and x ∈ Ω, such that

|(ix)^{α} − f_{k,α,0}(x)| ≤ \frac{∥x∥^{2}}{4k}  (6)

holds for all x ∈ Ω and all k ∈ IN. In what follows, we shall use the symbols c, c_{0}, c_{1}, … to denote generic constants which may differ by certain factors, and we just state the variables that the constants do not depend on. Such “constants” may, however, depend on other variables that occur in the context, and these are mainly α, Φ, and Ω.

Now we assert that all the functions f_{k,α,z}(x) are in the spaces F_{φ} for all radial basis functions that satisfy the basic assumptions (4).

Lemma. Let the given radial basis function satisfy (4). Then there are positive constants c_{1} and c_{2}, independent of k and z, such that

|f_{k,α,z}|_{φ}^{2} = ∫ \frac{|f̂_{k,α,z}(ω)|^{2}}{φ(ω)} dω ≤ c_{1} · k^{|α|−s_{0}/2}  (7)

holds for all z ∈ IR^{d} and all k ≥ c_{2}.

Proof: By elementary calculations,

\begin{aligned}f̂_{k,α,z}(ω) & = (2π)^{d}e^{i⟨z,ω⟩}D^{α}γ_{k}(ω) \\ & = (2π)^{d}e^{i⟨z,ω⟩} (\frac{k}{π} )^{d/2} D^{α}e^{−k∥ω∥^{2}} \end{aligned}

and we have to check these derivatives for k → ∞. Using Hermite polynomials H_{n}(x) defined here as H_{n}(x) = (−1)^{n}e^{x^{2}/2}D^{n}e^{−x^{2}/2}, we first treat the univariate case by

\frac{d^{n}}{dx^{n}} e^{−x^{2}/2} = (−1)^{n}H_{n}(x)e^{−x^{2}/2} = (\sqrt{2k})^{−n} \frac{d^{n}}{dω^{n}} e^{−kω^{2}}

with x = ω\sqrt{2k}, giving

D^{n}e^{−kω^{2}} = (\sqrt{2k})^{n}(−1)^{n}H_{n}(ω\sqrt{2k})e^{−kω^{2}}.

The multivariate case is a product of univariate cases, because

D^{α}e^{−k∥ω∥^{2}} = ∏_{j=1}^{d} \frac{∂^{α_{j}}}{∂ω{α_{j} \atop j}} e^{−kω_{j}^{2}}

= ∏_{j=1}^{d}(\sqrt{2k})^{α_{j}}(−1)^{α_{j}}H_{α_{j}}(ω_{j}\sqrt{2k})e^{−kω_{j}^{2}}

= (−\sqrt{2k})^{|α|}e^{−k∥ω∥^{2}}H_{α}(ω\sqrt{2k}),

where we defined a multivariate Hermite polynomial H_{α}(z) = ∏_{j=1}^{d} H_{α_{j}}(z_{j}) for z ∈ IR^{d}, α ∈ ZZ_{≥0}^{d}. Using (4), the integral of (7) on U is bounded by

k^{d+|α|} ∫_{U} ∥ω∥^{d+s_{0}} · e^{−2k∥ω∥^{2}} · |H_{α}(ω\sqrt{2k})|^{2}dω

up to constants which are independent of k and z. By transformation u = ω\sqrt{2k} we get the bound

k^{d+|α|}k^{−d−s_{0}/2} ∫_{∥u∥≤R\sqrt{2k}} ∥u∥^{d+s_{0}}e^{−∥u∥^{2}}|H_{α}(u)|^{2}du ≤ c · k^{|α|−s_{0}/2}.

On the complement of U and for 2k > γ we can use the bound |H_{n}(x)| ≤ c(n) · |x|^{n} for Hermite polynomials H_{n}(x) for large arguments, and we find

\begin{aligned}& k^{2|α|}k^{d} ∫_{R}^{∞} r^{d−1}e^{γr^{2}}r^{2|α|}e^{−2kr^{2}}dr \\ ≤ & c · \frac{k^{2|α|}k^{d}}{(2k − γ)^{(d+2|α|)/2}} ∫_{R\sqrt{2k−γ}}^{∞} s^{d−1+2|α|}e^{−s^{2}}ds \end{aligned}

which vanishes exponentially for k → ∞. Thus we finally have (7) with constants that depend only on α, Ω, d, φ, and R, but not on z or k.

The above argument proves that for s_{0} > 0 the delta functional is not continuous on F_{φ}, because one can pick any α with 0 ≤ |α| < s_{0}/2 to get

lim_{k→∞} \frac{|δ(f_{k,α,z})|}{∥f_{k,α,z}∥_{φ}} = ∞

for all z ∈ IR^{d} \ {0}. More generally,

Theorem 1. A functional of the form

ϵ_{x}f = f(x) − ∑_{j=1}^{N} u_{j}(x)f(x_{j})  (8)

is bounded on F_{φ}, iff it vanishes on all polynomials of degree less than s_{0}/2.

Proof: The functional is continuous if it vanishes on P_{q}^{d} (see [5], [9]). To prove the converse, assume ϵ_{x} to be bounded on F_{φ}. Then

lim_{k→∞} |ϵ_{x}f_{k,α,z}| ≤ c(x) · lim_{k→∞} ∥f_{k,α,z}∥_{φ} = 0

for all |α| < s_{0}/2 and all z. On the other hand, by the mollification and (6),

lim_{k→∞} ϵ_{x}f_{k,α,z} = (i(x + z))^{α} − ∑_{j} u_{j}(x)(i(x_{j} + z))^{α}

for all |α| < s_{0}/2 and all z. But this means that ϵ_{x} must vanish on all polynomials of degree < s_{0}/2 if it should be continuous on F_{φ}.

Theorem 1 throws some light on the admissible values of q. It is possible to define the space F_{φ} of (5) with values of q that do not satisfy q > s_{0}/2, but this inequality is required to let the usual radial basis function interpolants be contained in F_{φ}. At the same time q > s_{0}/2 ensures that functionals of the form (8) are continuous on F_{φ} if they are exact on P_{q}^{d}, and Theorem 1 proves a partial result for the converse assertion, making reproduction of low-degree polynomials a necessary condition for boundedness of the error functional. The consequence is that one cannot use values of q that do not satisfy q > s_{0}/2 without throwing the interpolants out of F_{φ} or making the error functional discontinuous on this space.

Due to Theorem 1, there are two situations:

1. For s_{0} ∉ 2IN, continuity of the error functional is equivalent to exactness on IP_{q}^{d} for the smallest q > s_{0}/2. This value of q yields the smallest possible error bound.

2. For s_{0} = 2k ∈ 2IN the necessary and sufficient conditions differ. Exactness on IP_{q}^{d} for any q > k is sufficient for continuity, and exactness on IP_{q}^{d} for all q ≤ k is necessary. This gap occurs for thin–plate splines ϕ(r) = r^{2k} log r.

Furthermore the proof of Theorem 1 now implies that interpolants from multiquadrics Φ(∥x∥) = \sqrt{c^{2} + ∥x∥^{2}} without addition of constant polynomials will not lead to a bounded error functional on the space F_{φ}.

On the space F_{φ} defined for thin-plate splines Φ(x) = ∥x∥^{3}, interpolation by multiquadrics will have an unbounded error, if linear polynomials are not added to the multiquadrics.

However, there can possibly exist larger and more interesting spaces than F_{φ} that still contain the interpolants and may lead to error estimates in weaker norms. These spaces could possibly allow smaller values of q, but they are still to be found.

§4. Approximative Polynomial Reproduction

We use results concerning the approximation order of radial basis function interpolants [4, 5, 8]. For a given radial basis function Φ with a generalized Fourier transform coinciding with φ on IR^{d}\{0} satisfying (4), a given compact set Ω ⊂ IR^{d} and a given constant ρ > 0 there are positive constants c_{3} and c_{4} such that for all finite sets X = {x_{1}, …, x_{N_{X}} } ⊂ IR^{d} of centers and all points x ∈ Ω with

h_{X,ρ}(x) := sup_{∥y−x∥_{2}≤ρ} min_{1≤j≤N_{X}}∥y − x_{j}∥_{2} ≤ c_{3}  (9)

the interpolation error of x is bounded by c_{4}·∥f∥_{φ}h_{X,ρ}^{s_{∞}/2}(x), where s_{∞} is determined by φ at infinity via (4). Here, the constants c_{3} and c_{4} are independent of x, X, and h_{X,ρ}(x). With this (7) yields

|f_{k,α,z}(x) − ∑_{j} u_{j}^{∗}(x)f_{k,α,z}(x_{j})|^{2} ≤ c · k^{|α|−s_{0}/2} · h_{X,ρ}^{s_{∞}}(x)  (10)

for all z ∈ IR^{d}, where u_{j}^{∗}(x) are optimal for Φ, and where c is independent of k, z, X, and x. This estimate interests us for |α| ≥ s_{0}/2, in order to prove something about the behaviour of optimal interpolation processes on high– degree polynomials. Defining s_{X,k,α}(x) := ∑_{j=1}^{N_{X}} u_{j}^{∗}(x)f_{k,α,0}(x_{j}) as a quasiinterpolant that is an interpolant to a mollification of the polynomial (ix)^{α}, we get

|f_{k,α,0}(x) − s_{X,k,α}(x)|^{2} ≤ c_{5}k^{|α|−s_{0}/2} · h_{X,ρ}^{s_{∞}}(x)

and add (6) to yield

|x^{α} − (−i)^{|α|}s_{X,k,α}(x)| ≤ \frac{c_{6}}{k} + c_{7}k^{|α|/2−s_{0}/4}h_{X,ρ}^{s_{∞}/2}(x)

with constants that are independent of k, X, and h_{X,ρ}(X).

This bounds the approximation error of a polynomial in the space of radial basis function interpolants. Special approximations of the monomial x^{α} are obtained here via interpolation of the test functions f_{k,α,0} that uniformly converge to x^{α} on compact sets at the rate O(k^{−1}) for k → ∞. The error depends on the centers X and their local density h_{X,ρ}(x) at x, of course, but also on the free parameter k that controls the mollification.

We now tie k to h_{X,ρ}(x) by making the two terms in the error bound approximately equal, and we are interested in the convergence rate of the error for centers X with densities h_{X,ρ}(x) tending to zero. This is done by choosing k ≈ h_{X,ρ}(x)^{−s_{∞}/(2+|α|−s_{0}/2)} which tends to infinity for h_{X,ρ}(x) tending to zero, whenever the inequalities |α| > s_{0}/2 − 2, s_{∞} > 0 hold. Under these conditions, we finally have

|x^{α} − (−i)^{|α|}s_{X,k,α}(x)| ≤ c_{7}(h_{X,ρ}(x))^{s_{∞}/(2+|α|−s_{0}/2)}  (11)

for all sets X of centers and all points x ∈ Ω satisfying (9), and with a constant c_{7} that is independent of x, X, and h_{X,ρ}(x). From (11) we can read off a number of statements:

Theorem 2. In the space F_{φ} generated by any of the usual radial basis functions Φ, polynomials are dense with respect to uniform convergence on compact sets. Polynomials of degree less than s_{0}/2 must necessarily be exactly reproduced to make error functionals of the form (8) continuous on F_{φ}. Polynomials of degree n ≥ s_{0}/2 can be approximated on compact sets by interpolants of certain mollifications of themselves with a guaranteed uniform convergence rate that depends on the radial basis function and on the local density h_{X,ρ} of the set X of centers. For h_{X,ρ}(x) → 0 the convergence order in terms of h_{X,ρ} is s_{∞}/(2 + n − s_{0}/2) with the following special implications:

1. For multiquadrics, inverse multiquadrics, and Gaussians with fixed parameters c, the convergence order is arbitrarily large.

2. For thin-plate splines with parameter β, the order is at least β/(2 + n − β/2).

3. For splines in Sobolew spaces of order k, the approximation order is at least k/(2 + n − k/2).

Proof: For the last three statements we simply use the appropriate values of s_{0} and s_{∞}.

§5. References

1. Buhmann, M. D., Multivariable interpolation using radial basis functions, Ph. D. dissertation, Cambridge 1989

2. Jones, D.S., The Theory of Generalised Functions, Cambridge University Press, 1982 .

3. Light, W. , Some aspects of radial basis function approximation, in: Approximation Theory, Spline Functions and Applications, S.P. Singh (ed.), Kluwer Adademic Publishers, Dordrecht, 1992.

4. Madych, W. R., and Nelson, S. A., Multivariate Interpolation and conditionally positive definite functions, Approx. Th. Appl. 4 (1988), 77–89

5. Madych, W. R., and Nelson, S. A., Multivariate Interpolation and conditionally positive definite functions II, Math. Comp. 54 (1990), 211–230

6. Madych, W. R., and Nelson, S. A., Bounds on Multivariate Polynomials and Exponential Error Estimates for Multiquadric Interpolation, in J. Approx. Th. 70 (1992) 94—113

7. Micchelli, C. A., Interpolation of scattered data: distance matrices and conditionally positive definite functions, Constr. Approx. 2 (1986), 11–22

8. Schaback, R., Comparison of radial basis function interpolants, in Multivariate Approximations: From CAGD to Wavelets, K. Jetter and F. Utereas (eds.), World Scientific Publishing Co., 1993

9. Wu, Zong-min, and R. Schaback, Local Error Estimates for Radial Basis Function Interpolation IMA J. of Numer. Anal., 13 (1993), 13–27

Help in proofreading was provided by Armin Iske and Marko Weinrich.