Seems like a bug in Matlab R2015a "sort" routine

a =[24     9    27    14; 1    40     9    33; 17    16    31    35];

[aSort, iSort] = sort(a, 2, 'descend');
% Those should be exactly the same

isequal(a(iSort), aSort)

But they are not!

Lets see what do we have:

aSort =
    27    24    14     9
    40    33     9     1
    35    31    17    16

iSort =
     3     1     4     2
     2     4     3     1
     4     3     1     2

Apparently iSort is totally messed up! It is not global indexing as it should be!!! It should span 1:numel(a), but instead is spans 1:size(a, 2)- only indexes along sorted dimension! Seems like a severe bug to me. One has to play with "ind2sub" and "sub2ind" to get it right...

Answer 1

[Guillaume]

The bug is in your interpretation of the indices returned by sort. As per the documentation of sort in R2015a:

[B,I] = sort(...) additionally returns a collection of index vectors for any of the previous syntaxes. I is the same size as A and describes the arrangement of the elements of A into B along the sorted dimension. For example, if A is a numeric vector, B = A(I).

It behaves exactly as documented. This may not be what you want but it is what it does. It's not particularly difficult to transform the column indices into linear indices:

iSort = sub2ind(size(A), repmat((1:size(A, 1))', 1, size(A, 2)), iSort);

[Nikolay S.]

O, pardon. I was positive it works for matrices, while it says clearly vector. My mistake. It would be nice to have it working for matrices as well... Thanks for the quick answer!

[Guillaume]

Nikolay S. comment moved from Answer to here:

Guillaume thanks for the solution. For the given case it works (though one may find the code a bit complicated). However, what it would be for N dimensional matrix, sorted along dimension K? I have failed to generate code for this general case.

[Guillaume]

Yes, a generic code is a lot more complicated. This would work:

A = randi([1 100], [10 20 15 5]); %demo data
k = 3; %demo data

[B, isort] = sort(A, k);
%build subscript indices arrays with ndgrid:
subs = cell(1, ndims(A)); %store all subscript arrays in a cell array
[subs{:}] = size(A); %size of each dimension
subs = cellfun(@(s) 1:s, subs, 'UniformOutput', false); %convert to vectors of 1:size(dim)
[subs{:}] = ndgrid(subs{:}); %convert to arrays of subscripts
subs{k} = isort; %replace subscript for kth dimension by isort
isort = sub2ind(size(A), subs{:});

isequal(A(isort), B) %should return 1