F Teytaud
email: fteytaud@lri.fr

O Teytaud Tao

Why one must use reweighting in Estimation Of Distribution Algorithms

Keywords: G, 1, 6 [Optimization]: Unconstrained optimization; G, 3 [Probability and statistics]: Multivariate statistics

We study the update of the distribution in Estimation of Distribution Algorithms, and show that a simple modification leads to unbiased estimates of the optimum. The simple modification (based on a proper reweighting of estimates) leads to a strongly improved behavior in front of premature convergence.

INTRODUCTION

The American election of 1936 is famous due to the error, in a very important part of American history [2], in the poll organized by the Literary Digest [1], in spite of a huge sampling -2 millions of questionnaires. This error was due to the absence of reweighting: the readers of the Literary Digest were much more often Republicans than Democrats, and were therefore much more likely to vote for Landon; whereas Roosevelt finally got 61 % of votes, the Literary Digest predicted a comfortable win for Landon. The good result was predicted by George Gallup, with only a much smaller sample (50 000 people).

Estimation of Distribution Algorithms are similar to polls: they are based on samplings. However, to the best of our knowledge, the reweighting has not been experimented or analyzed in this context. This paper is based on this idea.

This paper considers the adaptation of the distribution in Estimation of Distribution Algorithms (EDA). EDA include e.g. UMDA [START_REF] Mühlenbein | Evolutionary computation and Wright's equation[END_REF], Compact Genetic Algorithm [START_REF] Harik | The compact genetic algorithm[END_REF], Population-Based Incremental Learning [3], Relative Entropy [START_REF] Mühlenbein | The estimation of Dim[END_REF], Cross-Entropy [4] and Estimation of Multivariate Normal Algorithm (EMNA) [START_REF] Larranaga | Estimation of Distribution Algorithms[END_REF]. We will here focus on EMNA (cf Algorithm 1).

There are several advantages in EDA; simplicity, possible use of prior knowledge, compliance with mixed (continuous/discrete) spaces, applicability for multimodal optimization (with multimodal distributions calibrated by Expectation-Maximization algorithms); importantly, EDA are often nearly parameter free. However, an issue is premature convergence [START_REF] Shapiro | Drift and scaling in estimation of distribution algorithms[END_REF][START_REF] Grahl | The correlation-triggered adaptive variance scaling idea[END_REF][START_REF] Liu | Model learning and variance control in continuous edas using pca[END_REF][START_REF] Posík | Preventing premature convergence in a simple eda via global step size setting[END_REF]: for example, with a Gaussian EDA, if the initial point is too far from the target (formally, if the squared distance to the optimum divided by the initial variance is too large), then the EDA might diverge.

Algorithm 1 The EMNA algorithm.

Initialize σ ∈ R, x ∈ R N . while Halting criterion not fulfilled do for l = 1..λ do z l = σN l (0, Id) x l = x + z l y l = f (x l ) end for Sort the individuals by increasing fitness;

y (1) < y (2) < • • • < y (λ) . z avg = 1 µ µ X i=1 z (i) σ = v u u u t µ X i=1 ||z (i) -z avg || 2 µ×N x = x + z avg end while

VARIANCE REDUCTION

Random points have a lot of nice features, in particular the fact that it's easier to avoid bias with simple Monte-Carlo. How, several forms of variance reductions are possible without introducing any bias. These techniques are also termed importance sampling in recent literature, referring to either correcting an incorrect distribution (as in the present paper), or improving the variance by changing the distribution.

• Using symmetries of the problem provides huge improvements. This was very clearly shown with Buffon's needle. The problem in Buffon's needle consists in estimating the probability for a needle thrown at random on a table with parallel lines on it, to cross at least one of the lines (Buffon's needle can be used for approximating π, yet it's essentially a toy problem nowadays).

The interesting point is that this probability is more quickly estimated by launching something else than a needle: launch (several times) a cross, made of two needles joint at their middle, estimate the probability for this cross to meet at least one line, apply some algebra to the frequency, and you'll find a good estimate, much more precise than if you had simply launched a needle. Using such symmetries (using e.g. so-called antithetic variables) was exploited successfully in the field of EDA in [START_REF] Gelly | On the ultimate convergence rates for isotropic algorithms and the best choices among various forms of isotropy[END_REF], and in various other fields as welle.g. [START_REF] Salmerón | Importance sampling in bayesian networks using antithetic variables[END_REF] for Bayesian Networks. We will not take care of this specific techniques, which can be used simultaneously as (and independently of) the reweighting that we propose below.

• Stratification (related, but not equivalent to, quota sampling -in quota sampling the distribution in each "quota" is not specified whereas in stratification it is): choosing a partition of the domain, and sample independently each part of the partition; then, each point should be reweighted so that its weight is consistent with its real probability. Stratification can very strongly reduce the variance of the estimate. When you can't sample with the target distribution (as well as in EDA, when we use Gaussian sampling whereas we would like to sample uniformly -there's just no uniform distribution such that we are sure that for λ sufficiently large we will go sufficiently far away from the parent to cover the optimum), then reweighting is the key solution for removing the bias.

Discussions related to variance reduction can be found in [5,[START_REF] Gross | Fundamentals of Queueing Theory[END_REF]. Some interesting historical and other related elements can be found in [START_REF] Schruben | Pseudorandom number assignment in statistically designed simulation and distribution sampling experiments[END_REF][START_REF] Morgan | Elements of simulation[END_REF].

WHY EMNA IS NOT (ALWAYS) CON-SISTENT AND HOW TO MAKE IT CON-SISTENT

In an update step of EMNA, we generate x1, . . . , x λ , independently identically distributed (i.i.d). Typically the distribution is Gaussian, centered on x and with step-size σ. We consider this simple case, but the adaptation to other distributions is straightforward. y1, . . . , y λ are the fitness values of the xi; for some unknown fitness function f , yi = f (xi). We then consider the sorted population; x (1) , . . . , x (λ) have fitness y (1) ≤ y (2) ≤ • • • ≤ y (λ) and are a permutation of x1, . . . , xn. We will assume for the sake of clarity that all fitness values are distinct, but the result does not depend on this assumption.

Then we define

x = 1 µ µ X i=1 x (i) . (1)
x is the estimate of the location of the optimum; it's the center of the next distribution. With the notation Êλ (resp. Pλ ) for the empirical averaging (resp. frequency) operator for the λ points (x1, y1),. . . ,(x λ , y λ ), this is equivalent to

x = Êλ (x|y ≤ y (µ) ).
µ is usually equal to λ/4; we will consider, for simplicity, that λ = 4µ, i.e. that λ can be divided by 4, but this assumption is not central; µ/λ might just converge to some limit in ]0, 1[ as λ is large.

We now analyze conditions under which EDA are consistent.

Theorem 1 (Consistency of EDA). Consider x1, x2, x3, x4, . . . i.i.d with density d(x), and consider y1, y2, y3, y4, . . . such that the (xi, yi) are independently and identically distributed as (x, y) with y = f (x). We assume that:

1. The repartition function P (y ≥ t) is continuous at 1/4 and the 1 4 -quantile of y is well defined, i.e. there is one and only one t such that P (y

≤ t) = 1 4 . 2. If I ⊂ R is compact then x|y ∈ I has bounded range; 3. ∀t > 0, E λ x|y < t = arg min f .
Then, with µ → ∞ and λ = 4µ, x → arg min f .

Remarks:

• Condition 1 is stronger than necessary. The proof is simpler with this condition, which is enough to ensure that (i) there is a 1 4 quantile (ii) there are not so many points very close to this quantile. A careful reduction of this assumption is beyond the scope of this paper.

• Condition 2 in particular holds if the function is coercive. With coercive functions, for ||x|| sufficiently large, f (x) > sup I.

• Condition 3 is the important condition for this paper: it looks like a simple technical assumption, but it does not hold in many important cases, as discussed later.

Proof: Assumption 1 implies that the 1 4 quantile of y is well defined. The convergence of an empirical quantile (here y (µ) ) to the quantile (and the existence of this quantile) holds almost surely as soon as this quantile is well defined (classical corollary of Kolmogorov-Smirnov's theorem). So assumption 1 implies y (µ) → t;

(2) for the only t such that P (y ≤ t) = 1 4 . Consider some fixed 1 > ǫ > 0, assumption 1 implies that

Pλ (y ∈]t -ǫ, t + ǫ[) ≤ k(ǫ) (3)
for some function k(.) such that limǫ→0 k(ǫ) = 0. Proof of Eq. 3: We need k(ǫ) such that P (y ≤ t + ǫ) -P (y ≤ t -ǫ) ≤ k(ǫ); we can just set k(ǫ) = P (y ≤ t + ǫ) -P (y ≤ t -ǫ). By the continuity assumption, k(ǫ) goes to 0 as ǫ → 0. (proof of Eq. 3)

Let's define

I = [t -ǫ, t + ǫ]. Then,
• for λ sufficiently large, by Eq. 3,

P (y ∈ I) ≤ k(ǫ); (4)
• Thanks to assumption 2,

Êλ (x|y ∈ I) < k ′ (5)
for some constant k ′ independent of ǫ. Now, consider the following equation for x.

˛Ê λ x|y < y (µ) -Êx|y < t ˛≤ P (y ∈ I) ˛Ê λ (x|y ∈ I) -Êλ x|y < t (6)
Eqs. 4 and 5 plugged in Eq. 6 lead to:

Êλ x|y < y (µ) = Êλ x|y < yt + O(k(ǫ)). (7)
The limit of Eq. 7 for ǫ → 0 is Êx|y < yµ → E λ x|y < t, which in turn leads to

x → E λ x|y < t.

Eq. 8 and assumption 3 conclude the proof.

Optimum (center of the level set)

Estimate of the optimum (biased towards the center of the Gaussian)

Level set A possible solution consists in using a uniform distribution in a sufficiently large ball. However, Gaussian numbers are quite comfortable: they provide arbitrarily large values, whilst preserving small variance. How to have the best of both worlds ? The idea of reweighting consists in:

• sampling with your most comfortable distribution (here the Gaussian), with density c(.);

• be aware of the target distribution, i.e. the distribution with which there would be no bias; let t(.) be its density (here the uniform distribution);

• reweighting the points in order to correct the bias; the weight of a point x is the following ratio:

w(x) = t(x) c(x) .
The implementation for EDA is quite straightforward: replace the weight 1 by this ratio between the target density and the density used for sampling. We will define the corresponding algorithm and experiment it in the next section. Interestingly, this ratio is nearly constant among selected points when they are close to the center of the optimum, but not at all when the selected points are far from the optimum, i.e. precisely in the cases in which EDA are sensitive to premature convergence.

EXAMPLES AND COUNTER-EXAMPLES

Let's discuss the assumption of this theorem. Assumptions 1 and 2 are weak assumptions, implies by e.g. smoothness and coercivity assumptions. Assumption 3 is seemingly quite natural but indeed does not hold in many cases. This section is devoted to examples and counter-examples for assumption 3. Thanks to reweighting we can come back to the averages for the uniform distribution; so we have to see cases in which the center of level sets is (resp. is not) the optimum. Essentially, we will see that in many cases the problem is not solved; essentially we can deal with local convergence. As discussed in section 6, reweighting makes sense in many cases; however, the proof as made in this paper considers only consistency of the estimation of the optimum.

Quadratic functions

Let's first consider quadratic functions.

• The standard case of a quadratic positive definite function (Fig. 2) is well handled by the theorem. The optimum is at the center of ellipsoids and therefore assumption 3 of the theorem clearly holds for the uniform distribution -for a Gaussian distribution, we have to reweight the points in order to ensure consistency.

Typically, if the level sets are concentric ellipsoids, then E λ x|y < t is equal to the optimum when the distribution is uniform, but certainly not if the distribution is Gaussian -in the case of a Gaussian distribution, the points are much more densely distributed close to the center of the distribution (see Figure 1); therefore, the estimate is not consistent, and increasing λ does not solve the problem.

• The degenerated case (Figure 3) in which the level sets are infinite rectangles, with a line of optima, is not directly handled; we have only considered the case in which the optimum is a point and the level sets are bounded (assumption 2). We conjecture that the result should be nearly preserved but have no proof of it. • The case of a saddle point (Figure 4) is not well handled also -the level sets are not ellipsoids at all, and we have to estimate level sets and not just the location of the optimum with some x (see however section 6, first point).

Figure 4: Quadratic function with saddle point. This is a multimodal case. This is not handled by the theorem; we do not have to estimate the position of the optimum but the level set.

Ellipsoid level sets of non-quadratic functions

We have discussed quadratic functions; this is not equivalent to ellipsoid level sets. We point out that whenever the sampling is made thanks to the uniform distribution, the fact that all level sets are ellipsoids is not enough for ensuring that the center of a level set is the optimum. Consider for example the following fitness function:

f -1 (t) = {cos(0.9t) + t * cos(u) + sin(0.9t) + t * sin(u); u ∈ [0, 2π]} for t ≥ 0.
This fitness (see level sets in Fig. 5) has optimum in (1, 0), which is the center of almost no level set for the uniform distribution. For a fitness as on Fig. 5, correcting the sampling bias in order to "recover" the uniform sampling is not enough; i.e. the reweighting technique that we propose below is not proved consistent in that case. However, in spite of the lack of proof in the general case, we guess that our technique also avoids premature convergence in such cases. To the best of our knowledge, proofs of consistency in such nonquasi-convex functions are still very rare (see however [START_REF] Teytaud | Conditionning, halting criteria and choosing lambda[END_REF]).

Other cases; why we try to get to much with the consistency theorem

We have shown several counter-examples. However, these counter-examples are cases in which x does not converge to the optimum as λ → ∞, in one generation; this is not necessary for the convergence of the EDA. Therefore we require too much. For example, the two very similar Figures 6 and7 have probably very similar behaviors; but only the first one is handled, due to the dissymmetry of case 7.

EXPERIMENTAL RESULTS

This section is devoted to experimenting the convergence of the weighted version of EMNA in a framework in which the baseline version that not converge, namely an initial point far from the optimum and a small initial step-size. In this section we use Algorithm 2 (to be compared with Algorithm 1, the baseline EMNA).

Comparison on the sphere function with poor initialization: reweighting avoids premature convergence

First, algorithm 2 is now compared to Algorithm 1. The results are averaged over 11 runs in various dimensions. The initial point is (1, 1, . . . , 1), the number of generations is 25⌊D 3/2 ⌋ and the initial step-size is σ = 0.1. The population size is λ = 10D⌊ √ D⌋ in Fig. 8, λ = 10 × D 2 in Fig.

f (1) < f (2) < • • • < f (λ) . z avg = 1 P µ i w (i) µ X i=1 w (i) z (i) σ = v u u u t µ X i=1 w(i)||z (i) -z avg || 2 P µ i=1 w(i)×N
y = y + z avg end while 9, λ = 10 × D 3 in Fig. 10, with D the dimension. The "scores" are D log(x)/n, where n is the number of generations, x is the estimate of the optimum at the end of the run -this is the approximation of the normalized asymptotic convergence rate. 0 means a premature convergence, and a significantly negative number means a linear convergence; we have linear convergence in all cases with reweighting, and never without, as shown in Figure 10.

EMNA with reweighting avoids premature convergence even without centering assumption

We claimed earlier that assumption 3 could probably be relaxed, and in particular that even if fitness functions as in Figure 11 (all but the first plot) do not verify the assumption 3, it should lead to a linear convergence as for the "spherical" fitness function (first plot, in the same figure). Figure 12 shows that high values of x make the problem more difficult,

t + t × (1 + 0.2 * cos(8u)) × cos(u), t × (1 + 0.2 * cos(8u))sin(u); u ∈ [0, 2π]}.
Due to dissymmetry, the center of the level set if not the optimum and the condition of the theorem are not met. However, it is clear that the algorithm will converge in this case as well as we don't have to be exactly centered on the optimum. This will be discussed on a similar problem in section 5, see in particular Fig. 11.

but that the convergence holds anyway (Fig 12).

Looking closely at the evolution of σ with poor initialization

We will now investigate the evolution of σ in a restricted setting, namely only one offspring, in order to see if increasing the number of points solves the problem. Figure 13 clearly shows that increasing λ is useless for avoiding premature convergence when there's no reweighting, whereas Figure 14 shows that when weighting is used, premature convergence is consistently avoided when λ is sufficiently large. This shows the consistency of the weighted approach and the weakness of the unweighted approach.

DISCUSSION

We have proved a theorem of consistency of the estimation of optima in EDA thanks to reweighting. This proof has the following limitations:

• We have only shown that the algorithm is consistent for estimating the optimum when we know the distribution for which the optimum is the expected value in a level set (typically, uniform distribution if the optimum is at the center of the level sets). This shows for a Gaussian EDA, reweighting ensures that ellipsoid level sets are consistently estimated. This does not say anything for the approximate level sets for other distributions. However, this is only technical and the principle of reweighting can be extended far beyond this. This will be the subject of an extension of this work; the goal would be the proof of a property of type "with reweighting, level sets are consistently estimated for some metric". statistics could be used for improving the variance of estimates, e.g. bootstrap [START_REF] Efron | An introduction to the bootstrap[END_REF]6], confidence regions for M-statistics [START_REF] Van Der Vaart | Weak Convergence and Empirical Processes, With Applications to Statistics[END_REF]. This is not done in this paper. Interestingly, bootstrap can be used for estimating confidence intervals as well, and can be used in difficult cases (i.e. beyond the simple conditional expectation studied in this paper).

• Other classical tricks for improving samplings consist in using quasi-random sequences in the sampling [START_REF] Niederreiter | Random Number Generation and Quasi-Monte Carlo Methods[END_REF][START_REF] Owen | Quasi-Monte Carlo sampling[END_REF]; quasi-random sampling is an active area of research with strong recent improvements in large dimension [START_REF] Ecuyer | Recent Advances in Randomized Quasi-Monte Carlo Methods[END_REF][START_REF] Sloan | When are quasi-Monte Carlo algorithms efficient for high dimensional integrals[END_REF]. Such improvements have already been tested for mutations in evolutionary algorithms [START_REF] Teytaud | DCMA: yet another derandomization in covariance-matrix-adaptation[END_REF][START_REF] Teytaud | When does quasi-random work?[END_REF]. Also, antithetic variables are easy to use in EDA, see e.g. [START_REF] Gelly | On the ultimate convergence rates for isotropic algorithms and the best choices among various forms of isotropy[END_REF] • In this paper, we considered reweighting for correcting a bias (typically, Gaussian distribution instead of uniform distribution for estimating the mass center of an ellipsoidal level set used. In the case of integration, it has been pointed out that sometimes reweighting points or using them in a more complicated manner than simple weighted averages is more important than well distributing them [START_REF] Ecuyer | Recent Advances in Randomized Quasi-Monte Carlo Methods[END_REF]. For example, one can weight points proportionally to the probability of their Voronoi cell. With such a technique, [START_REF] Yakowitz | Weighted monte carlo integration[END_REF] points out that in dimension 1 and 2 respectively, the integration error can be reduced from 1/ √ n (random points) or log(n) d /n (quasirandom points) to O(1/n 4 ) and O(1/n 2 ) on twice differentiable functions. Unfortunately, such results only hold in small dimensions, and [START_REF] Yakowitz | Weighted monte carlo integration[END_REF] points out that possibly, in high dimension, naive Monte-Carlo methods have some form of optimality among various possible uses of (unbiased) random i.i.d samplings.

We think that the results can be extended essentially by taking into account that we do not have to ensure that the center of the level set is the optimum, but only that we reduce the size of the sampling whilst keeping a good density around the optimum. We should not loose in this process the following strong points of the approach presented in this paper:

• We don't need any convexity or quasi-convexity assumption. This is quite important as it is quite natural to take care of non quasi-convex fitness functions; to the best of our knowledge, only this paper and [START_REF] Teytaud | Conditionning, halting criteria and choosing lambda[END_REF] have considered non-quasi-convex fitness functions.

• The results have this advantage that the proposed modification is quite simple, and can be used in many cases: it is just a reweighting. Also, this is not only theory for theory (even if theory for theory is interesting in many cases); we can implement the modification and have immediate very clear improvements.

Figure 1 :

 1 Figure 1: Illustration of the bias induced by selection in the case of unweighted points and nonuniform (e.g. Gaussian) distribution.

Figure 2 :

 2 Figure 2: Standard cases of quadratic definite functions with or without bad conditioning. This case is handled by the theorem.

Figure 3 :

 3 Figure 3: Degenerated quadratic function; one eigenvalue is null.

Figure 5 :

 5 Figure 5: Level sets of a fitness function which does not satisfy assumption 3 whenever the level sets are ellipsoids.

Figure 6 :

 6 Figure 6: Here, the t level set is {t×(1+0.2 * cos(8u))× cos(u), t × (1 + 0.2 * cos(8u))sin(u); u ∈ [0, 2π]}. Thanks to symmetries, the assumption 3 is verified: the center of the level set is the optimum. Algorithm 2 The EMNA algorithm with weighted averages. Initialize σ ∈ R, y ∈ R N . while Halting criterion not fulfilled do for l = 1..λ do z l = σN l (0, Id) y l = y + z l f l = f (y l ) end for Let w(i) = 1/density(xi) with density the density of the distribution used for generating the offspring. Sort the indices by increasing fitness; f (1) < f (2) < • • • < f (λ) .

Figure 7 :

 7 Figure 7: Here, the t level set is { 1 10t + t × (1 + 0.2 * cos(8u)) × cos(u), t × (1 + 0.2 * cos(8u))sin(u); u ∈ [0, 2π]}.Due to dissymmetry, the center of the level set if not the optimum and the condition of the theorem are not met. However, it is clear that the algorithm will converge in this case as well as we don't have to be exactly centered on the optimum. This will be discussed on a similar problem in section 5, see in particular Fig.11.

  ). Other forms of weights can be Figure 10: Normalized convergence rates for various dimensionalities, λ = 10D 3 and p-value of significance. Results are averaged over 11 runs; the fitness is the sphere function and the initialization is poor; see text for more details.

	Dimension Score without Score with P-value
		reweighting	reweighting
	3	-0.0121133	-2.12388	0
	4	-0.0150117	-2.28472	0
	5	-0.0114129	-2.39826	0
	6	-0.0135526	-2.54386	0
	7	-0.0152274	-2.65001	0
	8	-0.0165787	-2.73935	0
	9	-0.0181852	-2.83664	0
	10	-0.0144745	-2.80969	0
	11	-0.0160375	-2.99463	0
	12	-0.0173464	-3.06692	0
	13	-0.0186479	-3.13602	0
	14	-0.0195956	-3.20573	0
	15	-0.020884	-3.2708	0
	16	-0.022177	-3.33649	0

Acknowledgements

REMOVED FOR DOUBLE BLIND

Figure 11: Fitness function with level set f -1 (t) = {(t(x+cos(u), tsin(u); u ∈ [0, 2π]} for x = 0, 0.2, 0.4, 0.6, 0.8 respectively. With x = 0, the third assumption of the theorem does not hold; however, results in Fig. 12 clearly show that EMNA with reweighting solves this fitness function even with poor initialization. We check the probability (percentage) that σ increases (therefore leading to convergence). We see that this probability does not increase but decreases to 0, consistently with theory.

On the other hand, we can see on Figure 14 that with reweighting the probability converges to 1 as λ becomes large. All probabilities are estimated on 30 runs.


F Teytaud
email: fteytaud@lri.fr

O Teytaud

On the parallel speed-up of Estimation of Multivariate Normal Algorithm and Evolution Strategies

Motivated by parallel optimization, we experiment EDA-like adaptation-rules in the case of λ large. The rule we use, essentially based on estimation of multivariate normal algorithm, is (i) compliant with all families of distributions for which a density estimation algorithm exists (ii) simple (iii) parameter-free (iv) better than current rules in this framework of λ large. The speed-up as a function of λ is consistent with theoretical bounds.

Introduction

Evolution Strategies (ES [START_REF] Rechenberg | Evolutionstrategie: Optimierung Technischer Systeme nach Prinzipien des Biologischen Evolution[END_REF]) are a robust optimization tool, known for its robustness (in particular in front of local minima) and simplicity. It is also known as suitable for parallel optimization, because it is population-based. However, it has been pointed out recently in [START_REF] Beyer | Covariance matrix adaptation revisited -the CMSA evolution strategy[END_REF] that usual step-size adaptation rules were far from being efficient for λ large, e.g. λ = 4N 2 where λ is the population size and N is the dimension of the search space.

The case of λ = 4N 2 as in [START_REF] Beyer | Covariance matrix adaptation revisited -the CMSA evolution strategy[END_REF] or λ >> 4N 2 is not purely theoretical. In our recent applications, we usually optimized on a cluster of 368 cores, and we recently organized an optimization on several clusters of a few hundreds cores on each cluster. With λ = 2000 cores, N = 22 satisfies λ = 4N 2 . Moreover, in many cases we have to launch several jobs simultaneously (more than the number of cores) in order to save up scheduling time, leading to λ much larger than the number of cores, and in robust optimization N ≤ 50 is far from being trivial.

In this paper, we: (i) describe the main step-size adaptation rules in the literature (section 2); (ii) experiment step-size adaptation rules for various values of λ (section 3); (iii) discuss the results in section 4. In all the paper, we assume that the reader is familiar with standard notations around ES (see e.g. [START_REF] Rechenberg | Evolutionstrategie: Optimierung Technischer Systeme nach Prinzipien des Biologischen Evolution[END_REF][START_REF] Beyer | The Theory of Evolution Strategies[END_REF] for more information) and we consider (µ/µ, λ)algorithms, i.e.: (i) at each iteration of the algorithm, λ points are generated according to some (usually but not necessarily) Gaussian distribution; (ii) the fitness values of these λ points are computed; (iii) then, the µ best points according to the fitness function are selected; averages are then w.r.t this subsample of the µ best points. All averages here are unweighted averages, but methods used in this paper can be applied in weighted cases. N (0, Id) and related notations (e.g. N i (0, Id)) denote standard multivariate Gaussian random variables with identity covariance matrix.

Methods

A central issue in ES is the adaptation of the distribution (step-size and beyond). The one-fifth rule has been successfully introduced in [START_REF] Rechenberg | Evolutionstrategie: Optimierung Technischer Systeme nach Prinzipien des Biologischen Evolution[END_REF] for this, and several other rules have been proposed later in the literature; some main rules are detailed in the rest of this section.

In this paper, we are considering minimization problems.

Cumulative

Step-size Adaptation (CSA)

Cumulative step-size adaptation has been proposed in [START_REF] Hansen | Completely derandomized self-adaptation in evolution strategies[END_REF]; essentially, this method compares the path followed by the algorithm to the path followed under random selection: if the followed path is larger, then the step size should increase. The detailed algorithm is presented in Algorithm 1.

Algorithm 1 Cumulative step-size adaptation.

Initialize σ ∈ R, y ∈ R N . while halting criterion not fulfilled do for i = 1..λ do si = Ni(0, Id) yi = y + σsi fi = f (yi) end for Sort the individuals by increasing fitness;

f (1) < f (2) < • • • < f (λ) . s avg = 1 µ P µ i=1 s (i) y = y + σs avg pσ = (1 -c)pσ + p µc(2 -c)s avg σ = σ exp[ ||pσ ||-χ N Dχ N ] end while where χ N is √ N ×(1-1 4.0×N + 1 21.0×N 2 )
. Following [START_REF] Hansen | Verallgemeinerte individuelle Schrittweitenregelung in der Evolutionsstrategie. Eine Untersuchung zur entstochastisierten,koordinatensystemunabhängigen Adaptation der Mutationsverteilung[END_REF], we choose c = 1 √ N and D = √ N . A main weakness of this formula is its moderate efficiency, for λ large, as pointed out in [START_REF] Beyer | Covariance matrix adaptation revisited -the CMSA evolution strategy[END_REF]. [START_REF] Beyer | Covariance matrix adaptation revisited -the CMSA evolution strategy[END_REF] therefore proposes mutative self-adaptation in order to improve the convergence rate as a function of λ; this mutative self-adaptation is presented below.

Mutative Self-Adaptation (SA)

Mutative self-adaptation has been proposed in [START_REF] Rechenberg | Evolutionstrategie: Optimierung Technischer Systeme nach Prinzipien des Biologischen Evolution[END_REF] and [START_REF] Schwefel | Adaptive Mechanismen in der biologischen Evolution und ihr Einfluss auf die Evolutionsgeschwindigkeit[END_REF], and extended in [START_REF] Beyer | Covariance matrix adaptation revisited -the CMSA evolution strategy[END_REF] for improving the convergence rate for λ large. The algorithm is simple and provides the state of the art results for λ large; it is presented in Algorithm 2.

Algorithm 2 Mutative self-adaptation. τ is equal to 1/ √ N ; other variants have been tested (e.g. 1/ √ 2N which is sometimes found in papers) without improvement in our case of λ large.

Initialize σ avg ∈ R, y ∈ R N . while Halting criterion not fulfilled do for i = 1..λ do σi = σ avg e τ N i (0,1) zi = σiNi(0, Id) yi = y + zi fi = f (yi) end for Sort the individuals by increasing fitness;

f (1) < f (2) < • • • < f (λ) . z avg = 1 µ P µ i=1 z (i) σ avg = 1 µ P µ i=1 σ (i) y = y + z avg end while

Statistical Step-size Adaptation (SSA)

We here propose simple step-size adaptation rules, inspired by Estimation of Distribution Algorithms (EDA), e.g. UMDA [START_REF] Mühlenbein | Evolutionary computation and Wright's equation[END_REF], Compact Genetic Algorithm [START_REF] Harik | The compact genetic algorithm[END_REF], Population-Based Incremental Learning [START_REF] Baluja | Population-based incremental learning: A method for integrating genetic search based function optimization and competitive learning[END_REF], Relative Entropy [START_REF] Mühlenbein | The estimation of distributions and the minimum relative entropy principle[END_REF], Cross-Entropy [START_REF] Cai | Cross entropy and adaptive variance scaling in continuous eda[END_REF] and Estimation of Multivariate Normal Algorithms (EMNA) [START_REF] Larranaga | Estimation of Distribution Algorithms[END_REF] (our main inspiration), which combine (i) the current distribution (possibly), (ii) statistical properties of selected points, into a new distribution. We show in this paper that forgetting the old estimate and only using the new points is a good idea in the case of λ large; in particular, premature convergence as pointed out in [START_REF] Shapiro | Drift and scaling in estimation of distribution algorithms[END_REF][START_REF] Grahl | The correlation-triggered adaptive variance scaling idea[END_REF][START_REF] Liu | Model learning and variance control in continuous edas using pca[END_REF][START_REF] Posík | Preventing premature convergence in a simple eda via global step size setting[END_REF] does not occur if λ >> 1 points are distributed on the search space with non-degenerated variance, and troubles around variance estimates for small sample size as in [START_REF] Donga | Unified eigen analysis on multivariate gaussian based estimation of distribution algorithms[END_REF] are by definition not relevant for us. Its advantages are as follows for λ large: (i) it's very simple and parameter free; the reduced number of parameters is an advantage of mutative self adaptation in front of cumulative step-size adaptation, but we show here that yet fewer parameters (0!) is possible and leads to better results, thanks to EMNA-like algorithms; (ii) we use it here in the case of a Gaussian, but it could easily be used also for any density estimation technique e.g. sums of Gaussians (i.e. no problem with multimodal distributions); (iii) it could easily be used for discrete spaces also. The algorithm is presented in Algorithm 3 for estimating only one step-size, but a straightforward extension is one step-size per axis or even a full covariance matrix (see section 3.3).

Experimental results

As we here focus on step-size adaptation-rules, we here do not consider complete Covariance-Matrix Adaptation (CMA) but only diagonal covariance matrix Algorithm 3 SSA. This is close to EMNA; we divide by µ×N in the adaptation rule (eq. defining σ) because the sum is over µ × N deviations (one per selected individual and per axis), and we add a trick against premature convergence. The constant K is here ∞ as this work is not devoted to premature convergence; however, we verified that the same convergence rates as those presented in this paper can be recovered with this more robust version.

Initialize σ ∈ R, y ∈ R N . while Halting criterion not fulfilled do for l = 1..λ do z l = σN l (0, Id) y l = y + z l f l = f (y l ) end for Sort the individuals by increasing fitness; f (1) < f (2) < • • • < f (λ) . z avg = 1 µ µ X i=1 z (i) if ||z avg || < Kσ then σ = v u u u t µ X i=1 ||z (i) -z avg || 2 µ×N else
σ = 2σ /** this avoids premature convergence in case of σ poorly chosen**/ end if y = y + z avg end while adaptation, but [START_REF] Ros | A simple modification in cma-es achieving linear time and space complexity[END_REF] has shown that this diagonal framework is indeed better in many (even non separable) cases, and we point out that all methods discussed in section 2 have anyway straightforward extensions in the framework of complete covariance matrix or covariance matrix by block. We first confirm results from [START_REF] Beyer | Covariance matrix adaptation revisited -the CMSA evolution strategy[END_REF] (superiority of mutative self-adaptation over cumulative step-size adaptation for λ large), in section 3.1; we then validate the statistical step-size adaptation rule in section 3.2. Table 1 presents the objective functions considered in this paper.

Validating SA for λ large

We here confirm results in [START_REF] Beyer | Covariance matrix adaptation revisited -the CMSA evolution strategy[END_REF], namely the moderate efficiency of CSA for λ large (at least under the various parameterizations we tested). Following [START_REF] Beyer | Covariance matrix adaptation revisited -the CMSA evolution strategy[END_REF], Figure 1 presents the numbers of iterations before a fixed halting criterion (f stop = 10 -10 ).

Validating the statistical step-size adaptation

We now present the main result of this paper, i.e. the efficiency of the statistical step-size adaptation rule. Results are presented in Fig. 2

f (y) = N X i=1 y 2 i Schwefel f (y) = N X i=1 ( i X j=1 yi) 2 Cigar f (y) = y 2 1 + 10000 × N X i=2 y 2 i
Table 1. Objective functions considered in this paper. µ = λ/4, µ = 1 respectively. Presented curves are N × log(||x -x * ||)/n as a function of λ, where: N is the dimensionality; x is the best point in the last λ offspring; x * is the optimum (x * = 0 for our test functions); n is the number of iterations before the halting criterion is met. Each run is performed until f (x) < 10e -50 . For the sake of statistical significance each point is the average of 300 independent runs.

Anisotropic case

The SSA algorithm can be adapted to the anisotropic case. We here consider the separable version, i.e. a diagonal covariance matrix for the Gaussian; this form of covariance matrix is intermediate between the full matrix and the matrix proportional to identity. It has been pointed out in [START_REF] Ros | A simple modification in cma-es achieving linear time and space complexity[END_REF] that this separable version is in many cases faster than the complete covariance matrix. We have in this case to determine one step-size per axis; see Algorithm 4.

Figure 5 presents the results of Algorithm 4. We can see that: (i) On the Schwefel function, the results are the same as in the case of the sphere function. The isotropic algorithms have, in this framework (non presented results), a result close to 0 (i.e. much worse). Therefore, the anisotropic step-size adaptation works for this moderately ill-conditioned function. (ii) On the Cigar function, the results are not exactly those of the sphere function, but almost; even on this much more ill-conditioned function, the anisotropic SSA works.

Discussion

First, we confirm the superiority of mutative self-adaptation over cumulative step-length adaptation, as already shown in [START_REF] Beyer | Covariance matrix adaptation revisited -the CMSA evolution strategy[END_REF], in the case of λ large. However, possibly, tuning CSA might lead to improved results for high values of λ; this point, beyond the scope of this paper, is not further analyzed here. Second, we experiment a simple statistical adaptation rule adapted from EMNA. This rule for step-size adaptation is (i) very simple (the most simple of all rules in section 2) (ii) fully portable to the full-covariance case (as well as the SA approach used in [START_REF] Beyer | Covariance matrix adaptation revisited -the CMSA evolution strategy[END_REF] for covariance matrix adaptation) (iii) computationally cheap (iv) highly intuitive (v) parameter-free (except K if we include the rule against premature convergence, which is only required if choosing σ sufficiently large is hard) (vi) efficient whenever λ is not known in advance as λ has no impact on the adaptation (important for fault-tolerance). It provides a speed-up of 100% (over mutative self-adaptation, which is already much better than CSA) on the sphere function for N = 3, decreasing to 33% for N = 10 and 25% for N = 100 (in the case λ = 150). The usual adaptation-rules use a combination of an old covariance matrix and a new estimate based on the current population -essentially, in SSA, we only keep the second term as λ large simplifies the problem. We point out that we only experimented the Gaussian case so that we can compare our results to standard rules, but our algorithm has the advantage that it can be adapted to all distributions for which the parameters can be estimated from a sample. Sums of Gaussians are a natural candidate for further work. Third, we show that the anisotropic version works in the sense that the convergence rate on the sphere was recovered with the Schwefel and the Cigar function.

Some by-products of our results are now pointed out, as a comparison with theoretical results in [START_REF] Teytaud | Lower bounds for evolution strategies using vcdimension[END_REF]: (i) The higher the dimensionality, the better the speedup for (µ/µ, λ)-algorithms; [START_REF] Beyer | The Theory of Evolutions Strategies[END_REF] has shown the linear speed-up as a function of λ as long as λ << N , and [START_REF] Teytaud | Lower bounds for evolution strategies using vcdimension[END_REF] has precised formally that the speed-up is linear until λ of the same order as the dimension. Roughly, a number of processors linear as a function of the dimension is efficient. This is visible on our experimental results. (ii) Also, (1, λ) algorithms (case µ = 1) have a less-than-linear (logarithmic) speed-up as a function of λ. This is visible in our experimental results. This is also consistent with [START_REF] Beyer | The Theory of Evolutions Strategies[END_REF] and [START_REF] Teytaud | Lower bounds for evolution strategies using vcdimension[END_REF]. (iii) We have proposed a very simple rule and got state-of-the-art results. This suggests that there is much to win by a refined work on the important case of λ large.

  , 3, 4 for µ = λ/2,

Fig. 1 .

 1 Fig.1. Performance of CSA and SA on the sphere function (number of iterations before f (x) < 10 -10 ) depending on the dimensionality. We confirm here the superiority of SA on cumulative step-size adaptation, at least for our version of SA and CSA. In these experiments, µ = 1 4 λ.

Fig. 2 .

 2 Fig.2. Log of distance to the optimum when the halting criterion (see text) is met, normalized by the dimension and the number of iterations. Results for µ = λ/2, on the sphere function. In all cases, SSA is the best rule for λ large.

Fig. 3 .

 3 Fig.3. Results for µ = λ/4 on the sphere function. All methods are better than for µ = λ/2. In all cases, SSA is the best rule for λ large.

Fig. 4 .

 4 Fig.4. Results for µ = 1 on the sphere function. SSA is not presented as it does not make sense for µ = 2. As shown by this figure (compared to Figs3 and 2), µ = 1 is quite weak for large dimension and the absence of SSA version for that case is therefore not a trouble.

Fig. 5 .

 5 Fig.5. We here present the results of the anisotropic version of the SSA algorithm (Algorithm 4). It has on ill-conditioned functions nearly the same performance as the isotropic version on the sphere, in particular in high dimension.

Acknowledgements We are grateful to H.-G. Beyer for useful email answers, and to A. Auger, N. Hansen, R. Ros for fruitful discussions. We also thank the MoGo-people (G. Chaslot, J.-B. Hoock, A. Rimmel) for interesting discussions around parallel applications of ES.

Algorithm 4 Anisotropic statistical step-size adaptation. We divide the average squared deviation by µ in the step-size adaptation rule (equation defining σ), instead of µ × N in the anisotropic case, because now the step-size is averaged over µ squared deviations only (one per selected individual). The subscript j refers to the axis.

while Halting criterion not fulfilled do for l = 1..λ do for ∈ {1, . . . , N } do z l,j = σjN l,j (0, 1) end for y l = y + z l f l = f (y l ) end for Sort the individuals by increasing fitness;

end for y = y + z avg end while


Merve Nur

Nicolas Amil

Christian Bredeche

Sylvain Gagné

Marc Gelly

Olivier Schoenauer

Teytaud

A Statistical Learning Perspective of Genetic Programming

This paper proposes a theoretical analysis of Genetic Programming (GP) from the perspective of statistical learning theory, a well grounded mathematical toolbox for machine learning. By computing the Vapnik-Chervonenkis dimension of the family of programs that can be inferred by a specific setting of GP, it is proved that a parsimonious fitness ensures universal consistency. This means that the empirical error minimization allows convergence to the best possible error when the number of test cases goes to infinity. However, it is also proved that the standard method consisting in putting a hard limit on the program size still results in programs of infinitely increasing size in function of their accuracy. It is also shown that cross-validation or hold-out for choosing the complexity level that optimizes the error rate in generalization also leads to bloat. So a more complicated modification of the fitness is proposed in order to avoid unnecessary bloat while nevertheless preserving universal consistency.

Most of the works cited earlier are in fact concerned with functional bloat, which is the most surprising, and the most disappointing kind of bloat. We will consider various levels of functional bloat: cases where

Introduction

This paper is about two important issues in Genetic Programming (GP), that is Universal Consistency (UC) and code bloat. UC consists in the convergence to the optimal error rate with regards to an unknown distribution of examples. A restricted version of UC is consistency, which focus on the convergence to the optimal error rate within a restricted search space. Both UC and consistency are well studied in the field of statistical learning theory. Despite their possible benefits, they have not been widely studied in the field of GP. Code bloat is the uncontrolled growth of program size that may occur in GP when relying on a variable length representation [START_REF] Koza | Genetic Programming: On the Programming of Computers by Means of Natural Selection[END_REF][START_REF] Langdon | The evolution of size in variable length representations[END_REF]. This has been identified as a key problem in GP for which there have been several empirical studies. However, very few theoretical studies addressed this issue directly. The work presented in this paper is intended to provide some theoretical insights on the bloat phenomenon and its link with UC in the context of GP-based learning taking a statistical learning theory perspective [START_REF] Vapnik | The Nature of Statistical Learning Theory[END_REF].

Statistical learning theory provides several theoretical tools to analyze some aspects of learning accuracy. Our main objective consists in performing both an in-depth analysis of bloat as well as providing appropriate solutions to avoid it. Section 2 shortly exposes issues of code bloat with GP. Section 3 and 4 present all the aforementioned results about code bloat avoidance and UC and propose a new approach ensuring both. Then, Section 5 provides some extensions of the previous theoretical results on the use of cross-validation and hold-out methodologies. Follows some experimental results in Section 6, illustrating the accuracy of the theoretical results. Section 7 finally concludes this paper with a discussion on the consequences of those theoretical results for GP practitioners and uncover some perspectives of work.

Code Bloat in GP

Due to length constraints, we do not introduce here some important theories around code bloat: introns, fitness causes bloat, and removal bias. The reader is refered to [START_REF] Banzhaf | Genetic Programming : an introduction[END_REF][START_REF] Blickle | Genetic programming and redundancy[END_REF][START_REF] Langdon | The evolution of size in variable length representations[END_REF][START_REF] Langdon | Fitness causes bloat: Mutation[END_REF][START_REF] Mcphee | Accurate replication in genetic programming[END_REF][START_REF] Nordin | Complexity compression and evolution[END_REF][START_REF] Soule | Exons and code growth in genetic programming[END_REF] for more informations around that. Some common solutions against bloat rely either on specific operators (e.g. size-fair crossover [START_REF] Langdon | Size fair and homologous tree genetic programming crossovers[END_REF], or different fair mutation [START_REF] Langdon | The evolution of size and shape[END_REF]), on some parsimony-based penalization of the fitness [START_REF] Soule | Effects of code growth and parsimony pressure on populations in genetic programming[END_REF] or on abrupt limitation of the program size such as the one originally used by Koza [START_REF] Koza | Genetic Programming: On the Programming of Computers by Means of Natural Selection[END_REF]. Also, some multi-objective approachs have been proposed [START_REF] Bleuler | Multiobjective genetic programming: Reducing bloat using SPEA2[END_REF][START_REF] Edwin | Reducing bloat and promoting diversity using multi-objective methods[END_REF][START_REF] Ekart | Maintaining the diversity of genetic programs[END_REF][START_REF] Luke | Lexicographic parsimony pressure[END_REF][START_REF] Silva | Dynamic maximum tree depth : A simple technique for avoiding bloat in tree-based GP[END_REF]. Some other more particular solutions have been proposed but are not widely used yet [START_REF] Ratle | Avoiding the bloat with probabilistic grammar-guided genetic programming[END_REF][START_REF] Zhang | Balancing accuracy and parsimony in genetic programming[END_REF]. Also, all proofs are removed due to length constraints. Readers familiar with mathematics like in e.g. [START_REF] Devroye | A Probabilistic Theory of Pattern Recognition[END_REF] should however be able to guess the main ideas.

Although code bloat is not clearly understood, it is yet possible to distinguish at least two kinds of code bloat. We first define structural bloat as the code bloat that necessarily takes place when no optimal solution can be approximated by a set of programs with bounded length. In such a situation, optimal solutions of increasing accuracy will also exhibit an increasing complexity (larger programs), as larger and larger code will be generated in order to better approximate the target function. This extreme case of structural bloat has also been demonstrated in [START_REF] Gustafson | Problem difficulty and code growth in genetic programming[END_REF]. The authors use some polynomial functions of increasing difficulty, and demonstrate that a precise fit can only be obtained through an increased bloat (see also [START_REF] Daida | What makes a problem GP-Hard? Analysis of a tunably difficult problem in genetic programming[END_REF] for related issues about problem complexity in GP). Another form of bloat is the functional bloat, which takes place when program length keeps on growing even though an optimal solution (of known complexity) does lie in the search space. In order to clarify this point, let us use a simple symbolic regression problem defined as follow: given a set S of test cases, the goal is to find a function f (here, a GP-tree) that minimizes the Mean Square Error (or MSE). If we intend to approximate a polynomial (e.g. 14 * x 2 with x ∈ [0, 1]), we may observe code bloat since it is possible to find arbitrarily long polynomials that gives the exact solution (e.g. 14x 2 +0 * x 3 +. . .), or sequences of polynomials of length growing to ∞ and accuracy converging to the optimal accuracy (e.g. P n (x) = 14x 2 + length of programs found by GP runs to infinity as the number of test cases runs to infinity whereas a bounded-length solution exists, and also cases where large programs are found with high probability by GP whereas a small program is optimal.

Another important issue is to study the convergence of the function given by GP toward the actual function used to generate the test cases, under some sufficient conditions and when the number of test cases goes to infinity. This property is known in statistical learning as Universal Consistency (UC). Note that this notion is slightly different from that of universal approximation, commonly referred in symbolic regression, where GP search using operators {+, * } is assumed to be able to approximate any continuous function. UC is rather concerned with the behavior of the algorithm when the number of test cases goes to infinity: the existence of a polynomial that approximates a given function at any arbitrary precision does not imply that any polynomial approximation built from a set of sample points will converge to that given function when the number of points goes to infinity. Or more precisely, UC can be stated informally as follows (a formal definition will be given later):

A GP setting corresponds to symbolic regression from examples if it takes as inputs a finite number of examples x 1 , . . . , x n with their associated labels y 1 , . . . , y n and outputs a program P n . Universal consistency holds if, when pairs (x 1 , y 1 ), . . . , (x n , y n ) (test cases) are identically independently distributed as the random variable (x, y), L(P n ) → L * where L(p) = P r(y = p(x)) and where L * = inf p measurable L(p). In all of this paper, Pr(.) denotes probabilities, as the traditional notation P (.) is used for programs.

Negative Results without Regularization and Resampling

Definition 1 precisely defines the programs space under examination. Theorem 1 evaluates its VC-dimension [START_REF] Vapnik | The Nature of Statistical Learning Theory[END_REF]. Many theorems, in the sequel, are based only on VC-dimensions and hold for other sets of programs as well.

It should be noted the mildness of the hypothesis behind our results. We consider any programs of bounded length, working with real variables, provided that the computation time is a priori bounded. Usual families of programs in GP verify this hypothesis and much stronger hypothesis. For example, usual treebased representations avoid loops and therefore all quantities that have to be bounded in lemma below (typically, number of times each operator is used) are bounded for trees of bounded depths. This is also true for direct acyclic graphs. We here deal with a very general case; much better constants can be derived for specific cases, without changing the fundamental results in the sequel of the paper.

Definition 1 (Set of programs studied). Let F (n, t, q, m, z) be the set of functions from R z-m towards {0, 1} which can be computed by a program with a maximum of n lines as follows:

(1) A run uses at most t operations. ( 2) Each line contains one operation among the followings:

-Operations α → exp(α) (at most q times); -Operations +, -, ×, and /; -Jumps conditioned on >, ≥, <, ≤, and =; -Output 0; -Output 1; -Labels for jumps; -Constants (at most m different); -Variables (at most z different, with z ≥ m).

We note F (n, t, q, m, z) as F for short when there is no ambiguity. The parameters n, t, q, m, z are then implicit.

The following property is central for our results.

Theorem 1 (Finite VC-dimension of the computing machine). Consider q ′ , t ′ and d ′ ≥ 0. Let F = F (n, t, q, m, z) be the set of programs described by Definition 1, where q ≤ q ′ , T (n, t, z) ≤ t ′ , and

1 + m ≤ d ′ . V Cdim(F ) ≤ t ′2 d ′ (d ′ + 19 log 2 (9d ′ )) ≤ (d ′ (q ′ + 1)) 2 + 11d ′ (q ′ + 1)(t ′ + log 2 (9d ′ (q ′ + 1))) If q = 0 (no exponential) then V Cdim(F ) ≤ 4d ′ (t ′ + 2).
Interpretation: The theorem demonstrates that interesting and natural families of programs have finite VC-dimension. Effective methods can associate a VCdimension to these families of programs.

We now consider how to use such results in order to ensure UC. First, we show why simple empirical risk minimization (i.e. minimizing the error observed without taking into account programs complexity) does not ensure consistency. More precisely, for some distribution of test cases and some i.i.d. (independent identically distributed) sequence of test cases {(x 1 , y 1 ), . . . , (x n , y n ), . . . }, there exists P 1 , . . . , P n , . . . such that ∀n ∈ N, ∀i ∈ {1, 2, . . . , n} P n (x i ) = y i , and however ∀n ∈ N Pr (P n (x) = y) = 0. This can be proved by considering that x is uniformly distributed in [0, 1] and y is a constant equal to 1. Then, consider P n , the program that compares its entry to x 1 , x 2 , . . . , x n , and outputs 1 if the entry is equal to x j for some j ≤ n, and otherwise outputs 0. With probability 1, this program output 0, whereas almost surely the desired output y is 1.

We therefore conclude that minimizing the empirical risk is not enough for ensuring any satisfactory form of consistency. Let's now show that structural risk minimization (i.e. taking into account a penalization for complex structures) can ensure UC and fast convergence when the solution can be written within finite complexity.

Theorem 2 (Universal consistency of genetic programming with structural risk minimization). Consider q k , t k , m k , n k , and z k increasing integer sequences. Define F k the set of programs with at most t k lines executed, z k variables, n k lines, q k exponentials, and

m k constants (F k = F (n k , t k , q k , m k , z k ) of Definition 1) and F = ∪ k F k . Then with q ′ k = q k , t ′ k = T (n k , t k , z k ), and
d ′ k = 1 + m k , define V k as: -If ∀k q k = 0, then V k = 4d ′ k (t ′ k + 2). -Otherwise, V k = (d ′ k (q ′ k + 1)) 2 + 11d ′ k (q ′ k + 1)(t ′ k + log 2 (9d ′ k (q ′ k + 1))).
Now given s test cases, consider P ∈ F minimizing L(P )+ 32 s V (P ) log(es), where V (P ) = V k where k is minimal such that P ∈ F k . Then, the generalization error, with probability 1, converges to L * ; moreover, if one optimal program belongs to F k , then for any s and ǫ such that V k log(es) ≤ sǫ 2 /512, the generalization error with s test cases is larger than L * + ǫ with probability at most ∆ exp(-sǫ 2 /128) + 8s

V k exp(-sǫ 2 /512) where ∆ = ∞ j=1 exp(-V j ).
Interpretation: This theorem shows that genetic programming for binary classification, provided that structural risk minimization is performed (i.e. if we optimize an ad hoc compromise between complexity of programs and accuracy on empirical data), is universally consistent and verifies some convergence rate properties.

We now prove the non-surprising fact that if it is possible to approximate the optimal function (the Bayesian classifier) without reaching it exactly, then the complexity of the program runs to infinity as soon as there is convergence of the generalization error to the optimal one.

Proposition 1 (Structural bloat in genetic programming).

Consider

F 1 ⊂ F 2 ⊂ F 3 ⊂ . . .
, where F V is a set of functions from X to {0, 1} with VC-dimension bounded by V . Consider (V (s)) s∈N a non decreasing sequence of integers and (P s ) s∈N a sequence of functions such that P s ∈ F V (s) .

Define L V = inf P ∈FV L(P ) and V (P ) = inf{V ; P ∈ F V } and suppose that

∀V L V > L * . Then, L(P s ) s→∞ -→ L * =⇒ V (P s ) s→∞ -→ ∞ .
Interpretation: This is structural bloat: if the space of programs approximates but does not contain the optimal function and cannot approximate it within bounded size, then bloat occurs. Note that for any F 1 , F 2 , . . . , the assumption ∀V L V > L * holds simultaneously for all V for many distributions, as we consider countable unions of families with finite VC-dimension (e.g. see [6, chap. 18]).

We now show that, even in cases in which an optimal short program exists, the usual procedure (known as the method of Sieves; see also [START_REF] Silva | Dynamic limits for bloat control: Variations on size and depth[END_REF]) defined below, consisting in defining a maximum VC-dimension depending upon the sample size and then using a family of functions accordingly, leads to bloat.

Theorem 3 (Bloat with the method of Sieves). Let F 1 , . . . , F k , . . . be nonempty sets of functions with finite VC-dimensions V 1 , . . . , V k , . . . , and let F = ∪ n F n . Then given s i.i.d. test cases, consider P ∈ F s minimizing the empirical risk L in F s .

From theorems about the method of Sieves, we already know that if V s = o(s/ log(s)) and V s → ∞, then Pr L( P ) ≤ L( P ) + ǫ(s, V s , δ) ≥ 1 -δ and almost surely L( P ) → inf P ∈F L(P ).

We now state that if V s → ∞, and noting V (P ) = min{V k ; P ∈ F k }, then ∀V 0 , δ 0 > 0, ∃ Pr, a distribution of probability on X and Y , such that ∃g ∈ F 1 such that L(g) = L * , and for s sufficiently large Pr V ( P ) ≤ V 0 ≤ δ 0 .

Interpretation: The result in particular implies that for any V 0 , there is a distribution of test cases such that ∃g; V (g) = V 1 and L(g) = L * , with probability 1, V ( P ) ≥ V 0 infinitely often as s increases. This shows that bloat can occur if we use only an abrupt limit on code size, if this limit depends upon the number of test cases ( a fortiori if there's no limit). Note that this result, proved thanks to a particular distribution, could indeed be proved for the whole class of classification problems for which the conditional probability of Y = 1 (conditionally to X) is equal to 1 2 in an open subset of the domain.

Universal Consistency without Bloat

In this section, we consider a more complicated case where the goal is to ensure UC, while simultaneously avoiding non-necessary bloat. This means that an optimal program does exist in a given family of functions and convergence towards the minimal error rate is performed without increasing the program complexity. This is achieved by: i) merging regularization and bounding of the VC-dimension, and ii) penalization of the complexity (i.e. length) of programs by a penalty term R(s, P ) = R(s)R ′ (P ) depending upon the sample size and the program. R(., .) is user-defined and the algorithm looks for a classifier with a small value of both R ′ and L. In the following, we study both the UC of this algorithm (i.e. L → L * ) and the no-bloat theorem (i.e. R ′ → R ′ (P * ) when P * exists). Note that the bound V s = o(log(s)) is much stronger than the usual limit used in the method of Sieves (see Theorem 3).

Theorem 4 (No-bloat theorem). Let F 1 , . . . , F k , . . . with finite VC-dimensions

V 1 , . . . , V k , . . . . Let F = ∪ n F n . Define V (P ) = V k with k = inf{t|P ∈ F t }.
Define L V = inf P ∈FV L(P ). Consider V s = o(log(s)) and V s → ∞. Consider also that Ps minimizes L(P ) = L(P )+R(s, P ) in F s , and assume that R(s, .) ≥ 0. Assume that sup P ∈FV s R(s, P ) = o(1). Then, L( Ps ) → inf P ∈F L(P ) almost surely. Note that for well chosen family of functions, inf P ∈F L(P ) = L * . Moreover, assume that ∃P * ∈ F V * L(P * ) = L * . With R(s, P ) = R(s)R ′ (P ) and with R ′ (s) = sup P ∈FV s R ′ (P ), we get the following results:

1. Non-asymptotic no-bloat theorem: For any δ ∈]0, 1], R ′ ( Ps ) ≤ R ′ (P * )+ (1/R(s))2ǫ(s, V s , δ) with probability at least 1 -δ. This result is in particular interesting for ǫ(s, V s , δ)/R(s) → 0. 2. Almost-sure no-bloat theorem: If for some α > 0, R(s)s (1-α)/2 = O(1), then almost surely R ′ ( Ps ) → R ′ (P * ) and if R ′ (P ) has discrete values (such as the number of instructions in P or many complexity measures for programs) then for s sufficiently large, R ′ ( Ps ) = R ′ (P * );

3. Convergence rate: For any δ ∈]0, 1], with probability at least 1 -δ, L Ps ≤ inf

P ∈FV s L (P ) + R (s) R ′ (s)
=o( 1) by hypothesis +2ǫ (s, V s , δ) , where ǫ (s, V, δ) = 4-log(δ/(4s 2V )) 2s-4

.

Interpretation: Combining a code limitation and a penalization leads to UC without bloat.

5 Some Negative Results with Subsampling: Hold-Out or

Cross-Validation

When one tries to learn a relation between x and y, the "true" cost function (typically the mean squared error of a given approximate relation, which is an expectation under some usually unknown law of probability) is generally not available. It is usually replaced by its empirical mean on a finite sample. Minimizing this empirical mean is natural, but this can be done over various families of functions (e.g. trees with depth 1, 2, 3, and so on). Choosing between these various levels is hard. Typically, the empirical mean decreases as the complexity is increased, but this decrease is not generally a decrease of the generalization error, as trees of larger depth have usually a very bad generalization error due to overfitting. Therefore, the problem is somewhat multi-objective: there is a conflict between the empirical error and the complexity level. This multi-objective optimization setting has been studied in [START_REF] Bleuler | Multiobjective genetic programming: Reducing bloat using SPEA2[END_REF][START_REF] Edwin | Reducing bloat and promoting diversity using multi-objective methods[END_REF][START_REF] Ekart | Maintaining the diversity of genetic programs[END_REF]. This section is devoted to hold-out and cross-validation as tools for UC without bloat. First, let's consider hold-out for choosing the complexity level. Consider X 0 , . . . , X N , Y 0 , . . . , Y N , 2(N + 1) samples (each of them consisting in n examples, i.e. X i = (X i,1 , X i,2 , . . . , X i,n ) and Y i = (Y i,1 , Y i,2 , . . . , Y i,n )), the X i 's being learning sets, the Y i 's being (hold-out) test sets. Consider that the function can be chosen in many complexity levels, F 0 ⊂ F 1 ⊂ F 2 ⊂ F 3 ⊂ . . . , where F 0 is non-empty and

F i = F i+1 . Note Lk (f ) the error rate of the function f in the set X k of examples: Lk (f ) = 1 n n i=1 l (f, X k,i ) where l (f, x) = 1 if f fails on
x and 0 otherwise. Define f k = arg min F k Lk (.). In hold-out, after the complete learning, the resulting classifier is

f k * (n) , where k * (n) = arg min k≤N (n) l k and l k = 1 n n i=1 l (f k , Y k,i ). In the sequel, we assume that f ∈ F k ⇒ 1 -f ∈ F k and that V Cdim(F k ) → ∞ as k → ∞.
The case with hold-out leads to different cases, namely:

Greedy case: all X k 's and Y k 's are independent; this means that we test separately each complexity level F k with different learning sets X k and test sets Y k .

Case with pairing: X 0 is independent of Y 0 , ∀k, X k = X 0 and ∀k, Y k = Y 0 ; this means that we use the same learning set for all complexity levels and the same test set for all complexity levels. This case is far more usual.

Theorem 5 (No bloat avoidance with greedy hold-out). Consider greedy hold-out for choosing between complexity levels 0, 1, . . . , N (n). If N (n) is a constant, then for some distribution of examples ∀k ∈ [0, N ], P (k

* (n) = k) → 1/(N + 1). If N (n) → ∞ as n → ∞,
then for some distribution of examples such that an optimal function lies in F 0 , greedy hold-out leads to k * (n) → ∞ as n → ∞ and therefore lim sup n→∞ k * (n) = ∞.

All the following results are in the general case of N a non decreasing function of n.

Proposition 2 (Bloat cannot be controlled by hold-out with pairing, first result). Consider the case with pairing. For arbitrarily large v, there exists a distribution with optimal function in F 0 such that lim inf n→∞ Pr(k

* (n) ≥ v) > 0.
Now, let's consider a distribution that depends on n. This is interesting, as it provides lower bounds on what can be guaranteed, for a given value of n, independently of the distribution. For technical reasons, and without loss of generality with renumbering of the F k , we assume that F v+1 has a VCdimension larger than F v . We can show that, with a distribution dependent on n, lim sup n→∞ k * (n) → ∞. This leads to this other negative theorem about the control of bloat by hold-out.

Proposition 3 (Bloat can not be controlled by hold-out with pairing, second result). lim sup n k * (n) = ∞, where the distribution depends on n but is always such that an optimal function lies in F 0 . This result above is in the setting of a distribution which depends on n; it is of course not interesting for modelizing the evolution of one particular problem as the number of examples increases, but it shows that no bound on k * (n) for n ≥ n 0 can be provided, whatever may be n 0 , for hold-out with pairing, unless the distribution of problems is taken into account.

Cross-Validation for the Control of Bloat. We now turn our attention to the case of cross-validation. We formalize N -folds cross-validation as follows:

f i k = arg min F k L(., X ′ i k ), k * = arg min 1 N N i=1 L(f i k , X i k ) X ′ i k = (X 1 k , X 2 k , . . . , X i-1 k , X i+1 k , X i+2 k , . . . , X N k ) for i ≤ N
where for any i and k, X i k is a sample of n points. Greedy cross-validation could be considered as in the case of hold-out above: all X i k could be independent. This leads to the same result (for some distribution, k * (n) → ∞) with roughly the same proof. We therefore only consider crossvalidation with pairing, i.e. ∀i, k, k ′ , X i k = X i k ′ . For short, we note X i k = X i .

Theorem 6 (Negative result on subsampling). Assume that F k has a VCdimension going to ∞ as k → ∞. One can not avoid bloat with only hold-out or cross-validation, in the sense that with paired hold-out, or greedy hold-out, or cross-validation, for any V , there exists some distribution for which almost surely, k * (n) > V infinitely often whereas an optimal function lies in F 0 .

Note that propositions above show in some cases stronger forms of bloat. If we consider greedy hold-out, hold out with pairing and cross-validation with pairing, then: (1) For some well-chosen distribution of examples, greedy holdout almost surely leads to (i

) k * (n) → ∞ if N → ∞ (ii) k * (n) asymptotically uniformly distributed in [[0, N ]] if N finite,
whereas an optimal function lies in F 0 (theorem 5). ( 2) Whatever may be V = V Cdim(F v ), for some well-chosen distribution, hold-out with pairing almost surely leads to k * (n) > V infinitely often whereas an optimal function lies in F 0 (proposition 2). (3) Whatever may be V = V Cdim(F v ), for some well-chosen distribution, cross-validation with pairing almost surely leads to k * (n) > V infinitely often whereas an optimal function lies in F 0 .

Experimental Results

Some theoretical elements presented in Sections 3 and 4 are verified experimentally in this section. The experimentation are conducted using Koza-style GP [START_REF] Koza | Genetic Programming: On the Programming of Computers by Means of Natural Selection[END_REF], with a problem setup similar to the classical symbolic regression example, modified for binary classification. This is covered by theoretical results above. The GP branches used are the addition, subtraction, multiplication, protected division, and if-less-than. This last branch takes four arguments, returning the third argument if the first argument is less than the second one, otherwise returning the fourth argument. The GP terminals are the x variable, and the 0 and 1 constants. The learning task consists in minimizing the error e(i) between the desired output y i = {-1, 1} and the obtained output ŷi of the tested GP tree for the x i input, as in the following: e(i) = max(1 -y i ŷi , 0). The fitness measure used in the experiments consists in minimizing the sum of the errors to which is added a complexity factor that approximate the VC-dimension of the

GP program: f = 1 s s i=1 e(i) + k t 2 log 2 (t)
s , where t is the number of nodes of the GP program tested, s is the number of test cases used for fitness evaluation, and k is a trade-off weight in the composition of the complexity penalization relatively to the accuracy term. The s test cases are distributed uniformly in x i ∈ [0, 1], with associated y i = {-1, 1}. For x i < 0.4, each y i are equal to 1 with probability 0.25 (so y i = -1 with probability 0.75), for x i ∈ [0.4, 0.6[, y i = 1 with probability 0.5, and for x i ≥ 0.6, y i = 1 with probability 0.75. Thus, the associated classifier with best generalization capabilities would return y * i = -1 for x i < 0.4, y * i = 1 for x i ≥ 0.6 and a random output for x i ∈ [0.4, 0.6[, with a minimal generalization error of 0.3. After the evolutions, each best-of-run classifier is thus evaluated by a fine sampling of the input space, with the generalization error evaluated as the difference between the output given by the tested best-of-run classifier and the output obtained by a classifier with best generalization capabilities. Five types of GP evolutions have been tested: i) no limitation on the tree size (no depth limit and complexity trade-off k = 0), ii) depth limitation on the tree size of 17 levels (complexity trade-off k = 0), iii) soft complexity penalty in the fitness (complexity trade-off k = 0.0001), iv) medium complexity penalty in the fitness (complexity trade-off k = 0.001), and v) important complexity penalty in the fitness (complexity trade-off k = 0.01). For the three last approaches, the depth limitation of 17 levels is still maintained. The selection method used is lexicographic parsimony pressure [START_REF] Luke | Lexicographic parsimony pressure[END_REF], that is regular tournament selection 4 participants, with the smallest participant taken in case of ties. Other GP parameters are: population of 1000 individuals; evolutions on 200 generations; crossover probability of 0.8; subtree, swap and shrink mutation of probability 0.05 each; and finally half-and-half initialization with maximal depth of 5. All the experiments have been implemented using the GP facilities of the Open BEAGLE (http://beagle.gel.ulaval.ca, [START_REF] Gagné | Genericity in evolutionary computation software tools: Principles and case study[END_REF]) C++ framework for evolutionary computations. The experiments have been conducted different number of test cases varying from s = 10 to s = 100 by steps of 10. One hundred evolutions is done for each combinations of approaches tested and number of test cases, for a total of 50 000 evolutions. Figure 1 shows the average generalization errors and tree size obtained for the different approaches in function of the number of test cases used for fitness evaluation. These results show that bloat occurs when no limitation of size occurs, even when lexicographic parsimony pressure is used (see curve No limit of Figure 1b), which validates Theorem 3. Then, as stated by Theorem 2, UC is achieved using moderate complexity penalization in the fitness measure, with a convergence toward optimal generalization error of 0.3 (see curve k=0.001 of Figure 1a). Third, as predicted by Theorem 4, increasing the penalization leads to both UC and no bloat (see curve k=0.01 of both Figures 1a and1b). Note that Theorem 3 asserts that this result cannot be achieved by a priori scaling of the complexity, and that Section 5 shows that this can not be achieved by cross-validation.

Conclusion

In this paper, we have proposed a theoretical study of two important issues in Genetic Programming (GP) known as Universal Consistency (UC) and code bloat. We have shown that the understanding of the bloat phenomenon in GP could benefit from classical results from statistical learning theory. The main limit of our work is that it deals only with the statistical elements of genetic programming (effect of noise) and not with the dynamics (the effect of bounded computational power). Application of theorems from learning theory has led to two original outcomes with both positive and negative results. Firstly, results on UC of GP: there is almost sure asymptotic convergence to the optimal error rate in the context of binary classification with GP with any of the classical forms of regularizations (from learning theory): the method of Sieves, or Structural Risk Minimization. Secondly, results on code bloat: i) if the ideal target function does not have a finite description then code bloat is unavoidable (structural bloat), and ii) code bloat can be avoided by simultaneously bounding the length of the programs with some ad hoc limit and using some parsimony pressure in the fitness function (functional bloat), i.e. by combining Structural Risk Minimization and Sieves. An important point is that all methods leading to no-bloat use a regularization term; in particular, cross-validation or hold-out methods do not reach no-bloat.

Fig. 1 .

 1 Fig. 1. Generalization errors and tree sizes observed for different size limitations. Figure (a) shows the average generalization errors observed, with apparently better results for the approaches where the fitness includes some parsimony pressure. Figure (b) shows the average tree sizes obtained, where important bloat is observed for the no limitation and maximum depth limitations

Acknowledgements This work was supported in part by the PASCAL Network of Excellence, and by postdoctoral fellowships from the ERCIM (Europe) and the FQRNT (Québec) to C. Gagné. We thank Bill Langdon for very helpful comments.