Positive constrained approximation via RBF-based partition of unity method

TL;DR

This paper introduces a local, positive-constrained RBF-based partition of unity method that improves accuracy in constructing positive interpolants, with applications demonstrated in population dynamics.

Contribution

It presents a novel local approach using PU techniques for positive RBF interpolation, avoiding restrictions to CSRBFs and optimizing constraints via error estimates.

Findings

Enhanced accuracy through local constraints

Effective in population dynamics applications

Flexible approach not limited to CSRBFs

Abstract

In this paper, we discuss the problem of constructing Radial Basis In this paper, we discuss the problem of constructing Radial Basis Function (RBF)-based Partition of Unity (PU) interpolants that are positive if data values are positive. More specifically, we compute positive local approximants by adding up several constraints to the interpolation conditions. This approach, considering a global approximation problem and Compactly Supported RBFs (CSRBFs), has been previously proposed. Here, the use of the PU technique enables us to intervene only locally and as a consequence to reach a better accuracy. This is also due to the fact that we select the optimal number of positive constraints by means of an a priori error estimate and we do not restrict to the use of CSRBFs. Numerical experiments and applications to population dynamics are provided to illustrate the effectiveness of the method in applied sciences.

Tables

Table 1. Table 1: The PC-PU Algorithm . Routine performing the PC-PU method.

INPUTS: $N$, number of data; ${𝒳}_N=\{(x_i,y_i),i=1,\dots ,N\}$, set of data points;                    ${ℱ}_N=\{f_i,i=1,\dots ,N\}$, set of data values; $d$, number of                    subdomains; $s$, number of evaluation points. OUTPUTS: ${𝒜}_s=\{\widehat{ℐ}({\tilde{x}}_i,{\tilde{y}}_i),i=1,\dots ,s\}$, set of approximated values. Step 1: A grid of evaluation points ${ℰ}_s=\{({\tilde{x}}_i,{\tilde{y}}_i),i=1,\dots ,s\}\subseteq \Omega$ is generated. Step 2: A grid of subdomain points ${𝒞}_d=\{({\overline{x}}_j,{\overline{y}}_j),j=1,\dots ,d\}\subseteq \Omega$ is                    constructed. Step 3: For each PU centre $({\overline{x}}_j,{\overline{y}}_j)$, $j=1,\dots ,d$, a subdomain, whose                    radius is given by (15), is constructed. Step 4: For each patch ${\Omega }_j$, $j=1,\dots ,d$,             Step 5: Find all data points ${𝒳}_{N_j}$ belonging to ${\Omega }_j$.             Step 6: Find all evaluation points ${ℰ}_{s_j}$ belonging to ${\Omega }_j$.             Step 7: Solve the unconstrained interpolation problem and a local                              interpolant $R_j$ is formed as in (2).             Step 8: If the local fit is positive ${\hat{N}}_j=0$ and ${\hat{R}}_j=R_j$, else                       Step 9: For ${\hat{N}}_j=1,\dots ,N_j$, compute the constraints as in (14)                                       and calculate (10).                              Consider ${\hat{N}}_j$ constraints, if ${\hat{N}}_j$ satisfies (13) and (6).             Step 10: Solve the constrained approximation problem and form a                              positive local approximant ${\hat{R}}_j$, see (5). Step 11: The interpolant (8) is formed by the weighted sum of the local fits.

Table 2. Table 2: MAEs and RMSEs computed with the Wendland’s C 2 superscript 𝐶 2 C^{2} kernel for f 1 subscript 𝑓 1 f_{1} .

$N$	method	MAE	RMSE
$300$	PU	$1.50E-01$	$1.52E-02$
	PC-PU	$1.50E-01$	$2.03E-02$
$1000$	PU	$7.36E-02$	$3.08E-03$
	PC-PU	$7.96E-02$	$6.44E-03$
$3500$	PU	$6.34E-02$	$1.47E-03$
	PC-PU	$8.40E-02$	$2.86E-03$
$8000$	PU	$2.43E-02$	$4.17E-04$
	PC-PU	$5.99E-02$	$1.03E-03$

Table 3. Table 3: MAEs and RMSEs computed with the IMQ C ∞ superscript 𝐶 C^{\infty} kernel for f 1 subscript 𝑓 1 f_{1} .

$N$	method	MAE	RMSE
$300$	PU	$1.39E-01$	$1.04E-02$
	PC-PU	$1.39E-01$	$1.44E-02$
$1000$	PU	$7.02E-02$	$2.88E-03$
	PC-PU	$7.02E-02$	$3.49E-03$
$3500$	PU	$5.89E-02$	$1.50E-03$
	PC-PU	$5.89E-02$	$1.66E-03$
$8000$	PU	$2.33E-02$	$3.46E-04$
	PC-PU	$2.33E-02$	$3.68E-04$

Table 4. Table 4: MAEs and RMSEs of the CSRBF-based global method and of the PC-PU approximant computed with the Wendland’s C 2 superscript 𝐶 2 C^{2} kernel for f 2 subscript 𝑓 2 f_{2} .

$N$	method	MAE	RMSE
$300$	CSRBF global	$3.36E-01$	$2.77E-01$
	PC-PU	$2.76E-01$	$1.91E-02$
$1000$	CSRBF global	$4.31E-01$	$2.59E-02$
	PC-PU	$8.84E-02$	$5.95E-03$
$3500$	CSRBF global	$1.69E-01$	$6.67E-02$
	PC-PU	$8.48E-02$	$2.61E-03$

Table 5. Table 5: MAEs and RMSEs of the Shepard’s method and of the PC-PU approximant computed with the IMQ C ∞ superscript 𝐶 C^{\infty} kernel for f 2 subscript 𝑓 2 f_{2} .

$N$	method	MAE	RMSE
$300$	Shepard	$8.38E-01$	$3.21E-01$
	PC-PU	$1.32E-01$	$1.48E-02$
$1000$	Shepard	$5.71E-01$	$7.75E-02$
	PC-PU	$8.62E-02$	$4.31E-03$
$3500$	Shepard	$1.87E-01$	$1.42E-02$
	PC-PU	$2.89E-02$	$9.73E-04$

Equations

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

$$\mathcal{I}(x_i,y_i)=\sum _{j=1}^d{R}_j(x_i,y_i)W_j(x_i,y_i)=\sum _{j\in I(x_i,y_i)}{f}_i{W}_j(x_i,y_i)=f_i,$$

$$R_j(x,y)=\sum _{k=1}^{N_j}{c}_k{\varphi }_{\epsilon }^{(k)}(x,y),(x,y)\in {\Omega }_j,$$

$${\varphi }_{\epsilon }^{(k)}(x,y)=\frac{1}{\sqrt{1+{\epsilon }^2{r}_k^2(x,y)}},$$

$${\varphi }_{\epsilon }^{(k)}(x,y)=(1-\epsilon r_k(x,y){)}_{+}^4(4\epsilon r_k(x,y)+1),$$

$$R_j(x_i,y_i)=f_i,i=1,\dots ,N_j,$$

$$A^{(j)}\mathbf{c}=\mathbf{f},$$

$$A_{ik}^{(j)}={\varphi }_{\epsilon }^{(k)}(x_i,y_i),i,k=1,\dots ,N_j.$$

$$h_{{\mathcal{X}}_N,\Omega }=\underset{(x,y)\in \Omega }{\sup }\left(\underset{(x_i,y_i)\in {\mathcal{X}}_N}{\min }{r}_i(x,y)\right).$$

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

$$({\hat{x}}_{N_j+1},{\hat{y}}_{N_j+1}),\dots ,({\hat{x}}_{N_j+{\hat{N}}_j},{\hat{y}}_{N_j+{\hat{N}}_j}),$$

$${\hat{R}}_j(x,y)=\sum _{k=1}^{N_j}{c}_k{\varphi }_{\epsilon }^{(k)}(x,y)+\sum _{\hat{k}=N_j+1}^{N_j+{\hat{N}}_j}{c}_{\hat{k}}{\hat{\varphi }}_{{\epsilon }_{\hat{k}}}^{(\hat{k})}(x,y),(x,y)\in {\Omega }_j,$$

$${\hat{R}}_j(x_i,y_i)=f_i,i=1,\dots ,N_j,c_i\ge 0,i=1,\dots ,N_j+{\hat{N}}_j,$$

$${\mathbf{a}}_k^T=-\{{\varphi }_{\epsilon }^{(k)}(x_i,y_i){\}}_{i=1}^{N_j},k=1,\dots ,N_j,$$

$${\mathbf{a}}_{\hat{k}}^T=-\{{\hat{\varphi }}_{{\epsilon }_{\hat{k}}}^{(\hat{k})}(x_i,y_i){\}}_{i=1}^{N_j}=(0,\dots ,0,-b_{\hat{k}},0,\dots ,0),\hat{k}=N_j+1,\dots ,2N_j,$$

$${\mathbf{a}}_{2N_j+1}^T=(f_1,\dots ,f_{N_j}),$$

$$\min {\left(\sum _{i=N_j+1}^{N_j+{\hat{N}}_j}{c}_i^2\right)}^{1/2},$$

$$\hat{\mathcal{I}}(x,y)=\sum _{j=1}^d\left(\sum _{k=1}^{N_j}{c}_k{\varphi }_{\epsilon }^{(k)}(x,y)+\sum _{\hat{k}=N_j+1}^{N_j+{\hat{N}}_j}{c}_{\hat{k}}{\hat{\varphi }}_{{\epsilon }_{\hat{k}}}^{(\hat{k})}(x,y)\right)W_j(x,y).$$

$$R_j^i(x,y)=\sum _{k=1,k\ne i}^{N_j}{c}_k{\varphi }_{\epsilon }^{(k)}(x,y),(x,y)\in {\Omega }_j,$$

$$e_i=f_i-R_j^i(x_i,y_i),$$

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

$$\left({\hat{e}}_1,\dots ,{\hat{e}}_{N_j+{\hat{N}}_j}\right)=\left(\frac{c_1}{{\left({\hat{A}}_{11}^{(j)}\right)}^{-1}},\dots ,\frac{c_{N_j+{\hat{N}}_j}}{{\left({\hat{A}}_{N_j+{\hat{N}}_j{N}_j+{\hat{N}}_j}^{(j)}\right)}^{-1}}\right),$$

$${\hat{A}}^{(j)}=\left(\begin{array}{ccc}{\varphi }_{\epsilon }^{(1)}(x_1,y_1)& \cdots & {\hat{\varphi }}_{{\epsilon }_{N_j+{\hat{N}}_j}}^{(N_j+{\hat{N}}_j)}(x_1,y_1)\\ ⋮& \ddots & ⋮\\ {\hat{\varphi }}_{{\epsilon }_{N_j+{\hat{N}}_j}}^{(N_j+{\hat{N}}_j)}(x_1,y_1)& \cdots & {\hat{\varphi }}_{{\epsilon }_{N_j+{\hat{N}}_j}}^{(N_j+{\hat{N}}_j)}({\hat{x}}_{N_j+{\hat{N}}_j},{\hat{y}}_{N_j+{\hat{N}}_j})\end{array}\right).$$

$${\hat{\mathbf{e}}}_j({\hat{N}}_j)=\left({\hat{e}}_1,\dots ,{\hat{e}}_{N_j+{\hat{N}}_j}\right).$$

$$||{\hat{\mathbf{e}}}_j({\hat{N}}_j)|{|}_{\infty }=\underset{k=1,\dots ,N_j}{\min }||{\hat{\mathbf{e}}}_j(k)|{|}_{\infty },$$

$$({\hat{x}}_k,{\hat{y}}_k)=({\bar{x}}_j+u_k\cos {\eta }_k,{\bar{y}}_j+u_k\sin {\eta }_k),$$

$$u_k=\delta \frac{\sqrt{k-1/2}}{\sqrt{{\hat{N}}_j-1/2}}\text{and}{\eta }_k=\frac{4k\pi }{1+\sqrt{5}}.$$

$$\delta =\sqrt{\frac{1}{d}}.$$

$$MAE=\underset{1\le i\le s}{\max }|f({\tilde{x}}_i,{\tilde{y}}_i)-\hat{\mathcal{I}}({\tilde{x}}_i,{\tilde{y}}_i)|,$$

$$RMSE=\sqrt{\frac{1}{s}\sum _{i=1}^s|f({\tilde{x}}_i,{\tilde{y}}_i)-\hat{\mathcal{I}}({\tilde{x}}_i,{\tilde{y}}_i){|}^2}.$$

$$f_1(x,y)=(x-0.5{)}^2+(y-0.4{)}^2,$$

Full text

Positive constrained approximation via RBF-based partition of unity method

Alessandra De Rossi and Emma Perracchione

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

Abstract. In this paper, we discuss the problem of constructing Radial Basis Function (RBF)-based Partition of Unity (PU) interpolants that are positive if data values are positive. More specifically, we compute positive local approximants by adding up several constraints to the interpolation conditions. This approach, considering a global approximation problem and Compactly Supported RBFs (CSRBFs), has been previously proposed in [39]. Here, the use of the PU technique enables us to intervene only locally and as a consequence to reach a better accuracy. This is also due to the fact that we select the optimal number of positive constraints by means of an a priori error estimate and we do not restrict to the use of CSRBFs. Numerical experiments and applications to population dynamics are provided to illustrate the effectiveness of the method in applied sciences.

1 Introduction

Given a set of samples, the scattered data interpolation problem consists in finding an approximating function that matches the given measurements at their corresponding locations. Furthermore, dealing with applications, we often have additional properties, such as the non-negativity of the measurements, which we wish to be preserved during the interpolation process. Note that, since such property is known as positivity-preserving property, to keep a common notation with existing literature, we use the term positive (instead of non-negative) function values or interpolants.

To preserve such property, mostly considering rational spline functions with $C^{1}$ or $C^{2}$ continuity, the recent research studies techniques which force the approximants to be positive. As example, the conditions under which the positivity of a cubic piece may be lost are investigated in [30]. Moreover, in order to preserve the positivity, the use of bicubic splines, coupled with a technique based on adding extra knots, has been investigated in [2, 8]. A positive fit is instead obtained by means of rational cubic splines in [21, 23]. The same authors also developed a positive surface construction scheme for positive scattered data arranged over triangular grids [22].

Note that all the above mentioned methods depend on a mesh. However, the positivity-preserving problem is also well-known in the field of meshfree or meshless methods. They include Shepard-type approximants [27, 29] and RBF interpolants [7, 37]. While the positivity-preserving problem has been widely investigated for the Modified Quadratic Shepard’s (MQS) approximant [1, 6], it remains a challenging computational issue for RBF interpolation. Indeed, even if such meshfree approach has been extensively studied in the recent years, especially focusing on the stability of the interpolant [13, 16], not a lot of effort has been addressed to construct positivity-preserving approximants. Such problem has been studied only in particular and well-known cases; for instance, in [33] it is analyzed for the thin plate spline. Other methods, which follow from the quasi-interpolation formula given in [40], have been proposed and effectively performed to construct positive approximants of positive samples [36].

Focusing on RBFs, our scope consists in preserving the positivity of the PU interpolant for a wider family of kernels. In [39] a global positive RBF approximant is constructed by adding up to the interpolation conditions several positive constraints and considering CSRBFs. Even if the optimal number of constraints is not investigated, the results are promising and show that this technique has a better accuracy than the Constrained MQS (CMQS) approximant. However, since a global interpolant is used, adding up other constraints to preserve the positivity implies that the shape of the curve/surface is consequently globally modified. As pointed out in [39], this might lead to a considerable decrease of the quality of the approximating function in comparison with the unconstrained CSRBF interpolation.

Thus here, in order to avoid such drawback, focusing on 2D data sets, the PU method is performed by imposing positive constraints on the local RBF interpolants. Such approach enables us to consider constrained interpolation problems only in those PU subdomains which do not preserve the required property. This leads to an accurate method compared with existing techniques [1, 6, 39]. Moreover, in contrast with [39], we can consider truly large data sets.

Specifically, in order to construct the Positive Constrained PU (PC-PU) approximant, following [39], we locally impose several positive constraints and we reduce to solve an optimization problem. The number of constraints is properly selected by means of an a priori error estimate. This is a fundamental step to maintain a good accuracy of the fit. Moreover, differently from [39], a wider family of RBFs (not only compactly supported) is considered. The main disadvantage of using CSRBFs is that they introduce large errors in the area in which the interpolant is negative. This is due to the fact that the shape of the fit is modified only within the support of the CSRBF and thus neighbouring points are not taken into account. Numerical evidence shows that the use of infinitely smooth globally defined RBFs leads to an improvement in this direction. Moreover, in order to stress the effectiveness of the proposed technique in applied sciences, we investigate several applications to biomathematics. Finally, comparisons with the unconstrained PU interpolant, with the Shepard’s method and with the one outlined in [39] are carried out.

The guidelines of the paper are as follows. In Section 2 we investigate the positivity of the PU approximant by considering extra positive constraints. The computational aspects of such method are analyzed in Section 3. Numerical experiments and applications to population dynamics are shown in Section 4 and 5, respectively. Finally, Section 6 deals with conclusions and future work.

2 Positive approximation of positive data values

In addition to the assigned interpolation conditions, our goal consists in considering several positive constraints. This allows to preserve the positivity of the measurements. Therefore, our approach turns out to be meaningful especially in applications, indeed in order to avoid the violation of biological or physical measurements, a positive fit is often necessary. Thus, before discussing the proposed method and without loosing generality, we consider an example of univariate interpolation to illustrate the scope of such numerical tool and motivate the reader.

Example 2.1** (Motivations and targets).**

One of the most common tumor, affecting mainly men over sixty years old, is the prostate cancer. Luckily, this disease has a very slow growth and a reliable biomarker after prostatectomy, the so-called Prostate Specific Antigen (PSA), suitable for the early diagnosis. In particular, if its value is larger than $0.2$ ng/mL, a relapse (a local or distal metastasis) occurs.

Clinical data of prostatectomized patients are available in [18]. They are used to investigate the evolution of the relapse via mathematical models [31]. Therefore, in order to validate such models, a data fitting results essential [25]. For such scope, one can reconstruct the PSA curve with the standard PU interpolant. Anyway, the positivity of the PSA values is not always preserved through the interpolation process, violating the biological constraint. In order to avoid this problem, we propose a novel technique, namely the PC-PU method, which allows to preserve the positivity property, see Figure 1.

2.1 The partition of unity method

Let us consider a bounded set $\Omega\subseteq\mathbb{R}^{2}$. Given a set of $N$ distinct data points (or data sites or nodes) ${\cal X}_{N}=\{(x_{i},y_{i}),i=1,\ldots,N\}\subseteq\Omega$ and a set of data values (or measurements or function values) ${\cal F}_{N}=\{f_{i}=f(x_{i},y_{i}),i=1,\ldots,N\}$, the PU method decomposes the domain $\Omega$ into $d$ subdomains $\Omega_{j}$, $j=1,\ldots,d$, such that $\Omega\subseteq\bigcup_{j=1}^{d}\Omega_{j}$ [4, 15, 24, 38]. In literature, the subdomains $\Omega_{j}$ are supposed to be circular patches of the same radius $\delta$ covering the domain $\Omega$.

According to [38], associated with these subdomains, we choose a $p$-stable partition of unity, i.e a family of nonnegative functions $\{W_{j}\}_{j=1}^{d}$, with $W_{j}\in C^{p}(\mathbb{R}^{2})$, such that

• i.

supp$(W_{j})\subseteq\Omega_{j}$,

• ii.

$\sum_{j=1}^{d}W_{j}(x,y)=1$ on $\Omega$,

• iii.

$||D^{\beta}W_{j}||_{L^{\infty}(\Omega_{j})}\leq\frac{C_{\beta}}{\xi_{j}^{|\beta|}},$ $\forall\beta\in\mathbb{N}^{2}:|\beta|\leq p,$ where $C_{\beta}>0$ is a constant and $\xi_{j}=\textrm{diam}(\Omega_{j})$.

Then, the global interpolant assumes the form

[TABLE]

where $R_{j}$ defines a local interpolant on each PU subdomain $\Omega_{j}$ and $W_{j}:\Omega_{j}\longrightarrow\mathbb{R}$.

Remark 2.1**.**

Since the functions $W_{j}$, $j=1,\ldots,d$, form a partition of unity, if the local fits $R_{j}$, $j=1,\ldots,d$, satisfy the interpolation conditions, then the global PU approximant inherits the interpolation property, i.e.

[TABLE]

where $I(x_{i},y_{i})=\{j/(x_{i},y_{i})\in\Omega_{j}\}.$

In what follows, the local interpolants $R_{j}$ are defined as linear combinations of RBFs centred at $(x_{k},y_{k})$, with $(x_{k},y_{k})\in{\cal X}_{N_{j}}={\cal X}_{N}\cap\Omega_{j}$, for $k=1,\ldots,N_{j}$, i.e.

[TABLE]

where $N_{j}$ is the number of nodes on $\Omega_{j}$ and $\phi_{\varepsilon}^{(k)}(x,y)$ is a RBF centred at $(x_{k},y_{k})$. A RBF depends on the so-called shape parameter $\varepsilon$ which governs the flatness of the function.

Among a large variety of RBFs, we can differentiate between compactly supported and globally defined RBFs. As examples of these two major classes, we consider the Wendland’s $C^{2}$ function and the Inverse MultiQuadric (IMQ). The latter is strictly positive definite and given by

[TABLE]

where $r_{k}^{2}(x,y)=(x-x_{k})^{2}+(y-y_{k})^{2}$ is the square of the Euclidean distance from the centre $(x_{k},y_{k})$. Moreover, in case of ill-conditioning of the interpolation matrix, it might be advantageous to work with locally supported functions. Indeed, they lead to sparse linear systems. Wendland found a class of RBFs which are smooth, locally supported and strictly positive definite. For example, the Wendland’s $C^{2}$ function is defined as

[TABLE]

where $(\cdot)_{+}$ denotes the truncated power function. Note that the shape parameter for a CSRBF identifies the support of the function.

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

[TABLE]

where $(x_{i},y_{i})\in{\cal X}_{N_{j}}$, and $f_{i}$ are the associated function values. Thus, the problem of finding the PU interpolant (1) requires to solve $d$ linear systems of the form

[TABLE]

where $\boldsymbol{c}=\left(c_{1},\ldots,c_{N_{j}}\right)^{T}$, $\boldsymbol{f}=\left(f_{1},\ldots,f_{N_{j}}\right)^{T}$ and the entries of the matrix $A^{(j)}\in\mathbb{R}^{N_{j}\times N_{j}}$ are

[TABLE]

If the considered RBF is strictly positive definite, the interpolant (1) is unique [15]. However, even if here we only focus on strictly positive definite RBFs, we remark that the uniqueness of the interpolant can be ensured also for the general case of strictly conditionally positive definite functions by modifying (2), see e.g. [37].

In order to be able to formulate error bounds, we need some further assumptions on the regularity of $\Omega_{j}$ and thus we give the following definitions.

Definition 2.1**.**

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

[TABLE]

Definition 2.2**.**

An open and bounded covering $\{\Omega_{j}\}_{j=1}^{d}$ is called regular for $(\Omega,{\cal X}_{N})$ if the following properties are satisfied [38]:

• i.

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

• ii.

there exist a constant $C_{r}>0$ and an angle $\theta\in(0,\pi/2)$ such that every subdomain $\Omega_{j}$ satisfies an interior cone condition (with angle $\theta$ and radius $C_{r}h_{{\cal X}_{N},\Omega}$),

• 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}$.

Thus, after defining the space $C_{\nu}^{p}(\mathbb{R}^{2})$ of all functions $f\in C^{p}$ whose derivatives of order $|\beta|=p$ satisfy $D^{\beta}f(x,y)={\cal O}((\sqrt{x^{2}+y^{2}})^{\nu})$ for $\sqrt{x^{2}+y^{2}}\longrightarrow 0$, we consider the following convergence result [15, 38].

Theorem 2.1**.**

Let $\Omega\subseteq\mathbb{R}^{2}$ be open and bounded and suppose that ${\cal X}_{N}=\{(x_{i},y_{i}),i=1,\ldots,N\}\subseteq\Omega$. Let us consider a strictly conditionally positive definite RBF which belongs to $C_{\nu}^{p}(\mathbb{R}^{2})$. Let $\{\Omega_{j}\}_{j=1}^{d}$ be a regular covering for $(\Omega,{\cal X}_{N})$ and let $\{W_{j}\}_{j=1}^{d}$ be $p$-stable for $\{\Omega_{j}\}_{j=1}^{d}$. Then the error between $f\in\mathscr{N}_{\phi}(\Omega)$, where $\mathscr{N}_{\phi}$ is the native space of the basis function, and its PU interpolant (1) can be bounded by

[TABLE]

for all $(x,y)\in\Omega$ and all $|\beta|\leq p/2$.

2.2 The positive constrained partition of unity approximant

If we compare the result reported in Theorem 2.1 with the global error estimate shown in [15, 37], we can see that the PU interpolant preserves the local approximation order for the global fit. In particular, the PU method can be thought as the Shepard’s method where $R_{j}$ are used instead of the data values $f_{j}$. Even if the classical Shepard’s approximant in its original formulation is overcome [29], it possesses a useful property; specifically, it lies within the range of the data. As a consequence, it is positive if the data values are positive [20]. The positivity-preserving property does not hold in its quadratic formulation nor for the PU approximant presented in the previous section.

In order to avoid such drawback for the PU method, we can directly act on the local RBF interpolants, following the strategy proposed in [39]. In such paper, a scheme devoted to construct a positive global fit is performed by defining several constraints. Extensive results show the good performances of such approach. In fact, the fit of positive samples is always positive, but, in comparison with the original unconstrained interpolation, a degrade of the quality of the approximation is observed. This is mainly due to the fact that a global method is considered and thus the shape of the surface is consequently globally modified (and not only in the area in which the interpolant is negative). Therefore, here we propose a new formulation for a positive fit considering (1).

Sufficient condition to have positive approximants on each subdomain $\Omega_{j}$ is that the coefficients $c_{k}$ of (2) are all positive. To such scope, following [39], at first we choose ${\hat{N}}_{j}$ added data

[TABLE]

on the subdomain $\Omega_{j}$. Then, the $j$-th approximation problem consists in finding a function ${\hat{R}}_{j}$ of the form

[TABLE]

such that

[TABLE]

where $\hat{\phi}_{\varepsilon_{{\hat{k}}}}^{({\hat{k}})}$ are CSRBFs.

Note that in (5) we consider different supports for the CSRBFs. In particular, if a constraint $(\hat{x}_{i},\hat{y}_{i})$ is added on $\Omega_{j}$ in a neighborhood of a point, namely $(x_{i},y_{i})$, we select $\varepsilon_{{\hat{i}}}$ such that only $(x_{i},y_{i})$ belongs to the support of the CSRBF. This choice is due to the fact that, doing in this way, at least when $\hat{N}_{j}=N_{j}$, the problem (5) subject to (6) admits solution. This is proved in [39] by using the Gordan’s Theorem [5]. A brief sketch of the proof will be given in what follows because important considerations arise. Let us define

[TABLE]

where $b_{{\hat{k}}}$ are positive real numbers. Then, since a vector $\boldsymbol{v}$ such that $\boldsymbol{a}_{i}^{T}\boldsymbol{v}<0$, $i=1,\ldots,2N_{j}+1$, does not exist, from Gordan’s Theorem, we know that there exist nonnegative real numbers $c_{1},\ldots,c_{2N_{j}+1}$, such that $\sum_{i=1}^{2N_{j}+1}c_{i}\boldsymbol{a}_{i}=0$ and $c_{2N_{j}+1}>0$.

In practical use, the aim is to add as few data as possible, therefore, according to [39], one can find the minimum of $\sum_{i=N_{j}+1}^{2N_{j}}\textrm{sign}(c_{i})$, such that (6) is satisfied. However, this approach does not guarantee an optimal solution in terms of accuracy. Here instead, with a technique described in the next section, we select the optimal number of added data $\hat{N}_{j}$ which yields maximal accuracy. Nevertheless, our aim consists in perturbing (2) with small quantities. A feasible way is to find

[TABLE]

such that (6) is satisfied.

The local approach here proposed enables us to modify the shape of the surface only if a negative local approximant is found, in fact if the $j$-th original local fit of the form (2) is positive, we do not need to add other data and in this case ${\hat{R}}_{j}=R_{j}$ on $\Omega_{j}$. Therefore, for each subdomain, after selecting a suitable number of constraints $\hat{N}_{j}$, which can also be 0, the PC-PU approximant assumes the form

[TABLE]

Remark 2.2**.**

In [39] the authors limit their attention to CSRBFs. Here instead we will couple them with globally defined RBFs. In fact, even if specific supports of compactly supported kernels must be associated to the added constraints (see the second term in the right-hand side of (5)), we do not have any restrictions on the first term in the right-hand side of (5). Thus, we can use different types of RBFs. In what follows, we will point out that coupling RBFs and CSRBFs leads to a benefit in terms of accuracy.

3 The PC-PU approximant: algorithm and suitable selection of positive constraints

This section is devoted to describe the PC-PU algorithm. In particular, we focus on the suitable choice of the number of positive constraints. Even if the technique here discussed is robust enough to work in any domain $\Omega\subseteq\mathbb{R}^{2}$, for simplicity we consider $\Omega=[0,1]^{2}$.

3.1 Selection of positive constraints

Acting as explained in the previous section, we can ensure the positivity of the PU approximant. However, depending on the number of positive constraints, this might lead to a low accuracy. In [39], for a global RBF-based interpolant, $h$ random data are selected as constraints if $f_{h}$ is the minimum among all values $f_{i}$, $i=1,\ldots,N$. This criterion does not guarantee a good accuracy neither the existence of a solution and thus we design an alternative approach enabling us to select a suitable number of positive constraints.

The proposed method is based on an a priori error estimate. To this aim, several schemes have already been developed. Precisely, we focus on the so-called cross validation algorithm, see [15, 19]. A variant of such method, known in literature as Leave One Out Cross Validation (LOOCV), is detailed in [26]. This approach is always used to find the optimal value of the shape parameter of the basis function. Here instead, for each PU subdomain we are interested in selecting a suitable number of constraints $\hat{N}_{j}$.

To avoid complexities, let us first consider an interpolation problem on $\Omega_{j}$ of the form (2). Moreover, let us define the $j$-th interpolant obtained leaving out the $i$-th data on $\Omega_{j}$ as

[TABLE]

and 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 $(e_{1},\ldots,e_{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.

This implementation is computationally expensive. In fact, the matrix inverse, which requires ${\cal O}(N_{j}^{3})$ operations, must be computed for each node. This leads to a total computational cost of ${\cal O}(N_{j}^{4})$ operations. Thus, following [15, 26], we simplify the computation to a single formula. Precisely, we calculate

[TABLE]

where $c_{i}$ is the $i$-th coefficient of the interpolant based on the full data set and $(A^{(j)}_{ii})^{-1}$ is the $i$-th diagonal element of the inverse of the corresponding local interpolation matrix.

In our case, in order to guarantee a positive fit, we deal with an augmented local problem and therefore, as error estimate, we compute the following quantity

[TABLE]

where the symmetric matrix $\hat{A}^{(j)}$ is defined as

[TABLE]

Moreover, in order to stress the dependence of the error on $\hat{N}_{j}$, we use the notation

[TABLE]

Note that in our case the coefficients are not determined by directly computing the solution of a linear system, but they are found out by solving (7), subject to (6). Thus, to be more precise, we should refer to this criterion as LOOCV-like method, but in order to keep common notations, we will go on calling it simply LOOCV. Indeed, we are able to fix a criterion which enables us to select a suitable number of positive constraints. Specifically, focusing on the maximum norm, we compute (10) for $\hat{N}_{j}=1,\ldots,N_{j}$. Thus, if on $\Omega_{j}$ a negative fit is observed, we add $\hat{N}_{j}$ positive constraints if

[TABLE]

and the fit is positive. It is easy to see that we automatically ensure that the conditions (6) are satisfied, indeed they are fulfilled at least for ${\hat{N}}_{j}=N_{j}$.

The ${\hat{N}}_{j}$ new added data can be placed randomly within the $j$-th patch, but numerically we observe that selecting well distributed points on $\Omega_{j}$ leads to a better accuracy. As a consequence, on the subdomain $\Omega_{j}$ of centre $(\bar{x}_{j},\bar{y}_{j})$ and radius $\delta$, we consider ${\hat{N}}_{j}$ positive constraints, distributed as the seeds on a sunflower head, i.e. defined as [32, 35]

[TABLE]

$k=N_{j}+1,\ldots,N_{j}+{\hat{N}}_{j}$, where

[TABLE]

We conclude this section with an illustrative figure. A 2D view of a partition of unity structure covering a set of scattered data in the unit square is shown in Figure 2 (left); in the right frame we plot a set of points computed with (14).

3.2 Description of the algorithm

To make simpler the presentation, we summarize in the PC-PU Algorithm the steps needed to compute the PC-PU approximant.

At first, given a set of scattered data in $\Omega=[0,1]^{2}$, we construct the PU subdomains. They are circular patches centred at a grid of points ${\cal C}_{d}=\{(\bar{x}_{i},\bar{y}_{i}),i=1,\ldots,d\}\subseteq\Omega$ of radius

[TABLE]

According to [15], we choose the number of subdomains such that $N/d\approx 4$. The PC-PU approximant is then computed at a grid of $s$ evaluation points ${\cal E}_{s}=\{(\tilde{x}_{i},\tilde{y}_{i}),i=1,\ldots,s\}$, see Steps 1-3 of the PC-PU Algorithm.

Once the partition of unity structure is built, a local interpolation problem needs to be solved for each PU subdomain. Specifically, in the $j$-th local approximation problem, only those data sites and evaluation points belonging to $\Omega_{j}$ are involved, see Steps 4-6 of the PC-PU Algorithm. To find such points the so-called kd-tree partitioning structures are commonly and widely used [3, 15]. However, here we use the efficient partitioning structure proposed in [11]. It is a multidimensional procedure which leads to a saving in terms of computational time with respect to the previous partitioning procedures proposed in [9, 10].

Once both the sets ${\cal X}_{N}$ and ${\cal E}_{s}$ are organized among the different patches, an interpolant of the form (2) is constructed, see Step 7 of the PC-PU Algorithm. Then, only if the local fit is negative we add $\hat{N}_{j}$ positive constraints, as explained in the previous subsection, see Steps 8-9 of the PC-PU Algorithm. In this way, for each PU subdomain a positive local RBF approximant is computed for each evaluation point. Finally, the global fit (8) is applied accumulating all the $\hat{R}_{j}$ and $W_{j}$, see Steps 10-11 of the PC-PU Algorithm.

4 Numerical experiments

This section is devoted to show, by means of extensive numerical simulations, the performances of the PC-PU approximant. In [39] the authors point out that the global CSRBF constrained method possesses a better approximation behaviour than the CMQS approximant [1, 6].

Here, we compare our PC-PU fit with the original Shepard’s method. In fact, for this method the positivity-preserving property holds, while it is lost in the modified quadratic version [20]. Then, comparisons with the global method proposed in [39] and with the classical PU interpolant will be carried out. Obviously, since we perform a PU approximation, large data sets are considered. On the opposite, a global interpolant, such as the one outlined in [39], cannot handle large sets.

Experiments are performed considering several sets of random nodes contained in the unit square $\Omega=[0,1]^{2}$, a grid of $d=\lfloor\sqrt{N}/2\rfloor^{2}$ subdomain centres and a grid of $s=80\times 80$ evaluation points.

In order to test the accuracy of the proposed method, we compute, for different kernels with different order of smoothness, the Maximum Absolute Error (MAE) and the Root Mean Square Error (RMSE) whose formulas are

[TABLE]

and

[TABLE]

The errors are computed using as test functions

[TABLE]

and

[TABLE]

Experiments are carried out considering the Wendland’s $C^{2}$ function for the added nodes in (5). Furthermore, since for the given interpolation conditions the choice is arbitrary, we use both the Wendland’s $C^{2}$ and the IMQ $C^{\infty}$ functions. For them, we have to fix the shape parameters. We remark that the results are affected by the choice of the $\varepsilon$, in fact small values of $\varepsilon$ lead to problems of instability, while for large values the approximate solution might be inaccurate. In particular, many researchers already worked on the problem of finding stable approximations when $\varepsilon$ tends to zero (see e.g. [13, 14, 16]). Thus, referring to such papers one can easily guess which are the safe values. Here, for instance, we fix the shape parameter $\varepsilon$ equal to $0.1$ and $1$ for the Wendland’s $C^{2}$ and the IMQ $C^{\infty}$ functions, respectively. Moreover, we point out that our aim is to compare the behaviour of the PC-PU method versus the classical PU approach. In this sense, such comparison is independent from the shape parameter, indeed we always register a decrease of the accuracy when the PC-PU method is used.

Tables 2 and 3 show a direct comparison between the classical PU interpolant, which leads to a negative fit, and the PC-PU approximant. We can observe a better behaviour of the PC-PU approximant when globally supported RBFs are used instead of CSRBFs. In fact, with the latter, large errors are introduced in the region where a negative fit is observed. Roughly speaking, using globally defined RBFs, the errors of the classical PU method and of the PC-PU approach are close to each other. On the opposite, by means of CSRBFs, the PC-PU approximant preserves the positivity property with an error which is about two times the one of the unconstrained interpolant. In order to have a graphical prove, refer to Figures 3 and 4 in which we consider $N=3500$ nodes and the test function $f_{1}$. From these figures we can note that, both with the use of CSRBFs and globally defined RBFs there is no smoothing effect when the solution approaches zero, i.e. when constraints are used.

Then, in Tables 4 and 5 we show the errors obtained with the Shepard’s method and the CSRBF-based global method proposed in [39]. For this method we consider the Wendland’s $C^{2}$ kernel as CSRBF. As evident from Figures 5 and 6, in which we consider $N=3500$ random nodes and the test function $f_{2}$, our proposed local scheme maintains a better accuracy than the other considered methods. Moreover, from Figure 5, we can again observe a better behaviour of the PC-PU method coupled with globally defined RBFs rather than with CSRBFs. This behaviour, is also due to the fact that with CSRBFs we need to add more constraints. In particular, if on $\Omega_{j}$ a negative fit is observed, the PC-PU technique selects ${\hat{N}}_{j}$ constraints, computed with (14), such that the positivity of the local interpolant is ensured. In general, we note that when globally defined RBFs are used the number ${\hat{N}}_{j}$ is so that ${\hat{N}}_{j}\ll N_{j}$. On the opposite with CSRBFs, ${\hat{N}}_{j}$ is usually closer to $N_{j}$. In other words, we need to add more constraints in case of CSRBFs and, as stressed in [39], this causes a decrease of the fit accuracy.

5 Application to population dynamics

In this section we point out how the PC-PU procedure can be useful in applied sciences. Among several applications which can be investigated, we focus on a mathematical model dealing with wild herbivores in forests.

Let $H$, $G$ and $T$ represent respectively the herbivores, grass and trees populations of the environment in consideration. The model we consider is a classical predator with two prey system [34], in which the resources are consumed following a concave response function, usually called the Beddington-De Angelis function [12]. Let $\alpha$ and $\beta$ be the inverse of the herbivores maximal consumption of grass and trees, respectively. Let $r_{1}$ and $r_{2}$ denote the grass and trees growth rates, $K_{1}$ and $K_{2}$ their respective carrying capacities, $\mu$ the metabolic rate of herbivores, $c$ and $g$ the half saturation constants, $e\leq 1$ and $f\leq 1$ the conversion factors of food into new herbivore biomass and $a$ and $b$ the daily feeding rates due to grass and trees, respectively, the model reads as follows [28]

[TABLE]

All parameters are nonnegative. In particular, $K_{1}$, $K_{2}$, $c$ and $g$ are measured in biomass, $e$, $f$, $\alpha$ and $\beta$ are pure numbers, $\mu$, $r_{1}$, $r_{2}$, $a$ and $b$ are rates.

For the study of critical points and their stability, refer to [34]. The equilibrium that play a role in this investigation is the coexistence equilibrium point $E^{*}=(H^{*},G^{*},T^{*})$. It can be assessed only via numerical simulations. We test the model considering the natural park of the Dolomiti Bellunesi located in the Northern Italy. To estimate the parameters, we refer to [17, 34]. Basing our considerations on data related to the number of herbivores, the extension of the park and following tables providing the estimation of the annual net primary production of several environments shown in [17], we can fix $\mu=0.03$, $r_{1}=0.01$, $r_{2}=0.0006$, $\alpha=0.05^{-1}$, $\beta=8$, $e=0.605$, $f=0.001$, $K_{1}={3469640.64}$, $K_{2}={15695993.39}$, $c={101862.16}$ and $g={1001229580.18}$. Moreover, also the initial conditions $H(0)=268.750$, $G(0)=2313093.76$ and $T(0)=1046399.56$ are known.

By means of numerical simulations carried out in [28], we can observe that the herbivore population level appears to be very sensitive to small perturbations of several parameters and under the high risk of extinction. The results indicate that the parameters most affecting the system’s final configuration are $\alpha$ and $\mu$. Such consideration follows from the surface plotted in Figure 7; the surface shows the value attained by the herbivores in function of the parameters $\alpha$ and $\mu$. The dot represents the situation in the present ecosystem conditions. As we can graphically note from the surface, for small perturbations on the current values of $\alpha$ and $\mu$, the herbivores can extinguish, i.e. $H=0$.

In Figure 7 the surfaces are reconstructed with the standard PU interpolant and with the PC-PU approximant, left and right respectively. Even if it is not evident because of the large scale of the $z$-axis, the surface reconstructed with the PU interpolant reaches its minimum at $-1$, violating the biological constraint $H\geq 0$. From this consideration, the importance of having a robust tool enabling us to fit positive data values with a positive approximant follows.

6 Concluding remarks

In this paper we presented a robust technique devoted to preserve the positivity of the fit via a local RBF-based method. Extensive numerical simulations have been carried out to show the effectiveness of the method. Moreover, in order to motivate the importance of a meshless method that do not violate biological or physical constraints, we investigated an application in biomathematics.

Work in progress consists in extending the proposed tool in higher dimensions. Furthermore, in case of ill-conditioning, since in this context CSRBFs are less performing than globally defined RBFs, further investigations in the use of stable bases [13, 16] are needed.