Convergence of Euclidean Radial Basis Approximation on Spheres

J. Levesley

Department of Mathematics and Computer Science University of Leicester Leicester LE1 7RH

              jl1@mcs.le.ac.uk

August 11, 2000

Abstract

In this paper we investigate the convergence of radial basis interpolation when all of the data lie on a sphere. Here, we use strictly conditionally positive definite radial functions defined in the ambient Euclidean space. These are strictly conditionally positive definite of the same order when restricted to the sphere. The analysis of v. Golitschek and Light [9] can then be used to give an error estimate for radial approximation in Euclidean space when the data is restricted to the sphere.

1 Introduction

Let S^{d} be the unit sphere in IR^{d+1}, Π_{n} denote the polynomials of degree n in IR^{d+1} and P_{n} denote the polynomials on S^{d}. We will think of P_{n} as the restrictions of polynomials from Π_{n} to the sphere. Let juxtaposition xx_{i} be used to denote the usual inner product of vectors in IR^{d+1}.

A univariate function ϕ : IR_{+} → IR is said to be strictly conditionally positive definite of order m on IR^{d+1} if, for all N ∈ IN, and points x_{1}, x_{2}, …, x_{N} ∈ IR^{d+1}, whenever

∑_{i=1}^{N} c_{i}π(x_{i}) = 0

for all π ∈ Π_{m−1},

∑_{i,j=1}^{N} c_{i}c_{j}ϕ(∥x_{i} − x_{j}∥) > 0.

Similarly, ψ : [−1, 1] → IR is said to be strictly conditionally positive definite of order m on S^{d} if, for all N ∈ IN, and points x_{1}, x_{2}, …, x_{N} ∈ S^{d}, whenever

∑_{i=1}^{N} c_{i}p(x_{i}) = 0

for all p ∈ P_{m−1},

∑_{i,j=1}^{N} c_{i}c_{j}ψ(x_{i}x_{j}) > 0.

Let us now fix our points x_{1}, x_{2}, …, x_{N} ∈ S^{d}. We assume that we have a function f ∈ C(S^{d}) for which we know the values f_{i} = f(x_{i}), i = 1, …, N, and we wish to reconstruct the function f using interpolation. The usual procedure in this context is to use a strictly conditionally positive definite function ψ of some order, m say, on the sphere, and construct an approximation

r(x) = ∑_{i=1}^{N} α_{i}ψ(xx_{i}) + p(x), x ∈ S^{d},  (1)

where p ∈ P_{m−1}. We then compute the coefficients α_{1}, …, α_{N} and the polynomial p with the interpolation conditions

r(x_{j}) = f(x_{j}), j = 1, …, N,

and

∑_{i=1}^{N} α_{i}q(x_{i}) = 0, for all q ∈ P_{m−1}.

Building on the analysis presented in Levesley et. al. [2], von Golitschek and Light [9] give error estimates of the form

|f(x) − r(x)| ≤ C∥f∥_{ψ}h^{k}, x ∈ s^{d},

where h is a measure of the spacing of the points on the sphere, and k depends on the smoothness of the function ψ. Here, ∥f∥_{ψ} is the so called native space seminorm of f, this native space and hence the target function f, depending intimately on the original choice of ψ.

In this paper we use this error estimate to compute the convergence rate of the following interpolation strategy. Let ϕ : IR → IR be strictly conditionally positive definite of order m in Euclidean space, and approximate f ∈ C(S^{d}) using

s(x) = ∑_{i=1}^{N} α_{i}ϕ(∥x − x_{i}∥) + π(x),  (2)

where ∥ · ∥ is the Euclidean distance and π ∈ Π_{m−1}. We compute the coefficients α_{i}, i = 1, …, N, and the polynomial π with the interpolation conditions

s(x_{j}) = f(x_{j}), j = 1, …, N,

and

∑_{i=1}^{N} α_{i}ν(x_{i}) = 0, for all ν ∈ Π_{m−1}.

In Section 2 we show how each conditionally positive definite function in Euclidean space gives rise to a conditionally positive definite function on the sphere. We will show how to write the interpolation problem in the Euclidean space so that the corresponding interpolation problem on the sphere fits the framework in which the results of [9] may be applied. In Section 4 we give decay rates for Gegenbauer coefficients of the restrictions of the polyharmonic splines ϕ(t) = t^{β} or ϕ(t) = t^{β} log t in order to determine the appropriate native spaces. We then apply the results of [9] to give convergence rates for interpolation of functions defined on the sphere using polyharmonic spline approximation in the ambient space.

2 Recasting the interpolation problem on the sphere

Each conditionally positive definite function ϕ in Euclidean space gives rise to a conditionally positive definite function ψ on the sphere in the following way. Since, for x, x_{i} ∈ S^{d},

∥x − x_{i}∥^{2} = ∥x∥^{2} + ∥x_{i}∥^{2} − 2xx_{i}
= 2 − 2xx_{i},

we see that, for x, x_{i} ∈ S^{d},

ψ(xx_{i}) = ϕ((2 − 2xx_{i})^{1/2}) = ϕ(∥x − x_{i}∥).  (3)

Now, let

∑_{i=1}^{N} c_{i}p(x_{i}) = 0

for all p ∈ P_{m−1}. Then, since P_{n} = Π_{n}|_{S^{d}}, and x_{i} ∈ S^{d}, i = 1, …, N,

∑_{i=1}^{N} c_{i}π(x_{i}) = 0

for all π ∈ Π_{m−1}. As ϕ is strictly conditionally positive definite of order m,

∑_{i,j=1}^{N} c_{i}c_{j}ϕ(∥x_{i} − x_{j}∥) > 0,

and so by (3)

∑_{i,j=1}^{N} c_{i}c_{j}ψ(x_{i}x_{j}) > 0.

Thus we see that ψ is strictly conditionally positive definite of order m on the sphere.

We sum these findings up in the following proposition:

Proposition 1 Let ϕ be strictly conditionally positive definite of order m on IR^{d+1}. Then, ψ(t) = ϕ((2 − 2t)^{1/2}) is strictly conditionally positive definite of order m on S^{d}.

3 Basic theory and convergence rates

Let

[f, g] = ∫_{S^{d}} fgdμ

denote the usual inner product of functions defined on the sphere. Here dμ is used to denote the surface measure on the sphere. Let H_{n} = P_{n} ⋂ P_{n−1}^{⊥} be the space of spherical harmonics polynomials of degree n, where orthogonality is with respect to the above inner product. Suppose that {Y_{n}^{1}, …, Y_{n}^{d_{n}}} is an orthonormal basis for H_{n}. Then every element of L_{1}(S^{d}) has an associated Fourier series expansion

f ∼ ∑_{n=0}^{∞} ∑_{j=1}^{d_{n}} a_{n}^{j}(f)Y_{n}^{j},

where

a_{n}^{j}(f) = ∫_{S^{d}} fY_{n}^{j}dμ.

Let

ψ = ∑_{k=0}^{∞} ψ_{k}P_{k}^{λ},

where P_{k}^{λ} are the Gegenbauer polynomials, orthogonal with respect to the weight (1−t^{2})^{λ−1/2}; see Szegö [8] for details. Since ψ is conditionally positive definite of order m we know that ψ_{k} ≥ 0 for all k ≥ m (see [4]) . We will assume that ψ_{k} > 0 for all such k (this assumption will hold for all restrictions of polyharmonic splines to the sphere) so that ψ is strictly conditionally positive definite of order m. Then, the native space induced by ψ is

H_{ψ} = {f ∈ L_{2}(S^{d}) : ∥f∥_{ψ}^{2} = ∑_{k=m}^{∞} ψ_{k}^{−1} ∑_{j=1}^{d_{k}}(a_{k}^{j}(f))^{2} < ∞} .

We may assume, without loss of generality, that ψ_{k} = 0 for 0 ≤ k < m. For suppose that this were not the case. Let us consider the new basis function

\overline{ψ} = ψ − ∑_{k=0}^{m−1} ψ_{k}P_{k}^{λ}.

Then the interpolation problem for this function would be to find a function of the form

\overline{r}(x) = ∑_{i=1}^{N} \overline{α}_{i}\overline{ψ}(xx_{i}) + \overline{p}(x), x ∈ S^{d},

where \overline{p} ∈ P_{m−1}. We then compute the coefficients \overline{α}_{1}, …, \overline{α}_{N} and the polynomial \overline{p} with the interpolation conditions

\overline{r}(x_{j}) = f(x_{j}), j = 1, …, N,

and

∑_{i=1}^{N} \overline{α}_{i}q(x_{i}) = 0, for all q ∈ P_{m−1}.

By the last set of conditions,

∑_{i=1}^{N} \overline{α}_{i}\overline{ψ}(xx_{i}) = ∑_{i=1}^{N} \overline{α}_{i}ψ(xx_{i}),

so that the conditions satisfied by \overline{α}_{1}, …, \overline{α}_{N} and \overline{p} are the same as those satisfied by α_{1}, …, α_{N} and p. Since the solution of the above interpolation equations is unique \overline{r} = r.

We require that the first m Gegenbauer coefficients of ψ be zero to use Corollary 3.3 in [9], though this result is not stated in the most convenient form for our purposes. A slight alteration in their proof in the style of Theorem 4.4 of [3] gives us

Theorem 1 For each h ∈ (0, 1), let X_{h} = {x_{1}, …, x_{N(h)}} ⊂ S^{d} be such that

max_{x∈S^{d}} min_{y∈X_{h}} cos^{−1}(xy) < h.

Given f ∈ H_{ψ} let r_{f} be the interpolant to f as given by (1) and the following interpolation conditions. Then the are numbers C and h_{0} > 0 such that

∥f − r_{f}∥_{∞,S^{d}} ≤ C\sqrt{inf_{q∈Π_{m−1}^{1}} max_{t∈[1−h^{2},1]} |ψ(t) − q(t)|}∥f∥_{ψ}
= C\sqrt{inf_{q∈Π_{m−1}^{1}} max_{t∈[0,h]} |ϕ(t) − q(t)|}∥f∥_{ψ}

for all h < h_{0}, where Π_{m−1}^{1} denotes the univariate polynomials of degree m − 1.

4 Applications

In this section we will apply the results of the previous section when ψ(t) = (1 − t)^{β} for β > 0 or ϕ(t) = (1 − t)^{β} log(1 − t) for β ∈ IN. The corresponding functions in Euclidean space, t^{2β} and t^{2β} log t, are strictly conditionally negative definite of order [β] + 1 ([·] denotes the integer part) if [β] is even, and are strictly conditionally positive definite of order |β| + 1 if [β] is odd; see Micchelli [5].

In [6] the Gegenbauer coefficients for ψ above are computed. We give these results in the following proposition.

Proposition 2 Let λ ≥ −1/2, and k ∈ IN ∪ {0}.

i) For β ≥ 0,

a_{k}^{β,λ} := ∫_{−1}^{1}(1 − t)^{β}(1 − t^{2})^{λ−1/2}P_{k}^{λ}(t) dt

= (−1)^{k}2^{β+2λ+2k} \frac{Γ(β + λ + 1)Γ(λ + k + \frac{1}{2})Γ(2λ + k)}{k! Γ(2λ + 2k)Γ(β + 2λ + k + 1)}

×\frac{Γ(β + 1)Γ(λ + k)}{Γ(β − k + 1)Γ(λ)}

ii) If β ∈ IN, k ≤ β,

b_{k}^{β,λ} := ∫_{−1}^{1}(1 − t)^{β} log(1 − t)(1 − t^{2})^{λ−1/2}P_{k}^{λ}(t) dt

= (−1)^{k} 2^{β+2λ+2k} \frac{Γ(β + λ + \frac{1}{2}) Γ(λ + k + \frac{1}{2}) Γ(2 λ + k)}{k! Γ(2 λ + 2 k) Γ(β + 2 λ + k + 1)}

×\frac{Γ(β + 1) Γ(λ + k)}{Γ(β − k + 1)Γ(λ)}

×(log 2 + Ψ(β + λ + 1) + Ψ(β + 1) − Ψ(β + 2λ + k + 1)

−ψ(β − k + 1)),

where Ψ = Γ′/Γ is the digamma function; see [1, 6.3].

If k > β,

b_{k}^{β,λ} = (−1)^{β+1} 2^{β+2λ+2k} \frac{Γ(β + λ + \frac{1}{2}) Γ(λ + k + \frac{1}{2})}{k! Γ(2 λ + 2 k)}

×\frac{Γ(2 λ + k)Γ(β + 1) Γ(λ + k)(k − β)!}{Γ(β + 2 λ + k + 1) Γ(λ)}.

The Γ function is positive for positive argument. If x ∈ (−2l + 1, −2l+2) for some l ∈ IN then Γ(x) < 0 whilst if x ∈ (−2l, −2l+1) for some l ∈ IN then Γ(x) > 0. By examining the factor (−1)^{k}Γ(β −k +1) in i) above we see that, for l ∈ IN,

\begin{matrix}a_{k}^{β,λ} < 0, & k > β, & β ∈ (2(l − 1), 2l − 1), \\ a_{k}^{β,λ} > 0, & k > β, & β ∈ (2l − 1, 2l). \end{matrix}

Similarly,

\begin{matrix}b_{k}^{β,λ} < 0, & k > β, & β even, \\ b_{k}^{β,λ} > 0, & k > β, & β odd. \end{matrix}

Thus, (−1)^{[β]+1}ψ is strictly conditionally positive definite of order [β] + 1 on the sphere, which, in conjunction with the remarks at the start of this section, confirms the conclusions of Proposition 1.

Using the asymptotic estimate

Γ(az + b) ∼ \sqrt{2π}e^{−az}(az)^{az+b+1/2}

(see [1, 6.1.39]) we see that

a_{k}^{β,λ}

∼ 2^{2k} \frac{e^{−k}k^{λ+k}e^{−k}k^{2λ+k−1/2}e^{−k}k^{λ+k−1/2}}{e^{−k}k^{k+1/2}e^{−2k}(2k)^{2λ+2k−1/2}e^{−k}k^{β+2λ+k+1/2}e^{k}k^{β−k+1/2}}
∼ k^{−2(β+1)}.

Using the same asymptotic expression we see that

b_{k}^{β,λ} ∼ k^{−2(β+1)}.

Hence, in either case, the Hilbert space of functions we approximate is

H_{ψ} = {f : ∑_{k=m}^{∞} k^{2(β+1)} ∑_{j=1}^{d_{k}}(a_{k}^{j}(f))^{2} < ∞} .

It is known (see e.g. Schaback [7]) that, for ϕ(t) = t^{2β}, β > 0, or ϕ(t) = t^{2β} log t with β ∈ IN, and m = [β] + 1,

inf_{q∈Π_{m−1}^{1}} max_{t∈[0,h]} |ϕ(t) − q(t)| ≤ Ch^{β},

for some constant C independent of h.

Combining this with Theorem 1 we have the following result for polyharmonic spline interpolation on the sphere:

Theorem 2 Let ϕ(t) = t^{2β}, β > 0, or ϕ(t) = t^{2β} log t for β ∈ IN. Let m = [β] + 1. For each h ∈ (0, 1), let X_{h} = {x_{1}, …, x_{N(h)}} ⊂ S^{d} be such that

max_{x∈S^{d}} min_{y∈X_{h}} cos^{−1}(xy) < h.

Given f ∈ H_{ψ} let s_{f} be the interpolant to f as given by (2) and the following interpolation conditions. Then the are numbers C and h_{0} > 0 such that

∥f − s_{f}∥_{∞,S^{d}} ≤ Ch^{β}∥f∥_{ψ}.

Remarks 1 It is reassuring that we recover the same order of convergence as for interpolation using, for example thin plate splines (ϕ(t) = t^{2} log t) in the plane, where we expect to get O(h) convergence for h-separated data.

References

[1] M. Abramowitz and I. A. Stegun. Handbook of Mathematical Functions, National Bureau of Standards. Dover Publications, 1964.

[2] J. Levesley, W. A. Light, D. L. Ragozin, and X. Sun. A simple approach to the variational theory for interpolation on spheres. In New Developments in Approximation Theory (Dortmund 1998), volume 132 of Inter. Ser. Numer. Math. Birkhauser, Basel, 1999.

[3] W. A. Light and H. Wayne. Spaces of distributions, interpolation by translates of a basis function and error estimates. Numer. Math. 81 (1999), 415–450.

[4] Z. Luo and J. Levesley Norm estimates for inverses of matrices arising from interpolation by strictly conditional positive definite kernels on spheres. Numer. Func. Anal. and Opt. 19 (1998), 849– 871.

[5] C. Micchelli. Interpolation of scattered data: Distance matrices and conditionally positive definite functions. Constr. Approx. 2 (1986), 11–22.

[6] C. Odell and J. Levesley. Evaluation of some integrals arising in approximation on the sphere. Preprint, 2000.

[7] R. Schaback. Reconstruction of Multivariate Functions from Scattered Data. University of Göttingen, 1997.

[8] G. Szegö. Orthogonal Polynomials, volume XXIII of American Mathematical Society Colloquium Publications. American Mathematical Society, 1959.

[9] M. von Golitschek and W. A. Light. Interpolation by polynomials and radial basis functions on spheres. Constr. Approx. to appear, 2000.