Speed up permsort by utilizing stability of the default sorting algorithm

[petvana]

This is a draft on how to speed up permsort. The idea is that it is not necessary to maintain the stability by Perm if a stable sorting algorithm is used. Thus PermUnstable is introduces with a simplify and faster lt(p::PermUnstable, a::Integer, b::Integer) implementation.

Notice the PR builds upon #47383 that has not been merged yet.

@LilithHafner Do you like this idea (as the original author of the changes)?

julia> using BenchmarkTools

#47383
julia> @btime sortperm(x) setup=(x=rand(1000));
  29.741 μs (5 allocations: 15.92 KiB)

# This PR
julia> @btime sortperm(x) setup=(x=rand(1000));
  23.349 μs (5 allocations: 15.92 KiB)

[StefanKarpinski]

This is a great improvement but I think the naming is confusing: I would have guessed that PermUnstable would be used together with an unstable sorting algorithm, not the other way around.

[petvana]

PermUnstable is really a bad name. I can imagine smth like FastPerm or PermFast, inspired by fast math. Currently. the CI fails and I would blame the reverse ordering in IEEEFloatOptimization:

julia/base/sort.jl

Line 625 in 812c917

scratch = _sort!(v, a.next, PermT(Reverse, ip), (;kw..., lo, hi=j))

[LilithHafner]

Thanks! sortperm has always been a weak link in Julia's sorting, performance-wise (#939).

sortperm and sortperm! allow the user to specify an algorithm and the docstring states "The permutation is guaranteed to be stable even if the sorting algorithm is unstable, meaning that indices of equal elements appear in ascending order." While Base only defines stable algorithms, some packages (e.g. SortingAlgorithms.jl) export unstable algorithms. We could

1. Only use the non-stabilizing Perm type for algorithms we recognize, optionally define a trait to let external algorithms opt into declaring stability, and support fastpaths in e.g. IEEESortOptimization for both Perm and PermUstable.
2. Use PkgEval to determine that this PR is okay as a "minor change". This is dicey because stability is subtle and may break things without being detected by tests.

I think it is possible to get substantially larger speedups by computing sortperm with a ZipVector-like data structure, but so long as PermUnstable is not exported, this change shouldn't interfere with that possibility. Indeed it synergizes well for cases when a ZipVector is not appropriate (small arrays, perhaps?).

[petvana]

@nanosoldier runbenchmarks("sort")

[nanosoldier]

Your job failed.

[LilithHafner]

What about sizes between 10 and 40 when we currently perform IEEEFloatOptimizations and then automatically dispatch to InsertionSort? I imagine the allocations could be expensive in that regime.

w.r.t nanosoldier, we know there is a problem, but we don't know what it is (#47795)

[petvana]

What about sizes between 10 and 40

Yes, there is a slight regression and we should address these allocations. (Btw, I would expect only a single allocation on master for sortperm and none for sortperm! if in-place InsertionSort is used.)

julia> @btime sortperm(x) setup=(x=rand(12)); # master
  531.895 ns (5 allocations: 368 bytes)

julia> @btime sortperm(x) setup=(x=rand(12)); # PR
  574.639 ns (7 allocations: 688 bytes)

[LilithHafner]

I would expect only a single allocation on master for sortperm and none for sortperm! if in-place InsertionSort is used.

There is some type instability in sortperm (separated from the sorting by function barriers) This accounts for 3 of the allocations. Also, the dispatch to InsertionSort for sizes smaller than 40 is only present for UIntMappable types and orders, and Perm orderings are not UIntMappable so we go through the new allocating QuickSort. This accounts for the 4th allocation, and one allocation is expected. Unfortunately, there is a substantial regression in 1.9 compared to 1.8 on small sortperm. This would be covered by BaseBenchmarks after JuliaCI/BaseBenchmarks.jl#305 & after it stops failing. I'll open an issue.

[petvana]

I'll evaluate PR properly once allocation regression (#47811) is fixed to measure the influence of extra allocations and optimize it.

[petvana]

@nanosoldier runbenchmarks("sort")

[nanosoldier]

Your benchmark job has completed - successfully executed benchmarks. A full report can be found here.

[petvana]

@nanosoldier runbenchmarks("sort", vs = ":master")

[nanosoldier]

Your benchmark job has completed - no performance regressions were detected. A full report can be found here.

[LilithHafner]

I benchmarked this with JuliaCI/BaseBenchmarks.jl#305 and found a minor regression for

@benchmark sort(x; by=x -> x isa Symbol ? (0, x) : (1, x)) setup=(x=[rand() < .5 ? randstring() : Symbol(randstring()) for _ in 1:30])

but several more substantial improvements to other benchmarks. It would be nice to know what causes the regression for non-concrete eltypes.

[LilithHafner]

It would be nice to have code re-use; the loop here is almost exactly the same as the loops in partition!

[petvana]

It would be nice. However, partition! left a pivot on the place and split to two parts with. Thus it seems challenging.

[LilithHafner]

Ideally, this should utilize the scratch space re-use system, though it this PR is a performance improvement without, then that isn't essential.

[petvana]

I've tried to re-use the scratch array, but not sure if i got your @getkw concept right.

[petvana]

@LilithHafner Btw, is scratch type-stable here, i.e., can compiler infer that scratch cannot be Nothing?

[LilithHafner]

I recall that removing the if statement entirely introduced dynamic dispatch. IICU, all of type inference is an implementation detail, so I'm not totally sure, but I believe that we need an if statement to force union splitting though it works well whether we use _sort!(..., nothing) or _sort!(..., scratch) because the compiler can determine the type of scratch at compile time..

[petvana]

@nanosoldier runbenchmarks("sort", vs = ":master")

[nanosoldier]

Your benchmark job has completed - possible performance regressions were detected. A full report can be found here.

[petvana]

The more I think about it, the more I'm convinced the whole Perm concept is inefficient because it generates an extreme number of cache misses. It can be easily beaten by one-liner sp for large arrays (without float-specific optimizations I believe).

julia> x = rand(10000000);

julia> sp(x) = sort(collect(a for a in zip(x, eachindex(x))), by = x -> x[1]) .|> x -> x[2]
sp (generic function with 1 method)

julia> @time sortperm(x); # This PR, already compiled
  1.389301 seconds (7 allocations: 152.588 MiB, 1.04% gc time)

julia> @time sp(x);
  1.182990 seconds (8 allocations: 534.058 MiB, 1.73% gc time)

julia> sp(x) == sortperm(x)
true

[LilithHafner]

+1 for

The more I think about it, the more I'm convinced the whole Perm concept is inefficient