# Rank-based Non-dominated Sorting

Bogdan Burlacu  University of Applied Sciences Upper Austia  Heuristic and
Evolutionary Algorithms Laboratory  Josef Ressel Centre for Symbolic
Regression

###### Abstract

Non-dominated sorting is a computational bottleneck in Pareto-based multi-
objective evolutionary algorithms (MOEAs) due to the runtime-intensive
comparison operations involved in establishing dominance relationships between
solution candidates. In this paper we introduce Rank Sort, a non-dominated
sorting approach exploiting sorting stability and ordinal information to avoid
expensive dominance comparisons in the rank assignment phase. Two algorithmic
variants are proposed: the first one, _RankOrdinal_ (RO), uses ordinal rank
comparisons in order to determine dominance and requires $O(N)$ space; the
second one, _RankIntersect_ (RS), uses set intersections and bit-level
parallelism and requires $O(N^{2})$ space. We demonstrate the efficiency of
the proposed methods in comparison with other state of the art algorithms in
empirical simulations using the NSGA2 algorithm as well as synthetic
benchmarks. The _RankIntersect_ algorithm is able to significantly outperform
the current state of the art offering up to 30% speed-up for many objectives.
C++ implementations are provided for all algorithms.

##  1 Introduction

Multiobjective optimization problems (MOP) consist of finding the best
compromise or trade-off between a collection of competing objectives.
Formally, a MOP can be defined as:

$\text{Minimize}\leavevmode\nobreak\ F(X)\text{ subject to }x\in X$ (1)

where $F:X\to Y$, $X\in\mathbb{R}^{n}$ is the _decision variable space_ (also
called a _feasible set_) and $Y\in\mathbb{R}^{m}$ is the _objective function
value space_.

Pareto optimality is an established formalism used to approach this type of
problem, based on the concept of dominance. Given two vectors
$y,z\in\mathbb{R}^{m}$, it defines the following dominance relationships:

equality $\displaystyle y=z\Leftrightarrow y_{i}=z_{i},\forall i=1,...,m$ (2)
strong dominance $\displaystyle y\prec z\Leftrightarrow y_{i}<z_{i},\forall
i=1,...,m$ (3) weak dominance $\displaystyle y\preceq z\Leftrightarrow
y_{i}\leq z_{i},\forall i=1,...,m\leavevmode\nobreak\
\text{and}\leavevmode\nobreak\ y\neq z$ (4)

Thus, if $y\preceq z$ then we say that _$y$ dominates $z$_ (_$z$ is dominated
by $y$_). If $y\npreceq z$ and $z\npreceq y$ then $y$ and $z$ are _mutually
non-dominated_. A _Pareto front_ represents a set of mutually non-dominated
solutions.

Note that $\prec$ is a _strict partial order_ and $\preceq$ is a _partial
order_ on the objective space $R^{m}$. This allows to more generally describe
Pareto relationships in terms of “predecessors” and “successors” as will be
needed later on when Rank Sort is introduced.

When optimizing multiple objectives, it is often the case that a solution
cannot be further optimized with respect to any objective without worsening at
least one of the others. The concept of _Pareto optimality_ represents an
intuitive interpretation of this situation. If $x^{*}$ is a _Pareto optimal_
solution, then there is no other solution $x\in X$ such that $x\preceq x^{*}$.

Evolutionary algorithms (EA) are one of the most popular methods for solving
multi-objective optimization problems (MOP). In Pareto-based MOEAs, non-
dominated sorting (NDS) refers to the process of dividing the population into
Pareto fronts which are subsequently used to guide parent selection.

In this paper we present an approach to non-dominated sorting which can
substantially reduce the runtime costs of the parent selection step. This is
motivated by the fact that, due to its transactional nature and the need to
carry rank information from one individual to the next, non-dominated sorting
does not lend itself well to parallelization and can substantially increase
the runtime of MOEAs.

The remainder of the paper is organized as follows: Section 2 gives an
overview of existing algorithms from the literature. Section 3 discusses the
proposed approach, Section 4 provides experimental results and Section 5 is
dedicated to conclusions.

##  2 Related work

Improvements in the asymptotic complexity of non-dominated sorting algorithms
are mainly driven by new techniques to reduce (or eliminate altogether) the
number of Pareto dominance calculations necessary to divide the population
into Pareto fronts.

We distinguish between two main categories of approaches: those that perform
direct dominance comparisons and those that perform sorting in order to
determine the Pareto ranking of the population.

###  2.1 Inference-based approaches

Inference-based methods avoid unnecessaneuripsry dominance comparisons by
exploiting the transitive property of dominance relationships (e.g., if
$a\prec b$ and $b\prec c$, then $a\prec c$) in order to improve the asymptotic
complexity and runtime performance.

__Deductive sort (DS)__ (McClymont and Keedwell, 2012) iterates over the
solutions in natural order and maintains a strict order of assessment and
comparison to ensure that all dominated solutions are discarded and not
incorrectly added to the wrong front. Once a front has been filled, all
solutions assigned to that front are ignored and the process is repeated until
all solutions are assigned to a front.

__Efficient non-dominated sorting__ (Zhang et al., 2015) avoids unnecessary
dominance comparisons by comparing a solution that needs to be assigned to a
front only with the other solutions that have already been assigned. Two
algorithmic variants are proposed, differing in the search strategy used to
find if any solution in a given front dominates the current solution. Two
search strategies are proposed: sequential search (ENS-SS) and binary search
(ENS-BS).

_Hierarchical Non-dominated Sorting (HS)_ (Bao et al., 2017) sorts all
solutions lexicographically and then uses successive rounds of comparison
between the first solution and the succeeding solutions to establish dominance
relationships. Dominated solutions are discarded in the current round of
comparison so they will not be compared again with the first solution or other
solutions which are non-dominated with the first solution. Once a solution is
ranked it is removed from the comparison and added to its respective front.

###  2.2 Sort-based approaches

In contrast with inference-based approaches where dominance is explicitly
considered, sort-based approaches rely on ordinal ranking information (in some
form or another) in order to determine dominance relationships between
solutions. This typically involves a more expensive preprocessing step where
the solutions are sorted w.r.t. every objective.

The most notable examples of sort-based approaches are Best Order Sort (BOS)
(Roy et al., 2016) and Merge Non-dominated Sort (MNDS) (Moreno et al., 2020).
Both approaches share the same principneuripsle of using a solution’s
_dominance set_ (i.e., the set of solutions dominating the current solution)
to update its rank. The dominance set is obtained as the intersection between
the sets of solutions ranked better than the current solution according to
each objective. The ranking is obtained by stable sorting – a key
characteristic which is also employed by our algorithm.

A key difference between BOS and MNDS is given by the actual implementation of
the set intersection procedure for obtaining the dominance set. BOS uses a
linked list to store a solution’s dominance set while MNSE uses a bitset. Both
algorithms are deeply integrated with merge sort which they use in the sorting
step and employ various other auxiliary data structures to keep track of
already ranked solutions or duplicates (MNDS).

__Best Order Sort (BOS)__ (Roy et al., 2016)111Implementation available at
https://github.com/proteekroy/Best-Order-Sort sorts the population according
to each objective preserving lexicographical order. The resulting ranking is
used to compute the dominance set for each solution. Finally, a solution’s
rank becomes the rank of the worst-ranked solution in the dominance set plus
one.

__Merge Non-dominated Sort (MNDS)__ (Moreno et al., 2020)222Available in the
open-source framework jMetal (https://github.com/jMetal/jMetal). Employs a
stable-sorting algorithm (merge sort) to rank solutions according to each
objective and uses the ranking to construct a solution’s dominance set as the
set of solutions that dominate the current solution in every objective. Then,
the rank of the current solution is updated according to the worst rank in the
dominance set. Prior to rank assignment, the algorithm removes duplicate
solutions from the population and keeps track of them in a separate list.
Finally, the ranks of the solutions are obtained based on the dominance sets
and duplicates are inserted again with their corresponding rank.

A summary of the runtime complexity and space requirements of the described
algorithms is given in Table 1.

Table 1: Summary of algorithm time and space complexity. For DS, the worst case is $O(MN^{3})$ cf. Mishra and Buzdalov Mishra and Buzdalov (2020). Algorithm | Best | Worst | Space
---|---|---|---
DS | $O(MN^{2})$ | $O(MN^{3})$ | $O(N)$
HS | $O(MN\sqrt{N})$ | $O(MN^{2})$ | $O(N)$
MNDS | $O(MN\log N)$ | $O(MN^{2})$ | $O(N^{2})$
ENS-SS | $O(MN\sqrt{N})$ | $O(MN^{2})$ | $O(1)$
ENS-BS | $O(MN\log N)$ | $O(MN^{2})$ | $O(1)$
RO | $O(MN\log N)$ | $O(MN^{2})$ | $O(N)$
RS | $O(MN\log N)$ | $O(MN^{2})$ | $O(N^{2})$

##  3 Rank-based Non-dominated Sorting

Similar to MNDS and BOS, the Rank Sort algorithm uses the concept of a
“dominance set” to compute solution ranks. However, while MNDS and BOS define
the dominance set as the set of predecessors of the current solution, Rank
Sort defines the dominance set as the set of _successors_ of the current
solution. Here, the notions of “predecessor” and “successor“ correspond to the
$\preceq$ and $\succeq$ Pareto relationships, respectively.

Considering the dominance set as the set of successors brings an important
advantage: during rank update, it is sufficient to consider _equally ranked
successors_ of the current solution (considering provisional ranks iteratively
computed by the algorithm). Consequently, the rank update becomes a simple
_increment operation_ of the successor’s rank. This insight simplifies the
rank update phase and reduces the number of operations such that the total
number of rank updates will be equal to the sum of ranks in the final Pareto
front assignment.

Due to this rank update rule, the algorithm’s behavior can be interpreted as
establishing each new Pareto front by moving dominated individuals from the
current Pareto front into the new one.

We present two variants that differ in their space complexity:

  * •

_RankOrdinal_ (RO) requires $O(N)$ space by avoiding set intersections.
Instead, it takes the smallest set of successors w.r.t. any objective
$k=1,...,M$ as the dominance set:

$D_{\text{ordinal}}(S_{i})=D_{k}(S_{i}),k=\operatorname*{arg\,min}_{l}|D_{l}(S_{i})|$
(5)

  * •

_RankIntersect_ (RS) requires $O(N^{2})$ space for storage of dominance sets
and employs bit-level parallelism to speed up the computation of set
intersections. The dominance set is defined as:

$D_{\text{intersect}}(S_{i})=\bigcap D_{k}(S_{i}),k=1,...,M$ (6)

Here, $D_{k}(S_{i})$ represents the set of successors w.r.t. objective $k$:

$neuripsD_{k}(S_{i})=\\{S_{j}\leavevmode\nobreak\ |\leavevmode\nobreak\
\text{rank}(S_{i})=\text{rank}(S_{j})\text{ and }S_{j}\succeq S_{i}\\}$ (7)

Node that Equation 7 should be taken in the context of iterative rank
assignment within the algorithm’s inner loop, where the ranks are not yet
final.

_RankIntersect_ is the faster algorithm overall however its space requirements
might make it less suitable for resource-constrained environments, where
_RankOrdinal_ might be preferable.

###  3.1 RankOrdinal Algorithm

_RankOrdinal_ (RO) uses the comparison of ordinal ranks to establish
dominance. Let $i\in 1,...,N$ be an index over the solutions $S_{i}$, let
$k\in 1,...,M$ be an index over the objectives, and let $F_{k}(i)$ be a
function returning the value of the $k$-th objective for solution $S_{i}$.

For every objective $k$, we compute the permutation vector $\mathbf{p}_{k}$
and ordinal rank vector $\mathbf{r}_{k}$:

$\displaystyle\mathbf{p}_{k}$
$\displaystyle=\big{[}p_{k}(1)\leavevmode\nobreak\ \ldots\leavevmode\nobreak\
p_{k}(N)\big{]}^{T}\leavevmode\nobreak\ \text{\small such
that}\leavevmode\nobreak\ F_{k}(i)\leq F_{k}(i+1)$ (8)
$\displaystyle\mathbf{r}_{k}$
$\displaystyle=\big{[}r_{k}(1)\leavevmode\nobreak\ \ldots\leavevmode\nobreak\
r_{k}(N)\big{]}^{T}\leavevmode\nobreak\ \text{such that}\leavevmode\nobreak\
r_{k}\left(p_{k}(i)\right)=i$ (9)

Let $\mathbf{P}=[\mathbf{p}_{k}^{T}]\in\mathbb{N}^{N\times M}$ and
$\mathbf{R}=[\mathbf{r}_{k}^{T}]\in\mathbb{N}^{N\times M}$ be the permutation
matrix and ordinal rank matrix, where $k=1,...,M$.

It is clear from Eq. 8 that for every $\mathbf{p}_{k}\in\mathbf{P}$, the
leftmost solution $S_{p_{k}(1)}$ is always non-dominated and in general,
$S_{p_{k}(i)}\nsucceq S_{p_{k}(j)}$ for every $i<j$. Conversely, a solution
$S_{p_{k}(i)}$ can only dominate other solutions that come after it in the
permutation. Therefore, it will be efficient to examine solutions left-to-
right in the partial permutation where $S_{p_{k}(i)}$ is farthest from the
left, such that there remain fewer successors to be examined.

By virtue of Eqs. 8, 9, the matrix $\mathbf{R}$ helps infer a dominance
relationship between two solutions:

$\displaystyle\text{$S_{i}$ non-dominated with $S_{j}$},\leavevmode\nobreak\
\text{if}\leavevmode\nobreak\ r_{k}(i)>r_{k}(j),\exists k=1,...,M$ (10)
$\displaystyle S_{i}\preceq S_{j},\leavevmode\nobreak\
\text{if}\leavevmode\nobreak\ r_{k}(i)<r_{k}(j),\forall k=1,...,M$ (11)

Non-dominated sorting algorithms assign solutions to their corresponding
Pareto fronts based on the _domination rank_ value (see for example (Deb et
al., 2002)). Testing for dominance in permutation order gives our approach the
advantage of a natural rank update mechanism.

$S_{i}\preceq S_{j}\text{ and
}\text{rank}(i)=\text{rank}(j)\implies\text{rank}(j)\leftarrow\text{rank}(i)+1$
(12)

When the rank of $S_{j}$ is updated, the rank of $S_{i}$ will have already
been updated at a previous iteration. The entire algorithm is described in
pseudocode in Algorithm 1. The outer loop iterates over solution indices given
by $\mathbf{p}_{1}$ (the first column of $\mathbf{P}$). The inner loop of the
algorithm will iterate over the smallest set of successors, as defined by 5.
The corresponding column of $\mathbf{P}$ and the start index for the partial
permutation are identified at line 6 in Algorithm 1. The rank assignment
procedure is performed at lines 9–11. The comparison of ranks $r_{k}$ is
actually an element-wise comparison between the columns of matrix $\mathbf{R}$
corresponding to the two solutions (fast in practice due to vectorization and
the fact that the elements are integers).

####  3.1.1 Solution equality

Note that due to its reliance on sorting, ordinal ranking is unable to detect
the situation when two solutions are equal in all objectives. To handle this
case it is necessary to extend for example Eqs. 11 and 12 to explicitly check
for equality. This can be cheaply implemented by e.g. hashing objective values
in the preprocessing phase. However, it is a better design choice to separate
duplicates handling logic from non-dominated sorting.

Input: $S_{1},...,S_{N}$

Output: Pareto fronts $F_{1},F_{2},...$

1 Function  _RankOrdinal(_$S_{1},...,S_{N}$_)_:

2 $\mathbf{P}\leftarrow$ StableSort(_$S_{1},...,S_{N}$_);

3 $\mathbf{R}\leftarrow$ calculate using Eq. 9;

4 rank $\leftarrow$ zero-initialized array of size $N$;

5 foreach  _$S_{i}\in\mathbf{p}_{1}$_ do

6 $k\leftarrow\operatorname*{arg\,max}_{l}r_{l}(i)$;

7 if  _$r_{k}(i)=N$_ then

continue;

// $S_{i}$ cannot dominate any other

8

9 foreach  _$S_{j}\in\mathbf{p}_{k}$_ do

10 if  _$\text{rank}(i)=\text{rank}(j)$ and $S_{i}\preceq S_{j}$ cf. Eq 11_
then

11 $\text{rank}(j)\leftarrow\text{rank}(j)+1$;

12

13

14

15 $\mathbf{F}\leftarrow$ list of fronts of size $\max_{i}\text{rank}(i$)+1;

16 foreach  _$i$ from $1$ to $N$_ do

17 append $i$ to front $F_{\text{rank}(i)}$;

18

19 return $\mathbf{F}$;

20

21

Algorithm 1 _RankOrdinal_

####  3.1.2 Correctness of the algorithm

The permutations of solutions w.r.t. every objective are generated by calling
a sorting procedure. It is important to note that the algorithm is correct
only when the sorting procedure is stable, as otherwise Eqs. 10, 11 would not
be reliable. This issue is apparent in certain extended definitions of
dominance such as $\varepsilon$-dominance (Laumanns et al., 2002) where in the
absence of stability there can be no guarantee about the order in which
$\varepsilon$-dominated solutions are added to their respective fronts.

Many non-dominated sorting algorithms explicitly make use of sorting in order
to avoid the cost of dominance comparisons in at least one dimension (Zhang et
al., 2015; Roy et al., 2016; Bao et al., 2017; Zhou et al., 2017; Xue et al.,
2020; Moreno et al., 2020). However, only a few explicitly consider stability
(Roy et al., 2016; Moreno et al., 2020). While this is almost never a problem
in empirical runs, it can be a potential source of inconsistency in MOEAs, for
example if two equal solution candidates are placed on the Pareto front(s) in
different relative order due to sorting instability. Therefore, in some
particular cases, use of an unstable stable sorting procedure may introduce an
undesirable dependency to an implementation detail that can cause unexplained
behavior or an inability to reproduce results (i.e., runs with a fixed seed).
This was also observed by (Buzdalov, 2018) who called this phenomenon “bug-
compatible” non-dominated sorting.

####  3.1.3 Computational complexity

Since the preprocessing part has a fixed cost (one lexicographic sorting step
followed by $M-1$ regular sorting steps, its complexity is $O(MN\log N$). For
the rank assignment part of the algorithm, we first consider the total number
of iterations performed. Regardless how little work is done per iteration
(e.g., when the conditions allow to skip over an already ranked solution),
this quantity is important in determining the final complexity of the
algorithm.

As shown in Algorithm 1, the rank assignment procedure consists of two nested
loops: the outer loop iterates over all the solutions $S_{i}$ in the order
given by $\mathbf{p}_{1}$, while the inner loop iterates over the successors
of $S_{i}$ in the permutation $\mathbf{p}_{k}$ where
$k=\operatorname*{arg\,max}_{l}r_{l}(i)$.

#### Worst case

The worst case occurs when the input is degenerate, namely when the points to
be sorted either belong to a single front or belong to individual fronts.

  * •

When the solutions are non-dominated and belong to a single Pareto front, this
corresponds to the unique situation where matrix $\mathbf{P}$ has the form:

$\mathbf{P}=\begin{bmatrix}1&2&...&N\\\ N&N-1&...&1\end{bmatrix}^{T}$ (13)

In this case, the size of dominance set $D_{\text{ordinal}}(S_{i})$ will be:

$|D_{\text{ordinal}}(S_{i})|=\min(i,N-i)$ (14)

Therefore, the total number of iterations will be:

$\text{total
iterations}=2\cdot\sum_{i=1}^{N/2}i=2\cdot\frac{\frac{N}{2}\left(\frac{N}{2}+1\right)}{2}\approx\frac{N^{2}}{4}$
(15)

  * •

When each solution $S_{i}$ belongs to its own unique front, this corresponds
to the unique situation where matrix $\mathbf{P}$ has the form:

$\mathbf{P}=\begin{bmatrix}1&2&...&N\\\ 1&2&...&N\end{bmatrix}^{T}$ (16)

In this case, the size of dominance set $D_{\text{ordinal}}(S_{i})$ will be:

$|D_{\text{ordinal}}(S_{i})|=N-i$ (17)

The total number of iterations in the worst case will be the sum:

$\text{total
iterations}=\sum_{i=1}^{N}i=\frac{N(N+1)}{2}\approx\frac{N^{2}}{2}$ (18)

Note that the results above hold for $M>2$ since adding more columns to
$\mathbf{P}$ would not influence the resulting fronts because of the
degenerate configuration of the first two columns. Therefore, _RankOrdinal_
has the worst-case complexity of $O(MN^{2})$.

#### Best case

In the following we make the argument that the best-case complexity of
_RankOrdinal_ is $O(MN\log N)$. We consider the problem of finding the Pareto
set of randomly selected points in the unit hypercube.

  * •

A full dominance check is always precluded by the condition
$\text{rank(i)}=\text{rank}(j)$ in Algorithm 1 line 10. Since the algorithm
works by incrementing the ranks of dominated solutions within the current
Pareto front (in order to “move” them to a new front), this condition will be
true a finite number of times, and this number is proportional to the size of
the Pareto fronts.

  * •

Pareto front size is bounded in expectation by generalized harmonic numbers
and thus scales in proportion to the natural logarithm (Yukish et al., 2004),
(Köppen et al., 2005). Therefore, the number of time the condition
$\text{rank}(i)=\text{rank}(j)$ is true will also be proportional to the
logarithm of $N$, bounding the dominance checks.

Additionally, it can be shown that the expected minimum of $M$ uniform i.i.d.
random variables is

$\mathbb{E}\left[\min(X_{i})\right]=\frac{1}{M+1}$ (19)

For random objective values sampled from the unit hypercube, the total number
of iterations (regardless of the amount of work done inside the inner loop and
without considering the rank equality condition) will be $\frac{N^{2}}{M+1}$.

Figure 1, determined empirically by taking 1000-run averages of inner loop
iterations and dominance checks of _RankOrdinal_ , shows that the number of
dominance comparisons grows much slower than the total number of iterations.

Figure 1: RankOrdinal total number of iterations and total number of rank
comparisons

####  3.1.4 Example

We consider the case $N=10,\leavevmode\nobreak\ M=2$ with the following
points:

$\displaystyle S_{1}$ $\displaystyle=\begin{bmatrix}0.79&0.35\end{bmatrix}$
$\displaystyle S_{2}$ $\displaystyle=\begin{bmatrix}0.40&0.71\end{bmatrix}$
$\displaystyle S_{3}$ $\displaystyle=\begin{bmatrix}0.15&0.014\end{bmatrix}$
$\displaystyle S_{4}$ $\displaystyle=\begin{bmatrix}0.46&0.82\end{bmatrix}$
$\displaystyle S_{5}$ $\displaystyle=\begin{bmatrix}0.28&0.98\end{bmatrix}$
$\displaystyle S_{6}$ $\displaystyle=\begin{bmatrix}0.31&0.74\end{bmatrix}$
$\displaystyle S_{7}$ $\displaystyle=\begin{bmatrix}0.82&0.52\end{bmatrix}$
$\displaystyle S_{8}$ $\displaystyle=\begin{bmatrix}0.84&0.19\end{bmatrix}$
$\displaystyle S_{9}$ $\displaystyle=\begin{bmatrix}0.85&0.78\end{bmatrix}$
$\displaystyle S_{10}$ $\displaystyle=\begin{bmatrix}0.96&0.83\end{bmatrix}$

The following permutations are obtained by stable-sorting according to the two
objectives:

$\displaystyle\mathbf{p}_{1}$
$\displaystyle=\begin{bmatrix}3&5&6&2&4&1&7&8&9&10\end{bmatrix}^{T}$
$\displaystyle\mathbf{p}_{2}$
$\displaystyle=\begin{bmatrix}3&8&1&7&2&6&9&4&10&5\end{bmatrix}^{T}$

The corresponding matrices $\mathbf{P}$ and $\mathbf{R}$ are:

$\displaystyle\mathbf{P}$
$\displaystyle=\begin{bmatrix}3&5&6&2&4&1&7&8&9&10\\\
3&8&1&7&2&6&9&4&10&5\end{bmatrix}^{T}$ $\displaystyle\mathbf{R}$
$\displaystyle=\begin{bmatrix}6&4&1&5&2&3&7&8&9&10\\\
3&5&1&8&10&6&4&2&7&9\end{bmatrix}^{T}$

The algorithm will iterate over $\mathbf{p}_{0}$ and assert dominance between
each solution $S_{i}$ (highlighted with a circle at each iteration) and its
successors. The successors $S_{j}$ that are examined in the inner loop are
highlighted with an underline. The conditions which apply at each inner loop
iteration are encoded as follows:

$S_{j}^{\star}$

$\to$ $S_{i}\preceq S_{j}$ determined after a full rank comparison.

$S_{j}^{\bullet}$

$\to$ $S_{i}\npreceq S_{j}$ determined after a full rank comparison.

$S_{j}^{\square}$

$\to$ $S_{j}$ skipped due to $\text{rank}(i)\neq\text{rank}(j)$ (Alg. 1 line
10)

The progress of the algorithm is illustrated in Figure 2(a). The resulting
fronts are:
$F_{1}=\\{3\\},F_{2}=\\{1,2,5,6,8\\},F_{3}=\\{4,7\\},F_{4}=\\{9\\},F_{5}=\\{10\\}$.

Figure 2: Step-by-step operation of the RankOrdinal and RankIntersect
algorithms

# | Outer iteration $S_{i}$ (circled), inner iteration $S_{j}$ (underlined) | Domination ranks
---|---|---
1 | 3 | 5 | 6 | 2 | 4 | 1 | 7 | 8 | 9 | 10 | 2 | 2 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2
| 3 |  8⋆ |  1⋆ |  7⋆ |  2⋆ |  6⋆ |  9⋆ |  4⋆ |  10⋆ |  5⋆ |  |  |  |  |  |  |  |  |  |
2 | 3 | 5 | 6 | 2 | 4 | 1 | 7 | 8 | 9 | 10 | 2 | 2 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2
| 3 | 8 | 1 | 7 | 2 | 6 | 9 | 4 | 10 | 5 |  |  |  |  |  |  |  |  |  |
3 | 3 | 5 | 6 | 2 | 4 | 1 | 7 | 8 | 9 | 10 | 2 | 2 | 1 | 3 | 2 | 2 | 2 | 2 | 3 | 3
| 3 | 8 | 1 | 7 | 2 | 6 |  9⋆ |  4⋆ |  10⋆ |  5∙ |  |  |  |  |  |  |  |  |  |
4 | 3 | 5 | 6 | 2 | 4 | 1 | 7 | 8 | 9 | 10 | 2 | 2 | 1 | 3 | 2 | 2 | 2 | 2 | 3 | 3
| 3 | 8 | 1 | 7 | 2 |  6∙ |  9□ |  4□ |  10□ |  5∙ |  |  |  |  |  |  |  |  |  |
5 | 3 | 5 | 6 | 2 | 4 | 1 | 7 | 8 | 9 | 10 | 2 | 2 | 1 | 3 | 2 | 2 | 2 | 2 | 3 | 4
| 3 | 8 | 1 | 7 | 2 | 6 | 9 | 4 |  10⋆ |  5□ |  |  |  |  |  |  |  |  |  |
6 | 3 | 5 | 6 | 2 | 4 | 1 |  7⋆ |  8∙ |  9□ |  10□ | 2 | 2 | 1 | 3 | 2 | 2 | 3 | 2 | 3 | 4
| 3 | 8 | 1 | 7 | 2 | 6 | 9 | 4 | 10 | 5 |  |  |  |  |  |  |  |  |  |
7 | 3 | 5 | 6 | 2 | 4 | 1 | 7 |  8□ |  9⋆ |  10□ | 2 | 2 | 1 | 3 | 2 | 2 | 3 | 2 | 4 | 4
| 3 | 8 | 1 | 7 | 2 | 6 | 9 | 4 | 10 | 5 |  |  |  |  |  |  |  |  |  |
8 | 3 | 5 | 6 | 2 | 4 | 1 | 7 | 8 |  9□ |  10□ | 2 | 2 | 1 | 3 | 2 | 2 | 3 | 2 | 4 | 4
| 3 | 8 | 1 | 7 | 2 | 6 | 9 | 4 | 10 | 5 |  |  |  |  |  |  |  |  |  |
9 | 3 | 5 | 6 | 2 | 4 | 1 | 7 | 8 | 9 |  10⋆ | 2 | 2 | 1 | 3 | 2 | 2 | 3 | 2 | 4 | 5
| 3 | 8 | 1 | 7 | 2 | 6 | 9 | 4 | 10 | 5 |  |  |  |  |  |  |  |  |  |
(a) RankOrdinal Example

# | Perm. index | Rank | Dominance set | Rank set | Domination ranks
---|---|---|---|---|---
1 | $p_{2}(1)=3$ | rank$(S_{3})$ = 1 | $D(S_{3})=\\{1,2,4,5,6,7,8,9,10\\}$ | $r_{1}=\\{3\\}$ | 2 | 2 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2
2 | $p_{2}(2)=8$ | rank$(S_{8})$ = 2 | $D(S_{8})=\\{9,10\\}$ | $r_{2}=\\{1,2,4,5,6,7,9,10\\}$ | 2 | 2 | 1 | 2 | 2 | 2 | 2 | 2 | 3 | 3
3 | $p_{2}(3)=1$ | rank$(S_{1})$ = 2 | $D(S_{1})=\\{7\\}$ | $r_{2}=\\{1,2,4,5,6,8\\}$ | 2 | 2 | 1 | 2 | 2 | 2 | 3 | 2 | 3 | 3
4 | $p_{2}(4)=7$ | rank$(S_{7})$ = 3 | $D(S_{7})=\\{9,10\\}$ | $r_{3}=\\{7\\}$ | 2 | 2 | 1 | 2 | 2 | 2 | 3 | 2 | 4 | 4
5 | $p_{2}(5)=2$ | rank$(S_{2})$ = 2 | $D(S_{2})=\\{4\\}$ | $r_{2}=\\{1,2,5,6,8\\}$ | 2 | 2 | 1 | 3 | 2 | 2 | 3 | 2 | 4 | 4
6 | $p_{2}(6)=6$ | rank$(S_{6})$ = 2 | $D(S_{6})=\emptyset$ | $r_{2}=\\{1,2,5,6,8\\}$ | 2 | 2 | 1 | 3 | 2 | 2 | 3 | 2 | 4 | 4
7 | $p_{2}(7)=9$ | rank$(S_{9})$ = 4 | $D(S_{9})=\\{10\\}$ | $r_{4}=\\{9\\}$ | 2 | 2 | 1 | 3 | 2 | 2 | 3 | 2 | 4 | 5
8 | $p_{2}(8)=4$ | rank$(S_{4})$ = 3 | $D(S_{4})=\emptyset$ | $r_{3}=\\{4,7\\}$ | 2 | 2 | 1 | 3 | 2 | 2 | 3 | 2 | 4 | 5
9 | $p_{2}(9)=10$ | rank$(S_{10})$ = 5 | $D(S_{10})=\emptyset$ | $r_{5}=\\{10\\}$ | 2 | 2 | 1 | 3 | 2 | 2 | 3 | 2 | 4 | 5
10 | $p_{2}(10)=5$ | rank$(S_{5})$ = 2 | $D(S_{5})=\emptyset$ | $r_{2}=\\{1,2,5,6,8\\}$ | 2 | 2 | 1 | 3 | 2 | 2 | 3 | 2 | 4 | 5
(b) RankIntersect Example

###  3.2 RankIntersect Algorithm

The _RankIntersect_ (RS) algorithm takes a solution’s dominance set as the
intersection of its objective-wise dominance sets:

$\displaystyle D(S_{i})=\bigcap_{k}D_{k}(S_{i})$ (20)

The intersections are efficiently computed by exploiting bit-level
parallelism. In a bitset, integer values are encoded as positions of the set
bits. This has the advantage of reducing set intersections to simple logical
$\land$ operations that are optimized in the hardware.

The bitset approach is also employed by the MNDS Moreno et al. (2020), however
MNDS considers a dominance set as the set of solutions (or equivalently, their
permutation indices) that dominate the solution under consideration and
updates the considered solution’s rank using the maximum rank in the dominance
set (the same rank update mechanism is also shared by best order sort Roy et
al. (2016)). This leads to inefficiencies due to the overlap in dominance sets
and the need to compute maximum ranks.

According to Equation 7, for each solution $S_{i}$, _RankIntersect_ updates
the ranks of successors $S_{j}$ where the current rank of $S_{j}$ equals the
rank of $S_{i}$ and $S_{j}\succeq S_{i}$.

To implement this idea, _RankIntersect_ maintains a separate collection of
sets (bitsets) corresponding to each dominance rank value (from 1 to …), such
that each set contains the permutation indices of individuals with that rank.
This rank set is updated along with the rank updates and used to reduce each
individual’s dominance set according to Equation 7. The resulting algorithm is
illustrated in Algorithm 2:

  * •

The first part of the algorithm up to line 11 is responsible for preprocessing
and initialization of data structures. A work bitset $b$ is used to initialize
dominance sets $B(i)$ for each individual $i$.

  * •

Line 12 marks the beginning of the algorithmic loop where objective-wise
dominance sets are intersected. When the last objective is reached, the
algorithm proceeds to rank the rank assignment phase.

  * •

Lines 19–20 perform the bitset intersections including the additional
intersection with the corresponding rank set.

  * •

The dominance set consisting of equally ranked successors is obtained at line
21. These successors will have their rank incremented at line 25. Therefore,
they are removed from the current rank set (line 22) and added to the next one
(line 23).

The resulting algorithm is functionally equivalent with  _RankOrdinal_ whose
correctness has already been discussed. Relying on bitset operations requires
extra storage space for the bitsets but offers a speed advantage as bit-level
parallelism can accelerate intersection operations with a factor equal to the
basic data block size (i.e., a 64-bit integer used as a basic data block can
store 64 permutation indices), which can be further increased by
vectorization.

Input: $S_{1},...,S_{N}$

Output: Pareto fronts $F_{1},F_{2},...$

1 Function  _RankIntersect(_$S_{1},...,S_{N}$_)_:

2 $\mathbf{P}\leftarrow$ StableSort(_$S_{1},...,S_{N}$_);

3 $\mathbf{B}\leftarrow$ list of bitsets associated to each $S_{i}$;

4 $\mathbf{K}\leftarrow$ list of bitsets associated to each dominance rank
$r$;

5 $\mathbf{K}(1)\leftarrow\text{new bitset}$;

set all bits in $\mathbf{K}(1)$;

// all solutions start with rank 1

6 $b\leftarrow$ auxiliary bitset used to compute the set intersections;

7 set all bits in $b$;

8 foreach  _i from 1 to N_ do

9 reset bit $i$ in $b$;

10 $B(i)\leftarrow b$;

11

12

13 foreach  _$k$ from $2$ to $M$_ do

14 set all bits in $b$;

15 foreach  _$i$ in $\mathbf{p}_{k}$_ do

16 reset bit $i$ in $b$;

17

18 if  _$k <M$_ then

19 $\mathbf{B}(i)\leftarrow\mathbf{B}(i)\land b$;

20

21 else

$r\leftarrow\mathbf{K}\left(\text{rank}(i)\right)$;

// bitset for current rank

$s\leftarrow\mathbf{K}\left(\text{rank}(i+1)\right)$;

// bitset for next rank

22 $v\leftarrow\mathbf{B}(i)\land b\land r$;

$r\leftarrow r\land\lnot v$;

// remove dominance set from $r$

$s\leftarrow s\lor v$;

// add dominance set to $s$

23

24 foreach  _$j$ in $v$_ do

25 $\text{rank}(j)=\text{rank}(j)+1$;

26

27

28

29

30

31 $\mathbf{F}\leftarrow$ list of fronts of size $\max_{i}\text{rank}(i$)+1;

32 foreach  _$i$ from $1$ to $N$_ do

33 append $i$ to front $F_{\text{rank}(i)}$;

34

35 return $\mathbf{F}$;

36

37

Algorithm 2 _RankIntersect_

####  3.2.1 Computational complexity

Since _RankIntersect_ does not perform any dominance checks we discuss
computational complexity in terms of the costs of computing the necessary set
intersections.

#### Worst case

Logical operations $\land,\lor$ over bitsets are performed by iterating over
$\lceil\frac{N}{64}\rceil$ basic data blocks which are internally stored in a
list by the bitset. This is due to the fact that a data block is represented
by a 64-bit unsigned integer. The outer loop at line 14 performs $N$
iterations, while the inner loop performs $\lceil\frac{N}{64}\rceil$
iterations to intersect the bitsets. Therefore the total number of iterations
will be in the order of $N^{2}$ which leads to $O(MN^{2})$ worst-case
complexity.

#### Best case

Since Algorithm 2 successively intersects bitsets going from one objective to
the next (lines 12, 17), it becomes possible to memorize the regions at the
beginning and end of the list where the data blocks have become zero after the
intersection and then skip them in the next iteration. This creates a short-
circuiting mechanism for early stopping in the case of non-dominance. When all
solutions are non-dominated, as all the dominance sets will be empty, the
remaining complexity is $O(MN\log N)$ given by the sorting phase.

We conclude that the overall complexity of _RankIntersect_ is $O(MN^{2})$ in
the worst case and $O(MN\log N)$ in the best case. However, the advantage of
_RankSort_ comes from the fact that its operations are very basic and can
exploit the characteristics of modern processors. For $M>2$, most runtime
effort is spent performing objective-wise dominance set intersections at line
17. However in practice, a relatively modern CPU can intersect four basic data
blocks at a time using AVX2 instructions (since a “wide” SIMD type can hold
four 64-bit unsigned integers). For large $N$, the runtime effort shifts
towards the sorting phase, suggesting that additional runtime benefits might
be achieved by parallelizing the sorting.

####  3.2.2 Example

We consider the same input from Section 3.1.4 consisting of $N=10$ points with
$M=2$. The matrix $\mathbf{P}$ was:

$\mathbf{P}=\begin{bmatrix}3&5&6&2&4&1&7&8&9&10\\\
3&8&1&7&2&6&9&4&10&5\end{bmatrix}^{T}$

The algorithm will iterate over $\mathbf{p}_{2}$ (the last column of
$\mathbf{P}$) in the order: 3, 8, 1, 7, 2, 6, 9, 4, 10, 5. Since initially all
individuals start at rank one, the first rank set $r_{1}$ will contain all the
permutation indices. Each subsequent rank set is initialized with $\emptyset$.
Then, the dominance set and rank sets are updated as shown in Figure 2(b). As
their ranks get updated, the permutation indices of the respective solutions
are removed from the current rank set and inserted into the next one. Note
that the dominance sets and rank sets shown in Figure 2(b) represent the final
values resulting after the assignments at lines 19 and 20, respectively, in
Algorithm 2.

##  4 Benchmark Results

The effectiveness of the proposed algorithms is measured both synthetically
and in a practical optimization setting with the help of the _Pagmo_ framework
(Biscani and Izzo, 2020), a C`++` library for massively parallel optimization
that provides a unified interface to optimization algorithms and problems.

We extended _Pagmo_ with C`++` implementations for the tested algorithms and
ran experiments with the DTLZ benchmark set (Deb et al., 2005) (problems DTLZ1
to DTLZ5) using the NSGA2 algorithm (Deb et al., 2002).

All tests were performed on a workstation containing an AMD Ryzen™ 5950X CPU
with 8Mb of L2 cache and 64Mb of L3 cache, and 64Gb of DDR4 3600Mhz-CL16
memory. The experiments were run on a single thread to minimize effects such
as resource contention or CPU pipeline stalls and to be able to acquire
accurate profiling results. The code was compiled on GNU/Linux using the GNU
C`++` compiler version 11.2.0 with compilation flags: -O3 -mavx2 -mfma.

Each tested algorithm was implemented according to descriptions and pseudocode
made available in their respective publications. A reasonable amount of effort
has been expended for profiling and tuning each implementation. In the case of
MNDS, the code has been ported from the original Java implementation in
jMetal.

Since not all algorithms handle duplicate solutions in the same way, the
procedure for handling duplicates in the population was abstracted away and
made common between implementations. This ensures consistency and
reproducibility among algorithms such that, regardless of the concrete non-
dominated sorting implementation, the same random seed leads to the same
results.

Due to runtime and space constraints, we first compare a larger selection of
algorithms, consisting of RO, RS, MNDS, ENS-BS, ENS-SS, DS and HS on the DTLZ1
test problem, using 5 repetitions of the NGSA2 algorithm. We report the
average elapsed time over the 5 algorithmic runs (Figure 3).

We then focus on a direct comparison between the three fastest algorithms, RO,
RS and MNDS, where we report the average elapsed time over all runs (20
repetitions) and all problems (DTLZ1 to DTLZ5).

The NSGA2 algorithm was parameterized as follows:

  * •

Population size: 500, 1000, 2500, 5000, 7500, 10000

  * •

Maximum generations: 100

  * •

Problem dimension: 25

  * •

Crossover rate: 95%

  * •

Crossover distribution rate $\eta_{c}$: 10

  * •

Mutation rate: 10%

  * •

Mutation distribution rate $\eta_{m}$: 50

The results indicate that bitset-based approaches RS and MNDS are the overall
fastest, followed by RO. In the case of two objectives, the Efficient Non-
dominated Sort family of algorithms (ENS-SS and ENS-BS) also performs very
well, being similarly fast to RS and MNDS and marginally ahead of RO. This is
a remarkable result especially considering their simplicity. However, as the
number of objectives increases ($M=3$, $M=5$) these two algorithms become
slower due to the overhead of direct dominance comparisons within their
respective search strategy. The elapsed time for each algorithm with varying
number of objectives is shown in Figure 3.

Figure 3: Pagmo DTLZ1 – elapsed time (seconds) averaged over five algorithmic
runs

We perform a detailed comparison of the fastest algorithms, RO, RS and MNDS,
on the DTLZ benchmark set consisting of problems DTLZ1 to DTLZ5. Each
algorithm was tested with the same random seeds over 20 repetitions for each
configuration. The number of objectives was varied from 1 to 10 with an
increment of 1 and from 10 to 20 with an increment of 2.

Figure 4 shows that _RankIntersect_ provides a consistent runtime benefit,
finishing ahead of MNDS in all cases. _RankOrdinal_ performs better than MNDS
for small populations but becomes slower as the population size increases
($N>5000$), exhibiting an unexpected performance decrease for $M<10$. This
aspect is not yet fully understood, but we suspect some inefficiency at the
level of the ordinal rank comparison, also considering the fact that the size
of the dominance set is inversely proportional with $M$.

Next, we tested the performance of all algorithms on a synthetic benchmark
where the objective values are uniformly sampled from the unit hypercube. The
results shown in Figure 5 again show that RO is slightly faster than MNDS for
lower population sizes but becomes slower as population size increases
($N>5000$). The same slowdown for $M<10$ is observed for RO on synthetic data.
The other algorithms are noticeably slower than the trio RO, RS, MNDS but
perform similarly well to each other. As the points are random, their
respective ranking strategies are dominated by the runtime cost of Pareto
dominance checks.

Figure 4: RO, RS, MNDS – Pagmo DTLZ benchmark set – elapsed time (seconds)
averaged over twenty runs and five problems (DTLZ1 to DTLZ5 Figure 5: RO, RS,
MNDS – Synthetic benchmark randomly uniform objective values

##  5 Conclusion

In this paper we introduced two simple and performant algorithms for non-
dominated sorting, one of the most runtime-intensive components of Pareto-
based MOEAs.

In contrast to similar approaches BOS or MNDS, Rank Sort defines a dominance
set as the set of solutions dominated by the current solution (and not
dominating the current solution). The main insight leading to the better
performance of the Rank Sort algorithms is that a solution’s rank needs only
be updated when the solution is found to be dominated by another solution of
the same rank. This further reduces the size of the dominance set that needs
to be examined.

The first algorithm, _RankOrdinal_ , uses the comparison of ordinal ranks to
establish dominance. Instead of performing set intersection operations to
compute dominance sets, it simply iterates over the smallest objective-wise
dominance set. _RankOrdinal_ is slower than _RankIntersect_ and MNDS but
requires only $O(N)$ space. It is also slightly slower than ENS-SS and ENS-BS
in the two-objective case, but outperforms them as the number of objectives
increases.

The second algorithm, _RankIntersect_ , uses set intersections to compute
dominance sets. These are efficiently implemented using bit-level parallelism,
which is the main factor in its performance.

The two algorithms have $O(MN\log N)$ best-case and $O(MN^{2})$ worst-case
asymptotic complexity, in line with other non-dominated sorting algorithms. As
already apparent from Algorithms 1 and 2, the proposed methods are simple and
easy to include in other frameworks. They use only basic data structures and
have low constant factors. Nevertheless, C`++` implementations are provided
for all the methods tested in this paper.

Future work will focus on a better understanding of the computational
complexity of the proposed methods and on exploring other algorithmic ideas.
For example, performance can likely be improved by hybridizing the two steps
of the algorithm: stable sort and rank update into a single hybrid step in
order to further minimize the number of operations.

It will also be interesting to profile each algorithm and express its
complexity in terms of low-level hardware events (executed instructions,
branches, branch misses, cache misses). This will serve to provide a unified
framework for the analysis of asymptotic complexity.

We also plan to provide efficient implementations of other state of the art
non-dominated sorting algorithms and to publish all the algorithms as a stand-
alone open-source library.

## References

  * Bao et al. [2017] Chunteng Bao, Lihong Xu, Erik D. Goodman, and Leilei Cao.  A novel non-dominated sorting algorithm for evolutionary multi-objective optimization.  _Journal of Computational Science_ , 23:31–43, 2017.  ISSN 1877-7503.  doi: https://doi.org/10.1016/j.jocs.2017.09.015.  URL https://www.sciencedirect.com/science/article/pii/S1877750317310530.
  * Biscani and Izzo [2020] Francesco Biscani and Dario Izzo.  A parallel global multiobjective framework for optimization: pagmo.  _Journal of Open Source Software_ , 5(53):2338, 2020.  doi: 10.21105/joss.02338.  URL https://doi.org/10.21105/joss.02338.
  * Buzdalov [2018] Maxim Buzdalov.  Generalized offline orthant search: One code for many problems in multiobjective optimization.  In _Proceedings of the Genetic and Evolutionary Computation Conference_ , GECCO ’18, page 593–600, New York, NY, USA, 2018. Association for Computing Machinery.  ISBN 9781450356183.  doi: 10.1145/3205455.3205469.  URL https://doi.org/10.1145/3205455.3205469.
  * Deb et al. [2002] Kalyanmoy Deb, Samir Agrawal, Amrit Pratap, and T. Meyarivan.  A fast and elitist multiobjective genetic algorithm: Nsga-ii.  _IEEE Trans. Evol. Comput_ , 6(2):182–197, 2002\.
  * Deb et al. [2005] Kalyanmoy Deb, Lothar Thiele, Marco Laumanns, and Eckart Zitzler.  Scalable test problems for evolutionary multiobjective optimization.  In _Evolutionary Multiobjective Optimization_ , 2005.
  * Köppen et al. [2005] Mario Köppen, Raul Vicente-Garcia, and Bertram Nickolay.  The pareto-box problem for the modelling of evolutionary multiobjective optimization algorithms.  In _Adaptive and Natural Computing Algorithms_ , 2005.
  * Laumanns et al. [2002] Marco Laumanns, Lothar Thiele, Kalyanmoy Deb, and Eckart Zitzler.  Combining Convergence and Diversity in Evolutionary Multiobjective Optimization.  _Evolutionary Computation_ , 10(3):263–282, 09 2002.  ISSN 1063-6560.  doi: 10.1162/106365602760234108.  URL https://doi.org/10.1162/106365602760234108.
  * McClymont and Keedwell [2012] Kent McClymont and Ed Keedwell.  Deductive Sort and Climbing Sort: New Methods for Non-Dominated Sorting.  _Evolutionary Computation_ , 20(1):1–26, 03 2012\.  ISSN 1063-6560.  doi: 10.1162/EVCO˙a˙00041.  URL https://doi.org/10.1162/EVCO_a_00041.
  * Mishra and Buzdalov [2020] Sumit Mishra and Maxim Buzdalov.  If unsure, shuffle: Deductive sort is $o(mn^{3})$, but $o(mn^{2})$ in expectation over input permutations.  In _Proceedings of the 2020 Genetic and Evolutionary Computation Conference_ , GECCO ’20, page 516–523, New York, NY, USA, 2020. Association for Computing Machinery.  ISBN 9781450371285.  doi: 10.1145/3377930.3390246.  URL https://doi.org/10.1145/3377930.3390246.
  * Moreno et al. [2020] Javier Moreno, Daniel Rodriguez, Antonio J. Nebro, and Jose A. Lozano.  Merge nondominated sorting algorithm for many-objective optimization.  _IEEE Transactions on Cybernetics_ , page 1–11, 2020.  ISSN 2168-2275.  doi: 10.1109/tcyb.2020.2968301.  URL http://dx.doi.org/10.1109/TCYB.2020.2968301.
  * Roy et al. [2016] Proteek Chandan Roy, Md. Monirul Islam, and Kalyanmoy Deb.  Best order sort: A new algorithm to non-dominated sorting for evolutionary multi-objective optimization.  In _Proceedings of the 2016 on Genetic and Evolutionary Computation Conference Companion_ , GECCO ’16 Companion, pages 1113–1120, New York, NY, USA, 2016. ACM.  ISBN 978-1-4503-4323-7.  doi: 10.1145/2908961.2931684.  URL http://doi.acm.org/10.1145/2908961.2931684.
  * Xue et al. [2020] Lingling Xue, Peng Zeng, and Haibin Yu.  Setnds: A set-based non-dominated sorting algorithm for multi-objective optimization problems.  _Applied Sciences_ , 10(19), 2020.  ISSN 2076-3417.  doi: 10.3390/app10196858.  URL https://www.mdpi.com/2076-3417/10/19/6858.
  * Yukish et al. [2004] Michael A. Yukish, W. Simpson, Mark Traband, Soundar R. T. Kumara, and Richard C. Benson.  Algorithms to identify pareto points in multi-dimensional data sets.  Technical report, Mechanical Engineering Dept., The Pennsylvania State University, State College, 2004.
  * Zhang et al. [2015] Xingyi Zhang, Ye Tian, Ran Cheng, and Yaochu Jin.  An efficient approach to nondominated sorting for evolutionary multiobjective optimization.  _IEEE Transactions on Evolutionary Computation_ , 19(2):201–213, 2015.  doi: 10.1109/TEVC.2014.2308305.
  * Zhou et al. [2017] Yuren Zhou, Zefeng Chen, and Jun Zhang.  Ranking vectors by means of the dominance degree matrix.  _IEEE Transactions on Evolutionary Computation_ , 21(1):34–51, 2017.  doi: 10.1109/TEVC.2016.2567648.


# Drude weight in presence of nonlocal potentials

Raffaele Resta  Istituto Officina dei Materiali IOM-CNR, Strada Costiera 11,
314151 Trieste, Italy  Donostia International Physics Center, 20018 San
Sebastián, Spain

(run through LaTeX on at 10:41)

###### Abstract

The nonlocal potential contributes an extra term to the velocity operator; I
show here that such term affects the formal expression of the Drude weight in
a nontrivial way. Notably, the present main result fixes a disturbing
discrepancy in the Dreyer-Coh-Stengel sum rule [Phys. Rev. Lett. 128, 095901
(2022)].

Since the early times of norm-conserving pseudopotentials in the late 1970s it
became clear that in presence of nonlocality the electronic velocity operator
$\hat{\bf v}$ differs from $\hat{\bf p}/m$ by an extra term, owing to the
basic definition $\hat{\bf v}=i[\hat{H},{\bf r}]/\hbar$.

Within band-structure theory one writes the Schrödinger Eq. as

$\hat{H}_{\bf k}|u_{j{\bf k}}\rangle=\epsilon_{j{\bf k}}|u_{j{\bf k}}\rangle,$
(1)

where the eigenstates $|u_{j{\bf k}}\rangle$ are lattice-periodical with band
energies $\epsilon_{j{\bf k}}$, and

$\hat{H}_{\bf k}={\rm e}^{-i{\bf k}\cdot{\bf r}}\hat{H}{\rm e}^{i{\bf
k}\cdot{\bf r}}=\frac{1}{2m}\left(\hat{\bf p}+\hbar{\bf
k}\right)^{2}+\hat{V}_{{\bf k}}.$ (2)

The velocity of a Bloch electron in band $j$ is then

${\bf v}_{j{\bf k}}=\frac{1}{\hbar}\frac{\partial\epsilon_{j{\bf
k}}}{\partial{\bf k}}=\langle u_{j{\bf k}}|\,\left(\frac{\hat{\bf p}+\hbar{\bf
k}}{m}+\frac{1}{\hbar}\frac{\partial\hat{V}_{{\bf k}}}{\partial{\bf
k}}\right)\,|u_{j{\bf k}}\rangle,$ (3)

and the macroscopic current density in dimension $d$ is (for double band
occupancy)

${\bf j}=-\frac{2e}{\hbar}\sum_{j}\int_{\rm BZ}\\!\\!\\!\\!\\![d{\bf
k}]f(\epsilon_{j{\bf k}})\frac{\partial\epsilon_{j{\bf k}}}{\partial{\bf k}},$
(4)

where $[d{\bf k}]=d{\bf k}/(2\pi)^{d}$ and $f(\epsilon)$ is the Fermi
occupation factor.

Little attention has been devoted so far to nonlocal effects on the Drude
weight. The customary expressions—derived in the local case—areAllen06 ;
rap157

$D^{(\rm local)}_{\alpha\beta}=\frac{2\pi e^{2}}{\hbar^{2}}\sum_{j}\int_{\rm
BZ}\\!\\!\\!\\!\\![d{\bf k}]f(\epsilon_{j{\bf
k}})\,\frac{\partial^{2}\epsilon_{j{\bf k}}}{\partial k_{\alpha}\partial
k_{\beta}},$ (5) $D^{(\rm local)}_{\alpha\beta}=-2\pi e^{2}\sum_{j}\int_{\rm
BZ}\\!\\!\\!\\!\\![d{\bf k}]f^{\prime}(\epsilon_{j{\bf
k}})\;v_{jk_{\alpha}}v_{jk_{\beta}};$ (6)

the Fermi-volume and the Fermi-surface expressions are related to each other
by an integration by parts; Eq. (6) also admits a semiclassical
derivation.semicl

To the best of the author’s knowledge, Eqs. (5) and (6) are adopted as they
stand in the existing literature, even when nonlocal potentials are
adopted;Dreyer22 ; Marchese23  in the following, I am going to show that the
above expressions need instead to be modified in the nonlocal case.

In presence of a macroscopic field ${\bf E}$ the vector potential contributes
to the velocity by the term $e{\bf A}/c$, linear in time. The time-derivative
of Eq. (4) is then, to lowest order in ${\bf E}$,

$\displaystyle\partial_{t}j_{\alpha}$ $\displaystyle=$
$\displaystyle-\frac{2e}{\hbar}\sum_{j}\int_{\rm BZ}\\!\\!\\!\\!\\![d{\bf
k}]f(\epsilon_{j{\bf k}})\frac{\partial^{2}\epsilon_{j{\bf k}}}{\partial
k_{\alpha}\partial A_{\beta}}\partial_{t}A_{\beta}$ (7) $\displaystyle=$
$\displaystyle\frac{2ec}{\hbar}\sum_{j}\int_{\rm BZ}\\!\\!\\!\\!\\![d{\bf
k}]f(\epsilon_{j{\bf k}})\frac{\partial^{2}\epsilon_{j{\bf k}}}{\partial
k_{\alpha}\partial A_{\beta}}E_{\beta};$

this clearly shows that the many-electron system undergoes free acceleration.

In the local case one may exploit

$\frac{\partial}{\partial{\bf A}}=\frac{e}{\hbar
c}\frac{\partial}{\partial{\bf k}},$ (8)

in order to obtain

$\displaystyle\partial_{t}j_{\alpha}$ $\displaystyle=$
$\displaystyle-\frac{2e^{2}}{\hbar^{2}}\sum_{j}\int_{\rm
BZ}\\!\\!\\!\\!\\![d{\bf k}]f(\epsilon_{j{\bf
k}})\frac{\partial^{2}\epsilon_{j{\bf k}}}{\partial k_{\alpha}\partial
k_{\beta}}E_{\beta}$ (9) $\displaystyle=$ $\displaystyle\frac{D^{(\rm
local)}_{\alpha\beta}}{\pi}E_{\beta}.$

The Fourier transform of Eq. (9) yields the familiar expression

$j_{\alpha}(\omega)=\sigma^{(\rm Drude)}_{\alpha\beta}(\omega)E_{\beta},$ (10)
$\sigma^{(\rm Drude)}_{\alpha\beta}(\omega)=D^{(\rm
local)}_{\alpha\beta}\left[\delta(\omega)+\frac{i}{\pi\omega}\right];$ (11)

a derivation can be found e.g. in the Appendix of Ref. rap160, .

In the nonlocal case Eq. (8) no longer holds; it is then expedient to define

$\displaystyle\tilde{\bf v}_{j{\bf k}}$ $\displaystyle=$
$\displaystyle\frac{1}{m}\langle u_{j{\bf k}}|\,\left(\hat{\bf p}+\hbar{\bf
k}\right)\,|u_{j{\bf k}}\rangle=\frac{c}{e}\langle u_{j{\bf
k}}|\,\frac{\partial\hat{H}_{{\bf k}}}{\partial{\bf A}}\,|u_{j{\bf k}}\rangle$
(12) $\displaystyle=$ $\displaystyle\frac{c}{e}\frac{\partial\epsilon_{j{\bf
k}}}{\partial{\bf A}}.$

Integrating by parts Eq. (7) one gets

$\displaystyle\partial_{t}j_{\alpha}$ $\displaystyle=$
$\displaystyle-\frac{2ec}{\hbar}\sum_{j}\int_{\rm BZ}\\!\\!\\!\\!\\![d{\bf
k}]f^{\prime}(\epsilon_{j{\bf k}})\frac{\partial\epsilon_{j{\bf k}}}{\partial
k_{\alpha}}\frac{\partial\epsilon_{j{\bf k}}}{\partial A_{\beta}}E_{\beta};$
(13) $\displaystyle=$ $\displaystyle\frac{D_{\alpha\beta}}{\pi}E_{\beta}.,$

where the Drude weight is now expressed as

$D_{\alpha\beta}=-2\pi e^{2}\sum_{j}\int_{\rm BZ}\\!\\!\\!\\!\\![d{\bf
k}]f^{\prime}(\epsilon_{j{\bf k}})\;v_{jk_{\alpha}}\tilde{v}_{jk_{\beta}},$
(14)

to be compared to Eq. (6).

The Drude weight measures the inverse inertia of the many-electron system. In
the special case of a flat potential its value is

$D^{(\rm free)}_{\alpha\beta}=\frac{\pi e^{2}n}{m}\delta_{\alpha\beta},$ (15)

where $n$ is the density of conducting electrons; in the general case one then
recasts $D_{\alpha\beta}$ as

$D_{\alpha\beta}=\frac{\pi e^{2}n^{*}_{\alpha\beta}}{m},$ (16)

where $n^{*}_{\alpha\beta}$ has the meaning of the effective electron density
contributing to the adiabatic current.Scalapino93 ; rap157

Remarkably, the improved expression of Eq. (14) fixes a disturbing discrepancy
in the sum rule obeyed by the Born effective charges in metals, recently found
by Dreyer, Coh, and Stengel.Dreyer22  Suppose that the sublattice $s$ is
displaced with velocity ${\bf v}_{s}$: the dimensionless Born tensor is by
definition (in either insulators or metals)

$Z^{*}_{s,\alpha\beta}=Z_{s}\delta_{\alpha\beta}+\frac{V_{\rm
cell}}{e}\frac{\partial j_{\alpha}}{\partial v_{s\beta}},$ (17)

where $eZ_{s}$ is the bare core charge and $V_{\rm cell}$ is the unit-cell
volume. When the whole lattice is rigidly translated with velocity ${\bf v}$,
the linearly-induced current density (due to electrons and nuclei) is

$j^{(\rm total)}_{\alpha}=\frac{e}{V_{\rm
cell}}\left(\sum_{s}Z^{*}_{s,\alpha\beta}\right)v_{\beta}.$ (18)

In the insulating case no current may flow for a rigid translation of the
whole solid, whence the basic relationship $\sum_{s}Z^{*}_{s,\alpha\beta}=0$,
called the acoustic sum rule. In the metallic case some current instead flows;
in the reference frame of the nuclei it is carried solely by the electrons,
all moving with velocity $-{\bf v}$, and whose effective density is
$n^{*}_{\alpha\beta}$. Therefore from Eq. (16)

$j^{(\rm
total)}_{\alpha}=e\,n^{*}_{\alpha\beta}v_{\beta}=\frac{mD_{\alpha\beta}}{\pi
e}v_{\beta}.;$ (19)

by comparing Eq. (19) to Eq. (18) one gets the Dreyer-Coh-Stengel sum rule:

$\frac{1}{V_{\rm cell}}\sum_{s}Z^{*}_{s,\alpha\beta}=\frac{m}{\pi
e^{2}}D_{\alpha\beta}.$ (20)

As said above, the sum rule is exact—even in the nonlocal case—provided Eq.
(14) is adopted for the Drude weight.

In insulators the Born effective charges manifest themselves in the phonon
spectra: they generate a depolarization field ${\bf E}$ which contributes to
the restoring force at the zone center. No such contribution may exist in
metals in the adiabatic limit, since macroscopic fields are screened therein
(Faraday-cage effect). A few calculations of Born charges in metals are
available;Dreyer22 ; Marchese23  their possible experimental
manifestations—e.g. in the case of a binary metal—are presently under
discussion.Hickox23

I thank Francesco Mauri for bringing this problem to my attention. Work
supported by the Office of Naval Research (USA) Grant No. N00014-20-1-2847.

## References

  * (1) P. B. Allen, in: Conceptual foundations of materials: A standard model for ground- and excited-state properties, S.G. Louie and M.L. Cohen, eds. (Elsevier, 2006), p. 139.
  * (2) R. Resta, J. Phys. Condens. Matter 30, 414001 (2018).
  * (3) N. W. Ashcroft and N. D. Mermin, Solid State Physics (Saunders, Philadelphia, 1976), Ch. 13.
  * (4) C. E. Dreyer, S. Coh, and M. Stengel, Phys. Rev. Lett. 128, 095901 (2022).
  * (5) G. Marchese et al., Nat. Phys. (2023). https://doi.org/10.1038/s41567-023-02203-3.
  * (6) G. Bellomia and R. Resta, Phys. Rev. B 102, 205123 (2020).
  * (7) D. J. Scalapino, S. R. White, and S. C. Zhang, Phys. Rev. B 47, 7995 (1993).
  * (8) D. Hickox-Young , D. Puggioni, and J. M. Rondinelli, Phys. Rev. Materials 7, 010301 (2023).