用 Rcpp 更慢地优化 R 目标函数，为啥？
Posted
技术标签:
【中文标题】用 Rcpp 更慢地优化 R 目标函数，为啥？【英文标题】：Optimizing R objective function with Rcpp slower, why?用 Rcpp 更慢地优化 R 目标函数，为什么？ 【发布时间】：2020-06-02 18:22:11 【问题描述】：我目前正在研究一种贝叶斯方法，该方法需要对每次迭代的多项式 logit 模型进行多个优化步骤。我正在使用 optim() 来执行这些优化，并且使用 R 编写了一个目标函数。分析表明 optim() 是主要瓶颈。
在四处挖掘之后，我发现了this question，他们建议用Rcpp 重新编码目标函数可以加快这个过程。我按照建议使用Rcpp 重新编码了我的目标函数，但它最终变慢了（大约慢了两倍！）。
这是我第一次使用Rcpp（或任何与 C++ 相关的东西），但我无法找到矢量化代码的方法。知道如何让它更快吗？
Tl;dr：Rcpp 中函数的当前实现不如矢量化 R 快；如何让它更快？
一个可重现的例子：
-
在
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);
')
-
比较它们的效率：
## 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
-
现在打电话给他们
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？）可以产生任何收益吗？使用 C++ 优化器重新编码 Rcpp 中的所有内容是否更好？
【问题讨论】：
【参考方案1】：一般来说，如果您能够使用向量化函数，您会发现它（几乎）与直接在 Rcpp 中运行代码一样快。这是因为 R 中的许多向量化函数（Base R 中的几乎所有向量化函数）都是用 C、Cpp 或 Fortran 编写的，因此通常没有什么收获。
也就是说，您的R 和Rcpp 代码都有改进。优化来自仔细研究代码，并删除不必要的步骤（内存分配、求和等）。
让我们从Rcpp 代码优化开始。
在您的情况下，主要优化是删除不必要的矩阵和向量计算。代码本质是
-
Shift 测试版
计算 exp(shift beta) [log-sum-exp] 总和的对数
使用 Obs 作为移动 beta 的索引并对所有概率求和
减去 log-sum-exp
利用这一观察，我们可以将您的代码减少到 2 个 for 循环。请注意，sum 只是另一个 for 循环（或多或少：for(i = 0; i < max; i++) sum += x ），因此避免求和可以进一步加快代码速度（在大多数情况下，这是不必要的优化！）。此外，您的输入 Obs 是一个整数向量，我们可以通过使用 IntegerVector 类型来进一步优化代码，以避免将 double 元素转换为 integer 值（归功于 Ralf Stubner 的回答）。
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 循环中不必要的计算。我也使用了 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 函数仍然比我优化的 Cpp 函数快得多（0.3 ms 或 4 %）。我最好的选择是 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
结果还是一样。
结论：
作为一个简短的结论，值得注意的是，这是一个示例，将代码转换为 Rcpp 并不值得麻烦。情况并非总是如此，但通常值得再次查看您的函数，看看您的代码中是否存在执行不必要计算的区域。尤其是在使用内置矢量化函数的情况下，通常不值得花时间将代码转换为 Rcpp。如果将for-loops 与不容易矢量化的代码一起使用以删除for循环，则通常可以看到很大的改进。
【讨论】：
您可以将Obs 视为IntegerVector 删除一些演员表。
在感谢您在回答中注意到这一点之前将其合并。它只是从我身边经过。我在@RalfStubner 的回答中对此表示赞赏。 :-)
正如您在这个玩具示例（仅拦截 mnl 模型）中注意到的那样，线性预测变量 (beta) 在观测值 Obs 上保持不变。如果我们有随时间变化的预测变量，则需要根据设计矩阵X 的值对每个Obs 进行隐式计算denom。话虽这么说，我已经在我的其余代码中实施了您的建议，并获得了一些非常好的收益:)。感谢@RalfStubner、@Oliver 和@thc 的富有洞察力的回复！现在继续我的下一个瓶颈！
很高兴我们能提供帮助。在更一般的情况下，在第二个for-loop 的每个步骤中计算减法面额，这将为您带来最大的收益。同样在更一般的情况下，我建议使用model.matrix(...) 来创建矩阵以供函数输入。【参考方案2】：
相对于 Ralf 和 Olivers 的答案，我能想到四个潜在优化。
（你应该接受他们的回答，但我只是想加我的 2 美分）。
1) 在单独的 C++ 文件中使用 // [[Rcpp::export(rng = false)]] 作为函数的注释头。这导致我的机器速度提高了约 80%。 （这是 4 个中最重要的建议）。
2) 尽可能首选cmath。 （在这种情况下，它似乎没有什么区别）。
3) 尽可能避免分配，例如不要将beta 转换为新向量。
4) 拉伸目标：使用SEXP 参数而不是Rcpp 向量。 （留给读者作为练习）。 Rcpp 向量是非常薄的包装器，但它们仍然是包装器，并且开销很小。
如果不是因为您在optim 中的紧密循环中调用该函数，这些建议并不重要。所以任何开销都非常重要。
替补：
microbenchmark("llmnl_int_R_v1" = optim(beta, llmnl_int, 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)),
"llmnl_int_C_v3" = optim(beta, llmnl_int_C_v3, Obs = mnl_sample,
n_cat = n_cat, method = "BFGS", hessian = F,
control = list(fnscale = -1)),
"llmnl_int_C_v4" = optim(beta, llmnl_int_C_v4, Obs = mnl_sample,
n_cat = n_cat, method = "BFGS", hessian = F,
control = list(fnscale = -1)),
times = 1000)
Unit: microseconds
expr min lq mean median uq max neval cld
llmnl_int_R_v1 9480.780 10662.3530 14126.6399 11359.8460 18505.6280 146823.430 1000 c
llmnl_int_R_v2 697.276 735.7735 1015.8217 768.5735 810.6235 11095.924 1000 b
llmnl_int_C_v2 997.828 1021.4720 1106.0968 1031.7905 1078.2835 11222.803 1000 b
llmnl_int_C_v3 284.519 295.7825 328.5890 304.0325 328.2015 9647.417 1000 a
llmnl_int_C_v4 245.650 256.9760 283.9071 266.3985 299.2090 1156.448 1000 a
v3 是 Oliver 对rng=false 的回答。 v4 包含建议 #2 和 #3。
功能：
#include <Rcpp.h>
#include <cmath>
using namespace Rcpp;
// [[Rcpp::export(rng = false)]]
double llmnl_int_C_v4(NumericVector beta, IntegerVector Obs, int n_cat)
int n_Obs = Obs.size();
//2: Calculate log sum only once:
// double expBetas_log_sum = log(sum(exp(betas)));
double expBetas_log_sum = 1.0; // std::exp(0)
for (int i = 1; i < n_cat; i++)
expBetas_log_sum += std::exp(beta[i-1]);
;
expBetas_log_sum = std::log(expBetas_log_sum);
double ll_sum = 0;
//3: Use n_Obs, to avoid calling Xby.size() every time
for (int i = 0; i < n_Obs; i++)
if(Obs[i] == 1L) continue;
ll_sum += beta[Obs[i]-2L];
;
//4: Use that we know denom is the same for all I:
ll_sum = ll_sum - expBetas_log_sum * n_Obs;
return ll_sum;
【讨论】：
【参考方案3】：使用以下观察可以使您的 C++ 函数更快。至少第一个也可以与您的 R 函数一起使用：
计算denom[i] 的方式对于每个i 都是相同的。因此使用double denom 并且只进行一次计算是有意义的。最后，我还考虑减去这个常用术语。
您的观察结果实际上是 R 端的整数向量，并且您也在 C++ 中将它们用作整数。使用IntegerVector 开头会导致不必要的大量转换。
您也可以在 C++ 中使用 IntegerVector 索引 NumericVector。我不确定这是否有助于提高性能，但它会使代码更短。
更多的变化与风格相关而不是性能。
结果：
double llmnl_int_C(NumericVector beta, IntegerVector 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];
;
double denom = log(sum(exp(betas)));
NumericVector Xby = betas[Obs - 1];
return sum(Xby) - n_Obs * denom;
对我来说，这个函数比你的 R 函数快大约十倍。
【讨论】：
感谢您的回答拉尔夫，没有发现输入类型。我已将其纳入我的回答中，并给予您信任。 :-)以上是关于用 Rcpp 更慢地优化 R 目标函数，为啥？的主要内容，如果未能解决你的问题，请参考以下文章
在构建 R 包时从另一个 Rcpp 函数调用 Rcpp 函数