Hermite Interpolation

with Radial Basis Functions on Spheres

Gregory E. Fasshauer

Abstract. We show how conditionally negative definite functions on spheres coupled with strictly completely monotone functions (or functions whose derivative is strictly completely monotone) can be used for Hermite interpolation. The classes of functions thus obtained have the advantage over the strictly positive definite functions studied in [17] that closed form representations (as opposed to series expansions) are readily available. Furthermore, our functions include the historically significant spherical multiquadrics. Numerical results are also presented.

AMS classification: 41A05, 41A63, 42A82. Key words and phrases: Spherical interpolation, Hermite interpolation, Radial basis functions.

1. Introduction

In 1975 R. Hardy mentioned the possibility of using multiquadric basis functions for Hermite interpolation (see [10], or the survey paper [11]). This problem, however, was not further investigated until the paper [29] by Wu appeared. Since then, the interest in this topic seems to have increased significantly. In [29] Wu deals with Hermite-Birkhoff interpolation in IR^{m}. His results are limited in the sense that one can have only one interpolation condition per data point (i.e., some linear combination of function value and derivatives). Sun [25] eliminates this restriction. He also deals with the Euclidean setting and gives results analogous to those of Micchelli [15] for ordinary interpolation. Narcowich and Ward [18] and Narcowich [17] develop an even more general theory of Hermite interpolation for conditionally positive definite (matrix-valued) functions in IR^{m}, and positive definite kernels on compact Riemannian manifolds (the sphere being a special case), respectively. Freeden [8, 9] studied a similar problem on the unit sphere in IR^{3} employing a slightly different type of linear functionals.

As in the Euclidean case, there are functions on spheres which are not positive definite, but only conditionally positive (negative) definite. These functions were characterized by Menegatto [13]. It is the goal of this paper to derive sufficient conditions for which Hermite interpolation based on these functions is well posed. We will see below that using this approach it is possible to obtain many closed form representations of basis functions which can be used in practice.

We introduce the interpolation problem and state our main theorems in Section 2. Section 3 contains some basic definitions needed in the further development. The central part of the paper is Section 4 which contains the proofs of our main results. In the last section we examine two first-order Hermite interpolation problems on the sphere.

2. Formulation of the Problem and Main Results

We will formulate all our results for S^{m}, the unit sphere in IR^{m+1}. However, we mention that an extension to S^{∞}, the unit sphere in ℓ_{2}, is possible using the functions cos^{k}(t) in expansion (2.3) below instead of the normalized Gegenbauer polynomials p_{k}^{λ}(cos t) (see e.g. [14, 23]).

We are interested in

Problem 2.1. Suppose we are given data d_{i} = L_{i}f, i = 1, …, n, generated by a function f ∈ C^{r}(S^{m}) and some linearly independent set L = {L_{1}, …, L_{n}} of linear functionals which are either bounded in the supremum norm on C^{r}(S^{m}) or can be approximated by such functionals. We want to find an interpolant of the form

s(x) = ∑_{j=1}^{n}c_{j}L_{j}^{y}g(δ(x, y)),  x, y ∈ S^{m},  (2.1)

satisfying

L_{i}s = d_{i}, i = 1, …, n.

Here g is a real-valued function in C^{r}[0, π], and δ(x, y) is the geodesic distance between two points x and y on the unit sphere S^{m}, i.e.,

δ(x, y) = arccos(x · y),

where x · y stands for the dot product of x and y in IR^{m+1}. In analogy to the Euclidean case we will frequently refer to g as a “radial function”. This is motivated by the fact that we only consider g in its composition with the geodesic distance, say G(ξ) = g(δ(x, ξ)), and thus for any fixed x ∈ S^{m} we have

δ(x, y) = δ(x, z) ⟹ G(y) = G(z), for all y, z ∈ S^{m}.

We write L^{y} in (2.1) to indicate that the linear functional acts on a function of the variable y. It should also be pointed out that the linear independence of the L_{i} is equivalent to the absence of redundancies in the given data.

The interpolation problem above will have a unique solution if the interpolation matrix A whose entries are given by

a_{ij} = L_{i}^{x}L_{j}^{y}g(δ(x, y)), i, j = 1, …, n,  (2.2)

is nonsingular.

Remark. For standard Hermite interpolation we can think of the L_{i} either as point evaluation functionals (yielding Lagrange interpolation conditions), or as functionals evaluating certain derivatives of a given function at specified points on the unit sphere S^{m}. However, we do allow the set L to contain more general functionals such as, e.g., local integrals. Furthermore, we stress that there is no assumption which requires that – if interpolation to derivative information is being performed – these derivatives must be in consecutive order.

For the case in which g ∘ δ is a C^{∞} strictly positive definite kernel on S^{m} × S^{m} Narcowich [17] showed in certain cases the nonsingularity of A by an argument involving the theory of distributions. We will be able to avoid the use of this sophisticated machinery and establish our proofs much in the flavor of those in [15]. The precise definition of strict complete monotonicity used in the following two theorems will be given in the next section.

Theorem 2.2. Let

h(t) = a_{0} + ∑_{k=1}^{∞}a_{k}[1 − p_{k}^{λ}(cos t)],  (2.3)

where λ = (m − 1)/2, all a_{k} ≥ 0 (and a_{k} > 0 for some odd k), ∑_{k=1}^{∞}a_{k} < ∞, p_{k}^{λ}(x) = P_{k}^{λ}(x)/P_{k}^{λ}(1), and P_{k}^{λ}(x) are the usual Gegenbauer polynomials (see e.g. [26]).

If H is a continuous and nonnegative function on [0, ∞), and H′ is strictly completely monotone on [0, ∞), then the matrix A in (2.2) formed by using the function g with g(t) = H(h(t)) is symmetric. Moreover,

(i) if L contains at least two point evaluation functionals, then A has one positive and n − 1 negative eigenvalues.

(ii) if a_{0} > 0 and L contains at least one point evaluation functional, then A has

one positive and n − 1 negative eigenvalues. (iii) if L consists only of functionals which annihilate constants, then A is negative

definite.

Theorem 2.3. If h is as in (2.3), and if H is strictly completely monotone on [0, ∞), then the matrix A formed by using g(t) = H(h(t)) in (2.2) is symmetric and positive definite.

Remark 1. Since the function H with H(t) = t^{α}, 0 < α < 1, has a strictly completely monotone derivative, this is a special case of Theorem 2.2. We point out that in the historic treatment of isometric imbeddings (Blumenthal and Schoenberg first studied imbeddings with metric transforms of the form t^{α}, see e.g. [21]), as well as in the discussion in [5] this function plays an important role.

Remark 2. Note that Theorem 2.2 does not cover the case where the generating set L consists only of two nonempty subsets; one containing derivative functionals, and the other one containing functionals which do not include point evaluation functionals.

3. Background

In the formulation of Theorems 2.2 and 2.3 we used the notion of a strictly completely monotone function.

Definition 3.1. A function F : [0, ∞) → IR which is in C[0, ∞) ∩ C^{∞}(0, ∞) and which satisfies

(−1)^{ℓ}F^{(ℓ)}(t) ≥ 0, ℓ = 0, 1, 2, …,

is called completely monotone on [0, ∞). If the inequality is strict, then F is strictly completely monotone.

The following characterization of a strictly completely monotone function in terms of Laplace-Stieltjes transforms can be found in [25], and is a modification of the well-known Bernstein-Widder theorem (see [28]).

Theorem 3.2. A function F : [0, ∞) → IR is strictly completely monotone on [0, ∞) if and only if F is of the form

F (t) = ∫_{0}^{∞} e^{−σt}dμ(σ),

where μ is a finite nonnegative Borel measure on [0, ∞) which satisfies

∫_{0}^{∞} dμ(σ) < ∞  and  ∫_{ϵ}^{∞} dμ(σ) > 0

for some positive constant ϵ.

Note that a strictly completely monotone function can be characterized as a completely monotone function which is not constant.

Schoenberg [22] gave the following characterization of a function whose derivative is strictly completely monotone.

Theorem 3.3. In order that F ∈ C[0, ∞) and F′ is strictly completely monotone on [0, ∞) it is necessary and sufficient that

F (t) − F (0) = ∫_{0}^{∞}[1 − e^{−σt}]\frac{dμ(σ)}{σ},

where μ is a finite nonnegative Borel measure on [0, ∞) which satisfies

∫_{1}^{∞} \frac{dμ(σ)}{σ} < ∞  and  ∫_{ϵ}^{∞} dμ(σ) > 0

for some positive constant ϵ.

In the discussion that follows we will not require the environment of a normed linear space. Usually it will be sufficient to work with the following “relaxed” metric spaces defined e.g. in [22].

Definition 3.4. A set S is called a semi-metric space, if S is provided with a distance function d : S × S → [0, ∞), such that for arbitrary points P, Q ∈ S

(i) d(P, Q) = d(Q, P), (ii) d(P, Q) = 0 ⟺ P = Q.

The following property (again due to Schoenberg) plays one of the central roles in our discussion and is a generalization of the classical notion of positive definite functions given by Mathias and Bochner (see e.g. [24]).

Definition 3.5. Let S be a semi-metric space with distance function d. A function F ∈ C[0, ∞) is called positive definite in S, if

∑_{i=1}^{n} ∑_{j=1}^{n}F (d(x_{i}, x_{j}))c_{i}c_{j} ≥ 0  (3.1)

for any n distinct points x_{i} ∈ S, i = 1, …, n, arbitrary real constants c_{i}, and all n = 2, 3, . . .. The function F is strictly positive definite if equality in (3.1) implies c_{1} = ⋯ = c_{n} = 0.

Remark 3. This distinction between strictness and nonstrictness is the essential modification of these classical ideas for their application in approximation theory.

4. Proofs

The key ingredients in the proofs of Theorems 2.2 and 2.3 are the following two lemmas.

Lemma 4.1. The function g given by g(t) = e^{−σh(t)}, with h as in (2.3) and σ > 0, is strictly positive definite on S^{m}.

Proof: This proof is based essentially on Schoenberg’s discussion of isometric imbeddings of (semi-)metric spaces into the Hilbert space ℓ_{2} (for a discussion of isometric imbeddings see e.g., [27]). He considered one semi-metric space S with distance function d and a new semi-metric space F (S) based on the same set S, but with distance function F ∘d. Schoenberg calls the function F the metric transform. For isometric imbeddings of S into ℓ_{2} he gave the following characterization (see [22]), where the simplicial class Π(S) denotes the class of all functions F for which F (S) is isometrically imbeddable in ℓ_{2}:

A function F (t) (F (0) = 0) belongs to the simplicial class Π(S) of a semimetric space S if and only if e^{−σF^{2}(t)} is positive definite in S for all σ > 0. In order to facilitate the following arguments, we let f(t) = h(t)−a_{0}. In [23] the

isometric imbeddings of S^{m} into ℓ_{2} are characterized as having metric transforms of the form F^{2}(t) = f(t) with a_{k} ≥ 0, k = 1, 2, … (we will make use of the fact that, for our purposes, at least one a_{k} with odd index must be positive below). Thus, we can immediately conclude that e^{−σf}, and therefore also e^{−σh}, is positive definite in S^{m}. It remains to establish strictness.

From Schoenberg’s work we know that the functions F with F^{2}(t) = f(t) are metric transforms associated with imbeddings of S^{m} into ℓ_{2}. Therefore, by the geometric interpretation of an isometric imbedding given by Menger and Schoenberg (see [20, 22]), we know that corresponding to the (distinct) data points x_{i}, i = 1, …, n, on S^{m} there exist (distinct) points y_{i} in Euclidean space IR^{N} (N = n − 1) such that their Euclidean distances equal the geodesic distances of the x’s (the points y_{i} are the vertices of simplices which motivates the terminology simplicial class above.). We use this characterization to interpret the values f(δ(x_{i}, x_{j})) as Euclidean distances ∥y_{i} − y_{j}∥^{2}. We can do this if we assume that at least one a_{k} with odd index in (2.3) is positive. This assumption on the a_{k} ensures that f(t) = 0 does not have two roots t = 0, and t = π violating property (ii) of Definition 3.4. Thus, for the rest of the proof we can assume that we are dealing with Euclidean distances ∥y_{i} − y_{j}∥^{2}.

It is now a standard argument in the radial basis function literature dealing with interpolation in Euclidean spaces (see e.g. [15]) to consider the Fourier transform

e^{−∥τ∥^{2}} = \frac{1}{(\sqrt{π})^{N}} ∫_{IR_{N}} e^{2\sqrt{−1}(τ·θ)}e^{−∥θ∥^{2}}dθ,

where τ and θ are vectors in IR^{N}. Then one sets τ = y_{i} − y_{j}, and shows that after multiplication by c_{i}c_{j} and summation over i and j it follows that c_{1} = ⋯ c_{n} = 0.

Clearly, our arguments also hold for

e^{−σh(δ(x_{i},x_{j}))} = e^{−σa_{0}}e^{−σf(δ(x_{i},x_{j}))} = e^{−σa_{0}}e^{−σ∥y_{i}−y_{j}∥^{2}}

with σ > 0, and thus the statement of the lemma follows.

Remark 1. We point out that the statement of Lemma 4.1 does not contradict the results of Narcowich (see e.g. [17, Theorems 3.6 and 4.2]) since Narcowich considers a notion of positive definiteness in the sense of distributions which is stronger than that of Definition 3.5.

Lemma 4.2. If the function g is of the form g(t) = e^{−σh(t)}, σ > 0, with h as in (2.3), then the matrix A defined in (2.2) is positive definite.

Proof: By Lemma 4.1 we know that g is (strictly) positive definite on S^{m}. Therefore g can be expressed as (see [17, 23])

g(t) = ∑_{k=0}^{∞} b_{k}P_{k,m+1}(x · y),  b_{k} ≥ 0,  (4.1)

where P_{k,m+1} are Legendre polynomials of degree k in m + 1 dimensions, and cos(t) = cos(arccos(x · y)) = x · y. It is known that the series above converges absolutely and uniformly for t ∈ [0, π] (see [23]). We use Legendre polynomials here instead of Gegenbauer polynomials, since they are the same up to a constant factor, and are somewhat easier to handle.

We split the proof into two parts. First we show that the matrix A is positive definite if the functionals in L are bounded in the supremum norm on C(S^{m}). Then we establish the validity of the lemma for such functionals which can be approximated by continuous linear functionals (this covers in particular those functionals which correspond to evaluation of derivatives).

If we are dealing with continuous linear functionals on (C(S^{m}), ∥ · ∥_{∞}) we can employ the Riesz representation theorem to regard the linear functionals in L as integrals of the form

L_{i}f = ∫_{S^{m}} f(x)dμ_{i}(x), i = 1, …, n,

where the μ_{i} are regular signed Borel measures on S^{m}. With this and (4.1) we can rewrite the entries of the matrix A as

L_{i}^{x}L_{j}^{y}e^{−σh(δ(x,y))} = ∫_{S^{m}} ∫_{S^{m}} ∑_{k=0}^{∞} b_{k}P_{k,m+1}(x · y)dμ_{i}(x)dμ_{j}(y), i, j = 1, …, n.

(4.2)

The other ingredient we use is the addition theorem for spherical harmonics which states that

∑_{l=1}^{N(m,k)} Y_{k,l}(x)Y_{k,l}(y) = \frac{N(m, k)}{ω_{m}}P_{k,m+1}(x · y),  (4.3)

where Y_{k,l}, l = 1, …, N(m, k), are the spherical harmonics of order k on S^{m}, N(m, k) = \frac{(2k+m−1)Γ(k+m−1)}{Γ(k+1)Γ(m)} is the dimension of the space of these functions, and

ω_{m} = \frac{2π^{(m+1)/2}}{Γ((m+1)/2)} is the surface area of S^{m} (see e.g. [16]).

Replacing P_{k,m+1} in (4.2) via (4.3) we arrive at

L_{i}^{x}L_{j}^{y}e^{−σh(δ(x,y))} = ∫_{S^{m}} ∫_{S^{m}} ∑_{k=0}^{∞} b_{k} \frac{ω_{m}}{N(m, k)} ∑_{l=1}^{N(m,k)} Y_{k,l}(x)Y_{k,l}(y)dμ_{i}(x)dμ_{j}(y).

Since the infinite series is absolutely and uniformly convergent as noted earlier we can interchange the infinite summation with the integration which yields

\begin{aligned}L_{i}^{x}L_{j}^{y}e^{−σh(δ(x,y))} & = ∑_{k=0}^{∞} b_{k} \frac{ω_{m}}{N(m, k)} ∫_{S^{m}} ∫_{S^{m}} ∑_{l=1}^{N(m,k)} Y_{k,l}(x)Y_{k,l}(y)dμ_{i}(x)dμ_{j}(y) \\ & = ∑_{k=0}^{∞} b_{k} \frac{ω_{m}}{N(m, k)} ∑_{l=1}^{N(m,k)} ∫_{S^{m}} Y_{k,l}(x)dμ_{i}(x) ∫_{S^{m}} Y_{k,l}(y)dμ_{j}(y). \end{aligned}

Now, multiplying by c_{i}c_{j} and summing over i and j, we obtain

c^{T} Ac = ∑_{i=1}^{n} ∑_{j=1}^{n} c_{i}c_{j} ∑_{k=0}^{∞} b_{k} \frac{ω_{m}}{N(m, k)} ∑_{l=1}^{N(m,k)} ∫_{S^{m}} Y_{k,l}(x)dμ_{i}(x) ∫_{S^{m}} Y_{k,l}(y)dμ_{j}(y)

= ∑_{k=0}^{∞} b_{k} \frac{ω_{m}}{N(m, k)} ∑_{l=1}^{N(m,k)} ∑_{i=1}^{n} ∑_{j=1}^{n} c_{i}c_{j} ∫_{S^{m}} Y_{k,l}(x)dμ_{i}(x) ∫_{S^{m}} Y_{k,l}(y)dμ_{j}(y)

= ∑_{k=0}^{∞} b_{k} \frac{ω_{m}}{N(m, k)} ∑_{l=1}^{N(m,k)} [∑_{i=1}^{n} c_{i} ∫_{S^{m}} Y_{k,l}(x)dμ_{i}(x)]^{2} ,

which is obviously nonnegative. Since certainly not all of the b_{k} are zero, we see that the only possibility for c^{T} Ac = 0 is

∑_{i=1}^{n} c_{i} ∫_{S^{m}} Y_{k,l}(x)dμ_{i}(x) = 0

for l = 1, …, N(m, k), and all k such that b_{k} > 0. Remembering that this is the same as

∑_{i=1}^{n} c_{i}L_{i}Y_{k,l} = 0,

and making use of the fact that L is a linearly independent set, it follows that c_{1} = ⋯ = c_{n} = 0. This establishes the claim of the lemma for the case of continuous linear functionals.

We illustrate the second part of the proof by assuming all functionals in L to be derivative functionals which can be written as L_{i} = lim_{ϵ→0} L_{i}^{ϵ}. The L_{i}^{ϵ} are defined by considering a regular parametric curve C in the neighborhood N_{i} ∩ S^{m} of x_{i}m, parametrized such that C(0) = x_{i}. The approximating functionals are then defined as

L_{i}^{ϵ}f = \frac{f(x_{i} + ϵ) − f(x_{i})}{ϵ},

where x_{i} + ϵ is used as symbolic notation for the point C(ϵ) in the neighborhood of x_{i}. The “direction” of the derivative is thus defined by the choice of C. In the last section of the paper a special choice of C (leading to partial derivatives with respect to the spherical coordinates on S^{2}) is employed.

For any fixed value of ϵ the same argumentation as above will yield

c^{T} A^{ϵ}c = ∑_{k=0}^{∞} b_{k} \frac{ω_{m}}{N(m, k)} ∑_{l=1}^{N(m,k)} [∑_{i=1}^{n} c_{i}L_{i}^{ϵ}Y_{k,l}]^{2} ≥ 0.

Here A^{ϵ} denotes the matrix corresponding to the use of the linear functionals L_{i}^{ϵ}. Now it is also clear that

c^{T} Ac = lim_{ϵ→0} c^{T} A^{ϵ}c ≥ 0.

As mentioned above, the infinite series (4.1) converges absolutely and uniformly for t ∈ [0, π], and therefore we can interchange the summation and the limit and arrive at

c^{T} Ac = ∑_{k=0}^{∞} b_{k} \frac{ω_{m}}{N(m, k)} ∑_{l=1}^{N(m,k)} [∑_{i=1}^{n} c_{i} lim_{ϵ→0} L_{i}^{ϵ}Y_{k,l}]^{2} ≥ 0.  (4.4)

Finally, as in the first part of the proof, it follows that the quadratic form (4.4) is positive since the functionals L_{i} = lim_{ϵ→0} L_{i}^{ϵ} were assumed to form a linearly independent set. Any other derivative functionals can be approximated in an analogous way, and thus the proof is complete.

Remark 2. Lemma 4.2 by itself gives a class of radial basis functions suitable for Hermite interpolation on S^{m}. In analogy to the Euclidean setting we refer to the functions e^{−σh} as spherical Gaussians. We point out that this allows for closed form representations of basis functions (since only one coefficient in the series (2.3) needs to be positive) which can be manipulated much more easily in practice than the long series expansions required for strictly positive definite functions. The same is true for the basis functions covered by Theorems 2.2 and 2.3.

Proof of Theorem 2.2. We will split the proof into three parts. First we show that the interpolation matrix A is symmetric, then we prove it is negative definite on a subspace of dimension n−1, and finally we give an argument for the eigenvalue structure of A and thus its nonsingularity.

a) Note that g is a continuous function of the geodesic distance, which itself is a function of the dot product, i.e., G(x, y) = g̃(x·y) = g(arccos(x·y)) is continuous on S^{m} ×S^{m}. Furthermore, by using the same approximation argument combined with the Riesz representation theorem as in the proof of Lemma 4.2 along with a version of Fubini’s theorem for finite measures (see e.g. [4, p. 193]) we can interchange the order in which we apply L^{x} and L^{y}, i.e.,

a_{ij} = L_{i}^{x}L_{j}^{y}g̃(x · y) = L_{j}^{y}L_{i}^{x}g̃(x · y).

Finally, we switch x and y, so that by the symmetry of the dot product we have

L_{j}^{y}L_{i}^{x}g̃(x · y) = L_{j}^{x}L_{i}^{y}g̃(y · x) = L_{j}^{x}L_{i}^{y}g̃(x · y) = a_{ji}.

b) We will now show that A is negative definite whenever

∑_{i=1}^{n} c_{i}L_{i}[1] = 0,  (4.5)

where c(≠ 0) ∈ IR^{n}.

If we apply Theorem 3.3 to the function H and substitute h(δ(x, y)) for t in the integral representation given there, the associated quadratic form is

\begin{aligned}∑_{i=1}^{n} ∑_{j=1}^{n}c_{i}c_{j}L_{i}^{x}L_{j}^{y}H(h(δ(x, y))) & − ∑_{i=1}^{n} ∑_{j=1}^{n}c_{i}c_{j}L_{i}^{x}L_{j}^{y}H(0) = \\ & ∫_{0}^{∞} ∑_{i=1}^{n} ∑_{j=1}^{n}c_{i}c_{j}L_{i}^{x}L_{j}^{y}[1 − e^{−σh(δ(x,y))}]\frac{dμ(σ)}{σ}. \end{aligned}

Here the interchange of the functionals with the integral is permitted by the ideas used in the proof of Lemma 4.2 (Riesz representation theorem, and approximation), Fubini’s theorem for finite measures, and the condition on the measure μ stated in Theorem 3.3 which guarantees the convergence of the integral.

We can now use condition (4.5) to get

∑_{i=1}^{n} ∑_{j=1}^{n}c_{i}c_{j}L_{i}^{x}L_{j}^{y}H(h(δ(x, y))) = − ∫_{0}^{∞} ∑_{i=1}^{n}∑_{j=1}^{n}c_{i}c_{j}L_{i}^{x}L_{j}^{y}e^{−σh(δ(x,y))} \frac{dμ(σ)}{σ}. (4.6)

Finally, Lemma 4.2 tells us that the integrand is positive, and the only possibility for (4.6) to be zero is that the measure μ has all its mass concentrated at the origin, which would contradict the assumption of strict complete monotonicity of the function H′.

c) To complete the proof we need to discuss four cases:

(i) L contains only point evaluation functionals, (ii) L contains at least two point evaluation functionals, (iii) L contains one point evaluation functional, and a_{0} > 0, (iv) L contains only functionals of the form L_{i}[1] = 0.

(i) In this case condition (4.5) simply becomes ∑_{i=1}^{n}c_{i} = 0, and thus A is conditionally (or almost) negative definite in the classical sense (see e.g. [3]). To determine the eigenvalue structure of A we recall that a matrix which is negative definite on a subspace of dimension n − 1 has at least n − 1 negative eigenvalues by the Courant-Fischer theorem (see e.g. [12]). If n ≥ 2 we can show that A must have also at least one positive eigenvalue by concluding from the definition of g that trace(A) ≥ 0. For the special case n = 1 it is necessary to have a_{0} > 0. Then g(0) > 0 and therefore A is nonsingular for n = 1 also (the only entry of A is positive).

(ii) From part b) of the proof we already know that A is negative definite on a subspace of dimension n−1 and therefore has at least n−1 negative eigenvalues. To show that A must have at least one positive eigenvalue in this case, we reorganize the matrix A into the following symmetric block form

A = ( \begin{matrix}A_{11} & A_{12} \\ A_{12}^{T} & A_{22} \end{matrix} ) ,  (4.7)

where A_{11} is a block (of size at least 2×2) which corresponds to Lagrange interpolation. Using the symmetry of A and the fact that A_{11} is invertible as shown in (i), we can perform a congruence transformation with the transformation matrix

P = ( \begin{matrix}I & −A_{11}^{−1}A_{12} \\ 0 & I \end{matrix} ) ,

and we then obtain

Ã = P^{T} AP = ( \begin{matrix}A_{11} & 0 \\ 0 & A_{22} − A_{12}^{T}A_{11}^{−1}A_{12} \end{matrix} ) .

Since it is well known from linear algebra (see e.g. [12]) that congruent symmetric matrices have the same inertia (an integer triple In(M) = (p, n, z) representing the number of positive, negative and zero eigenvalues of the symmetric matrix M), and it is also easily seen that the inertia of a block-diagonal matrix is the sum of the inertias of those diagonal blocks, it follows that In(Ã) ≥ In(A_{11}). Finally, since we just showed in (i) that A_{11} must have at least one positive eigenvalue, the same must be true for Ã, and by congruence for the Hermite matrix A.

(iii) The case n = 1 and a_{0} > 0 will work analogously based on the argumentation

at the end of (i). (iv) If all functionals in L annihilate constants, i.e., L_{i}[1] = 0, for all i = 1, …, n,

then condition (4.5) is void and it follows that A is negative definite.

Remark 3. If the elements of the interpolation matrix A are arranged as in (4.7), then efficient block decomposition algorithms can be used to invert A. In particular, the inertia argument above shows that the Schur complement, A_{22} − A_{12}^{T}A_{11}^{−1}A_{12}, is not only symmetric, but also negative definite.

Proof of Theorem 2.3. The symmetry of the interpolation matrix follows as in the proof of Theorem 2.2. We will now show that A is positive definite. Recall that the radial functions g are of the form g(t) = H(h(t)), where H is strictly completely monotone on [0, ∞). Consider the quadratic form

c^{T} Ac = ∑_{i=1}^{n} ∑_{j=1}^{n}c_{i}c_{j}L_{i}^{x}L_{j}^{y}H(h(δ(x, y)).

Using the Bernstein-Widder characterization of completely monotone functions (see Theorem 3.2), this becomes

c^{T} Ac = ∫_{0}^{∞} ∑_{i=1}^{n}∑_{j=1}^{n}c_{i}c_{j}L_{i}^{x}L_{j}^{y}e^{−σh(δ(x,y))}dμ(σ),

where an argument analogous to that used in the proof of Theorem 2.2 permits the change of the L’s with the integral. Now it follows directly from Lemma 4.2 and the fact that the strict complete monotonicity of H implies that the measure μ does not have all of its mass at the origin that A is positive definite.

5. An Example

In this section we discuss an example of first-order Hermite interpolation on the sphere S^{2}. We will consider the following situation:

Example 5.1. Assume we are given ν distinct points x_{i}, i = 1, …, ν, on the unit sphere S^{2} in IR^{3} represented in spherical coordinates

x = cos θ cos ϕ
y = sin θ cos ϕ
z = sin ϕ

with θ ∈ [0, 2π) and ϕ ∈ [−π/2, π/2], and a function f ∈ C^{1}(S^{2}), as well as associated data of the form

d_{i} = \{ \begin{matrix}f(x_{i}) & i = 1, …, ν, \\ \frac{∂f}{∂ϕ}(x_{i−ν}) & i = ν + 1, …, 2ν, \\ \frac{∂f}{∂θ} (x_{i−2ν}) & i = 2ν + 1, …, 3ν. \end{matrix}

To be specific, we could imagine using spherical multiquadrics (see e.g. [11]), i.e., basis functions of the form

g(t) = \sqrt{1 + R^{2} − 2R cos t},

where R(≠ 1) is a user-selectable parameter, to interpolate to this set of data.

It can easily be seen that the spherical multiquadrics satisfy the conditions of Theorem 2.2 if we let H(t) = t^{1/2}, a_{0} = (1 − R)^{2}, and a_{1} = 2R. Note that the condition R ≠ 1 (or equivalently a_{0} ≠ 0) ensures differentiability of g. The formulation of Example 5.1 means that in Problem 2.1 we would have n = 3ν functionals L_{i}, the first ν corresponding to point evaluations at the data points, and the others corresponding to evaluation of the partial derivatives with respect to the spherical coordinates at the same points.

Proposition 5.2. If the data points x_{i}, i = 1, …, ν, contain either the North or the South Pole, the interpolation matrix A for Example 5.1 is singular.

Proof: We can interpret the partial derivative \frac{∂f}{∂θ} as an inner product of the form ∇f · ϑ where ∇f is the gradient of f (in Cartesian coordinates), and ϑ = (− sin θ cos ϕ, cos θ cos ϕ, 0)^{T}. At the North Pole we have ϕ = π/2 and thus ϑ is the zero vector. Therefore, it is easily seen that a complete row (column) in the matrix A will be zero. The same holds for the South Pole with ϕ = −π/2.

Remark 1. It is immediate to see that the linear functionals chosen in the example above are linearly dependent if one of the poles is among the data points, since the functional corresponding to \frac{∂f}{∂θ} is the zero functional at the poles.

Remark 2. It is clear that the real trouble here is caused by the singularity in the parametrization chosen for the sphere. Thus, the above difficulties could be avoided using a combination of parametrizations which are regular at the poles. Another possibility would be the use of covariant derivatives.

Remark 3. Even though the interpolation matrix is singular for the situation described in Proposition 5.2, one can still solve the interpolation problem since the singularity of the matrix is consistent with zero data on the right-hand side (see the numerical experiment at the end of the paper).

If this approach does not seem appealing, it is possible to give a positive result concerning a slightly modified version of Example 5.1.

Proposition 5.3. If the data points x_{i}, i = 1, …, ν, are distinct, and if α_{i} and β_{i} are two linearly independent directions in the tangent plane to the sphere at x_{i}, then interpolation to data of the form

d_{i} = \{ \begin{matrix}f(x_{i}) & i = 1, …, ν, \\ D_{α_{i−ν}} f(x_{i−ν}) & i = ν + 1, …, 2ν, \\ D_{β_{i−2ν}} f(x_{i−2ν}) & i = 2ν + 1, …, 3ν, \end{matrix}

using any of the functions g classified in Theorems 2.2 and 2.3 is well posed. Here D_{α_{i}} and D_{β_{i}} are the usual directional derivatives in IR^{3}.

Proof: The hypotheses imply the linear independence of the functionals in L.

Remark 4. If we compare

φ = (− cos θ sin ϕ, − sin θ sin ϕ, cos ϕ)^{T} , ϑ = (− sin θ cos ϕ, cos θ cos ϕ, 0)^{T}

with

α_{i} = (− cos θ_{i} sin ϕ_{i}, − sin θ_{i} sin ϕ_{i}, cos ϕ_{i})^{T} , β_{i} = (− sin θ_{i}, cos θ_{i}, 0)^{T} ,

we see that Proposition 5.3 can be interpreted as a “scaled” version of Example 5.1. Note that in this case α_{i} and β_{i} are even orthonormal directions.

Remark 5. If one would like to use the Hermite basis functions introduced in this paper, but certain derivative data is not available, it is possible to use any of the well established methods for derivative generation. For a survey of these methods see e.g. [1].

We close by illustrating Example 5.1 numerically. To this end we use spherical multiquadrics with parameter R = 0.4, and choose the test function f as

f(x, y, z) = \frac{1 + 10xyz + x^{8} + e^{2y^{3}} + e^{2z^{2}}}{14} , (x, y, z) ∈ S^{2}.

We take 18 nearly equally spaced points on the sphere (obtained by subdividing a regular tetrahedron) and evaluate the Hermite interpolant on a finer grid of 3846

Hermite (part.) Hermite (dir.) Lagrange (18) Lagrange (54)
max. error  1.218017e-02  1.218017e-02 1.493461e-01 1.335387e-02
rms. error  1.233846e-05  1.233846e-05 4.226617e-03 1.251992e-05
cond(A)  2.970847e+19 3.508120e+03 1.443921e+03 1.438680e+10

Table 1. Numerical comparison.

points. The table below contains maximum and root-mean-square errors on this grid, as well as the approximate condition number of the interpolation matrices. In order to provide some insight into the advantages of Hermite basis functions we include the results for Hermite interpolation according to Proposition 5.3, as well as those for Lagrange interpolation at the 18 points and at 54 points which were generated by clustering 3 points at each of the 18 original locations.

As noted above, the interpolation matrix for Hermite interpolation based on the use of partial derivatives with respect to the spherical coordinates results in a singular matrix. Nevertheless, the errors for both Hermite interpolants are the same. The condition number observed for the Lagrange interpolant with clustered points shows that it is highly recommended to replace such clusters by one representative point, and to use Hermite basis functions there instead of the usual (Lagrange) radial basis functions. This idea was incorporated in an adaptive least squares fitting method in [6]. The basic algorithm (without the Hermite functions) is also described in [7].

References

1. Alfeld, P., Scattered data interpolation in three or more variables, in Mathematical Methods in Computer Aided Geometric Design, T. Lyche and L. Schumaker (eds), Academic Press, New York, 1989, 1–33.

2. Alfeld, P., M. Neamtu, and L. L. Schumaker, Fitting scattered data on spherelike surfaces using spherical splines, J. Comp. Appl. Math. 73 (1996), 5–43.

3. Donoghue, W. F., Monotone Matrix Functions and Analytic Continuation, Springer, New York, 1974.

4. Dunford, N. and J. T. Schwartz, Linear Operators Part I: General Theory, Interscience Publishers, New York, 1966.

5. Dyn, N., T. Goodman, and C. A. Micchelli, Positive powers of certain conditionally negative definite matrices, Indaga. Math. Proc. Ser.A. 89 (1986), 163–178.

6. Fasshauer, G. E., Radial Basis Functions on Spheres, dissertation, Vanderbilt University, 1995.

7. Fasshauer, G. E., Adaptive least squares fitting with radial basis functions on the sphere, in Mathematical Methods for Curves and Surfaces, M. Dæhlen, T. Lyche, and L. Schumaker (eds), Vanderbilt University Press, Nashville & London, 1995, 141–150.

8. Freeden, W., Spline methods in geodetic approximation problems, Math. Meth. Appl. Sci. 4 (1982), 382–396.

9. Freeden, W., A spline interpolation method for solving boundary value problems of potential theory from discretely given data, Num. Meth. Part. Diff. Eq. 3 (1987), 375–398.

10. Hardy, R. L., Research results in the application of multiquadric equations to surveying and mapping problems, Survg. Mapp. 35 (1975), 321–332.

11. Hardy, R. L., Theory and applications of the multiquadric-biharmonic method, Comput. Math. Appl. 19 (1990), 163–208.

12. Horn, R. A. and C. R. Johnson, Matrix Analysis, Cambridge University Press, Cambridge, 1985.

13. Menegatto, V. A., Interpolation on spherical domains, dissertation, University of Texas Austin, 1992.

14. Menegatto, V. A., Strictly positive definite kernels on the Hilbert sphere, Appl. Anal. 55 (1994), 91–101.

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

16. Müller, C., Spherical Harmonics, Springer Lecture Notes in Mathematics, Vol. 17, 1966.

17. Narcowich, F. J., Generalized Hermite interpolation and positive definite kernels on a Riemannian manifold, J. Math. Anal. Appl. 190 (1995), 165–193.

18. Narcowich, F. J. and J. D. Ward, Generalized Hermite interpolation via matrixvalued conditionally positive definite functions, Math. Comp. 63 (1994), 661– 687.

19. Powell, M. J. D., The theory of radial basis functions in 1990, in Advances in Numerical Analysis II: Wavelets, Subdivision, and Radial Basis Functions, W. Light (ed), Oxford University Press, Oxford, 1992, 105–210.

20. Schoenberg, I. J., On certain metric spaces arising from euclidean spaces by a change of metric and their imbedding in Hilbert space, Ann. of Math. 38 (1937), 787–793.

21. Schoenberg, I. J., Metric spaces and positive definite functions, Trans. Amer. Math. Soc. 44 (1938), 522–536.

22. Schoenberg, I. J., Metric spaces and completely monotone functions, Ann. of Math. 39 (1938), 811–841.

23. Schoenberg, I. J., Positive definite functions on spheres, Duke Math. J. 9 (1942), 96–108.

24. Stewart, J., Positive definite functions and generalizations, an historical survey, Rocky Mountain J. Math. 6 (1976), 409–434.

25. Sun, X., Scattered Hermite interpolation using radial basis functions, Linear Algebra Appl. 207 (1994), 135–146.

26. Szegő, G., Orthogonal Polynomials, Amer. Math. Soc. Coll. Publ., Vol. XXIII, 1959.

27. Wells, J. H. and R. L. Williams, Embeddings and Extensions in Analysis, Springer, New York, (1975).

28. Widder, D. V., The Laplace Transform, Princeton University Press, Princeton, 1941.

29. Wu, Z., Hermite-Birkhoff interpolation of scattered data by radial basis functions, Approx. Theory Appl. 8 (1992), 1–10.

Gregory E. Fasshauer Department of Computer Science and Applied Mathematics Illinois Institute of Technology Chicago, IL 60616

fass@amadeus.csam.iit.edu