Optymalizacja funkcji celu R z wolniejszym Rcpp, dlaczego?


16

Obecnie pracuję nad metodą bayesowską, która wymaga wielu kroków optymalizacji wielomianowego modelu logit dla każdej iteracji. Korzystam z Optim () do przeprowadzania tych optymalizacji, a funkcja celu napisana w R. Profilowanie wykazało, że Optim () jest głównym wąskim gardłem.

Po przekopaniu znalazłem pytanie, w którym sugerują, że przekodowanie funkcji celu Rcppmoże przyspieszyć proces. Postępowałem zgodnie z sugestią i przekodowałem swoją funkcję celu Rcpp, ale okazało się, że jest wolniejsza (około dwa razy wolniejsza!).

To był mój pierwszy raz Rcpp(lub cokolwiek związanego z C ++) i nie byłem w stanie znaleźć sposobu na wektoryzację kodu. Masz pomysł, jak to zrobić szybciej?

Tl; dr: Bieżąca implementacja funkcji w Rcpp nie jest tak szybka jak wektoryzacja R; jak zrobić to szybciej?

Odtwarzalny przykład :

1) Zdefiniuj funkcje celu w : Ri Rcpplog-prawdopodobieństwo prawdopodobieństwa modelu wielomianowego tylko do przechwytywania

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) Porównaj ich skuteczność:

## 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) Teraz dzwoniąc do nich 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

Byłem nieco zaskoczony, że wektoryzacja implementacji w języku R była szybsza. Wdrożenie bardziej wydajnej wersji w Rcpp (powiedzmy, z RcppArmadillo?) Może przynieść jakieś korzyści? Czy lepiej jest przekodować wszystko w Rcpp przy użyciu optymalizatora C ++?

PS: pierwsze publikowanie w Stackoverflow!

Odpowiedzi:


9

Ogólnie, jeśli jesteś w stanie używać funkcji wektoryzowanych, przekonasz się, że jest (prawie) tak szybki, jak uruchomienie kodu bezpośrednio w Rcpp. Wynika to z faktu, że wiele funkcji wektoryzowanych w R (prawie wszystkie funkcje wektoryzowane w Bazie R) są zapisane w C, Cpp lub Fortran i jako takie często niewiele można zyskać.

To powiedziawszy, istnieją ulepszenia do zdobycia zarówno w twoim, jak Ri Rcppkodzie. Optymalizacja wynika z dokładnego przestudiowania kodu i usunięcia niepotrzebnych kroków (przypisanie pamięci, sumy itp.).

Zacznijmy od Rcppoptymalizacji kodu.

W twoim przypadku główna optymalizacja polega na usunięciu niepotrzebnych obliczeń macierzy i wektorów. Kod jest w istocie

  1. Shift beta
  2. obliczyć log sumy exp (przesunięcie beta) [log-sum-exp]
  3. użyj Obs jako wskaźnika przesuniętej wersji beta i zsumuj wszystkie prawdopodobieństwa
  4. odejmij log-sum-exp

Korzystając z tej obserwacji, możemy zredukować twój kod do 2 pętli for. Zauważ, że sumjest to po prostu kolejna pętla for (mniej więcej:), for(i = 0; i < max; i++){ sum += x }więc unikanie sum może przyspieszyć kod (w większości przypadków jest to niepotrzebna optymalizacja!). Ponadto dane wejściowe Obssą wektorem liczb całkowitych, a my możemy dalej optymalizować kod za pomocą tego IntegerVectortypu, aby uniknąć rzucania doubleelementów na integerwartości (Podziękowania dla odpowiedzi Ralfa Stubnera).

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;
}')

Zauważ, że usunąłem sporo przydziałów pamięci i usunąłem niepotrzebne obliczenia w pętli for. Również użyłem tego denomsamego dla wszystkich iteracji i po prostu pomnożyłem dla końcowego wyniku.

Możemy przeprowadzić podobne optymalizacje w twoim kodzie R, co skutkuje poniższą funkcją:

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
}

Zwróć uwagę, że złożoność funkcji została drastycznie zmniejszona, dzięki czemu inni łatwiej ją czytają. Żeby się upewnić, że gdzieś w kodzie nie pomyliłem się, sprawdźmy, czy zwracają te same wyniki:

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

to ulga.

Wydajność:

Użyję mikrodruku do zilustrowania wydajności. Zoptymalizowane funkcje są szybkie, więc uruchomię te funkcje 1e5razy, aby zmniejszyć efekt śmieciarza

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

Tutaj widzimy taki sam wynik jak poprzednio. Teraz nowe funkcje są około 35 razy szybsze (R) i 40 razy szybsze (Cpp) w porównaniu do ich pierwszych przeciw-części. Co ciekawe, zoptymalizowana Rfunkcja jest nadal bardzo nieznacznie (0,3 ms lub 4%) szybsza niż moja zoptymalizowana Cppfunkcja. Moim najlepszym zakładem jest to, że z Rcpppakietu jest trochę narzutu , a jeśli to zostanie usunięte, dwa będą identyczne lub R.

Podobnie możemy sprawdzić wydajność za pomocą Optim.

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

Po raz kolejny wynik jest taki sam.

Wniosek:

Podsumowując, warto zauważyć, że jest to jeden przykład, w którym konwersja kodu do Rcpp nie jest tak naprawdę warta kłopotów. Nie zawsze tak jest, ale często warto ponownie przyjrzeć się funkcji, aby sprawdzić, czy istnieją obszary kodu, w których wykonywane są niepotrzebne obliczenia. Zwłaszcza w sytuacjach, w których używa się wbudowanych funkcji wektoryzowanych, często nie warto poświęcać czasu na konwersję kodu do Rcpp. Częściej widać wielkie ulepszenia, jeśli używa for-loopssię kodu, którego nie można łatwo wektoryzować w celu usunięcia pętli for.


1
Możesz traktować Obsjako IntegerVectorusunięcie niektórych rzutów.
Ralf Stubner

Właśnie to uwzględniłem, zanim podziękowałem za zauważenie tego w swojej odpowiedzi. Po prostu mi minęło. Uwierzyłem ci to w mojej odpowiedzi @RalfStubner. :-)
Oliver

2
Jak zauważyłeś na tym przykładzie zabawki (model mnlowy tylko do przechwytywania), predyktory liniowe ( beta) pozostają stałe w trakcie obserwacji Obs. Gdybyśmy mieli predyktory zmieniające się w czasie , konieczne byłoby niejawne obliczenie denomdla każdego z nich Obs, na podstawie wartości macierzy projektowej X. Biorąc to pod uwagę, wdrażam już twoje sugestie w pozostałej części mojego kodu z kilkoma naprawdę fajnymi korzyściami :). Dziękuję @RalfStubner, @Oliver i @thc za bardzo wnikliwe odpowiedzi! Teraz przechodzę do następnego wąskiego gardła!
smildiner

1
Cieszę się, że mogliśmy pomóc. W bardziej ogólnym przypadku obliczanie odejmowania denomu na każdym etapie drugiego kroku, for-loopco da ci największy zysk. Również w bardziej ogólnym przypadku sugerowałbym użycie model.matrix(...)do utworzenia macierzy dla danych wejściowych w swoich funkcjach.
Oliver

9

Funkcję C ++ można przyspieszyć, korzystając z następujących obserwacji. Przynajmniej pierwszy może być również używany z funkcją R:

  • Sposób obliczania denom[i]jest taki sam dla każdego i. Dlatego sensowne jest użycie a double denomi wykonanie tego obliczenia tylko raz. Na koniec odejmuję też odejmowanie tego wspólnego terminu.

  • Twoje obserwacje są w rzeczywistości wektorem liczb całkowitych po stronie R i używasz ich również jako liczb całkowitych w C ++. Używanie IntegerVectorna początku sprawia, że ​​wiele rzutów nie jest konieczne.

  • Możesz indeksować za NumericVectorpomocą również IntegerVectorw C ++. Nie jestem pewien, czy to poprawia wydajność, ale sprawia, że ​​kod jest nieco krótszy.

  • Więcej zmian, które są bardziej związane ze stylem niż z wydajnością.

Wynik:

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;
}

Dla mnie ta funkcja jest około dziesięć razy szybsza niż twoja funkcja R.


Dzięki za odpowiedź Ralph, nie zauważył typu danych wejściowych. Uwzględniłem to również w mojej odpowiedzi, dając ci kredyt. :-)
Oliver

7

Mogę wymyślić cztery potencjalne optymalizacje w stosunku do odpowiedzi Ralfa i Oliversa.

(Powinieneś zaakceptować ich odpowiedzi, ale chciałem tylko dodać 2 centy).

1) Użyj // [[Rcpp::export(rng = false)]]jako nagłówka komentarza do funkcji w osobnym pliku C ++. Prowadzi to do ~ 80% przyspieszenia na mojej maszynie. (To najważniejsza propozycja spośród 4).

2) Preferuj, cmathjeśli to możliwe. (W tym przypadku nie ma to znaczenia).

3) Unikaj alokacji, gdy to możliwe, np. Nie przechodź betana nowy wektor.

4) Cel SEXProzciągnięcia : użyj parametrów zamiast wektorów Rcpp. (Pozostawione jako ćwiczenie dla czytelnika). Wektory Rcpp są bardzo cienkimi opakowaniami, ale nadal są opakowaniami i istnieje niewielki narzut.

Te sugestie nie byłyby ważne, gdyby nie fakt, że wywołujesz funkcję w ciasnej pętli optim. Tak więc wszelkie koszty ogólne są bardzo ważne.

Ławka:

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 to odpowiedź Olivera na rng=false. Wersja 4 zawiera sugestie 2 i 3.

Funkcja:

#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;
}
Korzystając z naszej strony potwierdzasz, że przeczytałeś(-aś) i rozumiesz nasze zasady używania plików cookie i zasady ochrony prywatności.
Licensed under cc by-sa 3.0 with attribution required.