The most important and crucial step of EDAs is how to build the probabilitic model P (x) to express the promising solutions. In EDAs for global continuous optimization problem, the Gaussian distribution is a common one. Some other complex models, like Gaussian mixture, histogram etc., are also used [33J. Sun, Q. Zhang, and E.P. Tsang, "DE/EDA: A new evolutionary algorithm for global optimization", Inf. Sci., vol. 169, no. 3-4, pp. 249-262, 2005.
[http://dx.doi.org/10.1016/j.ins.2004.06.009] ]. It must be noted that, different dependency relations can appear between variables for large-scale optimization problems. Considering all the dependencies needs a complex probability model, which is hard to realize. Therefore, we specifically focus on the use of univariate probability model to reduce the amount of calculation. Even though the dependency between variables is very important for system modeling, the optimum solution is the most important issue for optimization problem. Hence, every variable is assumed independent of any variable in the probability model. The Gaussian model is used to model and estimate the distribution of promising solutions in every dimension of the problem.

In order to construct a Gaussian pdf model of the promising solutions, we should obtain the statistical information of promising solutions. Hence, statistical techniques have been extensively applied to the optimization problems. Fortunately, these parameters can be efficiently computed by the maximum likelihood estimations [31P.A. Bosman, and D Thierens, "Numerical optimization with real-valued estimation of distribution algorithms", In: Scalable Optimization via Probabilistic Modeling, Springer: Berlin, Heidelberg, 2006, pp. 91-120.
[http://dx.doi.org/10.1007/978-3-540-34954-9_5] ]. In the algorithm assuming full independence, every variable is assumed independent of any variable. That is, the probability distribution P(x₁, x₂, … , x_D) of the vector (x₁, x₂, … , x_D) of m variables is assumed to consist of a product of the distributions of individual variables:

(1)

where µ^k_i is the mean and σ^k_i is the standard deviation of k-th generation and i-th variable. D is the dimension size. This is very suitable for calculation. Different from the discrete EDAs, the number of parameters to be estimated does not grow exponentially with D.

The pdf for variables x_i is parameterized by the mean σ^k_i and the standard deviation , which is defined by

(2)

Therefore, the probability distribution P(x₁, x₂, … , x_D) of the vector (x₁, x₂, … , x_D) of m variables is

(3)

The parameters (µ^k_i, σ^k_i) can be estimated according to the selected best individuals. The parameters (μ_i, σ_i) can be updated every iteration.

The mean and standard deviation parameters of promising population can be computed adaptively by the maximum likelihood technique according to the selected promising solutions.

(4)

(5)

µ_i (k) is the mean of i-th variable in k-th iteration, NB is the selected individuals size. σ²_i (k) is the covariance of i-th variable in k-th iteration.

Fig. (2) Cartogram of sampling data.

The probability sampling is used to generate new individuals using the learned probabilistic models instead of crossover or mutation operators. The sampling method depends on the type of probabilistic model and the characteristics of the problem. For normal pdf problem, a conversion can be used in order to convert the normal pdf to a standard normal pdf.

Supposing,

(6)

The normal pdf about x is converted to a standard normal pdf about y.

(7)

The variable x can be calculated by

(8)

In the probability models, every variable (x₁, x₂, … , x_m) is assumed independent of any variable. The mean and standard deviation of variable x_i is μ_i and σ_i, when n→∞,

(9)

If x_i is the evenly distributed random number of [0, 1]

(10)

(11)

Therefore,

(12)

when n→∞, y→N(0, 1). We can select an appropriate n to generate a normal pdf for probability sampling. Fig. (2) shows the cartogram of sampling data in different n. From the figure we can see the sampling data follow the pdf better with the increasing of n.

The parameters of Gaussian model are always learnable in the process of optimization. Some iterative learning approaches are used in some literatures [34M. Nakao, T. Hiroyasu, M. Miki, H. Yokouchi, and M. Yoshimi, "Real-coded estimation of distribution algorithm by using probabilistic models with multiple learning rates", Proc. Comput. Sci., vol. 4, pp. 1244-1251, 2011.
[http://dx.doi.org/10.1016/j.procs.2011.04.134] -37M. Gallagher, M. Frean, and T Downs, "Real-valued evolutionary optimization using a flexible probability density estimator", In: Proceedings of the Genetic and Evolutionary Computation Conference, Orlando: USA, 1999, pp. 840-846.], and we can conclude as follows.

(13)

(14)

where α and β are a fixed weights of µ_i (k) and µ_i (k-1). The learning method depending on the class of models used, this step involves parametric or structural learning, also known as model fitting and model selection, respectively. This can improve the performance of EDAs, no matter how simple or complex the learning rule is. In order to update the probability model in this paper, we adopt a different learning rule for u and σ.u is the candidate of optimum solution. Therefore, we hope it updated in real-time. Hence, we adopt a fast learning rule for u (α=1, β=0). The standard deviation σ should be a large value to search widely at beginning [34M. Nakao, T. Hiroyasu, M. Miki, H. Yokouchi, and M. Yoshimi, "Real-coded estimation of distribution algorithm by using probabilistic models with multiple learning rates", Proc. Comput. Sci., vol. 4, pp. 1244-1251, 2011.
[http://dx.doi.org/10.1016/j.procs.2011.04.134] ]. Thus, the σ is set a squared size of 1/2 of domain. In order to update the P(x), the learning rule of σ is shown in equation (15).

(15)

w is the learning rate. A bigger w can provide better real time ability of σ. A smaller w can ensure better robust of σ. In the paper, the value of w is adaptive to the iterations.

Fig. (3) Population operation diagram.

(16)

where w_max and w_min are the maximal and minimal value of w. FEs is the current iterations, and Max_FEs is the maximal iterations.

Elitism strategy is an effective strategy to ensure the best individual(s) is selected as the next generation in EAs, because the best individual(s) maybe include the information of optimal solution. Therefore, elitism can improves the convergence performance of EAs in many cases [38C.W. Ahn, and R.S. Ramakrishna, "Elitism-based compact genetic algorithms", IEEE Trans. Evol. Comput., vol. 7, no. 4, pp. 367-385, 2003.
[http://dx.doi.org/10.1109/TEVC.2003.814633] ], and elitism has long been considered an effective method for improving the efficiency of EAs [39R.C. Purshouse, and P.J. Fleming, "Why use elitism and sharing in a multi-objective genetic algorithm", In: Proceedings of the Genetic and Evolutionary Computation Conference, New York: USA, 2002, pp. 520-527.]. This is achieved by simply copying the best individual(s) directly to the new generation [40I.J. Leno, S.S. Sankar, M.V. Raj, and S.G. Ponnambalam, "An elitist strategy genetic algorithm for integrated layout design", Int. J. Adv. Manuf. Technol., vol. 66, no. 9-12, pp. 1573-1589, 2013.]. However, the number of best individuals selected as the next generation must be handled properly and carefully otherwise may lead to premature convergence or can not improve the efficiency of algorithm. Fig. (3) is the process of new population generation. The elitism individuals will be selected as the new generation directly, and the best individuals are used to establish a probability model to generate other individuals of next generation.

(17)

where Elitism() is the operator to copy the best solution to Pop(k+1), and Sample() is the sampling function. N is the population size, NB is the number of best individuals selected to build probability model.

It is widely accepted that a local search procedure is efficient in improving the solutions generated by the EDA. Kinds of strategies are proposed to enhance the performance of EDA, including the combing with PSO [4C.W. Ahn, J. An, and J. Yoo, "Estimation of particle swarm distribution algorithms: Combining the benefits of PSO and EDAs", Inf. Sci., vol. 192, no. 1, pp. 109-119, 2012.
[http://dx.doi.org/10.1016/j.ins.2010.07.014] ] and DE [18Y. Wang, B. Li, and T. Weise, "Estimation of distribution and differential evolution cooperation for large scale economic load dispatch optimization of power systems", Inf. Sci., vol. 180, no. 12, pp. 2405-2420, 2010.
[http://dx.doi.org/10.1016/j.ins.2010.02.015] ], quantum EDA [13M.E. Platel, S. Schliebs, and N. Kasabov, "Quantum-inspired evolutionary algorithm: a multimodel EDA", IEEE Trans. Evol. Comput., vol. 13, no. 6, pp. 1218-1232, 2009.
[http://dx.doi.org/10.1109/TEVC.2008.2003010] ], etc. In this paper, we use a chaos perturbation as the local search strategy. The principle of perturbation is shown as Fig. (4).

Fig. (4) The principle of perturbation.

The running of chaos operator is conditional. In this paper, the chaos perturbation is running under the condition of slower convergence.

(18)

where is the i-th variable of j-th individual of k-th iteration. η is a perturbation coefficient. z_i is the chaotic variable, which can be generated by chaotic models. There are many chaotic models to generate the chaotic variables [41M.S. Tavazoei, and M. Haeri, "Comparison of different one-dimensional maps as chaotic search pattern in chaos optimization algorithms", Appl. Math. Comput., vol. 187, no. 2, pp. 1076-1085, 2007.
[http://dx.doi.org/10.1016/j.amc.2006.09.087] ], such as Logistic mapping, Cube mapping or infinite folding mapping. Logistic chaotic model is the most popular one, which folded within a limited number under a limited range [42Q. Xu, S. Wang, L. Zhang, and Y. Liang, "A novel chaos danger model immune algorithm", Commun. Nonlinear Sci. Numer. Simul., vol. 18, no. 11, pp. 3046-3060, 2013.
[http://dx.doi.org/10.1016/j.cnsns.2013.04.017] ]. The logistic model is shown as follows.

(19)

It is a typical chaotic system. is the control variable, and a definite time series can be generated by iteration for any .

The model of infinite folding mapping is shown as equation (20):

Fig. (5) Distribution property of logistic, cube and infinite folding mapping.

(20)

The cube mapping is shown as equation (21)

(21)

The Cube and infinite folding mapping do not need control variable u. The distribution of chaotic sequence in logistic mapping is asymmetric, which follows the property of Chebyshev. The three mappings have 10000 iterations, and the statistics the results of the distribution properties of the mappings are shown in Fig. (5). The distribution property of infinite folding mapping is better than other one. Therefore, the infinite folding mapping is adopted to generate the chaotic variables. Additionally, the infinite folding mapping is insensitivity to the initial value [42Q. Xu, S. Wang, L. Zhang, and Y. Liang, "A novel chaos danger model immune algorithm", Commun. Nonlinear Sci. Numer. Simul., vol. 18, no. 11, pp. 3046-3060, 2013.
[http://dx.doi.org/10.1016/j.cnsns.2013.04.017] ].

The above way to generate new solutions does not take into account the feasibility of the solutions. The individuals could be out of the domain due to the perturbation. Therefore, a repair procedure is needed if illegal individuals are constructed.

(22)

[lb_iub_i] is the domain of i-th variable.

With the design above, the procedure of the ALREEDA is illustrated as following.

Initialization: Set parameters:, Max_FEs, NP, NB, and generate population Pop(0).

While (stop criteria ?)

Generate w according FEs.

Evaluation: Calculate the fitness of all individuals, and store the elitism.

Statistical information obtaining: Select NB individuals to estimate the parameter of the probabilistic model, and update the parameter.

Probabilistic model building: According to estimated parameters, build the probabilistic model of each variable x_i

Probabilistic sampling:

Make use of the sampling technology sampling (NP-NB) individuals.

Elitism strategy: Combine the sampling individuals with elitism and generate new Pop(k).

Chaos perturbation (perturbation criteria ?):

In the iteration optimization algorithm, such as EAs, the maximal iteration number (Max_FEs) is an essential parameter [19W. Dong, T. Chen, P. Tino, and X. Yao, "Scaling up estimation of distribution algorithms for continuous optimization", IEEE Trans. Evol. Comput., vol. 17, no. 6, pp. 792-822, 2013.
[http://dx.doi.org/10.1109/TEVC.2013.2247404] ], while it is maybe varied. The Max_FEs is set to 3E+6 in this paper. For EDAs, the population size NP and the promising solution number NB are important except for Max_FEs. It is obvious that for an easy problem, a small value of NP is sufficient, but for difficult problems, a large value of NP is recommended in order to avoid trapping to a local optimum. The large NP may provide better optimization with larger calculation [16H. Karshenas, R. Santana, C. Bielza, P. Larra, and N. Aga, "Regularized continuous estimation of distribution algorithms", Appl. Soft Comput., vol. 13, no. 5, pp. 2412-2432, 2013.
[http://dx.doi.org/10.1016/j.asoc.2012.11.049] ]. However, it may vary from problem to problem. We pay attention to the performance of EDA instead of the population size. We provide some choices (100, 200, 300, 500, 1000) to obtain better performance for the ALREEDA. We have a comparisons of different population size, and select the best population size as the final decision of the algorithm on the problem with the given problem size. The NB is also an important parameter for the probabilistic model learning of EDA, and we also have a comparison to determine a proper NB. In the comparisons, the error recorded finally is the absolute margin between the fitness of the best solution found and the fitness of the global optimum. The semi-log graphs show log10 (error) vs. FES for the first 6 functions.

In Fig. (10), it is the testing result when the dimension D is 100 and has different NP and NB. We can see from the figure, the optimization results are similar for f₁, f₃, f₅, f₆ and f₇ when the population size is 100, 200, 300 or 500. According to f₂ and f₄, it is a proper choice that the population size is 500. Further more, the algorithm have a better performance when NP is 500 and NB is 50.

We also do some tests when the dimension D is 500. Firstly, we test the algorithm under different population size 100, 200 or 500. From Fig. (11), we can see the performance of the algorithm on f₁, f₃, f₅ and f₆ are similar regardless of population size. It only effect the optimization time. According to the performance of the algorithm on f₂, f₄ and f₇, the population size is selected as 500. 50 is a suitable value for NB according to the comprehensive performance of the algorithm on the benchmarks.

Fig. (10) Convergence graph for function 1-7 (D=100, different NP).

Due to the time reason, we didn’t do the population size test for D=1000. We testify the algorithm directly on NP=500 and NB=50. Fig. (12) shows the performance of the algorithm, and also have a comparison of different dimension size when NP=500 and NB=50.

According to Fig. (12), we also can conclude that f1, f5 and f6 are easy to solve for different dimension sizes. However, the performance is promising for f2 when D=100, the results are worse when D is 500 or 1000 with D=100. For f3 and f4, the error is increasing with the dimension.

From the above mentioned, the population size is selected as 500, and NB is 50 for dimension 100, 500 and 1000. We also can conclude that the ALREEDA is scale to the dimension size except for the f2.

We also do a small test to testify the effect of the elitism strategy. The optimization process partial enlarge graph of f7, f4, and f2 are shown in Fig. (13). f7 is the most obvious one. The optimization process is concussive without elitism strategy due to the complex f7 function, though the optimization is convergent finally. The stability of optimization process for f4 is better, and it only has a small shock. For f2, it has tiny fluctuation due to the characteristics of function. The convergence is smooth and steady when the elitism strategy is added to the algorithm.

In this section, we compare ALREEDA with LSEDA-gl. LSEDA-gl is a robust univariate EDA, which is proposed for large scale optimization and has good performance. In LSEDA-gl, an effective sampling under mixed Gaussian and Levy probability distribution is introduced to balance optimization and learning. And a restart mechanism is used to stop some variables shrinking dramatically to zero solely. The optimization results of 100-D, 500-D, and 1000-D are indicated in Table 1. The results of ALREEDA on f3, f4, and f7 outperform LSEDA-gl. The results of the two algorithms are similar on f1, f5 and f6. For f2, the two algorithms are similar on 100-D. However, the performance of LSEDA-gl is better than ALREEDA on 500-D and 1000-D. The better and comparable results are marked by bold type in Table 1.

Table 1

Error values for functions 1-7 with D=100, 500 and 1000.

The convergent graph on 1000-D of the two algorithms is shown in Fig. (14). From the figure we can see the convergent process of ALREEDA is gentle than LSEDA-gl. In comparison, ALREEDA outperforms LSEDA-gl on f3, f4 and f7 even though the convergence is kindly. However, the performance of ALREEDA is wore on f2 than LSEDA-gl on 1000-D.

In order to indicate the performance of the algorithm, we compare the algorithm with other state of the art algorithms, like the cooperative coevolution DE algorithm MLCC [10Z. Yang, K. Tang, and X Yao, "Multilevel cooperative coevolution for large scale optimization", In: IEEE Congress on Evolutionary Computation, CEC 2008, IEEE Press: Hong Kong, China, 2008, pp. 1663-1670. Available from: http://sci2s.ugr.es/sites/default/files/ files/TematicWebSites/EAMHCO/contributionsCEC08/yang08mcc.pdf], the dynamic multi-swarm particle swarm optimizer DMS-PSO [43S. Zhao, J.J. Liang, P.N. Suganthan, and M.F. Tasgetiren, "Dynamic multi-swarm particle swarm optimizer with local search for large scale global optimization", In: Proceedings on Evolutionary Computation, IEEE: Hong Kong, China, 2008, pp. 3845-3852.
[http://dx.doi.org/10.1109/CEC.2008.4631320] ], the efficiency population utilization strategy based PSO EPUS-PSO [44S. Hsieh, T. Sun, C. Liu, and S. Tsai, "Solving large scale global optimization using improved particle swarm optimizer", In: Evolutionary Computation, IEEE: Hong Kong, China, 2008, pp. 1777-1784.], the improved PSO DEwSAcc [45A. Zamuda, J. Brest, B. Boskovic, and V. Zumer, "Large scale global optimization using differential evolution with self-adaptation and cooperative co-evolution", In: IEEE Congress on Evolutionary Computation, IEEE: Hong Kong, 2008, pp. 3718-3725.
[http://dx.doi.org/10.1109/CEC.2008.4631301] ] and jDEdynNP-F [46J. Brest, A. Zamuda, B. Boskovic, M.S. Maucec, and V. Zumer, "High-dimensional real-parameter optimization using self-adaptive differential evolution algorithm with population size reduction", In: IEEE Congress on Evolutionary Computation, IEEE: Hong Kong, 2008, pp. 2032-2039.
[http://dx.doi.org/10.1109/CEC.2008.4631067] ], and the multiple trajectory search algorithm MTS [47L. Tseng, and C. Chen, "Multiple trajectory search for large scale global optimization", In: Evolutionary Computation, IEEE: Hong Kong, 2008, pp. 3052-3059.], which took part in the competition of CEC’08 on large scale optimization. According to the competition result, a comparison is shown in Table 2 on D=100, 500 and 1000. In the table, in order to express the results clearly, the better and comparable results are marked by bold type.

Fig. (11) Convergence graph for function 1-7 (D=500, different NP).

Fig. (12) Convergence graph for function 1-7 (different D, the same NB).

Fig. (13) The optimization process partial enlarge graph of f7, f4 and f2.

Fig. (14) Convergence graph for function 1-7 (D=1000).

Table 2

Error values for functions 1-7 with D=100, 500 and 1000.

ALREEDA has a good performance on f1, f5, f6 and f7, and it is comparable with the MLCC and JDEdynNP-F with dimensions 100, 500 and 1000. For function f2, ALREEDA has a good performance than other algorithms on D=100. However, it becomes worse suddenly, although the results are comparable with other algorithms on D=500 and 1000. This may be caused by the using of simple probabilistic model, and also f2 is a non-separable function. MTS has a better performance than other algorithms for f2, f3 and f4 are multimodal functions with separable and non-separable characteristic, which are a challenge for any algorithms. For function f3, MTS has a good performance. MLCC outperformed other algorithms with high precisions on f4, especially for D=1000. Although the ALREEDA is not the best solver for f3 and f4, it has the similar performance as others, and it scales well on D=100, 500 and 1000. f7 is a very complex function, the performance of ALREEDA is promising.