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.

```
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> 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> 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:

```
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> 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:

```
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> 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> 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)

```
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
```