Разница в производительности между 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.В примере standalone 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.