Ich arbeite derzeit an einer Bayes'schen Methode, die mehrere Schritte zur Optimierung eines multinomialen Logit-Modells pro Iteration erfordert. Ich verwende optim (), um diese Optimierungen durchzuführen, und eine in R geschriebene Zielfunktion. Eine Profilerstellung ergab, dass optim () der Hauptengpass ist.
Nachdem ich mich umgesehen hatte, fand ich diese Frage, in der vorgeschlagen wurde, dass das Neukodieren der Zielfunktion mit Rcppden Prozess beschleunigen könnte. Ich folgte dem Vorschlag und kodierte meine Zielfunktion mit Rcpp, aber es wurde langsamer (ungefähr zweimal langsamer!).
Dies war mein erstes Mal mit Rcpp(oder irgendetwas im Zusammenhang mit C ++) und ich konnte keinen Weg finden, den Code zu vektorisieren. Irgendeine Idee, wie man es schneller macht?
Tl; dr: Die derzeitige Implementierung der Funktion in Rcpp ist nicht so schnell wie die vektorisierte R; wie kann man es schneller machen?
Ein reproduzierbares Beispiel :
1) Definieren Sie Zielfunktionen in Rund Rcpp: Log-Wahrscheinlichkeit eines Intercept-Only-Multinomial-Modells
library(Rcpp)
library(microbenchmark)
llmnl_int <- function(beta, Obs, n_cat) {
n_Obs <- length(Obs)
Xint <- matrix(c(0, beta), byrow = T, ncol = n_cat, nrow = n_Obs)
ind <- cbind(c(1:n_Obs), Obs)
Xby <- Xint[ind]
Xint <- exp(Xint)
iota <- c(rep(1, (n_cat)))
denom <- log(Xint %*% iota)
return(sum(Xby - denom))
}
cppFunction('double llmnl_int_C(NumericVector beta, NumericVector Obs, int n_cat) {
int n_Obs = Obs.size();
NumericVector betas = (beta.size()+1);
for (int i = 1; i < n_cat; i++) {
betas[i] = beta[i-1];
};
NumericVector Xby = (n_Obs);
NumericMatrix Xint(n_Obs, n_cat);
NumericVector denom = (n_Obs);
for (int i = 0; i < Xby.size(); i++) {
Xint(i,_) = betas;
Xby[i] = Xint(i,Obs[i]-1.0);
Xint(i,_) = exp(Xint(i,_));
denom[i] = log(sum(Xint(i,_)));
};
return sum(Xby - denom);
}')
2) Vergleichen Sie ihre Effizienz:
## Draw sample from a multinomial distribution
set.seed(2020)
mnl_sample <- t(rmultinom(n = 1000,size = 1,prob = c(0.3, 0.4, 0.2, 0.1)))
mnl_sample <- apply(mnl_sample,1,function(r) which(r == 1))
## Benchmarking
microbenchmark("llmml_int" = llmnl_int(beta = c(4,2,1), Obs = mnl_sample, n_cat = 4),
"llmml_int_C" = llmnl_int_C(beta = c(4,2,1), Obs = mnl_sample, n_cat = 4),
times = 100)
## Results
# Unit: microseconds
# expr min lq mean median uq max neval
# llmnl_int 76.809 78.6615 81.9677 79.7485 82.8495 124.295 100
# llmnl_int_C 155.405 157.7790 161.7677 159.2200 161.5805 201.655 100
3) Rufen Sie sie jetzt an optim:
## Benchmarking with optim
microbenchmark("llmnl_int" = optim(c(4,2,1), llmnl_int, Obs = mnl_sample, n_cat = 4, method = "BFGS", hessian = T, control = list(fnscale = -1)),
"llmnl_int_C" = optim(c(4,2,1), llmnl_int_C, Obs = mnl_sample, n_cat = 4, method = "BFGS", hessian = T, control = list(fnscale = -1)),
times = 100)
## Results
# Unit: milliseconds
# expr min lq mean median uq max neval
# llmnl_int 12.49163 13.26338 15.74517 14.12413 18.35461 26.58235 100
# llmnl_int_C 25.57419 25.97413 28.05984 26.34231 30.44012 37.13442 100
Ich war etwas überrascht, dass die vektorisierte Implementierung in R schneller war. Die Implementierung einer effizienteren Version in Rcpp (z. B. mit RcppArmadillo?) Kann zu Gewinnen führen. Ist es eine bessere Idee, alles in Rcpp mit einem C ++ - Optimierer neu zu codieren?
PS: Erstes Posting bei Stackoverflow!
quelle
Obsals einIntegerVectorpaar Würfe zu entfernen.betableiben die linearen Prädiktoren ( ) über die Beobachtungen konstantObs. Wenn wir zeitvariable Prädiktoren hätten, wäre eine implizite Berechnungdenomfür jedenObserforderlich, basierend auf dem Wert einer EntwurfsmatrixX. Davon abgesehen setze ich Ihre Vorschläge für den Rest meines Codes bereits mit einigen wirklich schönen Gewinnen um :). Vielen Dank an @RalfStubner, @Oliver und @thc für Ihre sehr aufschlussreichen Antworten! Nun zu meinem nächsten Engpass!for-loopSie den größten Gewinn erzielen. Auch im allgemeineren Fall würde ich vorschlagenmodel.matrix(...), Ihre Matrix für die Eingabe in Ihre Funktionen zu erstellen.Ihre C ++ - Funktion kann mithilfe der folgenden Beobachtungen beschleunigt werden. Zumindest die erste könnte auch mit Ihrer R-Funktion verwendet werden:
Die Art und Weise, wie Sie berechnen,
denom[i]ist für alle gleichi. Es ist daher sinnvoll, a zu verwendendouble denomund diese Berechnung nur einmal durchzuführen. Am Ende ziehe ich auch das Subtrahieren dieses allgemeinen Begriffs heraus.Ihre Beobachtungen sind tatsächlich ein ganzzahliger Vektor auf der R-Seite, und Sie verwenden sie auch in C ++ als Ganzzahlen. Die Verwendung
IntegerVectorvon a macht zunächst viel Casting unnötig.Sie können a auch
NumericVectormit einemIntegerVectorin C ++ indizieren . Ich bin nicht sicher, ob dies die Leistung verbessert, aber es macht den Code etwas kürzer.Einige weitere Änderungen, die eher mit dem Stil als mit der Leistung zusammenhängen.
Ergebnis:
Für mich ist diese Funktion ungefähr zehnmal schneller als Ihre R-Funktion.
quelle
Ich kann mir vier mögliche Optimierungen für die Antworten von Ralf und Olivers vorstellen.
(Sie sollten ihre Antworten akzeptieren, aber ich wollte nur meine 2 Cent hinzufügen).
1) Verwenden Sie
// [[Rcpp::export(rng = false)]]als Kommentar-Header für die Funktion in einer separaten C ++ - Datei. Dies führt zu einer Beschleunigung von ca. 80% auf meinem Computer. (Dies ist der wichtigste Vorschlag unter den 4).2)
cmathWenn möglich bevorzugen . (In diesem Fall scheint es keinen Unterschied zu machen).3) Vermeiden Sie nach Möglichkeit die Zuordnung, z. B. wechseln Sie nicht
betain einen neuen Vektor.4) Stretch-Ziel: Verwenden Sie
SEXPParameter anstelle von Rcpp-Vektoren. (Als Übung dem Leser überlassen). Rcpp-Vektoren sind sehr dünne Wrapper, aber sie sind immer noch Wrapper und es gibt einen kleinen Overhead.Diese Vorschläge wären nicht wichtig, wenn Sie die Funktion nicht in einer engen Schleife aufrufen würden
optim. Jeder Overhead ist also sehr wichtig.Bank:
v3 ist Olivers Antwort mit
rng=false. v4 enthält die Vorschläge 2 und 3.Die Funktion:
quelle