{'fraction_non_alphanumeric': 0.03495651443216671, 'fraction_numerical': 0.008613083483887231, 'mean_word_length': 4.616328456819255, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Student's t Distribution based Estimation of Distribution Algorithms for Derivative-free Global Optimization

In this paper, we are concerned with a branch of evolutionary algorithms termed estimation of distribution (EDA), which has been successfully used to tackle derivative-free global optimization problems. For existent EDA algorithms, it is a common practice to use a Gaussian distribution or a mixture of Gaussian components to represent the statistical property of available promising solutions found so far. Observing that the Student's t distribution has heavier and longer tails than the Gaussian, which may be beneficial for exploring the solution space, we propose a novel EDA algorithm termed ESTDA, in which the Student's t distribution, rather than Gaussian, is employed. To address hard multimodal and deceptive problems, we extend ESTDA further by substituting a single Student's t distribution with a mixture of Student's t distributions. The resulting algorithm is named as estimation of mixture of Student's t distribution algorithm (EMSTDA). Both ESTDA and EMSTDA are evaluated through extensive and in-depth numerical experiments using over a dozen of benchmark objective functions. Empirical results demonstrate that the proposed algorithms provide remarkably better performance than their Gaussian counterparts.


I. INTRODUCTION
The estimation of distribution algorithm (EDA) is an evolutionary computation (EC) paradigm for tackling derivative-free global optimization problems [1]. EDAs have gained great success and attracted increasing attentions during the last decade [2]. Different from most of EC algorithms which are normally equipped with meta-heuristics inspired selection and variation operations, the EDA paradigm is characterized by its unique variation operation which uses probabilistic models to lead the search towards promising areas of the solution space. The basic assumption is that probabilistic models can be used for learning useful information of the search space from a set of solutions that have already been inspected from the problem structure and this information can be used to conduct more effective search. The EDA paradigm provides a mechanism to take advantage of the correlation structure to drive the search more efficiently. Henceforth, EDAs are able to provide insights about the structure of the search space. For most existent EC algorithms, their notions of variation and selection usually arise from the perspective of Darwinian evolution. These notions offer a direct challenge for quantitative analysis by the existing methods of computer science.
In contrast with such meta-heuristics based EC algorithms, EDAs are more theoretically appealing, since a wealth of tools and theories arising from statistics or machine learning can be used for EDA analysis. Convergence results on a class of EDAs have been given in [3].
All existent EDAs fall within the paradigm of model based or learning guided optimization [2]. In this regard, different EDAs can be distinguished by the class of probabilistic models used. An investigation on the boundaries of effectiveness of EDAs shows that the limits of EDAs are mainly imposed by the probabilistic model they rely on. The more complex the model, the greater ability it offers to capture possible interactions among the variables of the problem. However, increasing the complexity of the model usually leads to more computational cost. The EDA paradigm could admit any type of probabilistic model, among which the most popular model class is the Gaussian model. The simplest EDA algorithm just employs univariate Gaussian models, which regard all design variables to be independent with each other [4,5]. The simplicity of models makes such algorithms easy to implement, while their effectiveness halts when the design variables have strong interdependencies. To get around of this limitation, several multivariate Gaussian based EDAs (Gaussian-EDAs) have been proposed [4,6,7]. To represent variable linkages elaborately, Bayesian networks (BNs) are usually adopted in the framework of multivariate Gaussian-EDA, while the learning of the BN structure and parameters can be very time consuming [8,9]. To handle hard multimodal and deceptive problems in a more appropriate way, Gaussian mixture model (GMM) based EDAs (GMM-EDAs) have also been proposed [10,11]. It is worthwhile to mention that the more complex the model used, the more likely it encounters model overfitting, which may mislead the search procedure [12].
The Student's t distribution is symmetric and bell-shaped, like the Gaussian distribution, but has heavier and longer tails, meaning that it is more prone to producing values that fall far from its mean. Its heavier tail property renders the Student's t distribution a desirable choice to design Bayesian simulation techniques such as the adaptive importance sampling (AIS) algorithms [13][14][15]. Despite the purpose of AIS is totally different as compared with EDAs, AIS also falls within the iterative learning and sampling paradigm like EDAs. The Student's t distribution has also been widely used for robust statistical modeling [16]. Despite of great success achieved in fields of Bayesian simulation and robust modeling, the Student's t distribution has not yet gained any attention from the EC community.
In this paper, we are concerned with whether the heavier tail property of the Student's t distribution is beneficial for developing more efficient EDAs. To this end, we derive two Student's t distribution based algorithms within the paradigm of EDA. They are entitled as estimation of Student's t distribution algorithm (ESTDA) and estimation of mixture of Student's t distribution algorithm (EMSTDA), implying the usages of the Student's t distribution and the mixture of Student's t distributions, respectively.
We evaluate both algorithms using over a dozen of benchmark objective functions. The empirical results demonstrate that they indeed have better performance than their Gaussian counterparts. In literature, the closest to our work here are the posterior exploration based Sequential Monte Carlo (PE-SMC) algorithm [13] and the annealed adaptive mixture importance sampling algorithm [17], while both were developed within the framework of Sequential Monte Carlo (SMC), other than EDA. In addition, the faster Evolutionary programming (FEP) algorithm proposed in [18] is similar in the spirit of replacing the Gaussian distribution with an alternative that has heavier tails. The major difference between FEP and our work here lies in that the former falls within the framework of evolutionary strategy while the latter is developed within the EDA paradigm. In addition, the former employs the Cauchy distribution, other than the student's t or mixture of student's t distributions here, to generate new individuals. To the best of our knowledge, we make the first attempt here to bring Student's t probabilistic models into the field of population based EC.
The remainder of this paper is organized as follows. In Section II we revisit the Gaussian-EDAs. Besides the single Gaussian distribution model based EDA, the GMM-EDA is also presented. In Section III we describe the proposed Student's t model based EDAs in detail. Like in Section II, both the single distribution and the mixture distribution model based EDAs are presented. In Section IV, we present the performance evaluation results and finally in Section V, we conclude the paper.

II. GAUSSIAN DISTRIBUTION BASED EDA
In this section, we present a brief overview on the EDA algorithmic paradigm, the Gaussian-EDA and the GMM-EDA.
In this paper, we are concerned with the following derivative-free continuous global optimization problem where χ denotes the nonempty solution space defined in R n , and f : χ → R is a continuous real-valued function. The basic assumption here is that f is bounded on χ, which means ∃f l > −∞, fu < ∞ such that f l ≤ f (x) ≤ fu, ∀x ∈ χ. We denote the minimal function value as f * , i.e., there exists an

A. Generic EDA scheme
A generic EDA scheme to solve the above optimization problem is presented in Algorithm 1. In this paper, we only focus on the probabilistic modeling part in the EDA scheme, not intend to discuss the selection operations. In literature, the tournament or truncation selection are commonly used with the EDA algorithm. Here we adopt truncation selection for all algorithms under consideration. Specifically, we select, from {xi} N i=1 , M individuals which produce the minimum objective function values. The stopping criteria of Algorithm 1 is specified by the user based on the available budget, corresponding to the acceptable maximum number of iterations here.

Algorithm 1: A main scheme of EDA
Input: the probabilistic model p(x|θ) with parameter θ, the population size N , the selection size M , which satisfies M < N , and θ 1 , the initialized value of θ. Output:x and f (x), which represent estimates of x * and f (x * ), respectively. 1 Initialization: Set θ = θ 1 and then draw a random samplex from p(·|θ = θ 1 ). Set the iteration index k = 1; 2 while the stopping criterion is not met do 3 Sample N individuals from p(x|θ k ), denote them by x 1 , x 2 , . . . , x N ; if y min < f (x) then 7 Letx = x min and f (x) = y min ; Update θ k−1 , based on the selected individuals, to get θ k . See Subsections II-B, II-C, III-B and III-C for details on specific operations with respect to the Gaussian model, GMM, the Student's t and the mixture of Student's t models, respectively. Set k = k + 1.
As mentioned before, we employ a truncation selection strategy, which indicates that the selected individuals are those which produce minimum objective function values. Given the selected individuals {xi, yi} M j=1 , a maximum likelihood (ML) estimate of the parameters of the Gaussian distribution is shown to be: where T denotes transposition and all vectors involved are assumed to be column vectors both here and hereafter.

C. GMM-EDA
As its name indicates, the GMM-EDA uses a GMM to represent p(x|θ), namely , l is the index of the mixing components, w the probability mass of the mixing component which satisfies L l=1 w l = 1, and L denotes the total number of mixing components included in the mixture and is preset as a constant. Now we focus on how to update θ k−1 to get θ k based on the selected individuals {xj } M j=1 . We use an expectation-maximization (EM) procedure, as shown in Algorithm 2, to update the parameters of the mixture model. For more details on EM based parameter estimation for GMM, readers can refer to [19][20][21][22]. The stopping criterion can be determined by checking the changing rate of the GMM parameters. If θnew has not been changed a lot for a fixed number of successive iterations, then we terminate the iterative process of the EM procedure. Compared with a benchmark EM procedure, here we add a components deletion operation at the end of the Expectation step. We delete components with extremely small mixing weights, because their roles in the mixture are negligible, while at the Maximization step, they will occupy the same computing burden as the rest of components. This deletion operation is also beneficial for avoiding disturbing numerical issues such as singularity of covariance matrix. For the influence of this operation on the the number of survival mixing components, see Section IV-B.

III. STUDENT'S T DISTRIBUTION BASED EDA
In this section, we describe the proposed ESTDA and EMSTDA in detail. Both of them are derived based on the application of the Student's t distribution in the EDA framework as presented in Algorithm 1. To begin with, we give a brief introduction to the Student's t distribution.

A. Student's t distribution
Suppose that x is a d dimensional random variable that follows the multivariate Student's t distribution, denoted by S(·|µ, Σ, v), where µ denotes the mean, Σ a positive definite inner product matrix and v ∈ (0, ∞] is the degrees of freedom (DoF). Then the density function of x is: where denotes the Mahalanobis squared distance from x to µ with respect to Σ, A −1 denotes the inverse of A and Γ(·) denotes the gamma function. For ease of understanding, we give a graphical description of univariate Student's t pdfs corresponding to different DoFs in comparison with a standard Gaussian pdf in Fig.1. It is shown that the Student's t pdfs have heavier tails than the Gaussian pdf. The smaller the DoF v, the heavier the corresponding pdf's tails. From Fig.1, we see that the Student's t pdf is more likely to generate an individual further away from its mean than the Gaussian model due to its long flat tails. It leads to a higher probability of escaping from a local optimum or moving away from a plateau. We can also see that heavier tails lead to a smaller hill around the mean, and vice versa. A smaller hill around the mean indicates that the corresponding model tuning Algorithm 2: EM procedure to estimate parameters of the GMM within the EDA paradigm as presented in Algorithm 1 Input: the individuals {x j } M j=1 , the model's current parameter value θ k−1 = {w k−1,l , µ k−1,l , Σ k−1,l } L l=1 and the admittable smallest weight W of a mixing component (we set W = 0.02 as a default value for use in our experiments) Output: θ k 1 Initialization: Set θ old = θ k−1 . Set the iteration index i = 1; 2 while the stopping criterion is not met do 3 The Expectation step: where ǫ(l|x j ) denotes the probability of the event individual x j belonging to the lth mixing component, which is calculated as follows Delete mixing components whose mixing weights are smaller than W . Then update the value of L and increase the mixing weights of the remaining mixing components proportionally to guarantee that their summation equals 1;

4
The Maximization step: Set Set θ old = θ new and let k = k + 1.
operation spends less time in exploiting the local neighborhood and thus has a weaker fine-tuning ability. Hence heavier tails do not always bring advantages. This has been demonstrated by our numerical tests. See Fig.2d, which shows that for several two-dimensional (2D) benchmark test problems, lighter tails associated with v = 50 provide better convergence result than heavier tails corresponding to v = 5 and 10.
First, let's figure out, given the parameter value θ, how to sample an individual x from p(·|θ). It consists of two steps that is simulating a random draw τ from the gamma distribution, whose shape and scale parameters are set identically to v/2, and sampling x from the Gaussian distribution N (·|µ, Σ/τ ).
Now let's focus on, given a set of selected individuals {xj } M j=1 , how to update the parameters of the student's t distribution in a optimal manner in terms of maximizing likelihood. As mentioned above, in generation of {xj } M j=1 , a corresponding set of gamma distributed variables {τj } M j=1 is used. We can record these gamma variables and then use them to easily derive an ML estimate for parameters of the Student's t distribution [23,24]: the mixture of Student's t model can also be found in [13][14][15], where the mixture learning procedure is performed based on a weighted sample set. The EM operations presented in Algorithm 3 can be regarded as an application of the EM methods described in [13][14][15] to an equally weighted sample set. Similarly as in Algorithm 2, we add a components deletion operation at the end of Expectation step. This deletion operation were used in the adaptive annealed importance sampling [14] algorithm and the PE-SMC algorithm [13], which also adopt the EM procedure to do parameter estimation for a Student's t mixture model. By deleting components with extremely small mixing weights, we can avoid consuming computing burdens to update parameter values for negligible components. As presented before in Subsection II-C, this deletion operation is also useful to avoid common numerical issues such as singularity of covariance matrix. For the influence of this operation on the the number of survival mixing components, see Subsection IV-B.

IV. PERFORMANCE EVALUATION
In this section, we evaluate the presented ESTDA and EMSTDA using a number of benchmark objective functions, which are designed and commonly used in literature for testing and comparing optimization algorithms. See the Appendix Section for details on the involved functions. The Gaussian-EDA and the GMM-EDA, which are described in Sections II-B and II-C, respectively, are involved for comparison purpose.
Although scaling up EDAs to large scale problems has become one of the biggest challenges of the field [25], it is not our Algorithm 3: EM procedure to estimate parameters of a mixture of student's t model within the EDA paradigm as presented in Algorithm 1 Input: the individuals {x j } M j=1 , the model's current parameter value θ k−1 = {w k−1,l , µ k−1,l , Σ k−1,l } L l=1 and the admittable smallest weight W of a mixing component (we set W = 0.02 as a default value for use in our experiments) Output: θ k 1 Initialization: Set θ old = θ k−1 . Set the iteration index i = 1; 2 while the stopping criterion is not met do 3 The Expectation step: where ǫ(l|x j ) denotes the probability of the event individual x j belonging to the lth mixing component, which is calculated as follows ǫ(l|x j ) = w old,l S(x j |µ old,l , Σ old,l , v) L l=1 w old,l S(x j |µ old,l , Σ old,l , v) .
Delete mixing components whose mixing weights are smaller than W . Then update the value of L and increase the mixing weights of the remaining mixing components proportionally to guarantee that their summation equals 1;

4
The Maximization step:

B. A more thorough performance evaluation for ESTDA and EMSTDA
In this subsection, we perform a more thorough performance evaluation for ESTDA and EMSTDA. Both Gassuian-EDA and GMM-EDA are also involved for comparison purpose. Including GMM-EDA we hope to further demonstrate that the heavier tails of the Student's t distribution are beneficial for exploration of solution space via the EDA mechanism. By contrasting EMSTDA with ESTDA, we can inspect how much a mixture model can contribute to search the solution space via a Student's t based EDA mechanism.
To make a representative evaluation, we will consider in total 17 objective functions listed in Table I. For all but the Rastrigin and Michalewicz functions, we consider 2D cases. For the Rastrigin and Michalewicz functions, we also consider 5D and 10D cases. The population size per iteration N is set to be 1e3, 1e4, and 1e5 for 2D, 5D, 10D cases, respectively. The selection size M is fixed to be 0.2 × N for all cases involved. We didn't consider higher dimensional cases, for which a reliable probabilistic model building requires far more samples, while it is impractical to carry out so large numbers of samples in reasonable time scales. Definitions of all the involved objective functions are presented in the Appendix of the paper.
We fix the value of the DoF v to be 5 for all test functions, except that, for problems involving the Rastrigin function, v is set to be 50, since it gives better convergence behavior as reported in subsection IV-A. For every algorithm, the iterative EDA procedure terminates when the iteration number k in Algorithm 1 is bigger than 50. For EMSTDA and GMM-EDA, the iterative EM procedure terminates when the iteration number i in Algorithms 2 and 3 is bigger than 2. Each algorithm is run 30 times independently. Then we calculate the average and standard error of the converged fitness values. The results are given in Table I. For each test function, the best solution is marked with bold font. Then we count for each algorithm how many times it output I: Convergence results yielded from 30 independent runs of each algorithm on each test problem. a ± b in the table denotes that the average and standard error of the best fitness values obtained from 30 independent runs of the corresponding algorithm are a and b, respectively. The best solution for each problem is marked with bold font. a better solution as compared with all the rest of algorithms and put the result into Table II. We see that the EMSTDA get the highest score 7, followed by score 5 obtained by ESTDA, while their Gaussian counterparts only obtain relatively lower scores 0 and 2, respectively. This result further coincides with our argument that the Student's t distribution is more preferable than the Gaussian for use in designing EDAs. We also see that, between ESTDA and EMSTDA, the latter performs better than the former.
This demonstrates that using mixture models within a Student's t based EDA is likely to be able to bring additional advantages.
What slightly surprises us in Table II is that the GMM-EDA seems to perform worse than the Gaussian-EDA, while after a careful inspection on Table I To investigate the effect of the component deletion operator used in Algorithms 3 and 2, we record the number of survival mixing components per iteration of the ESTDA and EMSTDA, and calculate its mean averaged over 30 independent runs of each algorithm. The result for six typical objective functions is depicted in Fig. 3. We see that, for both algorithms, the mean number of survival components decreases as the algorithm iterates on. This phenomenon suggests that the mixture model mainly takes effect during the early exploration stage, while, when the algorithm goes into the latter exploitation phase, there usually remains much less components in the mixture model. Note that in some cases, e.g, the Cross-in-tray function case, more than one algorithms provide an identical best solution. For such cases, no algorithm will get a score in the statistics made here.

V. CONCLUDING REMARKS
In this paper we introduced the Student's t distribution, for the first time, into a generic paradigm of EDAs. Specifically, we developed two Student's t based EDAs, which we term as ESTDA and EMSTDA here. The justification for the proposed algorithms lies in the assumption that the heavier tails of the Student's t distribution are beneficial for an EDA type algorithm to explore the solution space. Experimental results demonstrate that for most cases under consideration, the Student's t based EDAs indeed performed remarkably better than their Gaussian counterparts implemented in the same algorithm framework. As a byproduct of our experimental study, we also demonstrated that mixture models can bring in additional advantages, for both Student's t and Gaussian based EDAs, in tackling global optimization problems.
All algorithms developed here are based on density models that use full covariance matrices. The reason is two-fold. One is for the ease of presentation of the basic idea, namely, replacing the commonly used Gaussian models by heavier tailed Student's

A. Ackley Function
This function is defined to be where a = 20, b = 0.2 and c = 2π. This function is characterized by a nearly flat outer region, and a large peak at the centre, thus it poses a risk for optimization algorithms to be trapped in one of its many local minima. The global minima is f (x ⋆ ) = 0, at x ⋆ = (0, . . . , 0).

D. Rastrigin Function
This function is evaluated on the hypercube xi ∈ [−5.12, 5.12], for all i = 1, . . . , d, with the form Its global minimum f (x ⋆ ) = 0 is located at x ⋆ = (0, . . . , 0). This function is highly multimodal, but locations of the local minima are regularly distributed.

Q. Rosenbrock Function
This function is also referred to as the Valley or Banana function. The global minimum lies in a narrow, parabolic spike.
However, even though this spike is easy to find, convergence to the minimum is difficult [26]. It has the following form and is evaluated on the hypercube xi ∈ [−5, 10], for all i = 1, . . . , d.