# Fast non dominated sorting

**URL:** <https://discourse.julialang.org/t/fast-non-dominated-sorting/24164>\
**Category:** Optimization (Mathematical)\
**Tags:** sort\
**Created:** [May 13, 2019, 7:16am UTC](https://discourse.julialang.org/t/fast-non-dominated-sorting/24164 "2019-05-13T07:16:03Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![Igor\_Douven](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/igor_douven/32/8433_2.png) [@Igor\_Douven](https://discourse.julialang.org/u/Igor_Douven)\
**Post date:** [May 13, 2019, 7:16am UTC](https://discourse.julialang.org/t/fast-non-dominated-sorting/24164/1 "2019-05-13T07:16:04Z")

</div>

I need a fast implementation of the so-called non dominated sorting algorithm. The best I could come up with so far is the following (the sorting is to be done on the basis of three criteria):

```julia
dominates(x, y) = all(x .<= y) && any(x .< y)

function nonDominatedSorting(arr::Array{Float64,2})
    fronts::Array{Array,1} = Array[]
    ind::Array{Int64,1} = collect(1:size(arr, 1))
    while !isempty(arr)
        s = size(arr, 1)
        red = dropdims(sum([dominates(arr[i,:], arr[j,:]) for i in 1:s, j in 1:s ], dims=1) .== 0, dims=1)
        a = 1:s
        sel::Array{Int64,1} = a[red]
        push!(fronts, ind[sel])
        da::Array{Int64,1} = deleteat!(collect(1:s), sel)
        ind = deleteat!(ind, sel)
        arr = arr[da, :]
    end
    return fronts
end

```

The timing is not too bad:

```julia
using BenchmarkTools

scores = rand(Float64, (50, 3));

@btime nonDominatedSorting($scores)
# 10.039 ms (31998 allocations: 21.81 MiB)

```

But I need this for simulations, so at the moment I use the faster

```julia
using RCall

function nds(arr::Array{Float64,2})
    @rput arr
    R"""
    library(nsga2R)
    ranking <- fastNonDominatedSorting(arr)
    """
    return @rget ranking
end

```

which gives

```julia
@btime nds($scores)
# 2.354 ms (227 allocations: 9.97 KiB)

```

That’s fine for what I want, but I’d still be happier to have a pure Julia implementation. So I’d be grateful for any suggestions for making the first version faster.

---

<div class="post-metadata">

**Author:** ![blegat](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/blegat/32/217090_2.png) [@blegat](https://discourse.julialang.org/u/blegat)\
**Post date:** [May 13, 2019, 7:36am UTC](https://discourse.julialang.org/t/fast-non-dominated-sorting/24164/2 "2019-05-13T07:36:59Z")

</div>

You might have more answers in [https://discourse.julialang.org/c/user/perf](https://discourse.julialang.org/c/user/perf) (reference this post, if you post there).  
Here is a first improvement:

```julia
julia> x = rand(Int, 100000);

julia> y = rand(Int, 100000);

julia> dominates(x, y) = all(x .<= y) && any(x .< y)

julia> dominates2(x, y) = all(i -> x[i] <= y[i], eachindex(x)) && any(i -> x[i] < y[i], eachindex(x))
dominates2 (generic function with 1 method)

julia> _time(f, x, y) = @time f(x, y) # When `@time` is called from global scope, it has a fixed cost in term of speed and allocation so we make this dummy _time function
_time (generic function with 1 method)

julia> _time(dominates, x, y)
  0.000118 seconds (3 allocations: 16.594 KiB)
false

julia> _time(dominates, x, y)
  0.000186 seconds (3 allocations: 16.594 KiB)
false

julia> _time(dominates2, x, y)
  0.000001 seconds
false

julia> _time(dominates2, x, y)
  0.000001 seconds
false

```

The reason for this improvement is that `dominates2` does create the intermediate array `x .<= y` so it does not allocate as much and it does not evaluate `x <= y` all the way and stops when it is false for one `i`.

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [May 13, 2019, 7:54am UTC](https://discourse.julialang.org/t/fast-non-dominated-sorting/24164/3 "2019-05-13T07:54:46Z")

</div>

> [@Igor\_Douven](#):
>
> arr[i,:]

You could try storing all data transposed so you can access `arr[:,i]` instead, as this is how it is stored in memory.

Also, this line  
`sum([dominates(arr[i,:], arr[j,:]) for i in 1:s, j in 1:s ], dims=1)`  
creates a matrix just to sum over it. Writing out this loop would probably save you both a lot of allocations and time.

---

<div class="post-metadata">

**Author:** ![Igor\_Douven](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/igor_douven/32/8433_2.png) [@Igor\_Douven](https://discourse.julialang.org/u/Igor_Douven)\
**Post date:** [May 13, 2019, 8:36am UTC](https://discourse.julialang.org/t/fast-non-dominated-sorting/24164/4 "2019-05-13T08:36:36Z")

</div>

Thank you, both @blegat and @baggepinnen – these suggestions helped a lot. I now have

```julia
dominates(x, y) = all(i -> x[i] <= y[i], eachindex(x)) && any(i -> x[i] < y[i], eachindex(x))

function nonDominatedSorting(arr0::Array{Float64,2})
    arr = arr0'
    fronts::Array{Array,1} = Array[]
    ind::Array{Int64,1} = collect(1:size(arr, 2))
    while !isempty(arr)
        s = size(arr, 2)
        red = Array{Bool,2}(undef, s, s)
        for i in 1:s, j in 1:s
            red[i, j] = dominates(arr[:, i], arr[:, j])
        end
        sel::Array{Int64,1} = collect(1:s)[dropdims(sum(red, dims=1) .== 0, dims=1)]
        push!(fronts, ind[sel])
        da::Array{Int64,1} = deleteat!(collect(1:s), sel)
        ind = deleteat!(ind, sel)
        arr = arr[:,da]
    end
    return fronts
end

```

which is a lot faster than the R function I used:

```julia
scores = rand(Float64, (50, 3))

@btime nonDominatedSorting($scores)
# 394.493 μs (8801 allocations: 987.38 KiB)

```

---

<div class="post-metadata">

**Author:** ![GunnarFarneback](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gunnarfarneback/32/1827_2.png) [@GunnarFarneback](https://discourse.julialang.org/u/GunnarFarneback)\
**Post date:** [May 13, 2019, 8:52am UTC](https://discourse.julialang.org/t/fast-non-dominated-sorting/24164/5 "2019-05-13T08:52:43Z")

</div>

This should give you another order of magnitude or so:

```julia
function dominates2(x, y)
    strict_inequality_found = false
    for i in eachindex(x)
        y[i] < x[i] && return false
        strict_inequality_found |= x[i] < y[i]
    end
    return strict_inequality_found
end

function nds2(arr)
    fronts = Vector{Int64}[]
    ind = collect(1:size(arr, 1))
    a = [arr[i, :] for i = 1:size(arr, 1)]
    while !isempty(a)
        s = length(a)
        red = [all(x -> !dominates2(x, y), a) for y in a]
        push!(fronts, ind[red])
        da = deleteat!(collect(1:s), red)
        ind = deleteat!(ind, red)
        a = a[da]
    end
    return fronts
end

```

Edit: You may want to check correctness, especially in the presence of equal values. Having a testset of tricky cases is good practice anyway.

---

<div class="post-metadata">

**Author:** ![Igor\_Douven](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/igor_douven/32/8433_2.png) [@Igor\_Douven](https://discourse.julialang.org/u/Igor_Douven)\
**Post date:** [May 13, 2019, 9:07am UTC](https://discourse.julialang.org/t/fast-non-dominated-sorting/24164/6 "2019-05-13T09:07:04Z")

</div>

Many thanks @GunnarFarneback – this is correct, as far as I have been able to check, and is indeed much faster:

```julia
@btime nds2($scores)
# 37.822 μs (266 allocations: 14.38 KiB)

```

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [May 13, 2019, 9:26am UTC](https://discourse.julialang.org/t/fast-non-dominated-sorting/24164/7 "2019-05-13T09:26:44Z")

</div>

> [@Igor\_Douven](#):
>
> arr = arr0’

Note that trsnaposing is lazy, this does not actually change the memory layout, so your `arr[:,i]` is still ineffective. If you do `arr = copy(arr0')` the data is copied and put into the correct memory layout.

---

<div class="post-metadata">

**Author:** ![Igor\_Douven](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/igor_douven/32/8433_2.png) [@Igor\_Douven](https://discourse.julialang.org/u/Igor_Douven)\
**Post date:** [May 13, 2019, 9:31am UTC](https://discourse.julialang.org/t/fast-non-dominated-sorting/24164/8 "2019-05-13T09:31:48Z")

</div>

Thanks, I didn’t know. This made the function about twice as fast, but it’s still slower than @GunnarFarneback’s `nds2`, so for now, I’ll go with that.

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [May 13, 2019, 11:10am UTC](https://discourse.julialang.org/t/fast-non-dominated-sorting/24164/9 "2019-05-13T11:10:13Z")

</div>

I’m not sure what this thing is doing, but you can get more speedup (\>2x over the one from @GunnarFarneback) by doing this:

```julia
function nds3(arr)
   fronts = Vector{Int64}[]
   ind = collect(axes(arr, 1))
   a = collect(eachrow(arr))
   while !isempty(a)
       red = [all(x -> !dominates2(x, y), a) for y in a]
       push!(fronts, ind[red])
       deleteat!(ind, red)
       deleteat!(a, red)
   end
   return fronts
end

```

---

<div class="post-metadata">

**Author:** ![Igor\_Douven](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/igor_douven/32/8433_2.png) [@Igor\_Douven](https://discourse.julialang.org/u/Igor_Douven)\
**Post date:** [May 13, 2019, 11:47am UTC](https://discourse.julialang.org/t/fast-non-dominated-sorting/24164/10 "2019-05-13T11:47:33Z")

</div>

That’s fantastic, many thanks!

What a nice community this is: within a couple of hours, the function got about 500 times faster thanks to the joint help of people on here!

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [May 13, 2019, 11:53am UTC](https://discourse.julialang.org/t/fast-non-dominated-sorting/24164/11 "2019-05-13T11:53:04Z")

</div>

BTW, if the `scores` matrix has a fixed, known and small size in the second dimension, you could restructure it as a vector of `SVector`s, to possibly get even further speedups.

---

<div class="post-metadata">

**Author:** ![Igor\_Douven](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/igor_douven/32/8433_2.png) [@Igor\_Douven](https://discourse.julialang.org/u/Igor_Douven)\
**Post date:** [May 13, 2019, 1:46pm UTC](https://discourse.julialang.org/t/fast-non-dominated-sorting/24164/12 "2019-05-13T13:46:26Z")

</div>

Thanks, the only way I can get this to work is by using

```julia
scores = @SArray rand(50, 3)

```

but then the timing actually gets worse.

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [May 13, 2019, 3:19pm UTC](https://discourse.julialang.org/t/fast-non-dominated-sorting/24164/13 "2019-05-13T15:19:28Z")

</div>

No I was thinking more along the lines of a `Vector` of `SVector`s, not a large `SMatrix`. StaticArrays are not particularly performant when they become large (\>100 elements, roughly). Also, you may have to change the code, I think.

---

<div class="post-metadata">

**Author:** ![GunnarFarneback](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gunnarfarneback/32/1827_2.png) [@GunnarFarneback](https://discourse.julialang.org/u/GunnarFarneback)\
**Post date:** [May 13, 2019, 3:31pm UTC](https://discourse.julialang.org/t/fast-non-dominated-sorting/24164/14 "2019-05-13T15:31:04Z")

</div>

Maybe something similar to this:

```julia
function nds4(arr)
   fronts = Vector{Int64}[]
   ind = collect(axes(arr, 1))
   a = SVector{3}.(eachrow(arr))
   while !isempty(a)
       red = [all(x -> !dominates2(x, y), a) for y in a]
       push!(fronts, ind[red])
       deleteat!(ind, red)
       deleteat!(a, red)
   end
   return fronts
end

```

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [May 13, 2019, 3:43pm UTC](https://discourse.julialang.org/t/fast-non-dominated-sorting/24164/15 "2019-05-13T15:43:01Z")

</div>

It seems we just broke 1000-fold speedup limit 😃

---

<div class="post-metadata">

**Author:** ![Igor\_Douven](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/igor_douven/32/8433_2.png) [@Igor\_Douven](https://discourse.julialang.org/u/Igor_Douven)\
**Post date:** [May 13, 2019, 5:06pm UTC](https://discourse.julialang.org/t/fast-non-dominated-sorting/24164/16 "2019-05-13T17:06:54Z")

</div>

Amazing! This gives

```julia
@btime nds4(scores)
# 5.510 μs (100 allocations: 7.30 KiB)

```

Many thanks to all!

---

<div class="post-metadata">

**Author:** ![braamvandyk](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/braamvandyk/32/5086_2.png) [@braamvandyk](https://discourse.julialang.org/u/braamvandyk)\
**Post date:** [May 13, 2019, 6:34pm UTC](https://discourse.julialang.org/t/fast-non-dominated-sorting/24164/17 "2019-05-13T18:34:52Z")

</div>

Julia. Come for the speed. Stay for the community.

---

<div class="post-metadata">

**Author:** ![yakir12](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yakir12/32/297_2.png) [@yakir12](https://discourse.julialang.org/u/yakir12)\
**Post date:** [May 13, 2019, 7:05pm UTC](https://discourse.julialang.org/t/fast-non-dominated-sorting/24164/18 "2019-05-13T19:05:16Z")

</div>

Is there a sorting package or some pace to post this version to…?

---

<div class="post-metadata">

**Author:** ![ffevotte](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ffevotte/32/6587_2.png) [@ffevotte](https://discourse.julialang.org/u/ffevotte)\
**Post date:** [May 14, 2019, 1:47pm UTC](https://discourse.julialang.org/t/fast-non-dominated-sorting/24164/19 "2019-05-14T13:47:05Z")

</div>

It looks like I’m arriving late here, but here are a few more thoughts:

I tend to think that dynamic operations like `deleteat!` and `push!` are costly in general. One might want to avoid them by performing the appropriate swaps in order to move the elements to the beginning of the array. The algorithm can then continue working on an `arr` array of unchanged length, but looking only at the end of it. It is also possible to examine the dominance on elements one by one, instead of computing the full `red` array storing whether each element is dominated.

Here is an implementation trying to take advantage of these observations.

```julia
function nds5(arr)
    fronts = Vector{Int64}[]
    ind = collect(axes(arr, 1))
    a = SVector{3}.(eachrow(arr))
    N = length(a)

    # All elements in arr[1:ii-1] are dominated
    ii = 1
    while ii<=N
        start = ii
        @inbounds for i in ii:N
            for j in start:N
                dominates2(a[j], a[i]) && @goto l1
            end

            # Swap values at indices `ii` and `i` in `ind` and `a`
            ind[ii], ind[i] = ind[i], ind[ii]
            a[ii], a[i] = a[i], a[ii]
            ii += 1

            @label l1
        end
        push!(fronts, ind[start:ii-1])
    end
    return fronts
end

```

Note that it is not strictly equivalent to the original functions, since the sorting is not stable : the same indices are output in a different order (whereas the original function would always have indices sorted in ascending order).

```julia
julia> nds4(scores)
9-element Array{Array{Int64,1},1}:
 [6, 42, 62, 76, 79, 91]                                                          
 [1, 14, 20, 23, 30, 37, 40, 50, 69, 75, 88]                                      
 [4, 10, 12, 18, 19, 22, 36, 48, 52, 61, 63, 64, 68, 83, 84, 87, 97, 98]          
 [2, 3, 11, 13, 21, 25, 26, 29, 38, 41 … 56, 57, 59, 60, 66, 71, 77, 80, 90, 96]
 [15, 16, 17, 35, 39, 43, 44, 45, 67, 70, 81, 82, 89, 94, 100]                    
 [5, 9, 24, 31, 32, 46, 49, 58, 65, 78, 85, 86]                                   
 [8, 27, 33, 34, 47, 53, 72, 73, 92]                                              
 [28, 74, 93, 95]                                                                 
 [7, 99]                                                                          

```

```julia
julia> nds5(scores)
9-element Array{Array{Int64,1},1}:
 [6, 42, 62, 76, 79, 91]                                                          
 [14, 20, 23, 30, 37, 40, 50, 69, 75, 88, 1]                                      
 [18, 19, 22, 10, 36, 12, 48, 52, 61, 63, 64, 68, 4, 83, 84, 87, 97, 98]          
 [11, 38, 41, 2, 13, 51, 25, 54, 55, 56 … 26, 3, 66, 29, 71, 21, 77, 80, 90, 96]
 [67, 39, 70, 15, 43, 44, 45, 81, 82, 16, 89, 17, 94, 35, 100]                    
 [9, 78, 5, 65, 31, 32, 85, 86, 24, 46, 49, 58]                                   
 [8, 33, 47, 27, 92, 53, 34, 72, 73]                                              
 [95, 28, 74, 93]                                                                 
 [99, 7]                                                                          

```

In the implementation above, it is interesting to observe that the `ind` array is shuffled in such a way that it contains all indices of the output partition _in the correct order_. All that is needed to reconstruct the output partition of `1:n` is the list of part boundaries to correctly split the `ind` array.

The following implementation relies on a custom `Partition` type, which allows manipulating partitions of `1:n` in the form of a 1D array of indices (which is a reordering of `1:n`) and an array of boundaries. This avoids manipulating arrays of arrays:

```julia
struct Partition
    ind :: Vector{Int64}
    bounds :: Vector{Int64}
end

Partition(n::Int64) = Partition(collect(1:n), [1])

function push_bound!(p::Partition, i::Int64)
    push!(p.bounds, i)
    p
end

Base.getindex(p::Partition, i::Int64) = if i==lastindex(p.bounds)
    view(p.ind, p.bounds[i]:lastindex(p.ind))
else
    view(p.ind, p.bounds[i]:p.bounds[i+1]-1)
end

function Base.show(io::IO, p::Partition)
    n = length(p.bounds)
    if p.bounds[n] == length(p.ind)+1
        n -= 1
    end
    println("$n-part Partition:")
    for i in 1:n
        println(" ", p[i])
    end
end

function nds6(arr)
    a = SVector{3}.(eachrow(arr))
    N = length(a)
    p = Partition(N)

    ii = 1
    while ii<=N
        start = ii
        @inbounds for i in ii:N
            for j in start:N
                dominates2(a[j], a[i]) && @goto l1
            end

            # Swap values at indices `ii` and `i` in `ind` and `a`
            p.ind[ii], p.ind[i] = p.ind[i], p.ind[ii]
            a[ii], a[i] = a[i], a[ii]
            ii += 1

            @label l1
        end
        push_bound!(p, ii)
    end
    return p
end

```

which yields the same results as `nds5`, except stored in another fashion:

```julia
julia> nds6(scores)
9-part Partition:
 [6, 42, 62, 76, 79, 91]
 [14, 20, 23, 30, 37, 40, 50, 69, 75, 88, 1]
 [18, 19, 22, 10, 36, 12, 48, 52, 61, 63, 64, 68, 4, 83, 84, 87, 97, 98]
 [11, 38, 41, 2, 13, 51, 25, 54, 55, 56, 57, 59, 60, 26, 3, 66, 29, 71, 21, 77, 80, 90, 96]
 [67, 39, 70, 15, 43, 44, 45, 81, 82, 16, 89, 17, 94, 35, 100]
 [9, 78, 5, 65, 31, 32, 85, 86, 24, 46, 49, 58]
 [8, 33, 47, 27, 92, 53, 34, 72, 73]
 [95, 28, 74, 93]
 [99, 7]

```

  

Now some benchmarks:

```julia
using BenchmarkTools

function bench(N)
    scores = rand(Float64, (N, 3));
    @assert nds1(scores) == nds2(scores) == nds3(scores) == nds4(scores)
    print("nds1"); @btime nds1($scores)
    print("nds2"); @btime nds2($scores)
    print("nds3"); @btime nds3($scores)
    print("nds4"); @btime nds4($scores)
    print("nds5"); @btime nds5($scores)
    print("nds6"); @btime nds6($scores)

    nothing
end

```

```julia
julia> bench(50)
nds1 12.494 ms (40026 allocations: 27.14 MiB)
nds2 42.062 μs (265 allocations: 14.36 KiB)
nds3 13.697 μs (240 allocations: 7.97 KiB)
nds4 5.633 μs (94 allocations: 7.00 KiB)
nds5 4.675 μs (64 allocations: 5.77 KiB)
nds6 3.211 μs (59 allocations: 4.94 KiB)

julia> bench(50)
nds1 14.268 ms (46021 allocations: 31.52 MiB)
nds2 48.169 μs (295 allocations: 15.48 KiB)
nds3 17.536 μs (266 allocations: 8.58 KiB)
nds4 6.547 μs (100 allocations: 7.30 KiB)
nds5 5.779 μs (65 allocations: 5.86 KiB)
nds6 3.303 μs (59 allocations: 4.94 KiB)

```

Most of the remaining allocations in `nds6` come from the Array-to-SArray conversion, which takes negligible time and is therefore probably not worth optimizing:

```julia
julia> f(x) = SVector{3}.(eachrow(x))
f (generic function with 1 method)

julia> @btime f($scores)
  693.521 ns (54 allocations: 4.20 KiB)

```

  

* * *

> **Full code (for anyone wanting to reproduce the results)**
>
> ```julia
> dominates(x, y) = all(x .<= y) && any(x .< y)
> 
> function dominates2(x, y)
> strict_inequality_found = false
> for i in eachindex(x)
> y[i] < x[i] && return false
> strict_inequality_found |= x[i] < y[i]
> end
> return strict_inequality_found
> end
> 
> function nds1(arr::Array{Float64,2})
> fronts::Array{Array,1} = Array[]
> ind::Array{Int64,1} = collect(1:size(arr, 1))
> while !isempty(arr)
> s = size(arr, 1)
> red = dropdims(sum([dominates(arr[i,:], arr[j,:]) for i in 1:s, j in 1:s ], dims=1) .== 0, dims=1)
> a = 1:s
> sel::Array{Int64,1} = a[red]
> push!(fronts, ind[sel])
> da::Array{Int64,1} = deleteat!(collect(1:s), sel)
> ind = deleteat!(ind, sel)
> arr = arr[da, :]
> end
> return fronts
> end
> 
> function nds2(arr)
> fronts = Vector{Int64}[]
> ind = collect(1:size(arr, 1))
> a = [arr[i, :] for i = 1:size(arr, 1)]
> while !isempty(a)
> s = length(a)
> red = [all(x -> !dominates2(x, y), a) for y in a]
> push!(fronts, ind[red])
> da = deleteat!(collect(1:s), red)
> ind = deleteat!(ind, red)
> a = a[da]
> end
> return fronts
> end
> 
> function nds3(arr)
> fronts = Vector{Int64}[]
> ind = collect(axes(arr, 1))
> a = collect(eachrow(arr))
> while !isempty(a)
> red = [all(x -> !dominates2(x, y), a) for y in a]
> push!(fronts, ind[red])
> deleteat!(ind, red)
> deleteat!(a, red)
> end
> return fronts
> end
> 
> using StaticArrays
> 
> function nds4(arr)
> fronts = Vector{Int64}[]
> ind = collect(axes(arr, 1))
> a = SVector{3}.(eachrow(arr))
> while !isempty(a)
> red = [all(x -> !dominates2(x, y), a) for y in a]
> push!(fronts, ind[red])
> deleteat!(ind, red)
> deleteat!(a, red)
> end
> return fronts
> end
> 
> function nds5(arr)
> fronts = Vector{Int64}[]
> ind = collect(axes(arr, 1))
> a = SVector{3}.(eachrow(arr))
> N = length(a)
> 
> ii = 1
> while ii<=N
> start = ii
> @inbounds for i in ii:N
> for j in start:N
> dominates2(a[j], a[i]) && @goto l1
> end
> 
> # Swap values at indices `ii` and `i` in `ind` and `a`
> ind[ii], ind[i] = ind[i], ind[ii]
> a[ii], a[i] = a[i], a[ii]
> ii += 1
> 
> @label l1
> end
> push!(fronts, ind[start:ii-1])
> end
> return fronts
> end
> 
> struct Partition
> ind :: Vector{Int64}
> bounds :: Vector{Int64}
> end
> 
> Partition(n::Int64) = Partition(collect(1:n), [1])
> 
> function push_bound!(p::Partition, i::Int64)
> push!(p.bounds, i)
> p
> end
> 
> Base.getindex(p::Partition, i::Int64) = if i==lastindex(p.bounds)
> view(p.ind, p.bounds[i]:lastindex(p.ind))
> else
> view(p.ind, p.bounds[i]:p.bounds[i+1]-1)
> end
> 
> function Base.show(io::IO, p::Partition)
> n = length(p.bounds)
> if p.bounds[n] == length(p.ind)+1
> n -= 1
> end
> println("$n-part Partition:")
> for i in 1:n
> println(" ", p[i])
> end
> end
> 
> function nds6(arr)
> a = SVector{3}.(eachrow(arr))
> N = length(a)
> p = Partition(N)
> 
> ii = 1
> while ii<=N
> start = ii
> @inbounds for i in ii:N
> for j in start:N
> dominates2(a[j], a[i]) && @goto l1
> end
> 
> # Swap values at indices `ii` and `i` in `ind` and `a`
> p.ind[ii], p.ind[i] = p.ind[i], p.ind[ii]
> a[ii], a[i] = a[i], a[ii]
> ii += 1
> 
> @label l1
> end
> push_bound!(p, ii)
> end
> return p
> end
> 
> using BenchmarkTools
> 
> function bench(N)
> scores = rand(Float64, (N, 3));
> @assert nds1(scores) == nds2(scores) == nds3(scores) == nds4(scores)
> print("nds1"); @btime nds1($scores)
> print("nds2"); @btime nds2($scores)
> print("nds3"); @btime nds3($scores)
> print("nds4"); @btime nds4($scores)
> print("nds5"); @btime nds5($scores)
> print("nds6"); @btime nds6($scores)
> 
> nothing
> end
> 
> ```

---

<div class="post-metadata">

**Author:** ![Igor\_Douven](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/igor_douven/32/8433_2.png) [@Igor\_Douven](https://discourse.julialang.org/u/Igor_Douven)\
**Post date:** [May 14, 2019, 3:10pm UTC](https://discourse.julialang.org/t/fast-non-dominated-sorting/24164/20 "2019-05-14T15:10:02Z")

</div>

The ordering of the indices is unimportant: the function calculates the so-called Pareto fronts, and within any front all solutions are equivalent. In other words, `nds5` would have done. I can’t now quite follow why, despite the extra work, `nds6` is still faster. In any case, fantastic to see that a further speedup was still possible.

In response to @yakir12, the non-dominated sorting algorithm is a key component in Deb’s NSGA-II algorithm, which is one of the most popular algorithms for agent-based optimization. At the moment, there are three Julia packages for this, but two seem to have been abandoned and one is, at least for my purposes, not flexible enough. For me, it’s much more useful to have the components of the NSGA-II algorithm available (another key component is the crowding sorting distance) and then to combine them in whichever way best suits my needs for the case at hand. I could imagine a package with fast versions of the necessary functions (like the one developed in this thread for non dominated sorting) to be useful for some people.

[Next page](https://discourse.julialang.org/t/fast-non-dominated-sorting/24164.md?page=2)