r - RcppArmadillo 和 Armadillo 之间的性能差异

Question

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

使用 X 的 1000000x100 矩阵运行大约需要 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)。

Answer 1

快速猜测:您的 native 程序使用加速 BLAS，而您的 R 构建则没有。

实际的“矩阵数学”是由 Armadillo 外包给 BLAS 库的。借助 RcppArmadillo，您可以获得 R 的构建目标。对于 native 程序，也许您会使用其他程序。这可能很简单，就像您的程序开始使用 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。