Positive constrained approximation via RBF-based partition of unity method

0.72 MB·

by  Alessandra De Rossi

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
7 Abstract. Inthispaper,wediscusstheproblemofconstructingRadialBa-
1 sisFunction(RBF)-basedPartitionofUnity(PU)interpolantsthatarepositive
0
if data values are positive. More specifically, we compute positive localapprox-
2
imants by adding up several constraints to the interpolation conditions. This

n
approach, considering a global approximation problem and Compactly Sup-
a
ported RBFs (CSRBFs), has been previouslyproposedin [39]. Here, the use of
J
the PU technique enables us to intervene only locally and as a consequence to
5
2 reach a better accuracy. This is also due to the fact that we select the optimal
numberofpositiveconstraintsbymeansofana priori errorestimateandwedo

] not restrict to the use of CSRBFs. Numerical experiments and applications to
A
population dynamics are provided to illustrate the effectiveness of the method
N
in applied sciences.
.
h
t
a 1 Introduction
m
[ Givenasetofsamples,thescattereddatainterpolationproblemconsistsinfind-

ing an approximating function that matches the given measurements at their
1
correspondinglocations. Furthermore,dealing with applications, we often have
v
0 additional properties, such as the non-negativity of the measurements, which
6 we wish to be preservedduring the interpolationprocess. Note that, since such
2 property is known as positivity-preserving property, to keep a common nota-
7
tion with existing literature, we use the term positive (instead of non-negative)
0
function values or interpolants.
.
1
To preservesuch property,mostly consideringrationalspline functions with
0
C1 or C2 continuity, the recent research studies techniques which force the
7
1 approximants to be positive. As example, the conditions under which the pos-
: itivity of a cubic piece may be lost are investigated in [30]. Moreover, in order
v
to preserve the positivity, the use of bicubic splines, coupled with a technique
i
X basedonaddingextraknots,hasbeeninvestigatedin[2,8]. Apositivefit isin-
r steadobtained by means of rationalcubic splines in [21, 23]. The same authors
a
alsodevelopedapositivesurfaceconstructionschemeforpositivescattereddata
arrangedover 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,
evenifsuchmeshfreeapproachhasbeenextensivelystudiedintherecentyears,
especially focusing on the stability of the interpolant [13, 16], not a lot of ef-
fort has been addressed to construct positivity-preserving approximants. Such
1
problemhasbeen studiedonly inparticularandwell-knowncases;for instance,
in[33]itisanalyzedforthethinplatespline. Othermethods,whichfollowfrom
thequasi-interpolationformulagivenin[40],havebeenproposedandeffectively
performed to construct positive approximants of positive samples [36].
Focusing onRBFs, our scope consists in preservingthe positivity ofthe PU
interpolantfor a wider family of kernels. In[39]a globalpositive RBFapproxi-
mantisconstructedbyaddinguptotheinterpolationconditionsseveralpositive
constraintsandconsideringCSRBFs. Eveniftheoptimalnumberofconstraints
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
thequalityoftheapproximatingfunctionincomparisonwiththeunconstrained
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,incontrastwith[39],we canconsidertrulylargedata sets.
Specifically,inordertoconstructthe PositiveConstrainedPU(PC-PU)ap-
proximant,following [39], we locally impose severalpositive constraints and we
reducetosolveanoptimizationproblem. Thenumberofconstraintsisproperly
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 dis-
advantage 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 direc-
tion. Moreover,in order to stress the effectiveness of the proposedtechnique in
appliedsciences,weinvestigateseveralapplicationstobiomathematics. Finally,
comparisonswiththeunconstrainedPUinterpolant,withtheShepard’smethod
and with the one outlined in [39] are carried out.
The guidelines of the paper are as follows. In Section 2 we investigate the
positivityofthePUapproximantbyconsideringextrapositiveconstraints. 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
Inadditiontothe assignedinterpolationconditions,ourgoalconsistsinconsid-
ering 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 physi-
cal measurements, a positive fit is often necessary. Thus, before discussing the
2
proposedmethodandwithoutloosinggenerality,weconsideranexampleofuni-
variateinterpolationtoillustratethe scope ofsuchnumericaltool andmotivate
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 inter-
polant. Anyway, thepositivity ofthePSAvalues is notalways 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.
10
20
8
mL) mL) 15
g/ 6 g/
n n
A ( A ( 10
S 4 S
P P
2 5
0 0
0 10 20 30 40 50 0 10 20 30 40 50
t (months) t (months)
Figure 1: Examples of curves fitting PSA values (plotted with blue dots). The
classical PU interpolant is plotted in orange and the PC-PU approximant in
blue.
2.1 The partition of unity method
LetusconsideraboundedsetΩ R2. GivenasetofN distinctdatapoints(or
⊆
data sites or nodes) = (x ,y ),i =1,...,N Ω and a set of data values
N i i
X { }⊆
(or measurements or function values) = f = f(x ,y ),i = 1,...,N , the
N i i i
F { }
PUmethoddecomposesthe domainΩ into dsubdomainsΩ , j =1,...,d, such
j
thatΩ d Ω [4,15,24,38]. Inliterature,thesubdomainsΩ aresupposed
⊆ j=1 j j
to be circuSlar patches of the same radius δ covering the domain Ω.
According to [38], associated with these subdomains, we choose a p-stable
partition of unity, i.e a family of nonnegative functions W d , with W
{ j}j=1 j ∈
Cp(R2), such that
i. supp(W ) Ω ,
j j
⊆
3
ii. d W (x,y)=1 on Ω,
j=1 j
P
iii. DβW Cβ , β N2 : β p, where C >0 is a constant and
|| j||L∞(Ωj) ≤ ξ|β| ∀ ∈ | |≤ β
j
ξ =diam(Ω ).
j j
Then, the global interpolant assumes the form
d
(x,y)= R (x,y)W (x,y), (x,y) Ω, (1)
j j
I ∈
Xj=1
whereR definesalocalinterpolantoneachPUsubdomainΩ andW :Ω
j j j j
R. −→
Remark 2.1. Since the functions W , j = 1,...,d, form a partition of unity,
j
if the local fits R , j = 1,...,d, satisfy the interpolation conditions, then the
j
global PU approximant inherits the interpolation property, i.e.
d
(x ,y )= R (x ,y )W (x ,y )= f W (x ,y )=f ,
i i j i i j i i i j i i i
I
Xj=1 j∈IX(xi,yi)
where I(x ,y )= j/(x ,y ) Ω .
i i i i j
{ ∈ }
Inwhatfollows,the localinterpolantsR aredefined aslinearcombinations
j
of RBFs centred at (x ,y ), with (x ,y ) = Ω , for k =1,...,N ,
k k k k ∈XNj XN ∩ j j
i.e.
Nj
R (x,y)= c φ(k)(x,y), (x,y) Ω , (2)
j k ε ∈ j
Xk=1
where N is the number of nodes on Ω and φ(k)(x,y) is a RBF centred at
j j ε
(x ,y ). A RBFdepends on the so-calledshape parameterε which governsthe
k k
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 C2 function and the Inverse MultiQuadric (IMQ).
The latter is strictly positive definite and given by
1
φ(k)(x,y)= , (3)
ε 1+ε2r2(x,y)
k
p
wherer2(x,y)=(x x )2+(y y )2isthesquareoftheEuclideandistancefrom
k − k − k
the centre (x ,y ). Moreover, in case of ill-conditioning of the interpolation
k k
matrix, it might be advantageous to work with locally supported functions.
Indeed, they lead to sparse linear systems. Wendland found a class of RBFs
whicharesmooth,locallysupportedandstrictlypositivedefinite. Forexample,
the Wendland’s C2 function is defined as
φ(k)(x,y)=(1 εr (x,y))4(4εr (x,y)+1), (4)
ε − k + k
where() denotesthetruncatedpowerfunction. Notethattheshapeparameter
+
·
for a CSRBF identifies the support of the function.
4
The coefficients c Nj in (2) are determined by enforcing the N local
{ k}k=1 j
interpolation conditions
R (x ,y )=f , i=1,...,N ,
j i i i j
where (x ,y ) , and f are the associated function values. Thus, the
i i ∈ XNj i
problem of finding the PU interpolant (1) requires to solve d linear systems of
the form
A(j)c=f,
T T
where c = c ,...,c , f = f ,...,f and the entries of the matrix
1 Nj 1 Nj
A(j) RNj×(cid:0)Nj are (cid:1) (cid:0) (cid:1)
∈
A(j) =φ(k)(x ,y ), i,k =1,...,N .
ik ε i i j
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,weremarkthattheuniquenessoftheinterpolantcanbe ensuredalsofor
thegeneralcaseofstrictlyconditionallypositivedefinitefunctionsbymodifying
(2), see e.g. [37].
Inordertobeabletoformulateerrorbounds,weneedsomefurtherassump-
tions on the regularity of Ω and thus we give the following definitions.
j
Definition 2.1. The fill distance, which is a measure of data distribution, is
given by
h = sup min r (x,y) .
XN,Ω (x,y)∈Ω(cid:18)(xi,yi)∈XN i (cid:19)
Definition 2.2. An open and bounded covering Ω d is called regular for
{ j}j=1
(Ω, ) if the following properties are satisfied [38]:
N
X
i. for each (x,y) Ω, the number of subdomains Ω with (x,y) Ω is
j j
∈ ∈
bounded by a global constant C,
ii. there exist a constant C > 0 and an angle θ (0,π/2) such that every
r
∈
subdomain Ω satisfies an interior cone condition (with angle θ and radius
j
C h ),
r XN,Ω
iii. the local fill distances h are uniformly bounded by the global fill dis-
XNj,Ωj
tance h .
XN,Ω
Thus, after defining the space Cp(R2) of all functions f Cp whose deriva-
ν ∈
tives of order β = p satisfy Dβf(x,y)= (( x2+y2)ν) for x2+y2 0,
| | O −→
we consider the following convergence result [1p5, 38]. p
Theorem 2.1. Let Ω R2 be open and bounded and suppose that =
N
⊆ X
(x ,y ),i = 1,...,N Ω. Let us consider a strictly conditionally posi-
i i
{tive definite RBF whic}h b⊆elongs to Cp(R2). Let Ω d be a regular covering
ν { j}j=1
for (Ω, ) and let W d be p-stable for Ω d . Then the error between
f NX(NΩ), where N{ j}isj=t1he native space o{f tjh}ej=b1asis function, and its PU
φ φ
∈
interpolant (1) can be bounded by
Dβf(x,y) Dβ (x,y) C′hp+2ν−|β| f ,
| − I |≤ XN,Ω | |Nφ(Ω)
for all (x,y) Ω and all β p/2.
∈ | |≤
5
2.2 The positive constrained partition of unity approxi-
mant
Ifwecomparethe resultreportedinTheorem2.1withtheglobalerrorestimate
shownin[15,37],wecanseethatthePUinterpolantpreservesthelocalapprox-
imation order for the global fit. In particular, the PU method can be thought
as the Shepard’s methodwhere R areusedinsteadofthe data valuesf . Even
j j
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-preservingpropertydoesnotholdinitsquadraticformulationnorfor
the PU approximant presented in the previous section.
In order to avoid such drawback for the PU method, we can directly act on
the localRBF interpolants, followingthe strategyproposedin [39]. In such pa-
per,aschemedevotedto constructapositiveglobalfitisperformedbydefining
several constraints. Extensive results show the good performances of such ap-
proach. Infact,thefitofpositivesamplesisalwayspositive,but,incomparison
with the original unconstrained interpolation, a degrade of the quality of the
approximationis observed. This is mainly due to the factthat a globalmethod
isconsideredandthustheshapeofthesurfaceisconsequentlygloballymodified
(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 Ω is
j
that the coefficients c of (2) are all positive. To such scope, following [39], at
k
first we choose Nˆ added data
j
(xˆ ,yˆ ),...,(xˆ ,yˆ ),
Nj+1 Nj+1 Nj+Nˆj Nj+Nˆj
onthesubdomainΩ . Then,thej-thapproximationproblemconsistsinfinding
j
a function Rˆ of the form
j
Nj Nj+Nˆj
Rˆ (x,y)= c φ(k)(x,y)+ c φˆ(kˆ)(x,y), (x,y) Ω , (5)
j k ε kˆ εkˆ ∈ j
Xk=1 kˆ=XNj+1
such that
Rˆ (x ,y )=f , i=1,...,N , c 0, i=1,...,N +Nˆ , (6)
j i i i j i j j
≥
where φˆ(kˆ) are CSRBFs.
εkˆ
Note that in (5) we consider different supports for the CSRBFs. In particu-
lar, if a constraint (xˆ ,yˆ) is added on Ω in a neighborhood of a point, namely
i i j
(x ,y ),weselectε suchthatonly(x ,y )belongstothesupportoftheCSRBF.
i i ˆi i i
This choice is due to the fact that, doing in this way, at least when Nˆ = N ,
j 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
aT = φ(k)(x ,y ) Nj , k =1,...,N ,
k −{ ε i i }i=1 j
aT = φˆ(kˆ)(x ,y ) Nj =(0,...,0, b ,0,...,0), kˆ =N +1,...,2N ,
kˆ −{ εkˆ i i }i=1 − kˆ j j
aT =(f ,...,f ),
2Nj+1 1 Nj
6
where b are positive real numbers. Then, since a vector v such that aTv <0,
kˆ i
i=1,...,2N +1, does not exist, fromGordan’s Theorem, we know that there
j
exist nonnegative real numbers c ,...,c , such that 2Nj+1c a = 0 and
1 2Nj+1 i=1 i i
c2Nj+1 >0. P
Inpracticaluse,theaimistoaddasfewdataaspossible,therefore,according
to[39],onecanfindtheminimumof 2Nj sign(c ),suchthat(6)issatisfied.
i=Nj+1 i
However, this approach does not gPuarantee 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 Nˆ which yields maximal accuracy.
j
Nevertheless,ouraimconsistsinperturbing(2)withsmall quantities. Afeasible
way is to find
1/2
Nj+Nˆj
min c2 , (7)
i
i=XNj+1 
such that (6) is satisfied.
The local approach here proposed enables us to modify the shape of the
surface only if a negative local approximantis found, in fact if the j-th original
localfitoftheform(2)ispositive,wedonotneedtoaddotherdataandinthis
case Rˆ = R on Ω . Therefore, for each subdomain, after selecting a suitable
j j j
numberofconstraintsNˆ ,whichcanalsobe0,thePC-PUapproximantassumes
j
the form
d Nj Nj+Nˆj
Iˆ(x,y)=  ckφ(εk)(x,y)+ ckˆφˆε(kˆkˆ)(x,y)Wj(x,y). (8)
Xj=1Xk=1 kˆ=XNj+1 
Remark2.2. In[39]theauthorslimittheirattentiontoCSRBFs. Hereinstead
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 suit-
able 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 Ω R2, for
⊆
simplicity we consider Ω=[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,
thismightleadtoalowaccuracy. In[39],foraglobalRBF-basedinterpolant,h
7
random data are selected as constraints if f is the minimum among all values
h
f , i=1,...,N. This criterion does not guarantee a good accuracy neither the
i
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-
calledcross validation algorithm,see[15,19]. Avariantofsuchmethod,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
ofthe basisfunction. Here instead,foreachPUsubdomainweareinterestedin
selecting a suitable number of constraints Nˆ .
j
To avoid complexities, let us first consider an interpolation problem on Ω
j
of the form (2). Moreover, let us define the j-th interpolant obtained leaving
out the i-th data on Ω as
j
Nj
Ri(x,y)= c φ(k)(x,y), (x,y) Ω ,
j k ε ∈ j
k=X1,k6=i
and let
e =f Ri(x ,y ),
i i− j i i
betheerroratthei-thpoint. Then,thequalityofthe localfitisdeterminedby
some norm of the vector of errors(e ,...,e )T, obtained by removing in turn
1 Nj
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 in-
verse, which requires (N3) operations, must be computed for each node.
O j
This leads to a total computational cost of (N4) operations. Thus, follow-
O j
ing [15, 26], we simplify the computation to a single formula. Precisely, we
calculate
c
i
e = , (9)
i −1
A(j)
ii
(cid:16) (cid:17)
where c is the i-th coefficient of the interpolant based on the full data set and
i
(A(j))−1 is the i-th diagonal element of the inverse of the corresponding local
ii
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 quan-
tity
eˆ ,...,eˆ = c1 ,..., cNj+Nˆj , (10)
(cid:16) 1 Nj+Nˆj(cid:17) Aˆ(j) −1 Aˆ(j) −1
(cid:16) 11(cid:17) (cid:16) Nj+NˆjNj+Nˆj(cid:17) 
where the symmetric matrix Aˆ(j) is defined as
φ(1)(x ,y ) φˆ(Nj+Nˆj)(x ,y )
 ε 1 1 ··· εNj+Nˆj 1 1 
Aˆ(j) = ... ... ... . (11)
 
φˆ(εNNjj++NNˆˆjj)(x1,y1) ··· φˆε(NNjj++NNˆˆjj)(xˆNj+Nˆj,yˆNj+Nˆj)
8
Moreover, in order to stress the dependence of the error on Nˆ , we use the
j
notation
eˆ (Nˆ )= eˆ ,...,eˆ . (12)
j j (cid:16) 1 Nj+Nˆj(cid:17)
Note that in our case the coefficients are not determined by directly computing
thesolutionofalinearsystem,buttheyarefoundoutbysolving(7),subjectto
(6). Thus, to be more precise, we should refer to this criterion as LOOCV-like
method, but in orderto keepcommonnotations,we will go oncalling 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 Nˆ = 1,...,N . Thus, if on Ω a negative fit is
j j j
observed, we add Nˆ positive constraints if
j
eˆ (Nˆ ) = min eˆ (k) , (13)
j j ∞ j ∞
|| || k=1,...,Nj|| ||
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 Nˆ =N .
j j
The Nˆ new added data can be placed randomly within the j-th patch, but
j
numerically we observe that selecting well distributed points on Ω leads to
j
a better accuracy. As a consequence, on the subdomain Ω of centre (x¯ ,y¯ )
j j j
and radius δ, we consider Nˆ positive constraints, distributed as the seeds on a
j
sunflower head, i.e. defined as [32, 35]
(xˆ ,yˆ )=(x¯ +u cosη ,y¯ +u sinη ), (14)
k k j k k j k k
k =N +1,...,N +Nˆ , where
j j j
k 1/2 4kπ
u =δ − and η = .
k pNˆ 1/2 k 1+√5
j
q −
Weconcludethissectionwithanillustrativefigure. A2Dviewofapartition
ofunity structurecoveringasetofscattereddatainthe unitsquareisshownin
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.
Atfirst,givenasetofscattereddatainΩ=[0,1]2,weconstructthePUsub-
domains. Theyarecircularpatchescentredatagridofpoints = (x¯ ,y¯),i=
d i i
C {
1,...,d Ω of radius
}⊆
1
δ = . (15)
rd
According to [15], we choose the number of subdomains such that N/d 4.
≈
The PC-PU approximant is then computed at a grid of s evaluation points
= (x˜ ,y˜),i=1,...,s , see Steps 1-3 of the PC-PU Algorithm.
s i i
E { }
Once the partition of unity structure is built, a local interpolation problem
needstobesolvedforeachPUsubdomain. Specifically,inthej-thlocalapprox-
imation problem, only those data sites and evaluation points belonging to Ω
j
9
1 0.5
0.8
0.4
0.6
y y 0.3
0.4
0.2
0.2
00 0.2 0.4 0.6 0.8 1 0.10.1 0.2 0.3 0.4 0.5
x x
Figure 2: Left: an illustrative example of PU subdomains covering the domain
Ω = [0,1]2. The red dots represent a set of scattered data and the blue circles
identify the PU subdomains. Right: an illustrative example of added points
computed with (14).
areinvolved,see Steps 4-6ofthePC-PU Algorithm. Tofindsuchpointsthe
so-called kd-tree partitioning structures are commonly and widely used [3, 15].
However,hereweusetheefficientpartitioningstructureproposedin[11]. Itisa
multidimensional procedure which leads to a saving in terms of computational
time with respect to the previous partitioning procedures proposed in [9, 10].
Onceboththesets and areorganizedamongthedifferentpatches,an
N s
X E
interpolantoftheform(2)isconstructed,see Step 7ofthePC-PU Algorithm.
Then,onlyifthelocalfitisnegativeweaddNˆ positiveconstraints,asexplained
j
in the previous subsection, see Steps 8-9 of the PC-PU Algorithm. In this
way,for eachPUsubdomaina positive localRBFapproximantis computedfor
each evaluation point. Finally, the global fit (8) is applied accumulating all the
Rˆ and W , see Steps 10-11 of the PC-PU Algorithm.
j j
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-preservingproperty holds, while it is lost in
the modified quadraticversion[20]. Then, comparisonswiththe globalmethod
proposed in [39] and with the classical PU interpolant will be carried out. Ob-
viously, 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 con-
tained in the unit square Ω=[0,1]2, a gridof d= √N/2 2 subdomaincentres
⌊ ⌋
and a grid of s=80 80 evaluation points.
×
10