Local Error Estimates for Radial Basis Function Interpolation of Scattered Data

Zong–min Wu and Robert Schaback

AMS classifications: 41A05, 41A63, 41A25, 65D05

Keywords: Multivariate interpolation, orders of convergence, Kriging, Multiquadrics, thin plate splines.

Abstract

Introducing a suitable variational formulation for the local error of scattered data interpolation by radial basis functions ϕ(r), the error can be bounded by a term depending on the Fourier transform of the interpolated function f and a certain “Kriging function”, which allows a formulation as an integral involving the Fourier transform of ϕ. The explicit construction of locally well–behaving admissible coefficient vectors makes the Kriging function bounded by some power of the local density h of data points. This leads to error estimates for interpolation of functions f whose Fourier transform f̂ is “dominated” by the nonnegative Fourier transform ψ̂ of ψ(x) = ϕ(∥x∥) in the sense ∫ |f̂|^{2}ψ̂^{−1}dt < ∞. Approximation orders are arbitrarily high for interpolation with Hardy multiquadrics, inverse multiquadrics and Gaussian kernels. This was also proven in recent papers by Madych and Nelson, using a reproducing kernel Hilbert space approach and requiring the same hypothesis as above on f̂, which limits the practical applicability of the results. This work uses a different and simpler analytic technique and allows to handle the cases of interpolation with ϕ(r) = r^{s} for s ∈ IR, s > 1, s ∉ 2IN, and ϕ(r) = r^{s} log r for s ∈ 2IN, which are shown to have accuracy O(h^{s/2}).

1 Introduction

Radial basis function interpolation to scattered data (x_{i}, f_{i}) ∈ IR^{n+1} for pairwise distinct points (“centres”) x_{1},…, x_{M} ∈ IR^{n} uses a function ϕ : IR_{≥0} → IR and the space IP_{q} of polynomials on IR^{n} with total order not exceeding q to construct the interpolant

s(x) := ∑_{i=1}^{M} a_{i}ϕ(∥x − x_{i}∥) + ∑_{i=1}^{Q} b_{i}p_{i}(x)  (1.1)

via the linear system

\begin{matrix}∑_{i=1}^{M} a_{i}ϕ(∥x_{j} − x_{i}∥) & + & ∑_{i=1}^{Q} b_{i}p_{i}(x_{j}) = f_{j}, & 1 ≤ j ≤ M \\ ∑_{j=1}^{M} a_{j}p_{i}(x_{j}) = 0, & & & 1 ≤ i ≤ Q = \binom{q+n−1}{n}, \end{matrix}  (1.2)

0File : usr/nam/rschaba/tex/fertig/wu rs 1 mod 2.tex, TEXed February 27, 1992, Status : 2nd revision

where p_{1},…, p_{Q} is a basis of IP_{q}. For a wide choice of functions ϕ and polynomial orders q, including the case q = Q = 0, the nonsingularity of the (M + Q) × (M + Q) system (1.2), written as

( \begin{matrix}A & P \\ P^{T} & 0 \end{matrix} ) ( \begin{matrix}a \\ b \end{matrix} ) = ( \begin{matrix}f \\ 0 \end{matrix} )  (1.3)

in matrix notation, has been established by Micchelli [8] and Powell [12]. Following [8], we assume F (r) = ϕ(\sqrt{r}) to be conditionally strictly positive definite of order q, which implies that A is positive definite on the subset of vectors u ∈ IR^{M} satisfying P^{T} u = 0. We handle the case of conditionally strictly negative definite functions ϕ(\sqrt{r}) by going over to the function −ϕ(\sqrt{r}). For q > 0 the additional condition

p(x_{j}) = 0, 1 ≤ j ≤ M for p ∈ IP_{q} implies p = 0  (1.4)

on q and the positions of data points is required. The numerical problems of (1.3) can be overcome by preconditioning methods (see Dyn, Levin and Rippa [5], and Dyn [4]), making the radial basis function approach a promising tool for multivariate interpolation.

We consider interpolation of values f_{i} = f(x_{i}) of a smooth function f on a domain Ω ⊆ IR^{n}. There are no further conditions on Ω. The sampling points x_{i} are allowed to be irregularly distributed over Ω under the restriction (1.4). For local error estimation, we measure the “density” of centres x_{j} near some x ∈ Ω by

h_{ρ}(x) := max_{y∈K_{ρ}(x)} min_{1≤j≤M} ∥y − x_{j}∥

for some fixed ρ > 0 , where K_{ρ}(x) = {y ∈ IR^{n} | ∥x − y∥ ≤ ρ}, using the Euclidean norm ∥.∥ (see also [9], [10], [11]).

Our goal is to prove local error bounds of the following form: Given constants ρ ∈ IR_{>0}, q ∈ IN_{≥0}, a radial basis function ϕ and a certain function space F_{ϕ} to be described below, we want to show the existence of positive constants k ≥ m, k, m ∈ IN, h_{0}, C ∈ IR such that for any distribution of points x_{j} ∈ IR^{n}, 1 ≤ j ≤ M, any function f from F_{ϕ} and any point x ∈ Ω with h_{ρ}(x) < h_{0} the inequality

|s^{(μ)}(x) − f^{(μ)}(x)| ≤ c_{f} · C · h_{ρ}^{k−|μ|}(x)  (1.5)

holds for the error and its μ–th derivatives for 0 ≤ |μ| ≤ m, where the constant c_{f} depends on f, ϕ, and F_{ϕ} only. Then we call (1.5) a local error bound of order k. Here and in the sequel we use the standard multi–index notation with |μ| := ∑_{i} μ_{i} for μ ∈ IN^{n}.

Convergence of interpolation on regular grids has been studied extensively by Buhmann (see e.g: [1] for a comprehensive treatment) and others. For the case of scattered data Duchon [3] treated the thin–plate spline case ϕ(r) = r^{2} log r, while Jackson [6] proved a general, but non–quantitative convergence result. The dissertation of Wu [13] of 1986 contained a rather general Hilbert space theory for Kriging and related radial basis function methods in the case of a single variable, including error estimates and convergence results. Recently, Madych and Nelson [9], [10] developed a reproducing kernel Hilbert space approach to get error bounds of arbitrarily high order for multiquadrics and inverse multiquadrics, i.e.: for ϕ(r) = (c^{2} + r^{2})^{s/2} with c ∈ IR_{>0}, s ∈ IR, s ∉ 2IN, s > −n and q > s/2 if s > 0.

This paper is based on the approach of Wu [13] and generalizes it to the multivariate case. There are some background connections to the work of Madych and Nelson which will be explained at the appropriate places. Like Madych and Nelson we require the restrictive condition ∫ |f̂|^{2}ψ̂^{−1}dt < ∞ for ψ(x) := ϕ(∥x∥) and show that k in (1.5) can be arbitrarily large for Gaussians ϕ(r) = exp(−cr^{2}), multiquadrics ϕ(r) = \sqrt{r^{2} + c^{2}}, and inverse multiquadrics ϕ(r) = 1/\sqrt{r^{2} + c^{2}}, c > 0. For these cases our error bounds are equivalent to those of Madych and Nelson. But our technique also yields the convergence order k = s/2 for m < s/2 < q in case of the radial basis functions ϕ(r) = r^{s} with s ∈ IR_{>0}, s ∉ 2IN, and ϕ(r) = r^{s} log r with s ∈ 2IN. Compared to the work of Madych and Nelson the methods of this paper are somewhat simpler and more direct, but do not provide or require additional information about the Hilbert space background. They generally imply that the order k in (1.5) is attained, if the Fourier transform ψ̂ satisfies

0 < ψ̂(t) ≤ O(∥t∥^{−n−2k})

for ∥t∥ → ∞.

2 Variational formulation

Since the equations (1.2) are solvable under the hypotheses of the preceding section, there is a Lagrange–type representation

s(x) = ∑_{j=1}^{M} f(x_{j})u_{j}(x), u_{j}(x_{i}) = δ_{ij}, 1 ≤ i, j ≤ M,  (2.1)

of the solution. Introducing vectors

\begin{aligned}R(x) & := & (ϕ(∥x − x_{1}∥),…, ϕ(∥x − x_{M}∥))^{T} \\ S(x) & := & (p_{1}(x),…, p_{Q}(x))^{T} , \end{aligned}

we assert

Theorem 2.2 The vector U(x) := (u_{1}(x),…, u_{M}(x))^{T} formed by the values of the Lagrange basis functions u_{1},…, u_{M} of (2.1) at x ∈ IR^{n} coincides with the solution U_{∗}(x) of the conditional minimization problem

min{U^{T} AU − 2U^{T} R(x) + ϕ(0) | U ∈ IR^{M}, P^{T} U = S(x)}.  (2.3)

Proof: If (2.3) is solved by Lagrange multiplier techniques, and if the solution is written in matrix form, there exists a vector V_{∗}(x) = (v_{1}(x),…, v_{Q}(x))^{T} of Lagrange multipliers such that

( \begin{matrix}A & P \\ P^{T} & 0 \end{matrix} ) ( \begin{matrix}U_{∗}(x) \\ V_{∗}(x) \end{matrix} ) = ( \begin{matrix}R(x) \\ S(x) \end{matrix} )  (2.4)

holds for the solution U_{∗}(x) of (2.3), which is unique since (2.4) has the same coefficient matrix as (1.3). For a single x = x_{j}, 1 ≤ j ≤ M, the right-hand side of (2.4) coincides with the j-th column of the coefficient matrix. Since the system is uniquely solvable, the vector U_{∗}(x_{j}) must coincide with U(x_{j}), the j-th unit vector. But as the components of U_{∗}(x) are of the form (1.1), we have U(x) = U_{∗}(x) for all x ∈ IR^{n}. □

Remarks: Cardinal interpolants like the u_{j} defined in (2.1) normally are constructed via the system (1.3) with data δ_{ij} for u_{j}. But (1.3) would yield the coefficients of u_{j}, not the vector U(x) of values u_{1}(x),…, u_{M}(x) at a fixed x ∈ IR^{n}. Furthermore, the value ϕ(0) in (2.3) will be required later, but could here be replaced by any other constant. While the approach of Madych and Nelson poses a variational problem in an infinite-dimensional Hilbert space, our approach uses the finite–dimensional conditional minimization problem (2.3).

We now do the same thing for derivatives. If ϕ is differentiable of order |μ| on (0, ∞) and of order 2|μ| around zero, and if we define ψ(x) = ϕ(||x||), the solution U_{∗}^{(μ)}(x) of the problem

min \{U^{T} AU − 2U^{T} R^{(μ)}(x) + ψ^{(2μ)}(0) | \begin{matrix}U ∈ IR^{M}, \\ P^{T} U = S^{(μ)}(x) \end{matrix} \}  (2.5)

uniquely exists and satisfies a nonsingular linear equation system like (2.4), which formally coincides with the system

( \begin{matrix}A & P \\ P^{T} & 0 \end{matrix} ) ( \begin{matrix}U^{(μ)}(x) \\ V^{(μ)}(x) \end{matrix} ) = ( \begin{matrix}R^{(μ)}(x) \\ S^{(μ)}(x) \end{matrix} )  (2.6)

obtained by differentiation of (2.4). Thus the derivatives U^{(μ)} of U(x) and V^{(μ)} of V (x) exist and satisfy (2.6), while U^{(μ)}(x) coincides with the solution U_{∗}^{(μ)}(x) of (2.5) by the same argument as above. It might be worth noting that (2.6) and (2.5) always imply polynomial reproduction up to order q in the sense

∑_{j=1}^{M} u_{j}^{(μ)}(x) · p(x_{j}) = p^{(μ)}(x)  (2.7)

for all p ∈ IP_{q}. The special choice of the additive constant in (2.5) will become apparent in (3.2); it will make the minimum value nonnegative, which is required for the following

Definition 2.8 Let F (r) = ϕ(\sqrt{r}) be conditionally strictly positive definite of order q, and assume ϕ ∈ C^{|μ|}(0, ∞), ϕ ∈ C^{2|μ|} around zero for μ ∈ IN_{≥0}^{n}. Then, for any distribution of centres x_{1},…, x_{M} satisfying (1.4), the nonnegative function function κ_{q}^{(μ)} defined by

(κ_{q}^{(μ)}(x))^{2} := min {U^{T} AU − 2U^{T} R^{(μ)}(x) + ψ^{(2μ)}(0) |U ∈ K_{q}^{(μ)}(x)}

with the set

K_{q}^{(μ)}(x) := {U = (u_{1},…, u_{M})^{T} ∈ IR^{M} |∑_{j=1}^{M} u_{j}p(x_{j}) = p^{(μ)}(x) for all p ∈ IP_{q} }  (2.9)

of admissible vectors is called the Kriging function at x.

As we shall see later, the Kriging function is the norm of the representer of the interpolation error functional on a reproducing-kernel Hilbert space. However, the next steps will directly prove the following facts:

a) The interpolation error can be bounded by the Kriging function,

b) The Kriging function can be expressed by a nonnegative integral.

Furthermore, section 5 will show that

c) there are admissible vectors from K_{q}^{(μ)}(x) for (2.5) which allow the integral (and thus the Kriging function) to be bounded by some power of h_{ρ}(x).

The first two steps will involve Fourier transforms which are not classically feasible unless we produce sufficiently “nice” admissible vectors U ∈ K_{q}^{(μ)}(x) for (2.5). From here on, we consider μ, q, and x as being fixed and do not always indicate dependence on these symbols in the notation.

3 Fourier transforms

If the hypotheses of Definition 2.8 hold and ψ(x) = ϕ(∥x∥) is an absolutely integrable function with a nicely behaving nonnegative Fourier transform ψ̂ satisfying

ψ(y) = ϕ(∥y∥) = \frac{1}{(2π)^{n}} ∫_{IR^{n}} e^{i⟨y,t⟩}ψ̂(t)dt, y ∈ IR^{n},  (3.1)

we can use identities like

∑_{j,k=1}^{M} w_{j}\overline{w_{k}}ϕ(∥x_{j} − x_{k}∥)

= \frac{1}{(2π)^{n}} ∫_{IR^{n}} ∑_{j,k=1}^{M} w_{j}\overline{w_{k}}e^{i⟨x_{j}−x_{k},t⟩}ψ̂(t)dt

= \frac{1}{(2π)^{n}} ∫_{IR^{n}} | ∑_{j=1}^{M} w_{j}e^{i⟨x_{j},t⟩}|^{2}ψ̂(t)dt

and

\begin{aligned}\frac{d^{μ}}{dx^{μ}} ϕ(∥x − x_{j}∥) & = & \frac{d^{μ}}{dx^{μ}} \frac{1}{(2π)^{n}} ∫_{IR^{n}} e^{i⟨x−x_{j},t⟩}ψ̂(t)dt \\ & = & \frac{1}{(2π)^{n}} ∫_{IR^{n}}(it)^{μ}e^{i⟨x−x_{j},t⟩}ψ̂(t)dt \end{aligned}

to get

\begin{aligned}& U^{T} AU − 2U^{T} R^{(μ)}(x) + ψ^{(2μ)}(0) \\ = & \frac{1}{(2π)^{n}} ∫_{IR^{n}} |∑_{j=1}^{M} u_{j}e^{i⟨x_{j},t⟩} − (it)^{μ}e^{i⟨x,t⟩}|^{2} ψ̂(t)dt \end{aligned}  (3.2)

for arbitrary U = (u_{1},…, u_{M}) ∈ IR^{M}, expressing the Kriging function via an integral. This explains the choice of ψ^{(2μ)}(0) as the additive constant in (2.5).

For most of the interesting radial basis functions ϕ, however, we have to use generalized Fourier transforms in (3.1) and (3.2). If (3.1) is interpreted as a generalized Fourier transform (e.g.: in the sense of [7]), the same interpretation applies to (3.2). Fortunately, we can circumvent these peculiarities, because the function

g_{U}(t) := ∑_{j=1}^{M} u_{j}e^{i⟨x_{j},t⟩} − (it)^{μ}e^{i⟨x,t⟩}  (3.3)

(for U, x_{j}, x, and μ fixed) has special properties which make the integral (3.2) well-defined in the classical sense for most of the interesting cases, provided that U = (u_{1},…, u_{M})^{T} lies in the admissible set K_{q}^{(μ)}(x) of (2.9). To prove this, we start with

Lemma 3.4 For every admissible vector U ∈ K_{q}^{(μ)}(x) the function g_{U}(t) has the property

|g_{U}(t)| ≤ \{ \begin{matrix}O(∥t∥^{q}) & for ∥t∥ → 0 \\ O(∥t∥^{|μ|}) & for ∥t∥ → ∞ \end{matrix} \} .  (3.5)

Proof. Let

e^{x} = p_{q}(x) + x^{q}r_{q}(x), p_{q} ∈ IP_{q}, |r_{q}(x)| ≤ e^{|x|}

be the Taylor expansion of e^{x} to order q, and let U ∈ K_{q}^{(μ)}(x) be admissible for (2.5). For q > 0 we use (2.9) to get

∑_{j=1}^{M} u_{j}p_{q}(i⟨x_{j} − x, t⟩) = (\frac{d^{μ}}{dy^{μ}} )_{y=x} p_{q}(i⟨y − x, t⟩)
= (it)^{μ}p_{q}^{(μ)}(0)
= \{ \begin{matrix}(it)^{μ} & 0 ≤ |μ| ≤ q − 1 \\ 0 & else \end{matrix} \} .

Then (3.3) yields

g_{U}(t)e^{−i⟨x,t⟩} = ∑_{j=1}^{M} u_{j}e^{i⟨x_{j}−x,t⟩} − (it)^{μ}

= ∑_{j=1}^{M} u_{j}p_{q}(i⟨x_{j} − x, t⟩) − (it)^{μ}

+∑_{j=1}^{M} u_{j}(i⟨x_{j} − x, t⟩)^{q}r_{q}(i⟨x_{j} − x, t⟩)

= ∑_{j=1}^{M} u_{j}(i⟨x_{j} − x, t⟩)^{q}r_{q}(i⟨x_{j} − x, t⟩)

for 0 ≤ |μ| ≤ q − 1. This proves (3.5), if (3.3) is used directly for ∥t∥ → ∞. The cases q = 0 and |μ| ≥ q now are easy.

With Lemma 3.4 and some additional assumptions on ψ̂ the integral in (3.2) can now be shown to exist classically:

Theorem 3.6 Let the generalized Fourier transform of ψ(x) = ϕ(∥x∥) exist and coincide with a continuous function ψ̂ on R^{n} \ {0} satisfying

0 < ψ̂(t) ≤ c \{ \begin{matrix}||t||^{−n−s_{0}} & for ∥t∥ → 0 \\ ||t||^{−n−s_{∞}} & for ∥t∥ → ∞ \end{matrix} \}  (3.7)

with constants c ∈ IR_{>0}, s_{0}, s_{∞} ∈ IR, where we additionally assume

2|μ| < s_{∞} and s_{0} < 2q.  (3.8)

Then for all U ∈ K_{q}^{(μ)}(x) we have (3.2) with a well-defined integral.

Proof: From Lemma 3.4, (3.8), and (3.7) we get the classical existence of the integral in (3.2). With (3.3) and the theorem on monotone convergence, (3.2) equals

\begin{matrix}lim_{m→∞} \frac{1}{(2π)^{n}} ∫_{R^{n}} |g_{U}(t)|^{2}e^{−||t||^{2}/m^{2}} ψ̂(t)dt \\ = lim_{m→∞} \frac{1}{(2π)^{n}} ∫_{IR^{n}} Ĝ_{m}(z)ψ(z)dz \end{matrix}  (3.9)

with the test functions

G_{m}(t) = |g_{U}(t)|^{2}e^{−||t||^{2}/m^{2}}

using the definition of the generalized Fourier transform for tempered distributions. Furthermore, the Fourier transform Ĝ_{m}(z) of G_{m} can be explicitly calculated up to a constant σ_{mn} as

m^{n}∑_{j,k} u_{j}u_{k}e^{−||z−(x_{j}−x_{k})||^{2}m^{2}/4}

− 2m^{n}(−1)^{n}∑_{k} u_{k}D^{μ}(e^{−||z−(x−x_{k})||^{2}m^{2}/4})
+ m^{n}D^{2μ}(e^{−||z||^{2}m^{2}/4}),

where D denotes differentiation with respect to z. Insertion into (3.9) yields (3.2), using the properties of the delta sequence σ_{mn}m^{n}e^{−||z||^{2}m^{2}/4}. □

The reproducing kernel Hilbert space approach of Madych and Nelson [9] [10] uses a different method of regularization of integrals of the form (3.2): they consider a space of test functions modulo L_{2}–orthogonality to polynomials and have to go all through a specific theory of distributions to make their variational problem well–defined. Note that in our approach there are no problems with the variational formulation; the specific integral representation needed elaboration. We now can write the Kriging function as an integral:

Theorem 3.10 Under the assumptions of Theorem 3.6 the Kriging function has the representation

(κ_{q}^{(μ)}(x))^{2} = min_{U∈K_{q}^{(μ)}(x)} \frac{1}{(2π)^{n}} ∫_{IR^{n}} |∑_{j=1}^{M} u_{j}e^{i⟨x_{j},t⟩} − (it)^{μ}e^{i⟨x,t⟩}|^{2} ψ̂(t)dt  (3.11)

with K_{μ,q}(x) defined as in (2.9), and the integral exists in the classical sense. □

4 Error bounds

If the data (x_{j}, f_{j}) stem from a smooth and absolutely integrable real-valued function f on IR^{n} with a nicely behaving Fourier transform f̂ satisfying

f(x) = \frac{1}{(2π)^{n}} ∫_{IR^{n}} e^{i⟨x,t⟩}f̂(t)dt, x ∈ IR^{n},

we can use the solution U^{(μ)}(x) = (u_{1}^{(μ)}(x),…, u_{M}^{(μ)}(x))^{T} of (2.6), which coincides with the derivatives of the Lagrange interpolation functions from (2.1), to find the error representation

|s^{(μ)}(x) − f^{(μ)}(x)|^{2} = |\frac{1}{(2π)^{n}} ∫_{IR^{n}} (∑_{j=1}^{M} u_{j}^{(μ)}(x)e^{i⟨x_{j},t⟩} − (it)^{μ}e^{i⟨x,t⟩}) f̂(t)dt|^{2}  (4.1)

of the interpolant s to f in the form (2.1), as far as it is feasible to take derivatives. Note that the bracketed function in (4.1) can be viewed as the Fourier transform of the representer of the error functional in some Hilbert space. It also is a special instance of a function of the form (3.3), but with optimal coefficients with respect to the minimization problem (2.5). The right-hand side of (4.1) is very similar to the Kriging function, and to relate the two we assume

c_{f}^{2} := \frac{1}{(2π)^{n}} ∫_{IR^{n}} |f̂(t)|^{2}(ψ̂(t))^{−1}dt < ∞  (4.2)

in addition to the hypotheses of Definition 2.8 and Theorem 3.6. Then, using the Cauchy– Schwarz inequality,

\begin{matrix}|s^{(μ)}(x) − f^{(μ)}(x)|^{2} ≤ \\ ≤ \frac{1}{(2π)^{n}} ∫_{IR^{n}} |∑_{j=1}^{M} u_{j}^{(μ)}(x)e^{i⟨x_{j},t⟩} − (it)^{μ}e^{i⟨x,t⟩}|^{2} ψ̂(t)dt · \frac{1}{(2π)^{n}} ∫_{IR^{n}} |f̂(t)|^{2}(ψ̂(t))^{−1}dt \\ = (κ_{q}^{(μ)}(x))^{2} · c_{f}^{2}, \end{matrix}  (4.3)

and the error of interpolation is pointwise bounded by the Kriging function.

This argument is feasible, whenever the Fourier transform f̂of f allows all intermediate integrals to exist classically. To cope with generalized Fourier transforms, we proceed as in the previous section:

Definition 4.4 A function f : IR^{n} → IR is dominated by a radial basis function ϕ satisfying (3.7) and (3.8) on R^{n} \ {0}, iff f has a generalized Fourier transform f̂ coinciding on R^{n} \ {0} with a continuous function satisfying (4.2) for ψ(x) = ϕ(∥x∥).

Remarks. The set F_{ϕ} of functions dominated by ϕ may be completed to form a Hilbert space with inner product

(f_{1}, f_{2})_{ϕ} = ∫_{IR^{n}} f̂_{1}(t)\overline{f̂_{2}(t)}(ψ̂(t))^{−1}dt

which was thoroughly studied by Madych and Nelson in [9], [10]. This is why the rather restrictive condition (4.2) occurs there, too, while (4.2) does not appear in [3] and [6]. We do not want to make direct use of Hilbert space properties here, but proceed directly to an estimate of κ_{q}^{(μ)}(x) in the next section. The application of the Cauchy-Schwarz inequality to (4.1) allows the error bound to be factored into a term c_{f} depending on f, but not on the data distribution, while the other factor consists of the Kriging function, which is independent of f but incorporates the sampling points. This is how the special form of (1.5) is obtained, and this is why we need (4.2). Of course, equation (4.3) shows that the Kriging function κ_{q}^{(μ)}(x) is nothing else than the norm of the representer of the error functional on Madych and Nelson’s Hilbert space. But knowledge of this fact does not improve the situation; a tight upper bound on κ_{q}^{(μ)}(x) is still to be constructed. This can be done by inserting special admissible vectors U ∈ K_{q}^{(μ)}(x) into (3.11). The next section will perform such a construction. We summarize the results of this section:

Theorem 4.5 If f is in the space F_{ϕ} of functions dominated by a radial basis function ϕ satisfying (3.8) and (3.7), then the interpolation error can be bounded by

|s^{(μ)}(x) − f^{(μ)}(x)| ≤ κ_{q}^{(μ)}(x) · c_{f},  (4.6)

where c_{f} is given by (4.2).

Proof: If g̃ is the bracketed function in (4.1), then g̃\sqrt{ψ̂} and f̂/\sqrt{ψ̂} are in L^{2}(IR^{n}). Thus the right-hand side of (4.1) is well-defined and (4.3) implies (4.6). □

5 Construction of admissible vectors for local error bounds

We start with a fundamental perturbation lemma which bridges the gap between regular and scattered data sets:

Lemma 5.1 For arbitrary constants ρ ∈ IR_{>0}, k ∈ IN there exist positive real constants h_{0}, c_{1}, c_{2} such that for any distribution of scattered centres x_{i} ∈ Ω, 1 ≤ i ≤ M and arbitrary points x ∈ Ω around which the local density of data points satisfies

h_{ρ}(x) := max_{y∈K_{ρ}(x)} min_{1≤i≤M} ∥y − x_{i}∥ ≤ h_{0}  (5.2)

and all μ ∈ IN_{≥0}^{n} with 0 ≤ |μ| ≤ k − 1 there is a vector Ũ := (ũ_{1}^{(μ)}(x),…, ũ_{M}^{(μ)}(x))^{T} in the admissible set K_{k}^{(μ)}(x) of (2.9) such that the inequalities

∥x_{j} − x∥ ≤ c_{1}h_{ρ}(x) for all j, 1 ≤ j ≤ M with ũ_{j}^{(μ)}(x) ≠ 0  (5.3)

∑_{j=1}^{M} |ũ_{j}^{(μ)}(x)| ≤ c_{2}h_{ρ}^{−|μ|}(x)  (5.4)

hold.

Proof. The matrix R_{k} = (α^{μ})_{0≤|μ|,|α|≤k−1} is a multidimensional form of the classical Vandermonde matrix; it is nonsingular because interpolation with polynomials of total order not exceeding k is possible on the set S_{k}^{n} := {α ∈ IN_{≥0}^{n} | 0 ≤ |α| ≤ k − 1}. Here and in the rest of this section the n-dimensional multi–indices α, β, μ, ν, and σ will always vary in S_{k}^{n}. When occurring as indices in matrices, we assume a fixed (e.g.: lexicographic) ordering, and the first index at the brackets in a matrix notation like B = (b_{μβ})_{μ,β} always indicates the row index. We use ρ, k, and n to define constants N, c_{1}, c_{2}, h_{0} by

\begin{aligned}(a) & N & := & max (1, 2(k − 1)(2k − 1)^{k}\binom{n+k−1}{n}n^{k/2}∥R_{k}^{−1}∥_{∞}) , \\ (b) & c_{1} & := & 1 + N(k − 1)\sqrt{n}, \\ (c) & c_{2} & := & 2(k − 1)!\binom{n+k−1}{n}^{2}∥R_{k}^{−1}∥_{∞}, \\ (d) & h_{0} & := & \frac{ρ}{N(k − 1)}, \end{aligned}  (5.5)

where ∥.∥_{∞} is the row-sum matrix norm. Now let the distribution of centres be dense enough to guarantee (5.2) for x ∈ IR^{n}. Then near all of the points x + Nhα with α ∈ S_{k}^{n}, which are in the ball K_{ρ}(x) := {y ∈ IR^{n} | ∥y − x∥ ≤ ρ} by (d) of (5.5), there must be points z_{α} ∈ X with

∥z_{α} − x − Nhα∥ ≤ h := h_{ρ}(x).  (5.6)

Now we get (5.3) for precisely these points z_{α}, if we use (b) of (5.5) for

\begin{aligned}∥z_{α} − x∥ & ≤ & h + ∥Nhα∥_{2} \\ & ≤ & h(1 + N(k − 1)\sqrt{n}) = c_{1}h. \end{aligned}  (5.7)

We now want to bound the difference of the two Vandermonde matrices V := ((α · N · h)^{μ})_{μ,α} and Z := ((z_{α} − x)^{μ})_{μ,α} with respect to regular and scattered data, respectively, using (5.6). By a straightforward expansion of the monomial t^{μ} we find

|t^{μ} − (t + y)^{μ}| ≤ |μ| ∥y∥_{∞} ( max_{0≤τ≤1}(∥t + τy∥_{∞}))^{|μ|−1}

and get the bound

|(z_{α} − x)^{μ} − (αNh)^{μ}| ≤ |μ|h((c_{1} + (k − 1)N)h)^{|μ|−1}  (5.8)

for the matrix elements of Z − V. The scaled regular Vandermonde matrix V = ((αNh)^{μ})_{μ,α} is expressible as

V = ((αNh)^{μ})_{μ,α} = (δ_{μα} · (Nh)^{|μ|})_{μ,α} · R_{k}  (5.9)

via the unscaled (integer) Vandermonde matrix R_{k}. We now use (5.8) and (5.9) for

\begin{aligned}∥V^{−1}(Z − V)∥_{∞} & = & ∥R_{k}^{−1} · (δ_{μα}(Nh)^{−|μ|})_{μ,α}(Z − V)∥_{∞} \\ & ≤ & ∥R_{k}^{−1}∥_{∞}max_{μ}(Nh)^{−|μ|}|μ|h^{|μ|}(c_{1} + (k − 1)N)^{|μ|−1}\binom{n + k − 1}{n} \\ & ≤ & ∥R_{k}^{−1}∥_{∞}\binom{n + k − 1}{n}\frac{1}{N} max_{μ} |μ| (\frac{c_{1}}{N} + k − 1)^{|μ|−1}. \end{aligned} (5.10)

With (a) and (b) of (5.5) we find

\frac{c_{1}}{N} + k − 1 = \frac{1}{N} + (k − 1)(1 + \sqrt{n})
≤ 1 + (k − 1)(1 + \sqrt{n})
≤ k + (k − 1)\sqrt{n}
≤ (2k − 1)\sqrt{n}

and continue (5.10) to get

∥V^{−1}(Z − V)∥_{∞} ≤ ∥R_{k}^{−1}∥_{∞} \binom{n + k − 1}{n}\frac{1}{N} · (k − 1)(2k − 1)^{k}n^{k/2} ≤ \frac{1}{2}.  (5.11)

By Neumann’s series, B := Z^{−1}V exists and its norm is bounded by 2. Some elementary matrix calculations give

(δ_{αμ}(Nh)^{|μ|})_{μ,α} = (ZBR_{k}^{−1})_{μ,α}

and we can define

ũ_{i}^{(μ)}(x) := \{ \begin{matrix}μ! (BR_{k}^{−1})_{α,μ} (Nh)^{−|μ|} & if x_{i} = z_{α} \\ 0 & otherwise \end{matrix} \} .

Note that the values ũ_{i}^{(μ)}(x) vanish except for those points x_{i} for which we have proven (5.7). This finishes (5.3). We now assert

∑_{j=1}^{M} ũ_{i}^{(μ)}(x)p(x_{j}) = p^{(μ)}(x)

for all p ∈ IP_{k}, which is equivalent to Ũ := (ũ_{1}^{(μ)}(x),…, ũ_{M}^{(μ)}(x))^{T} ∈ K_{k}^{(μ)}(x) of (2.9) and to

∑_{j=1}^{M} ũ_{i}^{(μ)}(x)p(x_{j} − x) = p^{(μ)}(0)  (5.12)

for all p ∈ IP_{k}. Representing an arbitrary p ∈ IP_{k} in the monomial basis as

p(y) = ∑_{|ν|<k} p_{ν}y^{ν}

we use the elements Z_{ν,α} = (z_{α} − x)^{ν} of the Vandermonde matrix Z to get

p(z_{α} − x) = ∑_{ν} p_{ν}(z_{α} − x)^{ν} = ∑_{ν} Z_{ν,α}p_{ν}

and prove (5.12) via

p^{(μ)}(0) = μ!p_{μ} = μ!(Nh)^{−|μ|} ∑_{|σ|<k} p_{σ}δ_{μ,σ}(Nh)^{|σ|}
= μ!(Nh)^{−|μ|}∑_{σ} p_{σ} (ZBR_{k}^{−1})_{σ,μ}
= μ!(Nh)^{−|μ|}∑_{α} (BR_{k}^{−1})_{α,μ}∑_{σ} Z_{σ,α}p_{σ}
= μ!(Nh)^{−|μ|}∑_{α} (BR_{k}^{−1})_{α,μ}p(z_{α} − x)

= ∑_{j=1}^{M} ũ_{j}^{(μ)}(x)p(x_{j} − x).

Finally, we establish (5.4) using (c) of (5.5) and

∑_{j=1}^{M} |ũ_{i}^{(μ)}(x)| ≤ μ!(Nh)^{−|μ|}∑_{α} | (BR_{k}^{−1})_{α,μ} |
≤ μ!(Nh)^{−|μ|}∥BR_{k}^{−1}∥_{1}
≤ (k − 1)!h^{−|μ|}\binom{n+k−1}{n}^{2}2∥R_{k}^{−1}∥_{∞},

where we had to go from the column-sum norm ∥.∥_{1} to the row-sum norm ∥.∥_{∞}. □

Now assume ρ ∈ IR_{>0} and an integer k ≥ q to be chosen for Lemma 5.1 such that (5.3) and (5.4) hold, whenever (5.2) is satisfied. Then the proof technique of Lemma 3.4 yields

\begin{aligned}|gŨ(t)e^{−i⟨x,t⟩}| & = & |∑_{j=1}^{M} ũ_{j}^{(μ)}e^{i⟨x_{j}−x,t⟩} − (it)^{μ}| \\ & ≤ & c_{2}h_{ρ}^{−|μ|}(x)c_{1}^{k}h_{ρ}^{k}(x)∥t∥^{k}e^{c_{1}h_{ρ}(x)∥t∥} \end{aligned}

for 0 ≤ |μ| < k, which specializes to

|gŨ(t)e^{−i⟨x,t⟩}| ≤ c_{2}c_{1}^{k}e^{c_{1}}∥t∥^{k}h_{ρ}^{k−|μ|}(x)

for h_{ρ}(x) · ∥t∥ ≤ 1. Direct application of (5.4) to gŨ implies

|gŨ(t)e^{−i⟨x,t⟩}| ≤ c_{2}h_{ρ}^{−|μ|}(x) + ∥t∥^{|μ|} ≤ (c_{2} + 1)∥t∥^{|μ|}.

for h_{ρ}(x) · ∥t∥ ≥ 1. Both bounds can be inserted into the Kriging function to get

\begin{aligned}(κ_{q}^{(μ)}(x))^{2} & ≤ & c_{2}^{2}c_{1}^{2k}e^{2c_{1}}h_{ρ}^{2k−2|μ|}(x)K_{1}(2k, ψ̂, h_{ρ}(x)) \\ & + & (1 + c_{2})^{2}K_{2}(2|μ|, ψ̂, h_{ρ}(x)) \end{aligned}

for 0 ≤ |μ| < min(k, s_{∞}/2) with the abbreviations

\begin{aligned}K_{1}(j, z, h) & := & \frac{1}{(2π)^{n}} ∫_{∥t∥≤1/h} ∥t∥^{j}z(t)dt, \\ K_{2}(j, z, h) & := & \frac{1}{(2π)^{n}} ∫_{∥t∥≥1/h} ∥t∥^{j}z(t)dt. \end{aligned}  (5.13)

To evaluate K_{1}(2k, ψ̂, h) for ψ̂ satisfying Theorem 3.6, we make use of (3.8) and s_{0} < 2q ≤ 2k to get

\begin{aligned}K_{1}(2k, ψ̂, h) & = & K_{1}(2k, ψ̂, h_{0}) + \frac{1}{(2π)^{n}} ∫_{1/h_{0}≤∥t∥≤1/h} ∥t∥^{2k}ψ̂(t)dt \\ & ≤ & O(1) + O(h^{s_{∞}−2k}) \end{aligned}

for h → 0. Likewise, we find

K_{2}(2|μ|, ψ̂, h) = O(h^{s_{∞}−2|μ|})

using 2|μ| < s_{∞} from (3.8). We summarize:

Theorem 5.14 Let ϕ satisfy the assumptions of Theorem 3.6, and let ρ ∈ IR_{>0} be given. Then there exist positive real constants h_{0} and C such that for any distribution of centres x_{i} ∈ IR^{n}, 1 ≤ i ≤ M, and any point x ∈ IR^{n} with (5.2) the Kriging function can be bounded by

κ_{q}^{(μ)}(x) ≤ Ch_{ρ}^{s_{∞}/2−|μ|}(x)  (5.15)

for 0 ≤ |μ| < s_{∞}/2.

Proof: Lemma 5.1 allows to pick any k ≥ s_{∞}/2 and to apply the above argument, where large values of k will lead to small values of h_{0} and large values of C. Then

\begin{aligned}(κ_{q}^{(μ)}(x))^{2} & ≤ & O(h_{ρ}^{2k−2|μ|}(x)) · (O(1) + O(h_{ρ}^{s_{∞}−2k}(x))) + O(h_{ρ}^{s_{∞}−2|μ|}(x)) \\ & ≤ & O(h_{ρ}^{2k−2|μ|}(x)) + O(h_{ρ}^{s_{∞}−2|μ|}(x)) \end{aligned}

yields (5.15). □

6 Special radial basis functions

All examples of this section are based on Theorems 4.5 and 5.14 when specialized to certain radial basis functions ϕ. We obtain specific error estimates of the form (1.5) valid for all functions dominated by ϕ in the sense of Definition 4.4. The attainable orders of the interpolation error will only depend on the constants s_{∞} for the radial basis functions in question.

In case of Gaussians

ϕ(r) = e^{−αr^{2}}, α ∈ IR_{>0}, r ∈ IR_{≥0}

the Fourier transform

ψ̂(y) = 2^{−n}(απ)^{−n/2}e^{−∥y∥^{2}/(4α)}

decays exponentially at infinity and we have arbitrarily large values of s_{∞}, giving

Theorem 6.1 For all functions f : IR^{n} → IR having generalized Fourier transforms f̂(t) satisfying

∫_{IR^{n}} ∥f̂(t)∥^{2}e^{∥t∥^{2}}dt < ∞

the error of interpolation by Gaussian radial basis functions has arbitrarily high local order in the sense of (1.5). □

In case of multiquadrics

ϕ(r) = (c^{2} + r^{2})^{s/2} s ∈ IR_{>−n}, s ∉ 2ZZ,

we have generalized Fourier transforms containing a Bessel function of the second kind and decaying exponentially at infinity like in the Gaussian case, but now with a singularity of type ∥y∥^{−n−s} at zero. Thus again s_{∞} is arbitrarily large, yielding

Theorem 6.2 For all functions f : IR^{n} → IR having generalized Fourier transforms f̂(t) satisfying

∫_{IR^{n}} |f̂(t)|^{2}∥t∥^{n+s}e^{∥t∥}dt < ∞

the error of interpolation by multiquadric radial basis functions

ϕ(r) = (c^{2} + r^{2})^{s/2} s ∈ IR_{>−n}, s ∉ 2ZZ

has unbounded local order in the sense of (1.5).

The above results for special radial basis functions were already obtained by Madych and Nelson [10] as a result of their reproducing kernel Hilbert space theory. Within the latter, the integral occurring in (3.2) plays an important part, as was pointed out by N. Dyn in her lucid survey [4] of Madych and Nelson’s work. There the integral is split at ∥t∥ = 1 and not at ∥t∥ = 1/h, which makes it impossible to handle the cases that follow.

The radial basis functions

\begin{aligned}ϕ(r) & = & r^{s}, & s ∈ IR_{>1}, s ∉ 2IN, \\ ϕ(r) & = & r^{s} log r, & s ∈ 2IN, \end{aligned}

have generalized Fourier transforms ψ̂(t) with behaviour like ∥t∥^{−n−s} on all of IR^{n}. This requires s_{∞} = s, q > s/2 and yields Theorem 6.3 For all functions f : IR^{n} → IR having generalized Fourier transforms satisfying

∫_{IR^{n}} |f̂(t)|^{2}∥t∥^{n+s}dt < ∞

the error of interpolation by radial basis functions

\begin{aligned}ϕ(r) & = & r^{s}, & s ∈ IR_{>1}, s ∉ 2IN, \\ ϕ(r) & = & r^{s} log r, & s ∈ 2IN, \end{aligned}

has local order s/2 in the sense of (1.5).

Remarks. The convergence rates obtained so far are quite different from those of Buhmann [1] for gridded data. Infinite orders, especially for the Gaussian kernel, are unexpected as are the rather low orders for the functions r^{s} and r^{s} log r. The gridded case is fundamentally different from the scattered data case in several aspects:

• Convergence orders, as given by our technique, are mainly dominated by the behaviour of the Fourier transform of ϕ near infinity. However, similar convergence orders are obtained, if multiquadric interpolation with a fixed constant c is used on a regular grid (see Buhmann and Dyn [2]).

• Our methods require the interpolated functions to have Fourier transforms which are “dominated” by ϕ̂ in the sense of (4.2). This requirement is not surprising, because a reasonable local convergence order cannot occur if the “high frequency part” of f exceeds that of ϕ. Interpolation on grids hides this phenomenon.

• The “dominance” requirement (4.2) imposes a rather hard restriction on f, which has no direct counterpart in the infinite grid case. Further research should try to weaken (4.2) at the expense of smaller (and possibly non–local) convergence orders.

References

[1] Buhmann, M.D., Multivariable Interpolation using Radial Basis Functions, Ph.D. dissertation, University of Cambridge, 1989.

[2] Buhmann, M.D., and N. Dyn, Error estimates for multiquadric interpolation, in Curves and Surfaces, P.-J. Laurent, A. LeMéhauté, L.L. Schumaker (eds.), Academic Press 1991

[3] Duchon,J., Sur L’erreur d’interpolation des fonctions de plusiers variables par les D^{m}− splines, RAIRO Analyse Numerique 12 (1978) 325-334.

[4] Dyn, N., Interpolation of scattered data by radial functions, in: Topics in Multivariate Approximation, C.K. Chui, L.L. Schumaker, and F.I. Utreras (eds.) Academic Press, 1987

[5] Dyn. N, Levin, D., and S. Rippa, Numerical procedures for surface fitting of scattered data by radial functions, SIAM J. Sci. Stat. Comput. 7, 639–659

[6] Jackson, I.R.H., An order of convergence for some radial basis functions, IMA J. of Numer. Anal. 1987 [7] Jones, D.S., The Theory of Generalised Functions, Cambridge University Press, Cambridge

1982

[8] Micchelli, C.A., Interpolation of scattered data: distance matrices and conditionally positive definite functions, Constructive Approximation 2, 11-22

[9] Madych, W.R., and Nelson, S.A., Multivariate Interpolation and conditionally positive definite functions, Approx. Theory and its Applications 4 (1988), 77–89

[10] Madych, W.R., and Nelson, S.A., Multivariate Interpolation and conditionally positive definite functions II, Math. of Computation 54 (1990), 211–230

[11] Madych, W.R., and Nelson, S.A., Error Bounds for Multiquadric Interpolation, in: Approximation Theory VI: Vol 2. C. K. Chui, L. L. Schumaker and J. D. Ward.(eds.) 413-416.

[12] Powell, M.J.D., Radial basis functions for multivariable interpolation: a review, in Numerical Analysis 1987, D.F. Griffiths and G.A. Watson (eds.), Longman Scientific & Technical (Harlow), 223–241

[13] Wu, Zong–min, Die Kriging–Methode zur Lösung mehrdimensionaler Interpolationsprobleme, Dissertation, Göttingen 1986

Author’s addresses:

Zong–min Wu, Dept. of Mathematics, Fudan University, Shanghai, 200433 People’s Republic of China

Robert Schaback, Institut für numerische und angewandte Mathematik, Universität Göttingen, Lotzestraße 16–18, D–3400 Göttingen, FRG