Optimiser la fonction objectif R avec Rcpp plus lentement, pourquoi?

16

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!

smildiner
la source

Réponses:

9

En général, si vous pouvez utiliser des fonctions vectorisées, vous constaterez qu'il est (presque) aussi rapide que l'exécution de votre code directement dans Rcpp. En effet, de nombreuses fonctions vectorisées dans R (presque toutes les fonctions vectorisées dans Base R) sont écrites en C, Cpp ou Fortran et en tant que telles, il y a souvent peu à gagner.

Cela dit, il y a des améliorations à gagner à la fois dans votre Ret dans le Rcppcode. L'optimisation provient de l'étude attentive du code et de la suppression des étapes inutiles (affectation de la mémoire, sommes, etc.).

Commençons par l' Rcppoptimisation du code.

Dans votre cas, l'optimisation principale consiste à supprimer les calculs de matrice et de vecteur inutiles. Le code est par essence

  1. Shift beta
  2. calculer le log de la somme de exp (shift beta) [log-sum-exp]
  3. utiliser Obs comme indice pour la bêta décalée et additionner toutes les probabilités
  4. soustraire le log-sum-exp

En utilisant cette observation, nous pouvons réduire votre code à 2 boucles pour. Notez qu'il sums'agit simplement d'une autre boucle for (plus ou moins for(i = 0; i < max; i++){ sum += x }:), donc éviter les sommes peut accélérer davantage le code (dans la plupart des cas, c'est une optimisation inutile!). De plus, votre entrée Obsest un vecteur entier, et nous pouvons optimiser davantage le code en utilisant le IntegerVectortype pour éviter de convertir les doubleéléments en integervaleurs (crédit à la réponse de Ralf Stubner).

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

Notez que j'ai supprimé un certain nombre d'allocations de mémoire et supprimé les calculs inutiles dans la boucle for. J'ai également utilisé denomle même pour toutes les itérations et simplement multiplié pour le résultat final.

Nous pouvons effectuer des optimisations similaires dans votre code R, ce qui se traduit par la fonction ci-dessous:

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
}

Notez que la complexité de la fonction a été considérablement réduite, ce qui la rend plus facile à lire pour les autres. Juste pour être sûr que je n'ai pas foiré le code quelque part, vérifions qu'ils retournent les mêmes résultats:

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

c'est un soulagement.

Performance:

Je vais utiliser un micro-benchmark pour illustrer la performance. Les fonctions optimisées sont rapides, donc je vais exécuter les fonctions 1e5fois pour réduire l'effet du garbage collector

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

Ici, nous voyons le même résultat qu'auparavant. Désormais, les nouvelles fonctions sont environ 35 fois plus rapides (R) et 40 fois plus rapides (Cpp) par rapport à leurs premières contreparties. Chose intéressante, la Rfonction optimisée est encore très légèrement (0,3 ms ou 4%) plus rapide que ma Cppfonction optimisée . Mon meilleur pari ici est qu'il y a des frais généraux dans le Rcpppackage, et si cela était supprimé, les deux seraient identiques ou le R.

De même, nous pouvons vérifier les performances à l'aide d'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

Encore une fois, le résultat est le même.

Conclusion:

En conclusion, il convient de noter qu'il s'agit d'un exemple, où la conversion de votre code en Rcpp ne vaut pas vraiment la peine. Ce n'est pas toujours le cas, mais il vaut souvent la peine de revoir votre fonction, pour voir s'il y a des zones de votre code, où des calculs inutiles sont effectués. Surtout dans les situations où l'on utilise des fonctions vectorisées intégrées, cela ne vaut souvent pas la peine de convertir du code en Rcpp. Plus souvent, on peut voir de grandes améliorations si l'on utilise for-loopsavec du code qui ne peut pas être facilement vectorisé afin de supprimer la boucle for.

Oliver
la source
1
Vous pouvez traiter Obscomme une IntegerVectorsuppression de certains lancers.
Ralf Stubner
Je viens de l'intégrer avant de vous remercier de l'avoir remarqué dans votre réponse. Il est simplement passé par moi. Je vous en donne le mérite dans ma réponse @RalfStubner. :-)
Oliver
2
Comme vous l'avez remarqué sur cet exemple de jouet (modèle mnl d'interception uniquement), les prédicteurs linéaires ( beta) restent constants au fil des observations Obs. Si nous avions des prédicteurs variant dans le temps, un calcul implicite de denomchacun Obsdeviendrait nécessaire, basé sur la valeur d'une matrice de conception X. 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!
smildiner
1
Je suis content que nous puissions vous aider. Dans le cas plus général, le calcul du soustrait du dénomination à chaque étape de la seconde for-loopvous donnera le plus grand gain. Dans le cas plus général, je suggère également d'utiliser model.matrix(...)pour créer votre matrice pour l'entrée dans vos fonctions.
Oliver
9

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 tous i. Il est donc logique d'utiliser a double 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 un IntegerVectoren 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:

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

Pour moi, cette fonction est environ dix fois plus rapide que votre fonction R.

Ralf Stubner
la source
Merci pour votre réponse Ralph, je n'ai pas repéré le type d'entrée. J'ai intégré cela dans ma réponse et je vous en remercie également. :-)
Oliver
7

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:

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 est la réponse d'Oliver rng=false. v4 est avec les suggestions # 2 et # 3 incluses.

La fonction:

#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
la source