Error Estimates and Condition Numbers for Radial Basis Function Interpolation

Robert Schaback

Institut für Numerische und Angewandte Mathematik Universität Göttingen D – 37083 Göttingen

Germany

Revised for

Advances in Computational Mathematics

October 6, 1994

\underline{Abstract}: For interpolation of scattered multivariate data by radial basis functions, an “uncertainty relation” between the attainable error and the condition of the interpolation matrices is proven. It states that the error and the condition number cannot both be kept small. Bounds on the Lebesgue constants are obtained as a byproduct. A variation of the Narcowich–Ward theory of upper bounds on the norm of the inverse of the interpolation matrix is presented in order to handle the whole set of radial basis functions that are currently in use.

1 Introduction

Interpolation by “radial” basis functions requires a function Φ : IR^{d} → IR, a space IP_{m}^{d} of d–variate polynomials of degree less than m, and interpolates data values y_{1},…, y_{N} ∈ IR at data locations (“centers”) x_{1},…, x_{N} ∈ IR^{d} by solving the system

\begin{aligned}∑_{j=1}^{N} α_{j}Φ(x_{j} − x_{k}) & + & ∑_{ℓ=1}^{Q} β_{ℓ}p_{ℓ}(x_{k}) & = & y_{k}, & 1 ≤ k ≤ N \\ ∑_{j=1}^{N} α_{j}p_{i}(x_{j}) & + & 0 & = & 0, & 1 ≤ i ≤ Q \end{aligned}  (1.1)

for a basis p_{1},…, p_{Q} of IP_{m}^{d}, where Q = ( \begin{matrix}m − 1 + d \\ d \end{matrix} ). See Table 1 for the most commonly

used examples. In self–evident matrix formulation the system (1.1) reads as

( \begin{matrix}A & P^{T} \\ P & 0 \end{matrix} ) ( \begin{matrix}α \\ β \end{matrix} ) = ( \begin{matrix}y \\ 0 \end{matrix} ) ,  (1.2)

and solvability is usually guaranteed by the requirements rank (P) = Q ≤ N and

λ∥α∥_{2}^{2} ≤ α^{T} Aα  (1.3)

for all α ∈ IR^{N} with P α = 0, where λ is a positive constant. The latter condition is called “conditional positive definiteness of order m” if it holds for a specific pair (m, Φ) and for arbitrary sets X = {x_{1},…, x_{N}} ⊂ IR^{d}. The condition rank P = Q ≤ N can be called “IP_{m}^{d}– nondegeneracy of X”, because on such sets polynomials from IP_{m}^{d} are uniquely determined by their values. Details can be found in review articles by M.J.D. Powell [14], N. Dyn [7] and M. Buhmann [5].

If m = 0, then IP_{m}^{d} and P disappear. In this case ∥A^{−1}∥_{2} ≤ λ^{−1} holds in the spectral norm. More generally, as we shall see in the next section, the quantity λ^{−1} controls the sensitivity of the solution vector α with respect to variations in the data vector y. Thus one is interested in lower bounds on λ that are as tight as possible. Such bounds were obtained by Ball [2], Narcowich and Ward [11], [12], [13], Ball, Sivakumar, and Ward [3], Baxter [4] and Sun [17], while lower bounds were supplied by Schaback [16]. If the data are values of a function f, one usually considers a fixed function space F and evaluates the error of the interpolant

s_{f} = ∑_{j=1}^{N} α_{j}Φ(· − x_{j}) + ∑_{l=1}^{Q} β_{ℓ}p_{ℓ}  (1.4)

defined by a solution of (1.1) with y_{k} = f(x_{k}). If F is defined via Φ itself in a natural way, the space F carries a specific seminorm |.|F and the bound for the error f(x) − s_{f}(x) takes the form

|f(x) − s_{f}(x)| ≤ |f|F · P (x),

where the power function P (x) just is the norm of the error functional on F evaluated at x. Note that P depends on x, X, Φ, and F, but not on f. For sake of completeness, we note that F contains all functions f : IR^{d} → R with

f(x) = (2π)^{−d} ∫_{IR^{d}} f̂(ω)e^{iω^{T} x}dω, x ∈ IR^{d}

where f̂ is in L_{1}(IR^{d}) and

|f|F^{2} := ∫_{IR^{d}} \frac{|f̂(ω)|^{2}}{φ(ω)} dω < ∞.

Here, φ denotes the d–variate generalized Fourier transform of Φ. It is assumed to be a positive continuous function on IR^{d} \ {0} satisfying a variational equation

∫_{IR^{d}} ∫_{IR^{d}} v(x)\overline{v(y)}Φ(x − y)dxdy = (2π)^{−d} ∫_{IR^{d}} φ(ω)|v(w)|^{2}dω  (1.5)

for all test functions v in the Schwartz space of tempered functions that additionally satisfy

v̂^{(α)}(0) = 0 for all α ∈ IN{d \atop ≥0}, |α| < m.

Note that φ coincides in the sense of Jones [8] with the generalized Fourier transform of Φ in the context of tempered distributions. This approach goes back to W. Madych and S.A. Nelson [9], where (1.5) explicitly occurs. Table 2 contains the functions φ corresponding to the various cases of Φ. If φ decays at least algebraically at infinity, e.g.:

0 < φ(ω) < C(1 + ∥ω∥)^{−k},

then F contains Sobolew space

W k 2 (IR^{d}) = {f | ∫_{IR^{d}} |f̂(ω)|^{2}(1 + ∥ω∥_{2})^{k}dω < ∞} .

However, if φ decays exponentially, e.g.: for Φ(x) = e^{−α∥x∥^{2}}, then F consists of C^{∞} functions. Further details can be found in W. Madych and S. Nelson [9], [10], N. Dyn [6], and R. Schaback [15].

Numerical observations and theoretical results have revealed that the error and the sensitivity, described by P (x) and λ^{−1}, seem to be intimately related. In particular, there is no case known where the error and the sensitivity are both reasonably small. There is a dichotomy: Either one goes for a small error and gets a bad sensitivity, or one wants a stable algorithm and has to take a comparably large error. This effect is reminiscent the Uncertainty Principle in quantum mechanics, and here it will take the very simple form

P^{2}(x) · λ^{−1}(x) ≥ 1  (1.6)

where λ(x) is defined via (1.3) for the matrix A_{x} that arises after adding an additional row and column for A of (1.2) for the location x = x_{0}. We shall prove (1.6) in section 2 and draw several conclusions.

To explain the latter, we have to introduce some additional notation. First, the known upper bounds for P (x) take the form

P^{2}(x) ≤ F (h(x))  (1.7)

where F : IR_{>0} → IR_{>0} is decreasing and

h(x) := max_{∥y−x∥_{2}≤ρ} min_{1≤j≤N} ∥y − x_{j}∥_{2}  (1.8)

measures the density of centers x_{j} from X around x. Note that h(x) depends on X and ρ > 0 (which is kept fixed), but not on Φ. Equation (1.7) holds uniformly for all points x and sets X such that

h(x) ≤ h_{0},

the positive constant h_{0} being dependent on d, m, ρ, and Φ. It requires no additional hypotheses on Φ exceeding those we made so far. The function F is dependent on the decay of the generalized Fourier transform φ of Φ at infinity. This is equivalent to being dependent on the smoothness of Φ. See the references in Table 1 for details concerning the construction of F.

Second, the known lower bounds for λ(x) take the analogous form

λ(x) ≥ G(q(x)),

where q(x) is the separation distance within the set X ∪ {x}, i.e.:

q(x) = \frac{1}{2} min_{0≤j<k≤N} ∥x_{j} − x_{k}∥  (1.9)

where we set x_{0} := x again. Table 1 contains references to some special instances; section 3 will provide a general theory based on previous work by Narcowich and Ward. No additional hypotheses on Φ are required. As h(x), the quantity q(x) depends on X and x, but there is a major difference between the two: if h(x) is small then there is a ball around x packed with centers from X of mutual distance at most h(x), while a small value of q(x) may possibly be attained for ∥x_{j} − x_{k}∥/2 with x_{j} and x_{k} far away from x. This difference is quite natural, because one cannot expect P (x) to be boundable from above in terms of q(x), nor can λ(x) in general be bounded from below by h(x). See Definition 2.2 below for a situation where such a bound exists.

Now our assertion (1.6) leads to a two–sided inclusion

G(q(x)) ≤ λ(x) ≤ P^{2}(x) ≤ F (h(x))

that serves to provide new upper bounds on λ(x) and new lower bounds on P (x). In addition, we now have an easy possibility to check how tight the bounds from the literature are in case of centers on the grid hZZ^{d}. Then q(x) = \frac{h}{2} + \sqrt{d}h is possible, and both λ(x) and P^{2}(x) should fit between

G(h) ≤ F (h\sqrt{d}),

and we shall see that F and G indeed differ only by a constant factor in case of thin–plate splines. In other cases, there are additional powers of δ that make a gap between

G(δ) = c · δ^{k} · F (δ) ≤ F (δ · \sqrt{d}),  δ → 0,  k > 0

which leaves room for further research.

For the reader’s convenience, we list the known examples for functions F and G in Table 1.

Φ(x) = ϕ(r), r = ∥x∥_{2}  F (h)  G(h)
r^{β}, β ∈ IR_{>0} \ 2IN  h^{β}  h^{β}
[2]: d = β = 1
thin–plate splines [18] [3], pg. 419: β ∈ (0, 2)
[13], § VI: β = m − d/2, d odd
(−1)^{1+β/2}r^{β} log r, β ∈ 2IN  h^{β}  h^{β}
thin–plate splines  [18]  [13], § VI: β = m − d/2, d even
(γ^{2} + r^{2})^{β/2}, β ∈ IR \ 2IN_{≥0}  e^{− \frac{δ}{h}}  h^{2}e^{− \frac{6}{h^{2}}} [11], pg. 90: γ = 1 = β, d = 2,
Multiquadrics  δ > 0. he^{− \frac{2d}{h}} [3], pgs. 422–423: γ = 1 = β
[10]  h^{β} exp(−12.76γd/h)
e^{−βr^{2}}, β > 0  e^{− \frac{δ}{h^{2}}}  h^{−d}e^{− \frac{γ}{h^{2}}} [13], pg. 90: β = 1
Gaussians  δ > 0 [10] h^{−d} exp(−40.71d^{2}/(βh^{2}))
\frac{2π^{d/2}}{Γ(k)} K_{k−d/2}(r)(r/2)^{k−d/2}  h^{2k−d} h^{2k−d}
2k > d,  as in [18]

Sobolev splines

Table 1: All entries are modulo factors that are independent of r and h, but possibly dependent on parameters of Φ. Unreferenced cases for G are treated in section 3, where we included the constants.

Another useful result of our analysis will be a bound on the Lagrange functions u_{1}(x),…, u_{N}(x) corresponding to interpolation by (1.1) on X = {x_{1},…, x_{N}}. In the above context, we get

1 + ∑_{j=1}^{N} u_{j}^{2}(x) ≤ \frac{P^{2}(x)}{λ(x)} ≤ \frac{F (h(x))}{G(q(x))} ,

and this serves to prove that the Lagrange functions cannot grow too badly in regions where there are sufficiently many regularly distributed centers.

2 Basic results

First we assert that λ^{−1} of (1.3) controls the sensitivity of the solution vector α ∈ IR^{N} of (1.2) with respect to perturbations of the data vector y ∈ IR^{N}. In fact,

α^{T} Aα = α^{T} y

follows from (1.2) and implies

∥α∥_{2} ≤ λ^{−1}∥y∥_{2},

if (1.3) holds. Similarly, an upper bound

α^{T} Aα ≤ Λ∥α∥_{2}^{2},  Λ > 0  (2.1)

of the quadratic form induced by A yields

Λ^{−1}∥y∥_{2} ≤ ∥α∥_{2}

and in case of a perturbation y + Δy of y that leads to a perturbation α + Δα of α, we get the condition–type estimate

\frac{∥Δα∥_{2}}{∥α∥_{2}} ≤ \frac{Λ}{λ} \frac{∥Δy∥_{2}}{∥y∥_{2}} .

Since the determination of the other solution vector β ∈ IR^{Q} of (1.2) can be viewed as a problem of polynomial interpolation, the main part of the sensitivity analysis of “radial” basis function interpolation problems consists in finding good bounds for λ and Λ.

We now turn to the function P (x). Due to the special choice of the space F and the seminorm |.|F, the function P^{2}(x) can be explicitly written as

\begin{aligned}P^{2}(x) = Φ(0) & − & 2∑_{j=1}^{N} u_{j}(x)Φ(x_{j} − x) \\ & + & ∑_{j=1}^{N} ∑_{k=1}^{N} u_{j}(x)u_{k}(x)Φ(x_{j} − x_{k}), \end{aligned}

where u_{1}(x),…, u_{N}(x) are the Lagrange basis functions for interpolation, i.e. (1.6) equals

s_{f} = ∑_{j=1}^{N} f(x_{j})u_{j}

(see e.g. Schaback [15]). We now set x_{0} := x and add a first row and column to the N × N matrix A in (1.2) to get a (N + 1) × (N + 1) matrix A_{x}. With

u_{x} := (1, −u_{1}(x),…, −u_{N}(x))^{T} ∈ IR^{N+1}

we then have

P^{2}(x) = u_{x}^{T} A_{x}u_{x}

and simply use (1.3) for A_{x} in the form

λ(x)∥γ∥_{2}^{2} ≤ γ^{T} A_{x}γ for all γ ∈ IR^{N+1}

to get

λ(x) (1 + ∑_{j=1}^{N} u_{j}^{2}(x)) ≤ P^{2}(x).  (2.2)

If we take (2.1) for A replaced by A_{x}, then

λ(x) ≤ \frac{P^{2}(x)}{1 + ∑_{j=1}^{N} u_{j}^{2}(x)} ≤ Λ(x).

We now can read off a series of statements:

Theorem 2.1 For any IP_{m}^{d}–nondegenerate set X = {x_{1},…, x_{N}} and all x ∈ IR^{d} \ X we have

1 ≤ 1 + ∑_{j=1}^{N} u_{j}^{2}(x) ≤ P^{2}(x)λ^{−1}(x).

If λ corresponds to A in (1.3), then

λ ≤ min {P_{j}^{2}(x_{j}) | X \ {x_{j}} is IP_{m}^{d}–nondegenerate, 1 ≤ j ≤ N}

where P_{j} is the power function for X \ {x_{j}}. □

In all practically relevant cases we have bounds of the form

\begin{aligned}P^{2}(x) & ≤ & F (h(x)) \\ λ(x) & ≥ & G(q(x)) \end{aligned}  (2.3)

with continuous and decreasing functions F and G for small arguments. The relation between q(x) and h(x), as defined in (1.9) and (1.8), is

2q(x) ≤ min_{1≤j≤N} ∥x − x_{j}∥ ≤ h(x),  (2.4)

but in order to relate the bounds in (2.3) we have to look at conditions under which h(x) can be bounded from above in terms of q(x).

Definition 2.2 We call a set of centers X = {x_{1},…, x_{N}} quasi–uniform for a compact set Ω ⊂ IR^{d}, if there are δ > 0 and h > 0 such that

a) h(x) ≤ h for all x ∈ Ω

b) 2q = min_{1≤j<k≤N} ∥x_{j} − x_{k}∥_{2} ≥ h · δ,

c) each set X \ {x_{j}}, 1 ≤ j ≤ N, is IP_{m}^{d}–nondegenerate.

Then we call h the density and δ the uniformity of X with respect to Ω.

Note that necessarily 0 < δ ≤ 1 due to (2.4), and that there are plenty of such sets if h is small enough.

Theorem 2.3 For quasi–uniform sets of centers with density h and quasi–uniformity δ we have two–sided bounds

G(\frac{1}{2} hε) ≤ λ(x) ≤ P^{2}(x) ≤ F (h)  (2.5)

for all x ∈ Ω with q(x) ≥ \frac{1}{2} h · ε. Furthermore, if (2.3) holds for arguments ≤ 2h in F and G, and if 2h ≤ ρ for ρ in (1.8), then

G(\frac{1}{2} hδ) ≤ λ ≤ F (2h).  (2.6)

Finally, F and G always satisfy

G(h) ≤ F (h\sqrt{d})

for sufficiently small arguments.

Proof: The bounds in (2.5) and the left part of (2.6) readily follow from the definition of quasi–uniformity, while for the right–hand part of (2.6) we have to prove

min_{1≤j≤N} h_{j}(x_{j}) ≤ 2h  (2.7)

if h_{j} is defined as

h_{j}(x) := sup_{∥y−x∥≤ρ} min_{k≠j^{1≤k≤N}} ∥y − x_{k}∥

like (1.8) after deletion of x_{j} from X = {x_{1},…, x_{N}}. For this, first take y ∈ IR^{d} with h < ∥y − x_{j}∥ ≤ ρ. Then there is some x_{k} ∈ X with k ≠ j and ∥y − x_{k}∥ ≤ h ≤ 2h, proving (2.7) in this case. If, however, y ∈ IR^{d} satisfies ∥y − x_{j}∥ ≤ h ≤ ρ, then we can find a z_{ε} ∈ IR^{d} with ∥z_{ε} − x_{j}∥ = h + ε for an arbitrarily small ε ∈ (0, h] such that y lies on the line between x_{j} and z_{ε}. Then there is some k ∈ {1,…, N}, k ≠ j, such that ∥z_{ε} − x_{k}∥ ≤ h, and consequently

∥y − x_{k}∥ ≤ ∥y − z_{ε}∥ + ∥z_{ε} − x_{k}∥ ≤ h + ε + h,

proving (2.7) for this case, too. The last assertion follows from the fact that for a uniform grid hZZ^{d} one has points x with

h(x) = h\frac{\sqrt{d}}{2} = \sqrt{d}q(x),

and thus

G(h/2) = G(q(x)) ≤ λ(x) ≤ P^{2}(x) ≤ F (h(x)) = F (h · \frac{\sqrt{d}}{2} ) .

□

The applications of these results can be read off Table 1. For thin–plate splines, the bounds in terms of F and G are of the same order, and the L_{2} Lebesgue function 1+∑_{j=1}^{N} u_{j}^{2}(x) will be bounded independent of N for x inside sets of sufficiently dense centers. Upper bounds of P^{2}(x) and lower bounds of λ(x) are best possible in terms of the order β. The results of the literature do not cover the full range of β, and this is why we add a section on lower bounds for λ(x).

The other cases show certain discrepancies between F (δ · \sqrt{d}) and G(δ) by factors δ^{k} that may or may not be consequences of insufficient proof techniques. But since in these cases the bounds decay exponentially anyway, one should rather look at the constants in the exponential in order to sharpen the bounds. As long as these constants are not known, one cannot say whether the L_{2} Lebesgue functions are exponentially increasing, decreasing or constant.

3 Lower bounds for λ

In this section we generalize the technique of Narcowich and Ward [11], [12], to provide Table 1 with a full set of examples for the G function. The main difference will be that we introduce Fourier transforms right from the start, which makes it much easier to treat large values of m, the order of conditional positive definiteness of Φ.

The starting point is that any function Φ of Table 1 satisfies a relation

∑_{j=1}^{N} ∑_{k=1}^{N} α_{j}α_{k}Φ(x_{j} − x_{k}) = (2π)^{−d} ∫_{IR^{d}} φ(ω) |∑_{j=1}^{N} α_{j}e^{ix_{j}^{T} ω}|^{2} dω  (3.1)

for all IP_{m}^{d}–nondegenerate sets X = {x_{1},…, x_{N}} and all vectors α ∈ IR^{N} satisfying the second set of equations in (1.1). Note that (3.1) can be derived from (1.5), as was shown by Madych and Nelson [9].

The left–hand side of (3.1) is the quantity α^{T} Aα that we want to bound from below, and we can do this by any minorant ψ on IR^{d} \ {0} of φ that satisfies

φ(ω) ≥ ψ(ω) ≥ 0 on IR^{d} \ {0}  (3.2)

and that itself leads to a similar quadratic form

∑_{j=1}^{N} ∑_{k=1}^{N} α_{j}α_{k}Ψ(x_{j} − x_{k}) = (2π)^{−d} ∫_{IR^{d}} ψ(ω) |∑_{j=1}^{N} α_{j}e^{ix_{j}^{T} ω}|^{2} dω  (3.3)

for another “radial” basis function Ψ and a weaker (or none) constraint on α ∈ IR^{N}. Furthermore, there should be an easy lower bound

α^{T} Bα ≥ λ∥α∥_{2}^{2}

for the left–hand side α^{T} Bα of (3.3). Then clearly for all α ∈ IR^{N} that are admissible,

α^{T} Aα ≥ α^{T} Bα ≥ λ∥α∥_{2}^{2},  (3.4)

as required. The basic trick of Narcowich and Ward now is to make B diagonally dominant, while ψ is obtained by chopping off φ appropriately. Before we proceed any further, here is the main result:

Theorem 3.1 With the function

φ_{0}(r) := inf_{∥ω∥_{2}≤2r} φ(ω)  (3.5)

we have

λ ≥ \frac{1}{2} \frac{φ_{0}(M)}{Γ (\frac{d}{2} + 1)} (\frac{M}{4\sqrt{π}} )^{d}  (3.6)

for any M > 0 satisfying

M ≥ \frac{12}{q} (\frac{πΓ^{2} (\frac{d}{2} + 1)}{9} )\frac{1}{d + 1}  (3.7)

or, a fortiori,

M ≥ \frac{6.38 d}{q} .  (3.8)

Proof: We start with any M > 0 and the characteristic function

χ_{M}(x) = \{ \begin{matrix}1 & ∥x∥_{2} ≤ M \\ 0 & else \end{matrix} \}

of the L_{2} ball with radius M. Then we define

ψ(ω) := ψ_{M}(ω) := \frac{φ_{0}(M)Γ (\frac{d}{2} + 1)}{2^{d}M^{d} π^{d/2}} (χ_{M} ∗ χ_{M})(ω)

and use the calculations of [13] to get (3.2) via

\begin{aligned}supp (ψ_{M}) & = & \{x ∈ IR^{d} : ∥x∥_{2} ≤ 2M\} =: C_{2M}(0) \\ ∥χ_{M} ∗ χ_{M}∥_{∞} & ≤ & vol(C_{M}(0)) = M^{d} \frac{2^{d}π^{d/2}}{Γ (\frac{d}{2} + 1)} . \end{aligned}

The radial basis function Ψ_{M} corresponding to ψ_{M} is obtained via the inverse Fourier transform

χ̌_{M}(x) = (\frac{M}{2π∥x∥_{2}} )^{d/2} J_{d/2}(M · ∥x∥_{2})

and the Convolution Theorem as

\begin{aligned}Ψ_{M}(x) & = & \frac{φ_{0}(M)Γ (\frac{d}{2} + 1)}{2^{d}M^{d} π^{d/2}} (χ_{M} ∗ χ_{M})^{∨}(x) \\ & = & \frac{φ_{0}(M)Γ (\frac{d}{2} + 1)}{2^{d}π^{d/2}} ∥x∥_{2}^{−d}J_{d/2}^{2}(M · ∥x∥_{2}) \end{aligned}

with J_{ν} denoting the Bessel function, satisfying

\begin{aligned}J_{d/2}^{2}(z) & ≤ & \frac{2^{d+2}}{πz} , z > 0 \\ lim_{z→0} z^{−d}J_{d/2}^{2}(z) & = & \frac{1}{2^{d}Γ^{2} (\frac{d}{2} + 1)} \end{aligned}  (3.9)

as was proven in [13]. The second formula yields

Ψ_{M}(0) = \frac{φ_{0}(M)}{Γ (\frac{d}{2} + 1)} (\frac{M}{4\sqrt{π}} )^{d}

and we assert diagonal dominance of the quadratic form in (3.3) by a suitable choice of M. We have

α^{T} Bα ≥ ∥α∥_{2}^{2} (Ψ_{M}(0) − max_{1≤j≤N} ∑_{k≠j^{k=1}}^{N} Ψ_{M}(x_{j} − x_{k}))

by Gerschgorin’s theorem, and the final bound will be of the form

λ ≥ \frac{1}{2} Ψ_{M}(0) = \frac{φ_{0}(M)}{2Γ (\frac{d}{2} + 1)} (\frac{M}{4\sqrt{π}} )^{d} ,  (3.10)

because we shall choose M such that

max_{1≤j≤N} ∑_{k≠j^{k=1}} Ψ_{M}(x_{j} − x_{k}) ≤ \frac{1}{2} Ψ_{M}(0).  (3.11)

This is done by a tricky summation argument of Narcowich and Ward [13] that proves (3.11) for M satisfying (3.7). It remains to show that (3.8) implies (3.7). We use a variation of Stirling’s formula in the form

Γ(1 + x) ≤ \sqrt{2πx}x^{x}e^{−x}e^{1/12x}, x > 0

to get

\frac{π}{9} Γ^{2} (\frac{d}{2} + 1) ≤ \frac{π^{2}}{9} d^{d+1}(2e)^{−d}e^{1/3d},

(\frac{π}{9} Γ^{2} (\frac{d}{2} + 1))\frac{1}{d + 1} ≤ d · (\frac{π^{2}}{9} )\frac{1}{d + 1} · (2e)^{− \frac{d}{d + 1}} e\frac{1}{3d(d + 1)}
≤ d\frac{π}{3\sqrt{2e}} · e^{1/6} ≤ d · 0.531

such that

M ≥ \frac{6.38}{q} d

is satisfactory for all cases. □

Note that this technique completely ignores the additional conditions on α that might lead to a larger lower bound. The advantage is that the result is fairly general and can be applied in all of the cases. We incorporated the results into Table 1 to supply a number of missing cases. This was done by application of Theorem 3.1 to the entries in Table 2. To keep the formulae short, we used (3.8) instead of (3.7), which would yield sharper, but much more complicated bounds. To treat multiquadrics as a specific example, we have to evaluate φ_{0}(M) via

φ_{1}(R) := inf_{0<r≤R} K_{ν}(r) · r^{−ν}, R > 0, ν = \frac{d + β}{2} .

We use equation 9.6.23 of Abramowitz and Stegun [1] to get

\frac{Γ (ν + \frac{1}{2} ) K_{ν}(r)}{\sqrt{π}(r/2)^{ν}} = ∫_{1}^{∞} e^{−rt}(t^{2} − 1)^{ν−1/2}dt

≥ ∫_{1}^{∞} e^{−rt}(t − 1)^{2ν−1}dt

= e^{−r} ∫_{0}^{∞} e^{−rs}s^{2ν−1}ds
= Γ(2ν)e^{−r} · r^{−2ν}

and

φ_{1}(R) ≥ \frac{Γ(2ν)}{Γ (ν + \frac{1}{2} )} \frac{\sqrt{π}}{2^{ν}} e^{−R} · R^{−2ν}.

By the doubling formula for the Γ function this can be simplified to

φ_{1}(R) ≥ 2^{ν−1}Γ(ν)R^{−2ν}e^{−R}.

Then we have

φ_{0}(M) = \frac{2π^{d/2}}{Γ (−\frac{β}{2} )} (2γ^{2})^{ν}φ_{1}(2γM)

≥ \frac{2π^{d/2}}{Γ (−\frac{β}{2} )} (2γ^{2})^{ν}2^{ν−1}Γ(ν)e^{−2γM}(2γM)^{−2ν}

= \frac{π^{d/2}Γ(ν)M^{−2ν}e^{−2γM}}{Γ (−\frac{β}{2} )}

and

λ ≥ \frac{1}{2^{2d+1}} · \frac{Γ (\frac{d + β}{2} )}{Γ (−\frac{β}{2} ) Γ (\frac{d}{2} + 1)} M^{−β} exp(−2γM)

as incorporated into Table 2 with M = 6.38d/q. For β < 0 and large values of q there is a better choice of M, but we leave this exotic case to the reader. Similar tricks can be done in the Gaussian and Sobolev cases.

Acknowledgement: The author thanks the referees for a series of helpful remarks and suggestions.

12

References

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

[2] K. Ball. Eigenvalues of Euclidean distance matrices. Journal of Approximation Theory, 68:74–82, 1992.

[3] K. Ball, N. Sivakumar, and J.D. Ward. On the sensitivity of radial basis interpolation to minimal distance separation. Journal of Approximation Theory, 8:401–426, 1992.

[4] B.J.C. Baxter, N. Sivakumar, and J.D. Ward. Regarding the p–norms of radial basis interpolation matrices. Technical Report 281, Department of Mathematics, Texas A & M University, 1992.

[5] M.D. Buhmann. New developments in the theory of radial basis function interpolation. In Multivariate Approximations: From CAGD to Wavelets, pages 35–75. K. Jetter and F.I. Utreras, editors; World Scientific, London, 1993.

[6] N. Dyn. Interpolation and approximation by radial and related functions. In C.K. Chui, L.L. Schumaker, and J.D. Ward, editors, Approximation Theory VI. Academic Press, Boston, 1989.

[7] N. Dyn, W.A. Light, and E.W. Cheney. Interpolation by piecewise–linear radial basis functions, I. Journal of Approximation Theory, 59:202–223, 1989.

[8] D.S. Jones. The Theory of Generalized Functions. Cambridge University Press, 1982.

[9] W.R. Madych and S.A. Nelson. Multivariate interpolation: a variational theory. Manuscript, 1983.

[10] W.R. Madych and S.A. Nelson. Bounds on multivariate polynomials and exponential error estimates for multiquadric interpolation. Journal of Approximation Theory, 70:94– 114, 1992.

[11] F.J. Narcowich and J.D. Ward. Norm of inverses and condition numbers for matrices associated with scattered data. Journal of Approximation Theory, 64:69–94, 1991.

[12] F.J. Narcowich and J.D. Ward. Norms of inverses for matrices associated with scattered data. In P.J. Laurent, A. Le Méhauté, and L.L. Schumaker, editors, Curves and Surfaces, pages 341–348. Academic Press, Boston, 1991.

[13] F.J. Narcowich and J.D. Ward. Norm estimates for the inverses of a general class of scattered–data radial–function interpolation matrices. Journal of Approximation Theory, 69:84–109, 1992.

[14] M.J.D. Powell. Univariate multiquadric interpolation: Some recent results. In P.J. Laurent, A. Le Méhauté, and L.L. Schumaker, editors, Curves and Surfaces, pages 371–382. Academic Press, 1991.

[15] R. Schaback. Comparison of radial basis function interpolants. In Multivariate Approximation. From CAGD to Wavelets, pages 293–305. K. Jetter and F. Utreras, editors; World Scientific, London, 1993.

[16] R. Schaback. Lower bounds for norms of inverses of interpolation matrices for radial basis functions. To appear in Journal of Approximation Theory, 1994.

[17] X. Sun. Norm estimates for inverses of Euclidean distance matrices. Journal of Approximation Theory, 70:339–347, 1992.

[18] Z. Wu and R. Schaback. Local error estimates for radial basis function interpolation of scattered data. IMA Journal of Numerical Analysis, 13:13–27, 1993.

Blatt bitte entfernen. Liegt nur bei, damit die Referenzierung fuer Tabelle 2 (die ja extra gedruckt wird) klappt.

Table 2: These table entries explicitly contain the relevant constants.