В целом, если вы можете использовать векторизованные функции, вы обнаружите, что это (почти) так же быстро, как запуск кода непосредственно в R cpp. Это связано с тем, что многие векторизованные функции в R (почти все векторизованные функции в Base R) написаны на C, Cpp или на Fortran, и поэтому зачастую получить мало что можно.
Тем не менее, есть улучшения для вашего кода R
и Rcpp
. Оптимизация происходит благодаря тщательному изучению кода и удалению ненужных шагов (назначение памяти, суммы и т. Д. c.).
Давайте начнем с оптимизации кода Rcpp
.
В вашем случае основной оптимизацией является удаление ненужных матричных и векторных вычислений. Код по сути:
- Сдвиг бета
- вычисление логарифма суммы exp (сдвиг бета) [log-sum-exp]
- использование Obs в качестве индекс для сдвинутой беты и сумма по всем вероятностям
- вычитают log-sum-exp
Используя это наблюдение, мы можем уменьшить ваш код до 2 for-циклов. Обратите внимание, что sum
- это просто другое значение для -l oop (более или менее: for(i = 0; i < max; i++){ sum += x }
), поэтому избегание сумм может еще больше ускорить код (в большинстве случаев это ненужная оптимизация!). Кроме того, ваш ввод Obs
является целочисленным вектором, и мы можем дополнительно оптимизировать код, используя тип IntegerVector
, чтобы избежать приведения элементов double
к значениям integer
(Кредит Ральфа Стубнера).
cppFunction('double llmnl_int_C_v2(NumericVector beta, IntegerVector Obs, int n_cat)
{
int n_Obs = Obs.size();
NumericVector betas = (beta.size()+1);
//1: shift beta
for (int i = 1; i < n_cat; i++) {
betas[i] = beta[i-1];
};
//2: Calculate log sum only once:
double expBetas_log_sum = log(sum(exp(betas)));
// pre allocate sum
double ll_sum = 0;
//3: Use n_Obs, to avoid calling Xby.size() every time
for (int i = 0; i < n_Obs; i++) {
ll_sum += betas(Obs[i] - 1.0) ;
};
//4: Use that we know denom is the same for all I:
ll_sum = ll_sum - expBetas_log_sum * n_Obs;
return ll_sum;
}')
Обратите внимание, что я удалил довольно много выделений памяти и удалил ненужные вычисления в for-l oop. Также я использовал, что denom
одинаково для всех итераций и просто умножено для конечного результата.
Мы можем выполнить аналогичные оптимизации в вашем R-коде, что приводит к следующей функции:
llmnl_int_R_v2 <- function(beta, Obs, n_cat) {
n_Obs <- length(Obs)
betas <- c(0, beta)
#note: denom = log(sum(exp(betas)))
sum(betas[Obs]) - log(sum(exp(betas))) * n_Obs
}
Обратите внимание, что сложность функции значительно уменьшена, что облегчает чтение другими пользователями. Просто чтобы быть уверенным, что я где-то не ошибся в коде, давайте проверим, чтобы они возвращали одинаковые результаты:
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))
beta = c(4,2,1)
Obs = mnl_sample
n_cat = 4
xr <- llmnl_int(beta = beta, Obs = mnl_sample, n_cat = n_cat)
xr2 <- llmnl_int_R_v2(beta = beta, Obs = mnl_sample, n_cat = n_cat)
xc <- llmnl_int_C(beta = beta, Obs = mnl_sample, n_cat = n_cat)
xc2 <- llmnl_int_C_v2(beta = beta, Obs = mnl_sample, n_cat = n_cat)
all.equal(c(xr, xr2), c(xc, xc2))
TRUE
хорошо, это облегчение.
Производительность:
Я буду использовать микробенчмарк для иллюстрации производительности. Оптимизированные функции работают быстро, поэтому я буду запускать функции 1e5
раза, чтобы уменьшить влияние сборщика мусора
microbenchmark("llmml_int_R" = llmnl_int(beta = beta, Obs = mnl_sample, n_cat = n_cat),
"llmml_int_C" = llmnl_int_C(beta = beta, Obs = mnl_sample, n_cat = n_cat),
"llmnl_int_R_v2" = llmnl_int_R_v2(beta = beta, Obs = mnl_sample, n_cat = n_cat),
"llmml_int_C_v2" = llmnl_int_C_v2(beta = beta, Obs = mnl_sample, n_cat = n_cat),
times = 1e5)
#Output:
#Unit: microseconds
# expr min lq mean median uq max neval
# llmml_int_R 202.701 206.801 288.219673 227.601 334.301 57368.902 1e+05
# llmml_int_C 250.101 252.802 342.190342 272.001 399.251 112459.601 1e+05
# llmnl_int_R_v2 4.800 5.601 8.930027 6.401 9.702 5232.001 1e+05
# llmml_int_C_v2 5.100 5.801 8.834646 6.700 10.101 7154.901 1e+05
Здесь мы видим тот же результат, что и раньше. Теперь новые функции примерно в 35 раз быстрее (R) и в 40 раз быстрее (Cpp) по сравнению с их первыми аналогами. Интересно, что оптимизированная функция R
на все еще немного (на 0,3 мс или 4%) быстрее, чем моя оптимизированная функция Cpp
. Моя лучшая ставка здесь заключается в том, что в пакете Rcpp
есть некоторые издержки, и если он будет удален, они будут идентичны или R.
Аналогичным образом мы можем проверить производительность с помощью Optim.
microbenchmark("llmnl_int" = optim(beta, llmnl_int, Obs = mnl_sample,
n_cat = n_cat, method = "BFGS", hessian = F,
control = list(fnscale = -1)),
"llmnl_int_C" = optim(beta, llmnl_int_C, Obs = mnl_sample,
n_cat = n_cat, method = "BFGS", hessian = F,
control = list(fnscale = -1)),
"llmnl_int_R_v2" = optim(beta, llmnl_int_R_v2, Obs = mnl_sample,
n_cat = n_cat, method = "BFGS", hessian = F,
control = list(fnscale = -1)),
"llmnl_int_C_v2" = optim(beta, llmnl_int_C_v2, Obs = mnl_sample,
n_cat = n_cat, method = "BFGS", hessian = F,
control = list(fnscale = -1)),
times = 1e3)
#Output:
#Unit: microseconds
# expr min lq mean median uq max neval
# llmnl_int 29541.301 53156.801 70304.446 76753.851 83528.101 196415.5 1000
# llmnl_int_C 36879.501 59981.901 83134.218 92419.551 100208.451 190099.1 1000
# llmnl_int_R_v2 667.802 1253.452 1962.875 1585.101 1984.151 22718.3 1000
# llmnl_int_C_v2 704.401 1248.200 1983.247 1671.151 2033.401 11540.3 1000
Еще раз результат тот же.
Вывод:
В качестве краткого заключения стоит отметить, что это один из примеров, где преобразование вашего кода в R cpp на самом деле не стоит беспокоиться. Это не всегда так, но часто стоит еще раз взглянуть на вашу функцию, чтобы увидеть, есть ли области вашего кода, где выполняются ненужные вычисления. Особенно в ситуациях, когда используются встроенные векторизованные функции, часто не стоит тратить время на преобразование кода в R cpp. Чаще всего можно увидеть большие улучшения, если использовать for-loops
с кодом, который нельзя легко векторизовать, чтобы удалить for-l oop.