Abstract

In this study, researchers describe and implement a parameter free adaptive algorithm for radial basis function interpolation of real valued bump functions. The method is based on thin plate splines which are known for their good approximation properties. Interpolation points are added and coarsened based on residuals computed on a finer discretization. Numerical examples in one dimension are presented to show the effectiveness of the method.

---

INTRODUCTION

If X is a finite set of scattered data points in R^d and if f = R^d→R is a function known only on X, then an important practical problem is that of finding a smooth function that interpolates f and provides to good approximation to f near X. To this end, interpolation by Radial Basis Functions (RBFs) has become a very powerful tool in multivariate approximation theory. In the past two decades, radial basis functions have been used extensively in fields like the numerical solution of partial differential equations like in Franke (1998) and Lorentz (2003).

They have also been applied in non-uniform sampling, mathematical finance, computer graphics and in optimization. The efficiency of the RBF interpolation on scattered and uniform data has been illustrated by several numerical examples and theoretical results by Wendland (2005) and Fasshauer (2007).

For problems with steep gradients, singularities and other topological changes, it is natural to use adaptive methods to reduce the error of the interpolation and approximation. Adaptive methods have been implemented for RBF interpolation by Behrens and Iske (2002), Sarra (2005) and Driscoll and Heryudono (2007) amongst others.

Although, several radial basis function like Gaussians and multiqudratics have very good approximation properties, they have a shape parameter and numerical experiments in Fasshauer (2007) and elsewhere show how the choice of the shape parameter has profound influence on both the accuracy and numerical stability of the solution of an interpolation problem. On the other hand, certain RBFs like radial powers and surface splines have the advantage of being shape parameter free. This means that one does not need to worry about the choice of a good shape parameter.

In Driscoll and Heryudono (2007), an adaptive residual subsampling method based on radial basis functions was proposed and implemented. In their research, they used multiquadrics which contain a shape parameter and they show how the value of the shape parameter affects the accuracy of their method. To this end, they had to adaptively modify the shape parameter during simulation.

In this study, the researchers will use a parameter free radial basis function, the thin plate splines (a type of surface splines), to implement an adaptive radial basis function interpolation scheme for bump functions and show how good results can be obtained without the need of adapting the shape parameter. A bump function is a function that is non-zero only on a compact subset of the interior of the domain of definition. The utility of thin plate splines interpolation in this context was highlighted by Gutzmer and Iske (1997) and Kaser (2003).

Radial basis function interpolation: Consider a data vector f = (f (x₁),...f (x_n))^T ∈ Rⁿ function values, obtained from a function f : R^d→R at a scattered finite point set X = {x₁,...x_n} ⊂ R^d, d≥1. Scattered data interpolation requires computing an appropriate interpolant s: R^d→R satisfying;

(1)

The radial basis function interpolation method utilizes a fixed radial function φ: [0, ∞] → R so that the interpolant s in Eq. 1 has the form:

(2)

Where ç -ç is the Euclidean norm on R^d and is the space of all polynomials in d variables of degree at most m-1 (order m), where m = m (φ) is known as the order of the basis function φ. Possible choices for φ along with their orders are shown in Table 1. When m = 0, the interpolant s in Eq. 2 has the form:

(3)

Using the interpolation conditions Eq. 1, the coefficients c = (c₁,..., c_n)^T of s in Eq. 3 can be obtained by solving the linear system:

Where for m>0, the interpolant s in Eq. 2 contains a nonzero polynomial part, yielding q additional degree of freedom where q is the dimension of the polynomial space . These additional degrees of freedom are usually eliminated using q vanishing moment conditions:

(4)

In total, this amounts to solving the (n + q)x(n + q) linear system:

(5)

Where:

Table 1:	Radial Basis Functions (RBFs), their parameters and order

for the coefficients of the polynomial part in Eq. 2. The important case is the thin plate spline (surface splines with k = 2). In this case, for one space dimension, we have φ (r) = r² log Υ with:

so that;

(6)

This requires us solving an (n + 2) x (n + 2) linear system to obtain the coefficients c₁,...c_n, d_0, d₁. The effectiveness of this interpolant for bump functions was shown by Bejancu (2001). We note, however that parameter free radial basis functions will not be able to achieve the spectral convergence rates that are possible with other basis functions such as Gaussians and multiquadrics. Moreover, the interpolation matrices are larger in this case because of the polynomial part (Fasshauer, 2007).

The adaptive interpolation method: We will use the residual subsampling of Driscoll and Heryudono (2007) coupled with the strategy of Kaser (2003) to provide a parameter free adaptive interpolation method.

A function f (x) defined on an interval (a, b) is picked and an initial discretization giving a set X_o of n_o equally spaced points is obtained. A thin plate spline interpolant s(x) is then obtained for the function f(x). The interpolation error is then computed at the set M_o of n_o^-1 midpoints of the subintervals. A midpoint m_i ∈ M_o becomes an interpolation point if the error at m_i exceeds a threshold value '_Υ and an interpolation x_i ∈ X_o is coarsened or removed is it lies between two points in M_o whose error is below a smaller threshold '_c.

This gives a new discretization X₁ of theinterval (a, b). This algorithm is repeated and further discretizations X₂, X₃, X₄,.... are generated until no further points are added or coarsened. For the case of thin plate spline interpolation, the local error estimate at x ε X is according to Wu and Schaback (1993) of the form:

(7)

Where C>0 is a constant depending on f and (for some radius p>0:

is the fill distance of X around x, with:

being the Euclidean distance between the point y and the set X. As shown by Wu and Schaback (1993), the reduction in the error (Eq. 7) around any interpolation point is accomplished by reducing the distance function;

in a neighbourhood of x_j. This shows that adaptivity can indeed reduce the interpolation error.

Numerical examples and discussions
Example 1: We will use a simple function f(x) = exp(–40x²) on the interval [–1, 1] with '_r = 1.5x10^-5 and 'α = 1x10^-6.. The results are shown in the Table 2 where cond(A) denotes the condition number of the interpolation matrix, E_∞ denotes the error in the infinity norm, na the number of points added, nc the number of points coarsened and n denotes the number of interpolation points.

Table 2:	Node distribution, error and condition number for Example 1

The method converges in 8 iterations and the errors are quite small. Moreover, the method is stable looking at the condition numbers of the interpolation matrices.

These results are similar to the multiquadric method while the conditioning of the interpolation matrix is even better. The plots for f(x), the node distribution at the final iteration and the errors are shown in Fig. 1a-c, respectively. We observe that the method adaptvely places more interpolation points close to where the function has sharp features.

We also note that when interpolating on 111 equally spaced points the error will be 1.1656x10^-4 with the condition number 2.6166x10⁶ which shows that the adaptive method gives superior results for the same number of interpolation points.

Example 2: We will consider the bump function given by the poduct (Bejancu, 2001):

Table 3:	Node distribution, error and condition number for Example 2

 

Fig. 1:	a) Plot of f (x) = exp(¯40x2), b) node distribution at final iteration, c) number of interpolation points vs. error

 

 

Fig. 2:	a) Plot of gp(x), b) node distribution at final iteratuion, c) number of interpolation points vs. error

 

(8)

Therefore, and the support of g_p is the interval:

We will still use '_Υ = 1.5x10^-6 and '_d = 1x10^-6 and use p = 3. The results are shown in Table 3. The method converges in 8 iterations and the errors and condition numbers of the interpolation matrix are reasonably small. The plots for gp^(x), the node distribution at the final iteration and the errors are shown Fig. 2a-c, respectively. We observe that the method adaptively places more interpolation points close to where the function has sharp features. Finally, when interpolating on 141 equally spaced points the error will be 4.9640x10^-5 with the condition number 1.0056x10⁶. Therefore, although the condition number of the adaptive method is a bit larger and the error is smaller. This once again proves the utility of the adaptive method.

CONCLUSION

We have shown that the adaptive radial basis function method provides good results for interpolation of bump functions and can be used for other function with sharp gradients. We also see that the absence of the parameter does not significantly affect our results.

How to cite this article:

Samuel T. Swem and Terhemen Aboiyar. A Parameter-Free Adaptive Method for Radial Basis Function Interpolation of Bump Functions.
DOI: https://doi.org/10.36478/jmmstat.2010.84.88
URL: https://www.makhillpublications.co/view-article/1994-5388/jmmstat.2010.84.88