A General Approach to Running Time Analysis of Multi-objective Evolutionary Algorithms

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).