Can't get sortperm to return the sorting order of columns of a matrix

General Usage
, ,
#1

I’m sorting the columns of each matrix in long vector of matrices (having the matrix columns in a canonical order makes it easy to identify approximately equal matrices). Here’s a minimal working example:

julia> t
3×2 transpose(::Matrix{Float64}) with eltype Float64:
 0.0  -0.5
 0.0  -0.5
 0.0   0.0

julia> sort(t,dims=2)
3×2 Matrix{Float64}:
 -0.5  0.0
 -0.5  0.0
  0.0  0.0

There are functions associated with each column of each matrix as well, and I need to permute these polynomials in the same way when the columns are sorted. But sortperm does not accept the same arguments that sort does. It seems the dims keyword is unsupported.

julia> sortperm(t,dims=2)
ERROR: MethodError: no method matching sortperm(::Transpose{Float64, Matrix{Float64}}; dims=2)
Closest candidates are:
  sortperm(::AbstractUnitRange) at /Applications/Julia-1.7.app/Contents/Resources/julia/share/julia/base/range.jl:1309 got unsupported keyword argument "dims"
  sortperm(::AbstractRange) at /Applications/Julia-1.7.app/Contents/Resources/julia/share/julia/base/range.jl:1310 got unsupported keyword argument "dims"
  sortperm(::StaticArrays.StaticVector; alg, lt, by, rev, order) at ~/.julia/packages/StaticArrays/58yy1/src/sort.jl:26 got unsupported keyword argument "dims"
  ...
Stacktrace:
 [1] top-level scope
   @ REPL[1417]:1

Am I missing something? Anyone have ideas for a simple workaround?

I would assume (naively?) that sort and sortperm are the same function under the hood and just return different parts of the sort results (sorted items vs. order of sort). EDIT: for the first time ever, I looked at the source code in Base. It seems that sortperm and sort are separate functions…still, I really need a work around for this…

julia> versioninfo()
Julia Version 1.7.1
Commit ac5cc99908 (2021-12-22 19:35 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin19.5.0)
  CPU: Apple M1 Pro
#2

sortperm is a generalized method for a 1-D sort. What we need to do then is obtain either an array of rows or an array of columns to sort over.

There are functions eachrow, eachcol, and eachslice that are generators. To turn these into arrays we need to collect them.

julia> A = rand(5, 5)
5×5 Matrix{Float64}:
 0.212887  0.505692  0.656134  0.23317   0.607276
 0.535299  0.36777   0.595655  0.168758  0.225478
 0.934146  0.98104   0.556343  0.473143  0.836663
 0.216165  0.253191  0.854917  0.165115  0.44384
 0.339784  0.940717  0.35101   0.169418  0.271189

julia> rowsperm = sortperm(collect(eachslice(A; dims = 1)))
5-element Vector{Int64}:
 1
 4
 5
 2
 3

julia> colsperm = sortperm(collect(eachslice(A; dims = 2)))
5-element Vector{Int64}:
 1
 4
 2
 5
 3

julia> @view A[rowsperm,:]
5×5 view(::Matrix{Float64}, [1, 4, 5, 2, 3], :) with eltype Float64:
 0.212887  0.505692  0.656134  0.23317   0.607276
 0.216165  0.253191  0.854917  0.165115  0.44384
 0.339784  0.940717  0.35101   0.169418  0.271189
 0.535299  0.36777   0.595655  0.168758  0.225478
 0.934146  0.98104   0.556343  0.473143  0.836663

julia> @view A[:,colsperm]
5×5 view(::Matrix{Float64}, :, [1, 4, 2, 5, 3]) with eltype Float64:
 0.212887  0.23317   0.505692  0.607276  0.656134
 0.535299  0.168758  0.36777   0.225478  0.595655
 0.934146  0.473143  0.98104   0.836663  0.556343
 0.216165  0.165115  0.253191  0.44384   0.854917
 0.339784  0.169418  0.940717  0.271189  0.35101
1 Like
#3

Also see Sortperm for matrix sorting in Julia lang - Stack Overflow

1 Like
#4

If you are wondering why sortperm does not have a dims parameter, the answer is that Julia will have a dims parameter in Julia version 1.9 since @pjentsch0 just implemented it:

It looks like sortslicesperm is on the way also due to the efforts of jgonzalez49

3 Likes
#5

I modified the example you referenced (https://stackoverflow.com/questions/68344823/sortperm-for-matrix-sorting-in-julia-lang so that the permutations would be non-trivial for a sort over columns and so that multiple components of the vectors would be needed to break ties. Your solution was just what I needed. Thanks.

julia> B =
       [ 2 4 4 4 5 ;
         1 2 2 3 5 ;
         1 2 3 3 3 ;
         1 2 2 5 6 ;
         1 3 4 4 4 ; ]|>transpose
5×5 transpose(::Matrix{Int64}) with eltype Int64:
 2  1  1  1  1
 4  2  2  2  3
 4  2  3  2  4
 4  3  3  5  4
 5  5  3  6  4

julia> sortperm(collect(eachslice(B,dims=2)))
5-element Vector{Int64}:
 2
 4
 3
 5
 1
1 Like
#6

Something I didn’t understand at first—a fundamental misunderstanding—sort(t,dims=2) does not preserve the columns necessarily. It sorts the columns of each row separately. What I really wanted was to sort the entire columns (preserving their contents).

julia> t
3×2 Matrix{Float64}:
 0.5  0.5
 0.5  0.0
 0.0  0.5

julia> sort(t,dims=2) # Doesn't preserve the columns
3×2 Matrix{Float64}:
 0.5  0.5
 0.0  0.5
 0.0  0.5

julia> sort(collect(eachslice(t,dims=2))) # Does preserve the columns
2-element Vector{SubArray{Float64, 1, Matrix{Float64}, Tuple{Base.Slice{Base.OneTo{Int64}}, Int64}, true}}:
 [0.5, 0.0, 0.5]
 [0.5, 0.5, 0.0]
1 Like
#7

Have you looked at sortslices()?

julia> sortslices(t, dims=2)
3×2 Matrix{Float64}:
 0.5  0.5
 0.0  0.5
 0.5  0.0
2 Likes
#8

Super helpful! Thanks. I wasn’t aware of sortslices. That makes my sorting syntax simpler (though I still need the sortperm too but @mkitti shared above that a perm version is coming).

#9

Per @Lilith , perhaps it is not coming:

#10

We don’t need sortslicesperm because we already have eachslice and sortperm, which @Gus_Hart put together to achieve the same result above:

As of 1.9.0, collect is unnecessary here because the result of eachslice is already indexable:

julia> @btime sortperm(eachslice(B,dims=2))
  316.099 ns (5 allocations: 176 bytes)
5-element Vector{Int64}:
 2
 4
 3
 5
 1

julia> @btime sortperm(collect(eachslice(B,dims=2)))
  522.126 ns (6 allocations: 432 bytes)
5-element Vector{Int64}:
 2
 4
 3
 5
 1
#11

It seems to be necessary at least with Julia 1.7.3. Which version are you using?

1 Like
#12

Good catch! I’m on 1.9.0-DEV.925. I imagine that collect became unnecessary after add Slices array type for eachslice/eachrow/eachcol by simonbyrne · Pull Request #32310 · JuliaLang/julia · GitHub, which would make 1.9.0-DEV the cutoff. I edited my previous post to reflect that.

2 Likes