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!
źródło
ObsjakoIntegerVectorusunięcie niektórych rzutów.beta) pozostają stałe w trakcie obserwacjiObs. Gdybyśmy mieli predyktory zmieniające się w czasie , konieczne byłoby niejawne obliczeniedenomdla każdego z nichObs, na podstawie wartości macierzy projektowejX. 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!for-loopco da ci największy zysk. Również w bardziej ogólnym przypadku sugerowałbym użyciemodel.matrix(...)do utworzenia macierzy dla danych wejściowych w swoich funkcjach.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żdegoi. Dlatego sensowne jest użycie adouble 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:
Dla mnie ta funkcja jest około dziesięć razy szybsza niż twoja funkcja R.
źródło
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:
v3 to odpowiedź Olivera na
rng=false. Wersja 4 zawiera sugestie 2 i 3.Funkcja:
źródło