RcppArmadillo与Armadillo之间的性能差异

Question

RcppArmadillo与Armadillo之间的性能差异

6

我正在试图理解使用RcppArmadillo编写的函数和使用Armadillo库编写的独立C++程序之间性能差异。例如，考虑以下简单函数，它使用传统教科书公式计算线性模型的系数。

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

using namespace Rcpp;
using namespace arma;

// [[Rcpp::export]]
void simpleLm(NumericMatrix Xr, NumericMatrix yr) {
   int n = Xr.nrow(), k = Xr.ncol();
   mat X(Xr.begin(), n, k, false);
   colvec y(yr.begin(), yr.nrow(), false);

   colvec coef = inv(X.t()*X)*X.t()*y;
}

这段代码在一个 1000000x100 的矩阵 X 上运行大约需要6秒钟。代码中的一些时间测量（未显示）表明所有的时间都花费在了 coef 计算上。

X <- matrix(rnorm(1000000*100), ncol=100)
y <- matrix(rep(1, 1000000))
system.time(simpleLm(X,y))

 user  system elapsed 
  6.028   0.009   6.040

现在考虑一个非常相似的C++函数，然后使用g++编译。

#include <iostream>
#include <armadillo>
#include <chrono>
#include <cstdlib>

using namespace std;
using namespace arma;

int main(int argc, char **argv) {
    int n = 1000000;
    mat X = randu<mat>(n,100);
    vec y = ones<vec>(n);

    chrono::steady_clock::time_point start = chrono::steady_clock::now();

    colvec coef = inv(X.t()*X)*X.t()*y;

    chrono::steady_clock::time_point end = chrono::steady_clock::now();

    chrono::duration<double, milli> diff = end - start;

    cout << diff.count() << endl;

    return 0;
}

在这里，coef 变量的计算只需要大约 0.5 秒钟，或者说只有使用 RcppArmadillo 时所需时间的 1/12。我使用的是 Mac OS X 10.9.2，R 3.1.0，Rcpp 0.11.1 和 RcppArmadillo 0.4.200.0。我使用 sourceCpp 函数编译了 Rcpp 示例。独立的 C++ 示例使用 Armadillo 4.200.0，并且我还使用 Homebrew 安装了 Mac 的 Fortran 编译器（brew install gfortran）。

- Abiel

你没有列出设置的优化标志：如果我没记错，R（因此也包括sourceCpp）默认为-O2，但你应该检查一下（在sourceCpp中尝试使用verbose=TRUE）。你应该确保编译独立的C++文件时也使用相同的优化级别。 - Kevin Ushey

是的 - R使用了在“configure；make；make install”运行时使用的任何内容，您可以通过“CXXFLAGS”和其他选项进行覆盖。优化不太可能导致Abiel在这里看到的数量级。 - Dirk Eddelbuettel

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Dirk Eddelbuettel · Accepted Answer

快速猜测：您的本地程序使用了加速的BLAS库，而您的R构建未使用。实际上，“矩阵数学”是由Armadillo委托给BLAS库完成的。使用RcppArmadillo时，您将获得R所构建的内容。对于本地程序，您可能使用其他库。这可能仅仅是因为您的程序可以使用Accelerate库，而R不能-我不太清楚，因为我不使用OS X。

但是，为了证明，在我的（i7，Linux）机器上，时间几乎相同。

首先，您的程序没有修改：

edd@max:/tmp$ g++ -std=c++11 -O3 -o abiel abiel.cpp -larmadillo -llapack
edd@max:/tmp$ ./abiel 
2454
edd@max:/tmp$

其次，您的程序需要包装成R可以调用的形式（见下文）：

R> library(Rcpp)
R> sourceCpp("/tmp/abielviaR.cpp")
R> abielDemo()
2354.41
[1] TRUE
R>

关于相同的内容。

abielviaR.cpp的代码如下。

#include <RcppArmadillo.h>
#include <chrono>

using namespace std;
using namespace arma;

// [[Rcpp::plugins(cpp11)]]
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
bool abielDemo() {
    int n = 1000000;
    mat X = randu<mat>(n,100);
    vec y = ones<vec>(n);

    chrono::steady_clock::time_point start = chrono::steady_clock::now();
    colvec coef = inv(X.t()*X)*X.t()*y;
    chrono::steady_clock::time_point end = chrono::steady_clock::now();
    chrono::duration<double, milli> diff = end - start;
    Rcpp::Rcout << diff.count() << endl;

    return true;
}

PS：你真的不应该使用（X'X)^(-1) X计算OLS。