В настоящее время я работаю над байесовским методом, который требует нескольких этапов оптимизации полиномиальной логит-модели на одну итерацию. Я использую optim () для выполнения этих оптимизаций, и целевая функция, написанная на R. Профилирование показало, что optim () является основным узким местом.
Покопавшись, я нашел вопрос, в котором они предполагают, что перекодирование целевой функции Rcppможет ускорить процесс. Я последовал предложению и перекодировал свою целевую функцию Rcpp, но в итоге он оказался медленнее (примерно в два раза медленнее!).
Это был мой первый раз Rcpp(или что-то связанное с C ++), и я не смог найти способ векторизации кода. Есть идеи, как сделать это быстрее?
Tl; dr: текущая реализация функции в Rcpp не такая быстрая, как векторизация R; как сделать это быстрее?
Воспроизводимый пример :
1) Определить целевые функции в Rи Rcpp: логарифмическая вероятность перехвата только полиномиальной модели
library(Rcpp)
library(microbenchmark)
llmnl_int <- function(beta, Obs, n_cat) {
n_Obs <- length(Obs)
Xint <- matrix(c(0, beta), byrow = T, ncol = n_cat, nrow = n_Obs)
ind <- cbind(c(1:n_Obs), Obs)
Xby <- Xint[ind]
Xint <- exp(Xint)
iota <- c(rep(1, (n_cat)))
denom <- log(Xint %*% iota)
return(sum(Xby - denom))
}
cppFunction('double llmnl_int_C(NumericVector beta, NumericVector Obs, int n_cat) {
int n_Obs = Obs.size();
NumericVector betas = (beta.size()+1);
for (int i = 1; i < n_cat; i++) {
betas[i] = beta[i-1];
};
NumericVector Xby = (n_Obs);
NumericMatrix Xint(n_Obs, n_cat);
NumericVector denom = (n_Obs);
for (int i = 0; i < Xby.size(); i++) {
Xint(i,_) = betas;
Xby[i] = Xint(i,Obs[i]-1.0);
Xint(i,_) = exp(Xint(i,_));
denom[i] = log(sum(Xint(i,_)));
};
return sum(Xby - denom);
}')
2) Сравните их эффективность:
## Draw sample from a multinomial distribution
set.seed(2020)
mnl_sample <- t(rmultinom(n = 1000,size = 1,prob = c(0.3, 0.4, 0.2, 0.1)))
mnl_sample <- apply(mnl_sample,1,function(r) which(r == 1))
## Benchmarking
microbenchmark("llmml_int" = llmnl_int(beta = c(4,2,1), Obs = mnl_sample, n_cat = 4),
"llmml_int_C" = llmnl_int_C(beta = c(4,2,1), Obs = mnl_sample, n_cat = 4),
times = 100)
## Results
# Unit: microseconds
# expr min lq mean median uq max neval
# llmnl_int 76.809 78.6615 81.9677 79.7485 82.8495 124.295 100
# llmnl_int_C 155.405 157.7790 161.7677 159.2200 161.5805 201.655 100
3) Теперь звоню им optim:
## Benchmarking with optim
microbenchmark("llmnl_int" = optim(c(4,2,1), llmnl_int, Obs = mnl_sample, n_cat = 4, method = "BFGS", hessian = T, control = list(fnscale = -1)),
"llmnl_int_C" = optim(c(4,2,1), llmnl_int_C, Obs = mnl_sample, n_cat = 4, method = "BFGS", hessian = T, control = list(fnscale = -1)),
times = 100)
## Results
# Unit: milliseconds
# expr min lq mean median uq max neval
# llmnl_int 12.49163 13.26338 15.74517 14.12413 18.35461 26.58235 100
# llmnl_int_C 25.57419 25.97413 28.05984 26.34231 30.44012 37.13442 100
Я был несколько удивлен, что векторизованная реализация в R оказалась быстрее. Реализация более эффективной версии в Rcpp (скажем, с RcppArmadillo?) Может принести какую-либо выгоду? Это лучшая идея, чтобы перекодировать все в Rcpp с помощью оптимизатора C ++?
PS: первая публикация в Stackoverflow!
источник
ObsкакIntegerVectorудаление некоторых бросков.beta) остаются постоянными в течение наблюденийObs. Если бы у нас были предикторы, изменяющиеся во времени, потребовался бы неявный расчетdenomдля каждого из них наObsоснове значения матрицы проектированияX. Тем не менее, я уже реализую ваши предложения в остальной части моего кода с некоторыми действительно хорошими достижениями :). Спасибо @RalfStubner, @Oliver и @thc за ваши очень проницательные ответы! Теперь переходим к следующему узкому месту!for-loopчто даст вам наибольшую выгоду. Также в более общем случае я бы предложил использоватьmodel.matrix(...)для создания матрицы для ввода в ваши функции.Ваша функция C ++ может быть сделана быстрее, используя следующие наблюдения. По крайней мере, первое также может быть использовано с вашей функцией R:
То, как вы рассчитываете,
denom[i]одинаково для всехi. Поэтому имеет смысл использоватьdouble denomи делать это вычисление только один раз. Я также исключаю вычитание этого общего термина в конце.Ваши наблюдения на самом деле являются целочисленным вектором на стороне R, и вы используете их также как целые числа в C ++. Использование
IntegerVectorдля начала делает ненужным много приведения.Вы также можете индексировать с
NumericVectorпомощьюIntegerVectorC ++. Я не уверен, если это помогает производительности, но это делает код немного короче.Еще несколько изменений, которые больше связаны со стилем, чем с производительностью.
Результат:
Для меня эта функция примерно в десять раз быстрее, чем ваша R-функция.
источник
Я могу вспомнить четыре возможных оптимизации ответов Ральфа и Оливера.
(Вы должны принять их ответы, но я просто хотел добавить свои 2 цента).
1) Использовать
// [[Rcpp::export(rng = false)]]в качестве заголовка комментария к функции в отдельном файле C ++. Это приводит к увеличению скорости моей машины на 80%. (Это самое важное предложение из 4).2) Предпочитаю,
cmathкогда это возможно. (В этом случае это не имеет значения).3) Избегайте выделения, когда это возможно, например, не переходите
betaв новый вектор.4) Растянуть цель: использовать
SEXPпараметры, а не векторы Rcpp. (Оставлено как упражнение для читателя). Векторы Rcpp - это очень тонкие оболочки, но они по-прежнему являются оболочками и имеют небольшие накладные расходы.Эти предложения не были бы важны, если бы не тот факт, что вы вызываете функцию в тесном цикле
optim. Поэтому любые накладные расходы очень важны.Скамейка:
v3 - это ответ Оливера
rng=false. v4 с предложениями № 2 и № 3 включены.Функция:
источник