The First Proven Performance Guarantees for the Non-Dominated Sorting Genetic Algorithm II (NSGA-II) on a Combinatorial Optimization Problem

1 Introduction

Many optimization problems consist of several conflicting objectives. In such a situation, it is not possible to compute a single optimal solution. The most common solution concept therefore is to compute a set of Pareto optima (solutions which cannot be improved in one objective without accepting a worsening in another objective) and let a decision maker select the final solution based on their preference.

Besides mathematical programming approaches, evolutionary algorithms (EAs) are the standard approach to multiobjective problems with many successful applications [Zhou et al., 2011]. EAs profit here from their general ability to work with sets of solutions ("populations"). The by far most prominent multi-objective evolutionary algorithm (MOEA) is the Non-Dominated Sorting Genetic Algorithm II (NSGA-II) proposed by Deb, Pratap, Agarwal, and Meyarivan [Deb et al., 2002] (with over 50,000 citations on Google Scholar).

While very successful in practice, this algorithm is only little understood from a fundamental perspective, giving the users little general advice on how to optimally employ this algorithm, e.g., how to set its parameters right. In fact, it was only at AAAI 2022 that the first mathematical runtime analysis of the NSGA-II was presented [Zheng et al., 2022], a work that was quickly followed up more runtime analyses on this algorithm (see the previous works section). All these works analyze the performance of the NSGA-II on simple benchmark problems, mostly multi-objective variants of the ONEMAX, LEADINGONES, and JUMP benchmarks wellstudied in the theory of single-objective randomized search heuristics [Neumann and Witt, 2010; Auger and Doerr, 2011; Jansen, 2013; Zhou et al., 2019; Doerr and Neumann, 2020].

In this work, we conduct the first mathematical runtime analysis of the NSGA-II on a classic combinatorial problem, namely the NP-complete bi-objective minimum spanning tree problem. In this problem, we are given an undirected graph with n vertices and m edges. In the basic single-objective version of the minimum spanning tree (MST) problem, we are also given non-negative integral edge weights and the task is to compute a minimum spanning tree. This problem is easily solved by classic algorithms. It has also been used to understand how EAs solve combinatorial optimization problems. As a first result in this direction, Neumann and Wegener [Neumann and Wegener, 2007] showed that the (1 + 1) EA computes a minimum spanning tree in an expected number of O(m2 log(nwmax)) iterations (and fitness evaluations). Here wmax denotes the maximum edge weight. Using a balanced mutation operator, this can be improved to O(mn log(nwmax)).

In the bi-objective variant of the problem, we are given two weight functions and the target is to compute the Pareto front of the problem of computing a spanning tree minimizing both weight functions. This problem is NP-complete, but it is possible to compute in polynomial time the extremal points of the Pareto front [Hamacher and Ruhe, 1994]. The first result on how MOEAs solve this problem is [Neumann, 2007]. It shows that the global SEMO algorithm, a multiobjective analogue of the basic (1 + 1) EA, computes the extremal points of the Pareto front in an expected number of O(m2nwmin(|F| + log(nwmax))) iterations, where F denotes the set of extremal points of the Pareto front and wmin denotes the minimum of the maximum edge weights of the two weight functions. As in [Neumann and Wegener, 2007], this guarantee improves by a factor of Ω(m/n) when using balanced mutation.

We note that [Neumann and Witt, 2022] propose the convex global SEMO algorithm and show that it can solve the bi-objective MST problem in polynomial time regardless of wmax. Given that this algorithm is very new and not yet established, we do not follow this line of research. [Neumann and Witt, 2022, Theorem 4] also implies a runtime bound for the classic global SEMO algorithm, but this becomes superior to ours only when the size of the Pareto front is at most nwmin/ℓ, where ℓ is a problem parameter which in general can only be estimated by Θ(m2 ). So this result appears to be an improvement only in special cases.

In this work, we conduct a mathematical runtime analysis of the NSGA-II on the bi-objective MST problem. This is the first runtime analysis of this algorithm on a combinatorial problem. Besides showing that such analyses are possible, it proves that also the NSGA-II can efficiently compute the extremal points of the Pareto front, and this in an expected number of O(m2 log(nwmax)) iterations, hence a number of O(Nm2 log(nwmax)) fitness evaluations, when the population size N is at least N ≥ 4nwmin. This results holds for various ways of generating the offspring population including the use of crossover. As in the previous works, we obtain a bound lower by a factor of Ω(m/n) when using balanced mutation.

We note that our runtime guarantees are smaller than the ones proven in [Neumann, 2007] for the global SEMO by essentially a factor of |F|, the number of extremal points, for which the only general upper bound is nwmin. This improvement stems from our observation that the NSGA-II makes progress towards the different extremal points in parallel, whereas the proof in [Neumann, 2007] assumed that these were found sequentially. We show the same improvement for the GSEMO, lowering Neumann's bound to O(m2nwmin log(nwmax)).

2 Previous Work

The multi-objective version of the minimum spanning tree problem, usually called multi-criteria minimum spanning tree problem (mc-MST), is an important combinatorial optimization problem with many applications in network design. We refer to Ehrgott [Ehrgott, 2005] for an extensive discussion of the problem and the different algorithmic approaches to it. Being NP-complete, many heuristic approaches have been developed [Arroyo et al., 2008], including many based on evolutionary algorithms [Knowles and Corne, 2000; Knowles and Corne, 2001; Bossek and Grimme, 2017; Parraga-Alava et al., 2017; Majumder et al., 2020].

In this work, we investigate how the NSGA-II solves the bi-objective MST problem. As in most theoretical works on randomized search heuristics, our aim is not so much finding the best possible algorithm to solve this problem (for this problem-specific algorithms will usually be superior), but we aim at understanding how a certain algorithm, here the NSGA-II, solves a certain problem. The broader aim is to understand which algorithms are suitable for which problems, what are the right parameter settings, and to detect possible short-comings and remedies for these. We refer to the textbooks [Auger and Doerr, 2011] for a broader introduction into the research field of mathematical analyses of randomized search heuristics and its achievements.

We briefly review the most relevant literature. The mathematical runtime analysis of EAs was started in the 1990s with analysis how very simple EAs such as the (1 + 1) EA optimize very simple benchmark problems such as ONEMAX or LEADINGONES [M¨uhlenbein, 1992; B¨ack, 1993; Rudolph, 1997; Droste et al., 2002]. The first runtime analyses of MOEAs followed a similar approach, estimating the runtime of multi-objective analogues of the (1+1) EA such as the SEMO or GSEMO on bi-objective analogues of ONEMAX and LEADINGONES [Laumanns et al., 2002; Giel, 2003; Thierens, 2003]. The analyses of single-objective EAs quickly progressed towards more complex problems such as shortest paths, maximum matchings, the partition problem, the MST problem, and many others [Neumann and Witt, 2010], or more complex algorithms such as the (1+λ) EA [Jansen et al., 2005], (µ+λ) EA [Witt, 2006], (µ+λ) EA [Antipov and Doerr, 2021], (1+ (λ, λ)) GA [Doerr et al., 2015], and non-elitist algorithms [Dang et al., 2021].

In contrast, due to the more difficult population dynamics of MOEAs, the progress was slower in multi-objective evolutionary computation (there is, however, a highly successful line of research on multi-objectivization, that is, solving a single-objective problem via MOEAs [Neumann and Wegener, 2006; Friedrich et al., 2010; Qian et al., 2019; Crawford, 2019]; this line of work is not very related to the research conducted in this work, though). There are only few mathematical results on simple MOEAs solving combinatorial optimization problems (to the best of our knowledge only for the MST [Neumann and Wegener, 2007; Neumann and Witt, 2022], shortest path problems [Horoba, 2009; Neumann and Theile, 2010], and the travelling salesman problem [Lai and Zhou, 2020]). Also, only relatively few results analyze more complex algorithms such as the (µ + 1) SIBEA [Brockhoff et al., 2008; Nguyen et al., 2015; Doerr et al., 2016], the MOEA/D [Li et al., 2016; Huang et al., 2019; Huang and Zhou, 2020], the (1 + (λ, λ)) GSEMO [Doerr et al., 2022], the NSGA-II [Zheng et al., 2022; Zheng and Doerr, 2022a; Zheng and Doerr, 2022b; Bian and Qian, 2022; Doerr and Qu, 2023a; Doerr and Qu, 2023b; Doerr and Qu, 2023c; Dang et al., 2023b; Dang et al., 2023a], the NSGA-III [Wietheger and Doerr, 2023], and the SMS-EMOA [Bian et al., 2023]. Very roughly speaking, the works on the NSGA-II show that this algorithm can find the Pareto front of simple bi-objective benchmark problems when the population size is chosen large enough (typically by a constant factor larger than the Pareto front). When the population size is only equal to the size of the Pareto front or when the number of objective is more than two, already on the simple ONEMINMAX benchmark the NSGA-II needs exponential time to find the full Pareto front. So far, no mathematical runtime analysis of the NSGA-II on a combinatorial optimization problem exists.

3 Preliminaries: Basic Notation

The Multi-Objective Minimum Spanning Tree problem is stated as follows. Given an input connected graph G = (V,E), and a weight function w: E -→ Nd on the edges of G, define the weight of any subgraph H of G denoted w(H) ∈ Nd as the sum of the weights of all edges present in H. We want to find all possible "optimum" subtree weight values w ∈ Nd , in the sense that no subtree of G has a weight value for which all coordinates are smaller that those of w, and there exists a spanning tree that has weight w. Here, we focus on the case where d = 2. The search space is the set of all subgraphs and is represented with S = {0, 1}m as a subgraph is a choice of edges. Let w = (w1, w2) : E → N2 be the weight function. For a search point s, we refer to its weight as the geometric point of N2 denoted by ps = w(s). We further define

• wmax i =max{wi(e), e ∈ E}, for i ∈ {1, 2},
• wmax =maxi∈{1,2} wmax i
• wmin =mini∈{1,2} wmax i
• wub = n2wmax.

As in [Neumann, 2007], the fitness of an individual s ∈ S is given by a vector f(s) = (f1(s), f2(s)) with

for i ∈ {1, 2}, and where wi(j) is the weight of edge ej with respect to the function wi. c(s) is the number of connected components in the graph described by s, and e(s) is the number of edges in this same graph. Note that for a spanning tree s, fi(s) = ∑j|sj=1 wi(j).

Definition 1 (Domination). For s, s′ ∈ S, we say that s dominates s′ and we note s ⪯ s′ if s Pareto-dominates s′ according to the fitness functions f1 and f2, i.e. if f1(s) ≤ f1(s′ ) and f2(s) ≤ f2(s′ ). We say that s strictly dominates s′ and we note s ≺ s′ if s ⪯ s′ and s and s′ have different fitness values.

Since the relation of domination only depends on the objective value, we also use "domination" to compare objective values rather than individuals.

Definition 2 (Pareto optimality). s ∈ S is called Pareto optimal if there is no search point s′ ∈ S that strictly dominates s. The set of all Pareto optimal search points Spareto is called the Pareto set. F = f(Spareto) is the set of all Pareto optimal objective vectors and is called the Pareto front.

The goal is to find for each q ∈ F of the considered objective function f an object s ∈ Spareto with f(s) = q. From now on, as in [Neumann, 2007], we denote by conv(F) the lower-left part of the convex hull of F. The reader may find an illustration of conv(F ) in appendix of the Arvix version of the paper.

Definition 3 (Extremal points). The extremal points of the Pareto front F are the vertices of the polygonal line forming conv(F) .

Note that for each spanning tree T on the convex hull there is a λ ∈ [0, 1] such that T is a minimum spanning tree with respect to the single weight function λw1 + (1 - λ)w2 (see e.g. [Knowles and Corne, 2001], [Neumann, 2007]).

Let q1, q2, . . . , qr be the extremal points sorted in increasing f1 value. Observe that q1 (resp. qr) realises by construction the minimum of w1 (resp. w2) in w(S).

Those points are interesting because they give a solution which is a 2-approximation of the Pareto front ([Neumann, 2007]).

3.1 Algorithms

Here, we describe the NSGA-II algorithm. Let N be the population size and n the chromosome size. We also give pseudocode for the global SEMO (GSEMO), since part of our results directly come from the study of this algorithm in [Neumann and Wegener, 2007], and since most of these results also apply to this algorithm.

The NSGA-II

The NSGA-II is an evolutionary algorithm, which means it maintains a population of N ∈ N solutions to an optimization problem, named individuals (here, subgraphs of G). This population evolves over multiple generations. In each iteration, the algorithm generates an offspring population from selected parents in the current population, using some reproduction and mutation mechanisms. Note however that these mechanisms are not intrinsic to the NSGA-II and have to be specified for the task at hand. From the combined parent and offspring population, the algorithm decides which individuals to keep according to their fitness. The NSGA-II is centered around two notions. First, rank, which defines how good an individual is in the current population, and second, crowding distance, which quantifies how much diversity it brings to the population.

Definition 4 (Rank). Let X be a population. We recursively define the rank of every individual in X. An individual that is not strictly dominated by any individual in X has rank 1. Then, if ranks 1, . . . , k are defined, an individual that is not strictly dominated by any individual of X, except those of

Algorithm 1: The NSGA-II.

• 1 Generate initial population P ∈ ({0, 1}m)N

2 repeat

• 3 Generate offspring population Q ∈ ({0, 1}m)N
• 4 Let R = P ∪Q
• 5 Sort Rwith fast-non-dominated-sort to get the sets Fi, i ∈ N of the individuals of rank i
• 6 Find icut =max{i | ∑i-1 k=0 |Fk| < N}
• 7 Calculate crowding distance of each individual in Ficut
• 8 Let˜ Ficut be the N - ∑icut-1 k=0 |Fk| individuals in Ficut with largest crowding distance, chosen at random in case of a tie
• 9 P = (⋃icut-1 i=0 Fi) ∪˜ Ficut

10 until forever;

rank 1, . . . , k, and that is dominated by at least one individual of rank k, has rank k + 1.

The ranks of individuals in a given population X can be computed in quadratic time using the fast-non-dominatedsort algorithm. The reader will find in appendix of the Arvix version of the paper the associated pseudocode.

We now introduce the notion of crowding distance.

Definition 5 (Crowding distance). Let X be a finite set, and f : X → R. Let x1, . . . , xN be the elements of X, sorted in increasing order of f values. Then, the crowding distance of each xi is defined as

Note that in our case, we are dealing with multiple fitness functions (f1, f2). The crowding distance of an element is then the sum of the crowding distances of this element for each function. Computing the crowding distance of each individual in a population of size N can be done naively with two sortings and O(N) subtractions and divisions. It can then be done in O(N log(N)), which is negligible compared to fast-non-dominated sort.

The pseudocode for NSGA-II can be found in Algorithm 1.

The Global Simple Evolutionary Multiobjective Optimizer (GSEMO)

GSEMO is also an evolutionary algorithm, studied on the mc-MST problem in [Neumann, 2007] and [Neumann and Wegener, 2007]. Its functioning is simpler than that of the NSGA-II. It generates one individual s at each generation, adding it to the population if it is not dominated by any other individual, and removes those that are dominated by s. A pseudocode for GSEMO, as described in [Neumann, 2007] can be found in Algorithm 2.

4 Analysis of the GSEMO and the NSGA-II on the Bi-Objective Minimum Spanning Tree Problem

In this section, we prove two main results on the expected runtime of the GSEMO and the NSGA-II on the Bi-Objective

Algorithm 2: GSEMO

• 1 Generate initial population P, which consists of a unique individual s ∈ {0, 1}n , chosen randomly.

2 repeat

• 3 Choose a random s ∈ P
• 4 Generate an offspring s′ from x, flipping each bit with probability 1n .
• 5 if no individual in P dominates s′ then 6
• Add s′ to P 7 Remove all individuals in P that s′ dominates

8 until forever;

Minimum Spanning Tree Problem.

To state our theorem on the NSGA-II in the strongest possible form (which will enable us to easily obtain results for several different algorithm variants in Section 5), we need to introduce two parameters which depend on the offspring generation mechanism which the NSGA-II works with.

For any population P containing no spanning tree, any s ∈ F1 and any position i ∈ {1, . . . , n} of a bit of value 0, let p1 be a lower bound (that does not depend on P, s, i and j) of the probability that there exists in the offspring a child generated with s as (one of) the parent(s), that differs from s on exactly bit i.

For any population P containing at least a spanning tree, any s ∈ F1, and any pair of bits of different value i, j ∈ {1, . . . , n}, let p2 be a lower bound (that does not depend on P, s, i and j) of the probability that there exists in the offspring a child generated with s as (one of) the parent(s), that differs from s on exactly bits i, j.

Theorem 6. The expected number of generations until the NSGA-II, working on the fitness function f with a population of size N ≥ 4((n - 1)wmin + 1) and with any offspring generation mechanism resulting in p1 and p2, constructs a population which includes a spanning tree for each extremal vector of conv(F ) is upper bounded by O( log n p1 + log(nwmax) p2 ).

For the readers' convenience, we note already here that in typical versions of the NSGA-II, we have p1 = Θ(1/m) and p2 = Θ(1/m2 ), recall that the lenghth of the bit-string representation is m for MST problems. Then the runtime bound above becomes O(m2 log(nwmax) generations. We refer to Section 5 for the details.

Our arguments developed in the proof of Theorem 6 easily yield the following runtime estimate for the GSEMO, which improves the previously known results asymptotically.

Theorem 7. The expected number of fitness evaluations until the GSEMO working on the fitness function f constructs a population, which includes a spanning tree for each extremal vector of conv(F) is upper bounded by O(m2nwmin log(nwmax)).

The proof is split into three parts. First, we prove that the NSGA-II and the GSEMO find spanning trees in negligible time. Secondly, we introduce an elitism property of the NSGA-II from the moment when the population includes spanning trees. Then, we bound the time taken by both NSGA-II and GSEMO to find the extremal points from this point.

4.1 Sampling the First Spanning Tree

Regarding the first phase of the optimization process, the proof of [Neumann, 2007, Lemma 5] bounds the expected time until the population of the GSEMO contains at least one spanning tree from above, giving the following lemma.

Lemma 8. The GSEMO working on the fitness function f constructs a population with at least one spanning tree in expected time O(m log n).

We argue that a similar bound holds for the NSGA-II.

Lemma 9. The NSGA-II working on the fitness function f constructs a population with at least one spanning tree in expected time O( log n p1 ).

The proof is essentially an adaptation of [Neumann, 2007, Lemma 5], the details of which are in the appendix of the Arvix version of this paper. We first bound the time before the population contains a connected graph: to do so, we observe that for a subgraph H with l ≥ 2 connected components, there are at least l - 1 edges of G that decrease l. This allows us to use p1 to give a lower bound for the probability that these specific edges are added, and thus that connected components merge fast enough. A similar argument is used for the second phase, to show that excessive edges are deleted at roughly the same speed.

4.2 An ElitismProperty of the NSGA-II

One of the observations that make the study of the GSEMO easier is that it has an elitism property: an individual will not disappear from the population unless it is replaced by a dominating one. Here, we introduce a lemma, inspired by [Zheng and Doerr, 2022a], which shows that under a certain condition on the population size, the NSGA-II has a similar property.

Lemma 10. Let P be a population such that |P| > 4((n - 1)wmin +1) and having at least one spanning tree. Let P ′ be the next population. For each individual s in P, there is an individual s′ in P ′ such that s′ ⪯ s.

The proposed proof for this lemma requires to introduce the notion of incomparable set, which is used to bound the size of F1.

Definition 11 (Incomparable set). A ⊆ S is an incomparable set if there is no pair of individuals s, s′ ∈ A, such that s′ ≺ s.

We now make the following observation regarding the size of incomparable sets.

Lemma 12. Let P be an incomparable set containing at least one spanning tree. Then the number of objective values of P is bounded by (n - 1)wmin + 1.

Proof. Since f is designed such that a spanning tree strictly dominates every non-spanningtree graph, P consists of spanning trees only. Among two spanning trees of P with the same f1 value, one of them necessarily dominates the other. P being an incomparable set, this does only happens if the two spanning trees have the same objective value. Thus, for a fixed value of f1, there is at most one point in f(P) having this f1 value. Since a spanning tree has n-1 edges, and since by definition 0 ≤ w1(e) ≤ wmax 1 for any edge e, we deduce that, for any spanning tree s, 0 ≤ f1(s) ≤ (n - 1)wmax 1 . It follows that

where ST is the set of elements of S representing a spanning tree constructed from G. We do the same with f2, which gives |w(P)| ≤ (n - 1)wmin + 1, this bound being the minimum of the two bounds previously obtained.

The subsequent formal proof for Lemma 10 is rather long, and may be found in appendix of the Arvix version of the paper. However, we give a sketch of the principal arguments: we observe that, by definition of the fronts, the lemma requires to be proven only for individuals of the first front. Then, we show that, for a given fitness value p, there are at most 4 individuals that have fitness p and positive crowding distance. Finally, we use Lemma 12 to conclude that, because the NSGA-II will first choose individuals with positive crowding distance, the size of the population is such that the new generation will represent all fitness values of the first front.

4.3 Sampling the Extremal Points of the Pareto Front

This section gives upper bounds for the expected runtimes of the two algorithms, assuming the population contains at least a spanning tree, which completes the proof of Theorem 6 and 7.

We build our analysis on the fact that the extremal points are the unique minimums for specific linear combinations of the weights w1 and w2. The following lemma gives the combination for qi, when 1 ≤ i ≤ r. To state it more easily, we introduce the points q0 = q1 + (0, 1) and qr+1 = qr + (1, 0). These points are chosen such that for 1 ≤ i ≤ r, w(S) lies entirely on one of the two half planes delimited by the lines qjqj+1 and qj-1qj

Lemma 13. Let i ∈ {1, . . . , r}. For all objective values p, let

Then, for a given objective value p = w(s) of some individual s, we have di(p) ≥ 0 and di(p) = 0 ⇔ p = qi.

Proof. As shown in Figure 1, the set of all objective values is contained in the intersection of the upper-half planes Hj of the lines qjqj+1, 0 ≤ j ≤ r. For any such j, let

This vector is normal to the line qjqj+1. Also, w(S) lies entirely in one of the half planes delimited by qjqj+1, that vj is pointing towards. Thus, for all p ∈ R2 , p ∈ Hj if and only if p · vj ≥ q · vj for any point q lying on qjqj+1, with equality if and only if p ∈ (qj, qj+1). Now let us fix i ∈ {1, . . . , r}. We know that qi lies on both (qi, qi+1) and (qi-1, qi), so, for all p ∈ w(S), p · vi ≥ qi · vi and p · vi-1 ≥ qi · vi-1, with equality if and only if p ∈ (qi, qi+1) or p ∈ (qi-1, qi), respectively. By summing these two inequalities, we get (p - qi) · (vi + vi-1) ≥ 0, that is di(p) ≥ 0, with equality if and only if p ∈ (qi, qi+1)∩(qi, qi-1). Since qi and qi+1 are extremal points, these two lines cannot be parallel, hence (qi, qi+1) ∩ (qiqi-1) = {q}, which concludes the proof.

We suppose from now on that the population P contains a spanning tree, which is true after an expected number of O(m log n) iterations by Lemmas 8 and 9.

To prove the claimed upper bound, we use multiplicative drift analysis [Doerr et al., 2012] on a quantity derived from these functions. A reminder for this technique may be found in appendix of the Arvix version of the paper.

To use the multiplicative drift theorem, we introduce a potential that has a multiplicative drift. For a given population P, let

Lemma 14. Let P(t) be the population of the NSGA-II or the GSEMO algorithm after t iterations. Then di(P(t)) for any 1 ≤ i ≤ r, and d(P (t)) are nonnegative integers and are non-increasing with respect to t.

Proof. Let pi(t) = argminp∈w(P(t)) di(p). Using Lemma 10 in the case of the NSGA-II, and by the definition in the case of the GSEMO, there exists an element p′i ∈ w(P(t + 1)) such that p′i ⪯ pi. Hence,

so di(P(t)) is non-increasing. Lemma 13 gives that di(P(t)) is nonnegative. Summing di for 1 ≤ i ≤ r gives the same two properties for d(P(t)).

In order to apply drift analysis, we need another lemma given by [Neumann and Wegener, 2007].

Lemma 15 ([Neumann and Wegener, 2007], Lemma 2). Let ˜w: S -→ R be any function which is a linear combination with nonnegative coefficients of w1 and w2. Let s be a search point describing a spanning tree T . Let ˜wopt be the minimum value taken by ˜w on spanning trees. Then there exists a set of n 2-bit flips resulting on new spanning trees, such that the average weight decrease of these flips is at least ˜w(s)-˜wopt n

We now prove the following lemma.

Lemma 16. The expected number of generations until the NSGA-II (resp. GSEMO), working on the fitness function f with a population of size N ≥ 4((n - 1)wmin + 1) and with any offspring generation mechanism (resp. its intrinsic mechanism), constructs a population P such that d(P) = 0 is upper bounded by O( log(nwmax) p2 ) (resp. O(m2wmin log(nwmax))).

Proof. For any i ∈ {1, . . . , r} and any iteration t, let si,t be an individual such that di(si,t) = di(P(t)). Then, let p2,(i,t) be a lower bound over all pairs of bit positions of the probability that there exists, in the offspring, a child generated with s(i,t) as (one of) the parent(s) that differs from s(i,t) on exactly this pair of bits. For the NSGA-II algorithm, by definition of p2, we have p2,(i,t) ≥ p2.

For the GSEMO algorithm, if there is at least one spanning tree in the population, the size of the population is upperbounded by (n - 1)wmin + 1 by Lemma 12. Therefore,

1

We also denote p2 = (n-1)wminem2 when analysing GSEMO, such that we have p2,(i,t) ≥ p2 for both algorithms. We now apply drift analysis on d(P(t)). Indeed, for i ∈ {1, . . . , r}, di(p) is a linear combination of w1 and w2 with non-negative coefficients minimum dopt i = 0.

Focusing only on the set of n 2-bits flips given by Lemma 15 on a search point being one of argmins∈P di(s), and noticing that, since di is non-increasing (Lemma 14), all other bit flips contribute positively to the drift, we have

Using this estimate, we show the following lower bound for the drift of our potential.

Let T = min{t ∈ N | d(P(t)) = 0}. Since d(P) is an non-negative integer, the minimum strictly positive value of d(P) is 1. Since each edge has a weight at most 2w2max , d(P) ≤ ∑ri=1 2mw2max ≤ 2rn2w2max . Finally, since the convex hull is an incomparable set, we have r ≤ (n - 1)wmin + 1 by Lemma 12. Hence we obtain

Using multiplicative drift theorem, we obtain

We then conclude the proof by plugging in the value of p2 for GSEMO, which does only one fitness evaluation per iteration.

Combining the last lemma with Lemma 13, and noticing that NSGA-II and GSEMO do respectively N and 1 fitness evaluation(s) per generation, we conclude the proof of theorems 6 and 7.

5 Study of Offspring GenerationMechanisms of the NSGA-II

This section is dedicated to demonstrating applications of Theorem 6 on the study of the influence of particular offspring generation mechanisms on the expected runtime of the NSGA-II before finding the extremal points.

We start with a very rudimentary mechanism, composed of fair selection (selecting every individual s in P as a parent), and standard bitwise mutation with constant c, 0 ≤ c ≤ m (mutating every bit of a given parent with probability cm ). Then we can define the lower bounds on the success probabilities as follows.

Theorem 6 then gives the following corollary.

Corollary 17. The expected time until the NSGA-II, working on the fitness function s with a population of size N ≥ 4((n1)wmin + 1) and using standard bitwise mutation and fair selection, constructs a population which includes a spanning tree for each extremal vector of conv(F) is upper bounded by O(m2 log(nwmax)) generations, and

fitness evaluations.

With note that the same result would hold when replacing bitwise mutation with fast mutation [Doerr et al., 2017], which proved to be advantageous in some multi-objective problems [Doerr and Zheng, 2021].

Theorem 6 also ensures that adding crossover with constant probability 0 < q < 1 does not worsen the asymptotic complexity of the NSGA-II. To be more precise, let us consider a randommonoparental reproduction schemeM, which, from an individual s ∈ S, generates an offspring M(s). From this operator, and a given crossover operator C we derive a biparental reproduction scheme M′ , which, from a pair of individuals s, s′ , generates, with probability q, the children C(s, s′ ) and C(s′ , s), and with probability (1 - q), the children M(s), and M(s′ ). Now, consider an offspring generation scheme which puts the selected parents into pairs and appliesM′ to these pairs. Finally, note that if p1, p2 and p′1 , p′2 are the pairs of probabilities defined in 6 for the offspring generation mechanism using M or M′ , then we have trivially p′1 ≥ (1 - q)p1 and p′2 ≥ (1 - q)p2. This gives the following corollary.

Corollary 18. For a given NSGA-II implementation, replacing M by M′ , as described above (i.e. adding crossover with constant probability) does not worsen the asymptotic expected runtime (up to a multiplicative constant).

Finally, to demonstrate the relevance of the p2 parameter, let us introduce a problem-specific "balanced" mutation operator, described in [Neumann and Wegener, 2007]. For any tree individual s ∈ S the mutation operator flips each 1-bit with probability 1 n-1 , and each 0-bit with probability 1 m-n+1 , and applies standard bitwise mutation otherwise. We use fair selection for this example. By definition, the value of p1 remains unchanged compared to standard bitwise mutation, that is, O( 1m ). However, one can show that p2 is increased from O( 1 m2 ) to O( 1 (m-n)n ) (see appendix of the Arvix version of the paper). Theorem 6 automatically gives:

Corollary 19. The expected runtime of the NSGA-II, working on the fitness function f, with |P| > 4((n - 1)wmin + 1), using fair selection and balanced mutation is O((m - n)n log(nwmax) + m log n) generations and

fitness evaluations.

6 Conclusion

In this first mathematical runtime analysis of the NSGA-II on a combinatorial optimization problem, we provided a general approach to proving runtime guarantees for MOEAs solving the bi-objective MST problem. For the global SEMO, this gave a bounds lower than the previously known ones by a factor of Ω(|F|). More interestingly, we could prove the same performance guarantees for the much more complex NSGA-II. Our result applies to several variants of the NSGA-II, including some that use crossover. As for the simple global SEMO algorithm, we obtain better guarantees when employing a balanced mutation operator (which supports the general belief that analyses on simple toy algorithms can nevertheless give useful hint for the use of more complex algorithms).

Overall, this work indicate that mathematical runtime analyses for the NSGA-II are possible also for combinatorial optimization problems. In this first work, we mostly concentrated on proving performance guarantees at all. For future work, it would be interesting to derive more insights on how to optimally use the NSGA-II on the particular problem (we only saw that balanced mutation is preferable). That such results are possible in principle is again indicated by the simpler works on artificial benchmarks, where, e.g., [Doerr and Qu, 2023a] gave some indications on the right mutation rate. Clearly, runtime analyses of the NSGA-II on other combinatorial optimization problems would also be desirable to put this research direction on a broader basis.

Acknowledgments

This work was supported by a public grant as part of the Investissements d'avenir project, reference ANR-11-LABX0056-LMH, LabEx LMH and a fellowship via the International Exchange Program of ´Ecole Polytechnique.