Diferença de desempenho entre RcppArmadillo e Armadillo
Pergunta
Estou tentando entender a diferença de desempenho entre uma função escrita em RcppArmadillo e outra escrita em um programa C++ independente usando a biblioteca Armadillo.Por exemplo, considere a seguinte função simples que calcula os coeficientes de um modelo linear usando a fórmula tradicional dos livros didáticos.
// [[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;
}
Isso leva cerca de 6 segundos para ser executado com um 1000000x100
matriz para X
.Alguns horários no código (não mostrados) indicaram que todo o tempo é gasto no coef
Cálculo.
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
Agora considere uma função muito semelhante escrita em C++ que é então compilada com 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;
}
Aqui o cálculo do coef
variável leva apenas cerca de 0,5 segundos, ou apenas 1/12 do tempo como quando feito com RcppArmadillo.
Estou usando o Mac OS X 10.9.2 com R 3.1.0, Rcpp 0.11.1 e RcppArmadillo 0.4.200.0.Compilei o exemplo Rcpp usando a função sourceCpp.O exemplo C++ autônomo usa Armadillo 4.200.0, e também instalei o compilador Fortran para Mac usando Homebrew (brew install gfortran
).
Solução
Adivinhação rápida:seu programa nativo usa BLAS acelerado, sua compilação R não.
A verdadeira "matemática matricial" é transferida por Armadillo para a biblioteca BLAS.Com o RcppArmadillo, você obtém aquilo contra o qual o R foi construído.Com um programa nativo, talvez você use outra coisa.Poderia ser tão simples quanto o seu programa usar as bibliotecas Accelerate, enquanto R não - eu realmente não sei, pois não uso o OS X.
Mas, para demonstrar, na minha máquina (i7, Linux), os tempos são quase idênticos.
Primeiro, seu programa, inalterado:
edd@max:/tmp$ g++ -std=c++11 -O3 -o abiel abiel.cpp -larmadillo -llapack
edd@max:/tmp$ ./abiel
2454
edd@max:/tmp$
Segundo, seu programa está encapsulado em algo que R pode chamar (veja abaixo):
R> library(Rcpp)
R> sourceCpp("/tmp/abielviaR.cpp")
R> abielDemo()
2354.41
[1] TRUE
R>
Aproximadamente o mesmo.
O código de abielviaR.cpp
segue.
#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 Você realmente não deveria calcular OLS via (X'X)^(-1) X.