Optimal Centers’ Allocation in Smoothing or Interpolating with Radial Basis Functions

mathematics Article Optimal Centers’ Allocation in Smoothing or Interpolating with Radial Basis Functions Pedro González-Rodelas 1,*,† , Hasan M. H. Idais 2,†, Mohammed Yasin 3,† and Miguel Pasadas 1,†   Citation: González-Rodelas, P.; Idais, H.M.H.; Yasin, M.; Pasadas, M. Optimal Centers’ Allocation in Smoothing or Interpolating with Radial Basis Functions. Mathematics 2022,10, 59. https://doi.org/ 10.3390/math10010059 Academic Editors: Sara Remogna, Domingo Barrera and María José Ibáñez Received: 20 November 2021 Accepted: 15 December 2021 Published: 24 December 2021 Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Copyright: © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). 1Departamento Matemática Aplicada, Universidad de Granada, 18071 Granada, Spain; [email protected] 2Department of Mathematics and Statistics, Arab American University, Jenin P240, Palestine; [email protected] 3Department of Mathematics, An-Najah National University, Nablus P400, Palestine; [email protected] *Correspondence: [email protected] † These authors contributed equally to this work. Abstract: Function interpolation and approximation are classical problems of vital importance in many science/engineering areas and communities. In this paper, we propose a powerful methodology for the optimal placement of centers, when approximating or interpolating a curve or surface to a data set, using a base of functions of radial type. In fact, we chose a radial basis function under tension (RBFT), depending on a positive parameter, that also provides a convenient way to control the behavior of the corresponding interpolation or approximation method. We, therefore, propose a new technique, based on multi-objective genetic algorithms, to optimize both the number of centers of the base of radial functions and their optimal placement. To achieve this goal, we use a methodology based on an appropriate modification of a non-dominated genetic classification algorithm (of type NSGA-II). In our approach, the additional goal of maintaining the number of centers as small as possible was also taken into consideration. The good behavior and efficiency of the algorithm presented were tested using different experimental results, at least for functions of one independent variable. Keywords: approximation; interpolation; RBFs; centers’ allocation; MOGA; NSGA-II algorithms 1. Introduction In the last decades, radial basis function (RBF) methods have been developed in a way that contributed to their great spread and emergence, mainly for their simplicity and for having been lauded for the ability to solve multivariate scattered data approximation and interpolation problems (see for example [ 1 – 3 ], and many other references included therein). The main advantages of these methods are that they do not need a mesh or triangulation, that they are simple to implement and dimension independent, and that no staircasing or polygonization for boundaries is required [4,5]. Many authors have presented several methodologies for approximating or fitting a surface to data (see for example [ 6 , 7 ]), such as the methodology that is based on the polynomial interpolation theorem of Lagrange [ 8 , 9 ] and another important approach, that is based on orthogonal polynomials [10]. Also, some seminal works on dynamic node adaptive strategy, or knot selection placement for sparse reduced data, were published in [11,12]. In [ 13 ], the methodologies for single objective problems (SOP) have been presented for selection and optimization of centers, mainly based on techniques for the selection of center points through non-dominated sorting, where promising points are selected as centers. Also, some successful hierarchical genetic algorihtm approaches for curve fitting with B-splines were applied in [14]. On the other hand, the authors in [ 15 ] developed techniques to perform centers’ selection and estimation simultaneously by a componentwise boosting algorithm. Moreover, Mathematics 2022,10, 59. https://doi.org/10.3390/math10010059 https://www.mdpi.com/journal/mathematics Mathematics 2022,10, 59 2 of 14 the centers of RBFs are selected among the vertices of the Voronoi diagram of the sample data point in [ 16 ]. The methodology in [ 17 , 18 ] was to optimize the shape parameter of radial basis functions for the solution of ordinary differential equations. The authors in [ 19 ] used an approach based on genetic algorithms to optimize centers and radial basis function neural networks (RBFNNs). However, they used one type of neural network where the basis function is of Gaussian type. Research in neural networks has also attracted many computer and mathematical scientists over the last years and covers a broad spectrum of different subjects and real-world applications, such as image processing, pattern recognition and optimization in general (the reader can also consult the new upcoming book of the authors Grienggrai Rajchakit, Praveen Agarwal and Sriraman Ramalingam [ 20 ] for more information about the multiple applications and stability properties of these type of neural networks and many other mathematical aspects of them). In this paper, a new methodology for optimal placement of random centers, for approximating or fitting a curve to data, using radial basis functions is developed. In fact, we chose a radial basis function under tension (RBFT) that depends on a positive parameter, that also provides a convenient way to control the behavior of the corresponding interpolation or approximation method. A new technique is presented to optimize both the number of radial basis centers and their optimal placement, using multi-objective genetic algorithms. More precisely, we used a new approach based on an appropriate modification of a goal-based non-dominated genetic classification algorithm (of type NSGA-II) and tested the good behavior and efficiency of the algorithm using different experimental results, at least for functions of one independent variable. However, the procedure and theory can be developed in the same way for functions of any number of variables, so we present it in this more general setting. The key points concerning this goal-based NSGA-II version, compared with the original NSGA, are the following (see for example [21,22]): • A fast, non-dominated sorting mechanism, grouping the entire population into a hierarchy of sub-populations, based on the Pareto dominance ordering. • Inside each sub-group, the similarity between members is evaluated on the Pareto front; to promote a diverse front of non-dominated solutions, the resulting groups and similarity measures are used. • For that, an adequate crowding distance is an appropriate mechanism of ranking members of a front which are dominating, or dominated by others. • The appropriate consideration of a certain elitist mechanism, in order to improve the convergence properties of the algorithm. • A niching operator without parameters is adopted, to maintain a certain level of diversity among the possible solutions. • One or several decision maker’s (DM) preferences within the evolutionary process are provided, in order to focus into a particular region of the Pareto front, where the more interesting features are present. • In this way, a so-called region of interest (ROI) appears; it is the set of non-dominated solutions that are preferred by the DM. • The DM can provide his/her preferences before (a priori), after (a posteriori), or during (interactively) the MOEA run. • Therefore, some trade-off between objectives specifies that the gain of a unit on one target is worth the degradation on some others and vice versa. • In our case, keeping the number of centers as small as possible while controlling some specific standard errors is our main goal. Then, these ranking mechanisms are used with the corresponding genetic selection operators (usually of Tournament type) to create the population of possible solutions of the next generations. Regarding implementation issues, we developed a MATLAB code through a suitable adaptation of some previous subroutines for the NSGA-II algorithm by Mostapha Kalami Mathematics 2022,10, 59 3 of 14 Heris in [ 23 ] (see also [ 24 ]); we have also implemented the corresponding approximation or interpolation procedures, using the specified RBFs according with [25]. 2. Proposed Methodology Let Ω be an open bounded connected nonempty subset of Rd (usually d= 1, 2, 3), with Lipschitz-continuous boundary. We also use the classical notation Hk(Ω) to denote the usual Sobolev space of all distributions u , whose partial derivatives, up to order k , are in the classical Lebesgue space L2(Ω). Let m> 1 be a positive integer (usually small) and let Πm−1(Rd) denote the space of polynomials on Rd of total degree at most m− 1, whose dimension is denoted by d(m) and let nq1, . . . , qd(m)o be the standard basis of Πm−1(Rd) . Let us give an arbitrary finite set {b1, . . . , bM}⊂Ω⊂Rd of distinct approximation points and a set of real values U=f≡{f1, . . . , fM} , that usually comes from the evaluation of a certain function f:Ω→R . We also use a set of centers {c1, . . . , cN}⊂Rd and, for each of them, the translation of an appropriate radial function Φ(·−ci),i=1, . . . , N, as in (2). The main goal is to approximate, in the best possible way, the points {(b1,f1), . . . , (bM,fM)}⊂Rd+1 by means of an adequate linear combination of radial basis functions (RBFs). 2.1. Spaces of Radial Functions We define the radial basis functions that is used in the sequel. We consider the following function (for ε∈R+,t≥0): φε(t) = −1 2ε3(e−ε√t+ε√t)(1) that is the fundamental solution, in the distribution sense, of the linear fourth-order differential equation d4φε dt4−ε2d2φε dt2=δ, where δ is the Dirac’s measure at the origin. The real parameter ε is called a tension parameter. In our case, we use radial basis functions that are the translations of the radial function (now, for x∈Rd) Φε(x) = φε(hxi2 d) = −1 2ε3e−ε<x>d+ε<x>d(2) where <·>ddenote the Euclidean norm on Rd. This type of radial basis functions under tension (RBFTs) were introduced by A. Bouhamidi and A. Le Mehauté [ 26 ] in 2004. They depend on a positive parameter that provides a convenient way of controlling the behavior of the corresponding interpolation or approximation method. Moreover, error estimates of the corresponding interpolation method [ 26 ] are given in terms of Sobolev semi-norms, that also help in the convergence of the associated procedure and provide attractive geometrical tension effects. 2.2. Smoothing Radial Basis Approximations Therefore, the idea is to use a finite-dimensional space HN , whose basis is of the type HN≡spannΦε(·−c1), . . . , Φε(·−cN),q1(·), . . . , qd(m)(·)o. Mathematics 2022,10, 59 4 of 14 In this case, the problem to be resolved can be formulated as follows: Given a set of real values and an approximation data set {(bi,fi):i=1, . . . , M} , we still want to obtain σN∈HNsuch that σN(bi)≈fi,i=1, . . . , M by means of the following minimization problem. Find σN∈HNsuch that J(σN)≤ J(υ),∀υ∈HN(3) where τ∈(0, ∞)and J(v) = M ∑ i=1 (fi−υ(bi))2+τ|υ|2 2 where |υ|2 2=∑ |α|=2ZΩ(∂αυ(x) )2dx with |α|=α1+. . . +αd,∀α= (α1, . . . , αd)∈Nd,∂αυ(x1, . . . , xd) = ∂|α|υ ∂xα1 1. . . ∂xαd d . Thus, σN(x) = N+d(m) ∑ l=1 βlwl(x), where wl(·) = Φε(·−cl),l=1, . . . , N ql−N(·),l=N+1, . . . , N+d(m); and β1 , . . . , βN+d(m)∈R are the control coefficients (i.e., the problem’s unknowns), obtained as the solution of the variational problem formulated in the following Theorem (see for example [27,28]). Theorem 1. The minimization problem (3) has a unique solution, that is also the only solution of the following variational problem. Find σN∈HNsuch that, for all υ∈HN, M ∑ i=1 σN(bi)υ(bi) + τ(σN,υ)2= M ∑ i=1 υ(bi)fi(4) where (u,v)2=∑ |α|=2ZΩ(∂αu(x)∂αυ(x) )dx Proof. The expression [[υ]] := M ∑ i=1 υ(bi)2+τ(v,v)2!1 2 constitutes a norm on HN equivalent to its usual norm. Then, as a consequence, the continuous symmetric bilinear form a:HN×HN−→ R Mathematics 2022,10, 59 5 of 14 defined by a(u,υ) = M ∑ i=1 u(bi)υ(bi) + τ(u,υ)2 is HN-elliptic. Besides, ψ:HN−→ R, so ψ(υ):= M ∑ i=1 fiυ(bi) is clearly a continuous linear form on this space. Then, applying the Lax–Milgram lemma (see [29] for example), there exists a unique function σN∈HNsuch that a(σN,υ) = ψ(υ),∀υ∈HN; that is, (4) holds. Moreover, σN is the unique function of HN that minimizes the functional e J(υ) = 21 2a(υ,υ)−ψ(υ), concluding that σNalso minimizes the functional J. By linearity, we can reduce the problem (4) to the following linear system: AA>+τRβ=Af where A= (wi(bj)),i=1, . . . , N+d(m),j=1, . . . , M. R= ∑ |α|=2ZΩ∂αwi∂αwjdx i,j=1,...,N+d(m) β=(βi)i=1,...,N+d(m)∈RN+d(m),f≡(fi)i=1,...,M Remark 1. When M=N and bi≡ci for i= 1, . . . N , we could also consider an interpolating problem at the centers, where the interpolating function ˆ σNhas explicitly the following form: ˆ σN(x) = N ∑ i=1 λiΦε(x−ci) + d(m) ∑ j=1 βjqj(x),x∈Rd, (5) where the coefficients λ= (λ1 , . . . , λN)> and β= (β1 , . . . , βd(m))> are now the solution of the following linear system: AεM M>Oλ β=f 0(6) where Aε= (Φε(ci−cj))1≤i,j≤N is an N×N matrix, M= (qj(ci)) 1≤i≤N 1≤j≤d(m) is an N×d(m) matrix, M> denotes the transpose of M , O is the d(m)×d(m) zero matrix, f≡(f1, . . . , fN)>∈RN and 0is the zero vector of Rd(m). Take into account that we are considering the interpolating conditions ˆ σN(ci) = fi,i=1, . . . , N; together with the constraints (the usual orthogonality conditions to the space of polynomials Πm−1(Rd) when we do not have the strict positiveness of the radial kernel function, to assure the Πm−1(Rd)-unisolvency of the centers; see for example [1]) Mathematics 2022,10, 59 6 of 14 N ∑ i=1 λiqj(ci) = 0, j=1, . . . , d(m). 3. Methodology Overview At this stage, it is important to define some ideas about multi-objective optimization problems and describe fundamental issues for NSGA-II-type algorithms (see [ 22 ]). The non-dominated sorting genetic algorithm is a multiple objective optimization algorithm, whose objective is to improve the adaptive fit of a population of candidate solutions to a Pareto front, constrained by a set of objective functions G≡(g1 , . . . , gl , . . . , gL) , each one defined in an appropriate domain. Amongst the first algorithms to explicitly exert certain tendency towards the discovery of nondominated solutions, we could mention, for example, Fonseca and Fleming’s MOGA [ 30 ]; it uses fitness sharing amongst solutions of the same rank, coupled along a proportional adaptive selection to help promoting diversity. The authors Srinivas and Deb’s [ 31 ] original nondominated sorting genetic algorithm (NSGA) works in a similar way, but assigns fitness based on dividing the population into a certain number of fronts of equal domination. To achieve this, the algorithm iteratively seeks all the nondominated points in the population that still have not been labeled as belonging to the current front and increments the front count, repeating until all solutions have been finally labeled. Therefore, each point in a given front obtains, as its raw score, the count of all solutions in inferior fronts. After this, Deb [ 22 ] and his coworkers, proposed a revised NSGA-II algorithm, which still uses the idea of non-dominated fronts, but also incorporates the following important changes: • A crowding distance metric is defined for each point as the average side length of the cuboid defined by its nearest neighbors in the same front. The larger the value, the fewer solutions reside in the vicinity of the point. • An appropriate survivor selection strategy, where the new population is obtained by accepting the individuals from progressively inferior fronts until it is full, so that, if not all the individuals in the last front considered can be accepted, they are chosen on the basis of their crowding distance. • Parent selection uses a modified tournament operator that considers, first, the dominance rank, then crowding distance as the second criterium. We can also see that, for the appropriate use of any NSGA-II procedure, it is necessary to describe the following fundamental issues: (1) Population initialization: the size of the population depends on the problem range. (2) Non-dominated sort. Before talking about this stage we define the following: •pand P, the individual and the main population, respectively. •Spis the set of all individuals dominated by p. •npis the number of individuals that pdominates. •prank is the rank that an individual phas, depending on the front assigned to it. •Qis the set for storing temporarily the individuals in the (i+1)-th front. (We recall that, in this context of multi-objective optimization, we say that p dominates q , or, equivalently, that q is dominated by p , if and only if p is no worse than q in all the objective functions, but it is strictly better in at least one of them.) We describe the non-dominated procedure as follows: ?For all individuals, ∀p∈P: ?Initialize Sp=∅and np=0 ?∀q∈P−{p}if pdominates q, then Sp=Sp∪{q} else, if qdominates p, then np=np+1 ?if np=0 , then F1=F1∪{p}and prank =1. Mathematics 2022,10, 59 7 of 14 ?Begin initializing the new front, from the first one (with i=1). ?Next, for every front Fi, while Fi6=∅, Q=∅ ∀p∈ Fiand ∀q∈Sp nq=nq−1 if nq=0, then Q=Q∪{q}and qrank =i+1 ?i=i+1 (augment the front counter by one). ?At this stage, the new front is Fi=Q. (3) Crowding distance: It is very important to compare between individuals belonging to the same front Fi , so that each individual has a crowding distance that estimates the density of solutions surrounding each member in this front. The method of calculating the crowding distance is summarized below. For all considered fronts Fi and every objective function gl , l= 1, . . . , L , we perform the following: ? We previously sort all of the members of the front Fi , according to the lth objective function sort(Fi , l)≡ {p[l] i,1 , . . . , p[l] i,Ni} , with Ni≡card(Fi) ; (where gl(p[l] i,1)≡gmin l and gl(p[l]) i,Ni≡gmax l ). Now, we use this arrangement to calculate the crowding distance for each individual in Fiwith the following passages. ?First, initialize the distances of each member of the front to zero. ? Namely, ∀p[l] i,j∈ Fi , with j= 1, . . . , Ni . Let us define and compute the following: ?d[l] i,1 =d[l] i,Ni=∞(to assure that these points are taken for next generations). ?d[l] i,j≡d[l] p[l] i,j (Fi) = 0, j=2, . . . , Ni−1. ?For k=2 : (Ni−1), d[l] i,k=d[l] i,k+gl(p[l] i,k+1)−gl(p[l] i,k−1) glmax −glmin (4) Selection: In order to select individuals, we use a crowded comparison operator that depends on a binary tournament selection. The operation depends on Fi and prank ; we say that p≺cqif prank <qrank or if p and q ∈ Fi then d[l] p(Fi)>d[l] q(Fi) (5) Genetic operator: Now, we describe the two parts of the simulated binary crossover and that of polynomial mutation, both performed componentwise: − Binary crossover: The following expressions represent the simulated binary crossover (for each component k ) of two contestants in the tournament, p1 and p2: c1,k=1 2((1−βk)p1,k+ (1+βk)p2,k) c2,k=1 2((1−βk)p1,k+ (1+βk)p2,k) Mathematics 2022,10, 59 8 of 14 with βk≥ 0 and ci,k represents the kth component to the ith child. However, the density of any random sample is generated by p(β) = (1 2(ηc+1)βηc, if 0 ≤β≤1 1 2(ηc+1)1 βηc+2, if β>1 where ηc is defined as the distribution index of crossover. Let u be a uniform sample random number included in (0, 1); then, we have β(u) =    (2u)1 ηc+1, if u<0.5 1 2(1−u)] 1 ηc+1 , otherwise − Polynomial Mutation: First, we define some important notions, such as the following: ◦ci,kis a child component. ◦pi,kis a parent component. ◦rkis a uniform sample random in (0, 1). ◦δkis a mutation distribution index. ◦[plower i,k,pupper i,k]is the box constraints imposed over the kth-component. Then, the polynomial mutation is defined componentwise by ci,k=pi,k+ (pi,kupper −pi,klower)δk where δk=   (2rk)1 ηl+1−1, , if rk<0 1−(2(1−rk))1 ηl+1, , if rk≥0 (6) Recombination and selection: The offspring population is combined with the current generation population and selection is performed to set the individuals of the next generation; the elitism is ensured, because the best individuals are added sequentially to the population. The selection of parents for the next generation is also based on their crowding distance, by selecting the individual at random, but taking into account their tournament score, then choosing the best individuals out of that set to be parents. (7) Selection of an appropriate region of interest (ROI), according with DM’s preferences: In this case, the preference of the user is always to maintain the number of centers as small as possible, if the degree of approximation remains within a preset range. In this manner, the proposed procedure attempts to find a set of solutions in the neighborhood of the corresponding Pareto optimal solution, so that the decisionmaker can have a better idea of the desired region. Different forms of fitness (or objective functions) can be used in a NSGA-II procedure, but the main goal, here, is to minimize some usual errors, such as EC and/or EL between the original function and the approximating or interpolating function, constructed from each population of random centers. At this step, it is also important to consider two error estimations that are appropriate normalizations of a discrete version of the usual norms in C(Ω)and L2(Ω), respectively, and they are given by the following expressions: EC= max i=1,...,K|f(ai)−σ(ai)| max i=1,...,K|f(ai)|(7) Mathematics 2022,10, 59 9 of 14 EL=v u u u u u u u t K ∑ i=1 (f(ai)−σ(ai))2 K ∑ i=1 (f(ai))2 (8) where f∈ C2(Ω) is a given function, σ≡σN is the approximating or σ≡ˆ σN interpolating RBF function associated with the given data set and {a1, . . . , aK}⊂Ω is another scattered random point test data set (TDS) where the errors are computed. However, before running any MOGA algorithm, there are always some key parameters that should be defined. The most important ones are shown in Table 1below. Table 1. Parameters and functions used by MOGA in the simulations results. Parameters of the MOGA Value Number of generations 20 Population size 40 Crossover function Binary crossover Selection function Binary tournament Crossover fraction 0.9 Pareto fraction 0.4 Mutation function Polynomial mutation Mutation rate 0.01 Centers’ deletion tolerance 0.5 ×10−2 Fitness’ functions Centers number and ECor/and EL Therefore, the most important goal of these methods will always be to try to improve the adaptive fit of a population of possible solutions to a Pareto front constrained by these often conflicting objectives. In this case, these objective functions are to minimize at least one or several normalized discrete versions of some approximate errors; however, the intention to maintain the number of centers as small as possible is also taken into account, so that we allow the procedure to delete some of the obtained centers, when they become sufficiently close to each other and the same level of approximation can be retained. 4. Simulation Examples The objective of this study is to analyze the performance of this procedure for different types of functions, optimizing the centers placement using radial basis functions by the multi-objective genetic algorithms (MOGA) already used and explained in more detail in [ 32 ], but using cubic spline basis functions instead of RBFs. Different experiments were carried out and we show, in this section, the results for the approximation of several functions, also presenting the evolution of the optimal distribution of points, together with the related Pareto fronts. Similar results could be obtained for the case of interpolating RBFs, as explained in the Remark 1. In order to analyze the behavior of the radial basis approximations, a TDS with a large (100) number of points was used, with a population number of 40 individuals, until 20 generations, in each one of the examples below. Concerning the shape parameter, it is also true that it can affect the approximation accuracy of the RBF and the overall computational stability of the procedure (see for example [ 33 , 34 ]), but it is not the goal of our study for the moment; therefore, a fixed value of ε= 0.1 was used without problems in all our computations. The same applies to the other smoothness parameter of the RBFs approximating procedure, τ , where a constant fixed value of 1 was used in all our numerical experiments.