Optimal selection of local approximants in RBF-PU interpolation

TL;DR

This paper introduces a method for optimally selecting local data set sizes and shape parameters in RBF-PU interpolation, enhancing accuracy for large scattered data problems with non-uniform density.

Contribution

It proposes a novel approach to adaptively choose local data points and shape parameters by minimizing an a priori error estimate, improving interpolation accuracy.

Findings

High accuracy achieved with non-homogeneous data density

Effective selection of local data set size and shape parameter

Numerical experiments validate the method's robustness

Abstract

The Partition of Unity (PU) method, performed with local Radial Basis Function (RBF) approximants, has been proved to be an effective tool for solving large scattered data interpolation problems. However, in order to achieve a good accuracy, the question about how many points we have to consider on each local subdomain, i.e. how large can be the local data sets, needs to be answered. Moreover, it is well-known that also the shape parameter affects the accuracy of the local RBF approximants and, as a consequence, of the PU interpolant. Thus here, both the shape parameter used to fit the local problems and the size of the associated linear systems are supposed to vary among the subdomains. They are selected by minimizing an a priori error estimate. As evident from extensive numerical experiments and applications provided in the paper, the proposed method turns out to be extremely accurate also when data with non-homogeneous density are considered.

Tables

Table 1. Table 1: CPU times (in seconds) obtained by running the sorting-based procedure ( t S − P ​ S subscript 𝑡 𝑆 𝑃 𝑆 t_{S-PS} ) and the integer-based one ( t I − P ​ S subscript 𝑡 𝐼 𝑃 𝑆 t_{I-PS} ).

$N$	$25000$	$50000$	$100000$	$200000$
$t_{I-PS}$	$5.13$	$10.68$	$21.99$	$45.00$
$t_{S-PS}$	$5.21$	$12.40$	$28.77$	$71.55$

Table 2. Table 2: RMSEs and MAEs computed on Halton points and obtained by using the IMQ as local RBF interpolant for f 1 subscript 𝑓 1 f_{1} .

$N$	method	RMSE	MAE
$289$	PU	$3.64E-03$	$5.66E-02$
	BLOOCV-PU	$1.03E-05$	$2.36E-04$
$1089$	PU	$7.57E-04$	$1.52E-02$
	BLOOCV-PU	$2.88E-06$	$7.89E-05$
$4225$	PU	$3.88E-04$	$1.01E-02$
	BLOOCV-PU	$3.84E-07$	$1.39E-05$
$16641$	PU	$8.27E-04$	$3.27E-02$
	BLOOCV-PU	$9.67E-08$	$3.15E-06$
$66049$	PU	$1.08E-05$	$1.09E-04$
	BLOOCV-PU	$2.68E-08$	$6.80E-07$

Table 3. Table 3: RMSEs and MAEs computed on Halton points and obtained by using the IMQ as local RBF interpolant for f 2 subscript 𝑓 2 f_{2} .

$N$	method	RMSE	MAE
$289$	PU	$2.59E-02$	$4.30E-01$
	BLOOCV-PU	$1.32E-02$	$2.76E-01$
$1089$	PU	$3.51E-03$	$6.20E-02$
	BLOOCV-PU	$2.11E-04$	$8.93E-03$
$4225$	PU	$8.63E-04$	$2.00E-02$
	BLOOCV-PU	$3.88E-06$	$1.12E-04$
$16641$	PU	$4.07E-04$	$1.18E-02$
	BLOOCV-PU	$8.26E-08$	$2.80E-06$
$66049$	PU	$1.23E-04$	$4.19E-03$
	BLOOCV-PU	$5.10E-08$	$1.76E-06$

Table 4. Table 4: RMSEs and MAEs computed on Halton points and obtained by applying the BLOOCV-PU method with the IMQ as local RBF interpolant and optimal ε j subscript 𝜀 𝑗 \varepsilon_{j} (but δ j = δ subscript 𝛿 𝑗 𝛿 \delta_{j}=\delta ) for f 1 subscript 𝑓 1 f_{1} .

N	RMSE	MAE
$289$	$3.00E-03$	$3.35E-02$
$1089$	$8.88E-04$	$1.25E-02$
$4225$	$2.48E-04$	$7.50E-03$
$16641$	$1.11E-04$	$2.49E-03$
$66049$	$8.64E-06$	$1.61E-04$

Table 5. Table 5: Fill and separation distances of different sets of Halton data and non-conformal points.

$N$	data set	$h_{{𝒳}_N,\Omega }$	$q_{{𝒳}_N}$	$h_{{𝒳}_N,\Omega }$ / $q_{{𝒳}_N}$
$289$	Halton	$7.46E-02$	$1.03E-02$	$7.20E+00$
	non-conformal	$1.87E-01$	$4.90E-03$	$3.81E+01$
$1089$	Halton	$3.93E-02$	$4.33E-03$	$9.07E+00$
	non-conformal	$1.39E-01$	$1.31E-03$	$1.06E+02$
$4225$	Halton	$2.19E-02$	$2.19E-03$	$9.97E+00$
	non-conformal	$9.94E-02$	$6.84E-04$	$1.45E+02$

Table 6. Table 6: RMSEs and MAEs computed on non-conformal points and obtained by using the Wendland C 6 superscript 𝐶 6 C^{6} as local RBF interpolant for f 1 subscript 𝑓 1 f_{1} .

$N$	method	RMSE	MAE
$289$	PU	$3.28E-02$	$1.90E-01$
	BLOOCV-PU	$3.64E-03$	$4.15E-02$
$1089$	PU	$1.12E-02$	$2.01E-01$
	BLOOCV-PU	$5.40E-04$	$9.11E-03$
$4225$	PU	$1.44E-02$	$2.23E-01$
	BLOOCV-PU	$1.24E-04$	$3.34E-03$
$16641$	PU	$1.12E-02$	$1.92E-01$
	BLOOCV-PU	$3.21E-05$	$7.05E-04$
$66049$	PU	$1.26E-02$	$2.45E-01$
	BLOOCV-PU	$1.14E-05$	$3.70E-04$

Table 7. Table 7: RMSEs and MAEs computed on non-conformal points and obtained by using the Wendland C 6 superscript 𝐶 6 C^{6} as local RBF interpolant for f 2 subscript 𝑓 2 f_{2} .

$N$	method	RMSE	MAE
$289$	PU	$5.30E-02$	$5.00E-01$
	BLOOCV-PU	$3.47E-02$	$3.13E-01$
$1089$	PU	$3.90E-02$	$5.00E-01$
	BLOOCV-PU	$7.11E-03$	$8.38E-02$
$4225$	PU	$4.63E-02$	$4.99E-01$
	BLOOCV-PU	$2.39E-03$	$4.77E-02$
$16641$	PU	$4.34E-02$	$5.00E-01$
	BLOOCV-PU	$7.51E-04$	$8.10E-03$
$66049$	PU	$4.03E-02$	$5.00E-01$
	BLOOCV-PU	$8.28E-05$	$1.27E-03$

Table 8. Table 8: RMSEs and MAEs computed on the glacier and black forest data sets and obtained by using the BLOOCV-PU method with the Matérn C 2 superscript 𝐶 2 C^{2} as local RBF interpolant.

Data set	RMSE	MAE
Glacier	$0.65$ m	$3.31$ m
Forest	$5.73$ m	$26.0$ m

Equations

$$\sum _{j=1}^d{W}_j(\mathbf{x})=1,\mathbf{x}\in \Omega ,$$

$${\Vert D^{\mathbf{\beta }}{W}_j\Vert }_{L^{\infty }({\Omega }_j)}\le \frac{C_{\mathbf{\beta }}}{{\left({\sup }_{\mathbf{x},\mathbf{y}\in {\Omega }_j}{\Vert \mathbf{x}-\mathbf{y}\Vert }_2\right)}^{|\mathbf{\beta }|}},j=1,\dots ,d.$$

$$W_j(\mathbf{x})=\frac{{\bar{W}}_j(\mathbf{x})}{\sum _{k=1}^d{\bar{W}}_k(\mathbf{x})},j=1,\dots ,d,$$

$$\mathcal{I}\left(\mathbf{x}\right)=\sum _{j=1}^d{R}_j\left(\mathbf{x}\right)W_j\left(\mathbf{x}\right),\mathbf{x}\in \Omega .$$

$$R_j(\mathbf{x})=\sum _{k=1}^{N_j}{c}_k^j\varphi (||\mathbf{x}-{\mathbf{x}}_k^j|{|}_2),$$

$$R_j({\mathbf{x}}_i^j)=f_i^j,i=1,\dots ,N_j,$$

$$A_j{\mathbf{c}}_j={\mathbf{f}}_j,$$

$$(A_j{)}_{ik}=\varphi (||{\mathbf{x}}_i^j-{\mathbf{x}}_k^j|{|}_2),i,k=1,\dots ,N_j.$$

$$l=\left(\begin{array}{c}L-1+M\\ L-1\end{array}\right).$$

$$\varphi (r)={\left(1+(\epsilon r{)}^2\right)}^{-1/2},$$

$$\varphi (r)=e^{-\epsilon r}(1+\epsilon r).$$

$$\varphi (r)={\left(1-\epsilon r\right)}_{+}^4\left(4\epsilon r+1\right),$$

$$\varphi (r)={\left(1-\epsilon r\right)}_{+}^8\left(32(\epsilon r{)}^3+25(\epsilon r{)}^2+8\epsilon r+1\right),$$

$$q_{{\mathcal{X}}_N}=\frac{1}{2}\underset{i\ne k}{\min }{\Vert {\mathbf{x}}_i-{\mathbf{x}}_k\Vert }_2.$$

$$h_{{\mathcal{X}}_N,\Omega }=\underset{\mathbf{x}\in \Omega }{\sup }\left(\underset{{\mathbf{x}}_k\in {\mathcal{X}}_N}{\min }{\Vert \mathbf{x}-{\mathbf{x}}_k\Vert }_2\right).$$

$$q_{{\mathcal{X}}_N}\le h_{{\mathcal{X}}_N,\Omega }\le C_{qu}{q}_{{\mathcal{X}}_N}.$$

$$C(\mathbf{x},\mathbf{\xi }(\mathbf{x}),\theta ,\gamma )=\{\mathbf{x}+\lambda \mathbf{y}:\mathbf{y}\in {\mathbb{R}}^M,||\mathbf{y}|{|}_2=1,{\mathbf{y}}^T\mathbf{\xi }(\mathbf{x})\ge cos(\theta ),\lambda \in [0,\gamma ]\},$$

$$|D^{\beta }f(\mathbf{x})-D^{\beta }\mathcal{I}(\mathbf{x})|\le C^{{}^{\prime }}{h}_{{\mathcal{X}}_N,\Omega }^{(k+\nu )/2-|\beta |}|f{|}_{{\mathcal{N}}_{\varphi }(\Omega )},$$

$$\left|f(\mathbf{x})-\mathcal{I}(\mathbf{x})\right|\le \sum _{j=1}^d\left|f_j(\mathbf{x})-R_j(\mathbf{x})\right|W_j(\mathbf{x})\le \underset{j=1,\dots ,d}{\max }{\Vert f_j-R_j\Vert }_{L_{\infty }({\Omega }_j)},$$

$$R_j^{(i)}(\mathbf{x})=\sum _{k=1,k\ne i}^{N_j}{c}_k^j\varphi (||\mathbf{x}-{\mathbf{x}}_k^j|{|}_2),$$

$$e_i^j=f_i^j-R_j^{(i)}({\mathbf{x}}_i^j),$$

$$e_i^j=\frac{c_i^j}{{\left(A_j\right)}_{ii}^{-1}},$$

$$\left(e_1^j,\dots ,e_{N_j}^j\right)=\left(\frac{c_1^j}{{\left(A_j\right)}_{11}^{-1}},\dots ,\frac{c_{N_j}^j}{{\left(A_j\right)}_{N_j{N}_j}^{-1}}\right).$$

$${\mathbf{e}}_j\left({\delta }_{j_p},{\epsilon }_{j_q}\right)=\left(e_1^j\left({\delta }_{j_p},{\epsilon }_{j_q}\right),\dots ,e_{N_j}^j\left({\delta }_{j_p},{\epsilon }_{j_q}\right)\right).$$

$$E_j=\left(\begin{array}{ccc}||{\mathbf{e}}_j({\delta }_{j_1},{\epsilon }_{j_1})|{|}_{\infty }& \cdots & ||{\mathbf{e}}_j({\delta }_{j_1},{\epsilon }_{j_Q})|{|}_{\infty }\\ ⋮& \ddots & ⋮\\ ||{\mathbf{e}}_j({\delta }_{j_P},{\epsilon }_{j_1})|{|}_{\infty }& \cdots & ||{\mathbf{e}}_j({\delta }_{j_P},{\epsilon }_{j_Q})|{|}_{\infty }\end{array}\right).$$

$$||{\mathbf{e}}_j({\delta }_j,{\epsilon }_j)|{|}_{\infty }=\underset{p=1,\dots ,P}{\min }\left(\underset{q=1,\dots ,Q}{\min }(E_j{)}_{pq}\right).$$

$$\tilde{\mathcal{I}}(\mathbf{x})=\sum _{j=1}^d{\tilde{R}}_j(\mathbf{x})W_j(\mathbf{x}),\mathbf{x}\in \Omega ,$$

$${\tilde{R}}_j(\mathbf{x})=\sum _{k=1}^{{\tilde{N}}_j}{c}_k^j{\varphi }_{{\epsilon }_j}(||\mathbf{x}-{\mathbf{x}}_k^j|{|}_2),$$

$${\mathcal{R}}_M=\prod _{m=1}^M\left[\underset{i=1,\dots ,N}{\min }{x}_{im},\underset{i=1,\dots ,N}{\max }{x}_{im}\right],$$

$$l_{box}=\underset{m=1,\dots ,M}{\max }\left(\underset{i=1,\dots ,N}{\max }{x}_{im}\right)-\underset{m=1,\dots ,M}{\min }\left(\underset{i=1,\dots ,N}{\min }{x}_{im}\right).$$

Full text

Optimal selection of local approximants in RBF-PU interpolation

Roberto Cavorettoa, Alessandra De Rossia, Emma Perracchioneb

a Department of Mathematics “G. Peano”, University of Torino, via Carlo Alberto 10, I–10123 Torino, Italy

b Department of Mathematics “T. Levi-Civita”, University of Padova, via Trieste 23, I–35121 Padova, Italy

Abstract. The Partition of Unity (PU) method, performed with local Radial Basis Function (RBF) approximants, has been proved to be an effective tool for solving large scattered data interpolation problems. However, in order to achieve a good accuracy, the question about how many points we have to consider on each local subdomain, i.e. how large can be the local data sets, needs to be answered. Moreover, it is well-known that also the shape parameter affects the accuracy of the local RBF approximants and, as a consequence, of the PU interpolant. Thus here, both the shape parameter used to fit the local problems and the size of the associated linear systems are supposed to vary among the subdomains. They are selected by minimizing an a priori error estimate. As evident from extensive numerical experiments and applications provided in the paper, the proposed method turns out to be extremely accurate also when data with non-homogeneous density are considered.

1 Introduction

Given a set of multivariate data, we aim at finding a function that accurately fits such points. This problem is rather common in many applied sciences, such as in physics, biology, geophysics and Earth’s topography. Moreover, dealing with applications, one often faces the problem of approximating large and irregular data sets, i.e. data which are far from being uniform or quasi-uniform. In these cases, since problems as lack of information, i.e. low density of points, or ill-conditioning, i.e. high density of data, arise, the fitting process becomes a challenging computational issue.

Because of the above mentioned problems, recently, the approximation theory has driven its attention on local techniques. Specifically, the approximation of irregularly distributed data via local schemes has gained much attention in both meshfree and mesh-dependent methods. For the latter, the problem results particularly hard and the choice of the mesh turns out to be crucial. As example, in [10], in order to build the local approximating fits, the authors consider spline functions on a uniform triangulation with $C^{1}$ or $C^{2}$ continuity. While, in [4], polyharmonic splines are effectively used to fit irregular and truly large data sets.

Another major class of techniques, which are known in literature as meshfree or meshless methods, includes RBF approximants [5, 25]. They obviously take advantage of being independent from the mesh and, as a consequence, they are easy to implement in any dimension. Indeed, a local hybrid approximation technique for data with non-homogeneous density, obtained by means of both splines and RBFs, is presented in [11]. But, since bivariate spline functions are involved, the method again strongly depends on the mesh. To avoid this drawback, we focus on purely meshless methods. In this context, the scattered data problem of huge and irregular sets of points is usually performed by means of least squares approximation [31, 36]. Here instead, our fitting criterion consists in exactly matching the measurements at their corresponding locations, i.e. we focus on interpolation. To this aim, the PU method performed with local RBF interpolants turns out to be really meaningful [3, 17].

The basic idea of the PU technique consists in decomposing the domain into several subdomains or patches forming a covering of the original domain. When the PU method is applied in the context of interpolation, except for particular cases [30], such subdomains are always supposed to be hyperspheres of a fixed size [7, 17, 24, 32]. But, in case of irregular data this might lead to inaccurate approximations. In [30] variable subdomains are used for an ad hoc 2D problem in finance; specifically, even if data are not quasi-uniform, they have a precise and well-known structure. Thus, the PU subdomains are constructed following exactly their distribution.

Our aim consists instead in developing a method which enables us to select, independently from the node distribution, suitable sizes of the different PU subdomains. Furthermore, we also take into account the critical choice of the shape parameter of the basis function. In fact, it can greatly influence the accuracy of final fit.

To such scope, we compute subsequent a priori error estimates depending on both the shape parameter and the size of the PU subdomain. Then, for each patch we select the optimal couple of values, i.e. the subdomain size and the shape parameter, used to solve the local interpolation problem. The error estimates are found out via a modified Leave One Out Cross Validation (LOOCV) scheme [18, 22, 29]. More precisely, since our problem depends on two quantities, for each patch we perform a Bivariate LOOCV (BLOOCV). The resulting method, named BLOOCV-PU, turns out to be extremely accurate compared with the classical PU technique. This improvement, in terms of accuracy, becomes particularly meaningful when data with non-homogeneous density are considered.

The complexity of the algorithm is also taken into account. Specifically, the computational issue consisting in finding all the points belonging to a given subdomain is performed with the use of a novel data structure, the so-called Integer-based Partitioning Structure (I-PS). It leads to a considerable saving in terms of computational time with respect to the most advanced searching techniques [7, 8].

Numerical experiments show the good performances of the BLOOCV-PU meth- od in case of quasi-uniform data, and underline the benefits of such a flexible approach with irregular points. Moreover, we investigate two applications with real world data, including a benchmark glacier data set and points with highly varying densities describing a terrain.

The guidelines of the paper are as follows. In Section 2, after briefly reviewing the main theoretical features of the PU method, we introduce the BLOOCV-PU interpolant. The computational aspects of such algorithm and its complexity are described in Sections 3 and 4, respectively. Extensive numerical experiments and applications with real world data, carried out in Sections 5 and 6 respectively, are devoted to test the accuracy and the flexibility of the BLOOCV-PU approximant. Finally, in Section 7, we deal with conclusions and work in progress.

We point out that the Matlab software is made available to the scientific community in a downloadable free package:

http://hdl.handle.net/2318/1559094.

2 Formulation of the BLOOCV-PU interpolant

In Subsection 2.1 we first briefly review the main theoretical aspects concerning the PU interpolation and then in Subsection 2.2 we focus on the local selection of suitable sizes of the patches and shape parameters. In order to achieve such scope, for each local interpolation problem we need to determine an error estimate depending on both the size of the subdomain and the shape parameter of the basis function.

2.1 The partition of unity method

The scattered data interpolation problem consists in recovering a function $f:\Omega\to\mathbb{R}$, $\Omega$ being a bounded set in $\mathbb{R}^{M}$, using a set of samples of $f$ on $N$ distinct data points or nodes ${\cal X}_{N}=\{\boldsymbol{x}_{i},i=1,\ldots,N\}\subset\Omega$, namely $\boldsymbol{f}=(f_{1},\ldots,f_{N})^{T}$, $f_{i}=f(\boldsymbol{x}_{i})$, with $\boldsymbol{x}_{i}\in{\cal X}_{N}$. More precisely, we aim at solving the above mentioned problem for truly large data sets. To this scope, the PU method, first introduced in [3, 26], is a well-suited numerical tool.

The idea behind the PU method is to start with a partition of the open and bounded domain $\Omega$ into $d$ subdomains $\Omega_{j}$, such that $\Omega\subseteq\cup_{j=1}^{d}\Omega_{j}$, with some mild overlap among them [3, 17, 24, 26, 35]. Together with these subdomains, we need a family of compactly supported, non-negative, continuous functions $W_{j}$, $j=1,\ldots,d$, which form a $k$-stable partition of unity, i.e.

[TABLE]

and for every $\boldsymbol{\beta}\in\mathbb{N}^{M}$, with $|\boldsymbol{\beta}|\leq k$, there exists a constant $C_{\boldsymbol{\beta}}>0$ such that

[TABLE]

More precisely, in what follows we consider the so-called Shepard’s weights which are defined as

[TABLE]

where $\bar{W}_{j}$ are compactly supported functions forming a partition of unity.

Once we choose the partition of unity $\{W_{j}\}_{j=1}^{d}$, the global interpolant is formed by the weighted sum of $d$ local approximants $R_{j}$, i.e.

[TABLE]

In particular here $R_{j}$ denotes a RBF interpolant defined on a subdomain $\Omega_{j}$ of the form

[TABLE]

where $\phi:[0,\infty)\rightarrow\mathbb{R}$ is called RBF, $||\cdot||_{2}$ denotes the Euclidean norm, $N_{j}$ indicates the number of data points belonging to $\Omega_{j}$ and $\boldsymbol{x}_{k}^{j}\in{\cal X}_{N_{j}}={\cal X}_{N}\cap\Omega_{j}$, with $k=1,\ldots,N_{j}$. We observe that if the local fits $R_{j}$, $j=1,\ldots,d$, satisfy the interpolation conditions then the global PU approximant inherits the interpolation property. This trivially follows from the fact that the functions $W_{j}$, $j=1,\ldots,d$, form a partition of unity.

The coefficients $\{c_{k}^{j}\}_{k=1}^{N_{j}}$ in (2) are determined by imposing the interpolation conditions

[TABLE]

which lead to the problem of solving $d$ linear systems of the form

[TABLE]

where $\boldsymbol{c}_{j}=(c_{1}^{j},\ldots,c_{N_{j}}^{j})^{T}$, $\boldsymbol{f}_{j}=(f_{1}^{j},\ldots,f_{N_{j}}^{j})^{T}$ and $A_{j}$ is the local interpolation matrix whose entries are given by

[TABLE]

The system (4) admits a unique solution if and only if the matrix $A_{j}$ is non-singular. Micchelli gave conditions on $\phi$ which guarantee the non-singularity of $A_{j}$ [27]. In particular, these conditions are fulfilled if $\phi$ is a strictly positive definite RBF. Nevertheless, even if here for simplicity we only consider strictly positive definite functions, we point out that such conditions are more general. Precisely, if $\phi$ is strictly conditionally positive definite of order $L$, a unique solution to the interpolation problem is ensured by adding to the expansion (2) certain polynomials which form a basis for the $l$-dimensional space $\Pi_{L-1}^{M}$ of polynomials of total degree less than or equal to $L-1$ in $M$ variables, where

[TABLE]

Among a large variety of known RBFs, which are defined in function of a shape parameter, we restrict our attention on strictly positive definite functions. Furthermore, we can distinguish between compactly and globally supported RBFs. As examples of these two classes we consider the compactly supported Wendland $C^{2}$ and $C^{6}$ functions and the globally defined Matérn $C^{2}$ and Inverse MultiQuadric (IMQ) functions [17]. The latter is infinitely smooth and its analytic expression is

[TABLE]

where $r$ is the Euclidean norm and $\varepsilon$ is a positive shape parameter governing the flatness of the RBF. The family of Matérn functions is instead characterized by a finite regularity. As example, the Matérn $C^{2}$ is defined as

[TABLE]

Concerning Compactly Supported RBFs (CSRBFs), a well-known class of functions which are smooth, locally supported and strictly positive definite is the family of Wendland’s functions. For instance, the Wendland $C^{2}$ and $C^{6}$ functions are respectively defined as

[TABLE]

[TABLE]

where $(\cdot)_{+}$ denotes the truncated power function.

Now, in order to give error bounds, we define two common indicators of data regularity:

Definition 2.1**.**

The separation distance is given by

[TABLE]

The quantity $q_{{\cal X}_{N}}$ represents the radius of the largest ball that can be centered at every point in ${\cal X}_{N}$ such that no two balls overlap.

Definition 2.2**.**

The fill distance, which is a measure of data distribution, is given by

[TABLE]

Remark 2.1**.**

The distances (10) and (11) give an idea of the node distribution, i.e. how uniform data are. Indeed, a set of data is supposed to be quasi-uniform with respect to a constant $C_{qu}$ if

[TABLE]

More specifically, the definition of quasi-uniform points has to be seen in the context of more than one data set. The idea is to consider a sequence of such sets so that the domain $\Omega$ is more and more filled out. Then, points are said to be quasi-uniform if (12) is satisfied by all the sets in this sequence with the same constant $C_{qu}$ [34].

Furthermore, we need some technical considerations on the regularity of the covering $\{\Omega_{j}\}_{j=1}^{d}$ and thus we give the following definitions.

Definition 2.3**.**

A subdomain $\Omega_{j}\subseteq\mathbb{R}^{M}$ satisfies an interior cone condition if there exists an angle $\theta\in(0,\pi/2)$ and a radius $\gamma>0$ such that, for all $\boldsymbol{x}\in\Omega_{j}$, a unit vector $\boldsymbol{\xi}(\boldsymbol{x})$ exists such that the cone

[TABLE]

is contained in $\Omega_{j}$.

Definition 2.4**.**

Suppose that $\Omega\subseteq\mathbb{R}^{M}$ is bounded and ${\cal X}_{N}=\{\boldsymbol{x}_{i},i=1,\ldots,N\}\subseteq\Omega$ is given. An open and bounded covering $\{\Omega_{j}\}_{j=1}^{d}$ is called regular for $(\Omega,{\cal X}_{N})$, if the following properties are satisfied:

• i.

for each $\boldsymbol{x}\in\Omega$, the number of subdomains $\Omega_{j}$, with $\boldsymbol{x}\in\Omega_{j}$ is bounded by a global constant $C$,

• ii.

each subdomain $\Omega_{j}$ satisfies an interior cone condition,

• iii.

the local fill distances $h_{{\cal X}_{N_{j}},\Omega_{j}}$ are uniformly bounded by the global fill distance $h_{{\cal X}_{N},\Omega}$.

Letting $C_{\nu}^{k}(\mathbb{R}^{M})$ the space of all functions $f\in C^{k}$ whose derivatives of order $|\boldsymbol{\beta}|=k$ satisfy $D^{\boldsymbol{\beta}}f(\boldsymbol{x})={\cal O}(||\boldsymbol{x}||_{2}^{\nu})$ for $||\boldsymbol{x}||_{2}\longrightarrow 0$, we consider the following convergence result [17, 34]:

Theorem 2.1**.**

Let $\Omega\subseteq\mathbb{R}^{M}$ be open and bounded and assume that ${\cal X}_{N}=\{\boldsymbol{x}_{i},i=1,$ $\ldots,N\}\subseteq\Omega$. Let $\phi\in C_{\nu}^{k}(\mathbb{R}^{M})$ be a strictly conditionally positive definite function of order $L$. Let $\{\Omega_{j}\}_{j=1}^{d}$ be a regular covering for $(\Omega,{\cal X}_{N})$ and let $\{W_{j}\}_{j=1}^{d}$ be $k$-stable for $\{\Omega_{j}\}_{j=1}^{d}$. Then the error between $f\in{\cal N}_{\phi}(\Omega)$, where ${\cal N}_{\phi}$ is the native space of $\phi$, and its PU interpolant is bounded by

[TABLE]

for all $\boldsymbol{x}\in\Omega$ and all $|\beta|\leq k/2$.

Remark 2.2**.**

The first assumption in Definition 2.4 plays a crucial role also in the implementation of the PU method. In fact, such property leads to the requirement that the number of subdomains is proportional to the number of data [34].

2.2 Choosing suitable shape parameters and PU

subdomain sizes

Usually, the shape parameter $\varepsilon$ can greatly affect the accuracy of the resulting interpolant. Therefore, techniques allowing to select a predicted optimal shape parameter via error estimates have already been designed. Precisely, if the function is supposed to be known, the error can be exactly evaluated and thus the optimal shape parameter can be found without uncertainty. Otherwise, all the techniques based on error estimates give an approximated optimal value. Anyway, with abuse of notation, in what follows we will use the term optimal in the sense that such approximation of the optimal value is close to the one that can be found via trials and errors, for which the knowledge of the exact solution is supposed to be provided [17].

In the context of the PU method, aside from the value of the shape parameter, the size of the PU subdomains also plays a crucial role, especially when data with highly varying densities are considered. In literature, subdomains often consist of hyperspherical patches of the same radius $\delta$ [7, 17, 32]. Here, always considering hyperspherical patches, we propose a novel method that allows to suitably select both the radius $\delta_{j}$ and the shape parameter $\varepsilon_{j}$ for each PU subdomain $\Omega_{j}$, basing our considerations on an a priori error estimate.

We will focus on the so-called cross-validation algorithm, see [17, 21], properly modified for a bivariate optimization problem. The cross-validation scheme has been firstly introduced in [1, 22]. A variant of such method, known in literature as LOOCV, is detailed in [29]. Recent modifications of the cross-validation method can be found in [19], where LOOCV is interpreted in the context of PDEs, and in [33].

Such approaches are always used in order to find the optimal value of the shape parameter for a global interpolation problem. Here instead we are interested in selecting, for each PU subdomain, the optimal couple $(\delta_{j},\varepsilon_{j})$. Carefully choosing, for each hypersherical patch such couple, leads to an accurate computation of the PU interpolant. In fact, supposing to have a regular covering, if we compare the result reported in Theorem 2.1 with the global error estimate shown in [34], we can see that the PU method preserves the local approximation order for the global fit. Thus, the problem truly reduces in finding accurate local interpolants. In other words, if we improve the accuracy of the local fits, then also the one of the PU interpolant has benefits. This is even more evident from the following simple upper bound

[TABLE]

which shows that the PU approximation error is governed by the worst local error.

Let us consider an interpolation problem on $\Omega_{j}$ of the form (2) and, for a fixed $i\in\{1,\ldots,N_{j}\}$, let

[TABLE]

be the $j$-th interpolant obtained leaving out the $i$-th data on $\Omega_{j}$. Moreover let

[TABLE]

be the error at the $i$-th point. Then the quality of the local fit is determined by some norm of the vector of errors $\boldsymbol{e}_{j}=(e^{j}_{1},\ldots,e^{j}_{N_{j}})^{T},$ obtained by removing in turn one of the data points and comparing the resulting fit with the known value at the removed point. Following [17, 29], we can simplify the computation to a single formula by calculating

[TABLE]

where $c^{j}_{i}$ is the $i$-th coefficient of the RBF interpolant $R_{j}$ based on the full data set and $\left(A_{j}\right)_{ii}^{-1}$ is the $i$-th diagonal element of the inverse of the corresponding local interpolation matrix.

Precisely, in order to obtain an error estimate, we compute the following vector

[TABLE]

In order to select the optimal couple $(\delta_{j},\varepsilon_{j})$ for each PU subdomain, we compute (15) for several values of the radius $(\delta_{j_{1}},\ldots,\delta_{j_{P}})$ and of the shape parameter $(\varepsilon_{j_{1}},\ldots,\varepsilon_{j_{Q}})$.

In (15), the dependence of the errors from the cardinality of the PU subdomain $N_{j}$, i.e. from the PU radius, is evident. Moreover, in this work, to stress the dependence of (15) also from the shape parameter, for a fixed $p\in\{1,\ldots,P\}$ and a fixed $q\in\{1,\ldots,Q\}$, we will use the notation

[TABLE]

Thus, focusing on the maximum norm, we compute

[TABLE]

Note that (16) provides an error estimate for several values of the PU radius and of the shape parameter. Therefore the $j$-th local approximant is computed considering the couple $(\delta_{j},\varepsilon_{j})$ if

[TABLE]

In other words, the BLOOCV-PU interpolant assumes the form

[TABLE]

where, for each subdomain $\Omega_{j}$, $\tilde{R}_{j}$ is given by

[TABLE]

and $\tilde{N}_{j}$ indicates the number of points in $\Omega_{j}$ of radius $\delta_{j}$.

Observe that, consistently with Definition 2.4 and Remark 2.2, if the number of patches is proportional to $N$ and if the subdomains form a covering of $\Omega$, then such covering is also regular. This trivially follows from the fact that a hypersphere of radius $\delta_{j}$ always satisfies an interior cone condition with constants independent from the space dimension; precisely, $\gamma=\delta_{j}$ and $\theta=\pi/3$ [34]. Therefore, all the considerations made in the previous subsection also hold for the BLOOCV-PU interpolant.

This approach obviously leads to a benefit in terms of accuracy, especially when irregular data are considered. However, we have to point out that the computation of the inverse for each couple $(\delta_{j_{p}},\varepsilon_{j_{q}})$ is particularly costly for large $\delta_{j_{p}}$. Therefore, for each $\Omega_{j}$ we need to carefully choose the extreme values of the discrete searching range for the radius, i.e. the interval $[\delta_{j_{1}},\delta_{j_{P}}]$.

Precisely, for each PU subdomain, we have at first to fix the intervals $[\delta_{j_{1}},\delta_{j_{P}}]$ and $[\varepsilon_{j_{1}},\varepsilon_{j_{Q}}]$, used to find out $(\delta_{j},\varepsilon_{j})$. Many researchers already worked on the problem of finding suitable values for shape parameter in order to increase the accuracy and, at the same time, avoid problems of instability. Thus, one can easily guess how to select a good range for the shape parameter [6, 15, 17, 18, 20]. In other words, for what concerns the shape parameter the notation simplifies, since for each subdomain we can consider the same discrete values, namely $(\varepsilon_{1},\ldots,\varepsilon_{Q})$. On the opposite, fixing for all the subdomains the same discretization $(\delta_{1},\ldots,\delta_{P})$ can lead to inaccurate solutions. More specifically:

Problem 2.1**.**

Arbitrarily fixing, for all the subdomains, the same searching interval $[\delta_{1},\delta_{P}]$ can lead to the following issues:

the union of the PU subdomains might not form a covering of the domain;

• 2.

in regions characterized by a low density of points the interval $[\delta_{1},\delta_{P}]$ can be too small to avoid empty patches or subdomains containing very few points;

• 3.

in regions characterized by a high density of points, the interval $[\delta_{1},\delta_{P}]$ can be too large and, in this case, both complexity and ill-conditioning grow.

In the subsequent section we detail a feasible scheme useful to determine the interval $[\delta_{1},\delta_{P}]$, in which we can search for the $j$-th suitable radius, avoiding the above mentioned problems. To reach this aim, we first need an efficient partitioning structure, used to organize points among the different subdomains. Therefore, we propose a novel multidimensional procedure, built ad hoc for the PU method and independent from the problem geometry.

3 Feasible computation of the BLOOCV-PU

interpolant

As already pointed out in Subsection 2.2, the searching interval $[\delta_{1},\delta_{P}]$ must be properly selected. Essentially, in order to obtain reliable error estimates, we want to make sure of having a reasonable number of points on each patch. Such consideration suggests the use of a $K$-nearest neighbor procedure. As example in [13], suitable supports of CSRBFs have been selected detecting, via a triangulation, the $K$-nearest neighbor set [2]. This turns out to be expensive and moreover fixing an arbitrary $K$ does not guarantee a good approximation. Thus, we will use a similar procedure to [13] only to determine the initial reasonable searching range $[\delta_{1},\delta_{P}]$ for the radius of the $j$-th patch and then such interval will be used in the computation of (16).

The complexity needed to construct the BLOOCV-PU interpolant will be taken into account. Specifically, we will not perform a $K$-nearest neighbor procedure, but we will use the new multidimensional I-PS. It leads to a considerable saving in terms of computational time with respect to [7, 8] and unlike them it can be applied in any space dimension $M$ (and not only for $M=2,3$).

We will treat the problem in the most general setting. Thus let us consider a set of data ${\cal X}_{N}=\{\boldsymbol{x}_{i}\in\Omega,i=1,\ldots,N\}$, where $\Omega\subseteq\mathbb{R}^{M}$ is a simply connected region. In order to perform the BLOOCV-PU method, we need to define an hyperrectangle ${\cal R}_{M}$ containing the scattered data

[TABLE]

and the bounding box containing the nodes, i.e. the box of edge

[TABLE]

Then, consistently with Remark 2.2, we define the PU centres as a grid of $d^{{\cal R}_{M}}=\left(d_{PU}^{{\cal R}_{M}}\right)^{M}$ points on ${\cal R}_{M}$, where

[TABLE]

and $V_{\Omega}$ is the hypervolume of the simply connected region $\Omega$. Then, to make sure that patches form a covering of the domain, we can set the radii of the hyperspheres $\delta_{j}$ such that

[TABLE]

With (21) we solve the issue $1.$ outlined in Problem 2.1.

The initial number of subdomains $d^{{\cal R}_{M}}$ is later reduced by taking only those $d$ centres lying in $\Omega$. This step provides the set of PU centres ${\cal C}_{d}=\{\boldsymbol{\bar{x}}_{i},i=1,\ldots,d\}$, used to construct the PU interpolant. Note that, because of (21), the set ${\cal C}_{d}$ forms a covering for $\Omega$. Furthermore, in the same way we also define a set ${\cal E}_{s}=\{\boldsymbol{\tilde{x}}_{i},i=1,\ldots,s\}$ on $\Omega$, which is used to evaluate the unknown function via BLOOCV-PU interpolation. Then, in order to organize points into the different patches and consequently choose a suitable searching interval for the radius (see Subsection 3.2), we consider the partitioning structure described in Subsection 3.1.

3.1 Multidimensional integer based partitioning structure

To make simpler the presentation, we first consider hyperspherical patches all having the same radius [8]

[TABLE]

Then, in order to solve the local interpolation problems, we need to develop a procedure enabling us to store the points among the different PU subdomains. Such scheme must be independent from the problem geometry and work in any dimension, as kd-trees [2, 12, 17, 37]. These are effective and widely used numerical tools, but they are not specifically constructed for the PU method.

Thus, starting from the bivariate and trivariate procedures, that in the following we will call the Sorting-based Partitioning Structures (S-PSs) [7, 8], our aim is to build a multidimensional procedure which allows to consider variable radii $\delta_{j}$. As in the S-PS, we store the points into $q^{M}$ blocks, where

[TABLE]

More precisely, we number blocks from $1$ to $q^{M}$, starting from the subspace of dimension $M-1$, obtained projecting along the first coordinate and thus parallel to the remaining ones. In order to fix the idea, in a 2D context they are numbered from bottom to top, left to right.

Remark 3.1**.**

In bivariate interpolation blocks are generated by the intersection of two orthogonal strips. In multivariate problems blocks are generated by the intersection of $M$ hyperrectangles. In what follows with abuse of notation we will continue to call such hyperrectangles with the term strips.

Then, in order to store the points among the different patches the following computational issue, known as containing query, needs to be solved

• •

given a PU centre $\boldsymbol{\bar{x}}_{j}$, find the $k$-th block containing the centre.

Such problem can be easily solved taking into account that, given a PU centre $\boldsymbol{\bar{x}}_{j}$, if $k_{m}$ is the index of the strip parallel to the subspace of dimension $M-1$ generated by $x_{r}$, $r=1,\ldots,M$ and $r\neq m$, containing the $m$-th coordinate of $\boldsymbol{\bar{x}}_{j}$, then the index of the $k$-th block containing the subdomain centre is

[TABLE]

To find the indices $k_{m}$, $m=1,\ldots,M$, in (24), we use an integer-based procedure consisting in rounding off to an integer value. Specifically, for each PU centre $\boldsymbol{\bar{x}}_{j}=(\bar{x}_{j1},\ldots,\bar{x}_{jM})$, we have that

[TABLE]

Then, exactly the same procedure is adopted in order to store into the different blocks both scattered data and evaluation points, i.e. the I-PS assigns:

• i)

to each scattered point $\boldsymbol{x}_{i}$, $i=1,\ldots,N$, the index of the block in which it lies,

• ii)

to each evaluation point $\boldsymbol{\tilde{x}}_{i}$, $i=1,\ldots,s$, the index of the block in which it lies.

Moreover, always supposing to have a fixed radius and assuming that the $j$-th centre belongs to the $k$-th block, from (23) the fact that we search for the points lying in the $j$-th patch among those lying in the $k$-th block and in its $3^{M}-1$ neighboring blocks easily follows.

On the opposite, here the radius is supposed to be variable for each patch and thus, if the radius $\delta_{j}$ is such that

[TABLE]

given the centre $\boldsymbol{\bar{x}}_{j}$ we search for the neighboring points in the $k$-th block and in its $(3+2^{n})^{M}-1$ neighboring blocks.

In [7, 8] nodes and evaluation points are organized in blocks by using recursive calls to a sorting routine, while here this step is replaced by (24) and (25). Such approach enables us to improve three aspects of the partitioning structure presented in the last mentioned papers. Precisely, the I-PS:

works in any dimension, while the S-PS only works for $M=2,3$,

• 2)

allows to work with variable radii, while the S-PS strictly depends on a fixed size of the subdomains,

• 3)

reduces the complexity of the sorting-based storing procedure (see Section 4).

3.2 Selection of a searching interval for the PU radius

When we deal with quasi-uniform or grid data, the number of points in each subdomain of radius $\delta$ is about constant. On the opposite, in case of irregular nodes, the number of points lying in the different patches is far from being constant or, even worst, we can have empty subdomains. This consideration turns out to be useful to determine the lower bound of the interval $[\delta_{1},\delta_{P}]$ for the $j$-th subdomain. In fact, given the number $N$ of scattered data in $\Omega$ and its hypervolume $V_{\Omega}$, from a simple proportion we have that a suitable number $K$ of points belonging to $\Omega_{j}$ of radius (21) is

[TABLE]

where $B(\delta)$ is the hypervolume of the hypersphere of radius $\delta$, defined as in (22). The value found in (26) represents the number of points we expect on average on each patch supposing to have a uniform node distribution.

Therefore, given the $j$-th subdomain of radius $\delta_{j_{1}}=\delta$, we compute its cardinality via the I-PS. Then, if such cardinality is less than the one given by (26), $\delta_{j_{1}}$ is updated as follows

[TABLE]

where $0<t<1$. The procedure continues in this way until (26) is satisfied, i.e. $\delta_{j_{1}}$ is determined, with recursive calls to the I-SP, so that

[TABLE]

Acting in this way, we solve the computational issue $2.$ outlined in Problem 2.1, i.e. there are few enough points for each patch.

Then, in order to avoid also the third issue of Problem 2.1, the simplest strategy, which takes into account the density of points and turns out to be effective, is to choose $P$ discrete values in an interval of the form

[TABLE]

Roughly speaking, since the upper bound of the searching interval is proportional to the lower bound and since this lower bound is large only if the density of points is low, we effectively avoid problems arising from high density of points, i.e. systems are not too large and the ill-conditioning is kept under control.

We end this section with the illustrative Figure 1, devoted to show how the classical PU structure is modified by means of the BLOOCV-PU algorithm. In the left frame we plot the classical structure by choosing $\delta_{j}=\delta$, $j=1,\ldots,d$, while in the right frame we show the result of the BLOOCV-PU method.

4 Complexity analysis

In this section we analyze the computational cost of the proposed method. It will be pointed out that the main cost is due to the computation of the error estimates, while the I-PS turns out to be really cheap.

4.1 The BLOOCV local implementation

For each PU subdomain, several error estimates are calculated via (14). Such computation needs ${\cal O}(N_{j}^{3})$ operations. Thus, simplifying the calculation via (14) is the key step which enables us to maintain a reasonable complexity cost. Indeed, evaluating the error via (13) is computationally expensive. In particular, the matrix inverse, which requires ${\cal O}(N_{j}^{3})$ operations, must be computed for each node. This step needs a total computational cost of ${\cal O}(N_{j}^{4})$ operations, $j=1,\ldots,d$, but using (14), the complexity cost significantly decreases.

However, the complexity of the proposed algorithm is quite high. The error estimate (14) needs to be computed for each subdomain $\Omega_{j}$ and for each $\delta_{j_{i}}$, $i=1,\ldots,P$, and $\varepsilon_{k}$, $k=1,\ldots,Q$.

Remark 4.1**.**

The computation of the error estimate for the shape parameter can be slightly speeded up by using the Matlab routine CostEpsilon.m, proposed in [17]. However, the same approach cannot be performed for an optimal choice of the radius $\delta_{j}$. Therefore, in our free software package, for easiness of the Matlab user, we carried out the standard implementation as in (16).

4.2 The integer-based partitioning structure

The I-PS, after organizing the scattered data into the different blocks, given a subdomain $\Omega_{j}$ searches for all the points lying in $\Omega_{j}$ in a reduced number of blocks. Specifically, in order to store the scattered data among the different blocks, it makes use of an integer-based procedure that assigns to each node $N_{i}$, $i=1,\ldots,N$, the corresponding block. This step requires ${\cal O}(N)$ time. Then, we apply the optimized searching routine already used in [7, 8]. Such procedure is performed in a constant time (refer to [8] for further details).

Observe that the I-PS turns out to be more efficient than the S-PS; in fact the latter, to store the points among the different blocks, needs ${\cal O}(N\log N)$ operations. Table 1 and Figure 2 support our findings. Specifically, we consider in a 2D framework different sets of Halton data. Tests have been carried out on a Intel(R) Core(TM) i7 CPU 4712MQ 2.13 GHz processor.

4.3 Computation of the PU interpolant

The computation of the local interpolation problems consists in solving $d$ linear systems of size $N_{j}\times N_{j}$, with $N_{j}\ll N$. This step involves a computational cost of ${\cal O}(N^{3}_{j})$ time, $j=1,\ldots,d$. Since the number $d$ of subdomains is bounded by ${\cal O}(N)$, this leads to ${\cal O}(N)$ operations for solving all of them. Finally, in order to get the value of the global fit (1), we have to add up a constant number of local RBF interpolant. This requires ${\cal O}(1)$ time.

5 Numerical experiments

This section is devoted to show, by means of extensive numerical simulations the flexibility and the accuracy of the proposed BLOOCV-PU method. It is applied fixing the initial intervals for the radii as in (29), with $h=2$ and $P=6$. These values are chosen so that we ensure to have enough points on each subdomain and a sufficient number of radii to test the accuracy of the interpolants. Moreover, we will use as PU weights the Wendland $C^{2}$ function, see (8).

Tests are carried out considering the so-called product and valley functions [4, 28], respectively defined as:

[TABLE]

To point out the accuracy of the BLOOCV-PU interpolant, we will refer to the Maximum Absolute Error (MAE) and the Root Mean Square Error (RMSE), whose formulas are:

[TABLE]

where $\tilde{\boldsymbol{x}}_{i}$, $i=1,\ldots,s$, forms a grid of $40\times 40$ points in which the interpolant is sought.

Concerning the data sets used in our numerical experiments we take, as quasi-uniform points, the Halton data [17], see Figure 3 (left). Then, in order to test the method with particularly hard nodes, we consider points coming from a Schwarz-Christoffel transformation, refer to [14, 16, 23] for further details. More precisely, we focus on a special case of conformal map from the unit disk onto a polygon. Thus, at first we define nodes in the unit disk and then we map them into a chosen polygon. As example, in Figure 3 (right) we show the result of conformally mapping onto a simply connected region points on concentric circles. We consider such points because they are far from being quasi-uniform, same time they are constructed with a specific rule and thus tests are repeatable. Moreover, they simulate practical situations, as it will be evident in Section 6 when we will deal with data coming from real life.

The BLOOCV-PU approach will be compared with the classical PU method, i.e. the PU scheme is applied with a fixed radius and a fixed shape parameter [7, 17, 34]. Thus, such classical approach requires to fix these two parameters. Concerning the radius, the classical PU method is generally applied considering a fixed size of the patches as in (22), while for the shape parameter, in literature, the choice is almost arbitrary. In fact, even if techniques allowing to obtain stable approximations when $\varepsilon\rightarrow 0$ have already been developed [6, 20], there is not an a priori good value for the latter. We will later discuss this concept providing several tests.

In our numerical experiments, we expect the following behavior classes depending on the distribution of the data set:

• i.

with quasi-uniform points: BLOOCV-PU and classical PU both give accurate approximations,

• ii.

with non-conformal points: the classical PU method fails, while the BLOOCV-PU maintains a good accuracy.

5.1 Tests with quasi-uniform data

In this subsection we consider Halton points and, because of their regularity, we take a smooth RBF as local approximant, specifically the IMQ, see (6). Therefore, choosing $30$ values for the shape parameter in the interval $[0.1,10]$ is reasonable. Results, obtained by means of the BLOOCV-PU interpolant with the functions $f_{1}$ and $f_{2}$, are shown in Tables 2 and 3, respectively. We also report the errors of the classical PU method, obtained by fixing the radius as in (22) and the shape parameter $\varepsilon=0.6$. For a graphical representation of the distribution of the absolute error refer to Figure 4. Such figure shows that the error for both test functions is larger close to the boundary. Moreover, for $f_{2}$ it also increases in correspondence of the oscillations of the test function.

From Tables 2 and 3, we can easily note that the results obtained with the BLOOCV-PU interpolant are more accurate than the ones carried out with the standard PU. Anyway, in order to get such accurate approximation, we have to pay in terms of efficiency. For instance, with $289$ and $1089$ points the BLOOCV-PU can be computed in $5.48$ s and $20.08$ s, respectively, while the classical PU interpolation only requires $0.21$ s and $0.60$ s.

One may argue that the results of the classical PU method might truly improve by varying the shape parameter. This is trivially true, but at the same time automatically choosing a safe shape parameter, which gives reliable approximations, is one of the main advantages of the BLOOCV-PU method. Anyway, in order to clarify this concept, we report in Table 4 the results obtained by choosing in the classical PU algorithm the optimal shape parameter, but not the optimal radius, for each subdomain. Note that this is only a particular case of the BLOOCV-PU algorithm and thus no longer investigated.

5.2 Tests with non-conformal data points

Let us now turn into the more complex case of irregular data, which gives the results shown in Tables 6 and 7. We use as data distribution the one shown in Figure 3 (right). In order to have a better understanding of their distribution we report in Table 5 the two indicators of data regularity, i.e. the separation distance (10) and the fill distance (11), of the non-conformal points and we compare them with Halton data. The last column gives an idea of the quasi-uniformity constant, see Remark 2.1. As evident, points coming from a Schwarz-Christoffel transformation are really far from being quasi-uniform. In this case, the interpolation process results particularly hard. In fact, we remark that the ill-conditioning primarily grows due to the decrease of the separation distance.

Since points are not well-distributed and ill-conditioning is expected, we choose as local approximant a CSRBF, specifically the Wendland $C^{6}$, see (9).

The classical PU approach is applied by fixing the shape parameter $\varepsilon=0.5$ and the radius as in (22).

In this case the classical PU method with a fixed size of the subdomains gives inaccurate approximations. We verify with numerical experiments that this does not depend on the shape parameter; in fact, neither the use the optimal $\varepsilon$ changes the order of the approximation errors. The BLOOCV-PU reveals its robustness in both cases. It follows that such method with the related flexibility of the patches turns out to be meaningful when irregular data sets are considered.

6 Applications to real world data sets

This section is devoted to test the BLOOCV-PU scheme with two real world data sets. The first one we consider is the so-called glacier data set. It consists of $8345$ points representing digitized height contours of a glacier [9, 31, 36]. The difference between the highest and the lowest point is $800$ m.

As second example, we consider the so-called black forest data set [10, 11]. It consists of $15885$ points representing a terrain in the neighborhood of Freiburg, Germany. In this case, the difference between the maximal and minimal heights is $1214$ m.

A 2D view of the data sets is plotted in Figure 5 (top left to right, respectively). Because of the high variability of the points, in both cases we use as local approximant a RBF characterized by a finite regularity, such as the Matérn $C^{2}$, see (7).

Figure 5 (bottom) shows the reconstruction of the surfaces defined by the two data sets. It has been obtained evaluating the BLOOCV-PU interpolant on a grid of $80\times 80$ points.

To test the accuracy we use, as validation set, $90$ and $170$ points of the glacier and black forest data sets, respectively. Such data are plotted in red in Figure 5. The errors obtained in these cases are shown in Table 8.

7 Concluding remarks

In this paper we provided a robust tool enabling us to safely select, for each PU subdomain, both its size and the shape parameter. Numerical evidence and applications with real world measurements show that the proposed method accurately fits data with highly varying densities. Moreover, the BLOOCV-PU implementation is carried out with a new multidimensional searching procedure which has been proved to be extremely fast.

Work in progress consists in varying the shape of the PU subdomains, which here are supposed to be hyperspherical patches. This is not trivial since several requirements for the covering might be not easily satisfied.