Je travaille actuellement sur une méthode bayésienne qui nécessite plusieurs étapes d'optimisation d'un modèle logit multinomial par itération. J'utilise optim () pour effectuer ces optimisations, et une fonction objective écrite en R. Un profilage a révélé qu'optim () est le principal goulot d'étranglement.
Après avoir fouillé, j'ai trouvé cette question dans laquelle ils suggèrent que le recodage de la fonction objectif Rcpppourrait accélérer le processus. J'ai suivi la suggestion et recodé ma fonction d'objectif avec Rcpp, mais cela a fini par être plus lent (environ deux fois plus lent!).
C'était ma première fois avec Rcpp(ou quoi que ce soit lié à C ++) et je n'ai pas pu trouver un moyen de vectoriser le code. Une idée de comment l'accélérer?
Tl; dr: L'implémentation actuelle de la fonction dans Rcpp n'est pas aussi rapide que R vectorisée; comment le rendre plus rapide?
Un exemple reproductible :
1) Définir les fonctions objectives dans Ret Rcpp: log-vraisemblance d'un modèle multinomial d'interception uniquement
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) Comparez leur efficacité:
## 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) Maintenant, appelez-les 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
J'ai été quelque peu surpris que l'implémentation vectorisée dans R soit plus rapide. La mise en œuvre d'une version plus efficace dans Rcpp (par exemple, avec RcppArmadillo?) Peut produire des gains? Est-ce une meilleure idée de tout recoder dans Rcpp en utilisant un optimiseur C ++?
PS: première publication sur Stackoverflow!
la source
Obscomme uneIntegerVectorsuppression de certains lancers.beta) restent constants au fil des observationsObs. Si nous avions des prédicteurs variant dans le temps, un calcul implicite dedenomchacunObsdeviendrait nécessaire, basé sur la valeur d'une matrice de conceptionX. Cela étant dit, j'implémente déjà vos suggestions sur le reste de mon code avec de très bons gains :). Merci @RalfStubner, @Oliver et @thc pour vos réponses très perspicaces! Passons maintenant à mon prochain goulot d'étranglement!for-loopvous donnera le plus grand gain. Dans le cas plus général, je suggère également d'utilisermodel.matrix(...)pour créer votre matrice pour l'entrée dans vos fonctions.Votre fonction C ++ peut être rendue plus rapide en utilisant les observations suivantes. Au moins, le premier pourrait également être utilisé avec votre fonction R:
La façon dont vous calculez
denom[i]est la même pour tousi. Il est donc logique d'utiliser adouble denomet de faire ce calcul une seule fois. J'élimine également la soustraction de ce terme commun à la fin.Vos observations sont en fait un vecteur entier du côté R, et vous les utilisez également comme entiers en C ++. Utiliser un
IntegerVectorpour commencer rend beaucoup de casting inutile.Vous pouvez également indexer un en
NumericVectorutilisant unIntegerVectoren C ++. Je ne sais pas si cela améliore les performances, mais cela rend le code un peu plus court.Quelques changements plus liés au style qu'à la performance.
Résultat:
Pour moi, cette fonction est environ dix fois plus rapide que votre fonction R.
la source
Je peux penser à quatre optimisations potentielles par rapport aux réponses de Ralf et Olivers.
(Vous devriez accepter leurs réponses, mais je voulais juste ajouter mes 2 cents).
1) Utiliser
// [[Rcpp::export(rng = false)]]comme en-tête de commentaire pour la fonction dans un fichier C ++ séparé. Cela entraîne une accélération de ~ 80% sur ma machine. (C'est la suggestion la plus importante des 4).2) Préférez
cmathsi possible. (Dans ce cas, cela ne semble pas faire de différence).3) Évitez autant que possible l'allocation, par exemple ne passez pas
betaà un nouveau vecteur.4) Objectif d'étirement: utilisez des
SEXPparamètres plutôt que des vecteurs Rcpp. (Laissé comme exercice au lecteur). Les vecteurs Rcpp sont des wrappers très fins, mais ils sont toujours des wrappers et il y a une petite surcharge.Ces suggestions ne seraient pas importantes, sinon pour le fait que vous appeliez la fonction en boucle serrée
optim. Donc, tout frais généraux est très important.Banc:
v3 est la réponse d'Oliver
rng=false. v4 est avec les suggestions # 2 et # 3 incluses.La fonction:
la source