Différence de performances entre RcppArmadillo et Armadillo
Question
J'essaie de comprendre une différence de performances entre une fonction écrite en RcppArmadillo et une fonction écrite dans un programme C++ autonome utilisant la bibliothèque Armadillo.Par exemple, considérons la fonction simple suivante qui calcule les coefficients d'un modèle linéaire à l'aide de la formule traditionnelle du manuel.
// [[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;
}
Cela prend environ 6 secondes pour s'exécuter avec un 1000000x100
matrice pour X
.Certains timings dans le code (non représentés) indiquaient que tout le temps était passé sur le coef
calcul.
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
Considérons maintenant une fonction très similaire écrite en C++ qui est ensuite compilée avec 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;
}
Ici le calcul du coef
La variable ne prend qu'environ 0,5 seconde, soit seulement 1/12 du temps comme avec RcppArmadillo.
J'utilise Mac OS X 10.9.2 avec R 3.1.0, Rcpp 0.11.1 et RcppArmadillo 0.4.200.0.J'ai compilé l'exemple Rcpp en utilisant la fonction sourceCpp.L'exemple C++ autonome utilise Armadillo 4.200.0 et j'ai également installé le compilateur Fortran pour Mac à l'aide de Homebrew (brew install gfortran
).
La solution
Devinez rapidement :votre programme natif utilise BLAS accéléré, votre build R ne le fait pas.
Les "mathématiques matricielles" actuelles sont confiées par Armadillo à la bibliothèque BLAS.Avec RcppArmadillo, vous obtenez ce contre quoi R est construit.Avec un programme natif, vous utilisez peut-être autre chose.Cela pourrait être aussi simple que votre programme utilise les bibliothèques Accelerate alors que R ne le fait pas - je ne sais pas vraiment car je n'utilise pas OS X.
Mais pour démontrer, sur ma machine (i7, Linux), les temps sont presque identiques.
Tout d'abord, votre programme, inchangé :
edd@max:/tmp$ g++ -std=c++11 -O3 -o abiel abiel.cpp -larmadillo -llapack
edd@max:/tmp$ ./abiel
2454
edd@max:/tmp$
Deuxièmement, votre programme enveloppé dans quelque chose que R peut appeler (voir ci-dessous) :
R> library(Rcpp)
R> sourceCpp("/tmp/abielviaR.cpp")
R> abielDemo()
2354.41
[1] TRUE
R>
À peu près pareil.
Le code de abielviaR.cpp
suit.
#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 Vous ne devriez vraiment pas calculer OLS via (X'X)^(-1) X.