用 Rcpp 优化 R 目标函数变慢，为什么？

Question

用 Rcpp 优化 R 目标函数变慢，为什么？

smi*_*ner 20 optimization r rcpp

我目前正在研究一种贝叶斯方法，该方法需要在每次迭代中对多项式 logit 模型进行多个优化步骤。我正在使用 optim() 来执行这些优化，并使用 R 编写的目标函数。分析显示 optim() 是主要瓶颈。

在挖掘之后，我发现了这个问题，他们建议重新编码目标函数Rcpp可以加快进程。我遵循了建议并用重新编码了我的目标函数Rcpp，但它最终变慢了（大约慢了两倍！）。

这是我第一次使用Rcpp（或任何与 C++ 相关的东西），我无法找到一种对代码进行矢量化的方法。知道如何使它更快吗？

Tl; dr：当前 Rcpp 中函数的实现不如矢量化 R 快；如何让它更快？

一个可重现的例子：

在Rand 中定义目标函数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);
}')

Run Code Online (Sandbox Code Playgroud)

比较它们的效率：

## 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

Run Code Online (Sandbox Code Playgroud)

现在调用它们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

Run Code Online (Sandbox Code Playgroud)

我有点惊讶 R 中的矢量化实现速度更快。在 Rcpp 中实现更高效的版本（例如，使用 RcppArmadillo？）可以产生任何收益吗？使用 C++ 优化器在 Rcpp 中重新编码所有内容是否更好？

Answer 1

Oli*_*ver 13

一般来说，如果您能够使用矢量化函数，您会发现它（几乎）与直接在 Rcpp 中运行您的代码一样快。这是因为 R 中的许多向量化函数（几乎所有 Base R 中的向量化函数）都是用 C、Cpp 或 Fortran 编写的，因此通常没什么好处。

也就是说，您的代码R和Rcpp代码都有改进。优化来自仔细研究代码，并删除不必要的步骤（内存分配、求和等）。

让我们从Rcpp代码优化开始。

在您的情况下，主要优化是删除不必要的矩阵和向量计算。代码本质上是

班次测试版
计算 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;
}')

Run Code Online (Sandbox Code Playgroud)

请注意，我已经删除了很多内存分配并删除了 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
}

Run Code Online (Sandbox Code Playgroud)

请注意，该函数的复杂性已大大降低，使其他人更容易阅读。为了确保我没有在某处搞乱代码，让我们检查它们是否返回相同的结果：

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

Run Code Online (Sandbox Code Playgroud)

嗯，这是一种解脱。

表现：

我将使用微基准测试来说明性能。优化后的函数速度很快，所以我会运行这些函数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

Run Code Online (Sandbox Code Playgroud)

在这里，我们看到了与之前相同的结果。现在，与它们的第一个对应部分相比，新功能大约快 35 倍 (R) 和快 40 倍 (Cpp)。有趣的是，优化后的R函数仍然比我的优化Cpp函数快一点（0.3 毫秒或 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

Run Code Online (Sandbox Code Playgroud)

再一次，结果是一样的。

结论：

作为一个简短的结论，值得注意的是，这是一个示例，其中将您的代码转换为 Rcpp 并不真正值得麻烦。情况并非总是如此，但通常值得再次查看您的函数，以查看您的代码中是否存在执行不必要计算的区域。特别是在使用内置向量化函数的情况下，通常不值得花时间将代码转换为 Rcpp。如果使用for-loops无法轻松矢量化的代码以删除 for 循环，则通常可以看到很大的改进。

正如您在这个玩具示例（仅截距 mnl 模型）中注意到的那样，线性预测变量（“beta”）在观测值“Obs”上保持恒定。如果我们有随时间变化的预测变量，则需要根据设计矩阵“X”的值对每个“Obs”隐式计算“decom”。话虽这么说，我已经在我的代码的其余部分实现了您的建议，并取得了一些非常好的收益:)。感谢@RalfStubner、@Oliver 和@thc 的非常有见地的回复！现在继续我的下一个瓶颈！ (2认同)

Answer 2

thc*_*thc 12

我可以想到对 Ralf 和 Olivers 答案的四个潜在优化。

（你应该接受他们的答案，但我只想加上我的 2 美分）。

1)// [[Rcpp::export(rng = false)]]在单独的 C++ 文件中用作函数的注释头。这导致我机器上的速度提高了约 80%。（这是 4 条建议中最重要的一条）。

2)cmath可能时优先。（在这种情况下，它似乎没有什么区别）。

3) 尽可能避免分配，例如不要转移beta到新向量中。

4) Stretch 目标：使用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

Run Code Online (Sandbox Code Playgroud)

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;
}

Run Code Online (Sandbox Code Playgroud)

Answer 3

Ral*_*ner 11

使用以下观察可以使您的 C++ 函数更快。至少第一个也可以与您的 R 函数一起使用：

您的计算方式denom[i]对于每个i. 因此，使用 adouble denom并且只进行一次计算是有意义的。我还考虑到最后减去这个常用术语。
您的观察结果实际上是 R 端的整数向量，并且您也在 C++ 中将它们用作整数。使用 anIntegerVector开始不需要很多强制转换。
您也可以在 C++ 中NumericVector使用 an来索引 a IntegerVector。我不确定这是否有助于提高性能，但它使代码更短一些。
更多与风格而非性能相关的变化。

结果：

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;
}

Run Code Online (Sandbox Code Playgroud)

对我来说，这个函数大约比你的 R 函数快十倍。

归档时间：	6 年前
查看次数：	974 次
最近记录：	5 年，5 月前