RcppArmadillo 和 Armadillo 之间的性能差异

Posted 2023-02-16

【中文标题】RcppArmadillo 和 Armadillo 之间的性能差异

【英文标题】：Performance difference between RcppArmadillo and Armadillo

【发布时间】：2014-04-14 00:19:02

【问题描述】：

我试图了解用 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。

我正在使用带有 R 3.1.0、Rcpp 0.11.1 和 RcppArmadillo 0.4.200.0 的 Mac OS X 10.9.2。我使用 sourceCpp 函数编译了 Rcpp 示例。独立的 C++ 示例使用 Armadillo 4.200.0，我还使用 Homebrew (brew install gfortran) 为 Mac 安装了 Fortran 编译器。

【问题讨论】：

您没有列出优化标志集：如果我没记错的话，R（因此 sourceCpp）默认为 -O2，但您应该检查一下（在 sourceCpp 中尝试 verbose=TRUE）。您应该确保您正在编译具有相同优化级别的独立 C++ 文件。

是的--R 使用configure; make; make install 运行时使用的任何内容，您可以通过CXXFLAGS 和朋友覆盖它。优化不太可能导致 Abiel 在这里看到的数量级。

【参考方案1】：

快速猜测：您的本机程序使用加速 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。

【讨论】：

感谢 Dirk，您解决了问题。我changedR 使用 Apple 的 Accelerate 框架中包含的 BLAS，现在我的时间在两个版本之间匹配。