A First Runtime Analysis of NSGA-III on a Many-Objective Multimodal Problem: Provable Exponential Speedup via Stochastic Population Update

1 Introduction

Evolutionary multi-objective algorithms (EMOAs) use principles of nature to optimize functions with multiple conflicting objectives by finding a Pareto optimal set [Emmerich and Deutz, 2018]. These have been commonly applied to a wide range of optimization problems in practice [Zapotecas-Mart´ınez et al., 2023; Keller, 2017; Gunantara, 2018] which also include neural networks [Liu et al., 2020], bioinformatics [Handl et al., 2008], engineering [Sharma and Chahar, 2022] or various fields of artificial intelligence [Luukkonen et al., 2023; Zhang et al., 2021; Monteiro and Reynoso-Meza, 2023]. As it is typical for real world problems [Stewart et al., 2008], such problems of EMOAs on many-objective problems has quickly gained huge importance in many research fields. However, when the number of objectives increases, the Pareto front grows exponentially in the number of objectives and hence problems typically become more challenging. Additionally, it becomes more difficult to identify dependencies between single objectives. NSGA-II [Deb et al., 2002], the most prominent EMOA (∼53.000 citations), is able to optimize biobjective problems efficiently (see [Zheng et al., 2022] for a first rigorous analysis and for example [Dang et al., 2023b; Doerr and Qu, 2022] for further rigorous ones or [Deb et al., 2002] for empirical results), but loses performance if the number of objectives grows (see [Zheng and Doerr, 2024a] for rigorous results where there is a huge difference already between two and three objectives or [Khare et al., 2003; Purshouse and Fleming, 2007] for empirical studies). This behaviour comes from the crowding distance, the tie breaker in NSGA-II, which is based on sorting search points in each objective. In case of two objectives, a sorting of nondominated individuals with respect to the first objective induces a sorting with respect to the second one (in reverse order) and hence the crowding distance is a good measure for the closeness of an individual to its neighbors in the objective space. However, for problems with three or more objectives, such a correlation between different objectives does not necessarily exist and hence individuals may have crowding distance zero even if they are not close to others in the objective space. To overcome this problem, the NSGA-III algorithm was designed in [Deb and Jain, 2014] which uses a set of predefined reference points instead of the crowding distance. It can optimize a broad class of different benchmarks with at least four objectives efficiently [Yannibelli et al., 2020; Gu et al., 2022], demonstrating its success in practice (∼6.000 citations). However, theoretical understanding of its success is still in its early development and lags far behind its practical impact. To the best of our knowledge, there are only a few papers which address rigorous runtime analyses of this algorithm [Wietheger and Doerr, 2023; Wietheger and Doerr, 2024; Opris et al., 2024; Opris, 2025] and only [Wietheger and Doerr, 2024] investigated a multimodal problem like the d-OJZJ problem, where the algorithm has to cross a fitness valley to cover the entire Pareto front. But especially on functions with local optima it is important to understand the mechanics of NSGA-III and its variants to obtain valuable insights when and why it performs well since such many-objective problems occur very often in real world scenarios [Geng et al., 2023]. This may help practitioners to design improved versions of NSGA-III with enhanced performance.

In this paper we significantly extend the results from [Wietheger and Doerr, 2024] and provide a runtime analysis of NSGA-III on d-OJZJ with and without stochastic population update. Our three main contributions are detailed as follows. Firstly, we show an expected runtime bound on NSGA-III on d-OJZJ of O(nk+d/2 /µ + n ln(n)) generations, or O(nk+d/2+n ln(n)µ) fitness evaluations, for a constant number d of objectives and population size µ ≥ (2n/d + 1)d/2 , but µ ∈ 2O(n) and 2 ≤ k ≤ n/(2d). This is by a factor of min{µ/nd/2 , nk-1 / ln(n)} smaller than the bound of O(nk ) computed in [Wietheger and Doerr, 2024]. Remarkably, for µ ∈ ω(nd/2 )∩o(nk-1 / ln(n)), our bound is o(nk µ), and hence if d = 2, NSGA-III outperforms NSGA-II if also µ ∈ o(n2 /k2 ) (compare with [Doerr and Qu, 2023a] for the tight runtime bound of Θ(µnk ) of NSGA-II on OJZJ if n + 1 ≤ µ ∈ o(n2 /k2 )). We also derive a lower runtime bound of NSGA-III on OJZJ of Ω(nk+1 /µ) generations, or Ω(nk+1 ) fitness evaluations in expectation, for 2 ≤ k ≤ n/4 and µ ∈ O(nk-1 ) when the number of objectives is d = 2. This implies tightness if also n + 1 ≤ µ = O(nk / ln(n)). This is another novel approach, since, as far as we know, the only rigorously proven lower runtime bounds for NSGA-III are provided in [Opris, 2025] on an artificial benchmark quite different to d-OJZJ.

Secondly, to be able to derive these runtime bounds, we have to investigate the population dynamics of NSGA-III more carefully than previous works did. To this end, we use the cover number of an objective vector v, which is the number of individuals in the current population with fitness vector v. Outgoing from some general results applicable to any fitness function, we prove with probability 1 - e-Ω(n) that after O(n) generations, when there are only Pareto optimal individuals, the cover number of every Pareto optimal vector is at most cµ/nd/2 for a suitable constant c. This also holds in all future iterations if this property about Pareto optimality is not violated. Hence, all solutions are spread out quite evenly on a large fraction of the Pareto front. For example, if µ = Θ(nd/2 ), the cover number of every v is then only at most constant.

Thirdly, we investigate NSGA-III with stochastic population update where the next generation is not only formed by selecting the first ranked solutions resulting from non-dominated sorting, i.e. in a greedy, deterministic manner, but also by selecting some individuals chosen uniformly at random [Bian et al., 2023]. Hence, low-ranked, but promising solutions also have a certain chance to survive. This feature may lead to an exponential speedup in the runtime. This version of NSGA-III can optimize OJZJ in O(k(bn/k)k ) generations, where b > 1 is a suitable constant, extending results from [Bian et al., 2023] for NSGA-II and d = 2 to NSGA-III and many objectives (see also [Zheng and Doerr, 2024b] for similar results regarding the SMS-EMOA). This is by at least a factor of Ω((k/b)k-1nd/2 /µ) smaller compared to the upper bound above which is exponentially large if k is linear in n and µ = Θ(nd/2 ).

To achieve these, we also have to adapt the arguments from [Opris et al., 2024] about the protection of good solutions also to the case when stochastic population update is considered.

There are several theoretical runtime analyses showcasing the efficiency of NSGA-II on bi-objective problems. The first was conducted by [Zheng et al., 2022] on classical benchmark problems, followed by results on a multimodal problem [Doerr and Qu, 2022], about the usefulness of crossover [Dang et al., 2024a; Doerr and Qu, 2023b], noisy environments [Dang et al., 2023a], approximations of covering the Pareto front [Zheng and Doerr, 2022], lower bounds [Doerr and Qu, 2023a], trap functions [Dang et al., 2024b] and stochastic population update [Bian et al., 2023]. There are also results on combinatorial optimization problems like the minimum spanning tree problem [Cerf et al., 2023] or the subset selection problem [Deng et al., 2024]. However, rigorous runtime results in many-objective optimization on simple benchmark functions appeared only recently for the SMS-EMOA [Zheng and Doerr, 2024b], the SPEA2 [Ren et al., 2024], variants of the NSGA-II [Doerr et al., 2025] and the NSGA-III [Wietheger and Doerr, 2023; Opris et al., 2024]. For the NSGA-III, in [Wietheger and Doerr, 2023] the first runtime analysis on the 3-ONEMINMAX problem with p ≥ 21n divisions along each objective was conducted for defining the set of reference points. In [Opris et al., 2024] this result was generalized to more than three objectives and runtime analyses for the classical d-COUNTINGONESCOUNTINGZEROES and dLEADINGONESTRAILINGZEROES benchmarks [Laumanns et al., 2004] for any constant number d of objectives were provided where it is also necessary to reach the Pareto front. They could also reduce the number of required divisions by more than half. Finally, the first runtime analysis of NSGA-III on d-OJZJ is given in [Wietheger and Doerr, 2024].

2 Preliminaries

For a finite set A we denote by |A| its cardinality and by ln the logarithm to base e. For n ∈ N let [n] := {1, . . . , n}. The number of ones in a bit string x is denoted by |x|1 and the number of zeros by |x|0, respectively. For two random variables Y and Z on N0 we say that Z stochastically dominates Y if P(Z ≤ c) ≤ P(Y ≤ c) for every c ≥ 0. For a d-objective function f : {0, 1}n → Nd0 , x ↦→ (f1(x), . . . , fd(x)), let fmax := max{fj(x) | x ∈ {0, 1}n , j ∈ [d]} be the maximum possible objective value. When d = 2, f is also called bi-objective. For two search points x, y ∈ {0, 1}n , x weakly dominates y, written as x ⪰ y, if fi(x) ≥ fi(y) for all i ∈ [d] and x (strictly) dominates y, written as x ≻ y, if one inequality is strict. We call x and y incomparable if neither x ⪰ y nor y ⪰ x. A set I ⊆ {0, 1}n is a set of mutually incomparable solutions if all search points in I are incomparable. Each solution x not dominated by any other in {0, 1}n is called Pareto-optimal and we call f(x) non-dominated fitness value. The set of all non-dominated fitness values is called Pareto front. For a

1 Initialize P0 ∼ Unif(({0, 1}n )µ ) 2 for t := 0 to∞do 3 Initialize Qt := ∅ 4 for i = 1 to µ do 5 Sample s from Pt uniformly at random 6 Create r by standard bit mutation on s with mutation probability 1/n 7 Update Qt := Qt ∪ {r} 8 Set Rt := Pt ∪Qt 9 if a = 1 then 10 Update Rt by choosing ⌈3µ/2⌉ solutions from Rt uniformly at randomwithout replacement 11 Partition Rt into layers F 1 t , F 2 t , . . . , F k t of non-dominated solutions 12 if a = 0 then 13 Find i∗ ≥ 1 such that ∑i∗-1 i=1 |F i t | < µ and ∑i∗ i=1 |F i t | ≥ µ 14 else 15 Find i∗ ≥ 1 such that ∑i∗-1 i=1 |F i t | < ⌈µ/2⌉ and ∑i∗ i=1 |F i t | ≥ ⌈µ/2⌉ 16 Compute Yt = ⋃i∗-1 i=1 F i t 17 Choose ˜ F i∗ t ⊂ F i∗ t such that |Yt ∪ ˜ F i∗ t | = µ if a = 0 and |Yt ∪ ˜ F i∗ t | = ⌈µ/2⌉ if a = 1 with Algorithm 2 18 Create the next population Pt+1 := Yt ∪ ˜ F i∗ t if a = 0 and Pt+1 := Yt ∪ ˜ F i∗ t ∪ ((Pt ∪Qt) \ Rt) if a = 1

population Pt the cover number ct(v) of v ∈ Nd0 is the number of individuals x ∈ Pt with f(x) = v and we say that v is covered if its cover number is at least 1.

The NSGA-III algorithm, originated in [Deb and Jain, 2014], with or without stochastic population update is shown in Algorithm 1. Initially, a population of size µ is created by choosing µ individuals from {0, 1}n uniformly at random. Then in each iteration t, a multiset Qt of µ new offspring is created by µ times choosing an individual s ∈ Pt uniformly at random and applying standard bit mutation on s, i.e. each bit is flipped independently with probability 1/n.

During the survival selection, the parent and offspring populations Pt and Qt are merged into Rt. When stochastic population update is turned on (i.e. a = 1) then Rt is updated by choosing ⌈3µ/2⌉ individuals from Rt uniformly at random without replacement. Then Rt is partitioned into layers F 1 t+1 , F 2 t+1 , . . . using the non-dominated sorting algorithm [Deb et al., 2002] where F 1 t+1 consists of all nondominated individuals, and F i t+1 for i > 1 of individuals only dominated by those from F 1 t+1 , . . . , F i-1 t+1 . Then if a = 0 the critical rank i∗ with ∑i∗-1 i=1 |F t i | < µ and ∑i∗ i=1 |F t i | ≥ µ

Algorithm 2: Selection procedure utilizing a set Rp of reference points to maximize a function

• 1 Compute the normalisation f n of f
• 2 Associate each x ∈ Yt ∪ F i∗ t with its reference point rp(x) such that the distance between f n (x) and the line through the origin and rp(x) is minimized
• 3 For each r ∈ Rp, set ρr := |{x ∈ Yt | rp(x) = r}|
• 4 Initialize ˜ F i∗ t = ∅ and R′ :=Rp
• 5 while true do
• 6 Determine rmin ∈ R′ such that ρrmin is minimal (where ties are broken randomly)
• 7 Determine xrmin ∈ F i∗ t \ ˜ F i∗ t which is associated with rmin and minimizes the distance between the vectors f n (xrmin ) and rmin (where ties are broken randomly)

is determined (i.e. there are fewer than µ search points in Rt with a lower rank than i∗ , but at least µ search points with rank at most i∗ ). If a = 1, then every individual not selected for Rt survives. The number of such individuals is µ - ⌈3µ/2⌉ = ⌊µ/2⌋. Hence, in this case, the critical index i∗ refers only on µ - ⌊µ/2⌋ = ⌈µ/2⌉ individuals. All individuals with a lower rank than i∗ are included in Pt+1, while the remaining individuals are selected from F t i∗ using Algorithm 2. Hereby, a normalized objective function f n is computed and then each individual with rank at most i∗ is associated with reference points. For the first, we use the normalization procedure from [Wietheger and Doerr, 2023] which can be also used for maximization problems as shown in [Opris et al., 2024]. We omit detailed explanations as they are not needed for our purposes. For a d-objective function f : {0, 1}n → Nd0 , the normalized fitness vector f n (x) := (f n 1 (x), . . . , f n d (x)) of a search point x is computed as

for each j ∈ [d] where ynad := (ynad 1 , . . . , ynad d ) and ymin := (ymin 1 , . . . , ymin d ) from the objective space are called nadir and ideal points, respectively. Computing the nadir point is not trivial and we have ynad j ≥ εnad, and ymin j ≤ ynad j ≤ ymax j for every j ∈ [d] where εnad is a positive threshold set by the user (see [Blank et al., 2019] or [Wietheger and Doerr, 2023] for the details). Further, ymax j and ymin j are the maximum and minimum value in objective j from all search points seen so far (i.e. from P0,Q0, . . . , Pt,Qt). After computing the normalisation, each individual x is associated with the reference point rp(x) such that the distance between f n (x) and the line through the origin and rp(x) is minimal (see also Figure 1). Ties are broken deterministically, which means that two or more individuals which have the same smallest distance to two or more reference points are associated to the same reference point. We use the same set of reference points Rp as proposed in [Deb and Jain, 2014]. The points are defined as

where p ∈ N is a parameter one can choose according to the fitness function f. These are uniformly distributed on the simplex determined by the unit vectors (1, 0, . . . , 0)⊺ , (0, 1, . . . , 0)⊺ , . . . , (0, 0, . . . , 1)⊺ .

Then, one iterates through all the reference points where the reference point with the fewest associated individuals that are already selected for the next generation Pt+1 is chosen. A reference point is omitted if it only has associated individuals that are already selected for Pt+1 and ties are broken uniformly at random. Next, from the individuals associated to that reference point who have not yet been selected, the one closest to the chosen reference point is selected for the next generation, where ties are again broken uniformly at random. Once the required number of individuals is reached (i.e. if |Yt| + | ˜ F i∗ t | = µ when a = 0 and |Yt| + | ˜ F i∗ t | = ⌈µ/2⌉ otherwise) the selection ends.

3 Structural Results

The following result from [Opris et al., 2024] also adapted to the case when stochastic population update is turned on (i.e. a = 1) shows that if a population covers a fitness vector v with a first-ranked individual x then it is covered for all future generations as long as f(x) is non-dominated. In other words, the best solutions are protected.

Lemma 1. Consider NSGA-III on a d-objective function f with εnad ≥ fmax and a set Rp of reference points for p ∈ N

with p ≥ 2d3/2 fmax. Let Pt be its current population. Assume µ ≥ |I| if a = 0 or µ ≥ 2|I| otherwise for the population size µ where I is a maximum set of mutually incomparable solutions. Let ˜ F 1 t , ˜ F 2 t , . . . be the layers of non-dominated fitness vectors of Pt ∪Qt. Then for every x ∈ ˜ F 1 t there is an x′ ∈ Pt+1 with f(x′ ) = f(x).

Proof. The case a = 0 is Lemma 3.4 in [Opris et al., 2024]. So suppose that a = 1. Let x ∈ ˜ F 1 t . If x /∈ Rt then the claim holds (by Line 18 in Algorithm 1 we have (Pt ∪ Qt) \ Rt ⊂ Pt+1). Suppose that x ∈ Rt. Then x ∈ F 1 t where F 1 t is the first layer of non-dominated solutions from Line 11 in Algorithm 1 (since if x is non-dominated with respect to Pt ∪ Qt then it is also non-dominated with respect to the updated Rt from Line 10 in Algorithm 1). If |F 1 t | ≤ ⌈µ/2⌉ all individuals in F 1 t (including x) survive. Hence, suppose that |F 1 t | > ⌈µ/2⌉ and let s := |{f(x) | x ∈ F 1 t }| be the number of different fitness vectors of search points from F 1 t . Then by Lemma 3.3 in [Opris et al., 2024] individuals with distinct fitness vectors in F 1 t are associated with different reference points. Thus there are s reference points, each associated with at least one individual. Since ⌈µ/2⌉ ≥ |I| ≥ s, at least one individual x′ ∈ Rt associated with the same reference point as x survives and we have f(x) = f(x′ ).

The population dynamics of NSGA-III significantly differs to that of NSGA-II. In case of NSGA-II, which uses the crowding distance as second tie breaker, all rank one individuals with crowding distance zero are treated equally in the survival selection no matter how they are distributed. In contrast, NSGA-III (without stochastic population update) keeps uniform distributions of Pareto optimal individuals across the search space. This behaviour can be formally described using the cover number ct(v) of a fitness vector v, defined as the number of individuals in Pt with fitness v, captured in the following lemma.

Lemma 2. Assume the same conditions as in Lemma 1 where stochastic population update is disabled. The following properties hold.

1. Let v ∈ F ∗ and 0 ≤ a ≤ µ/|I|. If ct(v) ≥ a then ct+1(v) ≥ a.
2. Let v ∈ F ∗ and suppose that ct+1(v) < ct(v). Then ct+1(w) ≤ ct(v) for every w ∈ F ∗ .
3. Suppose that every x ∈ Pt is Pareto optimal. Then mt :=max{ct(v) | v ∈ F ∗ } does not increase.

Proof. (1): The NSGA-III iterates through all reference points, always preferring a reference point r with the fewest associated individuals chosen for Pt+1 so far (see Line 6 in Algorithm 2), and selecting an individual for Pt+1 associated to r. Since, by Lemma 3.3 in [Opris et al., 2024], two Pareto optimal search points with distinct fitness are associated to two different reference points, NSGA-III iterates at least a times through all reference points with at least a many associated individuals to find Pt+1. Hence, the cover number ct+1(v) of v with respect to Pt+1 is still at least a.

(2): There is a reference point r to that all x with f(x) = v are associated. Hence, NSGA-III iterates at most ct(v) - 1

times through r, since ct+1(v) < ct(v), and consequently, through all reference points at most ct(v) times (see Line 6 in Algorithm 2). Thus, for every w ∈ V , at most ct(v) individuals x with f (x) = w survive.

(3): By iterating through all reference points at most mt times, the algorithm finds µ individuals in F 1 t (since all Pareto optimal individuals are in F 1 t ) and the survival selection ends. Then ct+1(w) ≤ mt for every w ∈ F ∗ .

4 The Many-Objective Jump Function

In this section we define the d-ONEJUMPZEROJUMPk (dOJZJk) function, defined in [Zheng and Doerr, 2024b], as a d-objective version of the bi-objective JUMPk benchmark. The latter was introduced in [Doerr and Qu, 2022] to understand how MOEAs can handle functions with local optima, i.e. changes of size k ≥ 2 are necessary to cover the whole Pareto front. Fix d ∈ N divisible by two and let n be divisible by d/2. For a bit string x let x := (x1 , . . . , xd/2 ) where all xj are of equal length 2n/d. For odd j ∈ [d] define yj := x(j+1)/2 and yj := xj/2 for even j ∈ [d]. For 2 ≤ k ≤ 2n/d the d-OJZJk(x) = (f1(x), . . . , fd(x)) is defined as

if j ∈ [1, . . . , d] is even. Note that for odd j ∈ [d] the (j + 1)-th objective is structurally identical to the j-th one with the roles of ones and zeros reversed. For every objective there are 2n/d + 1 different values and the maximum possible value is k + 2n/d. In [Zheng and Doerr, 2024b] it is shown that the Pareto front F ∗ of d-OJZJk is {(a1, 2k + 2n/d - a1, . . . , ad/2, 2k + 2n/d - ad/2) | a1, . . . , ad/2 ∈ {k, 2k, 2k +1, . . . , 2n/d-1, 2n/d, 2n/d+k}} and has cardinality (2n/d-2k+3)d/2 for k ≤ n/d. Further, a maximum set of mutually incomparable solutions I (its cardinality can be larger than |F ∗ |, compare with [Zheng and Doerr, 2024b]) fulfills |F ∗ | ≤ |I| ≤ (2n/d + 1)d/2 since I does not contain two search points with the same fitness vector and hence |I| ≤ |f({0, 1}n )| = (2n/d + 1)d/2 . For d = 2 we have that every non-Pareto optimal individual is dominated by a Pareto optimal one which is different in the many-objective case.

Lemma 3. Suppose that d = 2 and let y be not Pareto optimal. Then every Pareto optimal x strictly dominates y. 1

To simplify the analysis, we assume that d is constant.

5 Population Dynamics of NSGA-III on OJZJ

In this section we investigate the population dynamics of the steady state NSGA-III on d-OJZJ and demonstrate that, in an expected linear number of generations, the population is evenly distributed across a subset V ⊂ F ∗ already covered by Pt. Furthermore, if the population Pt is Pareto optimal in generation t and all future ones, the distribution can be even more characterized by bounding the cover number of every v ∈ F ∗ by ⌈µ/|V |⌉ from above. Notably, for |V | = |F ∗ |, this is the best possible distribution of Pt across the whole Pareto front one can achieve. We also demonstrate that once such a distribution is reached, it is maintained for all future generations.

1 A long version of this paper can be found at https://arxiv.org/abs/2505.01256 where the reader can find proofs for all lemmas stated in the paper.

Lemma 4. Assume the same conditions as in Lemma 1 where stochastic population update is disabled. Suppose µ ∈ 2O(n) . Let V ⊂ F ∗ be covered by Pt. The following properties are satisfied after O(n) generations and all future ones with probability at least 1 - e-Ω(n) and in expectation.

1. If all populations P0, P1, . . . consist only of Pareto optimal individuals, we have that ct(v) ≤ ⌈µ/|V |⌉.

Proof. (1): Fix a Pareto optimal x with v := f(x) ∈ V . We show that ct(v) ≥ β with probability at least 1 - e-Ω(n) after O(n) generations. By Lemma 2(1), ct(v) ≤ β cannot decrease since a maximum set of mutually incomparable solutions I fulfills |I| ≤ (2n/d + 1)d/2 and hence, ct(v) ≤ β ≤ µ/|I|. We define two phases where the second phase only applies if β > n.

Phase 1. We have ct(v) ≥ ℓ := min{β, n}.

In the next lemma we determine the expected duration of each phase under the assumptions of Lemma 4.

Lemma 5. Both phases are finished in O(n) generations with probability at least 1 - e-Ω(n) and in expectation.

Now we obtain by a union bound on both phases and every v ∈ V that in O(n) generations the cover number of every v ∈ V is at least β with probability at least 1-e-Ω(n) (since |V | = O(nd/2 ) = eo(n) ). The bound on the expected number of generations follows by applying the same arguments for an additional period of O(n) generations and by the fact that (1 + o(1))O(n) = O(n) such periods are sufficient.

(2): With the same argument as in (1) we obtain with probability 1-e-Ω(n) that, after O(n) generations, the cover number of all v ∈ V is at least ⌈µ/|V |⌉ or one of these cover numbers has decreased at least one time when it was at most ⌈µ/|V |⌉. Suppose that the former happens. If ⌈µ/|V |⌉ > µ/|V | we see that ⌈µ/|V |⌉ · |V | > |V | · µ/|V | = µ, a contradiction. Hence, ⌈µ/|V |⌉ = µ/|V | and all µ individuals are completely evenly distributed on V . If the latter happens, we see with Lemma 2(2) that the cover number of all vectors v ∈ F ∗ is at most ⌈µ/|V |⌉. By Lemma 2(3) the maximum cover number cannot increase, proving Lemma 4.

6 Upper Runtime Bounds

In this section we prove upper runtime bounds for NSGA-III with stochastic population update (a = 1), as well as without (a = 0), on f := d-OJZJ. For the variant where a = 1 we follow the proof in [Bian et al., 2023]. For technical reasons, we define for a Pareto optimal fitness vector v = f(x) the string L(v) ∈ {0, 1, ⊥}d as follows. For j ∈ [d] let

• L(v)j = 1 if vj = 2n/d + k and j is odd (attained if x(j+1)/2 = 1n since then fj(x) = 2n/d + k),
• L(v)j = 0 if vj = 2n/d + k and j is even (attained if xj/2 = 0n since then fj(x) = 2n/d + k), and
• L(v)j = ⊥ if k ≤ vj ≤ 2n/d (attained if k ≤ |x(j+1)/2 |1 ≤ 2n/d - k or x = 02n/d if j is odd and k ≤ |xj/2 |1 ≤ 2n/d - k or x = 12n/d if j is even).

Note that for every Pareto optimal x there is w ∈ {0, 1,⊥}d with w = L(f (x)).

Theorem 6. Assume the same conditions as in Lemma 1 for f := d-OJZJk where d is constant, 2 ≤ k ≤ n/(2d) and µ ∈ 2O(n) . Then the expected number of generations until the whole Pareto front is covered is at most O(nk+d/2 /µ + n ln(n)) for NSGA-III if a = 0, and O(k( 12en k )k ) for NSGA-III if a = 1.

Proof. We fix v ∈ F ∗ with string L(v) ∈ {0, 1,⊥}d and show that all Pareto optimal w ∈ F ∗ with L(w) = L(v) are covered after O(nk+d/2 /µ + n ln(n)) generations for NSGA-III if a = 0 and after O(k(12en/k)k ) generations for NSGA-III if a = 1. Since there are at most 3d different strings L(v) ∈ {0, 1,⊥}d , the runtime follows since d = Θ(1). Let Ov := {i ∈ [d] | (L(v))i = 1} and Zv := {i ∈ [d] | (L(v))i = 0}. By a classical Chernoff bound, the probability is at most e-Ω(µn) that for all initialized x there is a block j ∈ [d/2] with |xj |1 /∈ {k, . . . , 2n/d - k}. Suppose that this happens. Then the probability is at least n-n to create any individual with mutation. Hence, the expected number of generations to create a Pareto optimal x with Of(x) = Zf(x) = ∅ (that is, L(f(x)) = {⊥}d ) is 1 + n-n e-Ω(µn) = 1 + o(1). Now, suppose that there is a Pareto optimal x (that is, L(f(x)) ∈ {0, 1,⊥}d ) with Of(x) ⊊ Ov or Zf(x) ⊊ Zv. Let i ∈ Ov \ Of(x). Note that i is odd. Now we determine the time until all Pareto optimal fitness vectors w with L(w) = L(f(x)) are covered, captured by the following lemma where the same assumptions as in Theorem 6 hold.

Lemma 7. Let z ∈ Pt be Pareto optimal. Then in 8dn ln(n) generations (i.e. 8µdn ln(n) fitness evaluations) all Pareto optimal v with L(v) = L(f(z)) are covered with probability at least 1 - n-d . The expected number of generations is O(n ln(n)).

To generate a Pareto optimal y with Zf(y) = Zf(x) and Of(y) = Of(x) ∪ {i} we consider both versions of NSGA-III separately.

For NSGA-III without stochastic population update an additional phase of expected O(n) generations ensures that every fitness vector w ∈ F ∗ with L(w) = L(f(x)) has cover number at least c := ⌊µ/(2n/d + 1)d/2 ⌋ (by Lemma 4(1)). Then, to generate y, one can choose z with L(f(z)) = L(f(x)) and |z(i+1)/2 |1 = 2n/d-k as a parent (prob. c/µ), and flip k specific bits while keeping the remaining bits unchanged (prob.

(1-1/n)n-k /nk ≥ 1/(enk )). Hence, in one generation this happens with probability at least

where the first inequality is due to Lemma 10 in [Badkobeh et al., 2015]. Hence, in total, we obtain O(nk /c+n ln(n)) = O(nk+d/2 /µ + n ln(n)) generations in expectation to generate y.

For NSGA-III with stochastic population update we consider a sequence of k successive generations and for ℓ ∈ [k] we call the ℓ-th generation successful if an individual z is created with 2n/d - k + ℓ ones in block (i + 1)/2 while zj = xj for j ∈ [d/2] \ {(i + 1)/2} and z is not removed. Suppose that the (ℓ - 1)-th generation is successful. Then, generation ℓ is successful if in one trial one chooses a parent p with |p(i+1)/2 |1 = 2n/d - k + ℓ - 1 and zj = xj for j ∈ [d/2] \ {(i + 1)/2}, flips a zero bit in the ((i + 1)/2)-th block while keeping the remaining bits unchanged and finally keeps the new created individual. With probability at least

a desired individual z is created in generation ℓ. Note that z is not Pareto-optimal for ℓ < k. However, µ - ⌈3µ/2⌉ individuals chosen uniformly at random survive (see Line 10 and Line 18 in Algorithm 1 for a = 1) and hence even a non Pareto optimal individual remains with probability at least (µ - ⌈3µ/2⌉)/µ = ⌊µ/2⌋/µ ≥ (µ/2 - 1)/µ = 1/2 - 1/µ. For n sufficiently large, we have 1/2 - 1/µ ≥ 1/3 and hence with probability at least (k - ℓ + 1)/(12n), generation ℓ is successful. If all k generations are successful, the desired y is created, which happens with probability at least

(since k! ≥ (k/e)k by Stirling's approximation). Hence, the expected number of generations to create the desired y is at most k(12en)k /kk .

In either case, by symmetry, the above bounds also apply to create a Pareto optimal y such that Zf(y) = Zf(x) ∪ {i} for i ∈ Zv \ Zf(x) and Of(y) = Of(x) else. Since Ov and Zv are both finite, we obtain that a Pareto optimal individual z with Of(z) = Ov and Zf(z) = Zv is created in expected O(nk+d/2 /µ + n ln(n)) generations when stochastic population update is turned off and O(k( 12en k )k ) it it is turned on. With Lemma 7 we obtain that all Pareto optimal v with L(v) = L(f(z)) are covered in expected O(n ln(n)) generations, proving the theorem.

We see that for µ ∈ Ω(nd/2 ) ∩ O(nk+d/2-1 / ln(n)), the derived runtime if a = 0 is O(nk+d/2 ) in terms of fitness evaluations and, hence, independent of the population size µ. An important consequence is that, if the number of objectives is 2, NSGA-III outperforms NSGA-II if also µ ∈ ω(nd/2 )∩o(n2 /k2 ) (compare with [Doerr and Qu, 2023a] for the tight runtime bound of Θ(µnk ) of NSGA-II on 2-OJZJ in terms of fitness evaluations if n + 1 ≤ µ ∈ o(n2 /k2 )).

7 Lower Runtime Bounds

In this section we prove sharp lower runtime bounds of NSGA-III without stochastic population update on OJZJ for the bi-objective case which are tight for many parameter settings. To this purpose, we show that with probability 1-e-Ω(n) after O(n) generations, the cover number of every v ∈ F ∗ can be bounded from above by a sufficiently small value, which also remains for all future iterations. In other words, the population remains well-distributed across F ∗ .

Lemma 8. Assume the same conditions as in Lemma 1 for f := d-OJZJ, a = 0, d is constant, 2 ≤ k ≤ n/(2d) and µ ∈ 2O(n) . Suppose that all P0, P1, . . . ever seen by NSGA-III consist only of Pareto optimal individuals. Then, with probability at least 1-e-Ω(n) the cover number of each Pareto optimal fitness vector is at most ⌈2d/2+1µ/|I|⌉ after O(n) generations.

However, in the case d > 2 it may happen that non-Pareto optimal search points with rank one are created which survive (for example if, outgoing from a population on the Pareto front without 02n/d or 12n/d in each block j ∈ [d/2], µ - 1 individuals are cloned and one non-Pareto optimal individual y is generated with yi = 12n/d in one block i ∈ [d/2] while 1 ≤ |yj |1 < k in another block j). Such non-Pareto optimal solutions can remain for a long time in the population and hence, Lemma 2(3) cannot be applied from scratch, complicating the analysis. Hence, we present only the case d = 2.

Theorem 9. Assume the same conditions as in Lemma 1 for f := d-OJZJ and a = 0 where d = 2, 2 ≤ k ≤ n/4 and µ ∈ O(nk-1 ) ∩ poly(n). Then the expected number of generations required for covering the entire Pareto front is at least Ω(nk+1 /µ).

Proof. By a classical Chernoff bound, with probability 1 - e-Ω(n) , every individual x fulfills |x|1 ∈ {k, . . . , n - k} after initialization. Suppose that this happens. Then there are only Pareto optimal individuals in P0 and all future populations P1, P2, . . . since a non-Pareto optimal individual is dominated by every Pareto optimal one (see Lemma 3). Since there are µ Pareto optimal individuals in Rt, such a y never survives (due to the non-dominated sorting procedure). Note that |I| = |F ∗ | (since |F ∗ | = n - 2k + 3 ≥ n/2 + 3, f attains n+1 different values in each objective and every nonPareto optimal individual is dominated by a Pareto optimal one). By Lemma 8 applied on d = 2, with probability at least 1 - e-Ω(n) the cover number of every v ∈ F ∗ is at most ⌈4µ/|I|⌉ after at most O(n) generations. Suppose that this happens. Let n be sufficiently large such that µ ≤ cnk-1 for a constant c > 0 and the event from Lemma 8 applies on at most δn generations for a further constant δ > 0. Further, suppose that 1/n+ 5n/|F ∗ | ≤ α for a suitable constant α > 0 (due to |F ∗ | ∈ Θ(n)). We consider a phase of at most δn generations (i.e. µδn trials) and show that no individual y ∈ {0n , 1n } is generated within µδn trials with at least constant probability. Since in one trial one has to flip at least k specific one bits or k specific zero bits to create y, the probability that this happens in µδn trials is at most 1-(1-2/nk )µδn ≤ 1-(1-2/nk )cδnk ≤ 1-(1/16)cδ =: b where we used (1 - 2/n)n ≥ ((1 - 1/n)n )2 ≥ 1/16. Note that b < 1 is a constant. Hence, we obtain with probability at least 1-b-e-Ω(n) = Ω(1) that there is a generation t where every v ∈ F ∗ has cover number at most ⌈4µ/|F ∗ |⌉, and every individual x ∈ Pt fulfills x /∈ {1n , 0n }. Suppose that this happens. We estimate the expected time to create x = 1n (which is the only Pareto optimal search point with fitness f1(x) = 2n/d + k) from above. Note that for a given k ≤ ℓ ≤ 2k - 1 there are at most ⌈4µ/|F ∗ |⌉ ≤ 5µ/|F ∗ | different individuals y ∈ Pt with Hamming distance H(1n , y) = ℓ to 1n . Note also it requires to flip at least 2k specific bits to create 1n from a y with H(1n , y) ≥ 2k. Since k ≥ 2, all these considerations yield a probability of at most

to create 1n in one trial. Hence, in one generation we obtain by a union bound on µ trials that 1n is created with probability at most αµ(1/n)k+1 . Hence, the expected number of generations to create 1n is at least (1 - b - e-Ω(n) )nk+1 /(αµ) = Ω(nk+1 /µ) which proves the theorem.

Combining Theorems 6 and 9, we obtain for |I| = |F ∗ | ≤ µ ∈ O(nk-1 ) ∩ poly(n), d = 2 and 2 ≤ k ≤ n/4 the tight runtime of Θ(nk+1 /µ) for NSGA-III on 2-OJZJ.

8 Conclusions

In this paper we provided new insights in the population dynamics of NSGA-III, demonstrating that it quickly spreads its solutions evenly across the whole Pareto front. Our methods were developed with a level of generality that makes them applicable to a broad range of problems. Subsequently, we derived upper runtime bounds for NSGA-III on the ONEJUMPZEROJUMP benchmark, and, for certain regimes of gap size k and population size µ, even tight runtime bounds when the number of objectives is two. Additionally, we showed that stochastic population update, where solutions for survival are not always selected deterministically, can lead to an exponential speedup in runtime. We hope that the techniques developed mark a significant step forward in understanding the behavior of NSGA-III, while also shedding light on its strengths and limitations. To deepen this understanding, future work could focus on the usefulness of stochastic population update on more complex problems with many local optima like the many-objective REALROYALROAD function, proposed by [Opris, 2025], where the algorithm has to cross a large fitness valley. Another research direction could be deriving rigorous lower bounds for NSGA-III on d-OJZJ for d > 2, or on further classical problems such as OMM, COCZ, and LOTZ. We also hope that the gained deeper theoretical understanding of the dynamics of NSGA-III on OJZJ will offer useful insights for practitioners, such as facilitating the creation of refined versions of the algorithm, particularly to enhance performance in real world scenarios where many local optima occur.

Acknowledgements

This work benefited from fruitful discussions at Dagstuhl Seminar 24271 "Theory of Randomized Optimization Heuristics".