Почему стандартная медианная функция R намного медленнее, чем простая альтернатива C++?

Как отметил @joran, ваш код очень специализирован, и, вообще говоря, менее обобщенные функции, алгоритмы и т. д. часто более эффективны. Взгляните на median.default:

median.default
# function (x, na.rm = FALSE) 
# {
#   if (is.factor(x) || is.data.frame(x)) 
#     stop("need numeric data")
#   if (length(names(x))) 
#     names(x) <- NULL
#   if (na.rm) 
#     x <- x[!is.na(x)]
#   else if (any(is.na(x))) 
#     return(x[FALSE][NA])
#   n <- length(x)
#   if (n == 0L) 
#     return(x[FALSE][NA])
#   half <- (n + 1L)%/%2L
#   if (n%%2L == 1L) 
#     sort(x, partial = half)[half]
#   else mean(sort(x, partial = half + 0L:1L)[half + 0L:1L])
# }

Существует несколько операций для учета возможности отсутствия значений, и они определенно повлияют на общее время выполнения функции. Поскольку ваша функция не повторяет это поведение, она может исключить кучу вычислений, но, следовательно, не даст такого же результата для векторов с отсутствующими значениями:

median(c(1, 2, NA))
#[1] NA

median2(c(1, 2, NA))
#[1] 2

Пара других факторов, которые, вероятно, не имеют такого эффекта, как обработка NA, но на них стоит обратить внимание:

• median вместе с несколькими функциями, которые он использует, являются дженериками S3, поэтому на отправку методов тратится небольшое количество времени.
• median будет работать не только с целочисленными и числовыми векторами; он также будет обрабатывать Date, POSIXt и, возможно, кучу других классов и правильно сохранять атрибуты:

median(Sys.Date() + 0:4)
#[1] "2016-01-15"

median(Sys.time() + (0:4) * 3600 * 24)
#[1] "2016-01-15 11:14:31 EST"

Редактировать: я должен упомянуть, что функция ниже вызовет сортировку исходного вектора, поскольку NumericVector являются прокси-объектами. Если вы хотите избежать этого, вы можете либо Rcpp::clone ввести входной вектор и работать с клоном, либо использовать исходную подпись (с std::vector<double>), которая неявно требует копии при преобразовании из SEXP в std::vector.

Также обратите внимание, что вы можете сэкономить немного больше времени, используя NumericVector вместо std::vector<double>:

#include <Rcpp.h>

// [[Rcpp::export]]
double cpp_med(Rcpp::NumericVector x){
  std::size_t size = x.size();
  std::sort(x.begin(), x.end());
  if (size  % 2 == 0) return (x[size / 2 - 1] + x[size / 2]) / 2.0;
  return x[size / 2];
}

microbenchmark::microbenchmark(
  median(x),
  median2(x),
  cpp_med(x),
  times = 200L
)
# Unit: microseconds
#       expr    min      lq      mean  median      uq     max neval
#  median(x) 74.787 81.6485 110.09870 92.5665 129.757 293.810   200
# median2(x)  6.474  7.9665  13.90126 11.0570  14.844 151.817   200
# cpp_med(x)  5.737  7.4285  11.25318  9.0270  13.405  52.184   200

Yakk поднял в комментариях выше важную мысль, которую также подробно изложил Джерри Коффин, о неэффективности выполнения полной сортировки. Вот переписывание с использованием std::nth_element, проверенное на гораздо большем векторе:

#include <Rcpp.h>

// [[Rcpp::export]]
double cpp_med2(Rcpp::NumericVector xx) {
  Rcpp::NumericVector x = Rcpp::clone(xx);
  std::size_t n = x.size() / 2;
  std::nth_element(x.begin(), x.begin() + n, x.end());

  if (x.size() % 2) return x[n]; 
  return (x[n] + *std::max_element(x.begin(), x.begin() + n)) / 2.;
}

set.seed(123)
xx <- rnorm(10e5)

all.equal(cpp_med2(xx), median(xx))
all.equal(median2(xx), median(xx))

microbenchmark::microbenchmark(
  cpp_med2(xx), median2(xx), 
  median(xx), times = 200L
)
# Unit: milliseconds
#         expr      min       lq     mean   median       uq       max neval
# cpp_med2(xx) 10.89060 11.34894 13.15313 12.72861 13.56161  33.92103   200
#  median2(xx) 84.29518 85.47184 88.57361 86.05363 87.70065 228.07301   200
#   median(xx) 46.18976 48.36627 58.77436 49.31659 53.46830 250.66939   200

[Это скорее развернутый комментарий, чем ответ на заданный вами вопрос.]

Даже ваш код может быть открыт для значительных улучшений. В частности, вы сортируете весь ввод, даже если вас интересуют только один или два элемента.

Вы можете изменить это с O(n log n) на O(n), используя std::nth_element вместо std::sort. В случае четного числа элементов вы обычно хотите использовать std::nth_element, чтобы найти элемент непосредственно перед серединой, а затем использовать std::min_element, чтобы найти следующий за ним элемент, но std::nth_element также разделяет входные элементы, поэтому std::min_element имеет только запускать элементы выше середины после nth_element, а не весь входной массив. То есть после nth_element получается такая ситуация:

Сложность std::nth_element "линейна в среднем", и (конечно) std::min_element тоже линейна, так что общая сложность линейна.

Итак, для простого случая (нечетное количество элементов) вы получите что-то вроде:

auto pos = x.begin() + x.size()/2;

std::nth_element(x.begin(), pos, x.end());
return *pos;

... и для более сложного случая (четное количество элементов):

std::nth_element(x.begin(), pos, x.end());
auto pos2 = std::min_element(pos+1, x.end());
return (*pos + *pos2) / 2.0;

Из здесь можно ожидать, что max_element( ForwardIt first, ForwardIt last ) обеспечивает максимальное значение от первого до последнего, но делая: return (x[n] + *std::max_element(x.begin(), x.begin() + n)) / 2. элемент x.begin() + n кажется исключенным из расчета. Почему такое несоответствие?

Например. cpp_med2({6, 2, 1, 5, 3, 4}) производит x={2, 1, 3, 4, 5, 6} где:

n = 3
*x[n] = 4
*x.begin() = 2
*(x.begin() + n) = 4
*std::max_element(x.begin(), x.begin() + n) = 3

Так что cpp_med2({6, 2, 1, 5, 3, 4}) возвращает (4+3)/2=3,5, что является правильной медианой. Но почему *std::max_element(x.begin(), x.begin() + n) равно 3, а не 4? функция фактически исключает последний элемент (4) в максимальном расчете.

РЕШЕНО (я думаю): в:

Находит наибольший элемент в диапазоне [first, last)

) закрытие означает, что последнее исключается из расчета. Это правильно?

С наилучшими пожеланиями

Я не уверен, какую «стандартную» реализацию вы бы имели в виду.

В любом случае: если бы он был, он, будучи частью стандартной библиотеки, определенно не мог бы изменять порядок элементов в векторе (как это делает ваша реализация), поэтому ему определенно пришлось бы работать над копией.

Создание этой копии потребует времени и ЦП (и значительной памяти), что повлияет на время выполнения.