Overview of Estimation of Distribution Algorithms (EDAs)

Abstract

1  Introduction

Estimation of Distribution Algorithms (EDAs) [ Mühlenbein and Paaß1996, Larrañaga and Lozano2001] constitute a new and promising stochastic heuristic for optimization that can be considered as member of the family of evolutionary computation heuristics. Due to their recent development, their application to classic and real optimization problems is still starting. This work presents a first approximation to a type of combinatorial optimization problems -searching for the best permutation- that have a big relevance in areas such as Operational Research and Machine Learning.

2  Estimation of Distribution Algorithms

2.1  Introduction

EDAs [ Mühlenbein and Paaß1996, Larrañaga and Lozano2001] are non-deterministic, stochastic heuristic search strategies that form part of the evolutionary computation approaches, where number of solutions or individuals are created every generation, evolving once and again until a satisfactory solution is achieved. In brief, the characteristic that most differentiates EDAs from other evolutionary search strategies such as GAs is that the evolution from a generation to the next one is done by estimating the probability distribution of the fittest individuals, and afterwards by sampling the induced model. This avoids the use of crossing or mutation operators, and the number of parameters that EDAs require is reduced considerably.

In EDAs, individuals are not said to contain genes, but variables which dependencies have to be analyzed. Also, while in other heuristics from evolutionary computation the interrelations between the different variables representing the individuals are kept in mind implicitly (e.g. building block hypothesis), in EDAs the interrelations are expressed explicitly through the joint probability distribution associated with the individuals selected at each iteration. The task of estimating the joint probability distribution associated with the database of the selected individuals from the previous generation constitutes the hardest work to perform, as this requires the adaptation of methods to learn models from data developed in the domain of probabilistic graphical models.

Figure 1 shows the pseudocode of EDA, in which we distinguish four main steps in this approach:

1. At the beginning, the first population D₀ of M individuals is generated, usually by assuming an uniform distribution (either discrete or continuous) on each variable, and evaluating each of the individuals.
2. Secondly, a number N (N £ M) of individuals are selected, usually the fittest.
3. Thirdly, the n-dimensional probabilistic model that better expresses the interdependencies between the n variables is induced.
4. Next, the new population of M new individuals is obtained by simulating the probability distribution learnt in the previous step.

Next, some notation is introduced. Let X = (X₁, ¼, X_n) be a set of random variables, and let x_i be a value of X_i, the i^th component of X. Let y = (x_i)_{Xi Î Y} be a value of Y Í X. Then, a probabilistic graphical model for X is a graphical factorization of the joint generalized probability density function, r(X = x) (or simply r(x)). The representation of this model is given by two components: a structure and a set of local generalized probability densities.

Regarding the structure of the model, the structure S for X is a directed acyclic graph (DAG) that describes a set of conditional independences between the variables on X. Pa_i^S represents the set of parents -variables from which an arrow is coming out in S- of the variable X_i in the probabilistic graphical model, the structure of which is given by S. The structure S for X assumes that X_i and its non descendents are independent given Pa_i^S, i = 2, ¼, n. Therefore, the factorization can be written as follows:

r(x) = r(x1, ¼, xn) = n Õ i=1  r(xi | paiS).

(1)

Furthermore, regarding the local generalized probability densities associated with the probabilistic graphical model, these are precisely the ones appearing in Equation 1.

A representation of the models of the characteristics described above assumes that the local generalized probability densities depend on a finite set of parameters q_S Î Q_S, and as a result the previous equation can be rewritten as follows:

r(x | qS) = n Õ i=1  r(xi | paiS, qi)

(2)

where q_S = (q₁, ¼, q_n).

After having defined both components of the probabilistic graphical model, the model itself will be represented by MO = (S,q_S).

2.2  EDAs in discrete domains

In the particular case where every variable X_i Î X is discrete, the probabilistic graphical model is called Bayesian network. If the variable X_i has r_i possible values, x_i¹, ¼, x_i^r_i, the local distribution, p (x_i | pa_i^j,S,q_i) is:

p (xik | paij,S, qi) = qxik | paij º qijk

(3)

where pa_i^{1, S}, ¼, pa_i^{q_i, S} denotes the values of Pa_i^S, that is the set of parents of the variable X_i in the structure S; q_i is the number of different possible instantiations of the parent variables of X_i. Thus, q_i = Õ_{Xg Î Pai} r_g. The local parameters are given by q_i = ((q_ijk)_k=1^r_i)_j=1^q_i). In other words, the parameter q_ijk represents the conditional probability that variable X_i takes its k^th value, knowing that the set of its parent variables take their j^th value. We assume that every q_ijk is greater than zero.

All the EDAs are classified depending on the maximum number of dependencies between variables that they accept (maximum number of parents that any variable can have in the probabilistic graphical model).

2.2.1  Without interdependencies

The Univariate Marginal Distribution Algorithm (UMDA) [ Mühlenbein1998] is a representative example of this category, which can be written as:

pl(x;ql) = n Õ i=1  pl(xi;qil)

(4)

where q_i^l = { q_ijk^l} is recalculated every generation by its maximum likelihood estimation, i.e. [^(q)]_ijk^l = [(N_ijk^l-1)/(N_ij^l-1)]. N_ijk^l is the number of cases on which the variable X_i takes the value x_i^k when its parents are on their j^th combination of values for the l^th generation, and N_ij^l-1=å_k N_ijk^l-1.

2.2.2  Pairwise dependencies

An example of this second category is the greedy algorithm called MIMIC (Mutual Information Maximization for Input Clustering) [ De Bonet et al. 1997]. The main idea in MIMIC is to describe the true mass joint probability as closely as possible by using only one univariate marginal probability and n - 1 pairwise conditional probability functions.

2.2.3  Multiple interdependencies

We will use EBNA (Estimation of Bayesian Network Algorithm) as an example of this category. The EBNA approach was introduced for the first time in [ Etxeberria and Larrañaga1999], where the authors use the Bayesian Information Criterion (BIC) as the score to evaluate the goodness of each structure found during the search.

The initial model MO₀ in EBNA, is formed by its structure S₀ which is an arc-less DAG and the local probability distributions given by the n unidimensional marginal probabilities p(X_i = x_i)=[1/(r_i)], i = 1, ¼, n -that is, MO₀ assigns the same probability to all individuals. The model of the first generation -MO₁- is learnt using Algorithm B [ Buntine91], while the rest of the models are learnt following a local search strategy that received the model of the previous generation as the initial structure.

2.2.4  Simulation in Bayesian networks

In EDAs, the simulation of Bayesian networks is used merely as a tool to generate new individuals for the next population based on the structure learned previously. The method used in this work is the Probabilistic Logic Sampling (PLS) [ Henrion88]. Following this method, the instantiations are done one variable at a time in a forward way, that is, a variable is not sampled until all its parents have already been so.

2.3  EDAs in continuous domains

In this section we introduce an example of the probabilistic graphical model paradigm that assumes the joint density function to be a multivariate Gaussian density.

The local density function for the i^th variable is computed as the linear-regression model

f (xi | paiS, qi) º N (xi; mi + å xj Î pai  bji (xj - mj), vi)

(5)

where N (x; m, s²) is a univariate normal distribution with mean m and variance s².

Local parameters are given by q_i = (m_i,b_i, v_i), where b_i = (b_1i, ¼, b_i-1i)^t is a column vector. Local parameters are as follows: m_i is the unconditional mean of X_i, v_i is the conditional variance of X_i given Pa_i, and b_ji is a linear coefficient that measures the strength of the relationship between X_j and X_i. A probabilistic graphical model built from these local density functions is known as a Gaussian network. Gaussian networks are of interest in continuous EDAs because the number of parameters needed to specify a multivariate Gaussian density is smaller.

Next, an analogous classification of continuous EDAs as for the discrete domain is done, in which these continuous EDAs are also classified depending on the number of dependencies they take into account.

2.3.1  Without dependencies

In this case, the joint density function is assumed to follow a n-dimensional normal distribution, and thus it is factorized as a product of n unidimensional and independent normal densities. Using the mathematical notation X º N (x;m, å), this assumption can be expressed as:

fN ( x; m, å ) = n Õ i=1  fN (xi; mi, si) = n Õ i=1  1 Ö   2p   si e-[1/2]([(xi - mi)/(si)])2. (6)

An example of continuous EDAs in this category is UMDA_c [ Larrañaga et al. 2000].

2.3.2  Bivariate dependencies

An example of this category is MIMIC_c^G [ Larrañaga et al. 2000], which is basically an adaptation of the MIMIC algorithm [ De Bonet et al. 1997] to the continuous domain.

2.3.3  Multiple dependencies

Algorithms in this section are approaches of EDAs for continuous domains in which there is no restriction in the learning of the density function every generation. An example of this category is EGNA_BGe (Estimation of Gaussian Network Algorithm) [ Larrañaga et al. 2000]. The method used to find the Gaussian network structure is a Bayesian score+search. In EGNA_BGe a local search is used to search for good structures.

2.3.4  Simulation of Gaussian networks

A general approach for sampling from multivariate normal distributions is known as the conditioning method, which generates instances of X by sampling X₁, then X₂ conditionaly to X₁, and so on. The simulation of a univariate normal distribution can be done with a simple method based on the sum of 12 uniform variables.

3  Example of the Griewangk problem

This is a minimization problem proposed in [Törn and Zilinskas 1989]. The fitness function is as follows:

F(x) = 1 + n å i=1  xi2 4000 - n Õ i=1  cos æ ç è xi Öi ö ÷ ø .

(2)

The range of all the components of the individual is -600\leqslant x_i \leqslant 600 , i = 1, ... , n, and the fittest individual corresponds to a value of 0, that only can be obtained when all the components of the individual are 0.

4.3  The evolution of the search

As an example to illustrate the difference in the way of reaching the final result for all the continuous EDAs, Figure 3, 4, 5, 6, and 7, show which is the behaviour of the UMDAc, MIMICc, Edge Exclusion, EMNA-global and EMNA-adaptative algorithms for the Griewangtk problem with a dimension of 2.

5  Conclusions

At the light of the results obtained in the fitness values, we can conclude the following: generally speaking, for less dimensions EMNA_global, and EGNA type algorithms perform better that the rest, but when increasing the dimensions some of the algorithms show a poorer performance as a result of the higher complexity to take into account (e.g. the case of EGNA_BIC). The EMNA_a algorithm showed a very poor behaviour for all these optimization problems, and its additional computation effort made impossible to apply it to the 50 dimension cases.

An important aspect to take into account is that the EMNA_global algorithm appears to be the method that more quickly approaches to the best results, although these results are not always the optima. Nevertheless, once this algorithm is nearby the optimum solution it requires more time than algorithms as EGNA_BGe or EGNA_ee to satisfy the ending criterion.

Depending on the problem the ES method showed better results than the continuous EDAs, but when the type of problem to optimize presents many local minima and maxima, continuous EDAs show a more appropriated behaviour. The main drawback for continuous EDAs in general is the computation time they require, but for some problems the results that can be obtained with them are not comparable to methods in the ES category.

References