Elements of bi-cubic polynomial natural spline interpolation for scattered data: Boundary conditions meet partition of unity technique

This paper investigates one kind of interpolation for scattered data by bi-cubic polynomial natural spline, in which the integral of square of partial derivative of two orders to x and to y for the interpolating function is minimal (with natural boundary conditions). Firstly, bi-cubic polynomial natural spline interpolations with four kinds of boundary conditions are studied. By the spline function methods of Hilbert space, their solutions are constructed as the sum of bilinear polynomials and piecewise bi-cubic polynomials. Some properties of the solutions are also studied. In fact, bi-cubic natural spline interpolation on a rectangular domain is a generalization of the cubic natural spline interpolation on an interval. Secondly, based on bi-cubic polynomial natural spline interpolations of four kinds of boundary conditions, and using partition of unity technique, a Partition of Unity Interpolation Element Method (PUIEM) for fitting scattered data is proposed. Numerical experiments show that the PUIEM is adaptive and outperforms state-of-the-art competitions, such as the thin plate spline interpolation and the bi-cubic polynomial natural spline interpolations for scattered data.


Introduction
In the classical problem of variational spline interpolation one looks for a C 2 curves that satisfies the problem minimize{ ∫ [a,b] |s ′′ (t)| 2 dt : s ∈ S, s(t i ) = y i , i = 1, 2, · · · , n}, where τ := {a = t 1 < · · · < t n = b} is a partition of interval [a, b], and S := {s(t) : s ′′ ∈ L 2 [a, b], s α (t) are absolutely constinuous on [a, b], α = 0, 1} is a functional space on [a, b]. It is well known that the solution to the above problem is a cubic natural spline that approximate the thin-beam spline of mechanics [16]. The above problem is also called cubic natural spline interpolation for one dimension scattered data [1,2]. It is also well known that to compute the spline, one solves a linear system of the order of the number of scattered data points within [a, b] that are being interpolated. The global solution has many nice approximate properties. Though multivariate tensor product B-splines interpolants on gridded data play a significant role in applications, they do not work on scattered data [3,8,4,6,7]. The tensor product B-spline representations increase the dimension of the spline space and can not deal directly with multivariate scattered data interpolation. How to construct directly multivariate splines interpolatory basis for scattered data has been yet a challenging work
We denote Z = ℜ N the N-dimensional Euclidean space, and we will mark any element in the space by z = (z 1 , z 2 , · · · , z N ). Let A : X → Z be the interpolatory operator defined by Au = (u(x 1 , y 1 ), · · · , u(x N , y N )).
We study the spline interpolation problem for scattered data in the Hilbert space H 2,2 (R), which is also called bi-cubic natural spline interpolation problem. Its solution is named bi-cubic natural spline. Problem: Look for a function σ(x, y) ∈ X such that [12] σ = arg min{∥T u∥ 2 : Au = z, u ∈ X, z = (z 1 , z 2 , · · · , z N )} (1) where ∥T u∥ 2 = ∫∫ R (u (2,2) (x, y)) 2 dxdy with natural boundary conditions as follows: In practice, there are four kinds of boundary conditions, which are ac, ad, bc and bd as above. In this paper they are named as type I, II, III, and IV respectively. 996 ELEMENTS OF BI-CUBIC POLYNOMIAL NATURAL SPLINE INTERPOLATION FOR SCATTERED DATA Denote by N (T ) and N (A) the null spaces of the operator T and A correspondingly, and by R(T ), R(A) their ranges. Let A z := {u ∈ X|Au = z, u satisf ies type I natural boundary conditions} be a set of u ∈ X satisfying both interpolatory conditions and type I natural boundary conditions. We consider the existence and uniqueness of the interpolatory solution with type I natural boundary conditions. It is proved easily the theory as follows:

Theorem 1
The null subspace of the operator T is The boundary conditions of Theory 2 in [13] is marked by type I, and Theory 2 in [13] is marked by Lemma 1.
Theorem 3 [12] If then the solution of bi-cubic natural spline interpolation for scattered data with the boundary conditions of type I does exist and is unique, and its explicit and compactly format expression as follows where c kl (k = 0, 1; l = 0, 1) and λ i (i = 1, · · · , N ) are real numbers obtained by the interpolatory conditions. Here g i (x, y) = G(x i , a; x)G(y i , c; y)(i = 1, 2, · · · , N ) are bi-cubic natural spline basis where

Theorem 4
If N (T ) ∩ N (A) = {Θ X }, then the solution of bi-cubic natural spline interpolation for scattered data with the boundary conditions of type IV does exist and is unique, and its explicit and compactly format expression as follows where c kl (k = 0, 1; l = 0, 1) and λ i (i = 1, · · · , N ) are real numbers obtained by the interpolatory conditions. Here Proof From the theorem 2 and 3, the solution of the bi-cubic natural spline interpolation with the boundary conditions of type IV does exist and is unique. By the spline theory of Hilbert spaces, it is well known that if σ(x, y) is the interpolatory spline function solution in Hilbert spaces, then for all v ∈ N (T ), we have ⟨T * T σ, v⟩ = 0, and there exist coefficients λ i and k i such that For all u ∈ X, by using Taylor Theorem at x = b firstly, and by using Taylor Theorem at y = d secondly, and from boundary conditions (2) and (4), it is easy to find that andG(x, y) satisfies the boundary conditions (2), thus So we have where g i (x, y) = G 1 (x i , b; x)G 1 (y i , d; y)(i = 1, 2, · · · , N ) are bi-cubic natural spline basis and The coefficients c j (r i , t)(j = 0, 1) are obtained by the boundary conditions (5).
Because g Similarly, G 1 (y i , d; y) has the same expression as G 1 (x i , b; x). So, we also know that all g i (x, y) for i = 1, 2, · · · , N are bi-cubic natural spline basis and Similarly, as for the boundary conditions of type II there exists the property as follows: then the solution of bi-cubic natural spline interpolation for scattered data with the boundary conditions of type II does exist and is unique, and its explicit and compactly format expression as follows where c kl (k = 0, 1; l = 0, 1) and λ i (i = 1, · · · , N ) are real numbers obtained by the interpolatory conditions.
are bi-cubic natural spline basis with type II boundary conditions.
Proof It is similar to the proof of the previous theorem. By using Taylor Theorem at x = a for all u ∈ X firstly, and secondly by using Taylor Theorem at y = d for u(a, y), u (1,0) (a, y) and u (2,0) (τ, y). Similarly, as for the boundary conditions of type III there exists the property as follows: then the solution of bi-cubic natural spline interpolation for scattered data with the boundary conditions of type III does exist and is unique, and its explicit and compactly format expression as follows where c kl (k = 0, 1; l = 0, 1) and λ i (i = 1, · · · , N ) are real numbers obtained by the interpolatory conditions. Here Proof It is similar to the proof of the previous theorem. By using Taylor Theorem at x = b for all u ∈ X firstly, and secondly by using Taylor Theorem at y = c for u(b, y), u (1,0) (b, y) and u (2,0) (τ, y).

Resolution of Bi-cubic Natural Spline Interpolation for Scattered Data
Now, for all given scattered data set {(x i , y i , z i )|i = 1, 2, · · · , N }, we can fit the scattered data by using bi-cubic natural spline interpolation function as By the interpolatory conditions, λ i (i = 1, · · · , N ) and c kl (k = 0, 1; l = 0, 1) can be decided by a linear system as follows where So, they include the interpolation conditions of type I, II, III, and type IV.

Theorem 7
The coefficient matrix of the linear system (7), which come from bi-cubic natural spline interpolation for scattered data, is symmetric.
Proof We only prove type II, and the others are similar. It is only needed to prove that A is symmetric. In order to do this, we denote Thus, a ij = α ij β ij .
In the general, suppose x i x j , then there is Suppose y j y i , then there is So, a ij = a ji . This means that A is symmetric.
The interpolatory function is achieved with the values of the coefficients λ i and c ij obtained by solving the linear system (7).

Properties of the Bi-cubic Natural Spline Interpolation for Scattered Data
For all given scattered data set {(x i , y i )|i = 1, 2, · · · , N }, project them to axis x and axis y and get a partition:

Theorem 8
The type I interpolation for scattered data is a bi-linear polynomial on , is a linear polynomial to y and a cubic polynomial to and is a linear polynomial to x and a cubic polynomial to y on Proof The proof will use many times this property: It is a bi-linear polynomial to x and y.
It is also a bi-linear polynomial to x and y.
It is a linear polynomial to y and a cubic polynomial to x.
Stat., Optim. Inf. Comput. Vol. 8, December 2020 It is a linear polynomial to x and a cubic polynomial to y.
We can directly verify that there exist the property as follow:

Theorem 9
The type I interpolation for scattered data is a linear polynomial to x and a cubic polynomial to y on λ i x µ i y ν i = 0, µ = 0, 1, ν = 0, 1} be a set of natural bi-cubic spline, and A z := {u ∈ X|Au = z, u satisf ies type I natural boundary conditions} be a set of u ∈ X satisfying both interpolatory conditions and type I natural boundary conditions.

Lemma 2
Let u ∈ X satisfy type I natural boundary conditions, and let σ(x, y) ∈ N S. Then Proof Using integration by parts with respect to x, we have From boundary conditions u (0,2) (a, y) = u (1,2) (a, y) = 0 and σ ∈ N S, we have  (4,2) dxdy. Similarly, using integration by parts with respect to y, from type I boundary conditions we have Thus ∥T σ∥ 2 = min u∈Az ∥T u∥ 2 .

This means
From the formula, it follows that ∥T σ∥ 2 = min u∈Az ∥T u∥ 2 .
Theorem 11 (Best Approximating Property) Suppose that u ∈ A z , σ ∈ N S ∩ A z , then for all s ∈ N S, such that If there exists anotherσ satisfies this property, that the difference between σ andσ is only a bi-linear polynomial.
Stat., Optim. Inf. Comput. Vol. 8, December 2020 Proof Using u − s replace u and σ − s replace σ in the Theorem 10, it follows that If there has anotherσ satisfying the property too, that's mean ∥T (u −σ) Then The other types have similar properties, and they also have similar proofs.
In the subsection, we apply the above constructed interpolation elements for scattered data. We fit the function u by the bi-cubic natural spline interpolations for Ex1,2,3,4. Ex1. In order to compare with the respective interpolation effects of four types, we produce 3 sets of simple data: they are sparser near x = b than near x = a in the rectangular domain R. The interpolatory results for max errors and average errors are listed in Table 1 to Table 3, the interpolatory surfaces are plotted in Figure 1 to Figure 3.
The experimental results show that the effects of scattered data fitting in type III and in type IV are slightly better than the ones in type I and type II.   Ex2. In order to compare with the respective interpolation effects of four types, we produce 2 sets of simple data: they are sparser near the corner x = b, y = d than the others in the rectangular domain R. The interpolatory results for max errors and average errors are listed in Table 4 and Table 5, the interpolatory surfaces are plotted in Figure  4 and Figure 5. The experimental results show that the effects of scattered data fitting in type IV are slightly better than the ones in other types.
From Ex1 and Ex2, the interpolation errors for scattered data mainly include internal error and boundary error. The latter can be reduced by choosing the respective interpolation type when the data in one corner is very sparse. Moreover, for the general scattered data, the Partition of Unity Interpolation Element Method (PUIEM) is proposed.

Partition of Unity Interpolation Element Method (PUIEM)
We test our adaptivity on the partition of unity interpolation element method (PUIEM) for scattered data. The basic idea of the partition of unity method consists of decomposing the original domain into several subdomains or patches forming a covering of it, then constructing a local interpolation element approximant on each subdomain. Such an approach, as evident in our numerical analysis, enables us to significantly reduce the computational   instability due to severe ill-conditioning generated by interpolation element methods. Originally, the partition of unity finite element method was introduced in [14,15] to solve PDEs. Recently, this method has reached popularity because it allows one to efficiently and accurately solve large scale interpolation problems, as well as differential equation ones. However, the problem of constructing an adaptive method is still open and does not seem to have received the proper attention in literature.
Partitioning R into k patches ω j such that ∪ k j=1 ω j ⊇ R with mild overlap among the patches, we can define the PUIEM by considering a partition of unity {α j } k j=1 subordinated to the covering {ω j } k j=1 such that The weight α j : ω j → R is a compactly supported, nonnegative function with supp(α j ) ⊆ ω j . So the global PUIEM approximant assumes the form The approximant f j is chosen from the bi-cubic natural spline interpolation elements with different boundary conditions. Ex4. In order to compare with the respective interpolation effects of four types and PUIEM, we generate 621 random data in the rectangular domain R. For similarity, we consider the four small rectangular α j , j = 1, 2, 3, 4 which cover R as shown in Figure 6. f j (x, y), j = 1, 2, 3, 4 are the bicubic natural spline inpolation elements with corresponding boundary conditions. The interpolatory results for max errors and average errors are listed in Table 6, the error surfaces of interpolation by four types and PUIEM are plotted in Figure 7. The experimental results show that the effects of scattered data fitting in PUIEM are better than four types. Ex5. In order to compare with the PUIEM and the thin plate spline (φ(r) = r 2 ln |r|, r = √ (x − x i ) 2 + (y − y i ) 2 ) interpolation, we produce 2 sets of simple data: the gridded data and the scattered data. The latter are produced by random functions belong to R = [0, 1] × [0, 1]. Denote sd = N p×q the measure to scattered data, then 1 N ≤ sd ≤ 1.  The interpolatory results for max errors and average errors are listed in Table 6 and Table 7, the interpolatory surfaces are plotted in Figure 7 and Figure 8.    Figure 9 show that the effects of scattered data fitting in the PUIEM are better than the ones in the thin plate spline.
In conclusion, the experimental results show that Our PUIEM is adaptive and superior to the four natural spline interpolation types and the thin plate spline in the interpolation results.

Conclusion
In this paper, One kind of interpolation for scattered data by bi-cubic polynomial natural spline with different boundary conditions is studied. By the spline function methods of Hilbert space, their solutions are constructed as the sum of bi-linear polynomials and piecewise bi-cubic polynomials. Many properties of the solutions are studied and it is concluded that bi-cubic natural spline interpolation on a rectangular domain is a generalization of the cubic natural spline interpolation on an interval. Moreover, Based on partition of unity technique and bi-cubic polynomial natural spline interpolations of four kinds of boundary conditions, Partition of Unity Interpolation Element Method (PUIEM) for fitting scattered data is proposed. The experiments show that the PUIEM is adaptive and outperforms the bi-cubic polynomial natural spline interpolations of four kinds and the thin plate spline in the interpolation results for scattered data.
In a general way, the effects of scattered data fitting depend mainly on their distribution (e.g., sparsity and irregular boundaries). When the data in one corner is very sparse, it is suggested to choose the PUIEM or the respective bicubic natural spline interpolation type to minimize the error. Of course, it is also suggested that to choose more interpolatory points uniformly for data collection. In the future, the refined partition of unity will be considered and the large scale real scattered data (e.g., geographical data [18]) will be tested.