Lower Bounds for Norms of Inverses of Interpolation Matrices for Radial Basis Functions

Robert Schaback

Abstract: Interpolation of scattered data at distinct points x_{1},…, x_{n} ∈ IR^{d} by linear combinations of translates Φ(∥x − x_{j}∥_{2}) of a radial basis function Φ : IR_{≥0} → IR requires the solution of a linear system with the n by n distance matrix A := (Φ(∥x_{i} − x_{j}∥_{2}). Recent results of Ball, Narcowich and Ward, using Laplace transform methods, provide upper bounds for ∥A^{−1}∥_{2}, while Ball, Sivakumar, and Ward constructed examples with regularly spaced points to get special lower bounds. This paper proves general lower bounds by application of results of classical approximation theory. The bounds increase with the smoothness of Φ. In most cases, they leave no more than a factor of n^{−2} to be gained by optimization of data placement, starting from regularly distributed data. This follows from comparison with results of Ball, Baxter, Sivakumar, and Ward for points on scaled integer lattices and supports the hypothesis that regularly spaced data are near-optimal, as far as the condition of the matrix A is concerned.

Keywords: Radial basis functions, multiquadrics, scattered data interpolation, lower bounds, condition numbers.

Classification: 65D05, 41A05, 41A63

1 Introduction

Let Φ : IR_{≥0} → IR be a scalar (“radial”) function, and let n distinct points (“centres”) x_{1},…, x_{n} ∈ IR^{d} be given, forming a set X := {x_{1},…, x_{n}}. As reported by Hardy [8] and Franke [7], interpolation of real values y_{i}, 1 ≤ i ≤ n, at the centres x_{i} by linear combinations

s(x) := ∑_{j=1}^{n} α_{j}Φ(∥x − x_{j}∥_{2}), α_{j} ∈ IR, x ∈ IR^{d}  (1.1)

of translates of Φ(∥.∥_{2}) can produce very good numerical results, if Φ is the “multiquadric” Φ_{c}(r) = (c^{2} + r^{2})^{1/2}, for instance. The interpolation problem for a function of the form (1.1) requires the solution of the linear system

∑_{j=1}^{n} α_{j}Φ(∥x_{i} − x_{j}∥_{2}) = y_{i}, 1 ≤ i ≤ n.  (1.2)

In some cases the n functions Φ(∥ · −x_{j}∥_{2}) are augmented by a basis of the space IP_{m}^{d} of all d-variate polynomials of total degree less than m. Let

μ := μ(m, d) := \binom{m − 1 + d}{d}, m ≥ 1, μ(0, d) := 0  (1.3)

0File : campaper.tex, TEXed May 18, 1994, Status : Revised for J. Approx. Th.

denote the dimension of IP_{m}^{d}, and let a basis be given by q_{1},…, q_{μ}. From now on, we suppress the dependence of μ on m and d to keep the notation simple. For m > 0 we additionally assume q_{1},…, q_{μ} to be linearly independent over the set X = {x_{1},…, x_{n}} of centres, which implies n ≥ μ. Then the (n + μ) by (n + μ) matrix

A := ( \begin{matrix}Φ(∥x_{i} − x_{j}∥_{2}) & q_{k}(x_{i}) \\ q_{k}(x_{j}) & 0 \end{matrix} )_{1≤i,j≤n, 1≤k≤μ}  (1.4)

which occurs in the generalization

\begin{aligned}∑_{j=1}^{n} α_{j}Φ(∥x_{i} − x_{j}∥_{2}) & + & ∑_{k=1}^{μ} β_{k}q_{k}(x_{i}) & = & y_{i}, & 1 ≤ i ≤ n \\ ∑_{j=1}^{n} α_{j}q_{k}(x_{j}) & + & 0 & = & 0, & 1 ≤ k ≤ μ \end{aligned}  (1.5)

of the system (1.2) is called a (generalized) distance matrix. It is nonsingular in the cases

Φ(r) = (c^{2} + r^{2})^{β}, β, c > 0, β ∉ 2ZZ, m ≥ 0

(Multiquadrics for β = 1/2)

Φ(r) = (c^{2} + r^{2})^{−β}, β, c > 0, m ≥ 0

(Inverse Multiquadrics for β = 1/2)

Φ(r) = r^{2β} log r,  β = m − d/2 > 0, d ∈ 2IN

(Thin-plate splines)

Φ(r) = log(c^{2} + r^{2}), c > 0, m ≥ 0
Φ(r) = r^{2β},  β = m − d/2, m > d/2, d ∈ 2IN − 1, or

β ∈ (0, 1), d ∈ 2IN − 1, m ≥ 0, or β = 1/2, n ≥ 2, m ≥ 0

Φ(r) = exp(−αr^{2}), α > 0, m ≥ 0, (Gaussians).

This follows from the work of Micchelli [9] and Schoenberg [15] (see also Dyn [5] and Powell [14] for highly useful surveys of known results). Furthermore, condition numbers of A were often observed to be quite large. Special preconditioning strategies for solving the system (1.5) were supplied by Dyn, Levin, and Rippa [6] in a variety of special cases.

Some very interesting upper bounds for the spectral norm ∥A^{−1}∥_{2} of A were given by Ball [1], and Narcowich and Ward [10] [11] [12], using Laplace transform methods. These bounds are in terms of the “separation distance”

h := min_{1≤i≠j≤n} ∥x_{i} − x_{j}∥_{2}  (1.6)

and hold for arbitrary placements of centres x_{j}, 1 ≤ j ≤ n. Recently, a paper by Ball, Sivakumar, and Ward [2] derived lower bounds for ∥A^{−1}∥_{2} for special regular choices of centres by a similar technique. Furthermore, Baxter [4] has investigated the case of centres on subsets of the integer lattice. Using the Toeplitz structure of A, he derived bounds for ∥A^{−1}∥_{2} for a wide class of radial basis functions including multiquadrics. His bounds are asymptotically optimal in the sense that their behaviour for n → ∞ is best possible up to a constant factor.

Since the existing lower bounds of ∥A^{−1}∥_{2} only hold for data on finite regular grids, the value of ∥A^{−1}∥_{2} might be decreased by irregular placements of centres.

This paper complements the work of Ball, Baxter, Narcowich, Sivakumar, and Ward by producing general lower bounds for ∥A^{−1}∥_{2}, independent of the distribution of centres and independent of the separation distance, via a completely different approach. Combined with the other results, they provide a bound for the possible gain by optimizing the placement of centres with respect to the condition of A. The findings of this paper are roughly of the form

∥A^{−1}∥_{2} ≥ \frac{1}{nE(2(d!n/2)^{1/d} − 1, K, Φ)},

where E(ℓ, K, Φ) denotes the error of best Chebyshev approximation by polynomials of degree less than ℓ to the function Φ(\sqrt{r}) on the interval [0, K^{2}] defined via the diameter

K := max_{1≤i,j≤n}∥x_{i} − x_{j}∥_{2}  (1.7)

of the data set in the Euclidean norm. This relates the bound to the smoothness of Φ and the dimension d of the space, and in a very remarkable way indeed: the bound tends to be smaller, if the smoothness of the radial basis functions is decreased. The results can be generalized to other matrix norms than the spectral norm to which the other bounds in the literature are confined.

2 Basic Results

We assume that the radial function Φ, the n ≥ 2 centres x_{j} ∈ X := {x_{1},…, x_{n}}, the space dimension d, and the polynomial order m are such that the matrix A in (1.4) is nonsingular. For an arbitrary polynomial p ∈ IP_{ℓ}^{1} we define the matrix

A_{p} := ( \begin{matrix}p(∥x_{i} − x_{j}∥_{2}^{2}) & q_{k}(x_{i}) \\ q_{k}(x_{j}) & 0 \end{matrix} )_{1≤i,j≤n, 1≤k≤μ}  (2.1)

as an approximation to A. Note that A and A_{p} differ only in their upper n by n submatrices with entries Φ(∥x_{i} −x_{j}∥_{2}) and p(∥x_{i} −x_{j}∥_{2}^{2}), respectively. Thus a good approximation of Φ(\sqrt{r}) by a polynomial p ∈ IP_{ℓ}^{1} on the set T of all real values t_{ij} := ∥x_{i} − x_{j}∥_{2}^{2} will produce a good approximation of A by A_{p}.

We now pick two arbitrary norms ∥.∥_{r} and ∥.∥_{s} on IR^{n+μ} and define the usual operator norm

∥B∥_{r,s} := sup { ∥Bx∥_{r}/∥x∥_{s} | x ∈ IR^{n+μ} \ {0} } .  (2.2)

By the theorem of Weierstrass, there will be some p ∈ IP_{ℓ}^{1} for ℓ large enough such that

∥A − A_{p}∥_{s,r} · ∥A^{−1}∥_{r,s} < 1  (2.3)

holds. We now assert that (2.3) is a sufficient condition for A_{p} to be nonsingular. In fact, for arbitrary vectors x, y with A_{p}x = y we find

\begin{aligned}∥x∥_{r} & = & ∥A^{−1}Ax∥_{r} = ∥A^{−1}(y + (A − A_{p})x)∥_{r} \\ & ≤ & ∥A^{−1}∥_{r,s}∥y∥_{s} + ∥A^{−1}∥_{r,s}∥A − A_{p}∥_{s,r}∥x∥_{r}, \end{aligned}

and reordering yields a bound for ∥A_{p}^{−1}∥_{r,s}.

Now, whenever A_{p} is nonsingular, the n + μ polynomials

\begin{matrix}p(∥ · −x_{j}∥_{2}^{2}), & 1 ≤ j ≤ n = |X| \\ q_{k}(·), & 1 ≤ k ≤ μ \end{matrix}  (2.4)

span the n-dimensional space C(X), because the n first rows or columns of A_{p} are linearly independent. For any univariate polynomial p ∈ IP_{ℓ}^{1} we can define a subspace

P_{p} := {∑_{k=1}^{μ}β_{k}q_{k} + ∑_{t=1}^{n}γ_{t}p(∥ · −x_{t}∥_{2}^{2}) : ∑_{t=1}^{n}γ_{t}q_{k}(x_{t}) = 0, 1 ≤ k ≤ μ}  (2.5)

of the space IP_{2ℓ−1}^{d}, and whenever A_{p} is nonsingular, we have dim P_{p} ≥ n, because the polynomials in (2.4) span the n-dimensional space C(X). If we define

μ^{∗}(ℓ) := μ^{∗}(X, ℓ, m, d) := max_{p∈IP_{ℓ}^{1}} dim P_{p},  (2.6)

we conclude that

μ^{∗}(ℓ) ≥ |X| = n  (2.7)

holds, if ℓ is large enough to satisfy (2.3) for some p ∈ IP_{ℓ}^{1}.

We now turn this argument upside down. Clearly the function μ^{∗}(ℓ) is (weakly) monotonic in ℓ and has the obvious bounds

μ ≤ μ^{∗}(ℓ) ≤ |X| + μ.

Our argument above implies existence of the maximum in the definition of

ℓ^{∗} := ℓ^{∗}(X, m, d) := max {ℓ ≥ 0 | μ^{∗}(ℓ) < |X|} .

This yields

μ^{∗}(ℓ^{∗}) < |X| = n,  (2.8)

which will be needed in the proof of

Theorem 2.1 The inverse of A, if it exists, satisfies

∥A^{−1}∥_{r,s}^{−1} ≤ inf {∥A − A_{p}∥_{s,r} | p ∈ IP_{ℓ^{∗}(X,m,d)}^{1}}  (2.9)

where A_{p} is defined by (2.1).

Proof: Take any polynomial p ∈ IP_{ℓ^{∗}}^{1} and assume A_{p} to be nonsingular. Then the n + μ polynomials (2.4) span the space C(X). Since they are in the space P_{p} occurring in (2.6), the inequality (2.8) is violated. Thus A_{p} must be singular for all p ∈ IP_{ℓ^{∗}}^{1}. But then the inequality (2.3) cannot hold because it would imply the nonsingularity of A_{p}, as was shown above. This proves the assertion. □

Theorem 2.1 relates lower bounds for ∥A^{−1}∥_{r,s} to a somewhat peculiar matrix-valued approximation problem. The error matrix is zero except for the upper left n by n submatrix of A − A_{p} with entries

Φ(∥x_{i} − x_{j}∥_{2}) − p(∥x_{i} − x_{j}∥_{2}^{2}), 1 ≤ i, j ≤ n.

Thus it involves an approximation of Φ(\sqrt{r}) by polynomials of order at most ℓ^{∗} on the set T of all real values t_{ij} := ∥x_{i} − x_{j}∥_{2}^{2} such that the function values are arranged in n × n matrix form and such that the approximation error is measured via the matrix norm ∥.∥_{s,r} from (2.2). The actual form of the approximation problem is thus determined by the matrix norm chosen. For the spectral norm ∥A − A_{p}∥_{2,2} we get a rather nasty approximation problem, but for the norm

∥B∥_{∞,1} = max_{1≤i,j≤n}|b_{ij}|  (2.10)

we can use Chebyshev approximation on T. Here, we interpreted the norms ∥.∥_{r} and ∥.∥_{s} of (2.2) as the usual L_{r} and L_{s} norms, but this is by no means mandatory. Studying linear approximation problems for peculiar norms like the spectral matrix norm may be of independent interest in approximation theory.

We want to give good asymptotic lower bounds for ∥A^{−1}∥_{r,s} in case n → ∞, and we want to compare the bounds with those of Ball, Baxter, Narcowich, Sivakumar, and Ward for the spectral norm ∥A^{−1}∥_{2,2}. We saw the latter to be rather inconvenient for our approach while good asymptotics are mainly available for Chebyshev approximation on real intervals. Therefore we try to get as far as we can with best Chebyshev approximation to the function Φ(\sqrt{r}) on the real interval [0, K^{2}] by polynomials of order ℓ. Note that this is equivalent to a best Chebyshev approximation of order 2ℓ − 1 to the function Φ(|r|) on [−K, +K], which, by uniqueness and symmetry, must be a polynomial of maximal order ℓ in r^{2}. Then we denote the error by

E(ℓ) := E(ℓ, K, Φ) := min_{p∈IP_{ℓ}^{1}} max_{0≤r≤K^{2}}|Φ(\sqrt{r}) − p(r)|
= min_{p∈IP_{ℓ}^{1}}∥Φ(|r|) − p(r^{2})∥_{∞,[−K,K]},

writing the Chebyshev norm of continuous functions on a real interval [a, b] by ∥.∥_{∞,[a,b]}.

Theorem 2.2 The inverse of A, if it exists, satisfies

∥A^{−1}∥_{1,∞} ≥ \frac{1}{E(ℓ^{∗})}  (2.11)

and

∥A^{−1}∥_{2,2} ≥ \frac{1}{n · E(ℓ^{∗})}.  (2.12)

Proof: Inequality (2.11) readily follows from Theorem 2.1 and (2.10). Then we conclude (2.12) simply from (2.11) and

∥B∥_{1,∞} ≤ n · ∥B∥_{2,2}

for arbitrary n by n matrices B. □

Unfortunately, our approach yields a factor of n^{−1} in the bound for the spectral norm; the comparisons at the end of the paper will indicate that we often seem to miss the actual behaviour of ∥A^{−1}∥_{2,2} for regular data asymptotically by just this factor. But there appear to be no other handy links between the spectral norm for matrices and the Chebyshev norm for the matrix entries. The factor does not arise if we use (2.11) or replace E in (2.12) by the (unknown) error of best approximation by polynomials in the spectral norm.

Theorem 2.2 relates lower bounds of ∥A^{−1}∥ to the error of best Chebyshev approximations to Φ(|r|) by polynomials of order ℓ^{∗} in the variable r^{2}, and we shall see below that ℓ^{∗} behaves approximately like at least n^{1/d} for large n. Altogether, our lower bounds will grow astronomically with n, if Φ(|r|) can be extended to an entire function in the complex plane C (e.g.: for Gaussians). They will be of polynomial growth, if Φ(|r|) has only finitely many continuous derivatives, which is the case for thin–plate splines and positive non-even rational powers of r. However, they will still grow exponentially for all Φ(|r|) which are analytic in C around [0, K], e.g.: for multiquadrics and inverse multiquadrics, if c is fixed. To get around this, c must tend to zero for n tending to infinity, moving the singularity of Φ towards the real axis. The details are worked out in the rest of the paper, but due to classical Jackson–Bernstein theorems the lower bounds for ∥A^{−1}∥ are directly related to the smoothness of Φ: they get larger, if the smoothness increases.

3 Auxiliary Results

The remaining task now is to evaluate E(ℓ^{∗}(X, m, d), K, Φ) for the functions listed in the introduction. Since the quantities E(ℓ, K, Φ) can be estimated using classical results of approximation theory (this will be done in the following section), we first prove something about μ^{∗}(ℓ) and ℓ^{∗} as defined in the beginning of the preceding section. For this, we suppose X, m, and d to be fixed, and we first look for an upper bound for μ^{∗} from (2.6).

To derive simple results, the crude bound

μ^{∗}(ℓ) ≤ dim IP_{2ℓ−1}^{d} = μ(2ℓ − 1, d)

for 2ℓ > m suffices. It follows easily from (2.5) and uses the general definition of μ(·, d) from (1.3). But there is an improvement of this bound:

Lemma 3.1 Let p be a univariate polynomial of order ≤ ℓ, where ℓ ≥ m+1. Then the subspace P_{p}, as defined in (2.5), is a subspace of IP_{2ℓ−1−m}^{d} and has dimension at most

dim IP_{ℓ}^{d} + dim IP_{ℓ−1}^{d} − dim IP_{m}^{d} = μ(ℓ, d) + μ(ℓ − 1, d) − μ(m, d).  (3.1)

Furthermore,

μ^{∗}(ℓ) ≤ μ(ℓ, d) + μ(ℓ − 1, d) − μ(m, d).

Proof: By expansion of the translate

p_{y}(x) := p(∥x − y∥_{2}^{2}) = ∑_{k=0}^{ℓ−1} α_{k}∥x − y∥_{2}^{2k}

of p into sums of products of simpler terms we get the representation

\begin{aligned}p_{y}(x) & = & ∑_{k=0}^{ℓ−1}α_{k}(x^{T} x − 2x^{T} y + y^{T} y)^{k} \\ & = & ∑_{j=0}^{ℓ−1}(−2x^{T} y)^{j} ∑_{i=0}^{ℓ−j−1}(x^{T} x)^{i} ∑_{k=i+j}^{ℓ−1}α_{k} ( \begin{matrix}k \\ j \end{matrix} ) ( \begin{matrix}k − j \\ k − j − i \end{matrix} ) (y^{T} y)^{k−j−i}. \end{aligned}

Thus an arbitrary function taken from P_{p} in (2.5) can be written as

∑_{k=1}^{μ}β_{k}q_{k}(x) + ∑_{t=1}^{n}γ_{t}p(∥x − x_{t}∥^{2})

= ∑_{k=1}^{μ}β_{k}q_{k}(x)

+ ∑_{t=1}^{n}γ_{t}∑_{j=0}^{ℓ−1}(−2x^{T} x_{t})^{j} ∑_{i=0}^{ℓ−j−1}(x^{T} x)^{i} ∑_{k=i+j}^{ℓ−1} ( \begin{matrix}k \\ j \end{matrix} ) ( \begin{matrix}k − j \\ k − j − i \end{matrix} ) α_{k}(x_{t}^{T} x_{t})^{k−j−i}

= ∑_{k=1}^{μ}β_{k}q_{k}(x)

+ ∑_{j=0}^{ℓ−1} ∑_{i=0}^{ℓ−j−1}(x^{T} x)^{i} ∑_{k=i+j}^{ℓ−1} ( \begin{matrix}k \\ j \end{matrix} ) ( \begin{matrix}k − j \\ k − j − i \end{matrix} ) α_{k}∑_{t=1}^{n}γ_{t}(−2x^{T} x_{t})^{j}(x_{t}^{T} x_{t})^{k−j−i}

and we note that the sum over t yields zero whenever j + 2(k − i − j) ≤ m − 1. Thus the range of i and j can be restricted to m + j + 2i ≤ 2ℓ − 2, and the function will be contained in the subspace

Q = IP_{m}^{d} + span \{(x^{T} y)^{j}(x^{T} x)^{i} | \begin{matrix}0 ≤ j ≤ ℓ − 1, 0 ≤ i ≤ ℓ − 1 − j \\ m + j + 2i ≤ 2ℓ − 2, y ∈ IR^{d} \end{matrix} \}

of IP_{2ℓ−1−m}^{d}. But the number of spanning functions can be further reduced. We assert that Q lies in the subspace

Q̃ = IP_{ℓ}^{d} + span {(x^{T} y)^{j}(x^{T} x)^{ℓ−1−j} | m ≤ j ≤ ℓ − 2, y ∈ IR^{d}} .  (3.2)

The terms (x^{T} x)^{i}(x^{T} y)^{j} for 2i + j ≤ ℓ − 1 lie directly in IP_{ℓ}^{d}. Those with i + j = ℓ − 1 that are not in IP_{ℓ}^{d} will have j ≤ ℓ − 2, and the restriction m + j + 2i ≤ 2ℓ − 2 implies j ≥ m. Thus these terms clearly lie in Q̃. For the remaining cases we can define r := ℓ − 1 − i − j > 0 and split the terms as

(x^{T} x)^{i}(x^{T} y)^{j} = (x^{T} x)^{i−r}(x^{T} x)^{r}(x^{T} y)^{j}.

We now use that each factor (x^{T} x)^{r}(x^{T} y)^{j}, being a homogeneous polynomial in IP_{2r+1+j}^{d}, can be represented via terms of type (x^{T} y)^{2r+j}. This leads to terms of the form (x^{T} x)^{i−r}(x^{T} y)^{2r+j}, and these are in Q̃, as an easy check of the exponents reveals.

We now have to bound the dimension of Q̃ from above. If H_{j}^{d} denotes the space of homogeneous polynomials on IR^{d} of order j, then the space IP_{ℓ}^{d} is representable as

IP_{ℓ}^{d} = span {(x^{T} y)^{j} | 0 ≤ j ≤ ℓ − 1, y ∈ IR^{d}}
= span {H_{j}^{d} | 1 ≤ j ≤ ℓ } ,

(3.3)

and we shall consider the excess of Q̃ over IP_{ℓ}^{d}. The representation (3.2) can be rearranged as

Q̃ = IP_{ℓ}^{d} + span {H_{j+1}^{d}(x^{T} x)^{ℓ−1−j} | m ≤ j ≤ ℓ − 2} ,

and we can repeat the argument used for (3.3) to conclude that the dimension of the second part is bounded from above by dim IP_{ℓ−1}^{d} − dim IP_{m}^{d}. This proves (3.1), and the second assertion of the lemma is obvious because Q̃ contains IP_{m}^{d} for m ≤ ℓ. □

If we define

ℓ^{∗∗} := ℓ^{∗∗}(n, m, d) := max {ℓ ≥ 0 | μ(ℓ, d) + μ(ℓ − 1, d) − μ(m, d) < n} ,  (3.4)

we have

m ≤ ℓ^{∗∗} ≤ ℓ^{∗}.  (3.5)

Note that ℓ^{∗∗}(n, m, d) is monotonic with respect to n. For ℓ → ∞ and d fixed, clearly

μ(ℓ, d) = \frac{ℓ^{d}}{d!} + O(ℓ^{d−1}),

which implies the asymptotic behaviour

ℓ^{∗∗}(n, m, d) = (\frac{d!}{2}n)^{1/d} + O(1)  (3.6)

for n = |X| → ∞ and fixed values of m and d. This will be sufficient for later use, because (3.5) allows us to work with ℓ^{∗∗}(n, m, d) and (3.6) instead of ℓ^{∗}(X, m, d) for large n. Since we are not interested in cases with very small values of n, we pay no further attention to the unimportant restrictions n ≥ max(2, μ(m, d))and ℓ ≥ m + 1.

In the one-dimensional case we easily get

\begin{aligned}μ^{∗}(X, ℓ, m, 1) & = & 2ℓ − 1 for 0 ≤ m ≤ 2ℓ − 1 \\ ℓ^{∗}(X, m, 1) & = & ⌊\frac{|X|+m}{2}⌋, \end{aligned}

while for d = 2 we find

\begin{aligned}μ^{∗}(X, ℓ, m, 2) & ≤ & ℓ^{2} for 0 ≤ m ≤ ℓ \\ ℓ^{∗∗}(|X|, m, 2) & = & ⌊\sqrt{|X| − 1 + \frac{m(m+1)}{2}}⌋. \end{aligned}

4 Approximation Orders

To get bounds for E(ℓ, K, Φ) we apply some classical results of approximation theory. We start with

Theorem 4.1 Let f ∈ C[a, b] be a real-valued function with a holomorphic extension to an ellipse in the complex plane with foci a and b, and let f have a singularity at the boundary of the ellipse. If ρ = R(b − a)/2 is the sum of the two half-axes of the ellipse, then

\overline{lim}_{ℓ→∞} _{ℓ−1} \sqrt{inf_{p∈IP_{ℓ}^{1}}∥f − p∥_{∞,[a,b]}} = \frac{1}{R}.

If, in addition, the absolute value of the real part ℜ(f) of f is bounded by 1 on the ellipse, and if we set a = −1, b = +1 for simplicity, we have

inf_{p∈IP_{ℓ}^{1}}∥f − p∥_{∞,[−1,+1]} ≤ \frac{8}{πρ^{ℓ}} .  (4.1)

Proof: The first part is due to Bernstein and can be found in [13], p. 194. The refined statement is proven in [16], p. 203. □

To cope with functions like Φ(r) = r^{2β} we apply

Theorem 4.2 For any function f ∈ C^{k}[a, b] with a k-th derivative in the class Lip_{M}α we have

inf_{p∈IP_{ℓ}^{1}}∥f − p∥_{∞,[a,b]} ≤ \frac{c_{k}(b − a)^{k+α}}{ℓ^{k+α}} M  (4.2)

where

c_{k} = 12\frac{6^{k}k^{k}}{k!} (\frac{k + 1}{2})^{α} .

If the k-th derivative is continuous, but not contained in any Lipschitz class, then still

inf_{p∈IP_{ℓ}^{1}}∥f − p∥_{∞,[a,b]} ≤ \frac{c}{ℓ^{k}}  (4.3)

with a constant c that depends on f, k, and [a, b], but not on ℓ.

Proof: This is due to Jackson, see e.g.: [13], p.128.

5 Application to Multiquadrics

We first treat the case of functions Φ_{c}(r) = f(c^{2} + r^{2}) on [−1, 1] for c > 0, where f is analytic in C except for a singularity at the origin. The regularity ellipse of Φ_{c} then has half-axes of length \sqrt{1 + c^{2}} and c, yielding ρ = c + \sqrt{1 + c^{2}}. If |ℜΦ_{c}| is bounded by C_{c} on the ellipse, then (4.1) implies

E(ℓ, 1, f(c^{2} + r^{2})) ≤ \frac{8}{πρ^{2ℓ−1}} C_{c},

if we approximate Φ_{c}(r) = f(c^{2} + r^{2}) by p ∈ IP_{2ℓ−1}^{1} and write the result as a polynomial in r^{2} of order at most ℓ. For functions Φ_{c}(r) = (c^{2} + r^{2})^{β} on [−1, 1] with β > 0 we find

E(ℓ, 1, (c^{2} + r^{2})^{β}) ≤ \frac{8(1 + 2c^{2})^{β}}{πρ^{2ℓ−1}}

and the general case of Φ_{c}(r) = (c^{2} + r^{2})^{β} on [−K, +K] is easily recovered as

E(ℓ, K, (c^{2} + r^{2})^{β}) ≤ K^{2β} \frac{8(1 + 2γ^{2})^{β}}{πρ^{2ℓ−1}}

with the scaled quantities

γ = c/K, ρ = γ + \sqrt{1 + γ^{2}},

indicating that the relative size γ of c and K is crucial.

For β < 0 the real part of Φ_{c}(r) = (c^{2} + r^{2})^{β} is unbounded on the regularity ellipse. Thus we have to use the weaker form of Theorem 4.1 to get

E(ℓ, K, (c^{2} + r^{2})^{β}) ≤ \frac{C}{ρ^{2ℓ−1}}

on [−K, +K] with a constant C that will depend on c, K, γ, and β, but not on ℓ. Similar estimates hold for the radial basis function Φ_{c}(r) = log(c^{2} + r^{2}).

The value ρ^{j} can be bounded using

ρ = γ + \sqrt{1 + γ^{2}} ≥ 1 + γ ≥ (1 − γ/2)^{−1}

for 0 ≤ γ ≤ 1 to get

ρ^{j} ≥ exp(jγ/2).

Inserting everything into (2.12) we get the lower bound

∥A^{−1}∥_{2,2} ≥ \frac{π exp((2ℓ^{∗∗} − 1)γ/2)}{8nK^{2β}} \{ \begin{matrix}(1 + 2γ^{2})^{−β} & for & β > 0 \\ C(c, K, γ, β) & for & β < 0 \end{matrix} \}  (5.1)

with exponential growth for fixed γ ∈ [0, 1], where ℓ^{∗∗} is given by (3.4). Note that (3.6) implies

∥A^{−1}∥_{2,2} ≥ O (\frac{1}{n} exp (\frac{γ}{2} (2 (\frac{d!}{2}n)^{1/d} − 1)))  (5.2)

for n → ∞ and fixed γ, d. The theory for infinite grids (see Buhmann [3]) varies c proportional to the minimal distances of centres. For a fixed finite domain in IR^{d} with centres forming a dense subset, this strategy is resembled by letting γ vary like n^{−1/d}. Then the behaviour for n → ∞ is

∥A^{−1}∥_{2,2} ≥ O (\frac{1}{n})

for β > 0, because the exponential function in (5.2) gets a constant argument due to cancellation.

6 Application to other radial basis functions

For Φ(r) = r^{2β} with 2β = 2p + 1, p ∈ IN_{≥0} we again approximate Φ(|r|) on [−K, K] by polynomials in IP_{2ℓ−1}^{1} and take advantage of the symmetry of the best approximation. Since the derivatives up to order 2p = 2β − 1 are continuous and the k := 2p-th derivative is in a Lipschitz class Lip_{C}α with α = 1 and C := (2β)! we can invoke Theorem 4.2 to get

E(ℓ, K, Φ) ≤ \frac{c_{2p}(2K)^{2β}}{(2ℓ − 1)^{2β}} (2β)!

and the lower bound

∥A^{−1}∥_{2,2} ≥ \frac{(2ℓ^{∗∗} − 1)^{2β}}{nc_{2p}(2K)^{2β}(2β)!}

with the asymptotic behaviour

∥A^{−1}∥_{2,2} ≥ O (n^{\frac{2β}{d} −1})

for n → ∞.

A similar analysis holds for thin-plate splines of the form Φ(r) = r^{2β} log r with β ∈ IN_{≥1}. Derivatives up to order 2β − 1 are continuous, and for ℓ^{∗∗} ≥ β there is no problem with the polynomials that arise when taking derivatives. However, the highest continuous derivative is not contained in any Lipschitz class. Application of the restricted statement of Theorem 4.2 then implies

∥A^{−1}∥_{2,2} ≥ C(β, K)\frac{(2ℓ^{∗∗} − 1)^{2β−1}}{n} = O (n^{\frac{2β−1}{d} −1})

for n → ∞ where C(β, K) is a constant.

For Gaussians Φ(r) = exp(−α^{2}r^{2}) we can use the Taylor expansion to get

E(ℓ, K, Φ) ≤ \frac{(αK)^{ℓ}}{ℓ!} ,

providing a disastrous lower bound on ∥A^{−1}∥_{2,2} for fixed values of αK if n is large.

Remark. A variation of our approach would be to approximate Φ(\sqrt{r}) on [h^{2}, K^{2}] for the separation distance (1.6) as used by Ball, Narcowich, and Ward. This would make any of the classical radial basis functions analytic in a neighbourhood of [h^{2}, K^{2}], leading to exponential decrease of the approximation error E(ℓ, [h^{2}, K^{2}], Φ) with respect to polynomials of order ℓ for ℓ → ∞. The basic estimate has the same form as (5.2), but a detailed analysis reveals that the exponential term of this bound gets a constant argument because of K ≥ n^{1/d}h and ℓ ≈ n^{1/d}. This corresponds to the cancellation in (5.2) when γ or c are scaled to decrease with the minimal separation distance.

For illustration, consider the function (r^{2})^{β} = t^{β} for t = r^{2} ∈ [h^{2}, K^{2}]. After rescaling, it coincides with the multiquadric

t^{β} = (h^{2} + s\frac{K^{2} − h^{2}}{K^{2}} )^{β} = (\frac{K^{2} − h^{2}}{K^{2}} )^{β} (\frac{h^{2}K^{2}}{K^{2} − h^{2}} + s)^{β}

with exponent β and c = hK/\sqrt{K^{2} − h^{2}} for s ∈ [0, K^{2}]. Thus the behaviour of Φ(r) = r^{2β}, as far as our lower bounds are concerned, is roughly the same as for multiquadrics with exponent β and a scaling of c ≈ h ≈ γ ≈ n^{−1/d}, where the exponential in the bound (5.2) gets a constant argument due to cancellation. For multiquadrics themselves, the effect of introducing h just acts like a corresponding increase of c and does not yield any improvement.

7 Comparison with other bounds

In general, the upper bounds by Ball [1] and Narcowich, and Ward [10][11][12] provide estimates of ∥A^{−1}∥_{2,2} from above in terms of the separation distance h of (1.6). The corresponding lower bounds of Ball, Sivakumar, and Ward [2] hold for a specific regularly distributed set of centres and thus act as strict lower bounds for the worst possible data set with prescribed separation distance h.

In contrast to this, our lower bounds provide best–case estimates of ∥A^{−1}∥_{2,2} from below, because they hold for every distribution of the data points, including the best possible choice, if the latter should exist. These bounds must necessarily be smaller than the lower bounds of Ball, Sivakumar, and Ward, the difference being leeway for optimizing the placement of centres.

We start with a comparison of results for multiquadrics on increasing sets of centres with separation distance 1 in IR^{d}. In this case, the optimal bounds of Baxter [4] for ∥A^{−1}∥_{2,2} must lie between our lower bounds and the upper bounds of Ball, Narcowich, and Ward. Whenever Baxter’s bounds coincide with the latter, the conclusion is that regular data asymptotically realize the worst possible distribution with separation distance 1. The difference between Baxter’s bounds and ours may possibly be used for optimization of placements of centres, because we do not make assumptions on the separation distance or regularity of distribution.

The optimal bound of Baxter for ∥A^{−1}∥_{2,2} in the case of centres on integer grids takes the form

∥A^{−1}∥_{2,2} ≤ \frac{π}{4c} \frac{1}{e^{−cπ} + \frac{1}{3}e^{−3cπ} + \frac{1}{5}e^{−5cπ}…} ≈ \frac{π}{4c}e^{cπ}

for d = 1 and for multiquadrics Φ_{c}(r) = (c^{2} + r^{2})^{1/2}. The same value arises as the precise limit of ∥A^{−1}∥_{2,2} for n → ∞, when the n integer points i ∈ IR with 0 ≤ i ≤ n − 1 are taken.

Ball, Sivakumar, and Ward [2] get

∥A^{−1}∥_{2,2} ≥ C e^{c\sqrt{d}}

for n → ∞ and the same data distribution, where the constant C does not depend on c and d. The worst–case upper bound of Narcowich and Ward [10] is constant, too, for n → ∞ and data with separation distance 1, while the dependence on c is

∥A^{−1}∥_{2,2} ≤ Ce^{4cd}.

To compare with our results, we consider arbitrary distributions of centres in [0, n] and get

∥A^{−1}∥_{2,2} ≥ \frac{π}{8n^{2}} exp (\frac{n − 1}{n} \frac{c}{2}) (1 + 2\frac{c^{2}}{n^{2}} )^{−1/2}  (7.1)

from (5.1) with exponential behaviour for c → ∞ with n fixed. Thus no other choice of centres can get rid of this exponential increase of ∥A^{−1}∥_{2,2} with c.

The variation of (7.1) with n for fixed c is off from Baxter’s optimal bound for equidistant data by only a factor of n^{−2}. Note that a factor n^{−1} may be due to our special combination of matrix norms, and that there must be at least some leeway to optimize placements of centres, which is clearly bounded by gaining a factor of n^{−2}.

We now compare our bounds with those of Ball [1] and Narcowich and Ward [10] [11] [12] and Ball, Sivakumar, and Ward [2] for irregular centres. Since these results are in terms of the separation distance which does not enter explicitly into our results, we have to make sure that the scaling is fair. Thus we can either keep the separation distance fixed and let the centres spread out into all of IR^{d} when their number n tends to infinity, or we can consider large data sets of centres contained in the unit cube [0, 1]^{d} of IR^{d}, letting the separation distance tend to zero when the number of centres tends to infinity. We choose the latter possibility because it is somewhat more related to possible applications. The diameter K of the sets of centres we consider will thus always be bounded by \sqrt{d}, and the Euclidean separation distance h for n centres in [0, 1]^{d} will be at most

h ≤ \sqrt{d} (n^{1/d} − 1)^{−1} ≈ \sqrt{d} n^{−1/d},  (7.2)

as can be easily shown by summing the volumes of sufficiently small disjoint cubes around each centre. In the following we shall simply use h = \sqrt{d} n^{−1/d} because we are mainly interested in the case n → ∞.

For multiquadrics Φ_{c}(r) = (c^{2} + r^{2})^{1/2} the results of Narcowich and Ward [10], as refined in [2] by Ball, Sivakumar, and Ward, yield

∥A^{−1}∥_{2,2} ≤ C(d)\frac{1}{h} exp (2d\frac{c}{h})

with a constant C(d) not depending on c, h, and n. Due to (7.2) this estimate is always worse than

∥A^{−1}∥_{2,2} ≤ C(d) n^{1/d} exp (2cn^{1/d}d)

if arbitrary placements of n centres in [0, 1]^{d} are allowed. Regular distribution of centres on a scaled integer lattice will yield h = n^{−1/d} and

C(d) n^{1/d} exp (\frac{1}{2}cn^{1/d}\sqrt{d}) ≤ ∥A^{−1}∥_{2,2}

according to Ball, Sivakumar, and Ward [2], which has the same asymptotic behaviour for n → ∞ as the best possible upper bound. Thus, as far as upper bounds in terms of the separation distance are concerned, the approach of Narcowich and Ward gives a best possible result for n regular data with n → ∞.

Our approach allows arbitrary centres in [0, 1]^{d} and proves

C_{1}(d)\frac{1}{n\sqrt{d + 2c^{2}}} exp (c\frac{2(d!n/2)^{1/d} + C_{2}(d)}{2\sqrt{d}}) ≤ ∥A^{−1}∥_{2,2}

with suitable constants depending on d only. In the limit n → ∞, the exponential increase cannot be overcome by optimized placement of centres; there is only a factor of at most O(n^{−1−1/d}) to be gained.

The dependence on c is exponential in all cases, and the exponential behaviour is eliminated, if c ≈ n^{−1/d} is chosen. Then, up to constants C_{1}, C_{2} not depending on n, we have

C_{1}n^{−1} ≤ ∥A^{−1}∥_{2,2} ≤ C_{2} n^{1/d}

as the possible variation of ∥A^{−1}∥_{2,2} with the centres, optimizing from a regular distribution. Note that the factor n^{−1} occurring above was introduced by solving the matrix approximation problem of section 2 in the “wrong” norm.

For inverse multiquadrics Φ_{c}(r) = (c^{2} + r^{2})^{−1/2} the lower bound of Ball, Sivakumar, and Ward [2] for regular centres is

C(d) c exp (\frac{1}{2}cn^{1/d}\sqrt{d}) ≤ ∥A^{−1}∥_{2,2},

while we get

C(d)\frac{1}{nc} exp (\frac{c}{2\sqrt{d}} (2 (\frac{d!}{2}n)^{1/d} − 1)) ≤ ∥A^{−1}∥_{2,2}

for arbitrary centres, which is smaller by only a factor O(n^{−1}) for n → ∞.

We now consider thin–plate splines with Φ(r) = r^{2} log r with d = 2, m = 2. Here, the best possible form of the upper bound by Narcowich and Ward [12] is

∥A^{−1}∥_{2,2} ≤ C(d)n^{3},

while we get

C(d)n^{−1/2} ≤ ∥A^{−1}∥_{2,2}.

Numerical experiments indicate that regular data have ∥A^{−1}∥_{2,2} = O(n^{1/2}) for n → ∞. Thus we conjecture that our bound is off by a factor of at most O(n^{−1}) from the actual behaviour of ∥A^{−1}∥_{2,2} for regular data.

For Φ(r) = r^{2β} with β ∈ (0, 1), the best possible upper bound by Narcowich and Ward [11] is

∥A^{−1}∥_{2,2} ≤ C_{1}(d, β)n^{2β/d},

and regularly distributed centres yield

C_{2}(d, β)n^{2β/d} ≤ ∥A^{−1}∥_{2,2},

as was shown by Ball, Sivakumar, and Ward [2]. Our general bound is

C_{3}(d, β)n^{2β/d−1} ≤ ∥A^{−1}∥_{2,2}

for arbitrary centres and β = 1/2, containing a O(n^{−1}) factor again.

8 Conclusions

The results of this paper support the hypothesis that regular placement of centres is a good strategy as far as minimization of the condition number of the matrix A is concerned. The theoretically possible gain by optimization of placement of centres is not more than a factor of O(n^{−2}) or O(n^{−1}) for n → ∞ in all cases, and the proof technique indicates that the factor O(n^{−1}) arises for technical reasons only. Possibly the bounds on the Euclidean norm of A^{−1} in the literature can be generalized to hold also for the norm ∥A^{−1}∥_{1,∞}, and then the factor O(n^{−1}) can be eliminated.

Acknowledgement. The author thanks the two referees for a number of useful suggestions that helped to improve the presentation.

References

[1] Ball, K. M., Invertibility of Euclidean distance matrices and radial basis interpolation, CAT Report 201, Texas A& M University, College Station 1989

[2] Ball, K. M., Sivakumar, N., and J. D. Ward, On the sensitivity of radial basis interpolation with respect to minimal separation distance, CAT Report 240, Texas A& M University, College Station 1989 [3] Buhmann, M.D., Multivariable Interpolation using Radial Basis Functions, Ph.D. disser-

tation, University of Cambridge, 1989.

[4] Baxter, B.J.C., Norm estimates for inverses of Toeplitz distance matrices, preprint, Dept. of Applied Math. and Theor. Physics, Univ. of Cambridge, 1991

[5] 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

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

[7] R. Franke (1982). Scattered data interpolation: tests of some methods. Math. Comp., Vol. 38, pp. 181–200

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

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

[10] Narcowich, F.J., and J.D. Ward, Norms of Inverses and Condition Numbers for Matrices Associated with Scattered Data, J. Approx. Th. 64, 69-94, 1991

[11] Narcowich, F.J., and J.D. Ward, Norm Estimates for the Inverses of a General Class of Scattered Data Interpolation Matrices, Preprint, Texas A & M University, 1991

[12] Narcowich, F.J., and J.D. Ward, Norm Estimates for Inverses of Scattered Data Matrices Associated with Completely Monotonic Radial Functions, Preprint, Texas A & M University, 1991

[13] Natanson, I.P., Constructive Function Theory, Vol. I, Frederick Ungar Publishing Co. New York, 1964

[14] 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

[15] Schoenberg, I. J., Metric spaces and completely monotone functions, Ann. of Math. 39 (1938) 811-841

[16] Schönhage, A., Approximationstheorie, Berlin-New York : Gruyter & Co., 1971

Author’s address:

Robert Schaback, Institut für Numerische und Angewandte Mathematik Universität Göttingen Lotzestr. 16–18 D – 3400 Göttingen