为什么Boost的实现比R慢5-10倍？

Question

我正在开发一个应用程序，经常计算正则化不完全贝塔函数。该应用程序使用C++编写，并调用R::pbeta()。当我尝试将应用程序多线程化时，一些来自R::pbeta()的警告消息导致堆栈溢出。

因此，我转向了boost::math::ibeta()。一切都很顺利，直到我测量速度为止。下面的C++文件whyIsBoostSlower.cpp使用R::pbeta()或boost::math::ibeta()实现了正则化不完全贝塔函数。

// [[Rcpp::plugins(cpp17)]]
#include <boost/math/special_functions/beta.hpp>
// [[Rcpp::depends(BH)]]
#include <Rcpp.h>
using namespace Rcpp;


// Compute the regularized incomplete Beta function.
// [[Rcpp::export]]
NumericVector RIBF(NumericVector q, NumericVector a, NumericVector b, 
                   bool useboost = false)
{
  NumericVector rst(q.size());
  for (int i = 0, iend = q.size(); i < iend; ++i)
  {
    if (useboost) rst[i] = boost::math::ibeta( a[i], b[i], q[i] );
    else          rst[i] = R::pbeta( q[i], a[i], b[i], 1, 0 );
  }
  return rst;
}

在R中，我们通过对随机数调用函数300000次来测量速度。

Rcpp::sourceCpp("whyIsBoostSlower.cpp")


set.seed(123)
N = 300000L
q = runif(N) # Generate quantiles.
a = runif(N, 0, 10) # Generate a in (0, 10)
b = runif(N, 0, 10) # Generate b in (0, 10)


# Use R's pbeta(). This function calls a C wrapper of toms708.c:
#   https://svn.r-project.org/R/trunk/src/nmath/toms708.c
system.time({ Rrst = RIBF(q, a, b, useboost = F) })
# Windows 10 (seconds):
# user  system elapsed 
# 0.11    0.00    0.11 

# RedHat Linux:
# user  system elapsed 
# 0.097   0.000   0.097 


# Use Boost's implementation, which also depends on TOMS 708 by their claim:
#  https://www.boost.org/doc/libs/1_41_0/libs/math/doc/sf_and_dist/html/math_toolkit/special/sf_beta/ibeta_function.html
system.time({ boostRst = RIBF(q, a, b, useboost = T) })
# Windows 10:
# user  system elapsed 
# 0.52    0.00    0.52 

# RedHat Linux:
# user  system elapsed 
# 0.988   0.001   0.993 


range(Rrst - boostRst)
# -1.221245e-15  1.165734e-15

为了复现这个例子，需要安装R和包Rcpp。在Windows上，还需要安装包含GCC分发的Rtools。优化标志默认为-O2。 R::pbeta()和boost::math::ibeta()都基于ACM TOMS 708，但是boost::math::ibeta()在Windows上慢5倍，在Linux上慢10倍。

我认为这可能与在boost::math::ibeta()中设置Policy参数有关，但是如何设置呢？

谢谢！

FYI，R::pbeta()在R-4.2.3/src/nmath/pbeta.c中定义。 R::pbeta()调用了bratio()，该函数在R-4.2.3/src/nmath/toms708.c中定义，即https://svn.r-project.org/R/trunk/src/nmath/toms708.c。其中的代码是TOMS 708 Fortran代码的C翻译版本。这个翻译是由R的核心团队完成的。

相比之下，Boost声称“这个实现是基于1992年DiDonato和Morris在ACM上发表的《算法708：不完全贝塔函数比率的有效数字计算》。”（boost::math::ibeta()）

- user2961927

---

实际上，请提供与您的问题相匹配的测量和数据。 - πάντα ῥεῖ

你有没有检查启用优化的C++发布版本？你的测试集足够大吗？换句话说，没有上下文的测量是毫无意义的。而且你有没有检查你使用的库是否是线程安全的？ - Pepijn Kramer

@PepijnKramer 我修改了问题，使得结果可以轻松复现。 - user2961927

你的运行时间是否可以由Rcpp接口主导？ - Eugene

@Eugene 作为Rcpp，它确保适当的状态保存和重置RNG，包装tryCatch等，因此会有一些非常小的开销，但这是微不足道的成本，通常不会影响任何实际运行时成本（除非你测量病态的空函数）。另外请注意，在这里我们有接近3000个与Rcpp标签相关的问题（尽管我努力获取无意义的问题或重复的问题并进行删除）。还有一些与Boost和R的BH包交集的问题。 - Dirk Eddelbuettel

Answer 1

由于存在Beta、不完整Beta以及Beta分布（可以使用R::pbeta()函数进行计算），我希望确保我们在比较时是公平的。

因此，这里有一个修改过的版本的代码，为了简单起见，使用了两个不同的函数，并与GSL进行了比较，并进行了正式的基准测试调用：

代码

// [[Rcpp::depends(BH)]]
#include <boost/math/special_functions/beta.hpp>

// this also ensure linking with the GSL
// [[Rcpp::depends(RcppGSL)]]
#include <gsl/gsl_sf_gamma.h>

#include <Rcpp.h>
using namespace Rcpp;


// [[Rcpp::export]]
NumericVector bfR(NumericVector a, NumericVector b) {
    int n = a.size();
    NumericVector rst(n);
    for (int i = 0; i<n; i++) {
        rst[i] = R::beta(a[i], b[i]);
    }
    return rst;
}

// [[Rcpp::export]]
NumericVector bfB(NumericVector a, NumericVector b) {
    int n = a.size();
    NumericVector rst(n);
    for (int i = 0; i<n; i++) {
        rst[i] = boost::math::beta( a[i], b[i] );
    }
    return rst;
}

// [[Rcpp::export]]
NumericVector bfG(NumericVector a, NumericVector b) {
    int n = a.size();
    NumericVector rst(n);
    for (int i = 0; i<n; i++) {
        rst[i] = gsl_sf_beta( a[i], b[i] );
    }
    return rst;
}


/*** R

set.seed(123)
N <- 300000L
a <- runif(N, 0, 10) # Generate a in (0, 10)
b <- runif(N, 0, 10) # Generate b in (0, 10)
summary(bfR(a,b) - bfB(a,b))
summary(bfR(a,b) - bfG(a,b))
microbenchmark::microbenchmark(R = bfR(a, b), Boost = bfB(a, b), GSL = bfG(a, b), times=10)

*/

输出

当我们使用Rcpp::sourceCpp()时，'marked' R 代码部分也会被执行：

> Rcpp::sourceCpp("answer.cpp")

> set.seed(123)

> N <- 300000L

> a <- runif(N, 0, 10) # Generate a in (0, 10)

> b <- runif(N, 0, 10) # Generate b in (0, 10)

> summary(bfR(a,b) - bfB(a,b))
     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
-3.64e-12  0.00e+00  0.00e+00 -5.00e-17  0.00e+00  1.36e-12 

> summary(bfR(a,b) - bfG(a,b))
     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
-1.18e-11 -4.00e-17  0.00e+00  2.10e-16  0.00e+00  1.09e-11 

> microbenchmark::microbenchmark(R = bfR(a, b), Boost = bfB(a, b), GSL = bfG(a, b), times=10)
Unit: milliseconds
  expr      min       lq     mean   median       uq      max neval cld
     R  44.9314  45.2773  46.9782  46.1237  49.0056  50.6273    10 a  
 Boost 166.0146 167.2552 171.0441 169.5741 175.0108 180.6520    10  b 
   GSL  58.3259  58.5101  61.0364  59.6556  62.4862  67.5316    10   c
> 

目前我只能猜测Boost可能要经历一些额外的步骤，或者在抽象层面上承担了一些成本，因此在准确性上输给了R和GSL。（我注意到在Boost的文档页面上，结果（准确性）与GNU GSL以及R进行了比较：https://www.boost.org/doc/libs/1_82_0/libs/math/doc/html/math_toolkit/sf_beta/beta_function.html）

- Dirk Eddelbuettel

---

谢谢Dirk！链接很有帮助。所以我猜可能是因为Boost运行更长的近似级数，以将误差率降低到比R更低的水平。 - user2961927

没问题。可能值得澄清一下你想要的是什么：贝塔分布，还是函数，或者其他什么东西？ - Dirk Eddelbuettel

我之前是在寻找正式的不完全Beta函数，这恰好是R中的Beta CDF。你的代码测量了“纯粹”的Beta函数的速度，但我认为Boost速度较慢的原因是相同的。对于我的工作来说，Boost提供的极高精度并不值得那5-10倍的减速。我仍在寻找一个线程安全的解决方案，可能会采用“数值计算”中实现的“教科书算法”。 - user2961927