The Shape Parameter in the Shifted Surface Spline—An Easily Accessible Approach

: In this paper, we present an easily accessible approach to ﬁnding a suitable shape parameter in the shifted surface spline for function interpolation. We aim at helping more readers, including mathematicians and non-mathematicians, to use our method to solve practical problems. Hence, some highly complicated mathematical theorems and deﬁnitions are avoided. The major requirement, as in our previous approach, that the data points should be evenly spaced in the domain is also relaxed. This means that the data points are purely scattered without restrictions. The drawback is that the shape parameter thus obtained is not exactly the same as the theoretically predicted optimal value, which can always be achieved by using our previous rigorous approach. However, experiments show that the gap is quite small and the ﬁnal interpolation errors thus obtained are satisfactory.


Introduction
In the field of radial basis functions (RBF), one of the most frequently used basis functions is the shifted surface spline, defined by if n is even, (1) where x ∈ R n , c > 0, |x| stands for the Euclidean norm of x, log denotes the natural logarithm, m denotes the smallest integer not less than m, and λ, c are constants determined by the user.The latter are called shape parameters and greatly influence the quality of the approximation.From the viewpoint of computation, for even dimensions, λ should be as small as possible.Hence, we will let it be two in this paper whenever possible, especially in the experiment.As for c, its choice is a challenge.If this problem remains unsolved, the power of this radial function will be greatly limited.
The function h(x) defined in Equation ( 1) is obtained from the fundamental solution of the iterated Laplacian by the shifting |x| → (|x| 2 + c 2 ) 1/2 , with c > 0 (Dyn [1] and Yoon [2,3]).Its generalized Fourier transform is ĥ where l(λ, n) is a constant depending on λ, n (see Luh [4]), and Kν (t) := t ν K ν (t), K ν being the modified Bessel function of the second kind (Abramowitz et al. [5]).This Fourier transform will be needed in our main theorem.Moreover, it shows that, from the viewpoint of the Fourier transform, the well-known multiquadrics are simply a type of shifted surface spline, because they share the same Fourier transform as in Equation (2).For odd dimensions, h(x) is generally called the multiquadric.For multiquadrics, the choice of the shape parameter c has been dealt with by the author ( [6]).In this paper, we will focus on even dimensions.Hereafter, whenever referring to Equation (1), we will mean h(x) for even dimensions only.
One of the main advantages of the RBF approach for function interpolation is that it is mesh-free.For any set of scattered data points (x 1 , f (x 1 )), . . ., (x N , f (x N )), where x 1 , . . . ,x N are arbitrary points in R n and f (x 1 ), . . ., f (x N ) are arbitrary real numbers, one can always find a unique interpolator of the form interpolating these data points, where c 1 , . . ., c N are constants to be determined and p is a polynomial of degree ≤ m − 1.The only requirement for the data points is that x 1 , . . ., x N should be polynomially nondegenerate; that is, the only polynomial that vanishes at x 1 , . . ., x N is the zero polynomial.The detailed theory can be seen in Wendland [7] and Buhmann [8].When using Equation (1), we already have a theoretically rigorous and practically useful method to choose the optimal shape parameter c contained in h, as introduced in Luh [4].However, there is a severe restriction in [4] that the interpolation centers x 1 , . . ., x N should be evenly spaced in a simplex.This, to some extent, means that meshes are still needed.Moreover, the function domain must be a simplex.In this paper, we completely discard the two requirements and develop a new approach for finding a suitable c.

Materials and Methods
In order to make this paper self-contained, without overly referring to the complicated theory in [4], next, we recall its main results as in the following subsection.

The Orthodox Approach
Two definitions governing the distribution of the interpolation points are needed.Definition 1.Let E be an n-dimensional simplex in R n with vertices v 1 , . . ., v n+1 (see, e.g., [9]).For any point x ∈ E, its barycentric coordinates are the non-negative numbers λ 1 , . . ., λ n+1 satisfying Now, we need a parameter, called the degree, which indicates the amount of interpolation points in the domain.The larger the degree is, the more points are used.Definition 2. For any natural number k, the corresponding domain points of E of degree k are the evenly spaced points with barycentric coordinates As pointed out in Luh [4], the number of such points is equal to the dimension of P n k , the space of n-variate polynomials of degree less than or equal to k.In this subsection, we use N to denote this number.Thus, N := dimP n k .Two constants involved in the main theorem should be defined here.

Definition 3.
Let h be as in Equation (1).The constants ρ and ∆ 0 are defined as follows.
The approximated functions lie in the so-called native space C h,m .We refer to Luh, Madych and Nelson [10][11][12] and Wendland [7] for its definition.For such functions, we have the following main theorem describing the upper bound for the interpolation error.

Theorem 1.
Let h be as in Equation (1).For any positive number b 0 , there exist positive constants δ 0 , c 1 , C, ω, 0 < ω < 1, completely determined by h and b 0 , such that, for any n-dimensional simplex Q 0 of diameter b 0 , any f ∈ C h,m , and any 0 < δ ≤ δ 0 , there exists a number r satisfying the property that 1 3C ≤ r ≤ b 0 , and, for any n-dimensional simplex Q of diameter r, Q ⊆ Q 0 , there is an interpolating function s(•) as defined in Equation ( 3) such that for all x in Q, where C is defined by where ρ and c appeared in Definition 3 and Equation ( 1), respectively.The function s(•) interpolates f at x 1 , . . ., x N , which are evenly spaced points of degree k − 1 on Q as defined in Definition 2, with k = r δ .Here, f h is the h-norm of f in the native space.The numbers δ 0 , c 1 and ω are given by δ 0 := 1 3C(m+1) , where m appeared in Equation (1), where λ is as in Equation ( 1), l(λ, n) appeared in Equation ( 2), α n is the volume of the unit ball in R n , and ∆ 0 was defined in Definition 3; ω := 2 3 1 3C .

Remark 1.
The key to understanding this result is to grasp the parameter δ.Although δ is not explicitly defined, it is, in concept, similar to the well-known fill distance, which describes the spacing of the data points.In other words, the smaller δ is, the more data points are used in the domain.
The upper bound provided in Formula (4) of Theorem 1 still cannot be used directly to choose a suitable shape parameter c because the dependence of f h on c is not transparent.We must further introduce a function space, E σ , which is a subset of the native space C h,m .Definition 4. For any σ > 0, define where f denotes the classical Fourier transform of f and the integral is the general integral.For each f ∈ E σ , its norm is .
In the preceding definition, For all f ∈ E σ and δ > 0, if the parameter b 0 is fixed, the essential part in the upper bound of Formula (4) depending on c will become where c 0 := 24ρ(m + 1)δ and c 1 := 12ρb 0 .If b 0 is not fixed, we have The optimal value of c minimizes the function MN.Although there is a restriction c ≥ c 0 , it is harmless since c 0 is usually very close to zero.In order to offer the reader an insight into this approach, we present in Figures 1 and 2 two MN curves for Cases 1 and 2 in Luh [4].Note that with the same fill distance δ, the MN function values for n = 4 are larger than those of n = 2.It shows that for higher dimensions, more data points are needed in order to obtain a sharper MN curve.In fact, MN curves are very informative, as we shall see later.Graph of the MN function with ∆ 0.006

The Unorthodox Approach
The reader may be frustrated by the complexity of Theorem 1 and the definition of MN(c).In fact, this is not a serious issue, as will be explained in this subsection.As the inventor of this theory, the author knows that some restrictions can be relaxed, and this will greatly reduce the effort required.In Theorem 1, the parameter b 0 is used simply to control the size of the interpolation domain Q.Hence, one can try to let it be the diameter of the domain only.The shape of the domain is also not important.One can let it be of any shape, e.g., a cube.The parameter δ can be considered to be the widely used fill distance defined by where Q is the domain and X = {x 1 , . . ., x N } is the set of sample points.In fact, one can even drop the restriction of the fill distance and consider that δ is merely a positive parameter inversely proportional to the amount of data points.Then, we completely discard the severe restriction that the sample points be evenly spaced in a simplex, as defined by Definition 2. In other words, the sample points are allowed to be arbitrarily scattered in the domain.The definition of the MN functions remains unchanged as in Equations ( 5) and ( 6).All these form the so-called unorthodox approach.
In this paper, the shape parameter c is chosen by the unorthodox approach.Based on the principles introduced in Section 2.1 and the strategies introduced in this subsection, namely Section 2.2, we have a new set of criteria of choosing c, as introduced in the next two subsections.

b 0 Fixed
For f ∈ E σ , there are two cases.
Reason: In this case, MN(c) may not be monotonic on both [c 0 , c 1 ] and [c 1 , ∞).

b 0 Not Fixed
If the domain is unbounded and hence the sample points can be located in a bounded geometric object of arbitrarily large diameter b 0 , then MN(c) defined by Equation (6) applies.Since such MN(c) is not a monotonic function, one should find c * that minimizes MN(c) for c ∈ [c 0 , ∞) by using appropriate software, such as Mathematica or Matlab.Note that we never decrease b 0 because it will worsen the error bound.

Results
We perform a two-dimensional experiment in the unit square Ω = {(x, y)| 0 ≤ x ≤ 1, 0 ≤ y ≤ 1}.It is shown in Madych [13] that a function f (x) = e − c x 2 belongs to E σ if σ > 2 c.The parameter σ greatly affects the behavior of the MN function MN(c).Here, we choose σ = 1/10 and c = 1/21.Hence, the approximated function is f (x) = e − x 2 /21 for x ∈ Ω.The radial function defined in Equation ( 1) is chosen to have λ = 2. Thus, m = 1 + λ/2 = 2 and the interpolating function s is of the form Equation ( 3), where p is a two-dimensional polynomial of degree m − 1 = 1.Then, we adopt the unorthodox approach to find a suitable shape parameter c.
The first step is to analyze the MN curves.The constant ρ appearing in Equation ( 5) is only one by Definition 3, and the domain diameter b 0 is √ 2. With these understandings, we can now present the needed MN curves, as in Figures 3-7.MN curve for ¡=0.02MN curve for ¡=0.01These figures clearly show that as δ decreases, the optimal values of c move to 12b 0 = 16.97 ≈ 17 and are fixed there.It strongly suggests that one should choose c = 17 to make the interpolation.This method of choosing c is very reliable for the orthodox approach, as we saw in Luh [4], where the theoretically predicted optimal value coincides exactly with the experimentally optimal one.However, for the unorthodox approach, it is expected to lose some accuracy due to the removal of some severe restrictions.
In order to measure the quality of the approximation, we adopt the root-mean-square error defined by where (x i , y i ), i = 1, . . ., M, are test points evenly spaced in the domain Ω.The number M is always chosen to be larger than the number of the randomly spaced data points (interpolating centers) N. In this section, we use N t and N d to denote M and N, respectively.We point out that the randomly spaced points are generated by the Mathematica command RandomPoint.The condition number of the interpolating matrix for establishing s is denoted by COND.As is well known, the condition number for RBF interpolation may be quite large.In order to cope with the problem of ill-conditioning, the arbitrarily precise computer software Mathematica is used.There are always enough effective digits for the calculations.For example, if COND = 1E20, we adopt at least forty effective digits to the right of the decimal point for each step of the calculation.Hence, the results obtained are reliable.As for the crucial parameter δ contained in the MN function, we regard it as simply an indicator of the number of data points and it is inversely proportional to the data amount.When plotting the MN curves, the δs can be arbitrarily chosen, as long as their corresponding curves show clearly the trend of the optimal c values.The experimental results are presented in Tables 1-6, where the optimal value of c is marked by the symbol * .In this experiment, the most time-consuming step in solving the linear system for each c value took less than one second, even though forty effective digits were adopted for each calculation.

Discussion
We emphasize that our theory is reliable only when enough data points are used, and the choice of the shape parameter is determined by a series of MN curves, not only one, as stated in the third paragraph of Section 3.Moreover, our prediction may not be so accurate as in the orthodox rigorous approach of Luh [4].These two traits can be clearly seen in these tables.However, whenever enough data points are used, our predictions are always near the experimentally optimal values and the RMSs thus obtained are satisfactory.
Note that as the number of data points increases, the RMS does seem to improve whenever the optimal c is chosen.However, if c is fixed, say c = 120 or 210, the RMS does not improve as N d increases.This shows that increasing the number of data points is meaningful only when the shape parameter is well chosen.It may challenge the traditional concept of the convergence order based on the number of data points.

Conclusions
As a final conclusion, we believe that the approach presented in this paper is feasible and has greatly reduced the effort required in using the highly complicated Theorem 1 needed for the orthodox approach presented in [4].For pure mathematicians, who insist on absolute theoretical rigor, the approach in [4] is recommended.For applied mathematicians and non-mathematicians, who emphasize practical usefulness, we strongly suggest the method presented in this paper because it is much simpler and more useful.