Find k-th smallest index in an MCMC procedure I am repeatedly (100,000+ times) running the following function to find the kth order statistic. I ran a profiler and noticed this function was most expensive. I am looking for a way to optimize this code for speed. The reason I think it can be done faster is because I do not want find the kth largest value, but the index of the kth largest value. However, now I'm effectively looping twice through the vector x, which seems inefficient to me.
The code below is an example, in practice I use R and Rcpp to run the code. Any ideas?

#include &lt;iostream&gt;
#include "C:\c++\armadillo-7.950.1\include\armadillo"

using namespace std;
using namespace arma;

int qSelectIdxC(arma::vec&amp; x, const int k) {

    // ARGUMENTS
    // x: vector to find k-th largest element in
    // k: k-th statistic to look up

    // safety copy since nth_element modifies in place
    arma::vec y(x.memptr(), x.n_elem);

    // partially sorts y.
    std::nth_element(y.begin(), y.begin() + k, y.end());

    // the k-th largest value
    const double kthValue = y(k);

    // find and return the index of the k-th largest value;
    int idxK = std::find(x.begin(), x.end(), kthValue) - x.begin();
    return idxK;

}

int main() {

    vec test0 = regspace&lt;vec&gt;(0, 10); // 0, 1, ..., 10
    int ans0 = qSelectIdxC(test0, 5); // returns 5

    vec test1(100, fill::randu);
    int ans1 = qSelectIdxC(test1, 50);

    cout &lt;&lt; "ans0" &lt;&lt; ans0;
    cout &lt;&lt; "ans1" &lt;&lt; ans1;

    return 0;
}
To find the kth smallest/largest value you do not need to sort all values! Simply keep track of the k smallest/largest encountered values in a (sorted) buffer. This will take your run time from O(nlog(n)) to O(nlog(k)). And will also avoid the unnecessary copy and iteration, simply keep the index with k.

The following pseudocode illustrates the algorithm

template&lt;typename It&gt;
It min_k(It first, It last, int k){
   auto cmp_it_values = [](It lt, It rt){ return *lt &lt; *rt;};
   auto max_copy = std::min&lt;long&gt;(k, std::distance(first, last));
   auto start_it = first;
   std::advance(start_it, max_copy);

    k++; // k == 0 has to return smallest one element.
    std::vector&lt;It&gt; k_smallest;
    k_smallest.reserve(k+1);
    for(auto it = first; it != start_it; ++it){
        k_smallest.push_back(it);
    }
    std::sort(k_smallest.begin(), k_smallest.end(), cmp_it_values);

    for(auto it = start_it; it != last; ++it){
        if(k_smallest.empty() || *it &lt; *k_smallest.back()){
            auto insertion_point = std::lower_bound(k_smallest.begin(), k_smallest.end(), 
                                                    it, cmp_it_values);
            k_smallest.insert(insertion_point, it);
            if(k_smallest.size() &gt; k){
                k_smallest.pop_back(); // Remove the largest value
            }
        }
    }    
    return k_smallest.back(); // The iterator to the min(k, n) smallest value 
}


The above has O(nklog(k)) worst case and O(n + log(k)) best case behaviour. You can improve the worst case by using a max heap for k_smallest instead of a vector. This would have the promised O(nlog(k)) run time. I leave that as an exercise to the reader ;)

And you should read this using namespace std; is bad practice.

You can see from the benchmark here: https://ideone.com/B81Hs4 that the code above is twice as fast as your original code for the given values (you need to test with typical values for your application to verify) and the heap version is even faster.