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

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.

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,2,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, 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 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 $\Omega$ be an open bounded connected nonempty subset of ${\mathbb{R}}^d$ (usually $d=1,2,3$), with Lipschitz-continuous boundary. We also use the classical notation $H^k(\Omega )$ to denote the usual Sobolev space of all distributions u, whose partial derivatives, up to order k, are in the classical Lebesgue space ${\mathcal{L}}^2(\Omega )$.

Let $m>1$ be a positive integer (usually small) and let ${\Pi }_{m-1}\left({\mathbb{R}}^d\right)$ denote the space of polynomials on ${\mathbb{R}}^d$ of total degree at most $m-1$, whose dimension is denoted by $d\left(m\right)$ and let $\left\{q_1,\dots ,q_{d\left(m\right)}\right\}$ be the standard basis of ${\Pi }_{m-1}\left({\mathbb{R}}^d\right)$. Let us give an arbitrary finite set $\left\{{\mathbf{b}}_1,\dots ,{\mathbf{b}}_M\right\}\subset \Omega \subset {\mathbb{R}}^d$ of distinct approximation points and a set of real values $U=\mathbf{f}\equiv \left\{f_1,\dots ,f_M\right\}$, that usually comes from the evaluation of a certain function $f:\Omega \to \mathbb{R}$. We also use a set of centers $\left\{{\mathbf{c}}_1,\dots ,{\mathbf{c}}_N\right\}\subset {\mathbb{R}}^d$ and, for each of them, the translation of an appropriate radial function $\Phi (·-{\mathbf{c}}_i),i=1,\dots ,N$, as in (2).

The main goal is to approximate, in the best possible way, the points

$$\left\{({\mathbf{b}}_1,f_1),\dots ,({\mathbf{b}}_M,f_M)\right\}\subset {\mathbb{R}}^{d+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 $\epsilon \in {\mathbb{R}}^{+},t\ge 0$):

$${\varphi }_{\epsilon }\left(t\right)=-\frac{1}{2{\epsilon }^3}(e^{-\epsilon \sqrt{t}}+\epsilon \sqrt{t})$$

(1)

that is the fundamental solution, in the distribution sense, of the linear fourth-order differential equation

$$\frac{d^4{\varphi }_{\epsilon }}{d{t}^4}-{\epsilon }^2\frac{d^2{\varphi }_{\epsilon }}{d{t}^2}=\delta ,$$

where $\delta$ is the Dirac’s measure at the origin. The real parameter $\epsilon$ is called a tension parameter.

In our case, we use radial basis functions that are the translations of the radial function (now, for $\mathbf{x}\in {\mathbb{R}}^d$)

$${\displaystyle {\Phi }_{\epsilon }\left(\mathbf{x}\right)={\varphi }_{\epsilon }\left({\langle \mathbf{x}\rangle }_d^2\right)=-\frac{1}{2{\epsilon }^3}\left(e^{-\epsilon <\mathbf{x}{>}_d}+\epsilon <\mathbf{x}{>}_d\right)}$$

(2)

where $<·{>}_d$ denote the Euclidean norm on ${\mathbb{R}}^d$.

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 $H^N$, whose basis is of the type

$$H^N\equiv span\left\{{\Phi }_{\epsilon }(·-{\mathbf{c}}_1),\dots ,{\Phi }_{\epsilon }(·-{\mathbf{c}}_N),q_1(·),\dots ,q_{d\left(m\right)}(·)\right\}.$$

In this case, the problem to be resolved can be formulated as follows: Given a set of real values and an approximation data set $\left\{({\mathbf{b}}_i,f_i):i=1,\dots ,M\right\},$ we still want to obtain ${\sigma }_N\in H^N$ such that

$${\sigma }_N\left({\mathbf{b}}_i\right)\approx f_i,i=1,\dots ,M$$

by means of the following minimization problem.

Find ${\sigma }_N\in H^N$ such that

$$\mathcal{J}\left({\sigma }_N\right)\le \mathcal{J}\left(\upsilon \right),\forall \upsilon \in H^N$$

(3)

where $\tau \in (0,\infty )$ and

$${\displaystyle \mathcal{J}\left(v\right)=\sum _{i=1}^M{\left(f_i-\upsilon \left({\mathbf{b}}_i\right)\right)}^2+\tau {\left|\upsilon \right|}_2^2}$$

where

$${\displaystyle {\left|\upsilon \right|}_2^2=\sum _{\left|\alpha \right|=2}{\int }_{\Omega }{\left({\partial }^{\alpha }\upsilon \left(\mathbf{x}\right)\right)}^{2}d\mathbf{x}}$$

with

$$\left|\alpha \right|={\alpha }_1+\dots +{\alpha }_d,\forall \alpha =({\alpha }_1,\dots ,{\alpha }_d)\in {\mathbb{N}}^d,{\partial }^{\alpha }\upsilon (x_1,\dots ,x_d)=\frac{{\partial }^{\left|\alpha \right|}\upsilon }{\partial x_1^{{\alpha }_1}\dots \partial x_d^{{\alpha }_d}}.$$

Thus,

$${\displaystyle {\sigma }_N\left(\mathbf{x}\right)=\sum _{l=1}^{N+d\left(m\right)}{\beta }_l{w}_l\left(\mathbf{x}\right),}$$

where

$$w_l(·)=\left\{\begin{array}{cc}{\Phi }_{\epsilon }(·-{\mathbf{c}}_l),\hfill & l=1,\dots ,N\hfill \end{array}{q}_{l-N}(·),l=N+1,\dots ,N+d\left(m\right)\right.;$$

and ${\beta }_1,\dots ,{\beta }_{N+d\left(m\right)}\in \mathbb{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 ${\sigma }_N\in H^N$ such that, for all $\upsilon \in H^N$,

$${\displaystyle \sum _{i=1}^M{\sigma }_N\left({\mathbf{b}}_i\right)\upsilon \left({\mathbf{b}}_i\right)+\tau {({\sigma }_N,\upsilon )}_2=\sum _{i=1}^M\upsilon \left({\mathbf{b}}_i\right)f_i}$$

(4)

where

$${\displaystyle {(u,v)}_2=\sum _{\left|\alpha \right|=2}{\int }_{\Omega }\left({\partial }^{\alpha }u\left(\mathbf{x}\right){\partial }^{\alpha }\upsilon \left(\mathbf{x}\right)\right)d\mathbf{x}}$$

Proof. 

The expression

$$\left[\left[\upsilon \right]\right]:={\left(\sum _{i=1}^M\upsilon {\left({\mathbf{b}}_i\right)}^2+\tau {(v,v)}_2\right)}^{\frac{1}{2}}$$

constitutes a norm on $H^N$ equivalent to its usual norm. Then, as a consequence, the continuous symmetric bilinear form

$$a:H^N\times H^N⟶\mathbb{R}$$

defined by

$${\displaystyle a(u,\upsilon )=\sum _{i=1}^{M}u\left({\mathbf{b}}_i\right)\upsilon \left({\mathbf{b}}_i\right)+\tau {(u,\upsilon )}_2}$$

is $H^N$-elliptic. Besides, $\psi :H^N⟶\mathbb{R}$, so

$${\displaystyle \psi \left(\upsilon \right):=\sum _{i=1}^M{f}_i\upsilon \left({\mathbf{b}}_i\right)}$$

is clearly a continuous linear form on this space.

Then, applying the Lax–Milgram lemma (see [29] for example), there exists a unique function ${\sigma }_N\in H^N$ such that

$$a({\sigma }_N,\upsilon )=\psi \left(\upsilon \right),\forall \upsilon \in H^N;$$

that is, (4) holds. Moreover, ${\sigma }_N$ is the unique function of $H^N$ that minimizes the functional

$$\tilde{\mathcal{J}}\left(\upsilon \right)=2\left(\frac{1}{2}a(\upsilon ,\upsilon )-\psi \left(\upsilon \right)\right),$$

concluding that ${\sigma }_N$ also minimizes the functional $\mathcal{J}$. □

By linearity, we can reduce the problem (4) to the following linear system:

$$\left(\mathbf{A}{\mathbf{A}}^{\top }+\tau \mathbf{R}\right)\beta =\mathbf{A}\mathbf{f}$$

where ${\displaystyle \mathbf{A}=\left(w_i\left({\mathbf{b}}_j\right)\right),i=1,\dots ,N+d\left(m\right),j=1,\dots ,M.}$

$$\mathbf{R}={\left({\displaystyle \sum _{\left|\alpha \right|=2}{\int }_{\Omega }{\partial }^{\alpha }{w}_i{\partial }^{\alpha }{w}_{j}d\mathbf{x}}\right)}_{\begin{array}{c}\\ i,j=1,\dots ,N+d\left(m\right)\end{array}}$$

$$\beta ={\left({\beta }_i\right)}_{i=1,\dots ,N+d\left(m\right)}\in {\mathbb{R}}^{N+d\left(m\right)},\mathbf{f}\equiv {\left(f_i\right)}_{i=1,\dots ,M}$$

Remark 1.

When $M=N$ and ${\mathbf{b}}_i\equiv {\mathbf{c}}_i$ for $i=1,\dots N$, we could also consider an interpolating problem at the centers, where the interpolating function ${\widehat{\sigma }}_N$ has explicitly the following form:

$${\widehat{\sigma }}_N\left(\mathbf{x}\right)=\sum _{i=1}^N{\lambda }_i{\Phi }_{\epsilon }(\mathbf{x}-{\mathbf{c}}_i)+\sum _{j=1}^{d\left(m\right)}{\beta }_j{q}_j\left(\mathbf{x}\right),\mathbf{x}\in {\mathbb{R}}^d,$$

(5)

where the coefficients $\lambda ={({\lambda }_1,\dots ,{\lambda }_N)}^{\top }$ and $\beta ={({\beta }_1,\dots ,{\beta }_{d\left(m\right)})}^{\top }$ are now the solution of the following linear system:

$$\begin{array}{c}\hfill \left(\begin{array}{cc}{\mathbf{A}}_{\epsilon }& \mathbf{M}\\ {\mathbf{M}}^{\top }& \mathbf{O}\end{array}\right)\left(\begin{array}{c}\lambda \\ \beta \end{array}\right)=\left(\begin{array}{c}\mathbf{f}\\ \mathbf{0}\end{array}\right)\end{array}$$

(6)

where ${\mathbf{A}}_{\epsilon }={\left({\Phi }_{\epsilon }({\mathbf{c}}_i-{\mathbf{c}}_j)\right)}_{1\le i,j\le N}$ is an $N\times N$ matrix, $\mathbf{M}={\left(q_j\left({\mathbf{c}}_i\right)\right)}_{\begin{array}{c}1\le i\le N\\ 1\le j\le d\left(m\right)\end{array}}$ is an $N\times d\left(m\right)$ matrix, ${\mathbf{M}}^{\top }$ denotes the transpose of $\mathbf{M},\mathbf{O}$ is the $d\left(m\right)\times d\left(m\right)$ zero matrix, $\mathbf{f}\equiv {(f_1,\dots ,f_N)}^{\top }\in {\mathbb{R}}^N$ and $\mathbf{0}$ is the zero vector of ${\mathbb{R}}^{d\left(m\right)}$.

Take into account that we are considering the interpolating conditions

$${\widehat{\sigma }}_N\left({\mathbf{c}}_i\right)=f_i,i=1,\dots ,N;$$

together with the constraints (the usual orthogonality conditions to the space of polynomials ${\Pi }_{m-1}\left({\mathbb{R}}^d\right)$ when we do not have the strict positiveness of the radial kernel function, to assure the ${\Pi }_{m-1}\left({\mathbb{R}}^d\right)$-unisolvency of the centers; see for example [1])

$$\sum _{i=1}^N{\lambda }_i{q}_j\left({\mathbf{c}}_i\right)=0,j=1,\dots ,d\left(m\right).$$

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\equiv (g_1,\dots ,g_l,\dots ,g_L)$, 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:

• p and P, the individual and the main population, respectively.
• $S_p$ is the set of all individuals dominated by p.
• $n_p$ is the number of individuals that p dominates.
• $p_{rank}$ is the rank that an individual p has, depending on the front assigned to it.
• Q is 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, $\forall p\in P$:

⋆

Initialize $S_p=\varnothing$ and $n_p=0$

⋆

$\forall q\in P-\left\{p\right\}$

$$\mathrm{if} p \mathrm{dominates} q, \mathrm{then} S_p=S_p\cup \left\{q\right\}$$

$$\mathrm{else}, \mathrm{if} q \mathrm{dominates} p, \mathrm{then} n_p=n_p+1$$

⋆

$\mathrm{if} n_p=0,\mathrm{then} {\mathcal{F}}_1={\mathcal{F}}_1\cup \left\{p\right\} \mathrm{and} p_{ rank }=1$.

⋆

Begin initializing the new front, from the first one (with $i=1$).

⋆

Next, for every front ${\mathcal{F}}_i$, while ${\mathcal{F}}_i\ne \varnothing$,

$$\mathcal{Q}=\varnothing$$

$$\forall p\in {\mathcal{F}}_i \mathrm{and} \forall q\in S_p$$

$$n_q=n_q-1$$

$$if n_q=0, \mathrm{then} \mathcal{Q}=\mathcal{Q}\cup \left\{q\right\} \mathrm{and} q_{ rank }=i+1$$

⋆

$i=i+1$ (augment the front counter by one).

⋆

At this stage, the new front is ${\mathcal{F}}_i=\mathcal{Q}$.

(3)

Crowding distance: It is very important to compare between individuals belonging to the same front ${\mathcal{F}}_i$, 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 ${\mathcal{F}}_i$ and every objective function $g_l$, $l=1,\dots ,L$, we perform the following:

⋆

We previously sort all of the members of the front ${\mathcal{F}}_i$, according to the lth objective function $sort({\mathcal{F}}_i,l)\equiv \{p_{i,1}^{\left[l\right]},\dots ,p_{i,N_i}^{\left[l\right]}\}, with N_i\equiv card\left({\mathcal{F}}_i\right)$; (where $g_l\left(p_{i,1}^{\left[l\right]}\right)\equiv g_l^{ min }$ and $g_l(p_{i,N_i}^{\left[l\right])}\equiv g_l^{max}$). Now, we use this arrangement to calculate the crowding distance for each individual in ${\mathcal{F}}_i$ with the following passages.

⋆

First, initialize the distances of each member of the front to zero.

⋆

Namely, $\forall p_{i,j}^{\left[l\right]}\in {\mathcal{F}}_i,withj=1,\dots ,N_i$. Let us define and compute the following:

⋆

$d_{i,1}^{\left[l\right]}=d_{i,N_i}^{\left[l\right]}=\infty$ (to assure that these points are taken for next generations).

⋆

$d_{i,j}^{\left[l\right]}\equiv d_{p_{i,j}^{\left[l\right]}}^{\left[l\right]}\left({\mathcal{F}}_i\right)=0,j=2,\dots ,N_i-1.$

⋆

For $k=2:(N_i-1)$,

$$d_{i,k}^{\left[l\right]}=d_{i,k}^{\left[l\right]}+\frac{g_l\left(p_{i,k+1}^{\left[l\right]}\right)-g_l\left(p_{i,k-1}^{\left[l\right]}\right)}{{g_l}^{max}-{g_l}^{min}}$$

(4)

Selection: In order to select individuals, we use a crowded comparison operator that depends on a binary tournament selection. The operation depends on ${\mathcal{F}}_i$ and $p_{rank}$; we say that $p{\prec }_{c}q$ if

$$p_{rank}<q_{rank}$$

$$orifpandq\in {\mathcal{F}}_i$$

$$then{d}_p^{\left[l\right]}\left({\mathcal{F}}_i\right)>d_q^{\left[l\right]}\left({\mathcal{F}}_i\right)$$

(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, $p_1$ and $p_2$:

$$c_{1,k}=\frac{1}{2}\left((1-{\beta }_k)p_{1,k}+(1+{\beta }_k)p_{2,k}\right)$$

$$c_{2,k}=\frac{1}{2}\left((1-{\beta }_k)p_{1,k}+(1+{\beta }_k)p_{2,k}\right)$$

with ${\beta }_k\ge 0$ and $c_{i,k}$ represents the kth component to the ith child. However, the density of any random sample is generated by

$$p\left(\beta \right)=\left\{\begin{array}{cc}\frac{1}{2}({\eta }_c+1){\beta }^{{\eta }_c},\hfill & \hfill \mathrm{if}0\le \beta \le 1\\ \frac{1}{2}({\eta }_c+1)\frac{1}{{\beta }^{{\eta }_c+2}},\hfill & \hfill \mathrm{if}\beta >1\end{array}\right.$$

where ${\eta }_c$ is defined as the distribution index of crossover. Let u be a uniform sample random number included in $(0,1)$; then, we have

$$\beta \left(u\right)=\left\{\begin{array}{cc}{\left(2u\right)}^{\frac{1}{{\eta }_c+1}}\hfill & \hfill ,\mathrm{if} u<0.5\\ \frac{1}{{2(1-u)]}^{\frac{1}{{\eta }_c+1}}}\hfill & \hfill ,\mathrm{otherwise}\end{array}\right.$$

−

Polynomial Mutation: First, we define some important notions, such as the following:

○

$c_{i,k}$ is a child component.

○

$p_{i,k}$ is a parent component.

○

$r_k$ is a uniform sample random in $(0,1)$.

○

${\delta }_k$ is a mutation distribution index.

○

$[p_{i,k}^{ lower },p_{i,k}^{ upper}]$ is the box constraints imposed over the kth-component.

Then, the polynomial mutation is defined componentwise by

$$c_{i,k}=p_{i,k}+({p_{i,k}}^{upper}-{p_{i,k}}^{lower}){\delta }_k$$

where

$${\delta }_k=\left\{\begin{array}{cc}{\left(2r_k\right)}^{\frac{1}{{\eta }_l+1}}-1,\hfill & \hfill ,if{r}_k<0\\ 1-{\left(2(1-r_k)\right)}^{\frac{1}{{\eta }_l+1}},\hfill & \hfill ,if{r}_k\ge 0\end{array}\right.$$

(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 decision-maker 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 $E_{\mathcal{C}}$ and/or $E_{\mathcal{L}}$ 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 $\mathcal{C}(\Omega )$ and ${\mathcal{L}}^2(\Omega )$, respectively, and they are given by the following expressions:

$$E_{\mathcal{C}}={\displaystyle \frac{{\displaystyle \underset{i=1,\dots ,K}{max}|f\left({\mathbf{a}}_i\right)-\sigma \left({\mathbf{a}}_i\right)|}}{{\displaystyle \underset{i=1,\dots ,K}{max}\left|f\left({\mathbf{a}}_i\right)\right|}}}$$

(7)

$${\displaystyle E_{\mathcal{L}}=\sqrt{\frac{{\displaystyle \sum _{i=1}^K{(f\left({\mathbf{a}}_i\right)-\sigma \left({\mathbf{a}}_i\right))}^2}}{{\displaystyle \sum _{i=1}^K{\left(f\left({\mathbf{a}}_i\right)\right)}^2}}}}$$

(8)

where $f\in {\mathcal{C}}^2(\Omega )$ is a given function, $\sigma \equiv {\sigma }_N$ is the approximating or $\sigma \equiv \widehat{{\sigma }_N}$ interpolating RBF function associated with the given data set and $\left\{a_1,\dots ,a_K\right\}\subset \Omega$ 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 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\times {10}^{-2}$
Fitness’ functions	Centers number and $E_{\mathcal{C}}$ or/and $E_{\mathcal{L}}$

below.

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 $\epsilon =0.1$ was used without problems in all our computations. The same applies to the other smoothness parameter of the RBFs approximating procedure, $\tau$, where a constant fixed value of 1 was used in all our numerical experiments.

One-Dimensional Examples

Example 1.

$f_1:[0,\pi ]⟶\mathbb{R}$

$$\begin{array}{cc}& f_1\left(x\right)=0.12+0.25{exp}^{-4{(x-\frac{\pi }{4})}^2}cos\left(2x\right)sin\left(2\pi x\right)\hfill \end{array}.$$

To test the behavior of the approximation for the presented methodology, performed by optimization of the centers’ placement of radial basis functions by a NSGA-II algorithm, Figure 1 shows the results of the approximation of function $f_1\left(x\right)$, using radial basis function approximations, while, in Figure 2, we can see how the $E_{\mathcal{L}}$ error decreased clearly when the number of centers increased. We can also see, in the left column of Figure 1, the radial basis approximation corresponding to the last centers’ distribution, whose evolution is also showed in the right column, with the increase in the number of interior centers for the function $f_1\left(x\right)$ with 7, 9 and 11 centers, respectively.

Example 2.

$f_2:[0,1]⟶\mathbb{R}$

$$\begin{array}{cc}& f_2\left(x\right)=0.5e^{-5{(x-\pi /2)}^2}sin\left(4\pi x\right)cos\left(4x\right)\hfill \end{array}.$$

For Example 2, the numerical approximation results for $f_2\left(x\right)$ and the evaluation of centers are shown in Figure 3; on the other hand, the two normalized discretised errors for $f_2\left(x\right)$ are also described in

Table 2. Two approximation discretization errors for $f_2\left(x\right)$.

Centers	${\mathit{E}}_{\mathcal{L}}$	${\mathit{E}}_{\mathcal{C}}$
2	$4.3978\times {10}^{-2}$	$1.8918\times {10}^{-3}$
4	$2.6664\times {10}^{-2}$	$6.7993\times {10}^{-4}$
6	$4.0691\times {10}^{-3}$	$1.9208\times {10}^{-4}$
8	$3.9539\times {10}^{-3}$	$1.7241\times {10}^{-4}$

.

We can see, in the left column of Figure 3, the radial basis approximation corresponding to the last centers’ distribution, whose evolution is also showed at the right column, with the increase in the number of interior centers for the function $f_2\left(x\right)$, using 2, 4 and 6 interior centers, respectively.

Example 3.

$f_3:[-1,1]⟶\mathbb{R}$

$$\begin{array}{cc}& f_3\left(x\right)=1-e^{-50\left|x\right|}\hfill \end{array}.$$

The graphical result approximating $f_3\left(x\right)$ with 20 interior centers, using 20 generations in the NSGA-II algorithm, using radial basis functions, is showed in Figure 4.

5. Conclusions and Future Work

Although there are some other approaches to select knots amongst a uniform partition of the interval (or domain) of definition of the function to be approximated, or interpolated, they are usually based on splines and the evolutionary procedure used to improve the final allocation of these knots usually relies on basic genetic algorithms, not very performing for general multi-objective problems, or even on self-adjusting or hybrid neural networks, with the requirement of the previous hard training process, etc. This paper presents another important and interesting possibility for optimizing both the number and the corresponding centers’ placement in the problem of radial basis function approximation of functions of one or several independent variables, where we can clearly see the effectiveness of the strategy for different types of functions. 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. In addition, the preliminary numerical results of this MOGA strategy, for centers’ placement in approximating functions of one variable, show that the placement of the centers in radial basis approximation function also has an important and significant impact on the behavior of the final results, because the ideal place or location of the centers is not known in advance. However, one important feature of this approach is that completely the same theoretical procedure can be developed for problems in more spatial dimensions, for the well-known meshless character of these type of functions. It is also well known that increasing the number of centers for the radial basis in the RBF approximation usually also increases the accuracy of the approximation, but only up to a certain level, where the computational effort is not worthy enough; therefore, in our procedure, we can admit the coalescence of very close centers, maintaining the same level of accuracy. In addition, in the interpolating case, this possible, excessive proximity of the centers has to be avoided, so as not to have excessive matrix conditioning problems with the linear systems to be solved. On the other hand, in classic allocation methods, the designer should choose a priori the number of centers that will be used and, many times, they are just equidistant uniform partitions of the domain. However, using such a MOGA procedure of the type considered, they can move freely inside the domain and the corresponding Pareto fronts can also be directly improved using a different or variable number of centers, at the same time as simulations are being performed.

An open problem would be the optimal election of both the shape parameter $\epsilon$ and the smoothness parameter $\tau$, since they can depend not only on the number and distribution of data points and the data vector, but even on the precision of the computations, even though it is well known that very small values give, in general, highly accurate results. Then, the idea would be to introduce, in a future work, these parameters also as possible optimization variables in the general MOGA algorithm. In addition, the possibility of using the basis of RBFs of compact support, such as the well-known Wendland’s RBFs, would be a good alternative to the RBFs used in this work.