On the Impact of Crossover in Many-Objective Optimization: A Runtime Analysis of NSGA-III

1 Introduction

Multi-objective evolutionary algorithms (MOEAs), such as the non-dominated sorting genetic algorithm II (NSGAII) [Deb et al., 2002] and its extension for many-objective optimization, NSGA-III [Deb and Jain, 2014], have been used in thousands of applications and together have received more than 60,000 citations. Applications span a wide range of domains, including machine learning [Bi-wu Zhu et al., 2025], bioengineering [Rashmi et al., 2025], and artificial intelligence [Luukkonen et al., 2023], where also many studies involve four or more objectives. The most prominent EMOA for optimizing bi-objective problems is NSGA-II [Deb et al., 2002] (see [Vijai and P., 2025] for empirical results, or [Zheng et al., 2022; Dang et al., 2023] for rigorous ones), however, which is not well suited for solving problems where the number of objective increases [Campos Ciro et al., 2016; Doerr et al., 2025]. The reason is, that in the bi-objective case, an ordering with respect to the first objective of nondominated solutions implies also an ordering with respect to the second, which makes the crowding distance, the second tie breaker of NSGA-II, effective. However, this relation breaks down already for three objectives, which indicates the inefficiency of NSGA-II in the many-objective setting. However, NSGA-III, a refinement of NSGA-II, uses reference points instead of crowding distance to ensure that the solution set will be well-distributed across the objective space in a very natural way (see [Deb and Jain, 2014]). It has been shown both theoretically [Wietheger and Doerr, 2023; Opris et al., 2024a] and empirically [Campos Ciro et al., 2016] that NSGA-III effectively optimizes many-objective problems, which is a challenging task, as the size of the Pareto front and the number of incomparable solutions grow exponentially with the number of objectives. However, these first rigorous runtime analyses of the state-of-the-art NSGA-III were published just a few years ago. As a result, its theoretical understanding still lags behind its practical achievements, although some progress has been made since then [Opris, 2025a; Opris, 2026a; Opris, 2026b]. There is still a significant gap in our theoretical understanding of how crossover contributes to many-objective optimization, despite its importance as a fundamental operator in evolutionary computation [Pavai and Geetha, 2016] which hinder our general understanding when and why MOEAs perform well. Beyond up to two papers [Opris, 2025c; Opris, 2026a], which theoretically demonstrate that crossover yields an exponential runtime speedup on handcrafted benchmark functions, we are not aware of any further theoretical results, in contrast to the bi-objective setting [Dang et al., 2023; Doerr and Qu, 2023a]. But further empirical results on many different multiobjective constrained, and unconstrained problems [Sharma et al., 2021], and on multi- and many-objective Knapsack problems [Ishibuchi et al., 2014] show the huge potential of crossover operators in many-objective optimization.

Our contribution: In this paper, we build on the considerations from [Opris, 2025a; Doerr and Qu, 2023a], and provide a theoretical runtime analysis of NSGA-III on the many-objective m-OJZJk benchmark with and without crossover for any number of objectives m, and show that NSGA-III without crossover needs O(m2n ln(n/m) + δmnk (1 + 2n/m)m/2 /µ) generations with probability 1 - o(1) - e-(δ-3)m to cover the whole Pareto front. Here, µ denotes the population size, n the problem size, k is a parameter specific to the problem at hand, and δ > 3 is an additional parameter. The expected number of generations is O(m2n ln(n/m) +mnk (1 + 2n/m)m/2 /µ). One sees that this number is asymptotically the same for a wide range of population sizes µ. This robustness stems from the observation that NSGA-III can retain many individuals with the same fitness vector due to how it associates solutions with and iterates over reference points. This aspect was not considered in [Wietheger and Doerr, 2024]. In particular, our analysis improves the runtime bound given there for population sizes asymptotically larger than the size of the Pareto front, and it even extends the analysis from [Opris, 2025a] to an arbitrary number of objectives. With uniform crossover, we show that with probability 1 - o(1) - e-(δ-3)m, the number of generations until the whole Pareto front is covered is at most O(kδmk+1nk ck µ/k!) for NSGA-III, and the expected number of generations is O(kmk+1nk ck µ/k!) for a suitable constant c > 0. Thus, one obtains a speedup of order Ω((k - 1)!(1 + 2n/m)m/2 /(µ2 ckmk )) in the expected runtime compared to the case without crossover. If k = Ω(n/ ln(n)), m = O(log(n)) and µ = O((1 + 2n/m)m/2 ), this speedup even becomes exponential. In a nutshell, our proof applies the arguments from [Doerr and Qu, 2023b] for the bi-objective case sequentially to all blocks. A key difference is that we also need high-probability guarantees for finding a single Pareto-optimal solution in order to obtain all of them in reasonable time in parallel. Finally, we complement our analysis by also providing a lower runtime bound of NSGA-III without crossover on 4-OJZJ of Ω(nk /µ) generations, which is by a factor of (k - 1)!/(µ2dk ) larger than the upper runtime bound derived with crossover for a constant d > 0. This factor is also exponential if k = Ω(n/ ln(n)) and µ = poly(n). Extending this lower bound to a larger number of objectives appears considerably more difficult, since the interactions between single objectives becomes more complex, and requires a much deeper understanding of the underlying population dynamics which extends the scope of this paper. However, we expect our results to extend similarly to other MOEAs such as GSEMO, SPEA2, SMS-EMOA, and variants of PAES-25 [Opris, 2025b].

Related work: In single-objective optimization, the benefits of crossover are much better understood than in the manyobjective setting. On pseudo-Boolean benchmark problems such as JUMPk, where a fitness valley of size k must be crossed, it has been rigorously shown that uniform crossover yields a speedup depending on k and the crossover probability [Jansen and Wegener, 2002; K¨otzing et al., 2011; Dang et al., 2017; Opris et al., 2024b; Opris and Antipov, 2026]. An exponential performance gap in the runtime was proven in [Jansen and Wegener, 2005] on a function REALROYALROAD, which is specifically designed for 1point crossover. These insights have been used to prove also advantages through crossover for combinatorial optimization problems like shortest paths [Doerr et al., 2012; Doerr et al., 2013], solving complex data clustering problems [Gaeuman and Sutton, 2025] or NP-hard graph problems [Sutton and Lee, 2024], or even in more complex search spaces like permutation spaces [Opris et al., 2025].

In multiobjective optimization, only a few variants of the global simple evolutionary multiobjective optimizer (GSEMO) with crossover have been studied [Qian et al., 2013; Qian et al., 2020; Doerr et al., 2022], and rigorous analyses of NSGA-II with crossover on classical benchmark problems [Doerr and Qu, 2023b] and multi-objective variants of REALROYALROAD [Dang et al., 2023] have been conducted. However, these results were only restricted on bi-objective problems. The theoretical analysis of NSGA-III only succeeded recently. Based on a rigorous analysis of GSEMO on classical benchmark functions [Laumanns et al., 2004], in [Wietheger and Doerr, 2023] the first runtime analysis of NSGA-III on the easy 3-ONEMINMAX problem was conducted. This was then generalized by [Opris et al., 2024a; Opris, 2025a] on more than three objectives, where also key structural insights into the working principles of NSGA-III are given. Also, other pseudo-Boolean functions, where one has also to reach the Pareto front at a first glance, have been analyzed there. Then, in [Wietheger and Doerr, 2024; Opris, 2025a] m-OJZJ has been analyzed, but without investigating crossover. Similar analyses have then also been conducted on other popular MOEAs like the SPEA-2 [Ren et al., 2024] and the SMS-EMOA [Zheng and Doerr, 2024]. First theoretical results which showcase that MOEAs, particularly NSGA-III and SPEA-2, are quite robust with respect to the chosen population size can be found in [Opris, 2026b; Doerr et al., 2026]. However, apart from [Opris, 2025c; Opris, 2026a], we are not aware of any theoretical results on whether or how crossover can be beneficial in many-objective optimization. Moreover, these results are limited to Royal Road functions designed for specific crossover operators.

2 Preliminaries

Notation: For a finite set A, we write |A| for its cardinality. For n ∈ N, define [n] := {1, . . . , n} and 1n the vector of length n with only ones, while 0n is the corresponding vector with only zeros. Given a bit string x ∈ {0, 1}n , let |x|1 and |x|0 denote the number of ones and zeros in x, respectively. For x, y ∈ {0, 1}n denote by H(x, y) = ∑ni=1 |xi - yi| the Hamming distance of x and y. We use ln to denote the natural logarithm. Let Y and Z be random variables taking values in N0. We say that Z stochastically dominates Y if Pr(Z ≤ c) ≤ Pr(Y ≤ c) for all c ≥ 0. Consider an m-objective function f : {0, 1}n → Nm0 , x ↦→ (f1(x), . . . , fm(x)), and let fmax := max{fj(x) | j ∈ [m], x ∈ {0, 1}m} the maximum possible value of an objective. For two search points x, y ∈ {0, 1}n , we say that x weakly dominates y, denoted x ⪰ y, if fi(x) ≥ fi(y) for all i ∈ [m]. If, in addition, at least one of these inequalities is strict, then x dominates y, written x ≻ y. If neither x ⪰ y nor y ⪰ x holds, the two points are called incomparable. We say that a subset S ⊆ {0, 1}n consists only of mutually incomparable solutions if every pair of elements in S is incomparable. A solution x is called Pareto-optimal if it is not dominated by any other search point in {0, 1}n , and we call the set {f(x) | x ∈ {0, 1}n is Pareto optimal} the Pareto front. For a population Pt and a fitness vector v ∈ Nm0 , the cover number ct(v) is defined as the number of individuals x ∈ Pt with f(x) = v. We say that v is covered if ct(v) ≥ 1.

The NSGA-III algorithm: The NSGA-III algorithm ( [Deb and Jain, 2014]) with crossover probability pc ∈ [0, 1) and Algorithm 1: NSGA-III on an m-objective function f with population size µ and crossover probability pc

1 Initialise P0 ∼ Unif(({0, 1}n )µ ) 2 for t := 0 to∞do 3 Initialise Qt := ∅ 4 for i = 1 to µ/2 do 5 Select a1, a2 from Pt uniformly at random 6 Select u ∼ Unif([0, 1]) 7 if u ≤ pc then 8 Create y1 by uniform crossover on a1, a2 9 Create y2 by uniform crossover on a1, a2 10 else 11 Create y1, y2 as copies from a1, a2 12 Create z1, z2 by standard bit mutation on y1, y2 13 Update Qt := Qt ∪ {z1, z2} 14 Set Rt := Pt ∪Qt 15 Partition Rt into layers F 1 t , F 2 t , . . . , F k t of non-dominated fitness vectors 16 Find i∗ ≥ 1 such that ∑i∗-1 i=1 |F i t | < µ and ∑i∗ i=1 |F i t | ≥ µ 17 Compute Yt = ⋃i∗-1 i=1 F i t 18 Choose ˜ F i∗ t ⊂ F i∗ t such that |Yt ∪ ˜ F i∗ t | = µ with Algorithm 2 19 Create the next population Pt+1 := Yt ∪ ˜ F t i∗

even population size µ is presented in Algorithm 1. Initially, a population P0 of size µ is generated by selecting µ individuals uniformly at random from {0, 1}n . In each generation t, an offspring population Qt of size µ is created by performing the following operations µ/2 times. First, two parents a1 and a2 are selected uniformly at random from Pt. Then, uniform crossover is applied to (a1, a2) two times with probability pc to produce two intermediate solutions y1 and y2. That is, for creating one solution, and for each position i ∈ [n] independently, the entry from a1 is taken with probability 1/2, and otherwise the entry from a2. If uniform crossover is not applied, y1 and y2 are exact copies of a1 and a2. Finally, two offspring z1 and z2 are generated by applying standard bit mutation to y1 and y2, that is, 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 and Rt is updated by partitioning Rt into layers F 1 t , F 2 t , . . . using the non-dominated sorting algorithm [Deb et al., 2002] where F 1 t consists of all non-dominated individuals, and F i t for i > 1 of individuals only dominated by those from F 1 t , . . . , F i-1 t . Then the critical rank i∗ with ∑i∗-1 i=1 |F i t | < µ and ∑i∗ i=1 |F i t | ≥ µ 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∗ ). All individuals with a lower rank than i∗ are included in Pt+1, while the remaining individuals are selected from F i∗ t using Algorithm2. Hereby, a normalized objective function f n is computed and then each individual with rank at most i∗ is associated with reference Algorithm 2: Selection procedure utilizing a set Rp of reference points for maximizing a function, including uniform selection.

• 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 minimised
• 3 For each r ∈ Rp, set ρr := |{x ∈ Yt | rp(x) = r}|
• 4 Initialise ˜ F i∗ t = ∅ and R′ :=Rp
• 5 while | ˜ F i∗ t | < µ/2 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 minimises the distance between the vectors f n (xrmin ) and rmin (where ties are broken randomly)

15 Select µ- |Yt| -µ/2 distinct individuals uniformly at random from F i∗ t \ ˜ F i∗ t , and add them to ˜ F i∗ t

• 16 return ˜ F i∗ t

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., 2024a]. We omit detailed explanations as they are not needed for our purposes. For an m-objective function f : {0, 1}n → Nm0 , the normalized fitness vector f n (x) := (f n 1 (x), . . . , f nm(x)) of a search point x is computed as

for each j ∈ [m] where ynad := (ynad 1 , . . . , ynad m ) and ymin := (ymin 1 , . . . , ymin m ) 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 ∈ [m] 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. 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, if the number of all individuals already chosen for the critical layer is at most µ/2, 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, or |F i∗ t | = µ/2, the selection ends. If after this selection procedure still |Yt| + | ˜ F i∗ t | < µ, then the remaining µ - |Yt| - | ˜ F i∗ t | = µ - |Yt| - µ/2 = µ/2 - |Yt| individuals from Fi∗ \ ˜ F i∗ t are chosen uniformly at random. Such a uniform selection strategy has been shown to be successful, since it enables movement on fitness neutral environments (called plateaus), which helps to build up and preserve population diversity [Opris et al., 2024b; Dang et al., 2017; Doerr and Qu, 2023a]. The following result can be formulated and proven as in [Opris, 2025a]. For completeness, we provide a proof in the appendix.

Lemma 1. Consider NSGA-III on an m-objective function f with Pareto front F, and assume that εnad ≥ fmax. Let Rp denote a set of reference points, with p ≥ 2m3/2 fmax. Let Pt be the population at iteration t. Let S be a maximum set of mutually incomparable solutions, let v ∈ F, and 0 ≤ α ≤ ⌊µ/(2|S|)⌋. Then if ct(v) ≥ α then also ct+1(v) ≥ α.

The many-objective m-OJZJk benchmark: This benchmark has been defined the first time in [Zheng and Doerr, 2024] and is defined as follows, where m-OJZJk := (f1, . . . , fm), m is even, and n is divisible by m/2 (see also [Doerr and Qu, 2022] for the bi-objective version). For 2 ≤ k ≤ 2n/m the m-OJZJk = (f1(x), . . . , fm(x)) is defined as

if j ∈ [1, . . . ,m] is odd, and

if j ∈ [1, . . . ,m] is even. We often call k gap size. For every objective there are 2n/m + 1 different values and fmax = k + 2n/m. The Pareto front F of m-OJZJk is {(ℓ1, 2k + 2n/m - ℓ1, . . . , ℓm/2, 2k + 2n/m - ℓm/2) | ℓ1, . . . , ℓm/2 ∈ {k, 2k, 2k+1, . . . , 2n/m-1, 2n/m, 2n/m+ k}}, has cardinality (2n/m-2k+3)m/2 for k ≤ n/m, and a maximum set S of mutually incomparable solutions satisfies |F| ≤ |S| ≤ (2n/m+ 1)m/2 (see [Zheng and Doerr, 2024] for proofs).

For v ∈ F we introduce the following notation: Denote by rv ∈ {-1, 0, 1}m/2 the (m/2)-dimensional vector with (rv)j = -1 if v2j-1 = k, (rv)j = 1 if v2j-1 = 2n/m+ k, and (rv)j = 0 if 2k ≤ v2j-1 ≤ 2n/m. All search points x with f(x) ∈ {v ∈ Nm0 | rv = s} for s ∈ {-1, 0, 1}m/2 satisfy xj = 0m/2 if (rv)j = -1, xj = 1m/2 if (rv)j = 1, and k ≤ |xj |1 ≤ 2n/m- k if (rv)j = 0 for all j ∈ [m/2].

3 An Upper BoundWithout Crossover

To compare the performance of NSGA-III with and without crossover, we first generalize a result from [Opris, 2025a] for NSGA-III on m-OJZJk to an arbitrary number of objectives. This generalization requires a much more refined proof, since for large m the number of Pareto-optimal points can be exponentially. Interestingly, the runtime bound from [Opris, 2025a] carries over directly to this setting in terms of generations. In addition, we provide an upper bound on the runtime that holds with high probability. Our results are formulated for arbitrary crossover probabilities pc ∈ [0, 1), as parts of the analysis are later reused for the crossover case. However, the analysis here relies only on mutation steps.

Theorem 2. Consider NSGA-III on f := m-OJZJk for 2 ≤ k ≤ n/(2m) + 1, (1 + 2n/m)m/2 ≤ µ/2 ≤ no(n) , crossover probability pc ∈ [0, 1), and a number m of objectives with 2 ≤ m ≤ n/2. Further, assume the same conditions as in Lemma 1, and let δ > 0. Then with probability at least 1 - e-(δ-3)m -O(1/n4 ) the number of generations until the whole Pareto front F is covered is at most

The expected number of generations is

Proof. At first we prove that, with probability at most e-Ω(n2 ) , there exists no Pareto optimal individual x with k ≤ |xj|1 ≤ 2n/m-k for all j ∈ [m/2] after initialization. Then, with probability 1 - e-Ω(n2 ) , a single generation suffices to create such a Pareto optimal individual since any individual x can be created with probability at least n-n , independently of whether crossover is executed or not. This implies that the expected number of generations for creating a Pareto optimal individual is at most 1 + nn · e-Ω(n2 ) = 1 + o(1).

Lemma 3. With probability at most e-Ω(n2 ) , there exists no Pareto optimal individual x with k ≤ |xj |1 ≤ 2n/m - k for all j ∈ [m/2] after initialization.

After a successful initialization, we explore search points x satisfying rf(x) = s for all possible s ∈ {-1, 0, 1}m/2 in parallel. Specifically, we determine a suitable number of generations such that, for any fixed s, a search point x with rf(x) = s is found with high probability. We then apply a union bound over all 3m/2 possible s, which shows that, still with high probability, a corresponding x with rf(x) = s is found for every s within that time.

We fix an s ∈ {-1, 0, 1}m/2 and define Step j as the process that starts when there is a point y ∈ Pt with (rf(y))i = si for all i < j and (rf(y))i = 0 for i ≥ j, and ends when a search point z is generated with (rf(z))i = (rf(y))i for all i ≤ j, and (rf(z))i = 0 for i > j. Note that Step 1 directly builds on a successful initialization (since then there exists a search point z ∈ Pt with rf(z) = 0m/2 ), and that the desired z is created after Step m/2. For Step j we fix such a y ∈ Pt. Then in all future iterations there is always a y′ ∈ Pt with f(y′ ) = f(y) by Lemma 1, and hence, we can never fall back in steps. Further, we define Bj := {v ∈ F | vi = fi(y) for every i ∈ [m] \ {2j -1, 2j} and v2j-1 ∈ {2k, . . . , n-2k}}. Note that f(y) ∈ Bj. Now we estimate the number of generations to finish Step j.

• At first we estimate the number Xj of generations until the complete set Bj is covered. Particularly, that for every v ∈ Bj , there is x∗ ∈ Pt with f(x∗ ) = v.
• After Xj generations, we estimate the number of generations Yj until every Pareto optimal fitness vector v ∈ Bj has a cover number of at least ⌊µ/(2(1 + 2n/m)m/2 )⌋.
• Finally, after Xj +Yj generations, estimate the number of generations Zj until a desired z is created by possibly crossing the fitness valley of size k if necessary.

Then the total number of generations until a desired x is created is stochastically dominated by ∑m/2 j=1 (Xj + Yj + Zj). We now derive tail bounds for the random variables X := ∑m/2 j=1 Xj, Y := ∑m/2 j=1 Yj, and Z := ∑m/2 j=1 Zj separately in three consecutive lemmas. We only provide some proof sketches due to space restrictions and similar ideas to the considerations from [Opris, 2025a]. Their full proofs can be found in the appendix.

Lemma 4. Fix β > 0. Then for j ∈ [m/2] and n sufficiently large,

Particulary,

We fix an uncovered Pareto-front fitness vector v and bound the number of generations until it is covered with high probability. A union bound over all v ∈ Bj then yields a highprobability bound for covering all such vectors, which extends to X via another union bound by multiplying with m/2. We next estimate Y .

Lemma 5. For j ∈ [m/2] we have

Particulary,

The idea behind proving that each fitness vector from Bj attains the desired cover number after Ω(n/(1 - pc)) generations is again based on a parallelization argument. We split the proof into two phases. In the first phase, we determine an expected number of generations such that the cover number of a fixed v ∈ Bj is at least Ω(n) with high probability, which follows from repeatedly cloning an individual x with f(x) = v. In the second phase, we show that the desired cover number ct(v) is reached by applying a classical Chernoff bound to the number of newly created individuals in a single generation, which is Ω(n) in expectation. By a union bound, we obtain the claim for all v ∈ Bj, and finally also for Y . Next, we extimate Z.

Lemma 6. For γ := ⌊µ/(2(1 + 2n/m)m/2 )⌋ and j ∈ [m/2] the random variable Zj is stochastically dominated by a geometrically distributed random variable Z ′j with success probability σk := (1-pc)γ/(2enk )

We estimate Zj pessimistically as the number of generations required to select an individual with 2n/m-k ones (or, symmetrically, 2n/m- k zeros) in block j, of which there are at least ⌊µ/(2(1+2n/m)m/2 )⌋ many, and then flip k specific bits in block j to cross the fitness valley. Again, we apply a union bound to estimate Z.

Finally, we are in a position to apply a union bound to the random variables X, Y and Z. Using Lemmas 4, 5, and 6, respectively, we obtain for every δ ≥ 1 and

the inequality

So after K(k,m, n, pc, δ) generations, a desired x with rf(x) = s is created with probability at most p(m, n, δ), after a successful initialization. A union bound on at most 3m/2 such possible x with different rf(x)-values shows that with probability at most (due to 2n/m ≥ 4)

that for each s ∈ {-1, 0, 1}m/2 there is x with rf(x) = s. Suppose that this happens. Then we cover the Pareto front. Fix an uncovered v ∈ F. In an analogous way to the argument above, the number of generations required to cover v is stochastically dominated by X. This is because there already exists a search point y ∈ Pt with rf(y) = rv and hence, for every block j ∈ [m/2] with (rv)j = 0, the time to create a search point z satisfying (f2j-1(z), f2j(z)) = (v2j-1, v2j), while not changing all remaining blocks, is stochastically dominated by Xj . By Lemma 4, we obtain

By a union bound over all such possible search points, we obtain that the whole Pareto front is covered in an additional amount of (4+16m) ·mn/2 · ln(2n/m)/(1-pc) generations with probability at most (1 + 2n/m)m/2 · (2n/m)-4m ≤ (2n/m)-2m. Hence, for p := p(m, n, δ), the entire Pareto front is covered within

generations with probability at least

where we also used that the functions h1 : ]0, n/2[ → R, x ↦→ (2n/x)-2x , and h2 : ]0, n/2[ → R, x ↦→ x/2 · (2n/x)-2x , are strictly monotone decreasing (see Lemma ?? in the appendix). The latter propability also includes the event that in the first generation a successful initialization happens.

Now it remains to estimate the expected time until the whole Pareto front is covered. Here we assume δ = 4. If, after a successful initialization, the Pareto front is not covered after O(m2n ln(n/m)/(1 - pc) + mnk (1 + 2n/m)m/2 /((1 - pc)µ)) generations (which happens with probability at most e-m + O(1/n4 )), then we repeat all the arguments from above. The expected number of such periods is 1/(1-e-m - O(1/n4 )) = O(1), concluding the proof of Theorem 2.

4 An Upper Bound with Crossover

Now we show that NSGA-III with crossover can be much more efficient when optimizing m-OJZJk than without, especially for large k. The reason is that a pair x, y with |xj |1 = |yj |1 = 2n/m-k and maximum Hamming distance H(xj , yj ) = 2k can be created via mutation within k generations, and then a recombination of those xj and yj leads to the all one string in block j with probability Ω(1/4k ). This probability is by a factor nk /4k larger than the probability to cross the fitness valley in block j via mutation.

Theorem 7. Consider NSGA-III on f := m-OJZJk for 2 ≤ k ≤ n/(2m)+1, (1+2n/m)m/2 ≤ µ/2 ≤ no(n) , crossover probability pc ∈ (0, 1), and a number m of objectives with 2 ≤ m ≤ n/2. Further, assume the same conditions as in Lemma 1, and let δ > 0. Then, with probability at least 1 - e-(δ-3)m - O(1/n4 ), the number of generations until the whole Pareto front F is covered is at most

The expected number of generations is

Proof. We closely follow the ideas used in the proof of Theorem 2. For the initialization, with probability 1 - e-Ω(n2 ) , a single generation suffices to initialize an individual x with k ≤ |xj |1 ≤ 2n/m-k for all j ∈ [m/2]. The expected number of generations is also 1+o(1). Now we define Step j and the corresponding set Bj as in the proof of Theorem 1. Let Xj denote the number of generations until, for every v ∈ Bj , there exists a search point x∗ ∈ Pt with f(x∗ ) = v. After these Xj generations, let Z cross j be the number of additional generations until a desired z is created. As in Lemma 4, we see that

Now we estimate Z. To estimate Zj , we consider k + 1 consecutive generations. In the first k generations, we create two suitable individuals y1, x2 with maximum Hamming distance in block j which can then be used for recombination in the (k + 1)-th generation.

Lemma 8. For j ∈ [m/2] the random variable Zj is stochastically dominated by (k + 1)Z cross j where Z cross j is a geometrically distributed random variable with success probability σk := 1 10k pc/(2·4k+1 µ) 1+pc/(2·4k+1 µ) ∏kℓ=1 (1-pc)ℓ/(2emn) 1+(1-pc)·ℓ/(2emn) . Additionally, for δ > 1,

Proof. As in the proof of Lemma 4, for the desired search point z to create, we may assume that zj = 1m/2 . Let s := rf(z). Consider a sequence of k + 1 generations as follows: Within the first k generations, one may create a suitable pair of individuals x1, x2 such that |xj1 |1 = |xj2 |1 = 2n/m-k and the Hamming distance in block j is maximized, particularly H(xj1 , xj2 ) = 2k, while the individuals coincide in all other blocks. Moreover, for all i < j, we have (rf(x1))i = (rf(x2))i = si and we have (rf(x1))i = (rf(x2))i = 0 for i > j. In the (k+1)-th generation, perform uniform crossover on x1 and x2, and do not flip any bit during mutation, to create a suitable search point z with zj = 12n/m, particularly (rf(z))i = si for all i ≤ j, and (rf(z))i = 0 for all i > j. So for generation ℓ ≤ k, assume that there are two individuals xℓ,1, xℓ,2 with |xjℓ,1 |1 = |xjℓ,2 |1 = 2n/m - k, H(xjℓ,1 , xjℓ,2 ) = 2(ℓ - 1), where for ℓ = 1, the individuals xℓ,1, xℓ,2 are the same. They also coincide in all the other blocks, where (rf(xℓ,1))i = si for i < j, and (rf(xℓ,1))i = 0 for i > j is satisfied. Then in one trial, one may choose x1 as a parent (prob. at least 1/µ), omit the crossover, and flip exactly two bits, namely, one from 1 to 0 and one from 0 to 1 in block j, both of which do not contribute to H(x1, x2). Note that there are k - (ℓ - 1) such zeros, and 2n/m - 2(ℓ - 1) such ones. Hence, this happens with probability at least (1-pc)(k-ℓ+1)/n·(2n/m-(ℓ-1)-k)/n·(1-1/n)n-2 ≥ (1-pc)(k-ℓ+1)/n·(2n/m-2k+1)/(en) ≥ (1-pc)(k-ℓ+ 1)/n· (n/m-1)/(en) ≥ (1-pc)(k-ℓ+1)/(2emn) =: qℓ, and in one generation with probability at least

These x1 and x2 then both survive with probability at least 1/10: If |F 1 t | ≤ µ where |F 1 t | is the number of individuals with rank one, then both survive with probability 1. Otherwise, µ < |F 1 t | ≤ 2µ, and they survive with probability at least (when both are not selected within the first µ/2 iterations of the reference point scheme in Algorithm 2) (µ/2 2 )/(|F 1 t |-µ/2 2 ) ≥ (µ/2 2 )/(3µ/2 2 ) ≥ 1/10 for n sufficiently large. Hence, in a sequence of k generations, the probability that there is such a pair x1, x2 of individuals is at least 1/10k ∏kℓ=1 (qℓ/2)/(1 + qℓ/2). For the subsequent (k + 1)- th generation, there are two such individuals x1, x2. Now, choose x1, x2 as parents (prob. at least 1/µ2 ), perform uniform crossover to create y′ with (rf(y′ ))i = si for all i ≤ j, particularly (y′ )i = 12n/m and omit mutation (prob. pc/4k · (1 - 1/n)n ≥ pc/4k+1 ). Hence, such a desired y′ is generated in this (k + 1)-th generation with probability at least

which survives by Lemma 1. Hence, the probability is at least σk = pcross 10k · ∏kℓ=1 qℓ/2 1+qℓ/2 to create the desired y′ within k+1 generations. Hence, the random variable Z = ∑m/2 j=1 Zj, is indeed stochastically dominated by (k + 1)Z ∗ , where Z ∗ = ∑m/2 j=1 Z ∗ j is an independent sum of geometrically distributed random variables Z ∗ j with success probability σk. Its expectation is m/(2σk). Hence, we have that

for all λ ≥ 0 where d := m/(2σ2 k ). For λ = 4mδ/σk, where δ ≥ 1 is arbitrary, we observe that Pr(Z ≥ (k + 1) · (1 + 8δ)m/(2σk)) ≤ e-δm proving the lemma.

Now we take a union bound to obtain with Lemmas 4 and 8 for W := X + Z

Hence, after O(nmln(n/m)/(1 - pc) + kδm40k µ/pc · (2emn)k /(k!(1-pc)k )) = O(kδmk+1nk (80e)k µ/(pck!(1pc)k )) generations, we see that a desired x with rf(x) = s is created with probability at most (2n/m)-4m + e-δm. A union bound over at most 3m/2 choices of x with different rf(x) values shows that, with probability at most e-(3-δ)m + O(1/n4 ), after this number of generations there exists for each s ∈ -1, 0, 1 an x with rf(x) = s. Conditioned on this event, the final covering step proceeds as in Theorem 2 and succeeds with probability at least 1 -O(1/n4 ) within an additional O(m2n ln(n/m)/(1 - pc)) generations. Hence, with probability at least 1-e-(δ-3)m -O(1/n4 ), the whole Pareto front is covered in O(kδmk+1nk (80e)k µ/(pck!(1 - pc)k )) generations. If this does not happen, we repeat the above arguments to obtain an expected number of O(kmk+1nk (80e)k µ/(pck!(1 - pc)k )) generations, where we use δ = 4, proving the theorem.

5 A Lower Bound for 4-OJZJ

To complement the upper bound for the crossover case, we also establish a lower bound for mutation only in the case m = 4 and for population sizes asymptotically smaller than nk by a factor of n2 .

Theorem 9. Consider NSGA-III on f := 4-OJZJk for 2 ≤ k ≤ n/8 + 1, µ ∈ O(nk-2 ), where crossover is turned off. Then the expected number of generations until the whole Pareto front is covered is at least Ω(nk /µ).

Proof. By a classical Chernoff bound, with probability 1 - e-Ω(n) = 1 - o(1), each individual x satisfies |xj |1 ∈ {k, . . . , n/2 - k} for both j ∈ {1, 2} after initialization. Suppose that this happens. Then in all future populations P1, P2, . . . there are only individuals x with |xj |1 ∈ {0, k, . . . , n/2 - k, n/2} for at least one j ∈ {1, 2}, since other individuals are dominated by every y with |yj |1 ∈ {k . . . , n/2 - k} for both j ∈ {1, 2}.

Lemma 10. With probability at least 1/81 - o(1), we have that Ct := {s ∈ {-1, 0, 1}2 | there is x ∈ Pt with rf(x) = s} ̸= {-1, 0, 1}2 after cn ln(n) generations for a constant c > 0, and Ct can only increase by flipping either k specific ones or zeros in a block j ∈ {1, 2}.

Assume the event from Lemma 10 occurs. Then flipping k specified ones or zeros within a block happens with probability at most 2n-k per trial, so the expected number of trials is at least (1 - o(1))nk /2 = Ω(nk ). This yields an expected number of Ω(nk /µ) generations.

6 Conclusion

In this paper, we derived upper runtime bounds for the widely used NSGA-III algorithm on m-OJZJk with an arbitrary number m of objectives, both with and without crossover, in expectation and with high probability. Notably, crossover can significantly speed up the runtime and even lead to exponential improvements in suitable parameter regimes, for example when the gap parameter k is large while the number of objectives and the population size are small. To complement these findings, we also established a lower runtime bound for the case m = 4. We hope these results improve the understanding of crossover in MOEAs, especially for escaping local optima and improving search efficiency. Future work could investigate lower bounds for the crossover case, and for m > 4 if crossover is turned off. The latter is more challenging due to population dynamics, and we still lack a clear understanding of how non-Pareto-optimal individuals behave in the population. For crossover, the challenge is that knowledge about objective vectors alone are insufficient. One has to know also about the distribution of the genotypes of individuals. Another direction is to explore crossover in more complex many-objective settings, such as combinatorial optimization problems like shortest paths, or permutation spaces, like the many-objective flow-shop scheduling problem.