Почему стандартная функция R-медианы намного медленнее, чем простая альтернатива C ++?
Я сделал следующую реализацию медианы в C++ и использовал его в R с помощью Rcpp:
// [[Rcpp::export]]
double median2(std::vector<double> x){
double median;
size_t size = x.size();
sort(x.begin(), x.end());
if (size % 2 == 0){
median = (x[size / 2 - 1] + x[size / 2]) / 2.0;
}
else {
median = x[size / 2];
}
return median;
}
Если впоследствии я сравниваю производительность со стандартной встроенной медианной функцией R, я получаю следующие результаты через microbenchmark
> x = rnorm(100)
> microbenchmark(median(x),median2(x))
Unit: microseconds
expr min lq mean median uq max neval
median(x) 25.469 26.990 34.96888 28.130 29.081 518.126 100
median2(x) 1.140 1.521 2.47486 1.901 2.281 47.897 100
Почему стандартная медианная функция намного медленнее? Это не то, что я ожидал …
Решение
Как отметил @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"Редактировать:
Я должен отметить, что функция ниже приведет к сортировке исходного вектора поскольку NumericVectors являются прокси-объектами. Если вы хотите избежать этого, вы можете либо 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
Якк поднял замечательное замечание в комментариях выше — также разработанных Джерри Коффином — о неэффективности выполнения полной сортировки. Вот переписать с помощью 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;
Я не уверен, на какую «стандартную» реализацию вы бы ссылались.
В любом случае: если бы он был, он, будучи частью стандартной библиотеки, безусловно, не имел бы права изменять порядок элементов в векторе (как это делает ваша реализация), поэтому он определенно должен был бы работать с копией.
Создание этой копии потребует времени и ресурсов процессора (и значительного объема памяти), что повлияет на время выполнения.