Theory of Estimation-of-Distribution Algorithms

Estimation-of-distribution algorithms (EDAs) are general metaheuristics used in optimization that represent a more recent alternative to classical approaches like evolutionary algorithms. In a nutshell, EDAs typically do not directly evolve populations of search points but build probabilistic models of promising solutions by repeatedly sampling and selecting points from the underlying search space. Recently, there has been made significant progress in the theoretical understanding of EDAs. This article provides an up-to-date overview of the most commonly analyzed EDAs and the most recent theoretical results in this area. In particular, emphasis is put on the runtime analysis of simple univariate EDAs, including a description of typical benchmark functions and tools for the analysis. Along the way, open problems and directions for future research are described.


Introduction
Optimization is one of the most important fields in computer science, with many problems being NP-hard and thus not necessarily easy to solve. Hence, heuristics play a major role, i.e., optimization algorithms that try to yield solutions of good quality in a reasonable amount of time. Research over the past decades has resulted in many good heuristics developed for classical NP-hard problems. Unfortunately, these heuristics are tailored with specific problems in mind and exploit certain problem-specific properties in order to save computation time. Thus, they cannot be used for problems that do not feature these specific properties.
One alternative to problem-specific heuristics are general-purpose heuristics. The information about the problem to optimize that these algorithms have access to is fairly limited, up to the point that they are only able to compare the quality of different solutions relatively. This has the advantage that the problem itself does not have to be formalized but only the quality of a solution, as the problem formalization is communicated implicitly via the quality measure to the algorithm. In turn, this results in great reusability of the algorithm for different problems.
One class of general-purpose heuristics are evolutionary algorithms (EAs) [34]. EAs are characterized by creating new solutions from already generated solutions. Oftentimes, many solutions are stored and only changed (evolved) locally, preferably discarding bad solutions and saving good ones. Such algorithms are EAs in the classical sense [60].
The concept of EAs can be broadened when being less restrictive about what is being evolved. A similar approach to changing solutions directly is to instead change the procedure that generates the solutions in the first place. Thus, a solution-generating mechanism is evolved. Algorithms following this approach are called estimation-ofdistribution algorithms (EDAs) [29,37,54,55]. They are not EAs in the classical sense but can be considered EAs in the broad sense, as just described.
EDAs have been used very successfully in real-world applications [29,37,54,55] and recently gathered momentum in the theory community analyzing EAs [10,20,22,36,38,63,66]. The aim of theoretically analyzing EAs is to provide guarantees for the algorithms and to gain insights into their behavior in order to optimize the algorithms themselves. Common guarantees include the expected time until an algorithm finds a solution of sufficient quality, the probability to do so after a certain time, or the fact that the algorithm is even able to find desired solutions.
With this article, we provide a state-of-the-art overview about the theoretical results on EDAs for discrete domains, as that is their main field of application. To the best of our knowledge, while continuous EDAs exist, no detailed theoretical analyses have been conducted so far. We present the most commonly investigated EDAs and give an outline of the history of their analyses, providing deep insights into some of the latest results. After reading, the reader should be familiar with EDAs in general, the current state of theoretical research, and common tools used for the analyses.
In Section 2, we go more into detail about how EDAs work, we introduce the scenario used in most theoretical papers, and we provide different ways of classifying EDAs, stating the most commonly analyzed algorithms. Further, we mention some tools that are often used when deriving results for EDAs. Then, in Section 3, we give a short overview over the most commonly considered objective functions. In Section 4, we discuss the historically older results of convergence analyses on EDAs. After that, in Section 5, we present more recent results on EDAs, which consider the actual runtime of an algorithm. We end this article in Section 6 with some conclusions and open problems.

Estimation-of-Distribution Algorithms
In general, EDAs are problem-agnostic optimization algorithms that store a probabilistic model over the solution space. This model is the core part of these algorithms. It implies a probability distribution over the solution space and is iteratively refined, using samples. Ideally, the model converges to a state that only produces optimal solutions.
Since EDAs make use of sample sets -so-called populations -, they are quite similar in this respect to EAs. However, the main difference is that EAs exclusively store a population and progress using solely this information: by varying samplesso-called individuals -from the population. Thus, they have quite a local view over the solution space and advance locally. In contrast, the probabilistic model of an EDA models most of the time the entire solution space. Updates to the model are done using the old model as well as a population. Hence, EDAs employ a more general view over the solution space than classical EAs.
The probabilistic model of an EDA is used as an implicit probability distribution over the solution space, instead of an explicit distribution. This is usually done by constraining the distributions that can be modeled and by factorizing them, i.e., by writing the distribution as a product of marginal probabilities. Hauschild and Pelikan [29] distinguish between many different classes of EDAs with respect to how strongly constrained the models are. An advantage of factorizing a distribution is to save a lot of memory, since an explicit distribution would necessitate to store a probability for each solution, which is not feasible. With a factorization, only the factors have to be stored in memory. However, even then it is possible for the model to grow to sizes exponential in the input [25].
As mentioned above, EDAs also use populations, like EAs, sampled from their probabilistic model, in order to update said model. It is up to the EDA to decide what to do with its population. However, all so far theoretically analyzed EDAs have in common that they always discard their population after every iteration, valuing the model higher than the population.
In the following, we first state the optimization domain for EDAs that we consider in this chapter. Then we discuss different classifications of EDAs and name various algorithms that fall into those classes. Last, we mention the tools that are commonly used in the current theoretical research of EDAs.

Scenario
As in the theory of EAs, theoretical analyses of EDAs mainly consider pseudo-Boolean optimization, i.e., optimization of a function f : {0, 1} n → R, often referred to as fitness function. Conventionally, the function value of a bit string x is called the fitness of x.
The aspect of an EDA being a general-purpose solver is modeled as a classical black-box setting, where the algorithm only gains problem-specific information from querying the fitness function by inputting bit strings and receiving their respective fitness. In this setting, mostly two different scenarios have been of major interest: Convergence Analyses In this historically older topic, EDAs have been analyzed with respect to the convergence of their probabilistic model, i.e., if they succeed at all in optimizing certain fitness functions. We discuss this scenario in more depth in Section 4.

Runtime Analyses
A more recent trend is the analysis of an EDA's runtime on certain functions. In this scenario, the focus is the number of queries needed until an optimum or a solution of sufficient quality is sampled, i.e., the first hitting time of an algorithm sampling such a solution. Although sampling a desired solution can happen by chance, the analyses usually entail that the probabilistic model of an EDA makes it very likely for said solution to be sampled again. Section 5 goes into detail about this topic.

Classifications of EDAs
Arguably, the most straightforward way of classifying EDAs is with respect to the power of their underlying probabilistic model: univariate algorithms only use a single variable in their model per problem variable 1 . In contrast, multivariate algorithms use more than a single variable to model a problem variable. Thus, univariate EDAs are not able to capture dependencies between problem variables, whereas multivariate EDAs are explicitly constructed to do so. Pelikan et al [55] give a more fine-grained classification of EDAs, differentiating multivariate EDAs even further with respect to how many dependencies can be captured by the underlying probabilistic model.
Note that the classification into univariate and multivariate EDAs does not constrain the populations at all.

Univariate Algorithms
When optimizing a pseudo-Boolean function, univariate EDAs assume independence of all of the n different bit positions to optimize. Under this assumption, every probability distribution can be factorized into a product of n different probabilities p i , collected in a vector p of length n. A bit string x is then sampled by choosing each bit x i to be 1 with probability p i and 0 otherwise. Since each p i determines how frequently, in expectation, a 1 is sampled at position i, we call these probabilities frequencies, following the common naming convention ( [22]). The vector p is then consequently called frequency vector.
n-Bernoulli-λ -EDA Although the class of univariate EDAs does not limit the populations of the algorithms in any way, the most commonly considered univariate EDAs discard their entire population after every iteration. Thus, from a theoretical point of view, a run of a univariate EDA can be modeled as a series (p (t) ) t∈N of frequency vectors over the number of iterations t. Usually, p (0) models the uniform distribution by satisfying that p (0) i = 1/2 for each i. Friedrich et al [22] capture this class of univariate EDAs in a framework called the n-Bernoulli-λ -EDA (Algorithm 1).
The n-Bernoulli-λ -EDA samples λ individuals each iteration and performs an update to its frequency vector, using the current frequency vector as well as all of the just sampled individuals and their respective fitness. The function performing this update is called the update scheme and fully characterizes the algorithm.
Note that we do not specify a termination criterion. In fact, determining what a good criterion is may vary between different use cases of the algorithm. When considering the expected runtime of these algorithms (Section 5), we are interested in the number of fitness function evaluations until an optimal solution is sampled for the first time.
In many EDAs, if a frequency is either 0 or 1, all bits sampled at the respective position will be 0 or 1, respectively, and the update scheme will not change the frequency Algorithm 1: n-Bernoulli-λ -EDA with a given update scheme ϕ, optimizing f x ← offspring sampled with respect to p (t) ; t ← t + 1; 10 until termination criterion met; anymore. To prevent this, the algorithm is usually modified such that each frequency is only allowed to take values in an interval [m, 1 − m] ⊂ [0, 1], where m ∈ (0, 1/2] is called a margin; the values m and 1 − m are called borders. Usually, a margin of 1/n is chosen [7,10,50]. In a scenario with a margin, Line 8 of Algorithm 1 can be modified as follows: We continue to give an overview of the most commonly theoretically analyzed univariate EDAs and show how they fit into the n-Bernoulli-λ -EDA framework. We present the algorithms without a margin although they are commonly analyzed with a margin of 1/n. Since many of the following examples do not make use of the entire population of size λ (the population size) but select a certain number µ (the effective population size) of individuals according to their fitness values, we denote the kth-best individual as x (k) , where 1 ≤ k ≤ µ; ties are broken uniformly at random. Thus, x (1) denotes an individual with best fitness.
UMDA The arguably easiest update scheme is given by the Univariate Marginal Distribution Algorithm (UMDA; Algorithm 2) [49]. It samples λ individuals each iteration, of which µ best are chosen. Then, each frequency p i is set to the relative frequency of 1s at position i in the set of the µ best individuals, regardless of the current frequency.
The update scheme of UMDA allows it to go from any valid frequency to any other in a single step if not stuck. Thus, the difference of two consecutive frequencies p can only be trivially bounded by roughly 1. We call such a difference the step size of the algorithm.
PBIL A variant of UMDA that has an adjustable step size is the Population-Based Incremental Learning algorithm (PBIL; Algorithm 3) [4]. A frequency is updated in Algorithm 2: UMDA with population size λ , effective population size µ, optimizing f x ← offspring sampled with respect to p (t) ; t ← t + 1; 10 until termination criterion met; a way similar to UMDA, but the new frequency is a convex combination with parameter ρ of the current frequency and the relative frequencies of 1s at that position. Thus, the step size is now bounded by ρ, and UMDA is a special case of PBIL with ρ = 1. x ← offspring sampled with respect to p (t) ; t ← t + 1; 10 until termination criterion met; MMAS ib Another important univariate EDA is Max-Min Ant System with iterationbest update (MMAS ib ; Algorithm 4) [50], which is a special case of PBIL when setting µ = 1, i.e., when we only consider the best individual in each iteration. MMAS ib also falls into the general class of Ant Colony Optimization (ACO) algorithms [16]. Although, ACO spans an entire research topic independent of EDAs and is typically not considered to be an EDA, the process of how it produces solutions iteratively can be viewed as refining a probabilistic model. Thus, we view ACO as an EDA.
ACO considers graphs whose edges are weighted with probabilities: so-called pheromones. Additionally, the algorithm uses agents -so-called ants -that traverse  Figure 1: The ACO graph for pseudo-Boolean optimization with n = 5 bits.
the graph and thus construct paths. At each vertex v, if a path should be extended, an ant chooses an edge with a certain probability with respect to the pheromones on all of the outgoing edges of v. After the data of all ants has been collected, all pheromones decrease (they evaporate) and some are increased afterward, usually the ones that are part of the best solutions constructed. When considering pseudo-Boolean optimization, a graph for ACO can be modeled as a multigraph with n + 1 vertices from 0 to n, each vertex having exactly two outgoing edges to its direct successor (except for vertex n; see Figure 1). One of these edges is interpreted as a 0, the other one as a 1. Each solution is constructed by letting an ant traverse the graph starting at 0 and ending at n. The corresponding edges are then interpreted as a bit string of length n. Note how the probability to choose an edge corresponding to a 1 is equal to an n-Bernoulli-λ -EDA's frequency for that respective position.
MMAS ib is a variant of the Max-Min Ant System algorithm [61] that only makes an update with respect to the path of the best ant of each iteration, using a classical update rule in ACO.
Algorithm 4: MMAS ib with population size λ and evaporation factor ρ, optimizing f 1 t ← 0; x ← offspring sampled with respect to p (t) ; t ← t + 1; 10 until termination criterion met; cGA An algorithm with a different approach is the Compact Genetic Algorithm (cGA; Algorithm 5) [27]. It samples exactly two individuals each iteration and compares their bit values component-wisely. If the bits at position i are the same, the frequency p i is left unchanged. Otherwise, the frequency is adjusted by ±1/K, where K is an algorithm-specific parameter, often referred to as population size, such that the probability of sampling the bit value of the fitter individual is more likely in the next iteration.
Algorithm 5: cGA with population size K, optimizing f x ← offspring sampled with respect to p (t) ; t ← t + 1; 10 until termination criterion met;

Multivariate Algorithms
The class of multivariate EDAs consists of all algorithms that can use multiple variables to model one problem variable and thus express dependencies. A compact representation of such dependencies can be modeled as a directed graph whose vertices are the variables and whose edges denote dependencies among the variables. For each vertex, the probability distribution conditional on all its adjacent vertices with an incoming edge (its parents) is stored. This results in a factorization of the problem space that respects the given dependencies. Multivariate EDAs can assume a certain dependency model and only learn the respective (conditional) probabilities of the factorization, or they additionally try to learn a model that fits well to the samples.
The Factorized Distribution Algorithm (FDA) [48] falls into the former category. It assumes a factorization according to a so-called additively decomposable function (ADF), i.e., a function that is a sum of multivariate subfunctions. For each set of variables per subfunction, FDA creates a metavariable and expresses the objective function (the ADF) with respect to those metavariables. In each iteration, it samples solutions with respect to the factorization, selects a subset of them, and estimates the conditional probabilities based on these samples. Note that FDA is a generalization of the update of UMDA and coincides with it if no dependencies between the problem variables exist.
Another approach also using metavariables is the Extended Compact Genetic Algorithm (ECGA) [28]. Different from FDA, a metavariable of ECGA represents multiple variables at once (i.e., it is assumed that such variables are strongly correlated). Each iteration, the algorithm starts by placing each problem variable into its own class. Then, it greedily merges two classes such that a certain metric (the so-called Bayesian information criterion) is maximized, using samples from the current model. If no further improvement can be made, the merging process stops and the algorithm uses the newly created model.
The easiest of the multivariate cases is the one where each variable can be at most dependent on one other variable, i.e., a bivariate setting, and the arguably easiest probabilistic model in such a setting is a path. This model is used in the Mutual-Information-Maximization Input Clustering (MIMIC) algorithm introduced by De Bonet et al [11]. The idea of the underlying model is to construct a path that minimizes the Kullback-Leibler divergence with respect to the bivariate setting, i.e., to find a permutation that can explain the sample data best. However, since there are n! possible permutations for n variables, the authors suggest a greedy approach that makes use of the empirical entropies, i.e., the entropies of the sample data. First, a variable with minimal entropy is chosen as the start vertex of the path. Then, the path is continued by choosing a node that has minimal conditional entropy with respect to the currently last vertex in the path.
The Bivariate Marginal Distribution Algorithm (BMDA) [52] uses a somewhat similar approach. However, it does not consider paths as its model for dependency graphs but rather a forest of rooted trees. In order to determine which variables are dependent on which other variables, the Pearson's chi-square statistics is used as an indicator. If the indicator is too low, the corresponding variables are considered independent. The forest is then created greedily very similarly to regular algorithms for maximum spanning trees: iteratively, a vertex is added to one of the trees that has maximum Pearson chi-square value.
The Bayesian Optimization Algorithm (BOA) [53] is a very general multivariate EDA and constructs an arbitrary dependency graph with respect to a metric of choice. If wanted, the degree of incoming edges, i.e., the number of dependent variables, can be limited. Pelikan et al [53] propose the Bayesian Dirichlet metric as one possibility to determine the quality of a dependency graph, and they state that the general problem of finding an optimal graph is NP-hard. Thus, they suggest greedy algorithms or heuristics for efficiently creating good graphs.

Other Classifications
Another approach to classify EDAs is not to differentiate them by how many dependencies they can model but by certain invariances that their probabilistic models may have.
One such classification stems from the theory of EAs and was introduced by Lehre and Witt [39]. The authors consider a new black-box complexity known as unbiased black-box complexity in order to prove tighter lower bounds for commonly analyzed EAs. This definition is so general that it applies to any black-box algorithm optimizing pseudo-Boolean functions, thus, including EDAs.
Unbiased black-box complexity considers black-box algorithms optimizing perturbations of the hypercube, where a perturbation is any isometric automorphism of the hypercube. 3 For example, cyclicly shifting a bit string by one position to the right and changing the value of the first bit in the result is an isometric automorphism.
Given a fitness function and a perturbed variant of it, a black-box algorithm is said to be unbiased if the queries to the black box in the perturbed setting are the same as the queries of the unperturbed setting when inverted with respect to the perturbation. Thus, an unbiased algorithm does not favor certain positions over other positions or 1s over 0s, or vice versa, i.e., it has no bias in this respect.
When considering general-purpose algorithms, unbiasedness is a nice property to have, as it certifies that the algorithm has no bias with respect to the encoding of the search space. However, when considering certain problems, different values may have a strict, different meaning, such that unbiasedness with respect to those values does not make sense.
All of the EDAs presented in Section 2.2.1 are unbiased when using uniform tiebreaking.
A seemingly similar but unrelated property that many EDAs feature is that their probabilistic model does not change, in expectation, if all samples have the same fitness, i.e., there is no signal from the fitness function. Friedrich et al [22] call this property balanced, with respect to the n-Bernoulli-λ -EDA. However, this property has been considered before already by Shapiro [58], albeit with different terminology.
Although balancedness seems beneficial at first glance, it actually leads to the probabilistic model converging to one of the corners of the hypercube [22,58]. This is a general problem of martingales, i.e., random processes that do not change in expectation, with a bounded range, which will eventually end up at the bounds of their range. This means that balancedness implies a bias toward outer regions of the hypercube, also called genetic drift [3], as this is an inherent drift due to the genotypes of the sampled population. In order to overcome this bias and optimize successfully, the drift due to selection introduced by the fitness function has to be larger than the genetic drift.
Different approaches have been suggested in order to prevent an EDA's probabilistic model from quickly converging to a corner of the hypercube. Shapiro [58] proposes to reject updates done to the probabilistic model with a probability equal to the ratio of going from one model to the other. This has the advantage that the resulting implicit distribution is the uniform distribution over the hypercube. However, the transition probabilities have to be known and computed in order to get the correct rejection probabilities. Another approach proposed by Shapiro [58] and also by Friedrich et al [22] is to introduce an artificial bias that counteracts the one introduced from the balancedness.
In the context of balancedness, Friedrich et al [22] introduce another concept, which they call stable. An n-Bernoulli-λ -EDA is stable if the limit distribution of each frequency, when no fitness signal is received, is unimodal with its maximum at 1/2. This means that a stable n-Bernoulli-λ -EDA has a bias toward the center of the hypercube. The authors show that this concept is mutually exclusive with an n-Bernoulli-λ -EDA being balanced, as those have a bias toward the corners of the hypercube. The stable property is similar to the concept of an EDA's limit distribution being the uniform distribution, as considered by Shapiro [58].

Tools for Analyzing EDAs Theoretically
Most of the theoretical results on EDAs consider univariate algorithms, as we explain in Section 5. Thus, tools that make use of independent events are commonly used. However, that does not limit the use of these tools to the univariate case. Especially drift analysis, which we present later in this section, can be applied in any setting.
Many proofs make use of classical probabilistic concentration bounds, such as Markov's inequality, Chebyshev's inequality, or, most importantly, Chernoff bounds [44]. The latter is used very frequently, since the sampling process of a univariate EDA is usually done independently of the other samples. Thus, the bound can be applied.
Since the theory on EDAs usually considers first hitting times, more specialized tools suited for that purpose are used, as well. One such tool is the Coupon Collector problem [45], which gives highly concentrated first hitting time results if a certain number of events with low probability have to occur to reach the target. For EDAs, this can be thought of as a certain number of factors of the probabilistic model being at the wrong end of their spectrum, thus, slowing down optimization, since they need to be changed for the optimization process to succeed.
Another tool for determining first hitting times, and the most prominent one when looking at theory on EAs and EDAs in general, is drift theory. It is loosely akin to the potential method in complexity theory. To apply drift theory, one needs to define a potential that maps the stochastic process into the reals. Then, the expected difference of two consecutive steps of the process is considered: the drift. This can be thought of as the expected velocity of the process. If the drift can be bounded, the expected hitting time of the process reaching a target is easily deducible, i.e., if there is a known bias in the process toward a certain direction, the first hitting time can easily be bounded.
We now state the three most commonly used drift theorems. The most general theorem with respect to the prerequisites of the process -the Additive Drift Theorem (Theorem 1) -was stated by He and Yao [30]. However, the ideas used date back to Wald's equation [65].
If there is a constant δ > 0 such that, for all t < T , And if there is a δ > 0 such that, for all t < T , The Additive Drift Theorem can be applied when the expected difference of two potentials is known. However, oftentimes it is easier to determine the expected difference conditional on the current potential, i.e., E[X t − X t+1 | X t ]. Due to the law of total expectation, a lower bound for the conditional expected value is also a lower bound for the unconditional one.
A theorem more suited to processes whose potential changes at least linearly with respect to the current potential is the following Multiplicative Drift Theorem.
Theorem 2 (Multiplicative Drift [14]). Let (X t ) T ∈N be nonnegative random variables over R, each with finite expectation, and let T = min{t | X t < 1}.
If there is a constant δ > 0 such that, for all t < T , The Multiplicative Drift Theorem is not well suited if the difference in potential is dependent on the current potential but not in a linear fashion. Such cases are covered by the following Variable Drift Theorem. However, note that all these theorems assume that the difference in potential does not increase when getting closer to the goal.
If there exists a monotonically increasing function h : R ≥0 → R ≥0 such that 1/h is integrable and, for all t < T , The drift theorems above have been formulated in a simple, easy-to-read form that covers the most typical scenarios in which they are applied. However, more general drift theorems can be obtained [40,41], e. g., to apply Theorem 1 in unbounded state spaces, Theorems 2 and 3 with respect to arbitrary minimum states s min > 0 in the definition of T instead of state 1, and to allow processes adapted to arbitrary stochastic filtrations instead of the natural one implicit in the formulations above. These generalizations come partly at the expense of more complicated theorem statements, and sometimes require some additional technical assumptions on the underlying stochastic process.

Common Fitness Functions
The most commonly analyzed pseudo-Boolean functions for EDAs are ONEMAX [46] and LEADINGONES [56]. However, other functions have also been analyzed [6,7], with BINVAL being the most prominent one of them [17,48].
ONEMAX counts the number of 1s in a bit string. Thus, the unique optimum is the all-1s bit string.
This function can be generalized to a class of functions, each having a target bit string awhich denotes the unique global optimum -and yielding the number of incorrectly set bits. Note that any unbiased algorithm, as introduced in Section 2.2.3, behaves on ONEMAX exactly as on the generalized version. The ONEMAX function class is used to analyze how well an EDA performs as a hill climber. The usual expected runtime of an EDA on this function is Θ(n log n) [36,38,63,66].
Whereas ONEMAX is oftentimes considered to be the easiest pseudo-Boolean function, BINVAL is said to be the hardest [17]. In contrast to ONEMAX, where all bits are equally weighted, BINVAL uses exponentially scaled weights on its bit positions: That means that BINVAL represents value of a bit string interpreted as a binary unsigned integer.
Since the sum of all powers of 2 up to an exponent j is less than 2 j , BINVAL can be interpreted as a lexicographic order on the hypercube, where lexicographically greater bit strings have a better fitness.
As for ONEMAX, in its general form, the global optimum of BINVAL is any bit string a, and the fitness of any bit string is the weight of the respective index if the bit value is the same as the one of a, and it is 0 otherwise.
LEADINGONES yields the number of consecutive 1s in a bit string, starting from the left.
As in ONEMAX, the unique global optimum is the all-1s bit string. In its general version, the function yields the number of consecutively correctly chosen bits with respect to a fixed permutation π and a target bit string a. LEADINGONES is used to analyze how an EDA copes with dependencies between the bits. The known expected runtime of certain EDAs on this function is O(n 2 ) [10], which is compliant with the usual upper bound of EAs on this function [1].

Convergence Analyses
The earliest theoretical studies of EDAs focused mostly on their convergence, similar in style to research that had been done for evolutionary algorithms in the 1990s [56,64]. More precisely, it is studied how the algorithm behaves in the limit t → ∞, i. e., if the algorithm is allowed to run for an arbitrary amount of time. If optimal solutions will be found in this limit, the algorithm is considered effective.
Almost all convergence analyses of EDAs consider univariate models. An early publication by Höhfeld and Rudolph [32] studies the vector of frequencies p (t) in PBIL by a Markov-chain model and rigorously proves that if µ = 1 < λ and ρ > 0, it in expectation will converge to a some solution This solution need not be an optimal one but may correspond to a local optimum in which the search process is lead in the very first steps. If the fitness function f is a linear pseudo-Boolean function, then in fact E[p to the optimal solution x * . This includes classical benchmark functions like ONEMAX. However, as pointed out by Shapiro [58], convergence in expectation does not imply that PBIL eventually will sample the optimum of such functions. In fact, genetic drift may lock frequencies to values that make it impossible to sample the optimum.
PBIL was also theoretically analyzed by González et al [26] using a dynamical systems model. Convergence of the model to local optima of the fitness functions is proven for µ = 1, and it is argued that the actual PBIL will resemble the model if ρ is chosen sufficiently close to 0. Hence, the approach does not make predictions for high learning rates ρ, in particular, it excludes the special case of ρ = 1 as used in UMDA.
Several subsequent works consider UMDA and its generalization FDA. Mühlenbein and Mahnig [47] also use an approach similar to dynamical systems theory to derive a quantitative statement on the behavior of the frequencies in FDA and UMDA over time. In fact, both fitness proportionate and the usual truncation selection (take the best µ out of λ individuals) are considered. Specifically for the classical UMDA on ONEMAX, they derive that roughly where I is the so-called selection intensity that is determined from the ratio µ/λ and can be thought of as being constant. By solving a differential equation, the formula can be turned into an approximation of the expected frequency at time t. Interestingly, Equation (4) resembles a rigorous statement on the drift of the frequencies that was recently proven in [66] and is crucial for upper bounds on the runtime; see a more detailed discussion in Section 5.2.1. A more comprehensive convergence study of the FDA is given by Mühlenbein and Mahnig [48]. As a general assumption, the FDA is instantiated with the correct decomposition of an additively decomposable function , where X j ⊂ {1, . . . , n}, into its subfunctions. Then the algorithm will compute a probabilistic model, comprising unconditional and conditional frequencies from the sampled search points. Strong results are obtained if a fitness-proportionate selection scheme called Boltzmann selection is used. Under some assumptions on the initial population, the algorithm will converge to a distribution that is uniform on the set of optimal solutions. The drawback of this result is that Boltzmann selection is computationally very expensive. For the usual truncation selection, results building on simplifying assumptions are obtained. Moreover, using infinite-population models, the paper derives quantitative statements similar to Equation (4) on the time for a frequency of UMDA to converge to its optimum value, regarding ONEMAX and BINVAL.
In the 2000s, rigorous convergence proofs of the FDA (including UMDA) with fitness-proportionate [72] and truncation selection [71] followed. To study the regions of convergence and their stability, this research was supplemented by a fixed-point analysis for UMDA and the FDA with 2-tournament selection in [70]. It turns out that the FDA, given an appropriate decomposition of a non-linear function, converges under milder assumptions on the starting population than UMDA. Roughly, this indicates that a multivariate model, as used in FDA, can be superior to a univariate model, to which UMDA is restricted. However, also the analyses in [70][71][72] make the assumption of an infinite population size, which is very common in early convergence analyses of nature-inspired algorithms [64]. Infinite populations simplify the analysis since certain stochastic effects leading to a deviation from the expected behavior, so-called fluctuations such as genetic drift, vanish under this assumption. Often this type of analysis is accompanied by experiments, which support the validity of the statements also for finite population sizes. Theoretically motivated research often demands rigorous statements that also hold for finite populations, see the following sections on runtime analysis.
A more recent work by Wu and Kolonko [68] presents a convergence analysis of a so-called generalized cross-entropy optimization algorithm. The algorithm generalizes PBIL by adding so-called feasibility information to elements of the search space. This information corresponds to heuristic information used in ACO [15]. It is shown for constant ρ and under different assumptions on the feasibility information that the algorithm may stagnate in sub-optimal points due to genetic drift. However, for a time-dependent update scheme, almost sure convergence to a set of solutions that may include optimal points is proven. Finally, an initial runtime analysis on LEADINGONES is presented. However, this specific result is superseded by more detailed analyses in a follow-up paper [69] discussed below.
To conclude this overview of convergence analyses, we mention a very recent work by Ollivier et al [51]. They introduce the Information-Geometric Optimization algorithm, which is a very general EDA framework derived from three invariance properties: invariance under the parameterization of the search space, invariance under the parameterization of the probabilistic model, and invariance under monotone transformations of the fitness function. This means that IGO does not care about the encoding of the search space, the probabilistic model, or absolute fitness values. The authors show that IGO results in a general EDA that encompasses PBIL and cGA when it is used on the discrete hypercube considering Poisson binomial distributions. Further, they consider a time-continuous infinite-population version of IGO, which they call IGO flow, in the setting of linear pseudo-Boolean optimization and prove that it always converges to the optimum if the probabilistic model is not ill-initialized, i.e., none of the probabilities are initialized such that sampling the optimum is impossible.

Runtime Analyses
In contrast to convergence results as described in Section 4, the focus of runtime analyses is the number of iterations until an algorithm samples a solution of sufficient quality for the first time, usually an optimum. Normally, the analyses consider both the expected number of iterations and concentration results.
In this section, we first give an in-depth overview about the history of runtime analyses on EDAs, ending with a very detailed discussion of the most recent results. These results are summarized in Table 1. Then, we consider noisy scenarios, i. e., scenarios where the fitness function is perturbed by some kind of noise, usually as an additive term to the original fitness. In this setting, every time a solution is evaluated, the noise is drawn again and independently from any prior noise, and the goal is to optimize the underlying unperturbed function despite the noise.

Early Results
We start with a discussion of the first publications addressing runtime aspects of EDAs, which date back to the early 2000s. Although some of the runtime bounds proven in these publications can now be improved with state-of-the-art methods, the analyses already point out typical scenarios and challenges in the runtime behavior of EDAs, in particular regarding genetic drift. Also they give insights into fundamental properties of EDAs that distinguish them from other nature-inspired algorithms like EAs.

First Steps Towards Runtime Analyses
As pointed out above, rigorous runtime analyes must avoid the infinite-population model and derive statements for populations of finite sizes. However, the finiteness comes at a cost: if very small populations are used, there is a high risk of genetic drift and premature convergence in suboptimal regions of the search space. In a series of works, Shapiro [57][58][59] addresses sources of genetic drift in EDAs, quantifies its impact, and proposes measures to avoid it. In [58], he points out that the probability distribution evolved by an EDA may converge to suboptimal points and, using a dynamical systems approach, determines √ n as minimum population size for UMDA to avoid genetic drift on the ONEMAX problem and even exponential ones for NEEDLE. Later, Sudholt and Witt [63] and Krejca and Witt [36] gave rigorous proofs of the fact that genetic drift can happen up to population sizes of O( √ n log n) in cGA and UMDA. Alternatively, for PBIL, the learning rate ρ may be reduced to counteract genetic drift. Using a dynamical systems approach, Shapiro [57] derives that the learning rate should be O(1/ √ n) and O(2 −n ) to avoid genetic drift on the ONEMAX and NEEDLE function, respectively.
In [59], Shapiro also gives a rigorous theorem on the speed at which genetic drift moves the probabilistic model belonging to a specific class of EDAs called SML-EDA (including UMDA) into suboptimal regions. Also, a rigorous bound Ω(2 n/2 / √ n) is determined for the population size required to make genetic drift on NEEDLE unlikely.
Finally, in [58, p. 115], early conjectures on the runtime of UMDA appear. More precisely, the paper experimentally determines a runtime of Θ(λ √ n) for UMDA on ONEMAX (given that λ is asymptotically larger than √ n). This bound was rigorously proven in [66]. However, it should be noted that Shapiro's UMDA slightly differs from the standard.

First Runtime Analyses
The first rigorous runtime analysis of an EDA was given by Droste [17]. He considers cGA without borders and proves the general lower bound Ω(K √ n) for its expected runtime on all linear functions. Using classical drift analysis and Chernoff bounds, Droste also proves the bound O(K √ n) for ONEMAX, using K = Ω(n 1/2+ε ), i. e., slightly above the threshold stated by Shapiro [58]. This bound becomes O(n 1+ε ) for the smallest K covered by his analysis. Finally, Droste argues that BINVAL is more difficult to optimize than ONEMAX and asymptotically most difficult within the class of linear functions by proving that cGA without borders takes time O(Kn) with at least constant probability on this function if K = Ω(n 1+ε ), and expected time at least Ω(Kn). The upper bound is O(n 2+ε ) for the smallest possible K allowed. However, the lower bound Ω(Kn) does not come with a minimum value for K.
The results for ONEMAX were recently refined by Sudholt and Witt [63], using more advanced tools. In particular, all of Droste's upper bounds apply Chernoff bounds to show that genetic drift is unlikely, more precisely, he shows that the probability of a frequency dropping below 1/3 during the optimization is superpolynomially small. Using a negative drift theorem, the upper bound is improved from O(n 1+ε ) to O(n log n) in [63]. See Section 5.2.1 for more details. Regarding BINVAL, a very recent analysis by Witt [67] proves Droste's conjecture that the function is harder to optimize than ONEMAX since the expected optimization time of cGA on BINVAL is Ω(n 2 ) no matter how K is chosen. The idea of the analysis is to show for all K = o(n) that genetic drift will lock many frequencies to 0 before the optimum can be found.
The results discussed in the previous two paragraphs are summarized in the following theorem.
In the years following Droste's seminal work, runtime analysis focused more on UMDA and variants thereof. The first runtime analysis of a UMDA variant was given by Chen et al [5], who study the LEADINGONES function and a modification called TRAPLEADINGONES. An expected optimization time of O(λ n) of UMDA on LEA-DINGONES is derived under a the so-called no-random-error assumption, which is similar to an infinite-population model and basically eliminates genetic drift. The authors also show that TRAPLEADINGONES, which starts out in the same way as LEA-DINGONES but requires an almost complete change of the probabilistic model for the EDA to reach the global optimum, yields expected exponential optimization time for UMDA, using 2-tournament selection instead of the usual truncation selection. Moreover, a generalization of UMDA similar to PBIL is considered, but it turns out that the strongest bounds apply for UMDA.
Strictly speaking, Chen et al [5] only derived runtime bounds for a model of UMDA. In subsequent work [7], they therefore supplement a rigorous proof of the fact that UMDA, when using appropriate borders for the frequencies, with high probability requires superpolynomial time to optimize TRAPLEADINGONES. Similarly to Droste's early work, Chernoff bounds are applied to show that frequencies do not deviate much from their expected behavior, i. e., do not exhibit strong genetic drift. For the Chernoff bounds to be sufficiently strong, unusually large population sizes such as λ = Ω(n 2+ε ) are required.
This approach was successfully picked up and extended in a more comprehensive journal publication [8]. Using again λ = Ω(n 2+ε ), the authors show that UMDA without borders optimizes LEADINGONES in time O(λ n) with overwhelming probability. Furthermore, the utility of appropriately set frequency borders is shown on a modification called BVLO, where the fitness landscapes requires the frequency of the last bit to be changed from one extremal value to the other one. Here UMDA with borders has expected polynomial runtime, whereas UMDA without borders with overwhelming probability will be stuck at non-optimal solutions. Finally, using similar proof techniques, in particular Chernoff bounds, Chen et al [6] presented a constructed example function called SUBSTRING, on which simple EDAs and simple evolutionary algorithms behave fundamentally different. More precisely, it is proven that the (1+1) EA with any mutation probability c/n, where c > 0 is constant, with overwhelming probability needs exponential time to find the optimum of the function, while UMDA using λ = Ω(n 2+ε ) and λ /µ = O(1) with very high probability finds the optimum in time O(λ n). Specifically, it is beneficial for the optimization that UMDA can sample search points with high variances as long as all frequencies are close to 1/2. The (1+1) EA always samples with low variance in the vicinity of the best-so-far solution, which is detrimental on the specific example function.

Recent Advances
Only very few runtime analyses of EDAs were published in the years 2010-2014, most notably [50,68]. Starting from 2015, this research area gained significant momentum again, see, e. g., [10,20,21]. We now discuss the latest results in runtime analysis of EDAs. They mostly consider the standard benchmark function for EAs: ONEMAX. Using and advancing the toolbox for the analysis, matching upper and lower bounds are proven, giving a tight runtime result that allows a direct comparison of their performance with other nature-inspired algorithms.

Upper Bounds for ONEMAX
Interestingly, early runtime analyses of EDAs focused more on variants of LEADING-ONES instead of ONEMAX, which is the most commonly considered example function in evolutionary computation. In fact, the first runtime analysis of UMDA on ONEMAX was not published until 2015 [10]. A possible explanation is that the hierarchical structure of LEADINGONES makes it more accessible to a runtime analysis than ONEMAX: if the best-so-far LEADINGONES-value is k and the frequencies of the first k bits all have attained their maximum value, it is likely to sample only 1s there, which is typically needed for an improvement of the best seen function value. In contrast, there is no direct relationship between the ONEMAX-value and frequencies at specific bits. Also, modern runtime analyses of UMDA [10,38] reveal that a proof of runtime bounds for LEADINGONES can be relatively short and simple once the case of ONEMAX has been understood.
Results for cGA and MMAS ib Before we describe the advances made in the runtime analysis of UMDA in more detail, we discuss the state of the art for the simpler EDAs MMAS ib and cGA. As mentioned above, Droste [17] showed that cGA typically optimizes ONEMAX in time O(n 1+ε ), using K = n 1/2+ε . His variant of cGA does not use any borders on the frequencies, which is why he uses a comparatively large K to make convergence of a frequency to 0 by genetic drift sufficiently unlikely. More recent analyses of cGA and also other EDAs like UMDA mostly impose borders {1/n, 1 − 1/n} on the frequencies, as motivated in Section 2.2.1. Using a more careful analysis of the stochastic behavior of frequencies, the classical O(n log n) runtime can be obtained, as shown in the following summary of theorems.
Theorem 5 makes statements for two slightly different EDAs but the proofs of these statements follow roughly the same structure. Crucially, the effect of genetic drift is bounded: in the given time bound, e. g., O( √ nK) generations, the expected number of frequencies that drop below 1/3 is proven polynomially small, e. g., O(1/n 2 ). Such a statement is typically obtained from a negative drift theorem. Next, the drift of frequencies towards 1 induced by selection (the so-called bias) is analyzed. It turns out that this bias is at least proportionate to the sampling variance of the EDA: roughly each frequency p i increases by an expected amount of O p i (1 − p i )/(K ∑ n j=1 p j ) in each generation. An analysis of this variable drift, using the variable drift theorem, then gives the desired runtime bound. (As variable drift analysis was not available to Neumann et al [50], a unified and simpler proof of the statement for MMAS ib is given in [63].) In the unlikely event that a frequency has reached the wrong border 1/n due to genetic drift, an event of probability Ω(1/n) is sufficient to lift the frequency again, which is absorbed in the total runtime due to the low expected number of such bad frequencies.
1st Phase Transition Around √ n log n Theorem 5 requires K ≥ c √ n log n. Recent research reveals that cGA in fact exhibits a phase transition in the regime Θ( √ n log n), similarly for MMAS ib . If K ≤ c ′ √ n log n for a sufficiently small constant c ′ > 0, then genetic drift will outweigh the drift due to selection such that a significant number of frequencies will drop to the lower border. In this case, classical arguments on coupon collector processes show that the runtime must be at least Ω(n log n); see more arguments below in Section 5.2.2 on lower bounds. There are no upper bounds on the runtime of cGA and MMAS ib in the regime corresponding to K ≤ c ′ √ n log n, but it is conjectured that bounds resembling the existing ones for UMDA (see Theorems 6 and 7 below) can be obtained if K ∈ [c 1 log n, c 2 √ n log n] for appropriate constants c 1 , c 2 > 0.
Results for UMDA We complete this discussion of upper bounds by a review of recent advancements for UMDA. As mentioned, Dang and Lehre [10] were the first to prove upper bounds for UMDA on ONEMAX. If λ ≥ c log n for a sufficiently large constant c > 0 and λ ≥ 13eµ/(1 − c ′ ) for an arbitrarily small constant c ′ > 0, then the expected runtime of UMDA on ONEMAX is O(nλ log λ ). Hence, plugging in the smallest value of λ allowed in the statement, the bound is O(n log n log log n), i. e., slightly above the O(n log n)-bound discussed above w. r. t. cGA and MMAS ib . Dang and Lehre use a powerful proof technique to obtain their bound. Interestingly, the so-called level based theorem [9], which was originally developed for the analysis of population-based evolutionary algorithms, can be applied in this context. It is shown how the truncation selection of the best µ out of λ individuals leads to a reasonable chance of improving the best-so-far ONEMAX value and allows to satisfy the other conditions of the level-based theorem under certain parameter settings. As a side-result, using the same proof technique, also the bound O(nλ + n 2 ) w. r. t. the LEADINGONES function is obtained. Somewhat unusual, these proofs mostly consider populations instead of analyzing the values of single frequencies. For this to work, it is necessary that a frequency vector more or less unambiguously can be translated back into the population from which it was computed. This is possible in UMDA but not even in the slight generalization PBIL, where a frequency vector depends on a history of previous populations.
Obviously, the proof of the above-mentioned O(n log n log log n) bound immediately raised the question whether this was the best possible runtime of UMDA on ONE-MAX. Recently, two independent improvements of the bound were presented. The first one due to Lehre and Nguyen [38] builds on a refinement of the level-based analysis, carefully using properties of the Poisson-binomial distribution, and is summarized by the following theorem. We emphasize that UMDA always refers to the algorithm with borders 1/n and 1 − 1/n on the frequencies, i. e., Algorithm 2 extended by a step that narrows all frequencies down to the interval [1/n, 1 − 1/n]. Hence, Theorem 6 proves that the runtime of UMDA is O(n log n) for an appropriate choice of the parameters. This is tight due to the recent lower bound Ω(n log n) discussed below in Section 5.2.2. Interestingly, the set of appropriate choices for the O(n log n) behavior is confined to λ = Θ(log n), which corresponds to a parameter choice below the above-mentioned phase transition, i. e., a choice where the algorithm exhibits severe genetic drift. Also, the theorem includes a limit on µ, which is exactly in the regime of the phase transition. For greater values of µ and λ , Witt [66] independently derived runtime bounds, see the following theorem; it also includes the regime covered by Lehre and Nguyen [38], albeit with an assumption on the ratio λ /µ.

Theorem 7 ([66]).
1. Let λ = (1 + β )µ for an arbitrary constant β > 0 and let µ ≥ c √ n log n for some sufficiently large constant c > 0. Then the optimization time of UMDA, both with and without borders, on ONEMAX is bounded from above by O(λ √ n) with probability Ω(1). For UMDA with borders, also the expected optimization time is bounded in this way.
2. Let λ = (1 + β )µ for an arbitrary constant β > 0 and µ ≥ c log n for a sufficiently large constant c > 0 as well as µ = o(n). Then the expected optimization time of UMDA with borders on ONEMAX is O(λ n). For UMDA without borders, it is infinite with high probability if µ < c ′ √ n log n for a sufficiently small constant c ′ > 0.
The two statements of Theorem 7 reflect the above-mentioned phase transition. For µ ≥ c √ n log n, as demanded in the first statement, the behavior is similar to the one underlying Theorem 5 w. r. t. cGA and MMAS ib . Frequencies move smoothly towards the upper border and it is unlikely that frequencies exhibit genetic drift towards smaller values than 1/3. Hence, it is unlikely as well that UMDA without borders gets stuck with frequencies being at 0. The runtime O(n log n) is obtained for λ = c √ n log n for an appropriately large constant c > 0.
The second statement of Theorem 7 applies to a case where genetic drift is likely, but frequencies that have hit the lower border 1/n have a reasonable chance to recover in the given time span, which is O(nλ ) instead of only O( √ nλ ) now. In fact, the analysis carefully considers the drift of frequencies from the lower towards the upper border and analyzes the probability that a frequency leaves its upper border again. To do so, a very careful analysis of the bias introduced by selecting the best µ individuals is required. Without such selection, a single frequency would correspond to a so-called martingale, but due to selection there is a small drift upwards, similarly to what we described w. r. t. cGA above. Hence, the proof of Theorem 7 also gives insights into the stochastic process described by single frequencies. It is more involved than the one for cGA since UMDA can change frequencies globally instead of only by ±1/K. The runtime O(n log n) can be obtained again, this time for λ = c √ n log n. It is worth pointing out that Theorems 6 and 7 make non-overlapping statements. Theorem 7 also applies to λ above the phase transition and describes a transition of O(nλ ) to O(n √ λ ) in the runtime. However, it crucially assumes λ = (1 + Θ(1))µ in both statements, an assumption that was also useful in earlier analyses of EDAs [59] but restricts generality of the statements. In contrast, Theorem 6 applies to settings like µ = 1, λ = c log n and shows the O(n log n) bound also for this somewhat extreme choice of parameters.
We conclude this discussion of upper bounds by summarizing a recent study by Wu et al [69], who present the first runtime analysis of PBIL (called cross-entropy method (CE) in their paper). Using µ = n 1+ε log n for some constant ε > 0 and λ = ω(µ), they obtain that the runtime of PBIL on ONEMAX is O(λ n 1/2+ε/3 /ρ) with overwhelming probability. Hence, if ρ = Ω(1), including the special case ρ = 1 where PBIL collapses to UMDA, a runtime bound of O(n 3/2+(4/3)ε log n) holds, i. e., slightly above n 3/2 . In light of the detailed analyses of UMDA presented above, one may conjecture that this bound is not tight even if ρ < 1 is used, i. e., PBIL actually uses its learning approach to include solutions from several previous generations in the probabilistic model. In addition to that, a bound of the type O(n 2+ε ) on LEADINGONES is obtained if ρ = Ω(1), µ = n ε/2 and λ = Ω(n 1+ε ). Technically, Wu et al [69] use concentration bounds such as Chernoff bounds to bound the effect of genetic drift as well as anti-concentration results, in particular for the Poisson-binomial distribution, to obtain their statements. All bounds hold with high probability only since PBIL is formulated without borders. Probably, using a more detailed analysis of genetic drift and applying modern drift theorems, the bound for LEADINGONES can be improved to an expected O(n 2 ) runtime for all ρ = Ω(1), provided that the classical borders {1/n, 1 − 1/n} are used.

Lower Bounds for ONEMAX
Deriving lower bounds on the runtime of EDAs is often more challenging than deriving upper bounds. Roughly, most existing approaches show that the probabilistic model is not sufficiently adjusted towards the set of optimal solutions within a given time span. A relatively straightforward approach relates the runtime to the strength of updates of the algorithm. With respect to simple univariate algorithms like cGA and UMDA, one can show that frequencies do not increase more than by 1/K (with probability 1) resp. O(1/µ) (in expectation, assuming λ = (1 + Θ(1))µ) in a step. This naturally leads to a lower bound of Ω(K) resp. Ω(µ) on the runtime on ONEMAX. However, the bound is weak as it pessimistically assumes that each generation changes frequencies in the right direction. More detailed analyses reveal that cGA in the early phases of the optimization process only has a probability of O(1/ √ n) of performing a step where the two offspring differ in less than two bits, i. e., the probability that the outcome of a certain bit is relevant for selection is then only O(1/ √ n) [63]. Similar results can be obtained w. r. t. UMDA [36]. Thus, each bit only moves up to an expected amount of O(1/(K √ n)) resp. O(1/(µ √ n)) per generation. Then a drift analysis translates this into the lower bounds Ω(K √ n) resp. Ω(µ √ n) that appear in the following theorems. The first bound was already known for cGA without borders from Droste's work [17].
Theorem 8 (cf. [63]). The optimization time of cGA (with borders) with K ≤ poly(n) on ONEMAX is Ω(K √ n + n logn) with high probability and in expectation.

Theorem 9 ([36]
). Let λ = (1 + β )µ for some constant β > 0 and λ ≤ poly(n). Then the expected optimization time of UMDA on ONEMAX is Ω(µ √ n + n log n) (both with and without borders). [63] also state Theorem 8 in an analogous fashion for MMAS ib , with the parameter K replaced by 1/ρ. As its working principle is rather similar to cGA, we do not discuss MMAS ib further in this section.

Sudholt and Witt
The lower bounds Ω(K √ n) and Ω(µ √ n) we illustrated so far are very weak if K resp. µ are small. In fact, they can be even worse than the bounds Ω(n/logn) that follow from black-box complexity [18]. Until 2016, it was not clear whether the runtime of these simple EDAs also was bounded by Ω(n log n) or whether they could possibly optimize ONEMAX in o(n logn) time and hence be faster than simple evolutionary algorithms. A negative answer was given by the two above theorems, both of which also contain an Ω(n log n) term.
The proof of the bound Ω(n log n) is technically demanding. It relies on the following strategy: 1. Show that with high probability several frequencies, e. g., √ n many, reach the lower border before the optimum is sampled. This requires a detailed analysis of the stochastic behavior of several dependent, single frequencies instead of considering merely the sum P t : i of the frequencies, whose stochastic behavior is already quite well understood and can relatively easily be analyzed by drift analysis, as sketched in the paragraph following Theorem 5. In fact, in the detailed analysis of single frequencies, it is even required to show that some frequencies walk to the lower border while most other frequencies do not move up too far to the upper border; otherwise one cannot rule out with sufficiently high probability that the optimum is sampled in the meantime.
2. Once polynomially many frequencies have reached the lower border 1/n, a socalled coupon collector effect arises. A relatively straighforward generalization of the coupon collector theorem [44,45] to the case that still polynomially many bits have to be corrected, where a correction is made with probability at most 1/n, yields the following statement: Assume cGA reaches a situation where at least Ω(n ε ) frequencies attain the lower border 1/n. Then with high probability and in expectation, the remaining optimization time is Ω(n log n). The underlying modification of the coupon collector theorem may be called folklore in probability theory, but interesting for its own sake: collecting the last n ε coupons takes asymptotically the same time as collecting them all.
A major effort is required to flesh out the behavior sketched in Item 1 above. Roughly speaking, it is exploited that frequencies behave similarly to a martingale and can walk to the lower border due to genetic drift. However, the effect of genetic drift is dependent on many factors. When all frequencies have reached a border, genetic drift is much less pronounced than in situations where many frequencies are close to the medium value 1/2 (which is initially the case). To handle this dependency on time, it is shown that some frequencies move unusually fast, which means faster than the expected time, to the lower border while the majority of the frequencies is still at a medium value. More precisely, the proofs approximate the hitting time of the lower border by a normally distributed random variable, which is not sharply concentrated around this mean and exhibits exactly the desired reasonable probability of deviating from the mean. Additionally, the drift analysis features a novel use of potential functions that smooth out the variances of the movements of frequencies, which would be place-dependent and not applicable to the approximation by a normal distribution otherwise.
2nd Phase Transition Around log n Not much research has been done on very small values of the population size λ and K in UMDA and cGA, corresponding to very big ρ in MMAS ib . Neumann et al [50] give an exponential bound on the runtime of MMAS ib if ρ ≥ c/ log n, indicating a second phase transition in behavior around log n. Roughly speaking, if the set of possible values for a frequency becomes less than log n, then the scale is too coarse for the probabilistic model to adjust slowly towards the set of optimal solutions. For example, even after a frequency has reached its maximum 1 − 1/n once, an unlucky step may lead to a drastic decline in frequency which on average cannot be recovered in polynomial time. It is conjectured that cGA and UMDA will not optimize ONEMAX in polynomial time either if K ≤ c log n resp. λ ≤ c log n for a small constant c > 0.
Major Open Problems Even if we ignore the values below log n corresponding to the second phase transition just mentioned, the lower bounds given in Theorem 8 and 9 still do not give a complete picture of the runtime of the algorithms on ONEMAX. For example, for µ being in the medium regime between the phase transitions, i. e., µ being both ω(log n) and o( √ n log n), it is not clear whether a lower bound of the kind Ω(µn) (which would match the upper bound given above in Theorem 7) or any other runtime being ω(n log n) holds. It is an open problem to prove tight bounds on the runtime of the simple EDAs in this medium regime. As usual, we expect analyses to be harder for UMDA than for cGA, as the former algorithm can change frequencies in a global way, while the latter only changes locally by ±1/K. Some progress on the way to tight bounds has been made very recently by Lengler et al. [42], who prove a lower bound of Ω(K 1/3 n + n logn) on the expected optimization time of the cGA if K = O(n 1/2 /(log n log log n)). Hence, the expected optimization time will be Ω(n 7/6 /(log n log log n)) for K = O(n 1/2 /(log n log log n)), while it is bounded from above by O(n log n) for K = cn 1/2 if c is chosen as a sufficiently large constant. Hence, the runtime seems to depend in a multimodal way on K. Nevertheless, this still remains a conjecture since there are no upper bounds on the runtime of the cGA for K = o(n 1/2 ); there are only upper bounds for the UMDA if λ = o(n 1/2 ) that support this conjecture.
A summary of proven upper and conjectured bounds for the runtime of UMDA on ONEMAX is displayed in Figure 2. We believe that similar results hold for cGA and MMAS ib , with λ replaced by K and 1/ρ, respectively.
We have carried out experiments for UMDA on ONEMAX to gain some empirical insights into the relationship between λ and the runtime. The algorithm was implemented in the C programming language using the WELL512a random number generator. The problem size was set to n = 2000, λ was increased from 14 to 350 in steps of size 2, µ was set to λ /2, and, due to the high variance of the runs especially for small λ , an average was taken over 3000 runs for every setting of λ . The left-hand side of Figure 3 demonstrates that the runtime in fact shows a multimodal dependency on λ . Starting out from very high values, it takes a minimum at λ ≈ 20 and then increases again up to λ ≈ 70. Thereafter it falls again up to λ ≈ 280 and finally increases rather steeply for the rest of the range. The right-hand side also illustrates that the number of times the lower border is hit seems to decrease exponentially with λ . The phase transition where the behavior of frequencies turns from chaotic into stable is empirically located somewhere between 250 and 300.

New Advances in Tackling Genetic Drift
Genetic drift slows down optimization because it basically adds a random signal to the objective function. One reason this impacts the algorithms is their myopic behavior: they have to perform an update to their frequencies based only on information from the current iteration. Especially if this sample size is small, like for cGA or MMAS ib , the information gained during a single iteration may be too small to perform a sensible decision with respect to the update. In order to counteract such ill-informed updates, Doerr and Krejca [12] propose a new EDA that tries to reduce the number of incorrect frequency updates by not only relying on the information from a single iteration but from multiple previous iterations. Their significance-based cGA (sig-cGA) stores a frequency vector, like an n-Bernoulliλ -EDA, but additionally also stores a history H i for each bit position i. Each iteration, only two offspring are sampled, and the bits of the better individual are saved in the respective histories. Then the algorithm checks for each history whether a significance occurs, that is, whether the number of 1s or 0s saved is drastically more than expected when assuming that each 1 occurs with probability p i . The level of confidence can be regulated by a parameter called ε. If a significance of 1s is detected at a position, the respective frequency is set to 1 − 1/n, if a significance of 0s is detected, the frequency is set to 1/n; otherwise it is left unchanged. Overall, the algorithm only uses three different frequency values: 1 − 1/n, 1/2, and 1/n, where 1/2 is only used as a starting value -if a frequency once takes a value different from 1/2, it never returns.
This significance-based approach allows sig-cGA at the beginning of an optimization to keep frequencies at 1/2 until there is statistical proof that another value would be more beneficial. Thus, it can be thought of as an algorithm that is both balanced and stable. 4 The usefulness of this approach was shown by proving that this algorithm optimizes ONEMAX and LEADINGONES both in time O(n log n) in expectation and with high probability, which has not been proven for any other EA or EDA before [12]. Table 1: Expected run times (number of fitness evaluations) of various EDAs until they first find an optimum for the three functions ONEMAX (eq. (1)), LEADINGONES (eq. (3)), and BINVAL (eq. (2)).

Noisy Settings
In real-world optimization, the evaluation of a solution often involves a degree of uncertainty due to inaccuracies in the evaluation process. We call this uncertainty in the fitness noise. Since EDAs as general-purpose heuristics build on these inaccurate information, it is interesting to analyze how they perform when facing noise. Most EDA scenarios with noise consider ACO variants on single-destination shortestpath problems, mostly not in the context of EDAs at all [13,19,33,62]. However, some results analyze pseudo-Boolean optimization [23,24].

Combinatorial Optimization
Horoba and Sudholt [33] consider an acyclic weighted graph and are interested in finding a shortest path from each vertex to a single given destination. The noise is modeled by drawing a random nonnegative value η per edge weight w, possibly dependent on the edge, and determine its new weight w ′ by w ′ = w(1 + η). Thus, depending on the distribution of η, large weights increase more than small weights. The algorithm of interest is an ACO variant. It constructs paths from each node to the destination, using the perturbed weight and choosing an edge with a probability relative to its pheromone value with respect to the pheromones of all competing edges. The algorithm compares each constructed path with the currently best-so-far solution per node without re-evaluation. That means that the best-so-far solutions as well as their possibly perturbed weights are stored and used for look-up.
The authors provide instances on which the algorithm does not find a desired approximation within polynomial time with high probability. This is due to the best-sofar solution not being re-evaluated. Thus, if a non-optimal path is evaluated to be very good by chance, it will get reinforced many times, making it more unlikely to sample other paths that, additionally, have to be evaluated even better. However, the authors also prove that optimization will succeed if the noise follows the same distribution for every edge.
Sudholt and Thyssen [62] extend the results from Horoba and Sudholt [33] by considering a vaster range of noise distributions, showing how long it takes to approximate optimal solutions or even when optimization succeeds.
Doerr et al [13] consider a similar scenario to the one analyzed by Horoba and Sudholt [33], the difference being that the weights of the graphs are purely random, i.e., there is no ground-truth to rely on. This setting makes it harder to define what an optimal solution actually is.
The authors first consider a multigraph consisting of two nodes with multiple edges between those nodes. They call an edge preferred if its probability of being shorter than any other edge from the same vertex is at least 1/2 + δ , where δ > 0 is a constant, and they state how this scenario relates to armed-bandit settings. Using the same ACO algorithm as Horoba and Sudholt [33] but re-evaluating the best-so-far solution each iteration, the authors give an upper bound on the expected time until the pheromone on the preferred edge is maximal. They then provide examples of weight distributions that result in an edge being preferred. The paper concludes with a more general graph setting that assumes that there exists an inductively defined set of edges S, starting at a given node, such that each edge extending paths using edges from S is preferred. If S is a tree, the authors give an upper bound for the expected time until the considered ACO variant maximizes the pheromones on all of the edges in S.
Feldmann and Kötzing [19] analyze the same setting as Doerr et al [13] but investigate another ACO variant: MMAS-fp. This algorithm does not store best-so-far solutions but always makes an update with respect to the current samples; however, the update is done with respect to each sample's fitness. Thus, good solutions yield larger changes in the update than bad solutions. The authors explain the difference in this approach with respect to the works of Horoba and Sudholt [33] and Doerr et al [13] that MMAS-fp optimizes paths that are shortest in expectation. They prove this claim by providing upper bounds on the expected number of iterations until MMAS-fp finds expected shortest paths in graphs where, for each node, the difference between the expected lengths of different outgoing edges can be lower-bounded by a value δ > 0, which influences the runtime.

Pseudo-Boolean Optimization
Friedrich et al [23] (conference version [21]) also consider MMAS-fp, just like Feldmann and Kötzing [19], but in the setting of optimizing linear pseudo-Boolean functions. The noise is mostly modeled as Gaussian additive posterior noise, i.e., when evaluating the fitness of an individual, a normally distributed random variable is added to the fitness, every time anew and independently. The authors show that MMAS-fp does, what they call, scale gracefully in this scenario. That means that for every polynomially bounded variance of the noise, there is a configuration of MMAS-fp such that the runtime is polynomially bounded as well. Since the runtime results hold with high probability, by performing an uninformed binary search using restarts, the correct variance of a problem with Gaussian noise can be guessed correctly within polynomial time. Thus, MMAS-fp can be modified such that the runtime is, with high probability, polynomial if the variance of the noise is.
Additionally, the authors extend their results to other posterior noise than Gaussian. Further, they consider a prior noise model where, before evaluating the fitness of an individual, a uniformly randomly chosen bit is flipped. In both of these settings, they prove that the algorithm scales gracefully.
Friedrich et al [24] (conference version [20]) also consider cGA under the Gaussian additive posterior noise model. As for MMAS-fp, they prove that the algorithm scales gracefully. Further, they show that the (µ + 1) EA, a commonly analyzed EA, does not scale gracefully. Both results of Friedrich et al [23,24] suggest that EDAs are inherently more tolerant to noise than standard EAs, as the EDAs did not need to be modified to cope with noise, except for choosing correct parameters. The authors also compare the restart version of cGA with an approach that uses resampling in order to basically remove the noise in the fitness, as described by Akimoto et al [2]. Since the number of resamples is closely tied to the noise's variance, the cGA variant using restarts instead of resampling emerges victorious.

Conclusions and Open Problems
We have given an overview of the state of the art in the theory of discrete EDAs, where the most recent research surpasses convergence analyses and instead deals with the runtime of especially simple univariate EDAs like cGA, UMDA, and PBIL. In this domain, increasingly precise results have been obtained with respect to well-established benchmark problems like ONEMAX, but as we have emphasized in this article, there are several open problems even for this simple problem. In particular, a complete picture of the runtime of the simple EDAs depending on their parameters is still missing. We think that further results for benchmark functions will give insight into the right choice of specific EDAs, including the choice of parameters such as the population size and the borders on the frequencies depending on problem characteristics. We also expect that this research will lead to runtime results and advice on the choice of algorithms and parameters w. r. t. more practically relevant combinatorial optimization problems. Here in particular, noisy settings or, more generally, optimization under uncertainty seem to represent scenarios where EDAs can outperform classical evolutionary algorithms. Also, the combinatorial structure may favor the application of multivariate EDAs, a type of EDAs for which almost no theoretical results exist yet.