Maintaining Diversity Provably Helps in Evolutionary Multimodal Optimization

1 Introduction

In the real world, there exist a class of optimization problems where multiple (local) optimal solutions in the solution space correspond to a single point in the objective space, such as in truss structure optimization [Reintjes, 2022], space mission design [Schutze et al., 2011], rocket engine design [Kudo et al., 2011], and functional brain imaging [Sebag et al., 2005]. Such problems belong to multi-modal optimization problems (MMOPs). Note that there are usually two types of MMOPs [Deb and Saha, 2010; Preuss et al., 2021]: one with multiple local optima (typically with different objective function values) and the other with multiple (local) optimal solutions having identical objective function value. In this work, MMOPs we refer to is the latter.

Evolutionary algorithms (EAs) 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 generating new offspring solutions and replacing inferior ones. EAs, due to their population-based nature and the ability of performing the search globally, have become a popular tool to solve MMOPs [Das et al., 2011; Cheng et al., 2018; Tian et al., 2021; Pan et al., 2023]. Since 1990's, to deal with MMOPs, many ideas have been proposed to improve EAs' ability, including using niching technique [P´etrowski, 1996; Shir, 2012; Preuss, 2015], designing novel reproduction operators and population update methods [Kennedy, 2010; Storn and Price, 1997; Liang et al., 2024], and transforming an MMOP into a multiobjective optimization problem [Wessing et al., 2013; Cheng et al., 2018; Liu et al., 2022].

In contrast to algorithm design of EAs, theoretical studies (e.g., running time complexity analyses) in the area are relatively underdeveloped. This is mainly because sophisticated behaviors of EAs can make theoretical analysis rather difficult. Nevertheless, since the 2000s, there is increasingly interest in rigorously analyzing EAs. Some studies worked on running time analysis tools, such as drift analysis [He and Yao, 2001; Oliveto and Witt, 2011; Doerr et al., 2012], fitness level method [Wegener, 2002; Sudholt, 2013; Dang and Lehre, 2015], and switch analysis [Yu et al., 2015; Yu and Qian, 2015; Qian et al., 2016; Bian et al., 2018]. Some other studies devoted themselves to the analysis of various EAs for solving synthetic or combinatorial optimization problems [Neumann and Witt, 2010; Auger and Doerr, 2011; Zhou et al., 2019; Doerr and Neumann, 2020]. These theoretical results can help understand the mechanisms of EAs and inspire the design of more efficient EAs in practice.

In this paper, we consider running time analysis of EAs for MMOPs. We analytically show that when solving MMOPs, considering the diversity of the solutions in the solution space can be beneficial for the search of EAs. Specifically, we propose a simple method to deal with multiple (local) optimal solutions via comparing their crowdedness in the solution space. We consider both single-objective and multiobjective optimization cases, and incorporate the proposed method into a classical single-objective EA, (µ+1)-GA, and two well-established multi-objective EAs (MOEAs), NSGAII and SMS-EMOA.

For the single-objective case, we prove that the expected running time of (µ + 1)-GA using the proposed method for solving (i.e., finding the global optimum) Jump [Droste et al., 2002], a widely studied single-objective problem, is O(µ2 4k + µn log n + n√k(µ log µ + log n)), where n is the problem size, and k ≤ n/4, a parameter of Jump. Since the expected running time of the original (µ + 1)-GA for solving Jump is O((40eµ(µ + 1))k ( 1 µ-1 )k-1 ) 10e 9 nk (k-1)! + n√k(µ logµ + logn) = O(µ√k(40e2µn/k)k ) [Doerr and Qu, 2023c], our method can bring a polynomial acceleration when k is a constant, e.g., k = 4, and bring an exponential acceleration when k is large, e.g., k = n1/4 .

For the multi-objective case, we prove that the expected running time of NSGA-II and SMS-EMOA (with crossover) using the proposed method for solving (i.e., finding the whole Pareto front) OneJumpZeroJump [Doerr and Zheng, 2021], a widely studied bi-objective problem, is both O(µ2 4k + µn logn), where µ is the population size, n is the problem size, and k ≤ n/4, a parameter of OneJumpZeroJump. Note that the expected running time of the original NSGA-II [Doerr and Qu, 2023c] and SMS-EMOA (analyzed in this paper) for solving OneJumpZeroJump is O(µ2√k(Cn/k)k ), where C is a constant and k = o(√n). Thus, our method can also bring a polynomial acceleration when k is a constant, and bring an exponential acceleration when k is large, e.g., k = n1/4 . The main reason for the acceleration resulting from our method is that it can explicitly preserve solutions with high diversity, and thus can help enhance the exploration by working with crossover. Experiments are also conducted to validate the theoretical results. In addition, it is worth mentioning that since the expected running time of SMS-EMOA (without crossover) for solving OneJumpZeroJump is O(µnk ) [Bian et al., 2023], our results also contribute to the theoretical understanding of the effectiveness of using crossover in SMS-EMOA.

2 Preliminaries

In this section, we first introduce multimodal optimization and the considered benchmark problems, Jump and OneJumpZeroJump, followed by the description of the studied algorithms, (µ + 1)-GA, NSGA-II and SMS-EMOA.

2.1 Multimodal Optimization

Multimodal optimization refers to the optimization scenario that multiple (local) optimal solutions in the solution space have the same objective value, which arises in many realworld applications, e.g., flow shop scheduling [Basseur et al., 2002], rocket engine optimization [Kudo et al., 2011], and architecture design [Tian et al., 2021]. In this paper, we study two pseudo-Boolean (i.e., the solution space X = {0, 1}n ) multi-modal optimization problems (MMOPs), Jump and OneJumpZeroJump, which have been widely used in EAs' theoretical analyses [Droste et al., 2002; Doerr and Zheng, 2021; Bian et al., 2023; Doerr and Qu, 2023a; Doerr and Qu, 2023b; Doerr and Qu, 2023c; Lu et al., 2024].

The Jump problem as presented in Definition 1, is to maximize the number of 1-bits of a solution, except for a valley around 1n (the solution with all 1-bits) where the number of 1-bits should be minimized. Its optimal solution is 1n with function value n + k. We can see that 1n is global optimal, and any solution with (n-k) 1-bits is local optimal with objective value n; thus, the Jump problem is multi-modal, since multiple local optimal solutions correspond to an objective value.

Definition 1 ([Droste et al., 2002]). The Jump problem is to find an n bits binary string which maximizes

where k ∈ Z ∧ 2 ≤ k < n, and |x|1 denotes the number of 1-bits in x.

The OneJumpZeroJump problem is a bi-objective counterpart of the Jump problem. For multi-objective optimization, several objective functions (which are usually conflicting) need to be optimized simultaneously, and thus there does not exist a canonical complete order in the solution space X. And usually the Pareto domination relation is used (Definition 2) to compare solutions. A solution is Pareto optimal if it is not dominated by any other solution in X, and the set of objective vectors of all the Pareto optimal solutions is called the Pareto front. The goal of multi-objective optimization is to find the Pareto front or its good approximation.

Definition 2. Let f = (f1, f2, . . . , fm) : X → Rm be the objective vector. For two solutions x and y ∈ X :

• x weakly dominates y (denoted as x ⪰ y) if for any 1 ≤ i ≤ m, fi(x) ≥ fi(y);
• x dominates y (denoted as x ≻ y) if x ⪰ y and fi(x) > fi(y) for some i;
• x and y are incomparable if neither x ⪰ y nor y ⪰ x.

As presented in Definition 3 below, the first objective of OneJumpZeroJump is the same as the Jump problem, while the second objective is isomorphic to the first one, with the roles of 1-bits and 0-bits exchanged. The Pareto front of the OneJumpZeroJump problem is {(a, n + 2k - a) | a ∈ [2k..n] ∪ {k, n + k}}, whose size is n - 2k + 3, and the Pareto optimal solution corresponding to (a, n + 2k - a), a ∈ [2k..n] ∪ {k, n + k}, is any solution with (a - k) 1bits. Note that we use [l..r] (where l, r ∈ Z, l ≤ r) to denote the set {l, l +1, . . . , r} of integers throughout the paper. Since any solution with (a - k) 1-bits (where a ∈ [2k..n]) is Pareto optimal with objective vector (a, n + 2k - a), the OneJumpZeroJump problem is multimodal.

Definition 3 ([Doerr and Zheng, 2021]). The OneJumpZeroJump problem is to find n bits binary strings which maximize

where k ∈ Z ∧ 2 ≤ k < n/2, and |x|1 and |x|0 denote the number of 1-bits and 0-bits in x, respectively.

2.2 Evolutionary Algorithms

The (µ + 1)-GA algorithm is a classical genetic algorithm for solving single-objective optimization problems. The algorithm starts from an initial population of µ random solutions. In each iteration, it selects a solution x from the population P randomly as the parent solution. Then, with probability pc, it selects another solution y and applies uniform crossover on x and y to generate an offspring solution x′ ; otherwise, x′ is set as the copy of x. The uniform crossover operator exchanges each bit of two solutions independently with probability 1/2. Note that uniform crossover actually produces two solutions, but we only pick the first one. Afterwards, bitwise mutation is applied, which flips each bit of a solution independently with probability 1/n, on x′ to generate one offspring solution. Then, one solution with the worst objective value is removed. The complete procedure is provided in the supplementary material due to space limitation.

Now we introduce two well-established MOEAs, NSGAII [Deb et al., 2002] and SMS-EMOA [Beume et al., 2007]. The NSGA-II algorithm in Algorithm 1 adopts a (µ + µ) steady state mode. It starts from an initial population of µ random solutions (line 1). In each iteration, it selects µ solutions from the current population to form the parent population Q (line 4). Then, for each pair of the solutions in Q, uniform crossover and bit-wise mutation operators are applied sequentially to generate two offspring solutions x′′ and y′′ (lines 6-12), where the uniform crossover operator is applied with probability pc. After µ offspring solutions have been generated in P ′ , the solutions in P ∪ P ′ are partitioned into non-dominated sets R1, . . . , Rv (line 14), where R1 contains all the non-dominated solutions in P∪P ′ , and Ri (i ≥ 2) contains all the non-dominated solutions in (P ∪ P ′ ) \ ∪i-1 j=1Rj . Then, the solutions in R1, R2, . . . , Rv are added into the next population, until the population size exceeds µ (lines 15-18). For the critical set Ri whose inclusion makes the population size larger than µ, the crowding distance is computed for each of the contained solutions (line 19). Crowding distance reflects the level of crowdedness of solutions in the population. For each objective fj, 1 ≤ j ≤ m, the solutions in Ri are sorted according to their objective values in ascending order, and assume the sorted list is x1 , x2 , . . . , xk . Then, the crowding distance of the solution xl with respect to fj is set to ∞ if l ∈ {1, k} and (fj(xl+1 ) - fj(xl-1 ))/(fj(xk ) - fj(x1 )) otherwise. The final crowding distance of a solution is the sum of the crowding distance with respect to each objective. Finally, the solutions in Ri are selected to fill the remaining population slots where the solutions with larger crowding distance are preferred (line 20).

SMS-EMOA shares a similar framework with (µ + 1)- GA, and the main difference is that SMS-EMOA uses nondominated sorting and hypervolume indicator to evaluate the quality of a solution and update the population. Specifically, after the offspring solution x′′ is generated, the solutions in P ∪ {x′′ } are first partitioned into non-dominated sets R1, . . . ,Rv, and then one solution z ∈ Rv that minimizes ∆r(x,Rv) is removed, where ∆r(x,Rv) = HVr(Rv) - HVr(Rv \ {x}). Note that HVr(X) = Λ( ∪x∈X {f ′ ∈ Rm | ∀1 ≤ i ≤ m : ri ≤ f ′ i ≤ fi(x)}) denotes the hyper-

Algorithm 1 NSGA-II

Input: objective functions f1, f2 . . . , fm, population size µ,

• probability pc of using crossover Output: µ solutions from {0, 1}n 1: P ← µ solutions uniformly and randomly selected from {0,1}n with replacement; 2: while criterion is not met do 3: let P ′ = ∅; 4: generate a parent population Q of size µ; 5: for each pair of the parent solutions x and y in Q do 6: sample u from uniform distribution over [0, 1]; 7: if u < pc then 8: apply uniform crossover on x and y to generate two solutions x′ and y′ 9: else 10: set x′ and y′ as copies of x and y, respectively 11: end if 12: apply bit-wise mutation on x′ and y′ to generate x′′ and y′′ , respectively, and add x′′ and y′′ into P ′ 13: end for 14: partition P ∪P ′ into non-dominated sets R1, . . . ,Rv; 15: let P = ∅, i = 1; 16: while |P ∪Ri| < µ do 17: P = P ∪Ri, i = i +1 18: end while 19: assign each solution in Ri with a crowding distance; 20: sort the solutions inRi in ascending order by crowding distance, and add the last µ - |P | solutions into P 21: end while 22: return P

volume of a solution set X with respect to a reference point r ∈ Rm (satisfying ∀1 ≤ i ≤ m, ri ≤ minx∈X fi(x)), i.e., the volume of the objective space between the reference point and the objective vectors of the solution set, where Λ denotes the Lebesgue measure. A larger hypervolume implies a better approximation in terms of convergence and diversity.

Note that in this paper, we consider different methods to select the parent solutions in line 4 of Algorithm 1, i.e., fair selection which selects each solution in P once (the order of solutions is random), uniform selection which selects parent solutions from P independently and uniformly at random for µ times, and binary tournament selection which first picks two solutions randomly from the population P with replacement and then selects a better one for µ times. Furthermore, for all the three algorithms studied in this paper, the probability pc of using crossover belongs to the interval [Ω(1), 1 -Ω(1)].

3 Proposed DiversityMaintenance Method

In this section, we introduce the proposed method that is used to maintain the diversity of solutions in the population update procedure of (µ + 1)-GA, NSGA-II and SMS-EMOA. The general idea of our method is very simple - we take into account the (Hamming) distance of the solutions into the population update procedure of the algorithm when the solutions have the same objective value, such that the distant solutions in the solution space can be preserved.

First, let us consider (µ + 1)-GA. In the original population update procedure of (µ + 1)-GA, the solution with the minimal objective value is removed, with the ties broken uniformly. That is, the solutions with the minimal objective value are treated equally, regardless of their structure in the solution space. In order to take advantage of the genotype information of the solutions, our method will further consider the diversity of the solutions with the minimal objective value, and preserve the solutions with larger diversity. Specifically, the procedure of selecting solution z to remove is changed as follows (note that H(·, ·) denotes the Hamming distance of two solutions):

let S = {z | f(z) = minx∈P∪{x′′ } f(x)}; if |S| ≤ 2 then let z be a solution randomly selected from S else let (˜x, ˜ y) = arg maxx,y∈S,x̸=y H(x, y); let z be a solution randomly selected from S\{˜x, ˜ y} end if

Now we consider NSGA-II. When computing the crowding distance for the solutions in the critical non-dominated set, i.e., Ri in line 19 of Algorithm 1, the solutions with the same objective vector must be sorted together. Then, the solutions placed in the first or the last position among these solutions will be assigned a crowding distance larger than 0, while the other solutions will only be with a crowding distance of 0. In this case, the solutions (corresponding to the same objective vector) placed in the boundary positions will be preferred by the measure of crowding distance. In the original NSGAII, the solutions (corresponding to the same objective vector) are deemed as identical and their relative positions in the calculation of crowding distance are actually assigned randomly, which may harm the efficiency of the algorithm, as we will show later.

As can be seen, here the main idea of our method is to modify the solutions' position in the crowding distance calculation, so that the distant solutions in the solution space are preferred. Specifically, after the solutions in the critical non-dominated set Ri are sorted according to some objective fj in ascending order (here assuming that the sorted list is x1 , x2 , . . . , xk ), we reorder the list as follows:

let G = {fj(xl ) | l ∈ [1..k]}; for g ∈ G do let Sg = {l ∈ [1..k] | fj(xl ) = g}; let (˜a, ˜b) = arg maxa,b∈Sg ,a̸=b H(xa , xb ); let ˆl = min{l ∈ Sg}, ˆ r = max{r ∈ Sg}; swap the positions of x˜a and xˆl ; swap the positions of x˜b and xˆ r end for

That is, for the solutions with the same objective value, the pair of solutions with the largest Hamming distance are put in the first and last positions among these solutions. The crowding distance is then calculated based on the reordered list.

Finally, let us consider SMS-EMOA. In the original population update procedure of SMS-EMOA, the solution in the last non-dominated set Rv with the least ∆ value is removed. However, when several solutions in Rv have the same objective vector, they will all have a zero ∆ value, implying that one of them will be removed randomly. To make use of their crowdedness in the solution space, we first select the most crowded objective vector (i.e., the objective vector corresponding to the most solutions), and then randomly remove one of those corresponding solutions excluding the two solutions having the largest Hamming distance. Specifically, the procedure of selecting solution z that minimizes ∆r(x,Rv) will be changed to:

let G = {f(x) | x ∈ Rv}; let g∗ = argmaxg∈G |{x ∈ Rv | f(x) = g}|; let S = {x ∈ Rv | f(x) = g∗ }; if |S| > 2 then let (˜x, ˜ y) = arg maxx,y∈S,x̸=y H(x, y); let z be a solution randomly selected from S\{˜x, ˜ y}; else let z = arg minx∈Rv ∆r(x, Rv) end if

4 (µ+ 1)-GA on Jump

In this section, we show the expected running time of (µ+1)- GA for solving the Jump problem. Note that the running time of EAs is often measured by the number of fitness evaluations, the most time-consuming step in the evolutionary process. We prove in Theorem 1 that the expected number of fitness evaluations of (µ + 1)-GA using the proposed diversity maintenance method for solving Jump is O(µ2 4k+µn log n+ n√k(µ logµ + logn)). The proof idea is to divide the optimization procedure into three phases, where the first phase aims at driving the population towards the local optimum, i.e., all the solutions in the population have (n - k) 1-bits, the second phase aims at finding two local optimal solutions with the largest Hamming distance 2k, and the third phase applies uniform crossover to these two distant local optimal solutions to find the global optimum.

Theorem 1. For (µ + 1)-GA solving Jump with k ≤ n/4, if using the diversity maintenance method, and a population size µ such that µ ≥ 2, then the expected number of fitness evaluations for finding the global optimal solution 1n is O(µ2 4k + µn log n + n√k(µ log µ + log n)).

Proof. We divide the optimization procedure into three phases. The first phase starts after initialization and finishes until all the solutions in the population have (n - k) 1-bits; the second phase starts after the first phase, and finishes until the largest Hamming distance between solutions in the population reaches 2k; the third phase starts after the second phase, and finishes when the global optimal solution 1n is found. By Lemma 1 in [Dang et al., 2018], we can directly derive that the expected number of fitness evaluations of the first phase is O(n√k(µ logµ + logn)). Note that after the first phase finishes, all the solutions in the population will always have (n - k) 1-bits (except that 1n is found), because once a solution with the number of 1-bits in [0..n - k - 1] ∪ [n - k + 1..n - 1] is generated, it will be removed in the population update procedure.

Next, we consider the second phase. Let Jmax = maxx,y∈P H(x, y) denote the maximal Hamming distance between two solutions in the population P , where H(·, ·) denotes the Hamming distance of two solutions. We show that Jmax will not decrease, unless 1n is found. Assume that x and y are a pair of solutions which have the maximal Hamming distance. In each iteration, x or y can be removed only if the newly generated solution z has exactly (n - k) 1-bits and survives in the population update procedure (note that we can pessimistically assume that 1n has not been found), implying that all the solutions in P ∪ {z} have (n - k) 1-bits. According to the diversity maintenance method introduced in Section 3, there must exist a pair of solutions in P ∪{z} with Hamming distance at least Jmax, which will not be removed, implying that Jmax cannot decrease.

We then show that Jmax can increase to 2k in O(µn log n) expected number of fitness evaluations. Assume that currently Jmax = 2j < 2k, and (x, y) is a pair of corresponding solutions, i.e., H(x, y) = 2j. Let D(x, y) = {i ∈ [1..n] | xi = yi} denote the set of positions which have the identical bit for x and y. Suppose that the number of positions in D(x, y) that have 1-bits is l, then we have |x| + |y| = 2l+2j = 2n-2k, implying that l = n-k-j. Then, we can derive that the number of positions in D(x, y) that have 0-bits is n-2j - l = k - j. In each iteration, the probability of selecting x as a parent solution is 1/µ. In the reproduction procedure, a solution z with (n-k) 1-bits and H(y, z) = 2j+2 can be generated from x if crossover is not performed (whose probability is 1-pc), and a 1-bit together with a 0-bit from the positions in D(x, y) are flipped by bit-wise mutation (whose probability is ((n-k-j)(k-j)/n2 ) · (1-1/n)n-2 ). Thus, the probability of generating z is at least

This implies that the expected number of fitness evaluations for increasing Jmax to 2k is at most

where the inequality holds because n-k-j ≥ n-2k+1 ≥ n/2 for j ≤ k - 1 and k ≤ n/4, and the equality holds because pc ∈ [Ω(1), 1-Ω(1)], and the k-th Harmonic number Hk = ∑kj=1 1/j = O(log n).

Finally, we consider the third phase. When Jmax increases to 2k, there must exist two solutions x∗ and y∗ such that |x∗ |1 = |y∗ |1 = n - k and H(x∗ , y∗ ) = 2k. By selecting x∗ and y∗ as a pair of parent solutions, exchanging all the 2k different bits by uniform crossover, and flipping none of bits in bit-wise mutation, the solution 1n can be generated, whose probability is at least

Thus, the expected number of fitness evaluations of the third phase, i.e., finding 1n , is O(µ2 22k ).

Combining the three phases, the total expected number of fitness evaluations is O(n√k(µ logµ + logn)) + O(µn log n) + O(µ2 22k ), which leads to the theorem.

The expected number of fitness evaluations of the original (µ + 1)-GA for solving Jump has been shown to be O(µ√k(40e2µn/k)k ) when k = o(√n) [Doerr and Qu, 2023c], and O(µn logµ + nk /µ + nk-1 logµ) when k is a constant [Dang et al., 2018]. Thus, our result in Theorem 1 shows that when k is large, e.g., k = n1/4 , the expected running time can be reduced by a factor of Ω((10e2µn/k)k /(µn logn)), which is exponential; when k is a constant, the expected running time can be reduced by a factor of Ω(max(1, nk-1 /(µ2 logn))), which is polynomial. The main reason for the acceleration is that the proposed method prefers solutions with large Hamming distance, making the population quickly find two solutions with sufficiently large Hamming distance and then allowing crossover to generate the optimal solution.

Note that Dang et al. [2016] also proposed several types of diversity maintenance mechanisms for (µ+1)-GA solving the Jump problem, including mechanisms similar to ours, e.g., total Hamming distance mechanism which maximizes the sum of the Hamming distances between all solutions. Based on the running time bounds, their mechanisms and our proposed mechanism have own advantages. For example, when k and µ are fairly large (e.g., k = √n and µ = √n/3), the expected number of fitness evaluations of (µ + 1)-GA using the total Hamming distance mechanism and our proposed mechanism is O(4√n ) and O(n4√n ), respectively, implying that the total Hamming distance mechanism is faster by a factor of Θ(n). For small k and large µ (e.g., k = Θ(1) and µ = Θ(n)), however, our proposed mechanism achieves a better running time O(n2 log n), compared to O(n3 log n) for the total Hamming distance mechanism. Similar observation also holds for the other mechanisms.

5 NSGA-II on OneJumpZeroJump

In the previous section, we have shown that the proposed method can make (µ + 1)-GA faster on Jump, while in this section, we show that similar acceleration can be achieved in multi-objective scenario. In particular, we prove in Theorem 2 that the expected number of fitness evaluations of NSGA-II using the proposed diversity maintenance method for solving OneJumpZeroJump is O(µ2 4k + µn log n). The proof idea is to divide the optimization procedure into three phases: 1) to find the inner part F ∗ I = {(a, n + 2k - a) | a ∈ [2k..n]} of the Pareto front; 2) to find two Pareto optimal solutions with (n - k) (or k) 1-bits that have the same objective vector (n, 2k) (or (2k, n)) but have Hamming distance 2k; 3) to find the remaining two extreme vectors in the Pareto front, i.e., {(k, n + k), (n + k, k)}, corresponding to the two Pareto optimal solutions 0n and 1n , which can be accomplished by applying uniform crossover to those two distant solutions found in the last phase.

Theorem 2. For NSGA-II solving OneJumpZeroJump with k ≤ n/4, if using the diversity maintenance method, and a population size µ such that µ ≥ 4(n - 2k + 3), then the expected number of fitness evaluations for finding the Pareto front is O(µ2 4k + µn log n).

Proof. We divide the optimization procedure into three phases: the first phase starts after initialization and finishes until the inner part F ∗ I of the Pareto front is entirely covered by the population P; the second phase starts after the first phase and finishes until the maximal Hamming distance of the solutions which have the objective vector (n, 2k) (i.e., have (n - k) 1-bits) reaches 2k; the third phase starts after the second phase and finishes when the extreme Pareto optimal solution 1n is found. Note that the analysis of the last two phases for finding the other extreme Pareto optimal solution 0n is similar. By Lemma 4 in [Doerr and Qu, 2023a], we can directly derive that the expected number of fitness evaluations of the first phase is O(µn log n).

Then, we consider the second phase. Let Jmax = maxx,y∈P:f(x)=f(y)=(n,2k) H(x, y) denote the maximal Hamming distance of the solutions which have the objective vector (n, k). We will show that Jmax will not decrease after using the proposed diversity maintenance method. Because (n, k) is a point in the Pareto front, any corresponding solution (which has (n - k) 1-bits) must have rank 1, i.e., belong to R1 in the non-dominated sorting procedure. If |R1| ≤ µ, all the solutions in R1 will be maintained in the next population, implying the claim holds. If |R1| > µ, the crowding distance of the solutions in R1 needs to be computed. When the solutions in R1 are sorted according to some objective, the solutions corresponding to the objective vector (n, k) will be resorted using our proposed method, and a pair of solutions (denoted as x and y) with the largest Hamming distance will be put in the first or the last position among these solutions, thus having a crowding distance larger than 0. Since (n, k) has been obtained, the solutions in R1 must be Pareto optimal. Note that the size of the Pareto front is n-2k +3, thus the number of different objective vectors of the solutions in R1 is at most n-2k+3. For any objective vector, the corresponding solutions will have crowding distance larger than 0 only if they are put in the first or the last position among these solutions, when they are sorted according to some objective. As OneJumpZeroJump has two objectives, at most four solutions corresponding to one objective vector can have a crowding distance larger than 0, implying that at most 4(n-2k+3) solutions in R1 can have a positive crowding distance. Thus, x and y are among the best 4(n-2k+3) solutions in R1. As the population size µ ≥ 4(n-2k+3), they must be included in the next population, implying that Jmax will not decrease.

Next, we show that Jmax can increase to 2k in O(µn log n) expected number of fitness evaluations. The proof is similar to the argument in the proof of Theorem 1, and the main difference is that NSGA-II uses three different parent selection strategies, and produces µ solutions instead of one in each generation. For NSGA-II using uniform selection, Eq. (1) also holds here. But since in each generation, µ/2 pairs of parent solutions will be selected for reproduction, the probability of increasing Jmax is at least

which hold by 1 + a ≤ ea for any a ∈ R, and pc ∈ [Ω(1), 1 -Ω(1)]. For fair selection, each solution in the current population will be selected once; thus by Eq. (1), the probability of increasing Jmax in each generation is at least (1-pc)(k-j)(n-k-j)/(en2 ) = Ω((k-j)(n-k-j)/n2 ). For binary tournament selection, the above equation also holds, but the detailed analysis is provided in the supplementary material due to space limitation. Then, similar to Eq. (2), we can derive that the total expected number of generations for increasing Jmax to 2k is at most O(n log n).

Finally, we consider the third phase. For uniform selection, Eq. (3) directly applies here. Thus, the probability of finding 1n in each generation is at least

For fair selection, the probability that two solutions are paired together in each generation is at least 1/µ , because we can assume that the position of one solution is fixed and it is sufficient to put the other solution beside it. Thus by Eq. (3), the probability of finding 1n in each generation is at least Ω(1/(µ4k )). For binary tournament selection, the above equation also holds, and the detailed analysis is provided in the supplementary material. Thus, the expected number of generations for finding 1n is O(µ4k ).

Combining the three phases, the total expected number of generation is O(µ4k + n logn), implying that the expected number of fitness evaluations is O(µ2 4k +µn log n), because each generation of NSGA-II requires to evaluate µ offspring solutions. Thus, the theorem holds.

The expected number of fitness evaluations of the original NSGA-II for solving OneJumpZeroJump has been shown to be O(µ2√k(Cn/k)k ) if k = o(√n) [Doerr and Qu, 2023c], where C is a constant. Thus, our result in Theorem 2 shows that the expected running time can be reduced by a factor of Ω(√k(Cn/(4k))k / log n), which is polynomial for a constant k and exponential for large k, e.g., k = n1/4 . The main reason for the acceleration is similar to that of (µ+ 1)- GA, i.e., the proposed diversity maintenance method prefers distant solutions in the solution space, and thus enhances the exploration ability of the algorithm by working with the crossover operator.

6 SMS-EMOA on OneJumpZeroJump

In this section, we consider SMS-EMOA on the OneJumpZeroJump problem. We prove in Theorem 3 that the expected number of fitness evaluations of the original SMSEMOA is O(µ2√k(Cn/k)k ), where C is a constant, and then prove in Theorem 4 that the running time reduces to O(µ2 4k + µn logn) if using the proposed diversity maintenance method. Their proofs are similar to that of Theorem 2. That is, we divide the optimization procedure into three phases, where the first phase aims at finding the inner part F ∗ I of the Pareto front, the second phase aims at finding two Pareto optimal solutions with (n - k) (or k) 1-bits that have Hamming distance 2k, and the last phase aims at finding the remaining two extreme Pareto optimal solutions 1n and 0n by applying uniform crossover to those two distant solutions found in the last phase. The detailed proofs are provided in the supplementary material due to space limitation.

Theorem 3. For SMS-EMOA solving OneJumpZeroJump with k = o(√n), if using a population size µ such that µ ≥ 2(n-2k+3), then the expected number of fitness evaluations for finding the Pareto front is O(µ2√k(Cn/k)k ), where C is a constant.

Theorem 4. For SMS-EMOA solving OneJumpZeroJump with k ≤ n/4, if using the diversity maintenance method, and a population size µ such that µ ≥ 2(n-2k+3), then the expected number of fitness evaluations for finding the Pareto front is O(µ2 4k + µn log n).

By comparing Theorems 3 and 4, we can find that by using the diversity maintenance method, the expected number of fitness evaluations of SMS-EMOA for solving OneJumpZeroJump can be reduced by a factor of Ω(√k(Cn/(4k))k / logn), which is polynomial for a constant k and exponential for large k, e.g., k = n1/4 . The main reason for the acceleration is just similar to that of NSGAII. That is, for Pareto optimal solutions with (n - k) (or k) 1-bits, the proposed method always maintains the two ones with the largest Hamming distance in the population update procedure, which enables the algorithm to quickly increase the largest Hamming distance until reaching the maximum 2k, i.e., finishing the second phase.

7 Experiments

In the previous sections, we have proved that using the proposed diversity maintenance method in (µ+1)-GA, NSGA-II and SMS-EMOA can bring a significant acceleration for solving Jump and OneJumpZeroJump. However, as only upper bounds on the running time of the original algorithms have been derived, we conduct experiments to examine their actual performance to complement the theoretical results.

Specifically, we set the problem size n of Jump and OneJumpZeroJump from 10 to 30, with a step of 5, and set the parameter k of the two problems to 4. For (µ+1)-GA solving Jump, the population size µ is set to 2, as suggested in Theorem 1; while for NSGA-II and SMS-EMOA solving OneJumpZeroJump, the population size µ is set to 4(n-2k +3) and 2(n - 2k + 3), respectively, as suggested in Theorems 2 and 4. For each n, we run an algorithm 1000 times independently, and record the average number of fitness evaluations until the optimal solution (for Jump) or the Pareto front (for OneJumpZeroJump) is found. We can observe from Figure 1 that using the proposed method can bring a clear acceleration.

8 Conclusion

This paper gives a theoretically study for EAs solving MMOPs, a class of optimization problems with wide applications. Considering the characteristic of MMOPs that multiple (local) optimal solutions in the solution space correspond to a single point in the objective space, we propose to prefer those solutions diverse in the solution space when they are equally good in the objective space. We show that such a method, working with crossover, can benefit the evolutionary search via rigorous running time analysis. Specifically, we prove that for (µ + 1)-GA solving the widely studied single-objective problem, Jump, as well as NSGA-II and SMS-EMOA solving the widely studied bi-objective problem, OneJumpZeroJump, the proposed method can lead to polynomial or even exponential acceleration on the expected running time, which is also verified by the experiments. Note that diversity has been theoretically shown to be helpful for evolutionary optimization in various scenarios [Friedrich et al., 2009; Qian et al., 2013; Doerr et al., 2016; Osuna et al., 2020; Sudholt, 2020; Doerr et al., 2023; Qian et al., 2024]. This work theoretically shows the benefit of maintaining the diversity in the solution space for EAs solving MMOPs, thus contributing to this line of research. We hope our results can be beneficial for the design of better EAs for solving MMOPs, and may also encourage the diversity exploration in the solution space for multi-objective EAs, which is often ignored.

Acknowledgments

This work was supported by the National Science and Technology Major Project (2022ZD0116600) and National Science Foundation of China (62276124). Chao Qian is the corresponding author. The supplementary is available at arXiv.