nth smallest element

[Hermann Schwarting]

Hi,

I need to find the nth smallest element of a large matrix. The goal is 
to create a "threshold matrix" where the smallest 10% are 1 and the 
rest is 0. The following works:

% X is a matrix, typical sizes are 2000x2000 to 3000x3000
limit = ceil( numel(X) / 10 );
thresh = sort( X(:) )(limit);
X = X < thresh;

But this takes quite long for large matrices (e.g. if part of a loop). 
A full sort is not needed, so my question is if there is a more 
efficient way to find the nth smallest element?


I already tried to use the nth_element() function from the STL in 
an .oct file. Is it possible in principle? I don’t have experience 
with .oct files, but I threw together the following and it seems to 
work:

#include <octave/oct.h>
#include <algorithm>

DEFUN_DLD( nth_element, args, , "nth smallest element")
{
        int nargs = args.length();
        if( nargs != 2 ) {
                print_usage();
                return octave_value_list();
        }
        Matrix A( args(0).matrix_value() );
        octave_idx_type n( args(1).int_value() );
        if( n < 1 ){ n = 1; }
        if( n > A.numel() ){ n = A.numel(); }
        if( error_state ){
                print_usage ();
                return octave_value_list();
        }

        std::nth_element( &A(0), &A.elem(n-1), &A.elem(A.numel()-1) );
        return octave_value( A(n-1) );
}

In fact I would be surprised if the one line should do it type safe 
and for the general case. What do you think?


Here are some timing measurements:

  size        sort (s)  nth_element (s)
  1474x1474   0.83875   0.069808
  1748x1748   1.23738   0.110414
  1943x1943   1.59727   0.155705
  2191x2191   1.95997   0.175527
  2532x2532   2.78808   0.233461
  2560x2560   2.77560   0.227554
  2967x2967   4.13168   0.268378
  3187x3187   4.38392   0.350276
  3353x3353   4.82643   0.618831
  3531x3531   5.95325   0.610461
  3683x3683   5.66996   0.571548


Kind regards,
Hermann