Quality-Diversity Algorithms Can Provably Be Helpful for Optimization

1 Introduction

Search algorithms, e.g., simulated annealing, Evolutionary Algorithms (EAs), and Bayesian optimization, are traditionally to find a single highest-performing solution to the given problem. However, for many real-world tasks, e.g., in reinforcement learning [Eysenbach et al., 2018; Chalumeau et al., 2023], robotics [Cully et al., 2015; Salehi et al., 2022], and human-AI coordination [Lupu et al., 2021; Strouse et al., 2021], there is need for finding a set of high-performing solutions that exhibit varied behaviors. A behavior space can be characterized by dimensions of variation of interest to a user. For example, if searching for robot gaits that can run fast, the dimensions of interest can be the portion of time that each leg is in contact with the ground. Illuminating high-performing solutions of each area of the entire behavior space reveals relationships between dimensions of interest and performance, which can help improve the robustness in complex environments, e.g., help robots rapidly adapt to damage (i.e., discover a behavior compensating for the damage) [Cully et al., 2015].

Quality-Diversity (QD) algorithms [Pugh et al., 2016; Cully and Demiris, 2018] are a subset of EAs, and have emerged as a potent optimization paradigm for the challenging task of finding a set of high-quality and diverse solutions. They follow a typical evolutionary process, which iteratively select some parent solutions from the population (also called archive in the QD area), apply variation operators (e.g., crossover and mutation) to produce offspring solutions, and use these offspring solutions to update the population. Lehman and Stanley [2011b] proposed the QD algorithm, Novelty Search and Local Competition (NSLC), by maximizing a local competition objective (e.g., the number of nearest neighbors of a solution worse than itself) and a novelty objective (e.g., the average distance of a solution to its nearest neighbors) simultaneously. The Multi-Objective Landscape Exploration (MOLE) algorithm [Clune et al., 2013] is similar to NSLC, except that a global performance objective, instead of the local competition objective, is considered.

Mouret and Clune [2015] proposed the Multi-dimensional Archive of Phenotypic Elites (MAP-Elites) algorithm to achieve the goal of finding a set of high-quality and diverse solutions more straightforwardly. MAP-Elites maintains an archive by discretizing the behavior space into cells and storing at most one solution in each cell. The generated solutions compete only when their behaviors belong to the same cell, and MAP-Elites tries to fill in the cells with as high-quality solutions as possible. MAP-Elites performs better than NSLC and MOLE empirically [Mouret and Clune, 2015], and has become the most popular QD algorithm. Many advanced variants have been proposed to improve its sample efficiency, e.g., refining parent selection strategies [Cully and Demiris, 2018; Sfikas et al., 2021; Wang et al., 2022; Wang et al., 2023] and variation operators [Nilsson and Cully, 2021; Tjanaka et al., 2022; Faldor et al., 2023], improving archive structure [Vassiliades et al., 2018; Grillotti and Cully, 2022], using surrogate models [Lim et al., 2022; Bhatt et al., 2022], and applying cooperative coevolution [Xue et al., 2024].

They have achieved encouraging performance in various reinforcement learning tasks, e.g., exploration [Conti et al., 2018; Ecoffet et al., 2021], robust training [Kumar et al., 2020], and environment generation [Fontaine et al., 2021; Zhang et al., 2023].

Besides generating a diverse set of high-performing solutions, extensive experiments (e.g., on evolving locomoting virtual creatures [Lehman and Stanley, 2011b], producing modular neural networks [Mouret and Clune, 2015], designing soft robots [Mouret and Clune, 2015], encouraging exploration of agents in deceptive [Conti et al., 2018] and sparse [Colas et al., 2020] reward environments, and solving traveling thief problems [Nikfarjam et al., 2022a]) have shown that QD algorithms can find a better overall solution than traditional search algorithms which focus on finding the single highest-performing solution. That is, if only looking at the best achieved objective value, QD algorithms are still better. Even when searching for behavioral novelty only while ignoring the objective, the empirical results on maze navigation and biped walking show that it significantly outperforms objective-based search [Lehman and Stanley, 2011a]. The better optimization ability of QD algorithms may be due to the better exploration of the search space, which helps avoid local optima and thus find better objective values.

Compared to their wide application, the theoretical analysis of QD algorithms is far behind. To the best of our knowledge, there have existed only two pieces of works on runtime analysis of QD algorithms. Nikfarjam et al. [2022b] analyzed the MAP-Elites algorithm for solving the knapsack problem. They proved that by using carefully designed twodimensional behavior spaces, MAP-Elites can simulate dynamic programming behaviors to find an optimal solution within expected pseudo-polynomial time. Bossek and Sudholt [2023] further proved that for maximizing monotone submodular functions with a size constraint, MAP-Elites using the number of items contained by a set as the behavior descriptor can achieve a (1 - 1/e)-approximation ratio in expected polynomial time; by using the number of connected components as the behavior descriptor, MAP-Elites can solve the minimum spanning tree problem in expected polynomial time. They also analyzed MAP-Elites for solving any pseudoBoolean problem f : {0, 1}n → R. The number of 1-bits of a solution is used as the behavior descriptor, and the i-th cell of the behavior space stores the best found solution with the number of 1-bits belonging to [(i - 1)k, ik - 1], where k (dividing n + 1) is a granularity parameter. They proved the upper bound on the expected runtime of MAP-Elites until covering all cells. For the specific pseudo-Boolean problem OneMax which aims to maximize the number of 1-bits of a solution, they proved that the derived bound is tight, even for the time to find the optimal solution within each cell. There are also some preliminary studies trying to understand the behavior of novelty search from a theoretical perspective. Doncieux et al. [2019] hypothesized that novelty search asymptotically behaves like uniform sampling in the behavior space, and Wiegand [2024] showed that novelty search will converge in the archive space by utilizing k nearest neighbor theory.

In this paper, we try to theoretically examine the fundamental question: whether QD algorithms can be helpful for optimization, i.e., finding better objective values? For this purpose, we compare the popular QD algorithm MAP-Elites with (µ + 1)-EA via rigorous runtime analysis. (µ + 1)-EA maintains a population of µ solutions, and focuses on finding better objective values by iteratively generating an offspring solution via mutation and using it to replace the worst solution in the population if it is better. For fair comparison, the population size µ of (µ + 1)-EA is set to the number of cells of the archive of MAP-Elites, and all their parameters are set to the same. The only difference is that for MAP-Elites, the generated offspring solution only competes with the solution in the same cell of the archive, while for (µ + 1)-EA, the offspring competes with all solutions in the population, and replaces the worst one if it is better. Thus, their performance gap can be purely attributed to the QD algorithms' framework of searching for a set of high-quality and diverse solutions.

We consider two NP-hard problem classes with wide applications. For monotone approximately submodular maximization with a size constraint [Das and Kempe, 2011], we prove that MAP-Elites achieves the optimal polynomialtime approximation ratio of 1 - e-γ [Harshaw et al., 2019], where γ ∈ [0, 1] measures the closeness of a set function f to submodularity (Theorem 1). When f is submodular, γ = 1, and the approximation ratio of MAP-Elites is 1 - 1/e (Corollary 1); while there exists an instance where (µ + 1)-EA requires expected exponential time to achieve a nearly 1/2-approximation ratio (Theorem 2). For the set cover problem with m elements [Feige, 1998], MAP-Elites achieves the asymptotically optimal approximation ratio of lnm + 1 in expected polynomial time (Theorem 3), while (µ + 1)-EA requires expected exponential time to achieve an approximation ratio smaller than 2m+1 /m on some instance (Theorem 4). These results clearly show the advantage of MAP-Elites over (µ + 1)-EA, providing theoretical justification for that QD algorithms can be helpful for optimization. The analyses disclose that the simultaneous search of QD algorithms for high-performing solutions at diverse regions in the behavior space plays a key role, which provides stepping stones to high-performing solutions in a region that may be difficult to be discovered if searching only in that region (which may get trapped in local optima). Due to space limitation, we only give the proof sketch of the theorems in the paper, and the proof details are given in the supplementary.

2 Preliminaries

2.1 MAP-Elites Algorithm

MAP-Elites is the most popular QD algorithm [Mouret and Clune, 2015; Cully et al., 2015]. Given an objective function f : X → R to be maximized, MAP-Elites tries to find a set of high-quality (i.e., having large f values) and diverse solutions, instead of a single highest-performing solution (i.e., a solution having the largest f value). Specifically, a user chooses m dimensions of variation of interest that define a behavior space, and MAP-Elites aims to create a map of highperforming solutions at each point in the behavior space. The behavior space can be characterized by a behavior descriptor function b : X → Rm. For a solution x ∈ X, b(x) returns an m-dimensional vector describing x's behaviors.

Algorithm 1 MAP-Elites Algorithm

Input: an objective function f : X → R, and a behavior descriptor function b : X → Rm

Parameter: #initial solutions I

• 1: Initialize an empty archive M ← ∅, and let iter = 0;
• 2: repeat
• 3: if iter < I then
• 4: Choose x′ fromX uniformly at random
• 5: else
• 6: Choose x fromM uniformly at random;
• 7: Create x′ by mutating x
• 8: end if
• 9: if M(b(x′ )) = ∅ or f(x′ ) > f(M(b(x′ ))) then 10: M(b(x′ )) = x′
• 11: end if
• 12: iter = iter + 1
• 13: until some criterion is met
• 14: returnM

In our theoretical analysis, we use the simple default version of MAP-Elites [Mouret and Clune, 2015], as presented in Algorithm 1. For each behavior descriptor vector b(·), an archive M uses a cell (denoted as M(b(·))) to store the highest-performing solution found so far. Note that multiple solutions may have the same behavior descriptor vector. If no solution is found for a cell in the archiveM, the cell is empty. The archive M is initialized to be empty in line 1. The first I solutions are randomly generated (lines 3-4). After that, a solution x is randomly chosen from the archive M (line 6), and mutated to generate an offspring solution x′ (line 7) in each iteration. Each generated solution x′ will be used to update the corresponding cell M(b(x′ )), as shown in lines 9-11. If the cell is empty (i.e., M(b(x′ )) = ∅), x′ will be directly placed there. If x′ is higher-performing than the current occupant of the cell (i.e., f(x′ ) > f(M(b(x′ )))), the occupant will be replaced by x′ . Note that according to user preference or available computational resources, multiple behavior descriptor vectors may be merged into one cell in practice.

2.2 (µ+1)-EA

(µ+1)-EA is a typical EA focusing on finding solutions with larger f values, which has been widely used in theoretical analyses of EAs, e.g., [Jansen and Wegener, 2001; Witt, 2006; Storch, 2008]. As presented in Algorithm 2, (µ+1)-EA maintains a set of µ solutions, i.e., employs a population size of µ. It first initializes a population P by randomly generating µ solutions in line 1. Then, in each iteration, it generates one new solution x′ by mutating a solution x randomly chosen from the population P (lines 3-4), and uses x′ to replace the worst solution in P if x′ is better (lines 5-8).

2.3 Algorithm Settings &Performance Measure

To theoretically show that QD algorithms can be helpful for optimization (i.e., finding a single, best solution), we will compare MAP-Elites with (µ+1)-EA. As the problems studied in Sections 3 and 4 have the solution space X = {0, 1}n (where n denotes the problem size), we employ the commonly used bit-wise mutation operator in Algorithms 1 and 2,

Algorithm 2 (µ + 1)-EA

Input: an objective function f : X →R Parameter: population size µ

• 1: Initialize a population P by choosing µ solutions fromX uniformly at random;
• 2: repeat
• 3: Choose x from P uniformly at random;
• 4: Create x′ by mutating x;
• 5: Let y = argminy∈P f(y); ties are broken randomly;
• 6: if f(x′ ) > f(y) then
• 7: P ← (P \ {y}) ∪ {x′ }
• 8: end if
• 9: until some criterion is met
• 10: return P

which flips each bit of a solution x with probability 1/n to generate a new solution x′ . For fairness, the population size µ of (µ + 1)-EA is set to the number of cells in the archive M; as (µ+1)-EA generates µ initial solutions randomly, the number of randomly generated initial solutions (i.e., the parameter I) of MAP-Elites is set equal to µ. These settings are to make that the performance gap between MAP-Elites and (µ+1)-EA is purely due to the mechanism of finding a set of high-quality and diverse solutions of MAP-Elites.

To be consistent with the original setting of MAP-Elites in [Mouret and Clune, 2015], we have used strict survivor selection, i.e., f(x′ ) > f(M(b(x′ ))) in line 9 of Algorithm 1 and f(x′ ) > f(y) in line 6 of Algorithm 2. However, our analyses will not be affected if allowing accepting the same good solutions, i.e., replacing > with ≥. We consider maximization problems in Algorithms 1 and 2. For minimization, we only need slight modification, i.e., replacing > with < in both algorithms, and changing y = argminy∈P f(y) in line 5 of Algorithm 2 to y = arg maxy∈P f(y).

To measure the performance of QD algorithms, the three metrics, i.e., optimization, coverage and QD-Score, are often used [Pugh et al., 2016; Cully and Demiris, 2018]. Optimization corresponds to the largest objective value of solutions in the archive, i.e., the objective value of the highest-performing solution found by the algorithm. Coverage is the number of non-empty cells, i.e., the total number of solutions, in the archive. QD-Score is equal to the total sum of objective values across all solutions in the archive, which reflects both the quality and diversity of the found solutions.

Here, we focus on the optimization metric, serving for our purpose of showing that QD algorithms can be helpful for optimization, i.e., finding better objective values. Specifically, we will analyze the expected runtime complexity of MAPElites or (µ + 1)-EA, which counts the number of objective evaluations (the most time-consuming step in the evolutionary process) until finding an α-approximate solution, and has been a leading theoretical aspect for randomized search heuristics [Neumann and Witt, 2010; Auger and Doerr, 2011; Zhou et al., 2019; Doerr and Neumann, 2020]. Let OPT denote the optimal function value. An α-approximate solution x means f(x) ≥ α · OPT (where α ≤ 1) for maximization, and f (x) ≤ α · OPT (where α ≥ 1) for minimization.

3 Monotone Approximately Submodular Maximization with a Size Constraint

Let R and R+ denote the set of reals and non-negative reals, respectively. Given a finite set V = {v1, . . . , vn}, a set function f : 2V → R is monotone if ∀X ⊆ Y, f(X) ≤ f(Y ). We assume w.l.o.g. that monotone functions are normalized, i.e., f(∅) = 0. A set function f is submodular [Nemhauser et al., 1978] if ∀X ⊆ Y ⊆ V, f(Y ) - f(X) ≤ ∑v∈Y \X (f (X ∪ {v}) - f(X)). The submodularity ratio in Definition 1 measures how close a general set function f is to submodularity. When f is monotone, it holds: ∀X, l : 0 ≤ γX,l(f) ≤ 1; f is submodular iff ∀X, l : γX,l(f) = 1.

Definition 1 (Submodularity Ratio [Das and Kempe, 2011]). Let f : 2V → R be a set function. The submodularity ratio of f with respect to a set X ⊆ V and a parameter l ≥ 1 is

As presented in Definition 2, the studied problem is to find a subset with at most k items such that a given monotone and approximately submodular objective function is maximized. It is NP-hard, and has many applications, such as maximum coverage [Feige, 1998], influence maximization [Kempe et al., 2003], sensor placement [Krause et al., 2008], sparse regression [Das and Kempe, 2011], and unsupervised feature selection [Feng et al., 2019], just to name a few.

Definition 2 (Monotone Approximately Submodular Function Maximization with a Size Constraint). Given a monotone and approximately submodular function f : 2V → R+ and a budget k, the goal is to find a subset X ⊆ V such that

In the following analysis, we equivalently reformulate the above problem as an unconstrained maximization problem by setting f(X) = -1 when |X| > k, i.e., X does not satisfy the size constraint. Note that a subset X of V can be naturally represented by a Boolean vector x ∈ {0, 1}n , where the i-th bit xi = 1 if vi ∈ X, and xi = 0 otherwise. Throughout the paper, we will not distinguish x ∈ {0, 1}n and its corresponding subset {vi ∈ V | xi = 1} for convenience.

3.1 Analysis of MAP-Elites

To apply MAP-Elites in Algorithm 1 to maximize monotone approximately submodular functions with a size constraint, we use the number of 1-bits (denoted as |x|, which equals ∑ni=1 xi) of a solution x ∈ {0, 1}n as the behavior descriptor. Thus, the archive M maintained by MAP-Elites contains n + 1 cells, where for any 0 ≤ i ≤ n, the i-th cell stores the best found solution with i 1-bits. In each iteration of MAPElites, if an offspring solution x′ with j 1-bits is generated, the j-th cell in the archive M will be examined: if the cell is empty or the existing solution in the cell is worse than x′ , the j-th cell will then be occupied by x′ .

We prove in Theorem 1 that after running O(n2 (log n+k)) expected time, MAP-Elites achieves an approximation ratio of 1 - e-γmin , which reaches the optimal polynomial-time approximation ratio [Harshaw et al., 2019]. Inspired from the analysis of GSEMO (a multi-objective EA with similar procedure to MAP-Elites) [Friedrich and Neumann, 2015; Qian et al., 2015b; Qian et al., 2019], the proof can be accomplished by following the behavior of the greedy algorithm [Das and Kempe, 2011]. Note that the parameter I = n+1 implies that MAP-Elites generates n+1 initial solutions randomly. The value of I does not affect the analysis, and such a setting is just for fairness, because the population size µ of (µ + 1)-EA will be set to the number of cells (i.e., n + 1) contained by the archive M.

Theorem 1. For maximizing a monotone approximately submodular function f with a size constraint k, the expected runtime of MAP-Elites with the parameter I = n+1, until finding a solution x with |x| = k and f(x) ≥ (1 - e-γmin ) · OPT, is O(n2 (log n + k)), where γmin = minx:|x|=k-1 γx,k, and γx,k is the submodularity ratio of f w.r.t. x and k.

Proof Sketch. We divide the optimization process of MAPElites into two phases: (1) starts after initialization and finishes after finding the special solution 0 (i.e., the empty set); (2) starts after phase (1) and finishes after finding a solution with the desired approximation guarantee.

For phase (1), we consider the minimum number of 1-bits of the solutions in the archive M, denoted as Jmin. We show that Jmin will not increase, but can decrease in each iteration by selecting the solution x corresponding to Jmin and flipping only one 1-bit in lines 6-7 of Algorithm 1, which happens with probability (1/|M|) · (Jmin/n)(1 - 1/n)n-1 ≥ Jmin/(en(n + 1)), where |M| ≤ n + 1. Thus, the expected number of iterations required by phase (1) (i.e., making Jmin = 0) is at most ∑ni=1 en(n + 1)/i = O(n2 log n).

For phase (2), we consider Jmax, which is defined as the maximum value of j ∈ {0, 1, . . . , k} such that in the archive M, there is a solution x with |x| = j and f(x) ≥ (1 - (1 - γmin/k)j )·OPT. We only need to analyze the expected number of iterations until Jmax = k, which implies that ∃x ∈ M with |x| = k and f(x) ≥ (1 - (1 - γmin/k)k ) · OPT ≥ (1 - e-γmin ) · OPT, i.e., the desired approximation guarantee is reached. We show that Jmax will not decrease, but can increase in each iteration by selecting the solution x corresponding to Jmax and flipping only one specific 0-bit (according to Lemma 1 of [Qian et al., 2016]), which happens with probability (1/|M|) · (1/n)(1 - 1/n)n-1 ≥ 1/(en(n + 1)). Thus, the expected number of iterations required by phase (2) (i.e., making Jmax = k) is at most k · en(n + 1).

Summing up the expected time of the two phases leads to the theorem. Note that the additional evaluations for the n+1 initial solutions will not affect the derived time bound. □

For submodular f, ∀x : γx,k = 1, and γmin = 1. Thus, Corollary 1. For maximizing a monotone submodular function f with a size constraint k, the expected runtime of MAPElites with I = n + 1, until finding a solution x with |x| = k and f (x) ≥ (1 - 1/e) · OPT, is O(n2 (log n + k)).

Corollary 1 is consistent with Theorem 5.1 in [Bossek and Sudholt, 2023]. Note that the MAP-Elites algorithm they analyzed is slightly different from ours, which generates only one random initial solution (i.e., the parameter I = 1) and replaces the solution in a cell if the generated offspring belonging to the same cell is not worse (i.e., f(x′ ) > f(M(b(x′ )))

in line 9 of Algorithm 1 changes to f(x′ ) ≥ f(M(b(x′ )))). Thus, these slight differences do not affect the result.

3.2 Analysis of (µ+1)-EA

To apply (µ + 1)-EA in Algorithm 2 to maximize monotone approximately submodular functions with a size constraint, we set the population size µ to the number of cells contained by the archive M of MAP-Elites, which is n+1. As we have mentioned in Section 2.3, such a setting is for fairness, which makes the population size of (µ + 1)-EA equal to the largest archive size of MAP-Elites. In each iteration of (µ + 1)-EA, the generated offspring solution x′ will be used to replace the worst solution in the current population if x′ is better.

Next, we will give a concrete instance of the classic NPhard problem, maximum coverage [Feige, 1998] (which is an application of monotone submodular maximization with a size constraint), where (µ + 1)-EA fails to achieve a good approximation ratio in polynomial time. Given a family V of sets that cover a universe U of elements, the goal of maximum coverage as presented in Definition 3 is to select at most k sets from V to make their union maximal. The objective function f(x) = | ⋃i:xi=1 Si| is monotone and submodular.

Definition 3 (Maximum Coverage). Given a ground set U, a collection V = {S1, S2, . . . , Sn} of subsets of U, and a budget k, the goal is to find a subset x of V such that

The considered instance of maximum coverage as presented in Example 1 is from [Friedrich et al., 2010; Friedrich and Neumann, 2015], and can be obtained from the complete bipartite graph in Figure 1(a). The ground set U is given by the set of edges. Each vertex vi corresponds to one subset Si of U, which consists of the edges adjacent to the vertex vi. The goal is to select at most k = (1 + δ)n/3 subsets of U to make the number of contained edges maximal.

Example 1. The ground set U contains all the edges of the complete bipartite graph in Figure 1(a). Each set Si in V consists of the edges adjacent to the vertex vi, i.e., ∀1 ≤ i ≤ (1 + δ)n/3, Si = {(vi, v(1+δ)n/3+1), . . . , (vi, vn)}; ∀(1 + δ)n/3 + 1 ≤ i ≤ n, Si = {(vi, v1), . . . , (vi, v(1+δ)n/3)}, where (vi, vj) denotes the edge between the vertexes vi and vj. The budget k = (1 + δ)n/3, and δ is a small positive constant close to 0.

The two parts of the complete bipartite graph in Figure 1(a) contain (1 + δ)n/3 and (2 - δ)n/3 vertexes, respectively. As the budget k = (1 + δ)n/3 < (2 - δ)n/3, the unique optimal solution is x∗ = 1(1+δ)n/3 0(2-δ)n/3 (i.e., {S1, . . . , S(1+δ)n/3}), which covers all the edges, and has the objective value f(x∗ ) = (1 + δ)(2 - δ)n2 /9. If a solution

i.e., contains (1 + δ)n/3 sets from {S(1+δ)n/3+1, . . . , Sn}, it is local optimal, with f(xlocal) = (1 + δ)2n2 /9. Adding a set from {S1, . . . , S(1+δ)n/3} into xlocal will cover (2 - δ)n/3 - (1 + δ)n/3 = (1 - 2δ)n/3 more edges. If i sets from {S1, . . . , S(1+δ)n/3} have been added into xlocal, then

𝑣!

𝑣"

𝑣#$&!

𝑣#$

deleting one set from xlocal will decrease the number of edges by (1+δ)n/3-i. To make (1+δ)n/3-i ≤ (1-2δ)n/3, i should be at least δn. Thus, to improve the local optimal solution xlocal, at least δn sets from {S1, . . . , S(1+δ)n/3} need to be added into xlocal, and meanwhile at least δn sets from xlocal have to be deleted. The approximation ratio of xlocal is f (xlocal)/f (x∗ ) = (1 + δ)/(2 - δ).

We prove in Theorem 2 that (µ + 1)-EA requires at least exponential (w.r.t. n) expected runtime to achieve an approximation ratio larger than (1 + δ)/(2 - δ) on Example 1. The main proof idea is to show that there is some probability that the population of (µ + 1)-EA will get trapped in a local optimal solution xlocal, escaping from which requires at least 2δn bits to be flipped simultaneously.

Theorem 2. There exists an instance of monotone submodular maximization with a size constraint (i.e., the maximum coverage instance in Example 1), where the expected runtime of (µ+1)-EA with µ = n+1 for achieving an approximation ratio larger than (1+δ)/(2-δ) is at least exponential w.r.t. n, where δ is a small positive constant close to 0.

Proof Sketch. (µ+ 1)-EA generates n+ 1 initial solutions by sampling from {0, 1}n uniformly at random. We first consider the event A that among the generated n + 1 initial solutions, one is a local optimal solution xlocal, and the other n solutions do not satisfy the size constraint, i.e., contain more than k = (1 + δ)n/3 1-bits. It happens with probability P (A) ≥ (n + 1)(1 - o(1))/2n by the Chernoff bound.

Given the event A, we consider that in the first n iterations of (µ + 1)-EA, a local optimal solution is always selected from the population in line 3 of Algorithm 2, and no bits are flipped by the bit-wise mutation operator in line 4. Such an event (denoted as B) will make the population full of local optimal solutions, because the copied local optimal solution will replace one solution violating the size constraint (whose objective value is set to -1) in the population after each iteration. It happens with probability P(B) ≥ √2πn/(e(2e2 )n ).

To improve xlocal, it needs to flip at least δn of its (1+δ)n/3 leading 0-bits and δn of its (1+δ)n/3 1-bits, with probability at most ((1+δ)n/3 δn )((1+δ)n/3 δn ) ( 1n )2δn ≤ 1 2πδn ( e(1+δ) 3δn )2δn . This implies that the expected runtime to generate a solution better than xlocal is at least 2πδn(3δn/(e(1+ δ)))2δn .

Combining the above analysis and using the fact that the approximation ratio of xlocal is (1 + δ)/(2 - δ), we derive that the expected runtime of (µ + 1)-EA for achieving an approximation ratio larger than (1 + δ)/(2 - δ) is at least P(A) · P(B) · 2πδn(3δn/(e(1 + δ)))2δn ≥ ((n + 1)(1 - o(1))/2n ) · (√2πn/(e(2e2 )n )) · 2πδn(3δn/(e(1 + δ)))2δn , which is exponential w.r.t. n because δ is a small positive constant close to 0. Thus, the theorem holds. □

The comparison between Corollary 1 and Theorem 2 clearly shows the advantage of MAP-Elites over (µ + 1)- EA. That is, for monotone submodular maximization with a size constraint, MAP-Elites can always achieve the (1-1/e)- approximation ratio in expected polynomial time, while there exists an instance where (µ+1)-EA requires at least expected exponential time to achieve an approximation ratio of nearly 1/2. From the proofs, we can find that MAP-Elites explicitly maintains multiple cells of the behavior space, allowing it to achieve a good solution with k 1-bits by following a path across different cells (i.e., gradually adding the number of 1bits), while (µ + 1)-EA only focuses on the objective value, and may lose diversity in the population and get trapped in local optima (with k 1-bits). This discloses the benefit of simultaneously searching for high-performing solutions at diverse regions in the behavior space, which may provide stepping stones to high-performing solutions in a region that may be difficult to be discovered if searching only in that region.

4 Set Cover

Set cover [Feige, 1998] is a well-known NP-hard combinatorial optimization problem. Given a family V of sets (where each set has a corresponding positive weight) that cover a universe U of elements, the goal is to find a subset of V with the minimum weight such that all the elements of U are covered.

Definition 4 (Set Cover). Given a ground set U = {e1, . . . , em}, and a collection V = {S1, S2, . . . , Sn} of subsets of U with corresponding weights w : V → R+, the goal is to find a subset of V (represented by x ∈ {0, 1}n ) such that

Let wmax = max{wi | 1 ≤ i ≤ n} and wmin = min{wi | 1 ≤ i ≤ n} denote the maximum and minimum weights of a single set Si, respectively. In the following analysis, we equivalently reformulate the above problem as an unconstrained minimization problem by setting

where w(x) = ∑ni=1 wixi denotes the weight of x, c(x) = | ⋃i:xi=1 Si| denotes the number of elements covered by x, and λ > nwmax is a sufficiently large penalty. That is, a solution covering more elements is better; if covering the same number of elements, a solution with a smaller weight is better.

4.1 Analysis of MAP-Elites

To apply MAP-Elites in Algorithm 1 to solve the set cover problem, we use the number c(x) = | ⋃i:xi=1 Si| of covered elements of a solution x ∈ {0, 1}n as the behavior descriptor. Thus, the archive M maintained by MAP-Elites contains m+ 1 cells, where for any 0 ≤ i ≤ m, the i-th cell stores the best found solution covering i elements. In each iteration of MAPElites, if an offspring solution x′ with c(x′ ) = j is generated, the j-th cell inM will be examined: if the cell is empty or the existing solution in the cell is worse than x′ (i.e., has a larger weight than x′ ), the j-th cell will then be occupied by x′ .

We prove in Theorem 3 that after running O(mn(m + logn + log(wmax/wmin))) expected time, MAP-Elites achieves an approximation ratio of lnm+ 1, which has been optimal up to a constant factor, unless P=NP [Feige, 1998]. The proof is accomplished by following the behavior of the greedy algorithm [Chvatal, 1979], inspired from the analysis of GSEMO [Friedrich et al., 2010; Qian et al., 2015a]. The parameter I of MAP-Elites is set to m+ 1 (i.e., it generates m + 1 initial solutions randomly) for fairness, because the population size µ of (µ + 1)-EA will be set to the number of cells (i.e., m+ 1) contained by the archive M.

Theorem 3. For the set cover problem, the expected runtime of MAP-Elites with the parameter I = m+ 1, until finding a solution x with c(x) = m and f(x) ≤ (lnm + 1) · OPT, is O(mn(m + logn + log(wmax/wmin))), where c(x) = | ⋃i:xi=1 Si| denotes the number of elements covered by x.

Proof Sketch. We divide the optimization process of MAPElites into two phases. Phase (1) starts after initialization and finishes after finding the special solution 0 which does not cover any element. We consider the minimum weight of solutions in the archive M after running t iterations, denoted as Xt = min{w(x) | x ∈ M}. Note that Xt = 0 implies that the solution 0 has been found. By showing that the expected change of Xt after one iteration satisfies E(Xt - Xt+1 | Xt) ≥ Xt/(en(m+ 1)) and applying the multiplicative drift analysis theorem [Doerr et al., 2012], we derive that the expected runtime to make Xt = 0 (i.e., find the solution 0) is at most O(mn(log n + log(wmax/wmin))).

After that, phase (2) aims to find a solution with the desired approximation guarantee. Let Hi = ∑ij=1 1/j denote the ith harmonic number, and Rk = Hm - Hm-k. We consider Kmax, which is the maximum value of k such that ∃x ∈ M with c(x) = k and w(x) ≤ Rk ·OPT. Then, we only need to analyze the expected runtime until Kmax = m, implying having found a solution with c(x) = m (i.e., covering all the elements) and w(x) ≤ Rm ·OPT < (ln m+1) ·OPT. We show that Kmax will not decrease, but can increase by flipping a specific 0 bit (i.e., adding the set Si with the minimum ratio of its weight and the number of elements it newly covered) of the solution x corresponding to Kmax, which happens with probability (1/|M|)·(1/n)(1-1/n)n-1 ≥ 1/(en(m+1)), where |M| ≤m+1. Since it is sufficient to increase Kmax at most m times for making Kmax = m, the expected runtime of phase (2) is O(m2n). By combining the two phases, the total expected runtime is O(mn(m + log n + log(wmax/wmin))). Thus, the theorem holds. □

4.2 Analysis of (µ+1)-EA

Next, we will give a concrete instance of set cover, where (µ + 1)-EA fails to achieve a good approximation ratio in polynomial time. To apply (µ + 1)-EA in Algorithm 2 to solve the set cover problem, the population size µ is set to the number of cells (i.e., m+ 1) contained by the archive M of MAP-Elites for fair comparison. The considered instance is obtained from the complete bipartite graph in Figure 1(b).

Example 2. The ground set U contains all the n - 1 edges of the complete bipartite graph in Figure 1(b). Each set Si in V consists of the edges adjacent to the vertex vi, i.e., S1 = {(v1, v2), . . . , (v1, vn)}; ∀2 ≤ i ≤ n, Si = {(vi, v1)}, where (vi, vj) denotes the edge between the vertexes vi and vj . Their weights are: w1 = 2n ; ∀2 ≤ i ≤ n, wi = 1.

In this example, the size m of the ground set U is n - 1, wmax = 2n , and wmin = 1. The unique optimal solution is x∗ = 01n-1 (i.e., {S2, . . . , Sn}), which covers all the edges, and has the minimum weight w(x∗ ) = n - 1. We can also observe that xlocal = 10n-1 (i.e., containing only S1) is local optimal, with the weight w(xlocal) = 2n . Furthermore, xlocal has the smallest weight except x∗ , implying that it needs to flip all the n bits to improve xlocal. Its approximation ratio is w(xlocal)/w(x∗ ) = 2n /(n - 1) = 2m+1 /m. Note that if a solution does not cover all the edges, i.e., does not satisfy the constraint, its objective value will be penalized, as shown by the reformulated equivalent objective function in Eq. (1).

We prove in Theorem 4 that (µ + 1)-EA requires at least exponential (w.r.t. m, n, and log(wmax/wmin)) expected time to achieve an approximation ratio smaller than 2m+1 /m on Example 2. The proof idea is similar to that of Theorem 2, i.e., to show that there is some probability that the population will get trapped in the local optimum xlocal, escaping from which requires all the n bits to be flipped simultaneously.

Theorem 4. There is a set cover instance (i.e., Example 2), where the expected runtime of (µ + 1)-EA with µ = m + 1 for achieving an approximation ratio smaller than 2m+1 /m is at least exponential w.r.t. m, n, and log(wmax/wmin).

Proof Sketch. As the population size µ = m + 1 = n, (µ+1)-EA generates n initial solutions randomly by uniform sampling from {0, 1}n . We consider the event A that one initial solution is the local optimal solution xlocal, and the other n-1 initial solutions are all worse than xlocal, which happens with probability P (A) ≥ n(1-o(1))/2n . Given the event A, we then consider that in the first n-1 iterations of (µ+1)-EA, xlocal is always selected from the population, and no bits are flipped by bit-wise mutation. This event (denoted as B) will make the population consist of n copies of xlocal, and happens with probability P (B) ≥ √2π(n - 1)/(e(2e2 )n-1 ).

After that, it needs to flip all the n bits of xlocal to make an improvement, whose probability is 1/nn . Because the approximation ratio of xlocal is 2n /(n - 1), the expected runtime of (µ + 1)-EA for achieving an approximation ratio smaller than 2n /(n - 1) is at least P(A) · P(B) · nn ≥ (n(1-o(1))/2n ) · (√2π(n - 1)/(e(2e2 )n-1 )) · nn , which is exponential w.r.t. n. As m = n - 1 and log(wmax/wmin) = log(2n /1) = n in Example 2, the theorem holds. □

Theorems 3 and 4 show that MAP-Elites is superior to (µ+ 1)-EA for set cover, because MAP-Elites can always achieve the (lnm + 1)-approximation ratio in expected polynomial time, while (µ + 1)-EA requires expected exponential time to achieve an approximation ratio smaller than 2m+1 /m on some instance. As found in the last section, the simultaneous search for high-performing solutions at diverse regions in the behavior space makes MAP-Elites better. More specifically, MAP-Elites obtains a good solution covering all elements by following a path across different cells of the behavior space (i.e., gradually covering more elements), while (µ + 1)-EA may lose diversity in the population and get trapped in local optima in the space of solutions covering all elements.

5 Conclusion and Discussion

This paper contributes to providing theoretical justification for the benefit of QD algorithms, i.e., bringing better optimization. On the two NP-hard problems, monotone approximately submodular maximization with a size constraint, and set cover, we prove that MAP-Elites achieves the (asymptotically) optimal polynomial-time approximation ratio, while (µ + 1)-EA requires exponential expected time on some instances. The simultaneous search of QD algorithms for highperforming solutions with diverse behaviors brings better exploration of the search space, helping avoid local optima and find better solutions. In fact, we can find more theoretical evidences by comparing the results from previous separate works on QD and EA's runtime analysis, though not completely fair. For example, Bossek and Sudholt [2023] proved that for maximizing any monotone function over {0, 1}n , MAP-Elites that generates one initial solution and uses the number of 1-bits of a solution as the behavior descriptor can find an optimal solution in O(n2 logn) expected time, while Lengler and Zou [2021] proved that (µ + 1)-EA with µ0 ≤ µ ≤ n (where µ0 is some constant) needs superpolynomial time to optimize some monotone functions. In the future, we will theoretically verify the other benefit of QD algorithms that they can better create a map of high-performing solutions in the behavior space than separately searching for a high-performing solution in each cell of the behavior space, by analyzing the QD-Score metric, i.e., the sum of objective values across all the obtained solutions in the archive.

Many theoretical works are left to be done, e.g., comparing different QD frameworks (MAP-Elites vs. NSLC), and studying the influence of different components (e.g., parent selection and variation operators), which may inspire the design of better QD algorithms in practice. Another interesting work is to study the relationship between MAP-Elites and the corresponding multi-objective EA, GSEMO, as also indicated in [Bossek and Sudholt, 2023]. By treating the behavior descriptors as the extra objective functions to be optimized, GSEMO behaves somewhat similarly to MAP-Elites. The differences are: 1) MAP-Elites only compares solutions within the same cell (i.e., having the same behavior descriptor values) while GSEMO compares different behavior descriptor values based on the domination relationship; 2) MAPElites can control the granularity of the behavior space by setting the number of behavior descriptor values in a cell. The known results of MAP-Elites (including [Bossek and Sudholt, 2023] and ours) match that of GSEMO in [Neumann and Wegener, 2006; Friedrich et al., 2010; Friedrich and Neumann, 2015; Qian et al., 2019]. Thus, it would be interesting to show their performance gap on some problems, especially that MAP-Elites may overcome some difficulty of GSEMO due to its good diversity mechanism.

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 material is available at https://arxiv.org/abs/2401.10539.