如何使用 Rcpp 将 lgamma 应用于矩阵（速度会更快）？

Question

我想知道是否可以使用 Rcpp 将 lgamma 应用于大型矩阵的所有条目。我尝试使用矢量：

// lgammaRcpp.cpp
#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
NumericVector lgammaRcpp(NumericVector v){
    NumericVector out;
    out = lgamma(v);
    return(out);
}

我做了一个简单的微基准测试：

library("microbenchmark")
x <- round(runif(100000)+50000);
microbenchmark(
   lgammaRcpp(x),
   lgamma(x)
)

并且 Rcpp 稍快：

Unit: milliseconds
          expr      min       lq     mean   median       uq      max neval
 lgammaRcpp(x) 5.405556 5.416283 5.810254 5.436139 5.511993 8.650419   100    
     lgamma(x) 5.613717 5.628769 6.114942 5.644215 6.872677 9.947497   100

当我尝试使用 "NumericMatrix" 时，但是：

// [[Rcpp::export]]
NumericMatrix lgammaRcpp(NumericMatrix v){
    NumericMatrix out;
    out = lgamma(v);
    return(out);
}

有些错误我不明白，例如

/home/canghel/R/x86_64-pc-linux-gnu-library/3.4/Rcpp/include/Rcpp/vector   /Matrix.h:83:13: note: Rcpp::Matrix<RTYPE, StoragePolicy>& Rcpp::Matrix<RTYPE, StoragePolicy>::operator=(const Rcpp::Matrix<RTYPE, StoragePolicy>&) [with int RTYPE = 14; StoragePolicy = Rcpp::PreserveStorage]
 Matrix& operator=(const Matrix& other) {

我的问题是：1) 有没有办法修改我的函数以将 lgamma 应用于矩阵的所有条目？和 2) 是否值得，或者为 lgamma 函数调用的底层库是否与 C++ 和 R 相同？

Answer 1

Rcpp Sugar 倾向于 return 向量，除非另有说明。因此，在这种情况下，您将总是返回类型为 Numeric 的 Vector，例如NumericVector。在此处查看我关于不同糖功能的注释：https://github.com/coatless/rcpp-api

以下允许根据上述注释进行编译：

#include <Rcpp.h>

// [[Rcpp::export]]
NumericVector lgammaRcpp(NumericMatrix v) { 
    NumericVector out;
    out = lgamma(v);
    return(out);
}

由于使用的函数相同，因此您不太可能看到大幅加速。您的上述基准测试部分表明了这一点，可以通过查看 Rcpp Math defines 来验证。 现在，这并不是说没有好处。 特别是，这里的主要好处是如果您将例程完全封装在 C++ 中。在这种情况下，与从 C++.

R

Answer 2

使用 Rfast 包将 lgamma/digamma 等函数应用于矩阵似乎更好（即更快）。

library("microbenchmark");
library("RcppArmadillo");
library("Rfast");
sourceCpp("lgammaRcpp.cpp");

x <- matrix(round(runif(100000)+50000), 100, 1000);
microbenchmark(
    lgammaRcpp(x),
    lgamma(x),
    Rfast::Lgamma(x)
)
Unit: milliseconds
          expr      min       lq     mean   median       uq      max neval
lgammaRcppArma(x) 4.654526 4.919831 5.577843 5.413790 5.888895 9.258325   100
lgamma(x) 5.572671 5.840268 6.582007 6.131651 7.280895 8.779301   100
Rfast::Lgamma(x) 4.450824 4.588596 5.128323 4.791287 5.608678 6.865331   100

我有：

#include<RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]

// [[Rcpp::export]]
arma::mat lgammaRcpp(arma::mat m) {
    arma::mat out = lgamma(m);
    return(out);
}

如何使用 Rcpp 将 lgamma 应用于矩阵（速度会更快）？

How to apply lgamma to a matrix using Rcpp (and will it be faster)?

r

matrix

vectorization

rcpp