Positive Approximation and Interpolation Using Compactly Supported Radial Basis Functions

Abstract

We discuss the problem of constrained approximation and interpolation of scattered data by using compactly supported radial basis functions, subjected to the constraint of preserving positivity. The approaches are presented to compute positive approximation and interpolation by solving the two corresponding optimization problems. Numerical experiments are provided to illustrate that the proposed method is flexible.

1. Introduction

The problem of scattered data interpolation consists of constructing a function that interpolates data values which are known at some scattered points. However, we often have some additional information that we wish to confine to interpolation. For example, we know the quantity from which the data is sampled, positive, monotonic, or convex. Thus, it is important to construct a function which satisfies the underlying constraints.

In the past twenty years, the problem of positivity-preserving scattered data interpolation has been considered by many researchers. One class of methods for scattered data interpolation require the points to be triangulated mainly by using piecewise (rational) cubic interpolation and derived the conditions for preserving the shape of positive data. For related literature, we refer to the papers [1–6] and the references therein.

Another major class of methods for scattered data interpolation does not involve any prior triangulation, and can be viewed as meshless method. The two main types are radial basis functions (RBFs), which include multiquadrics and thin-plate splines, and so forth, and Shepard-type methods, which include the modified quadratic Shepard approach. Both types, RBF and Shepard, are widely used in practical fields.

However, there has been relatively little work done on the imposition of constraints for these meshless methods. For RBF, Utreras [7] considered how positivity can be imposed as a constraint in case of thin plate spline interpolation for two-dimensional data. For Shepard, Asim and Brodlie [8, 9] discussed the modified quadratic Shepard method, which interpolates scattered data of any dimensionality and can be constrained to preserve positivity and more general constraint.

In this paper, we show how radial basis functions method, which interpolates or approximates scattered data, can be constrained to preserve positivity.

2. Radial Basis Functions

A radial basis function (RBF) is a relatively simple multivariate function generated by a univariate function. Due to its simple form and good approximation behavior, the radial basis function approach has become an effective tool for multivariate scattered data interpolation during the last two decades [10–14].

For any given scattered data (𝑥_𝑗,𝑓_𝑗)∈ℝ^𝑠×ℝ,𝑗=1,…,𝑚, where points 𝑥₁,…,𝑥_𝑚∈ℝ^𝑠 are pairwise distinct, the so-called radial basis function interpolation is to use a function 𝜙∶ℝ₊→ℝ to construct the interpolant in the form of 𝑠(𝑥)=^𝑚_𝑗=1𝜆_𝑗𝜙(‖𝑥−𝑥_𝑗‖), (2.1) satisfying𝑠(𝑥_𝑖)=^𝑚_𝑗=1𝜆_𝑗𝜙(‖𝑥_𝑖−𝑥_𝑗‖)=𝑓_𝑖,𝑖=1,…,𝑚. (2.2)

The positive definiteness of 𝜙 guarantees that the above interpolation problem (2.2) possesses a unique solution and refers to 𝜙 as a classical radial basis function. If 𝜙 has compact support, then the positive definite linear system is sparse and reduces computational cost greatly. Thus, we bypass this problem by restricting 𝜙 to have compact support.

Compactly supported radial basis functions (CSRBFs) have only recently been constructed. Wu first constructed a broad variety of CSRBF [15]. Shortly after this, Wendland constructed these functions such that they possess the lowest degree among all CSRBFs which are positive definiteness for given space dimension and prescribed order of smoothness [16]. They are radial basis functions which are positive definite on ℝ^𝑠 for a given space dimension 𝑠 (𝑃𝐷_𝑠), belong to a prescribed smoothness class (𝐶^2𝑘), are compactly supported and easy to evaluate. Some examples of such radial basis functions are given in Table 1.

Table 1 

Some of Wendland's CSRBF 𝜙_𝑠,𝑘∈𝑃𝐷_𝑠∩𝐶^2𝑘.

It is a useful property and provides a good selection of Wendland's functions with respect to the order of continuity and the dimension of space. Thus, CSRBFs have become a popular tool for multivariate interpolation of large scattered data [17, 18].

In order to adapt the interpolation to scattered data of different densities, it is necessary to be able to scale the support of 𝜙. So from now on we assume that the radius 𝛼 of support of 𝜙 is one and replace 𝜙 by 𝜙_𝛼⋅(⋅)=𝜙(𝛼),for𝛼>0. (2.3)

Meanwhile, in order to achieve both the best possible approximation behavior and best possible stability with respect to the support of CSRBF, we adopt the strategy to choose the radius 𝛼 as done in paper [17] throughout this paper.

3. Positive Approximation for Positive Data

The problem we are addressing is the following. Given a set of 𝑚 scattered data 𝑥_𝑖,𝑥_𝑖∈ℝ^𝑠 with associated data values 𝑓_𝑖,𝑖=1,2,…,𝑚, where 𝑓_𝑖≥0, we seek an optimal approximating function 𝐹(𝑥) such that 𝐹(𝑥)≥0 in least-squares sense.

we first construct the function 𝐹(𝑥) as follows:𝐹(𝑥)=^𝑁_𝑖=1𝜆_𝑖𝜙_𝛼(‖𝑥−𝑐_𝑖‖), (3.1) where 𝜙_𝛼(‖𝑥−𝑐_𝑖‖) is a compactly supported radial basis function centered on 𝑐_𝑖 with radius 𝛼 of support. Usually, the number of 𝑁 is less than 𝑚 greatly and centers 𝑐_𝑖 are different from data points 𝑥_𝑖.

Thus, the so-called positive approximation for positive data is reduced to solve the following optimization problem min^𝑚_𝑗=1(𝐹(𝑥_𝑗)−𝑓_𝑗)²,subjectto𝐹(𝑥)≥0. (3.2)

Note that all 𝜙_𝛼(‖𝑥−𝑐_𝑖‖) are positive and 𝐹(𝑥) is a linear combination of 𝜙_𝛼(‖𝑥−𝑐_𝑖‖). So, we can restrict each 𝜆_𝑖 to be positive in order to guarantee the function to be everywhere positive. This sufficient condition is the kernel idea of our proposed method.

Therefore, the problem of positive approximation is transformed into a quadratic optimization problem subjected to linear constraintsmin^𝑚_𝑗=1(𝐹(𝑥_𝑗)−𝑓_𝑗)²,subjectto𝜆₁≥0,…,𝜆_𝑁≥0. (3.3)

Without loss of generality, we illustrate it with a very simple example in 1D.

Example 3.1. The set of 7 data points in Table 2 shows the velocity of wind. The velocity is inherently positive, and we therefore require the resulted approximating function to preserve this property. Here, we choose the Wendland's function as 𝜙(𝑟)=(1−𝑟)⁴₊(4𝑟+1). Figure 1 shows the positive approximation curve applied to the set of data for the velocity of wind using our proposed method (3.3) when 𝑐_𝑖 are chosen the same as 𝑥_𝑖. Table 3 shows the values of Error=∑⁷_𝑗=1(𝐹(𝑥_𝑗)−𝑓_𝑗)² when the number of centers 𝑁 is varied. Of course, the total error will become less if the centers 𝑐_𝑖 are chosen properly.

Table 2 

One dimensional data for the velocity of wind.

Table 3 

Total error for different parametric values of 𝑁.

Figure 1 

Positive approximation curve for data of wind velocity.

Remark 3.2. The proposed method has a good chance to work well if the data are coming from a function 𝑓 which is a convolution of the kernel with a positive function 𝑔, because then 𝑓 can be recovered well by an integration formula which has exactly the form of (3.1).

4. Positive Interpolation for Positive Data

The problem of positive interpolation we are addressing is the following. Given a set of 𝑚 scattered data 𝑥_𝑖,𝑥_𝑖∈ℝ^𝑠 with associated data values 𝑓_𝑖,𝑖=1,2,…,𝑚, where 𝑓_𝑖≥0, the so-called positive interpolation is to construct a function 𝐹(𝑥) such that 𝐹(𝑥_𝑖)=𝑓_𝑖,𝑖=1,…,𝑚 and 𝐹(𝑥)≥0, where 𝐹(𝑥) is a linear combination of the CSRBF.

The basic idea is also to restrict each combination coefficient 𝜆_𝑖 to be positive in order to guarantee that 𝐹(𝑥) is positive everywhere. Obviously, classical RBF interpolation (2.2) cannot guarantee that all 𝜆_𝑖 are positive.

One feasible way is to choose 𝑛 new added data 𝑦_𝑖,𝑖=𝑚+1,…,𝑚+𝑛 and then turn to find the following interpolant:𝐹(𝑥)=^𝑚_𝑖=1𝜆_𝑖𝜙_𝛼(‖𝑥−𝑥_𝑖‖)+^𝑚+𝑛_𝑗=𝑚+1𝜆_𝑗𝜙_𝛽𝑗(‖𝑥−𝑦_𝑗‖) (4.1) such that𝐹(𝑥_𝑖)=𝑓_𝑖𝜆,𝑖=1,2,…,𝑚,_𝑖≥0,𝑖=1,2,…,𝑚+𝑛. (4.2)

The above interpolant will naturally arise an inevitable problem: what makes the system of linear equations (4.1) have a unique solution 𝐹(𝑥) satisfying the linear constraints. In order to solve it in special cases, we require the following lemma.

Lemma 4.1 (see [19, Gordan's Theorem]). Let 𝐴₁,𝐴₂,…,𝐴_𝑚 be n-dimensional vectors; there does not exist a vector 𝑃 such that 𝐴^𝑇_𝑖𝑃<0 if and only if there exist nonnegative real numbers 𝜆₁,𝜆₂,…,𝜆_𝑚 such that ∑^𝑚_𝑖=1𝜆_𝑖𝐴_𝑖=0.

Undoubtedly, it is hard to give a general result which can determine whether there exists a solution when the added centers 𝑦_𝑗 and the radius 𝛽_𝑗 are chosen randomly. Thus, we discuss the following special case to choose centers 𝑦_𝑗 and the radius 𝛽_𝑗.

For each data 𝑥_𝑖,𝑖=1,…,𝑚, if we choose the new added data 𝑦_𝑖 in the neighborhood of 𝑥_𝑖 randomly and choose the radius 𝛽_𝑖 such that the influence domain of 𝜙_𝛽𝑖(‖𝑥−𝑦_𝑖‖) only contains the data 𝑥_𝑖, then for this case we have the following.

Theorem 4.2. If the new added data 𝑦_𝑗 and radius 𝛽_𝑗 are chosen as above, then there exists a solution subjected to the constraints (4.2).

Proof. Let 𝐴^𝑇_𝑗=−{𝜙_𝛽𝑗(‖𝑥_𝑖−𝑦_𝑗‖)}^𝑚_𝑖=1=(0,…,0,−𝑏_𝑗𝐴,0,…,0),𝑗=𝑚+1,…,2𝑚,_𝑇𝑗=−{𝜙_𝛼(‖𝑥_𝑖−𝑥_𝑗‖)}^𝑚_𝑖=1,𝑗=1,2,…,𝑚,𝐴^𝑇_2𝑚+1=(𝑓₁,…,𝑓_𝑚),∀𝑓_𝑖≥0, (4.3) where 𝑏_𝑗,𝑗=𝑚+1,…,2𝑚, are positive real numbers.
It is easy to see that there does not exist a vector 𝑃 such that 𝐴^𝑇_𝑖𝑃<0,𝑖=𝑚+1…,2𝑚+1. According to Gordon lemma, we know that there exist nonnegative real numbers 𝜆₁,𝜆₂,…,𝜆_2𝑚+1 such that ∑^2𝑚+1_𝑖=1𝜆_𝑖𝐴_𝑖=0 and 𝜆_2𝑚+1>0.
We naturally wish that there are many zero coefficients 𝜆_𝑗,𝑗=𝑚+1,…,𝑚+𝑛, in the formulas (4.1). That is to say, we add as few data as possible. For example, if the original interpolation function is positive, then we certainly need not add any new data. So from now on we hope to minimize the objective function sgn(𝜆_𝑚+1)+⋯+sgn(𝜆_2𝑚), where sgn(⋅) denotes the sign function.
Therefore, an approach is presented to compute the positive interpolation by solving the following optimization problem:min(sgn(𝜆_𝑚+1)+⋯+sgn(𝜆_2𝑚)),subjectto𝐹(𝑥_𝑖)=𝑓_𝑖𝜆,𝑖=1,2,…,𝑚,_𝑖≥0,𝑖=1,2,…,2𝑚. (4.4)
We illustrate it with an example using the same data given in Table 2 and make comparisons with several existing methods.

Example 4.3 (Example 3.1 is continued). In Figure 2, we show the original interpolation curve by using CSRBF without the constraint of positivity. Generally, the curve has good approximation behavior, but the curve goes beyond the range of the data values which makes physical nonsense and indeed the positivity is violated.

Figure 2 

One-dimensional data for the velocity of wind—CSRBF interpolation loses positivity.

In Figure 3, the positive curve applied to the same set of positive data using piecewise cubic method [8] is shown. Meanwhile, we reveal the positive curve using constrained modified quadratic Shepard method (CMQS) [9] in Figure 4.

Figure 3 

Positive interpolation curve for data of wind velocity using piecewise cubic.

Figure 4 

Positive interpolation curve for data of wind velocity using CMQS.

By contrast, in Figure 5, we show the positive interpolation curve applied to the same data of wind velocity using our proposed approach (4.1). From the experiment, the constructed curve possesses relatively good approximation behavior in contrast to piecewise cubic method and constrained modified quadratic Shepard method.

Figure 5 

Positive interpolation curve for positive data.

Remark 4.4. It is pointed out that we generally add corresponding new data 𝑦_𝑚+𝑘 if 𝑓_𝑘 is the minimal data among all value 𝑓_𝑖,𝑖=1,…,𝑚 instead of adding 𝑚 new data directly to solve the positive interpolation in practical use.

5. Numerical Example

The strong and sufficient condition that all 𝜆_𝑖 are positive may degrade the quality of the interpolating function, in comparison with the original unconstrained CSRBF interpolation. We evaluate the quality of the new interpolant by calculating the variances between the exact and calculated values on a set of test points by the following experiment in 2D.

Example 5.1. We illustrate it with the following function 𝑆(𝑥,𝑦) [9] which is defined as: 𝑆(𝑥,𝑦)=1,if(𝑦−𝑥)≥0.5,2(𝑦−𝑥),if0.5≥(𝑦−𝑥)≥0,cos(4𝜋𝑟)+121,if𝑟≤4,0,otherwise, (5.1) where 𝑟=(𝑥−1.5)²+(𝑦−0.5)².
The interpolated function 𝑆(𝑥,𝑦) is shown in Figure 6. Meanwhile, we show the positive surface generated by our proposed method at a random set of 100 points in Figure 7.
In order to measure the quality of new interpolant, we carry out the following experiment. Firstly, we construct two interpolants by constrained modified quadratic Shepard method and our proposed method, based on a series of data sets where the number of randomly chosen points increases. Secondly, we evaluate the interpolants on a grid of 25×25 points and calculate the variances between the exact and calculated values. The results are shown in Table 4.

Table 4 

Variances for two interpolants for 𝑆(𝑥,𝑦).

Figure 6 

Interpolated function 𝑆(𝑥,𝑦).

Figure 7 

Positive interpolating surface for 𝑆(𝑥,𝑦).

6. Conclusion

Positivity preserving interpolation is an interesting work. In this paper, we discuss the approximation and interpolation problem of scattered data by using CSRBF under the constraint of positivity. The approaches are presented to compute positive approximation and interpolation by solving the two corresponding optimization problems.

However, we have not discussed how to select the optimal number and position of the new added data 𝑦_𝑖 in order to both achieve the existence and better possible approximation behavior of interpolation. These problems would be attractive and remain to be our future work.

Acknowledgments

This work was supported by the Natural Science Foundation of Zhejiang Province (nos. Y7080068 and Y6090211), the Foundation of Department of Education of Zhejiang Province (nos. Y200802999 and Y200907579), and Tianyuan Foundation for Mathematics of China (no. 10926035). The authors are very grateful to the referees for their careful reading of the manuscript and many valuable comments and remarks, which greatly improved this paper.