Optimierung der R-Zielfunktion mit Rcpp langsamer, warum?

16

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!

smildiner
quelle

Antworten:

9

Wenn Sie in der Lage sind, vektorisierte Funktionen zu verwenden, werden Sie feststellen, dass diese (fast) so schnell sind, wie wenn Sie Ihren Code direkt in Rcpp ausführen. Dies liegt daran, dass viele vektorisierte Funktionen in R (fast alle vektorisierten Funktionen in Base R) in C, Cpp oder Fortran geschrieben sind und als solche oft wenig zu gewinnen ist.

Das heißt, es gibt Verbesserungen sowohl in Ihrem Rals auch in Ihrem RcppCode. Die Optimierung erfolgt durch sorgfältiges Studieren des Codes und Entfernen unnötiger Schritte (Speicherzuweisung, Summen usw.).

Beginnen wir mit der RcppCodeoptimierung.

In Ihrem Fall besteht die Hauptoptimierung darin, unnötige Matrix- und Vektorberechnungen zu entfernen. Der Code ist im Wesentlichen

  1. Shift Beta
  2. Berechnen Sie das Protokoll der Summe von exp (Shift Beta) [log-sum-exp]
  3. Verwenden Sie Obs als Index für das verschobene Beta und summieren Sie alle Wahrscheinlichkeiten
  4. subtrahieren Sie die log-sum-exp

Mit dieser Beobachtung können wir Ihren Code auf 2 for-Schleifen reduzieren. Beachten Sie, dass dies sumeinfach eine andere for-Schleife ist (mehr oder weniger for(i = 0; i < max; i++){ sum += x }:), sodass das Vermeiden der Summen den Code weiter beschleunigen kann (in den meisten Situationen ist dies eine unnötige Optimierung!). Außerdem ist Ihre Eingabe Obsein ganzzahliger Vektor, und wir können den Code weiter optimieren, indem wir den IntegerVectorTyp verwenden, um zu vermeiden, dass die doubleElemente in integerWerte umgewandelt werden (Dank an Ralf Stubners Antwort).

cppFunction('double llmnl_int_C_v2(NumericVector beta, IntegerVector Obs, int n_cat)
 {

    int n_Obs = Obs.size();

    NumericVector betas = (beta.size()+1);
    //1: shift beta
    for (int i = 1; i < n_cat; i++) {
        betas[i] = beta[i-1];
    };
    //2: Calculate log sum only once:
    double expBetas_log_sum = log(sum(exp(betas)));
    // pre allocate sum
    double ll_sum = 0;

    //3: Use n_Obs, to avoid calling Xby.size() every time 
    for (int i = 0; i < n_Obs; i++) {
        ll_sum += betas(Obs[i] - 1.0) ;
    };
    //4: Use that we know denom is the same for all I:
    ll_sum = ll_sum - expBetas_log_sum * n_Obs;
    return ll_sum;
}')

Beachten Sie, dass ich einige Speicherzuordnungen entfernt und unnötige Berechnungen in der for-Schleife entfernt habe. Außerdem habe ich verwendet, dass denomdies für alle Iterationen gleich ist und einfach für das Endergebnis multipliziert wird.

Wir können ähnliche Optimierungen in Ihrem R-Code durchführen, was zu der folgenden Funktion führt:

llmnl_int_R_v2 <- function(beta, Obs, n_cat) {
    n_Obs <- length(Obs)
    betas <- c(0, beta)
    #note: denom = log(sum(exp(betas)))
    sum(betas[Obs]) - log(sum(exp(betas))) * n_Obs
}

Beachten Sie, dass die Komplexität der Funktion drastisch reduziert wurde, damit andere sie leichter lesen können. Um sicherzugehen, dass ich den Code nicht irgendwo durcheinander gebracht habe, überprüfen wir, ob sie die gleichen Ergebnisse liefern:

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))

beta = c(4,2,1)
Obs = mnl_sample 
n_cat = 4
xr <- llmnl_int(beta = beta, Obs = mnl_sample, n_cat = n_cat)
xr2 <- llmnl_int_R_v2(beta = beta, Obs = mnl_sample, n_cat = n_cat)
xc <- llmnl_int_C(beta = beta, Obs = mnl_sample, n_cat = n_cat)
xc2 <- llmnl_int_C_v2(beta = beta, Obs = mnl_sample, n_cat = n_cat)
all.equal(c(xr, xr2), c(xc, xc2))
TRUE

Nun, das ist eine Erleichterung.

Performance:

Ich werde Microbenchmark verwenden, um die Leistung zu veranschaulichen. Die optimierten Funktionen sind schnell, daher werde ich die Funktionszeiten ausführen 1e5, um den Effekt des Garbage Collector zu reduzieren

microbenchmark("llmml_int_R" = llmnl_int(beta = beta, Obs = mnl_sample, n_cat = n_cat),
               "llmml_int_C" = llmnl_int_C(beta = beta, Obs = mnl_sample, n_cat = n_cat),
               "llmnl_int_R_v2" = llmnl_int_R_v2(beta = beta, Obs = mnl_sample, n_cat = n_cat),
               "llmml_int_C_v2" = llmnl_int_C_v2(beta = beta, Obs = mnl_sample, n_cat = n_cat),
               times = 1e5)
#Output:
#Unit: microseconds
#           expr     min      lq       mean  median      uq        max neval
#    llmml_int_R 202.701 206.801 288.219673 227.601 334.301  57368.902 1e+05
#    llmml_int_C 250.101 252.802 342.190342 272.001 399.251 112459.601 1e+05
# llmnl_int_R_v2   4.800   5.601   8.930027   6.401   9.702   5232.001 1e+05
# llmml_int_C_v2   5.100   5.801   8.834646   6.700  10.101   7154.901 1e+05

Hier sehen wir das gleiche Ergebnis wie zuvor. Jetzt sind die neuen Funktionen im Vergleich zu ihren ersten Gegenstücken ungefähr 35x schneller (R) und 40x schneller (Cpp). Interessanterweise ist die optimierte RFunktion immer noch geringfügig (0,3 ms oder 4%) schneller als meine optimierte CppFunktion. Meine beste Wette hier ist, dass das RcppPaket etwas Overhead hat , und wenn dies entfernt würde, wären die beiden identisch oder die R.

Ebenso können wir die Leistung mit Optim überprüfen.

microbenchmark("llmnl_int" = optim(beta, llmnl_int, Obs = mnl_sample, 
                                   n_cat = n_cat, method = "BFGS", hessian = F, 
                                   control = list(fnscale = -1)),
               "llmnl_int_C" = optim(beta, llmnl_int_C, Obs = mnl_sample, 
                                     n_cat = n_cat, method = "BFGS", hessian = F, 
                                     control = list(fnscale = -1)),
               "llmnl_int_R_v2" = optim(beta, llmnl_int_R_v2, Obs = mnl_sample, 
                                     n_cat = n_cat, method = "BFGS", hessian = F, 
                                     control = list(fnscale = -1)),
               "llmnl_int_C_v2" = optim(beta, llmnl_int_C_v2, Obs = mnl_sample, 
                                     n_cat = n_cat, method = "BFGS", hessian = F, 
                                     control = list(fnscale = -1)),
               times = 1e3)
#Output:
#Unit: microseconds
#           expr       min        lq      mean    median         uq      max neval
#      llmnl_int 29541.301 53156.801 70304.446 76753.851  83528.101 196415.5  1000
#    llmnl_int_C 36879.501 59981.901 83134.218 92419.551 100208.451 190099.1  1000
# llmnl_int_R_v2   667.802  1253.452  1962.875  1585.101   1984.151  22718.3  1000
# llmnl_int_C_v2   704.401  1248.200  1983.247  1671.151   2033.401  11540.3  1000

Das Ergebnis ist wieder dasselbe.

Fazit:

Als kurze Schlussfolgerung ist anzumerken, dass dies ein Beispiel ist, bei dem die Konvertierung Ihres Codes in Rcpp die Mühe nicht wirklich wert ist. Dies ist nicht immer der Fall, aber oft lohnt es sich, einen zweiten Blick auf Ihre Funktion zu werfen, um festzustellen, ob es Bereiche in Ihrem Code gibt, in denen unnötige Berechnungen durchgeführt werden. Insbesondere in Situationen, in denen vektorisierte Funktionen verwendet werden, lohnt es sich oft nicht, Code in Rcpp zu konvertieren. Häufiger kann man große Verbesserungen feststellen, wenn man for-loopsCode verwendet, der nicht einfach vektorisiert werden kann, um die for-Schleife zu entfernen.

Oliver
quelle
1
Sie können behandeln Obsals ein IntegerVectorpaar Würfe zu entfernen.
Ralf Stubner
Ich habe es nur aufgenommen, bevor ich Ihnen dafür gedankt habe, dass Sie dies in Ihrer Antwort bemerkt haben. Es ging einfach an mir vorbei. Ich habe Ihnen dies in meiner Antwort @RalfStubner gutgeschrieben. :-)
Oliver
2
Wie Sie an diesem Spielzeugbeispiel (Nur-Intercept-Mnl-Modell) bemerkt haben, betableiben die linearen Prädiktoren ( ) über die Beobachtungen konstant Obs. Wenn wir zeitvariable Prädiktoren hätten, wäre eine implizite Berechnung denomfür jeden Obserforderlich, basierend auf dem Wert einer Entwurfsmatrix X. 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!
Smildiner
1
Ich bin froh, dass wir helfen konnten. Im allgemeineren Fall wird die Berechnung des Subtraktions-Nennwerts in jedem Schritt der Sekunde berechnet, wodurch for-loopSie den größten Gewinn erzielen. Auch im allgemeineren Fall würde ich vorschlagen model.matrix(...), Ihre Matrix für die Eingabe in Ihre Funktionen zu erstellen.
Oliver
9

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 gleich i. Es ist daher sinnvoll, a zu verwenden double 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 einem IntegerVectorin 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:

double llmnl_int_C(NumericVector beta, IntegerVector 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];
    };

    double denom = log(sum(exp(betas)));
    NumericVector Xby = betas[Obs - 1];

    return sum(Xby) - n_Obs * denom;
}

Für mich ist diese Funktion ungefähr zehnmal schneller als Ihre R-Funktion.

Ralf Stubner
quelle
Vielen Dank für Ihre Antwort Ralph, hat den Eingabetyp nicht erkannt. Ich habe dies auch in meine Antwort aufgenommen und Ihnen die Ehre gegeben. :-)
Oliver
7

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:

microbenchmark("llmnl_int_R_v1" = optim(beta, llmnl_int, Obs = mnl_sample, 
                                      n_cat = n_cat, method = "BFGS", hessian = F, 
                                      control = list(fnscale = -1)),
             "llmnl_int_R_v2" = optim(beta, llmnl_int_R_v2, Obs = mnl_sample, 
                                      n_cat = n_cat, method = "BFGS", hessian = F, 
                                      control = list(fnscale = -1)),
             "llmnl_int_C_v2" = optim(beta, llmnl_int_C_v2, Obs = mnl_sample, 
                                      n_cat = n_cat, method = "BFGS", hessian = F, 
                                      control = list(fnscale = -1)),
             "llmnl_int_C_v3" = optim(beta, llmnl_int_C_v3, Obs = mnl_sample, 
                                      n_cat = n_cat, method = "BFGS", hessian = F, 
                                      control = list(fnscale = -1)),
             "llmnl_int_C_v4" = optim(beta, llmnl_int_C_v4, Obs = mnl_sample, 
                                      n_cat = n_cat, method = "BFGS", hessian = F, 
                                      control = list(fnscale = -1)),
             times = 1000)


Unit: microseconds
expr      min         lq       mean     median         uq        max neval cld
llmnl_int_R_v1 9480.780 10662.3530 14126.6399 11359.8460 18505.6280 146823.430  1000   c
llmnl_int_R_v2  697.276   735.7735  1015.8217   768.5735   810.6235  11095.924  1000  b 
llmnl_int_C_v2  997.828  1021.4720  1106.0968  1031.7905  1078.2835  11222.803  1000  b 
llmnl_int_C_v3  284.519   295.7825   328.5890   304.0325   328.2015   9647.417  1000 a  
llmnl_int_C_v4  245.650   256.9760   283.9071   266.3985   299.2090   1156.448  1000 a 

v3 ist Olivers Antwort mit rng=false. v4 enthält die Vorschläge 2 und 3.

Die Funktion:

#include <Rcpp.h>
#include <cmath>
using namespace Rcpp;

// [[Rcpp::export(rng = false)]]
double llmnl_int_C_v4(NumericVector beta, IntegerVector Obs, int n_cat) {

  int n_Obs = Obs.size();
  //2: Calculate log sum only once:
  // double expBetas_log_sum = log(sum(exp(betas)));
  double expBetas_log_sum = 1.0; // std::exp(0)
  for (int i = 1; i < n_cat; i++) {
    expBetas_log_sum += std::exp(beta[i-1]);
  };
  expBetas_log_sum = std::log(expBetas_log_sum);

  double ll_sum = 0;
  //3: Use n_Obs, to avoid calling Xby.size() every time 
  for (int i = 0; i < n_Obs; i++) {
    if(Obs[i] == 1L) continue;
    ll_sum += beta[Obs[i]-2L];
  };
  //4: Use that we know denom is the same for all I:
  ll_sum = ll_sum - expBetas_log_sum * n_Obs;
  return ll_sum;
}
thc
quelle