Разница в производительности между RcppArmadillo и Armadillo

Question

Я пытаюсь понять разницу в производительности между функцией, написанной на RcppArmadillo, и одной написанной в автономной программе на C++ с использованием библиотеки Armadillo. Например, рассмотрим следующую простую функцию, которая вычисляет коэффициенты для линейной модели с использованием традиционной формулы учебника.

// [[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; 
} 

Это занимает около 6 секунд, чтобы работать с 1000000x100 матрицей для X. Некоторые тайминги кода (не показаны) показали, что все время расходуется на расчет 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 секунд, или только 1/12 время как, когда сделано с RcppArmadillo.

Я использую Mac OS X 10.9.2 с R 3.1.0, Rcpp 0.11.1 и RcppArmadillo 0.4.200.0. Я скомпилировал пример Rcpp, используя функцию sourceCpp. В отдельном примере на C++ используется Armadillo 4.200.0, и я также установил компилятор Fortran для Mac с помощью Homebrew (brew install gfortran).

Answer 1

Быстрое предположение: ваша родная программа использует ускоренный BLAS, вы не создали R build.

Фактическая «матричная математика» обрабатывается Armadillo в библиотеке BLAS. С RcppArmadillo вы получаете то, к чему строится R. С родной программой, возможно, вы используете что-то еще. Это может быть так же просто, как ваша программа, чтобы использовать библиотеки ускорения, в то время как 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 Вы действительно не должны вычислять OLS через (X'X)^(- 1) X, хотя.