Leistungsunterschied zwischen RcppArmadillo und Gürteltier
Frage
Ich versuche, einen Leistungsunterschied zwischen einer in RcppArmadillo geschriebenen Funktion und einer in einem eigenständigen C ++ - Programm mit der Armadillo-Bibliothek geschriebenen Funktion zu verstehen.Betrachten Sie beispielsweise die folgende einfache Funktion, die die Koeffizienten für ein lineares Modell unter Verwendung der traditionellen Lehrbuchformel berechnet.
// [[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;
}
Dies dauert ungefähr 6 Sekunden, um mit einem zu laufen 1000000x100
matrix für X
.Einige Timings im Code (nicht gezeigt) zeigten an, dass die ganze Zeit für die ausgegeben wird coef
Berechnung.
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
Betrachten Sie nun eine sehr ähnliche Funktion, die in C ++ geschrieben ist und dann mit kompiliert wird 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;
}
Hier wird die Berechnung der coef
die Variable dauert nur etwa 0,5 Sekunden oder nur 1/12 der Zeit wie bei RcppArmadillo.
Ich verwende Mac OS X 10.9.2 mit R 3.1.0, Rcpp 0.11.1 und RcppArmadillo 0.4.200.0.Ich habe das Rcpp-Beispiel mit der Funktion sourceCpp kompiliert.Das eigenständige C ++ - Beispiel verwendet Armadillo 4.200.0, und ich habe auch den Fortran-Compiler für Mac mit Homebrew installiert (brew install gfortran
).
Lösung
Schnelle Vermutung:ihr natives Programm verwendet beschleunigtes BLAS, Ihr R-Build nicht.
Die eigentliche "Matrixmathematik" wird von Armadillo an die BLAS-Bibliothek weitergegeben.Mit RcppArmadillo bekommen Sie, worauf R aufgebaut ist.Mit einem nativen Programm verwenden Sie vielleicht etwas anderes.Es könnte so einfach sein, wie Ihr Programm die Beschleunigungsbibliotheken verwendet, während R dies nicht tut - ich weiß es nicht wirklich, da ich OS X nicht verwende.
Aber um zu demonstrieren, auf meinem (i7, Linux) Rechner sind die Zeiten nahezu identisch.
Zuerst Ihr Programm, unverändert:
edd@max:/tmp$ g++ -std=c++11 -O3 -o abiel abiel.cpp -larmadillo -llapack
edd@max:/tmp$ ./abiel
2454
edd@max:/tmp$
Zweitens ist Ihr Programm in etwas eingewickelt, das R aufrufen kann (siehe unten):
R> library(Rcpp)
R> sourceCpp("/tmp/abielviaR.cpp")
R> abielDemo()
2354.41
[1] TRUE
R>
Ungefähr gleich.
Der Kodex von abielviaR.cpp
folgen.
#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 Sie sollten OLS wirklich nicht über (X'X) ^ (-1) X berechnen.