In practical applications over a wide field of study one often faces the problem of reconstructing an unknown function f from a finite set of discrete data. These data consist of data sites X = {x1, …, xN} and data values fj = f(xj), 1 ≤ j ≤ N, and the reconstruction has to approximate the data values at the data sites. In other words, a function s is sought that either interpolates the data, i.e. that satisfies s(xj) = fj, 1 ≤ j ≤ N, or at least approximates the data, s(xj) ≈ fj. The latter case is in particular important if the data contain noise.

In many cases the data sites are scattered, i.e. they bear no regular structure at all, and there is a very large number of them, easily up to several million. In some applications, the data sites also exist in a space of very high dimensions. Hence, for a unifying approach methods have to be developed which are capable of meeting this situation. But before pursuing this any further let us have a closer look at some possible applications.

Surface reconstruction

Probably the most obvious application of scattered data interpolation and approximation is the reconstruction of a surface S. Here, it is crucial to distinguish between explicit and implicit surfaces.

Scattered data approximation is a recent, fast growing research area. It deals with the problem of reconstructing an unknown function from given scattered data. Naturally, it has many applications, such as terrain modeling, surface reconstruction, fluid-structure interaction, the numerical solution of partial differential equations, kernel learning, and parameter estimation, to name a few. Moreover, these applications come from such different fields as applied mathematics, computer science, geology, biology, engineering, and even business studies.

This book is designed to give a thorough, self-contained introduction to the field of multivariate scattered data approximation without neglecting the most recent results.

Having the above-mentioned applications in mind, it immediately follows that any competing method has to be capable of dealing with a very large number of data points in an arbitrary number of space dimensions, which might bear no regularity at all and which might even change position with time.

Hence, in my personal opinion a true scattered data method has to be meshless. This is an assumption that might be challenged but it will be the fundamental assumption throughout this book. Consequently, certain methods, that generally require a mesh, such as those using wavelets, multivariate splines, finite elements, box splines, etc. are immediately ruled out. This does not at all mean that such methods cannot sometimes be used successfully in the context of scattered data approximation; on the contrary, it just explains why these methods are not discussed in this book.

Up to now we have dealt only with the problem of recovering an unknown function from certain known function values. But sometimes it might be desirable to recover the function also from other types of data. For example, the function's derivatives might be known at certain points, but not the function itself. This becomes interesting if partial differential equations are considered.

In this chapter we deal with a more general problem than those we have discussed so far. Our approach includes in particular collocation and Galerkin methods for solving partial differential equations. But the methods we will derive are at the present time only able to compete with classical methods to a certain extent. In any case, whenever large data sets are considered one has to combine the methods introduced below with the fast-evaluation ideas of Chapter 15.

In this sense, this chapter should be seen as a unified introduction to a general class of recovery problems.

Optimal recovery in Hilbert spaces

We start by generalizing results from Chapters 11 and 13. In this section we restrict ourselves for simplicity to the Hilbert space setting, even though everything works in the case of semi-Hilbert function spaces also.

Let H be a Hilbert space and denote its dual by H*. Suppose that Λ = λ1, … λN ⊆ ℝ are certain given values. Then a generalized recovery problem would seek to find a function s ℝ H such that λj(s) = fj, 1 ≤ j ≤ N. We will call s a generalized interpolant.

Every book on numerical analysis has a chapter on univariate polynomial interpolation. For given data sites x1 < x2 < … < xN and function values f1, …, fN, there exists exactly one polynomial pf ∈ ΦN−1(ℝ1) that interpolates the data at the data sites. One of the interesting aspects of polynomial interpolation is that the space ΦN−1(ℝ1) depends neither on the data sites nor on the function values but only on the number of points. While the latter is intuitively clear, it is quite surprising that the space can also be chosen independently of the data sites. In approximation theory the concept whereby the space is independent of the data has been generalized and we will shortly review the necessary results. Unfortunately, it will turn out that the situation is not as favorable in the multivariate case.

The Mairhuber–Curtis theorem

Motivated by the univariate polynomials we give

Definition 2.1Suppose that Ω ⊆ ℝd contains at least N points. Let V ⊆ C(Ω) be an N-dimensional linear space. Then V is called a Haar space of dimension N on Ω if for arbitrary distinct points x1, …, xN ∈ Ω and arbitrary f1, …, fN ∈ ℝ there exists exactly one function s ∈ V with s(xi) = fi, 1 ≤ i ≤ N.

For our investigations on radial basis functions in the following chapters it is crucial to collect several results from different branches of mathematics. Hence, this chapter is a repository of such results. In particular, we are concerned with special functions such as Bessel functions and the Γ-function. We discuss the features of Fourier transforms and give an introduction to the aspects of measure theory that are relevant for our purposes. The reader who is not interested in these technical matters could skip this chapter and come back to it whenever necessary. Nonetheless, it is strongly recommended that one should at least have a look at the definition of Fourier transforms to become familiar with the notation we use. Because of the diversity of results presented here, we cannot give proofs in every case.

Bessel functions

Bessel functions will play an important role in what follows. Most of what we discuss here can be found in the fundamental book [187] by Watson.

The starting point for introducing Bessel functions is to remind the reader of the classical Γ-function and some of its features.

In numerical analysis, the concept of locally supported basis functions is of general importance. Several function spaces used for approximation possess locally supported basis functions. The most prominent examples in the one-dimensional case are the well-known B-splines. The general advantages of compactly supported basis functions are a sparse interpolation matrix on the one hand, and the possibility of a fast evaluation of the interpolant on the other.

Thus, it seems to be natural to look for locally supported functions also in the context of radial basis function interpolation and we will give an introduction to this field in this chapter.

At the outset, though, we want to point out one crucial difference from classical spline theory. While the support radius of the B-splines can be chosen proportional to the maximal distance between two neighboring centers, something similar will not lead to a convergent scheme in the theory of radial basis functions. The correct choice of the support radius is a very delicate question, which we will address in a later chapter on numerical methods.

General remarks

Gaussians, (inverse) multiquadrics, powers, and thin-plate splines share two joint features. They are all globally supported and are positive definite on every ℝdx. The truncated powers from Theorem 6.20, however, are compactly supported but are also restricted to a finite number of space dimensions.

We will see that the two features are connected. But let us first comment on conditionally positive definite functions.

So far w have been concerned with interpolation on an arbitrary domain Ω ⊂ ℝd. However, we have not used any topological information about Ω. Instead, we have employed only the fact that it is a subset of ℝd. As a matter of fact, without having more information on Ω this is the only way. But many applications provide us with additional information on the underlying domain. For example, problems coming from geology often relate to the entire earth, so that the unit sphere would be an appropriate model and the additional information should lead to a better approximant.

Hence, in this chapter, we want to give an introduction to the theory of scattered data interpolation on spheres and other compact manifolds by radial or zonal functions.

Spherical harmonics

Generally, functions on the sphere are expressed as Fourier series with respect to an orthonormal family called spherical harmonics. In this section we will review the results on these functions. Since this material is only necessary for the present chapter we did not incorporate it into Chapter 5. Moreover, we have to skip the proofs once again. The interested reader is referred to Müller's book [140].

The domain of interest is the d-variate unit sphere Sd−1 := {x ∈ ℝd : ∥x∥2 = 1} ⊆ ℝd. It has surface area

On Sd−1, we will employ the usual inner product

where dS(x) is given by the standard measure on the sphere.

The distance between two points is the geodesic distance, which is the length of the shorter part of the great circle joining x and y or, in other words, dist(x, y) = arccos(xT y).

So far, we have dealt with the following simple interpolation or approximation problem. An in general unknown function f is specified only at certain points X = {x1, …, xN}, and we are interested in recovering the function f on a region Ω that is well covered by the centers X. In a later chapter we will concentrate on more general problems. But let us stick to this particular one a little longer. Why should we use (conditionally) positive definite kernels for recovering f?

We have already learnt that recovering f is a difficult task and that radial basis functions are a powerful tool for doing this. In particular, they can be used (at least theoretically – we come back to the numerical treatment in a later chapter) with truly scattered data and in every dimension. Moreover, positive definite functions appeared quite naturally in the context of reproducing-kernel Hilbert spaces.

But this is not the end of the story. Interpolants based on (conditionally) positive definite kernels are optimal in several other ways and the present chapter is devoted to this subject.

Minimal properties of radial basis functions

Let us start with best approximation. We have seen that the native space NΦ(Ω) corresponding to a (conditionally) positive definite kernel Φ is an adequate function space. The interpolant sf,X is one candidate that uses the given information about f on X, but of course not the only one.

So far we have dealt mainly with the theoretical background of scattered data approximation and interpolation. Except for the moving least squares approximation we have not discussed an efficient implementation of our theory. But from the investigation of the condition number of straightforward interpolation matrices we know that special techniques have to be developed to get an efficient but still accurate approximation method. Several approaches will be discussed in the next chapter.

Crucial to all these methods is the choice of the underlying data structure. It is worth reflecting for a while on how the centers can be stored in a computer most successfully. Thus, in this chapter we will discuss different ways of representing the centers X = {x1, …, xN} ⊆ Ω ⊆ ℝd.

Let us start by collecting possible requests about the set of points X that the data structure should be able to answer efficiently.

The first question every user has to answer is whether all points should be kept in the main memory of the computer or whether the number of points is so large that it has to be stored on a hard disk or other external device. The difference is that in the latter case a reasonable ordering of the data points reduces the number of disk accesses, resulting in a dramatic reduction in run time. There is also an improvement in the first situation by a reasonable ordering, because of the cache of the computer, but it is not dramatic. We will concentrate on the situation where all points can be kept in the main memory of the computer.

The crucial point in local polynomial reproduction is the compact support of the basis functions uj. To be more precise, all supports have to be of the same size. The local support of the uj means that data points far away from the current point of interest x have no influence on the function value at x. This is often a reasonable assumption.

The last chapter did not answer the question how to construct families with local polynomial reproductions efficiently. The moving least squares method provided in this chapter forms an example of this.

Definition and characterization

Suppose again that discrete values of a function f are given at certain data sites X = {x1, …, xN} ⊆ Ω ⊆ ℝd. Throughout this chapter Ω is supposed to satisfy an interior cone condition with angle θ and radius r.

The idea of the moving least squares approximation is to solve for every point x a locally weighted least squares problem. This appears to be quite expensive at first sight, but it will turn out to be a very efficient method. Moreover, in many applications one is only interested in a few evaluations. For such applications the moving least squares approximation is even more attractive, because it is not necessary to set up and solve a large system.

The influence of the data points is governed by a weight function w :Ω × Ω → ℝ, which becomes smaller the further away its arguments are from each other.