Comparison of Radial Basis Function Interpolants

Robert Schaback

Abstract

This paper compares radial basis function interpolants on different spaces. The spaces are generated by other radial basis functions, and comparison is done via an explicit representation of the norm of the error functional. The results pose some new questions for further research.

§1. Introduction

We consider interpolation of real-valued functions f defined on a set Ω ⊆ IR^{d}, d ≥ 1. These functions are evaluated on a set X := {x_{1}, …, x_{N_{X}} } of N_{X} ≥ 1 pairwise distinct points x_{1}, …, x_{N_{X}} in Ω. If N ≥ 2, d ≥ 2 and Ω ⊆ IR^{d} are given with Ω containing at least an interior point, it is well known that there is no N-dimensional space of continuous functions on Ω that contains a unique interpolant for every f and every set X = {x_{1}, …, x_{N_{X}} } ⊂ Ω ⊆ IR^{d} consisting of N = N_{X} data points.

Thus the family of interpolants must necessarily depend on X. This can easily be achieved by using translates Φ(x − x_{j}) of a single continuous realvalued function Φ defined on IR^{d}, and further simplification is obtained by letting Φ be radially symmetric, i.e.:

Φ(x) := ϕ(∥x∥_{2})  (1)

with a continuous real-valued function ϕ on IR_{≥0} and the L_{2} norm ∥ · ∥_{2}. Interpolants s_{f} to f can then be constructed via the representation

s_{f}(x) = ∑_{j=1}^{N_{X}}α_{j}Φ(x − x_{j}),  (2)

Multivariate Approximation and Wavelets K. Jetter and F.I. Utreras (eds.) Copyright, World Scientific, Singapore, 1992, 1– 12. ISBN x-yy-zzzzz-w: All rights of reproduction in any form reserved. where the coefficients α_{1}, …, α_{N_{X}} ∈ IR solve the linear system

f(x_{k}) = ∑_{j=1}^{N_{X}}α_{j}Φ(x_{k} − x_{j}), 1 ≤ k ≤ N_{X},

provided that the symmetric N_{X} × N_{X} matrix

A_{Φ,X} := ( \begin{matrix}Φ(x_{1} − x_{1}) & … & Φ(x_{1} − x_{N_{X}}) \\ ⋮ & ... & ⋮ \\ Φ(x_{N_{X}} − x_{1}) & … & Φ(x_{N_{X}} − x_{N_{X}}) \end{matrix} )

is nonsingular. This is the simplest form of radial basis function interpolation, but for a variety of choices of Φ it is necessary to add polynomials to the interpolant (2).

So let P_{q}^{d} denote the space of d-variate polynomials of order not exceeding q, and let the polynomials p_{1}, …, p_{Q} be a basis of P_{q}^{d} in IR^{d}. The Q additional degrees of freedom of the extended representation

s_{f}(x) = ∑_{j=1}^{N_{X}}α_{j}Φ(x − x_{j}) + ∑_{ℓ=1}^{Q}β_{ℓ}p_{ℓ}(x)  (3)

are compensated by the Q additional equations

∑_{j=1}^{N_{X}}α_{j}p_{ℓ}(x_{j}) = 0, 1 ≤ ℓ ≤ Q.  (4)

With the matrix

P_{X}^{T} := ( \begin{matrix}p_{1}(x_{1}) & … & p_{1}(x_{N_{X}}) \\ ⋮ & ... & ⋮ \\ p_{Q}(x_{1}) & … & p_{Q}(x_{N_{X}}) \end{matrix} ) .

we can write the interpolation conditions

f(x_{k}) = ∑_{j=1}^{N_{X}}α_{j}Φ(x_{k} − x_{j}) + ∑_{ℓ=1}^{Q}β_{ℓ}p_{ℓ}(x_{k}), 1 ≤ k ≤ N_{X}

together with (4) as a linear system

( \begin{matrix}A_{Φ,X} & P_{X} \\ P_{X}^{T} & 0 \end{matrix} ) ( \begin{matrix}α \\ β \end{matrix} ) = ( \begin{matrix}f_{X} \\ 0 \end{matrix} ) ,  (5)

where the data from f form a vector f_{X} := (f(x_{1}), …, f(x_{N_{X}}))^{T}. Solvability of this system depends on two conditions. First, the matrix A_{Φ,X} should be nonsingular on the vectors α satisfying (4). Second, polynomials in P_{q}^{d} should be uniquely determined by their values on X, i.e.:

If p ∈ P_{q}^{d} satisfies p(x_{i}) = 0 for all x_{i} ∈ X then p = 0.  (6)

The discussion of the first condition is simplified if nonsingularity is replaced by positive definiteness:

Definition 1. A function Φ : IR^{d} → IR with Φ(−x) = Φ(x) is conditionally positive definite of order q 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 (4) the quadratic form

∑_{j,k=1}^{N_{X}}α_{j}α_{k}Φ(x_{j} − x_{k})

attains nonnegative values and vanishes only if α = 0.

In a fundamental paper C.A. Micchelli [14] related the conditional positive definiteness of functions Φ of the form (1) to complete monotonicity of derivatives of ϕ, and this technique allows to prove conditional positive definiteness for a variety of radial basis functions. We list a few examples:

Multiquadrics ϕ(r) = (c^{2} + r^{2})^{β/2} for β ∈ IR_{>−d} \ 2ZZ and 2q > β [8], Thin-plate splines ϕ(r) = r^{β} for β ∈ IR_{>0} \ 2ZZ and 2q > β [3, 4, 5], Thin-plate splines ϕ(r) = (−1)^{β/2+1}r^{β} log r for β ∈ 2IN, 2q > β [3, 4, 5], Gaussians ϕ(r) = e^{−αr^{2}} for α > 0 and q ≥ 0.

Our main purpose here is to study the error f(x) − s_{f}(x) of different radial basis function interpolants on different spaces. Curiously enough, each conditionally positive definite function Φ does not only define an interpolation method, but also defines an inner-product space F_{Φ} of functions. We describe the construction of such a space in the next section and introduce the exponentially decaying positive definite radial basis functions that generate the Sobolew spaces W_{2}^{k}(IR^{d}). Then we represent the norm of the error functional of a general linear quasi-interpolation method on such a space by a numerically accessible power function. Finally, we evaluate the power functions that arise from interpolation with a radial basis function Φ_{I} on a space F_{Φ_{S}} defined by a different radial basis function Φ_{S}. Inspection of the results leads to a number of open questions for further research.

§2. Spaces Generated By Radial Basis Functions

We assume the radial basis function Φ to be conditionally positive definite of order q on IR^{d} in the sense of Definition 1, and we now construct an associated function space F_{Φ} using ideas from Madych and Nelson [11].

Let Ω ⊆ IR^{d} be given, and let V be the set of all pairs (α, X) with the following properties:

X = {x_{1}, …, x_{N_{X}} } ⊆ Ω ⊆ IR^{d}, |X| = N_{X}

(α, X) satisfies (4), and

X satisfies (6).

To avoid pathological cases we assume Ω to be large enough to contain at least one set X satisfying (6), to make sure that V is non-empty. Then we define the set

F_{Φ} := P_{q}^{d} + {f_{α,X} | (α, X) ∈ V }  (7)

with functions

f_{α,X}(x) := ∑_{j=1}^{N_{X}}α_{j}Φ(x − x_{j}), (α, X) ∈ V  (8)

defined on all of IR^{d}. Note that V depends on Ω, and that this fact induces a subtle dependence of F_{Φ} on Ω which still is a mystery except for the case Ω = IR^{d}, treated in detail by Madych and Nelson in [12]. Since F_{Φ} contains all finite linear combinations of translates of Φ(x) with coefficients satisfying (4) and X satisfying (6), the space F_{Φ} and its closures under different topologies are very natural candidates to study radial basis function approximation and interpolation. Here, we avoid to take closures, because they turn out to be irrelevant to our purposes. We rather investigate (7) as it is via the following

Lemma 1. The sum in (7) is direct, and F_{Φ} is a vector space over IR. Furthermore, each function f ∈ F_{Φ} has a unique representation

f = p + f_{α,X} with (α, X) ∈ V, p ∈ P_{q}^{d}.

Proof: Assume that the function (8) is a polynomial p ∈ P_{q}^{d}. Then (4) implies

\begin{aligned}0 & = ∑_{k=1}^{N_{X}}α_{k}p(x_{k}) \\ & = ∑_{j,k=1}^{N_{X}}α_{j}α_{k}Φ(x_{j} − x_{k}) \end{aligned}

and α must vanish because Φ is conditionally positive definite of order q.

Let two functions f_{α,X} and f_{β,Y} with (α, X), (β, Y) ∈ V be given. Without loss of generality we assume

W := X ∪ Y =: {w_{1}, …, w_{N_{W}} }
⊇ Z := X ∩ Y =: {w_{1}, …, w_{N_{Z}} }
= {x_{1}, …, x_{N_{Z}} } = {y_{1}, …, y_{N_{Z}} },

where all of the sets X, Y, Z, and W contain pairwise distinct points. Then we can represent the sum of f_{α,X} and f_{β,Y} as

∑_{j=1}^{N_{Z}}(α_{j} + β_{j})Φ(x − w_{j}) + ∑_{j=N_{Z}+1}^{N_{X}}α_{j}Φ(x − x_{j}) + ∑_{j=N_{Z}+1}^{N_{Y}}β_{j}Φ(x − y_{j}).

On all polynomials p ∈ P_{q}^{d} we get

\begin{aligned}∑_{j=1}^{N_{Z}}(α_{j} + β_{j})p(w_{j}) + ∑_{j=N_{Z}+1}^{N_{X}}α_{j}p(x_{j}) + ∑_{j=N_{Z}+1}^{N_{Y}}β_{j}p(y_{j}) & = \\ ∑_{j=1}^{N_{Z}}α_{j}p(x_{j}) + ∑_{j=N_{Z}+1}^{N_{X}}α_{j}p(x_{j}) + ∑_{j=1}^{N_{Z}}β_{j}p(y_{j}) + ∑_{j=N_{Z}+1}^{N_{Y}}β_{j}p(y_{j}) & = 0, \end{aligned}

which proves that the sum of f_{α,X} and f_{β,Y} is representable in the form f_{γ,W} with a suitable vector γ ∈ IR^{N_{W}} such that (γ, W) ∈ V. The last assertion follows from a similar decomposition argument applied to the difference of two representations of a function f ∈ F_{Φ}.

On the space F_{Φ} in (7) we now define the bilinear form

(p + f_{α,X}, r + f_{β,Y})_{Φ} := ∑_{j=1}^{N_{X}}∑_{k=1}^{N_{Y}}α_{j}β_{k}Φ(x_{j} − y_{k}),  (9)

where p and r are arbitrary polynomials from P_{q}^{d}. Lemma 1 makes sure that this definition is consistent. Of course one can complete the pre-Hilbert space F_{Φ}/P_{q}^{d} into a rather interesting space of (generalized) functions on IR^{d}, but we refer the reader to Madych and Nelson [11,12,13] for details which are not relevant here.

§3. Radial Basis Functions for Sobolew Seminorms

We have seen that any conditionally positive definite function defines an inner product on a space of functions on IR^{d}. Conversely, one can ask for the radial basis functions that possibly generate a given inner product. The most prominent example would be the inner products generating Sobolew spaces W_{2}^{k}(IR^{d}) of functions having generalized derivatives up to order k on IR^{d}.

Theorem 1. For k > d/2 the inner product of Sobolew space W_{2}^{k}(IR^{d}) is given by (9) and (1) with the positive definite function

ϕ(r) = \frac{2π^{k}}{Γ(k)}K_{k−d/2}(2πr) · r^{k−d/2}  (10)

defined via the Macdonalds (or spherical Bessel) function K_{ν}.

Proof: The inner product of W_{2}^{k}(IR^{d}) can be written in the form

(f, g)_{W_{2}^{k}(IR^{d})} = ∫_{IR_{d}} f̂(ω)\overline{ĝ}(ω)(1 + ∥ω∥_{2}^{2})^{k}dω

using Fourier transforms (see Yosida [23], p. 155, eq. (30)). Now if

Φ̂(·)(ω) = (1 + ∥ω∥_{2}^{2})^{−k},  (11)

and if the inverse Fourier transform can be taken, then

(f_{α,X}, f_{β,Y})_{Φ} = ∑_{ℓ=1}^{N_{X}}∑_{j=1}^{N_{Y}}α_{ℓ}β_{j}Φ(x_{ℓ} − y_{j}),

= (2π)^{−d}∫_{IR_{d}} Φ̂(·)(ω)∑_{ℓ=1}^{N_{X}}α_{ℓ}e^{ix_{ℓ}·ω}∑_{j=1}^{N_{Y}}β_{j}e^{−iy_{j}·ω}dω

= (2π)^{−d}∫_{IR_{d}} f̂_{α,X}(ω)\overline{f̂}_{β,Y} (ω)(1 + ∥ω∥_{2}^{2})^{k}dω coincides with the inner product in Sobolew space. For k > d/2 the equation (11) has the solution (10). To prove this, combine p. 21 of Oberhettinger [16] with p. 158 of Stein and Weiss [19]. The radial basis function in (10) decays exponentially at infinity, behaves like a polynomial of degree 2k − d near zero, and is conditionally positive definite of order zero. See e.g.: Abramowitz and Stegun [1] for properties of Bessel functions. For spaces of odd dimension d, the radial basis function (10) can be calculated recursively in terms of monomials and the exponential function. In even dimensions the recursion starts at Bessel functions of index 0 and 1, which are readily available in any software package.

At this point we do not comment on the possible closures of F_{Φ} for various choices of Ω. An additional density argument will be needed for the proof of \overline{F_{Φ}} = W_{2}^{k}(Ω), which may be difficult for peculiar cases, e.g., Ω being finite, compact, or of dimension less that d.

§4. Error of Quasi-Interpolants on Radial Function Spaces

Now for any (λ, Y) ∈ V we can define a linear functional

φ_{λ,Y} (f) := ∑_{j=1}^{N_{Y}}λ_{j}f(y_{j})

that vanishes on P_{q}^{d} by (4) and satisfies

\begin{aligned}φ_{λ,Y} (f_{α,X}) & = ∑_{j=1}^{N_{Y}}∑_{k=1}^{N_{X}}λ_{j}α_{k}Φ(x_{k} − y_{j}) \\ & = (f_{λ,Y} , f_{α,X})_{Φ} \\ & = φ_{α,X}(f_{λ,Y}) \end{aligned}  (12)

for all (α, X) ∈ V. Thus φ_{λ,Y} is a continuous linear functional on F_{Φ} with respect to the seminorm | · |_{Φ} induced by the bilinear form (·, ·)_{Φ}. Its corresponding norm is

∥φ_{λ,Y} ∥_{Φ} = |f_{λ,Y} |_{Φ},

because its representer via the bilinear form is f_{λ,Y}. The square of this norm can be explicitly evaluated as

\begin{aligned}∥φ_{λ,Y} ∥_{Φ}^{2} & = |f_{λ,Y} |_{Φ}^{2} \\ & = ∑_{j=1}^{N_{Y}}∑_{k=1}^{N_{Y}}λ_{j}λ_{k}Φ(y_{k} − y_{j}). \end{aligned}  (13)

Now we specialize φ_{λ,Y} to be the error functional of some quasi-interpolant

s_{f,u,X}(x) := ∑_{i=1}^{N_{X}}u_{i}(x)f(x_{i})  (14)

that depends only on the data of a function f on a set X = {x_{1}, …, x_{N_{X}} } and which is exact on P_{q}^{d}. Here, the point x is fixed in Ω \ X and the values u_{1}(x), …, u_{N_{X}} (x) are just real numbers. We set Y := X ∪ {x} and λ := (1, −u_{1}(x), …, −u_{N_{X}} (x)) ∈ IR^{N_{X}+1} to get (λ, Y) ∈ V from the exactness of the quasi-interpolant (14) on P_{q}^{d}. Then we use (13) to derive the error bound

\begin{matrix}(f(x) − s_{f,u,X}(x))^{2} ≤ & |f|_{Φ}^{2} (Φ(0) − 2∑_{i=1}^{N_{X}}u_{i}(x)Φ(x − x_{i}) \\ + ∑_{i=1}^{N_{X}}∑_{j=1}^{N_{X}}u_{i}(x)u_{j}(x)Φ(x_{i} − x_{j})) \end{matrix}

for all f ∈ F_{Φ}. We call

\begin{aligned}P_{X,u,Φ}(x) := & (Φ(0) − 2∑_{i=1}^{N_{X}}u_{i}(x)Φ(x − x_{i}) \\ & + ∑_{i,j=1}^{N_{X}}u_{i}(x)u_{j}(x)Φ(x_{i} − x_{j}))^{1/2} \end{aligned}  (15)

the power function of the quasi-interpolant (14), because it precisely describes the quality of the quasi-interpolant at x. This is reminiscent of the notion of a power function of a statistical decision function.

Theorem 2. Let the radial basis function Φ be conditionally positive definite of order q on IR^{d} in the sense of Definition 1. Then any quasi-interpolant (14) that is exact on P_{q}^{d} satisfies

sup_{f∈F_{Φ}\P_{q}^{d}} \frac{|f(x) − s_{f,u,X}(x)|}{|f|_{Φ}} = P_{X,u,Φ}(x).  (16)

Discussion:

1. The sup in (16) can be extended to the Hilbert space completion of F_{Φ} without change of P_{X,u,Φ}(x). This is why we do not care about completions in this paper.

2. Though F_{Φ} will depend on the domain Ω, the right-hand side of (16) is independent of Ω, which is a rather startling fact at first sight. It can be explained by the observation that the sup is attained for the special function

f_{x}(y) := Φ(x − y) − ∑_{j=1}^{N_{X}}u_{i}(x)Φ(x_{i} − y)

which is in F_{Φ} and its completion whenever Y := X ∪ {x} is contained in Ω.

3. The expression (15) for the power function P_{X,u,Φ}(x) can be numerically evaluated at any x where the quasi-interpolant is defined. This allows a convenient numerical comparison between different quasi-interpolants on the same space F_{Φ}. We shall do this in section 6.

§5. Optimal Interpolants on Radial Function Spaces

Following Wu and Schaback [22] we now optimize the function (15) with respect to the N_{X} real variables u_{1}(x), …, u_{N_{X}} (x) under the constraints of exactness on P_{q}^{d}, which read as

p_{ℓ}(x) = ∑_{j=1}^{N_{X}}u_{j}(x)p_{ℓ}(x_{j}), 1 ≤ ℓ ≤ Q.

This yields a finite-dimensional convex optimization problem with linear constraints, and the solution vector u_{1}^{∗}(x), …, u_{N_{X}}^{∗} (x) with Lagrange multipliers v_{1}^{∗}(x), …, v_{Q}^{∗}(x) is characterized by the necessary and sufficient optimality conditions

( \begin{matrix}A_{Φ,X} & P_{X} \\ P_{X}^{T} & 0 \end{matrix} ) ( \begin{matrix}u^{∗}(x) \\ v^{∗}(x) \end{matrix} ) = ( \begin{matrix}R(x) \\ S(x) \end{matrix} )  (17)

with vectors

R(x) := (Φ(x − x_{1}), …, Φ(x − x_{N_{X}}))^{T} , S(x) := (p_{1}(x), …, p_{Q}(x))^{T}

and matrices as in (5). Thus (17) is uniquely solvable for all x ∈ IR^{d}, and since for x = x_{j} the right-hand side coincides with the j-th column of the coefficient matrix, we get u_{i}^{∗}(x_{j}) = δ_{ij}. This proves

Theorem 3. Among all linear quasi-interpolants of the form (14) that are exact on P_{q}^{d} and have data points in X, the radial basis function interpolation with centers in X is pointwise optimal with respect to minimization of the function (16).

Classical spline theory would minimize the seminorm induced by (9) under all interpolants in the space F_{Φ}. This variational approach was used in the radial basis function context by Duchon [3,4,5] and Madych and Nelson [11,12,13]. It is

a) a minimization of a “smoothness functional”

b) on an infinite-dimensional space

c) under all interpolants in the space, while the approach of Wu and Schaback [22] is

a) a minimization of the norm of the pointwise error functional

b) on a finite-dimensional space

c) under all quasi-interpolants that need not necessarily lie in the space. Under general conditions these two approaches are not equivalent. To

give a short account of the relation between the two problems, we include a simple proof of

Theorem 4. The optimal radial basis function interpolant in the sense of Theorem 3 is also optimal with respect to the minimization of the seminorm induced by (9) under all interpolants in the space F_{Φ}.

Proof: The optimal interpolant s_{f,X} to a function f with data on X which minimizes the seminorm induced by (9) is characterized by the usual interpolation conditions and the property

(s_{f,X}, g)_{Φ} = (s_{f,X}, f_{λ,Y})_{Φ} = 0

for all g = p + f_{λ,Y} ∈ F_{Φ} with (λ, Y) ∈ V that vanish on X. It suffices to show that this orthogonality condition is satisfied by the interpolant s_{f,X}^{∗} = f_{α^{∗},X} + p^{∗} of the form (3) with (α^{∗}, X) ∈ V. But with (12) we easily get

(s_{f,X}^{∗}, g)_{Φ} = (s_{f,X}^{∗}, f_{λ,Y})_{Φ}
= (f_{α^{∗},X}, f_{λ,Y})_{Φ}
= φ_{α^{∗},X}(f_{λ,Y})
= φ_{α^{∗},X}(g)
= 0.

We note that Laurent [10] has a beautiful result that guarantees equivalence of the two variational approaches considered here, provided that they both have solutions and that they are based on inner-product spaces. Going further back, the minimization of the norm of the representer of a linear functional on a Hilbert space dates seems to have been started by Golomb and Weinberger [7] and was carried forth by a series of others, including de Boor and Lynch [2], Sard [17], Larkin [9], and Dyn [6]. Somewhat related to this approach is the theory of optimal recovery, as given by Micchelli and Rivlin in [15].

§6. Numerical Results

We now perform numerical comparisons between different radial basis function interpolations on different spaces. Each space will be defined via a radial basis function Φ_{S}, and each interpolation will be carried out with a radial basis function Φ_{I}, where the corresponding orders of conditional positive definiteness are q_{S} and q_{I}, respectively.

We do not compare the interpolants themselves, but evaluate the functions (15). This will require q_{I} ≥ q_{S}, because the exactness order of the interpolant must at least equal the order of positive definiteness of the radial basis function defining the space (see the hypothesis of Theorem 2). Otherwise the interpolation error at x simply is not a continuous linear functional on the space in question. Tables 1 and 2 thus will have no entries in case q_{I} < q_{S}.

If Φ_{S} = Φ_{I}, there will be no better interpolant on F_{Φ_{S}} because of Theorem 3. We underlined these cases in Tables 1 and 2.

The numerical results were obtained in space dimension d = 1 for simplicity. The sets X contained 2n points in [−1, +1] with spacing h := 2/(2n − 1), and we set x = 0 to evaluate (15). Table 2 contains the actual values of (15) for n = 50 being fixed, while Table 1 gives approximate error orders p_{IS} for interpolation with Φ_{I} in F_{Φ_{S}} in the sense that (15) behaves like O(h^{p_{IS}}) for h → 0.

Spaces ⇒ Mq i.Mq r^{3}  W_{2}^{2}  r^{1}  W_{2}^{1}
Interpolations ⇓ q_{S} = 1 q_{S} = 0 q_{S} = 2 q_{S} = 0 q_{S} = 1 q_{S} = 0
Mq, q_{I} = 1 \underline{∞} ∞ - 1.5 0.5 0.5
Mq, q_{I} = 2 ∞  ∞ 1.5 1.5 0.5 0.5
i.Mq, q_{I} = 0 -  \underline{∞}  -  1.5  -  0.5
i.Mq, q_{I} = 1 ∞  ∞  -  1.5 0.5 0.5
i.Mq, q_{I} = 2 ∞  ∞ 1.5 1.5 0.5 0.5
r^{3}, q_{I} = 2 4.0  4.0 \underline{1.5} 1.5 0.5 0.5
W_{2}^{2}, q_{I} = 2 4.0  4.0 1.5 1.5 0.5 0.5
W_{2}^{2}, q_{I} = 0 -  4.0  -  \underline{1.5}  -  0.5
r^{1}, q_{I} = 1 2  2  -  1.5 \underline{0.5} 0.5
W_{2}^{1}, q_{I} = 1 2  2  -  1.5 0.5 0.5
W_{2}^{1}, q_{I} = 0 -  2  -  1.5  -  \underline{0.5}

Table 1. Estimated convergence orders.

Mq = Multiquadrics with β = 1/2
i.Mq = inverse Multiquadrics, β = −1/2
W_{2}^{k} = Sobolew function
Spaces ⇒ Mq  i.Mq  r^{3}  W_{2}^{2}  r^{1}  W_{2}^{1}
Interpolations ⇓ q_{S} = 1 q_{S} = 0 q_{S} = 2 q_{S} = 0 q_{S} = 1 q_{S} = 0
Mq, q_{I} = 1 \underline{0.002178} 0.061979 - 0.001675 0.155065 0.206437
Mq, q_{I} = 2 0.002178 0.061979 0.003082 0.001675 0.155065 0.206437
i.Mq, q_{I} = 0  -  \underline{0.061883}  -  0.001665  -  0.205494
i.Mq, q_{I} = 1 0.002182 0.061883  -  0.001665 0.154357 0.205494
i.Mq, q_{I} = 2 0.002182 0.061833 0.003064 0.001665 0.154357 0.205494
r^{3}, q_{I} = 2 0.002672 0.070031 \underline{0.002986} 0.001623 0.148943 0.198287
W_{2}^{2}, q_{I} = 2 0.002672 0.070038 0.002986 0.001623 0.148941 0.198284
W_{2}^{2}, q_{I} = 0  -  0.070038  -  \underline{0.001623}  -  0.198284
r^{1}, q_{I} = 1 0.010857 0.182091  -  0.002215 \underline{0.142857} 0.190178
W_{2}^{1}, q_{I} = 1 0.010893 0.182285  -  0.002230 0.148941 0.190178
W_{2}^{1}, q_{I} = 0  -  0.182285  -  0.002241  -  \underline{0.190178}

Table 2. Power function value at zero for 50 data points.

There are some remarkable observations to be made from the tables:

1. Quasi-optimal convergence: Nonoptimal interpolation processes with radial basis functions Φ_{I} on spaces defined by radial basis functions Φ_{S} ≠ Φ_{I} seem to achieve the order of the optimal interpolation on F_{Φ_{S}} if the optimal orders satisfy p_{II} ≥ p_{SS}. That is, radial basis function interpolants that are optimal on small and very smooth spaces apparently are quasi-optimal on larger and less smooth spaces.

2. Superconvergence: On even smoother spaces than their native space, a radial basis function interpolant may show an even better behavior. That is, nonoptimal interpolation processes with radial basis functions Φ_{I} on spaces defined by radial basis functions Φ_{S} ≠ Φ_{I} achieve a higher order than their optimal order on their basic space, if p_{II} < p_{SS}.

3. There is hardly any difference induced by variations of the polynomial orders q_{I} and q_{S}. The only visible deviation is in the last two entries of the W_{2}^{2} column of Table 2.

4. The values of (15) are astonishingly close to each other and to the optimal value.

The first two observations can be commented within the context of classical one-dimensional natural splines, which also solve a variational problem and which are a special case of the theory of this paper. On the Sobolew space H_{2}^{k}[−1, +1] with seminorm |f|_{k}^{2} := ∥f^{(k)}∥_{2}^{2}, the natural polynomial splines of order 2k are optimal and have the optimal error order k − 1/2, which is a saturation order (see Schumaker [18]). Natural splines of higher order 2n > 2k attain this order on H_{2}^{k}[−1, +1], too, and are quasi-optimal in the sense of the first observation. On the space C^{2k}[a, b], however, natural splines of order 2k will have an error order 2k in the interior of the domain, and this phenomenon was called superconvergence in the second observation. Since 2k is a saturation order on C^{2k}[a, b], no improvement is possible in this special case. These facts were proven by Swartz and Varga in [20] for the univariate spline case, but there is no proof so far for multivariate splines or general radial basis functions. This raises the question of saturation orders and saturation spaces for general cases of radial basis function interpolation. The underlined optimal orders p_{II} for Φ_{S} = Φ_{I} were theoretically proven by Madych and Nelson [11,12,13] and Wu and Schaback [22].

Of course, the above experiments were only done in the interior of the domain and for uniform meshes. Performance may be quite different if these assumptions are not satisfied. The classical spline case teaches us that the behavior near the boundary will not show superconvergence, and there is some theoretical investigation under way that suggests quasi-optimality to be dependent on asymptotically quasi-uniform meshes.

§7. References

1. Abramowitz, M., and I. A. Stegun, Handbook of Mathematical Functions, Dover, New York 1970

2. de Boor, C., and R. E. Lynch, On splines and their minimum properties, J. Math. Mech15 (1966), 953–969

3. Duchon, J., Interpolation des fonctions de deux variables suivant le principe de la flexion des plaques minces, RAIRO Anal. Num. 10 (1976), 5–12

4. Duchon, J., Splines minimizing rotation-invariant seminorms in Sobolev spaces, in Constructive Theory of Functions in Several Variables, W.

Schempp and K. Zeller (eds.), Springer, Berlin-Heidelberg-New York, 1977, 85–100

5. Duchon, J., Sur L’erreur d’interpolation des fonctions de plusieurs variables par les D^{m}− splines, RAIRO Anal. Num. 12 (1978), 325-334.

6. Dyn, N., On best nonlinear approximation in sign-monotone norms and in norms induced by inner products, SIAM J. Numer. Anal.16 (1979), 612–622

7. Golomb, W., and H. F. Weinberger, Optimal approximation and error bounds, in Numerical Approximation, R. E. Langer (ed.), Univ. of Wisconsin Press, 1959, 117-190

8. Hardy, R. L., Multiquadric equations of topography and other irregular surfaces, J. Geophysical Res. 76 (1971), 1905–1915

9. Larkin, F. M., Optimal Approximation in Hilbert Spaces with Reproducing Kernel Functions, Math. Comp. 24 (1970), 911–921

10. Laurent, P. J., Quadratic convex analysis and splines, in Methods of Functional Analysis in Approximation Theory, C. A. Micchelli, D. V. Pai, B. V. Limaye (eds.) ISNM 76, Birkhäuser, Basel, 1986, 17–43

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

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

13. 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.), Acadamic Press, New York, 413–416

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

15. Micchelli, C. A., and T. J. Rivlin, A survey of optimal recovery, in Optimal Estimation in Approximation Theory, C. A. Micchelli and T. J. Rivlin (eds.), Plenum Press, New York, 1977

16. Oberhettinger, F., Fourier Transforms of Distributions and Their Inverses, Academic Press, New York–London, 1973

17. Sard, A., Optimal Approximation, J. Functional Anal. 1 (1967), 222–224 18. Schumaker, L. L., Spline Functions: Basic Theory, Wiley, New York, 1981.

19. Stein, E. M., and G. Weiss, Introduction to Fourier Analysis on Euclidean Spaces, Princeton University Press, Princeton, 1971

20. Swartz, B. K., and R. S. Varga, Error bounds for spline and L-spline interpolation, J. Approx. Th. 6 (1972), 6–49

21. Wu, Zong–min, Die Kriging–Methode zur Lösung mehrdimensionaler Interpolationsprobleme, Dissertation, Göttingen 1986

22. Wu, Zong-min, and R. Schaback, Local Error Estimates for Radial Basis Function Interpolation IMA J. of Numer. Anal., to appear

23. Yosida, Functional Analysis, Springer, Berlin–Heidelberg–New York 1968

Help in proofreading was provided by Armin Iske and Marko Weinrich. Special thanks go to Nira Dyn, Pierre–Jean Laurent, and Larry Schumaker for pointing out some references.