Atualmente, estou trabalhando em um método bayesiano que requer várias etapas de otimização de um modelo de logit multinomial por iteração. Estou usando o optim () para realizar essas otimizações e uma função objetiva escrita em R. Uma criação de perfil revelou que o optim () é o principal gargalo.

Depois de pesquisar, encontrei esta pergunta na qual eles sugerem que a recodificação da função objetivo Rcpppoderia acelerar o processo. Eu segui a sugestão e recodifiquei minha função objetiva Rcpp, mas acabou sendo mais lenta (cerca de duas vezes mais lenta!).

Esta foi a minha primeira vez com Rcpp(ou qualquer coisa relacionada ao C ++) e não consegui encontrar uma maneira de vetorizar o código. Alguma idéia de como torná-lo mais rápido?

Tl; dr: A implementação atual da função no Rcpp não é tão rápida quanto o R vetorizado; como torná-lo mais rápido?

Um exemplo reproduzível :

1) Definir funções objetivas Re Rcpp: probabilidade logarítmica de um modelo multinomial de interceptação apenas

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) Compare sua eficiência:

## 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) Agora, chamando-os 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

Fiquei um pouco surpreso que a implementação vetorizada em R fosse mais rápida. Implementar uma versão mais eficiente no Rcpp (digamos, com o RcppArmadillo?) Pode gerar ganhos? É uma idéia melhor recodificar tudo no Rcpp usando um otimizador de C ++?

PS: primeira postagem no Stackoverflow!

Em geral, se você for capaz de usar funções vetorizadas, verá que é (quase) tão rápido quanto executar seu código diretamente no Rcpp. Isso ocorre porque muitas funções vetorizadas em R (quase todas as funções vetorizadas na Base R) são escritas em C, Cpp ou Fortran e, como tal, geralmente há pouco a ganhar.

Dito isto, existem melhorias a serem obtidas no código Re no seu Rcpp. A otimização vem do estudo cuidadoso do código e da remoção de etapas desnecessárias (atribuição de memória, somas etc.).

Vamos começar com a Rcppotimização do código.

No seu caso, a principal otimização é remover cálculos desnecessários de matriz e vetor. O código é essencialmente

1. Shift beta
2. calcular o log da soma de exp (shift beta) [log-sum-exp]
3. use Obs como um índice para o beta alternado e some todas as probabilidades
4. subtrair o log-sum-exp

Usando esta observação, podemos reduzir seu código para 2 for-loops. Observe que sumé simplesmente outro loop for (mais ou menos for(i = 0; i < max; i++){ sum += x }:), portanto, evitar as somas pode acelerar ainda mais o código (na maioria das situações, isso é otimização desnecessária!). Além disso, sua entrada Obsé um vetor inteiro, e podemos otimizar ainda mais o código usando o IntegerVectortipo para evitar converter os doubleelementos em integervalores (Crédito à resposta 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;
}')

Observe que removi algumas alocações de memória e removi cálculos desnecessários no loop for. Também usei denomo mesmo para todas as iterações e simplesmente multiplicado para o resultado final.

Podemos realizar otimizações semelhantes no seu código R, o que resulta na função abaixo:

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
}

Observe que a complexidade da função foi drasticamente reduzida, facilitando a leitura de outras pessoas. Apenas para ter certeza de que não estraguei o código em algum lugar, vamos verificar se eles retornam os mesmos resultados:

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

Bem, isso é um alívio.

Atuação:

Vou usar o microbenchmark para ilustrar o desempenho. As funções otimizadas são rápidas, portanto, executarei os 1e5tempos das funções para reduzir o efeito do coletor de lixo

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

Aqui vemos o mesmo resultado de antes. Agora, as novas funções são aproximadamente 35x mais rápidas (R) e 40x mais rápidas (Cpp) em comparação com suas primeiras contrapartes. Curiosamente, a Rfunção otimizada ainda é muito ligeiramente (0,3 ms ou 4%) mais rápida que a minha Cppfunção otimizada . Minha melhor aposta aqui é que há alguma sobrecarga no Rcpppacote e, se isso fosse removido, os dois seriam idênticos ou o R.

Da mesma forma, podemos verificar o desempenho usando o 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

Mais uma vez o resultado é o mesmo.

Conclusão:

Como uma conclusão curta, vale a pena notar que este é um exemplo, onde converter seu código em Rcpp realmente não vale a pena. Nem sempre é esse o caso, mas geralmente vale a pena examinar sua função de novo, para verificar se há áreas do seu código em que cálculos desnecessários são realizados. Especialmente em situações em que se usa funções vetorizadas de compilação, geralmente não vale a pena converter o código em Rcpp. Com mais frequência, pode-se ver grandes melhorias se for usado for-loopscom código que não pode ser facilmente vetorizado para remover o loop for.

Sua função C ++ pode ser mais rápida usando as seguintes observações. Pelo menos o primeiro também pode ser usado com sua função R:

• A maneira como você calcula denom[i] é a mesma para todos i. Portanto, faz sentido usar a double denome fazer esse cálculo apenas uma vez. Eu também fatorei a subtração desse termo comum no final.

• Suas observações são na verdade um vetor inteiro no lado R e você as está usando também como números inteiros em C ++. Usar um IntegerVectorpara começar faz com que muitas transmissões sejam desnecessárias.

• Você pode indexar um NumericVector usando umIntegerVector em C ++. Não tenho certeza se isso ajuda no desempenho, mas torna o código um pouco mais curto.

• Mais algumas mudanças relacionadas ao estilo do que ao desempenho.

Resultado:

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

Para mim, essa função é aproximadamente dez vezes mais rápida que a sua função R.

Eu posso pensar em quatro potenciais otimizações sobre as respostas de Ralf e Olivers.

(Você deve aceitar as respostas deles, mas eu só queria adicionar meus 2 centavos).

1) Use // [[Rcpp::export(rng = false)]] como cabeçalho de comentário para a função em um arquivo C ++ separado. Isso leva a uma velocidade de ~ 80% na minha máquina. (Esta é a sugestão mais importante das 4).

2) Preferir cmath quando possível. (Nesse caso, não parece fazer a diferença).

3) Evite alocação sempre que possível, por exemplo, não mude beta para um novo vetor.

4) meta de alongamento: use SEXP parâmetros em vez de vetores Rcpp. (Deixado como um exercício para o leitor). Os vetores Rcpp são invólucros muito finos, mas ainda são invólucros e há uma pequena sobrecarga.

Essas sugestões não seriam importantes, se não fosse pelo fato de você estar chamando a função em um loop restrito optim. Portanto, qualquer sobrecarga é muito importante.

Banco:

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 é a resposta de Oliver com rng=false. v4 está com as sugestões 2 e 3 incluídas.

A função:

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