--- title: "A General Approach to Running Time Analysis of Multi-objective Evolutionary Algorithms" authors: - "Chao Bian" - "Chao Qian" - "Ke Tang" doi: "10.24963/ijcai.2018/195" url: "https://ijcai.org/proceedings/2018/195" pdf: "https://ijcai.org/proceedings/2018/0195.pdf" venue: "Proceedings of the Twenty-Seventh International Joint Conference on Artificial Intelligence" year: 2018 pages: "1405-1411" track: "Main track" keywords: - "Heuristic Search and Game Playing: Combinatorial Search and Optimisation" - "Heuristic Search and Game Playing: Heuristic Search" --- # A General Approach to Running Time Analysis of Multi-objective Evolutionary Algorithms Chao Bian, Chao Qian, Ke Tang *Proceedings of the Twenty-Seventh International Joint Conference on Artificial Intelligence*, pages 1405-1411. DOI: [10.24963/ijcai.2018/195](https://doi.org/10.24963/ijcai.2018/195) ## Abstract Evolutionary algorithms (EAs) have been widely applied to solve multi-objective optimization problems. In contrast to great practical successes, their theoretical foundations are much less developed, even for the essential theoretical aspect, i.e., running time analysis. In this paper, we propose a general approach to estimating upper bounds on the expected running time of multi-objective EAs (MOEAs), and then apply it to diverse situations, including bi-objective and many-objective optimization as well as exact and approximate analysis. For some known asymptotic bounds, our analysis not only provides their leading constants, but also improves them asymptotically. Moreover, our results provide some theoretical justification for the good empirical performance of MOEAs in solving multi-objective combinatorial problems. ## Full text ## 1 Introduction Multi-objective optimization, which requires optimizing several objective functions simultaneously, arises in many areas. Since the objective functions are usually conflicting, the goal is to find a set of Pareto optimal solutions (or the Pareto front), which represent different optimal trade-offs between objectives. Evolutionary algorithms (EAs) [B¨ack, 1996] are a kind of randomized heuristic optimization algorithms, inspired by natural evolution. They maintain a set of solutions (called a population), and iteratively improve the population by genetic operators. Due to their population-based nature, EAs are popular for solving multi-objective optimization problems, and have been found well useful in many real-world applications [Coello Coello and Lamont, 2004]. However, due to their sophisticated behaviors of mimicking natural phenomena, the theoretical analysis of EAs is difficult. Much effort thus has been devoted to understanding the behavior of EAs from a theoretical point of view [Neumann and Witt, 2010; Auger and Doerr, 2011], but most of them focus on single-objective optimization. In fact, multi-objective EAs (MOEAs) are even more difficult to be analyzed owing to the hardness of multi-objective optimization. To the best of our knowledge, only a few pieces of casespecific studies have been reported on the running time analysis of MOEAs. The running time complexity, which measures how fast an algorithm solves an optimization problem, is one essential theoretical aspect. In [Giel, 2003], the GSEMO (a simple MOEA with a global mutation operator) was proved to find the Pareto front of LOTZ in O(n3 ) expected running time (where n is the problem size), while for another bi-objective problem COCZ, it needs O(n2 log n) expected time [Qian et al., 2013]. For mLOTZ and mCOCZ with m objectives, which are generalized from LOTZ and COCZ, both the expected running time of the SEMO (a counterpart of the GSEMO, but with a local mutation operator) were shown to be O(nm+1 ) [Laumanns et al., 2004]. More results on synthetic problems include [Friedrich et al., 2010; Giel and Lehre, 2010; Friedrich et al., 2011; Neumann, 2012; Doerr et al., 2013; 2016; Qian et al., 2016; Osuna et al., 2017]. The analysis on NP-hard multi-objective combinatorial problems has been only slightly touched. For biobjective minimum spanning trees, the GSEMO was proved able to find a 2-approximation of the Pareto front in expected pseudo-polynomial time [Neumann, 2007]; and for multiobjective shortest paths, a variant of the GSEMO can achieve an (1 + ϵ)-approximation in polynomial time [Horoba, 2009; Neumann and Theile, 2010]. Note that the analysis approaches employed in most of the above mentioned studies are case-specific, which cannot provide a general guidance to analyze the running time of a given MOEA solving a given problem. Meanwhile, ad hoc analyses starting from scratch are quite difficult. In this paper, we thus propose a general approach (Theorem 1) for estimating upper bounds on the running time of MOEAs. The idea is to contrast a given sophisticated MOEA process with an easilyanalyzable process, instead of directly analyzing the MOEA process. We apply the approach to diverse situations, including bi-objective and many-objective optimization as well as exact and approximate analysis. The theoretical results are: - For the GSEMO solving the two bi-objective problems LOTZ and COCZ, the expected time for finding the Pareto front is at most 6n3 and 3n2 logn (Theorems 2, 3), respectively, which are consistent with the known asymptotic bounds O(n3 ) [Giel, 2003] and O(n2 logn) [Qian et al., 2013], and further give the leading constants. - For the SEMO solving the scalable m-objective problem mCOCZ where m ≥ 4, the expected time for finding the Pareto front is O(nm) (Theorem 4), which is tighter than the known bound O(nm+1 ) [Laumanns et al., 2004] by a factor of n. Furthermore, a better upper bound O(n3 log n) is also derived for m = 4. - For optimizing two linear functions simultaneously, which often appears in multi-objective combinatorial problems such as multi-objective minimum spanning trees and multiobjective knapsacks, a variant of the GSEMO can find a good approximation of the Pareto front in polynomial time (Theorem 5). Our analysis thus provides some theoretical justification for the good empirical performance of MOEAs in solving multi-objective combinatorial problems [Zhou and Gen, 1999; Ishibuchi et al., 2015]. The rest of this paper starts with a section of preliminaries. Section 3 presents the proposed approach, which is applied in three subsequent sections to analyze MOEAs in different situations. Section 7 concludes the paper. ## 2 Preliminaries ## 2.1 Multi-objective Evolutionary Algorithms Multi-objective optimization requires simultaneously optimizing two or more objective functions, as shown in Definition 1. We consider maximization here, while minimization can be defined similarly. The objectives are usually conflicting, and thus there is no canonical complete order in the solution space S. The comparison between solutions relies on the domination relationship, as presented in Definition 2. A solution is Pareto optimal if there is no other solution in S that dominates it. The set of objective vectors of all the Pareto optimal solutions constitutes the Pareto front. The goal of multi-objective optimization is to find the Pareto front, that is, to find at least one corresponding solution for each objective vector in the Pareto front. Definition 1 (Multi-objective Optimization). Given a feasible solution space S and objective functions f1, f2, . . . , fm, multi-objective optimization can be formulated as Definition 2 (Domination). Let f = (f1, f2, . . . , fm) : S → Rm be the objective vector. For two solutions s and s′ ∈ S: 1. s weakly dominates s′ if ∀1 ≤ i ≤ m, fi(s) ≥ fi(s′ ), denoted as s ⪰ s′ ; 2. s dominates s′ if s ⪰ s′ and fi(s) > fi(s′ ) for some i, denoted as s ≻ s′ . The GSEMO algorithm is a simple MOEA for multiobjective optimization over the Boolean solution space S = {0, 1}n . As described in Algorithm 1, it randomly selects an initial solution, and then repeatedly tries to improve the population P. In each iteration, a solution uniformly selected from the current P is used to generate a new solution by bitwise mutation; then the newly generated solution is compared with the solutions in P, and only non-dominated solutions ## Algorithm 1 GSEMO Algorithm Given the solution space S = {0, 1}n and the objective vector f, the procedure: - 1: Choose s ∈ S uniformly at random. - 2: P ← {s}. - 3: repeat until some criterion is met - 4: Choose s from P uniformly at random. - 5: Create s′ by flipping each bit of s with prob. 1/n. - 6: if ∄z ∈ P such that z ≻ s′ - 7: P ← (P \ {z ∈ P | s′ ⪰ z}) ∪ {s′ }. are kept. The SEMO algorithm is the same as the GSEMO except that bit-wise mutation which searches globally is replaced by one-bit mutation which searches locally, i.e., line 5 of Algorithm 1 becomes "Create s′ by flipping a randomly chosen bit of s". These two algorithms explain the common structure of various MOEAs and are widely used in theoretical analyses [Laumanns et al., 2004; Friedrich et al., 2010; Qian et al., 2013]. We will also use them in case studies. The running time of a MOEA is usually measured by the number of calls to f (the most costly computational process) until finding the Pareto front or an approximation of the Pareto front. For any c ≤ 1, a set P of solutions is a c-approximation of the Pareto front if for each objective vector (f ∗ 1 , f ∗ 2 , . . . , f ∗m) in the Pareto front, there always exists a solution s ∈ P such that ∀i ∈ {1, 2, . . . ,m} : fi(s) ≥ c · f ∗ i . ## 2.2 Markov ChainModeling EAs often generate offspring solutions from their current solutions rather than the historical ones; thus, they can be modeled as Markov chains. Let X be the population space and X ∗ ⊆X be the target population space. If the goal is to find a c-approximation of the Pareto front, each population in X ∗ is a c-approximation of the Pareto front. Let ξt ∈ X be the population after t generations. Then, a MOEA can be described as a random sequence {ξ0, ξ1, ξ2, . . .}. Since ξi+1 can often be decided from ξi (i.e., P(ξi+1|ξi, . . . , ξ0) = P(ξi+1|ξi)), the random sequence forms a Markov chain {ξt}+∞ t=0 with state space X , denoted as "ξ ∈ X" for simplicity. Given a Markov chain ξ ∈ X and ξt0 = x, we define τ as a random variable such that τ = min{t ≥ 0 | ξt0+t ∈ X ∗ }. That is, τ is the number of steps needed to reach the target space for the first time. The mathematical expectation of τ, E[[τ|ξt0 = x]] = ∑∞i=0 i · P(τ = i|ξt0 = x), is called the conditional first hitting time (CFHT). If ξt0 is drawn from a distribution πt0 , the expectation of the CFHT over πt0 , E[[τ|ξt0 ∼ πt0 ]] = ∑x∈X πt0 (x)E[[τ|ξt0 = x]], is called the distribution-CFHT (DCFHT). For a Markov chain ξ ∈ X modeling a MOEA, E[[τ|ξ0 ∼ π0]] is just the expected number of iterations of the MOEA until reaching X ∗ . We always use E[[·]] to denote the expectation of a random variable. A Markov chain ξ ∈ X is absorbing, if ∀t ≥ 0 : P(ξt+1 ∈ X ∗ |ξt ∈ X ∗ ) = 1. Note that a Markov chain can be transformed to be absorbing by making it unchanged once finding a target state, which obviously does not affect its first hitting time. Lemma 1 is to compare the DCFHT of two absorbing Markov chains, which will be used in this paper. Lemma 1. [Yu et al., 2015] Given two absorbing Markov chains ξ ∈ X and ξ ′ ∈ Y with target spaces X ∗ and Y ∗ , respectively, let τ and τ ′ denote their hitting times, respectively, and let πt denote the distribution of ξt. Given a series of values {ρt ∈ R}+∞ t=0 with ρ = ∑+∞ t=0 ρt and a mapping φ : X → Y with ∀x ∈ X \ X ∗ : φ(x) /∈ Y ∗ , if for all t ≥ 0, where φ-1 (y) = {x ∈ X | φ(x) = y}, πφ t (u) = ∑x∈φ-1 (u) πt(x), then E[[τ |ξ0 ∼ π0]] ≤ E[[τ ′ |ξ ′0 ∼ πφ 0 ]] + ρ. ## 3 The Proposed Approach For a MOEA solving a multi-objective problem, we propose a general approach (i.e., Theorem 1) to analyze running time upper bounds for finding a c-approximation of the Pareto front. Note that c = 1 implies finding the Pareto front. The idea is to model the given MOEA process as a Markov chain ξ ∈ X , and compare it with an easy-to-analyze chain ξ ′ ∈ Y. From the condition Eq. (1) of Theorem 1, we can see that the long-term behavior of the given chain ξ ∈ X is waived, since E[[τ|ξt]] is not involved. In the comparison, we need a function hα,c to measure the goodness of x ∈ X. As shown in Definition 3, hα,c(x) is the sum of two terms: the maximal weighted sum of objective values of the solution in x, and the number of objective vectors in F ∗ which are c-approximated by x. It is easy to verify that g(F ∗ , x, c) ≤ |F ∗ |, and for any x which is a c-approximation of F ∗ , g(F ∗ , x, c) = |F ∗ |. When α and c are clear, we will write h for short. Definition 3. Given a Markov chain ξ ∈ X modeling a MOEA solving a multi-objective problem (f1, f2, . . . , fm) with the Pareto front F ∗ , let the target space X ∗ = {x ∈ X | x is a c-approximation of F ∗ }, where c ≤ 1. For a nonnegative real vector α = (α0, α1, ..., αm), we define a function hα,c : X → R as for any x ∈ X , where f0(s) = 1, g(F ∗ , x, c) = |{(f ∗ 1 , f ∗ 2 , . . . , f ∗m) ∈ F ∗ | ∃s ∈ x, ∀1 ≤ i ≤ m : fi(s) ≥ c · f ∗ i }| and | · | denotes the cardinality of a set. If ∀x ∈ X , hα,c(x) ∈ N0 and then hα,c is called a well-defined function. Theorem 1. Given a Markov chain ξ ∈X modeling a MOEA solving a multi-objective problem with the Pareto front F ∗ , let the target space X ∗ ={x∈X|x is a c-approximation of F ∗ }, where c ≤ 1. If there exists a well-defined function hα,c : X → N0 and a Markov chain ξ ′ ∈ Y = {0, 1}r with Y ∗ = {1r } satisfying that ∀x /∈ X ∗ , ∀t ≥ 0, where r = min{h(x)|x ∈ X ∗ }, [r] denotes {0, . . . , r} and δ<1, the expected number of iterations until the MOEA finding a c-approximation of F ∗ when starting from x0 is at most Proof. We use Lemma 1 to prove it by comparing the two Markov chains ξ ∈ X and ξ ′ ∈ Y. We construct the mapping function φ : X → Y as ∀x ∈ X, φ(x) = 1min{h(x),r} 0r-min{h(x),r} . Since h is well-defined, we have x ∈ X ∗ iff φ(x) ∈ Y ∗ = {1r }. Let πt denote the distribution of ξt. Then, for any t, where the first equality is by the absorbing property of ξ ∈ X and the definition of φ, the inequality is by Eq. (1), and the last equality is by the absorbing property of ξ ′ ∈ Y and the definition of πφ t in Lemma 1. Thus, the condition of Lemma 1 holds with ρt = δ · (1 - πt(X ∗ )). When the MOEA starts from x0, i.e., π0(x0) = 1, we get from Lemma 1 that By applying ∑+∞ t=0 (1 - πt(X ∗ )) = E[[τ|ξ0 = x0]] and φ(x0) = 1min{h(x0),r} 0r-min{h(x0),r} to the above inequality, the theorem holds. Using the above approach to derive running time upper bounds of a MOEA (modeled by a chain ξ ∈ X) for finding a c-approximation of the Pareto front, one needs to 1. design a well-defined function hα,c : X →N0 and construct a reference chain ξ ′ ∈Y={0, 1}r with Y ∗ ={1r }; 2. analyze the transition probabilities (i.e., P(ξt+1|ξt) and P(ξ ′t+1 |ξ ′t )) and the CFHT of the chain ξ ′ ∈Y (i.e., E[[τ ′ |ξ ′t ]]), and then examine Eq. (1) to get the value of δ; 3. divide E[[τ ′ |ξ ′0 ]] by 1-δ to get an upper bound on E[[τ |ξ0]]. Since the running time of a MOEA is measured by the number of fitness evaluations, the cost of each iteration should be considered. For example, the expected running time of the GSEMO (i.e., Algorithm 1) is just the expected number of iterations (i.e, E[[τ|ξ0]]), since we only need to evaluate the newly generated solution s′ in each iteration. ## 4 Bi-objective Analysis In this section, we apply the proposed approach to analyze the running time of the GSEMO solving two biobjective pseudo-Boolean problems LOTZ and COCZ, which are widely used in MOEAs' theoretical analyses. The derived running time bounds in Theorems 2 and 3 are consistent with the known asymptotic bounds O(n3 ) [Giel, 2003] and O(n2 log n) [Qian et al., 2013], and further give the leading constants. For LOTZ as presented in Definition 4, it is to maximize the number of leading 1-bits and trailing 0-bits of a solution s ∈ {0, 1}n simultaneously. The Pareto front is F ∗ = {(i, n - i) | 0 ≤ i ≤ n}, and the set of all the Pareto optimal solutions is S ∗ = {0n , 10n-1 , . . . , 1n-1 0, 1n }. Definition 4 (LOTZ [Giel, 2003]). The pseudo-Boolean function LOTZ: {0, 1}n → N20 is defined as where sj ∈ {0, 1} is the j-th bit of s. Theorem 2. For the GSEMO solving LOTZ, the expected running time for finding the Pareto front is at most 6n3 . Proof. We use Theorem 1 to prove it. Let ξ ∈X model the GSEMO solving LOTZ. From the algorithm procedure, we know that the solutions in any x ∈ X (i.e., in any possible population P) are incomparable. Let c=1 (i.e., the goal is to find the Pareto front F ∗ ). Then X ∗ = {S ∗ }. We design the function h by setting α=(0, 1, 1), i.e., Note that g(F ∗ , x, 1) in Definition 3 equals to |x ∩ S ∗ | here. For any x ∈X \ X ∗ , maxs∈x(f1(s)+f2(s)) ≤ n and |x ∩ S ∗ | < n+1, thus h(x) < 2n+1 = h(S ∗ ). It is easy to see that h(x) ∈N0. Thus, h is well-defined and r = 2n+1. For the Markov chain ξ ′ ∈ Y = {0, 1}r with Y ∗ = {1r }, we construct its transition behavior as follows: for any ξ ′t /∈ Y ∗ , ξ ′t+1 is generated by flipping the first 0-bit of ξ ′t with probability 1/r; otherwise, ξ ′t+1 = ξ ′t . Thus, for any ξ ′t ∈Y, E[[τ ′ |ξ ′t ]] = r · |ξ ′t |0, where | · |0 denotes the number of 0-bits. Then, we are to investigate Eq. (1). For any ξt = x /∈ X ∗ , assume that h(x) = k < r and let ˆ s ∈ arg maxs∈x(f1(s) + f2(s)). It is easy to verify that h(ξt+1) ≥ h(ξt) = k, since |x∩S ∗ | never decreases and if ˆ s is deleted, the newly included solution s′ must weakly dominate ˆ s (line 7 of Algorithm 1) and thus f1(s′ )+f2(s′ ) ≥ f1(ˆ s)+f2(ˆ s). We then show that h(ξt+1) ≥ h(ξt) + 1 = k + 1 with probability at least 1 en2 by considering two cases. (1) If |x ∩ S ∗ | = 0, ˆ s must not be Pareto optimal, and flipping only its first 0-bit will generate a new solution s′ with f1(s′ )+f2(s′ ) > f1(ˆ s)+f2(ˆ s) = k. (2) If |x∩S ∗ | > 0, flipping only the last 1-bit or the first 0-bit of a specific Pareto optimal solution in x ∩ S ∗ can generate a new Pareto optimal solution, and thus |ξt+1 ∩S ∗ | = |x∩S ∗ | +1. In both cases, h(ξt+1) ≥ h(x) + 1 = k + 1. Note that the probability of selecting a specific solution from x in line 4 of Algorithm 1 is 1 |x| and the probability of flipping only one specific bit in line 5 is 1n (1- 1n )n-1 ≥ 1 en . Since the solutions in x are incomparable, we have |x| ≤ n. Thus, our claim that P(h(ξt+1) ≥ k +1|ξt = x) ≥ 1 en2 holds. We then have where the first inequality is by min{h(ξt+1), r} = h(ξt+1), h(ξt+1) ≥ k and E[[τ ′ |ξ ′0 = 1i 0r-i ]] = r(r - i) decreasing with i, and the second is by P(h(ξt+1) ≥ k + 1|ξt = x) ≥ 1 en2 . Meanwhile, by the transition behavior of ξ ′ ∈ Y, we get By combining the above two formulas, the condition Eq. (1) of Theorem 1 holds with δ = 1- r en2 < 1- 2 en . Thus, we get that starting from any x0 ∈ X, the expected number of iterations until finding the Pareto front is at most E[[τ ′ |ξ ′0 = 1h(x0) 0r-h(x0) ]]/(1 - δ) ≤ E[[τ ′ |ξ ′0 = 0r ]] · en 2 = (2n+1)2 · en 2 ≤6n3 . Since only a new solution is evaluated in each iteration of the GSEMO, the theorem holds. For COCZ as presented in Definition 5, the first objective is to maximize the number of 1-bits, and the other is to maximize the number of 1-bits in the first half of a solution plus the number of 0-bits in the second half. The Pareto front is {(n/2+i, n-i) | 0 ≤ i ≤ n/2}, and the set of all the Pareto optimal solutions is S ∗ = {s ∈ {0, 1}n | ∑n/2 i=1 si = n/2}. Definition 5 (COCZ [Qian et al., 2013]). The pseudo-Boolean function COCZ: {0, 1}n → N20 is defined as where n is even. Theorem 3. For the GSEMO solving COCZ, the expected running time for finding the Pareto front is at most 3n2 log n. Proof. Let ξ ∈X model the GSEMO solving COCZ and c=1. Then, X ∗ ={{s0 , . . . , sn/2 }|si ∈S ∗ , ∑nj= n2 +1 sij = i}. We design the function h by setting α=(-n/4, 1/2, 1/2), i.e., It is easy to verify that h is well-defined and r=n+1. For the chain ξ ′ ∈Y ={0, 1}r with Y ∗ ={1r }, the transition behavior is constructed as follows. Suppose ξ ′t =s and ∑n/2 i=1 si = j. ξ ′t+1 is generated by flipping the first 0-bit of s with probability n/2-j r for j < n2 and |s|0 r for j = n2 ; otherwise, ξ ′t+1 = ξ ′t . Thus, E[[τ ′ |ξ ′0 = 1k 0r-k ]] - E[[τ ′ |ξ ′0 = 1k+1 0r-k-1 ]] = r n/2-k for k < n2 , and r n+1-k for n2 ≤ k < r. We then show that Eq. (1) holds with δ ≤1- 2 e(n+1) . For any ξt = x /∈ X ∗ , assume that h(x) = k < r and let ˆ s ∈ argmaxs∈x(f1(s)+f2(s)). It is easy to verify that h(ξt+1)≥ h(ξt)=k. We consider two cases. (1) If |x ∩ S ∗ |=0, ˆ s must not be Pareto optimal and k < n/2. Flipping only one 0bit in the first half of ˆ s will generate a new solution s′ with f1(s′ )+f2(s′ )=f1(ˆ s)+f2(ˆ s)+2. Thus, P(h(ξt+1)≥k+1)≥ n/2-k en|x| . By using the same analysis as Eqs. (2) and (3), we get where the inequality is by r = n + 1 and |x| ≤ n/2 + 1. (2) If |x ∩ S ∗ | > 0, we have k = n/2 + |x ∩ S ∗ | ≥ n/2 + 1. Let ˆ C ∈ argmaxC⊆x∩S ∗ |C| s.t. maxs∈C f1(s) - mins∈C f1(s) = |C| - 1. Let u = argmaxs∈ ˆ C f1(s) and v = arg mins∈ ˆ C f1(s). Then, flipping one 0-bit of u or flipping one 1-bit in the second half of v can generate a new Pareto optimal solution. We thus have P(h(ξt+1) = k+1) ≥ n-|u|1+|v|1-n/2 en|x| ≥ n/2+1-|x∩S ∗ | en|x| = n+1-k en|x| . Using the same analysis as case (1), we also get δ ≤ 1- 2 e(n+1) . By applying δ ≤ 1 - 2 e(n+1) and E[[τ ′ |ξ ′0 = 0r ]] ≤ 2r(log(n/2 + 1) + 1) to Theorem 1, the theorem holds. ## 5 Many-objective Analysis In this section, we apply the proposed approach to analyze the running time of the SEMO solving the scalable m-objective problem mCOCZ, where m ≥ 4. Note that as the number of objectives increases, the size of the Pareto front can exponentially increase and thus the problem becomes more complicated. Our derived running time bound in Theorem 4 is asymptotically tighter than the known bound O(nm+1 ) [Laumanns et al., 2004]. As presented in Definition 6, all the m objectives are cooperative in the first half of a solution (i.e., maximizing the number of 1-bits), and the second half is divided into m′ = m/2 blocks, each of which is to maximize the number of 1-bits and 0-bits simultaneously. The Pareto front is F ∗ = {( n2 +i1, n2 +n′-i1, . . . , n2 +im′ , n2 +n′-im′) | ∀1 ≤ j ≤ m′ : 0 ≤ ij ≤ n′ }, where n′ = n/m is the size of each block. The set of all the Pareto optimal solutions is S ∗ = {s ∈ {0, 1}n | ∑n/2 i=1 si = n/2}. It is easy to see that |F ∗ | = (n′ + 1)m′ and |S ∗ | = 2n/2 . Definition 6 (mCOCZ [Laumanns et al., 2004]). The pseudoBoolean function mCOCZ: {0, 1}n →Nm0 is defined as where m= 2 · m′ , n =m · n′ and m′ , n′ ∈ N0. Theorem 4. For the SEMO solving mCOCZ, the expected running time for finding the Pareto front is O(nm) for m > 4 and O(n3 log n) for m = 4. Proof. According to the intermediate result in Theorem 11 of [Laumanns et al., 2004], i.e., the number of mutations allocated to non-Pareto-optimal solutions is O(nm/2+1 logn), we only need to consider the number of mutations allocated to Pareto optimal solutions, and thus we can assume that after O(nm/2+1 logn) iterations in expectation, the population will always contain only Pareto optimal solutions. We then use Theorem 1 to analyze the expected number of iterations (denoted by E[T ]) until finding the Pareto front when starting from any set of Pareto optimal solutions. Let ξ ∈ X model the SEMO solving mCOCZ and c = 1. We design the function h by setting α = 0, i.e., h(x) = |x ∩ S ∗ |. It is easy to verify that h is well-defined and r = |F ∗ | =(n′ +1)m′ . For the chain ξ ′ ∈Y ={0, 1}r with Y ∗ = {1r }, we construct the transition behavior as follows: for any ξ ′t /∈Y ∗ , ξ ′t+1 becomes 1r with probability 1/(|ξ ′t |0)1-1/m′ ; otherwise, ξ ′t+1 =ξ ′t . Thus, E[[τ ′ |ξ ′t ]]=(|ξ ′t |0)1-1/m′ . We then show that Eq. (1) holds with δ ≤ 1- 1 2nr . For any ξt = x /∈ X ∗ , assume that h(x) = k < r. Note that for any t, ξt contains only incomparable Pareto optimal solutions. Since the number of Pareto optimal solutions in the population never decreases, h(ξt+1) ≥ h(ξt) = k. We then analyze Pg :=P(h(ξt+1) = k+1). For i ∈ {1, . . . ,m′ }, let Ai ={f2i-1(s)|s∈x} and ai =|Ai|. Note that 1≤ai ≤n′ +1 and f2i-1 corresponds to the i-th block in the second half of a solution, i.e., the bits in positions n2 +(i-1)·n′+1 ∼ n2 +i·n′ . Assume that a1 ≥ a2 ≥ · · · ≥ am′ . We consider two cases: (1) am′ ≤ n′ . For the m′ -th block, there is at least one missing f2m′-1 value. Thus, for any value of f1 in A1, there is at least one corresponding solution s in x, and flipping a proper bit in the m′ -th block of s can generate a Pareto optimal solution with a new value of f2m′-1; this implies that h(ξt+1) = k + 1. Since |x| = k, Pg ≥ a1 nk ≥ k1/m′ /(nk). (2) am′ = n′ +1. Let NF ∗ = F ∗ \{f(s) | s ∈ x} denote the set of objective vectors in F ∗ that have not been found. For 1≤i≤m′ , let Bi={f ∗ 2i-1 | (f ∗ 1 , . . . , f ∗m) ∈ NF ∗ } and bi = |Bi|. Suppose j ∈argmaxi bi. For any value z∈Bj, there is at least one corresponding solution s ∈ x with f2j-1(s) = z, and flipping a proper bit of s can generate a Pareto optimal solution s′ with f(s′ ) ∈ NF ∗ ; otherwise, all the objective vectors in F ∗ with f ∗ 2j-1 =z have been found, which contradicts with z∈Bj. Thus, Pg ≥ bj nk ≥ (r - k)1/m′ /(nk). Combining these two cases, Pg ≥ min{ k1/m′ nk , (r-k)1/m′ nk }. By using the same analysis as Eqs. (2-3), Eq. (1) holds with where the inequality is by (1- 1 r-k )1-1/m′ ≤1- 1-1/m′ r-k and m′ = m2 ≥2. If k≤ r2 , Pg≥ k1/m′ nk ≥ 1 nr1-1/m′ , thus δ≤1- 1 2nr . If k> r2 , Pg ≥ (r-k)1/m′ nk ≥ (r-k)1/m′ nr , thus δ≤1- 1 2nr . By applying δ ≤ 1 - 1 2nr and E[[τ ′ |ξ ′0 = 0r ]] ≤ r1-1/m′ to Theorem 1, we get that E[T] ≤ r1-1/m′ · 2nr = O(nm). Note that r = (n′ + 1)m′ and 2m′ = m. Thus, the expected running time for m > 4 is O(nm/2+1 log n+nm) = O(nm). We further derive a tighter upper bound for m=4 by using a different ξ ′ ∈Y and analyzing Pg more carefully. The transition behavior of ξ ′ is as follows: ξ ′t+1 is generated by flipping the first 0-bit of ξ ′t with probability |ξ ′0 |/r; otherwise, ξ ′t+1 = ξ ′t . Thus, E[[τ ′ |ξ ′t ]] = rH|ξ′t |0 . We then show that Eq. (1) holds with δ ≤ 1- 1n by re-analyzing Pg. Note that x⊆S ∗ . For any z ∈A2, let Cz ∈argmaxC⊆x |C| s.t. ∀s∈ C, f3(s)=z and maxs∈C f1(s)-mins∈C f1(s)=|C|-1. Let u=argmaxs∈Cz f1(s) and v=argmins∈Cz f1(s). Then, a new Pareto optimal solution can be generated by flipping one 0-bit in the first block of u or flipping one 1-bit in the first block of v , whose probability is n′-(f1(u)-n/2)+f1(v)-n/2 = nk n′+1-|Cz| nk . By summing over all z∈A2, we get that one part of Pg is at least (n′+1)a2-k nk . By considering z ∈ A1 and exchanging f1 and f3 in the analysis above, we can similarly get that the other part of Pg is at least (n′+1)a1-k nk . Thus, Pg ≥ (n′+1)(a1+a2)-2k nk ≥ 2(n′+1)√k-2k nk , where the second inequality is by a1a2 ≥ k. Then, the δ in Eq. (1) satisfies that where the last inequality is by k < r = (n′+1)2 . By applying δ ≤ 1 - 1n and E[[τ ′ |ξ ′0 = 0r ]] = rHr to Theorem 1, we get that E[T] = O(n3 logn). Thus, the expected running time for m = 4 is O(n3 log n + n3 log n) = O(n3 log n). ## 6 Approximate Analysis In this section, we consider a variant of the GSEMO solving the WOMM problem, a generalization of the previously studied OneMinMax problem [Giel and Lehre, 2010; Osuna et al., 2017]. As presented in Definition 7, it is to maximize the weighted sum of 1-bits and 0-bits of a solution at the same time. Note that optimizing two linear functions simultaneously often appears in multi-objective combinatorial problems such as multi-objective minimum spanning trees and multi-objective knapsacks. The Pareto front F ∗ can be exponentially large, e.g., if wi = 2i , |F ∗ | = 2n ; thus we analyze the running time until finding an approximation of F ∗ (i.e., approximate analysis), instead of that until finding F ∗ (i.e., exact analysis) in the previous two sections. Theorem 5 shows that a good approximation can be obtained in polynomial time, which is consistent with the good empirical performance of MOEAs in solving multi-objective combinatorial problems [Zhou and Gen, 1999; Ishibuchi et al., 2015]. Definition 7 (WOMM). The pseudo-Boolean function WOMM: {0, 1}n → R2 is defined as We introduce two useful techniques into the GSEMO, i.e., the initialization strategy [Qian et al., 2013] and the diversityselecting mechanism [Horoba, 2009]. The former uses the (1+1)-EA (a simple single-objective EA) to optimize each objective separately, and then uses the optimum of each objective as the initial population. The latter is used to keep a diverse population and then achieve a good spread over the Pareto front. The objective space is divided into boxes, and the box value of a solution s is defined as where l > 1. In each iteration, a non-empty box (i.e., a box which contains at least one solution in the population) is chosen uniformly at random (u.a.r.); then a solution in the chosen box is selected u.a.r. to generate a new solution. Let W = ∑ni=1 wi for simplicity. Lemma 2 gives an upper bound on the number of non-empty boxes. Lemma 2. For the GSEMO on WOMM, there are at most B = 2(logl(1 +W/w1) + 2) non-empty boxes in the population. Proof. Assume that in the population, there are m > B different non-empty boxes: (a1, b1), (a2, b2), . . . , (am, bm), where a1 ≤ a2 ≤ . . . ≤ am, and for i < j, bi > bj if ai = aj. Then, it must hold that b1 ≥ b2 ≥ . . . ≥ bm; otherwise, it implies that for some i, ai < ai+1 and bi < bi+1, which contradicts with the fact that the solutions in the population are incomparable. Thus, we get a1 - b1 ≤ a2 - b2 ≤ . . . ≤ am - bm. Note that ai, bi ∈ {0, 1, . . . ,L}, where L = ⌈logl(1 + W/w1)⌉, and there are at most 2L + 1 different values for ai - bi. Meanwhile, m>B>2L + 2. Thus, for some i, ai-bi = ai+1-bi+1; this implies that ai=ai+1, bi =bi+1, which contradicts with the assumption. Lemma 3. A = {0n , 10n-1 , 010n-2 , . . . , 0n-1 1} is a 1n - approximation of F ∗ . Furthermore, the approximation ratio is max{ δ 1+δ , 1n }, if ∀i,wi+1≥(1+ δ) · wi, where δ>0. Proof. We only need to show that for any s∈{0, 1}n , there exists a solution s′ ∈ A such that ∀i ∈ {1, 2} : fi(s′ ) ≥ fi(s)/n. Let k=max{i|si =1}. We consider two cases: (1) If k = 0, s = 0n ∈ A. (2) If 1 ≤ k ≤ n, f2(0k-1 10n-k ) = W-wk ≥f2(s); f1(0k-1 10n-k )/f1(s) ≥ wk/(∑ki=1 wi) ≥ 1/n, since wi increases with i. For the furthermore clause, (∑ki=1 wi)/wk ≤ 1-1/(1+δ)k 1-1/(1+δ) ≤ 1+δ δ , since wi+1≥(1+δ) ·wi; then, wk/(∑ki=1 wi)≥δ/(1+δ). Thus, the lemma holds. Theorem 5. For maximizing WOMM, the GSEMO with an initialization strategy and a diversity-selecting mechanism can find a 1n -approximation of the Pareto front in O(n2 (logl n+ logl(wn/w1))) expected running time, where l > 1. Furthermore, the approximation ratio is max{ δ 1+δ , 1n }, if ∀i,wi+1 ≥ (1 + δ) · wi, where δ > 0. Proof. Since the (1+1)-EA optimizes a linear function in O(n logn) expected time [Droste et al., 2002], the initialization costs O(n log n) time, and the population will contain 1n and 0n . We then use Theorem 1 to analyze the expected number of iterations until finding a 1n -approximation of F ∗ . Let ξ ∈X model the variant of the GSEMO solving WOMM and c = 1n . We design a well-defined function h by setting α=0, i.e., h(x)=g(F ∗ , x, 1n ). Note that r = |F ∗ |. For the chain ξ ′ ∈Y={0, 1}r with Y ∗ ={1r }, its transition behavior is constructed as follows: for any ξ ′t /∈ Y ∗ , ξ ′t+1 becomes 1r with probability 1/(1 + log |ξ ′t |0); otherwise, ξ ′t+1 = ξ ′t . We can easily derive that E[[τ ′ |ξ ′t ]] = 1 + log |ξ ′t |0 for ξ ′t /∈ Y ∗ . Then, we are to investigate Eq. (1). For any ξt =x /∈X ∗ , assume h(x)=k<r. It is easy to verify that h(ξt+1)≥h(ξt)= k. Note that the solution 0n will be always contained in the population, since no other solution can weakly dominate it. Furthermore, it has a unique box value bl(0n )=(0, ⌈logl(1+ W/w1)⌉), since for any s ̸= 0n , ⌈logl(1+f1(s)/w1)⌉ ≥ 1. Assume that there are j solutions (denoted by s1 , s2 , . . . , sj ) in A which are not weakly dominated by any solution in ξt. Note that j ≥1; otherwise, h(ξt)=r, which makes a contradiction. For 1 ≤ i ≤ j, si can be generated by choosing 0n from ξt and flipping a specific 0-bit, whose probability is at least 1B · 1n (1 - 1n )n-1 ≥ 1 Ben . Once si is generated, it will be added into the population, since no solution in ξt weakly dominates it. For 1≤i≤j, let ki=h(x ∪ {si }), then ki≥k. By using the same analysis as Eqs. (2) and (3), we have ≤ 1 + (1/(Ben)) · ∑ji=1 (k-ki)/(r-k) ≤ 1-1/(Ben), where the second inequality is by log t≤t-1 for t>0, and the last is by ∑ji=1 (ki-k)= ∑ji=1 (h(x∪{si })-h(x)) ≥ h(x∪ A)-h(x) = r-k. Thus, Eq. (1) holds with δ ≤ 1-1/(Ben). By applying δ≤1 - 1 Ben and E[[τ ′ |ξ ′0 = 0r ]] ≤1+n log 2 to Theorem 1, we get that the expected running time for finding a 1n -approximation of F ∗ is at most (1+n log 2)Ben = O(n2 (logl n+logl(wn/w1))). For the furthermore clause, the analysis still holds by setting c=max{ δ 1+δ , 1n }. ## 7 Conclusion In this paper, we propose a general approach for deriving upper bounds on the running time of MOEAs. The key is to reduce the analysis of a given MOEA process to that of an easyto-analyze process. We apply this approach to diverse situations, including bi-objective and many-objective optimization as well as exact and approximate analysis. Our analysis can give the leading constants as well as bring an asymptotic improvement for some known running time bounds, which displays the strength of the proposed approach. In addition, our results provide some theoretical justification for MOEAs well solving multi-objective combinatorial problems in practice. ## Acknowledgments This work was supported in part by the National Key Research and Development Program of China (2017YFC0804002), the NSFC (61603367, 61672478), the YESS (2016QNRC001), and the Science and Technology Innovation Committee Foundation of Shenzhen (ZDSYS201703031748284). ## References 1. [Auger and Doerr, 2011] A. Auger and B. Doerr. Theory of Randomized Search Heuristics: Foundations and Recent Developments. World Scientific, Singapore, 2011. 2. [Bäck, 1996] T. Bäck. Evolutionary Algorithms in Theory and Practice: Evolution Strategies, Evolutionary Pro- gramming, Genetic Algorithms. Oxford University Press, Oxford, UK, 1996. 3. [Coello Coello and Lamont, 2004] C. A. Coello Coello and G. B. Lamont. Applications of Multi-Objective Evolution- ary Algorithms. World Scientific, Singapore, 2004. 4. [Doerr et al., 2013] B. Doerr, B. Kodric, and M. Voigt. Lower bounds for the runtime of a global multi-objective evolutionary algorithm. In CEC, pages 432-439, 2013. 5. [Doerr et al., 2016] B. Doerr, W. Gao, and F. Neumann. Runtime analysis of evolutionary diversity maximization for ONEMINMAX. In GECCO, pages 557-564, 2016. 6. [Droste et al., 2002] S. Droste, T. Jansen, and I. Wegener. On the analysis of the (1+1) evolutionary algorithm. Theoret- ical Computer Science, 276(1-2):51-81, 2002. 7. [Friedrich et al., 2010] T. Friedrich, N. Hebbinghaus, and F. Neumann. Plateaus can be harder in multi-objective optimization. Theoretical Computer Science, 411(6):854- 864, 2010. 8. [Friedrich et al., 2011] T. Friedrich, C. Horoba, and F. Neu- mann. Illustration of fairness in evolutionary multi- objective optimization. Theoretical Computer Science, 412(17):1546-1556, 2011. 9. [Giel and Lehre, 2010] O. Giel and P. K. Lehre. On the ef- fect of populations in evolutionary multi-objective optimi- sation. Evolutionary Computation, 18(3):335-356, 2010. 10. [Giel, 2003] O. Giel. Expected runtimes of a simple multi- objective evolutionary algorithm. In CEC, pages 1918- 1925, 2003. 11. [Horoba, 2009] C. Horoba. Analysis of a simple evolution- ary algorithm for the multiobjective shortest path problem. In FOGA, pages 113-120, 2009. 12. [Ishibuchi et al., 2015] H. Ishibuchi, N. Akedo, and Y. No- jima. Behavior of multiobjective evolutionary algorithms on many-objective knapsack problems. IEEE Transactions on Evolutionary Computation, 19(2):264-283, 2015. 13. [Laumanns et al., 2004] M. Laumanns, L. Thiele, and E. Zit- zler. Running time analysis of multiobjective evolutionary algorithms on pseudo-Boolean functions. IEEE Transac- tions on Evolutionary Computation, 8(2):170-182, 2004. 14. [Neumann and Theile, 2010] F. Neumann and M. Theile. How crossover speeds up evolutionary algorithms for the multi-criteria all-pairs-shortest-path problem. In PPSN, pages 667-676, 2010. 15. [Neumann and Witt, 2010] F. Neumann and C. Witt. Bioin- spired Computation in Combinatorial Optimization: Al- gorithms and Their Computational Complexity. Springer- Verlag, Berlin, Germany, 2010. 16. [Neumann, 2007] F. Neumann. Expected runtimes of a sim- ple evolutionary algorithm for the multi-objective mini- mum spanning tree problem. European Journal of Op- erational Research, 181(3):1620-1629, 2007. 17. [Neumann, 2012] F. Neumann. Computational complex- ity analysis of multi-objective genetic programming. In GECCO, pages 799-806, 2012. 18. [Osuna et al., 2017] E. C. Osuna, W. Gao, F. Neumann, and D. Sudholt. Speeding up evolutionary multi-objective op- timisation through diversity-based parent selection. In GECCO, pages 553-560, 2017. 19. [Qian et al., 2013] C. Qian, Y. Yu, and Z.-H. Zhou. An anal- ysis on recombination in multi-objective evolutionary op- timization. Artificial Intelligence, 204:99-119, 2013. 20. [Qian et al., 2016] C. Qian, K. Tang, and Z.-H. Zhou. Se- lection hyper-heuristics can provably be helpful in evolu- tionary multi-objective optimization. In PPSN, pages 835- 846, 2016. 21. [Yu et al., 2015] Y. Yu, C. Qian, and Z.-H. Zhou. Switch analysis for running time analysis of evolutionary algo- rithms. IEEE Transactions on Evolutionary Computation, 19(6):777-792, 2015. 22. [Zhou and Gen, 1999] G. Zhou and M. Gen. Genetic algo- rithm approach on multi-criteria minimum spanning tree problem. European Journal of Operational Research, 114(1):141-152, 1999.