Optimizando la función objetivo de R con Rcpp más lento, ¿por qué?

16

Actualmente estoy trabajando en un método bayesiano que requiere múltiples pasos de optimización de un modelo logit multinomial por iteración. Estoy usando optim () para realizar esas optimizaciones, y una función objetivo escrita en R. Un perfil reveló que optim () es el principal cuello de botella.

Después de investigar, encontré esta pregunta en la que sugieren que recodificar la función objetivo Rcpppodría acelerar el proceso. Seguí la sugerencia y recodifiqué mi función objetivo Rcpp, pero terminó siendo más lenta (¡aproximadamente dos veces más lenta!).

Esta fue mi primera vez con Rcpp(o cualquier cosa relacionada con C ++) y no pude encontrar una manera de vectorizar el código. ¿Alguna idea de cómo hacerlo más rápido?

Tl; dr: La implementación actual de la función en Rcpp no ​​es tan rápida como la R vectorizada; ¿Cómo hacerlo más rápido?

Un ejemplo reproducible. :

1) Definir funciones objetivas en Ry Rcpp: log-verosimilitud de una intercepción solo modelo multinomial

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 su eficiencia:

## 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) Ahora llamándolos 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

Me sorprendió un poco que la implementación vectorizada en R fuera más rápida. ¿Implementar una versión más eficiente en Rcpp (por ejemplo, con RcppArmadillo?) Puede producir ganancias? ¿Es una mejor idea recodificar todo en Rcpp usando un optimizador de C ++?

PD: publicación por primera vez en Stackoverflow!

smildiner
fuente

Respuestas:

9

En general, si puede usar funciones vectorizadas, encontrará que es (casi) tan rápido como ejecutar su código directamente en Rcpp. Esto se debe a que muchas funciones vectorizadas en R (casi todas las funciones vectorizadas en Base R) están escritas en C, Cpp o Fortran y, como tal, a menudo hay poco que ganar.

Dicho esto, hay mejoras para ganar tanto en su Ry Rcppcódigo. La optimización proviene de estudiar cuidadosamente el código y eliminar pasos innecesarios (asignación de memoria, sumas, etc.).

Comencemos con la Rcppoptimización del código.

En su caso, la optimización principal es eliminar los cálculos innecesarios de matrices y vectores. El código es en esencia.

  1. Shift beta
  2. calcular el registro de la suma de exp (shift beta) [log-sum-exp]
  3. usa Obs como índice para la beta desplazada y suma todas las probabilidades
  4. reste el log-sum-exp

Usando esta observación podemos reducir su código a 2 for-loops. Tenga en cuenta que sumes simplemente otro bucle for (más o menos for(i = 0; i < max; i++){ sum += x }:), por lo que evitar las sumas puede acelerar los códigos (¡en la mayoría de las situaciones, esto es una optimización innecesaria!). Además, su entrada Obses un vector entero, y podemos optimizar aún más el código utilizando el IntegerVectortipo para evitar convertir los doubleelementos en integervalores (crédito a la respuesta 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;
}')

Tenga en cuenta que he eliminado bastantes asignaciones de memoria y eliminado cálculos innecesarios en el ciclo for. También he usado esodenom es lo mismo para todas las iteraciones y simplemente multiplicado por el resultado final.

Podemos realizar optimizaciones similares en su código R, lo que da como resultado la siguiente función:

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
}

Tenga en cuenta que la complejidad de la función se ha reducido drásticamente, lo que facilita la lectura de otros. Solo para asegurarme de que no he estropeado el código en alguna parte, verifiquemos que devuelvan los mismos 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

bueno eso es un alivio.

Actuación:

Usaré microbenchmark para ilustrar el rendimiento. Las funciones optimizadas son rápidas, por lo que ejecutaré los 1e5tiempos de las funciones para reducir el efecto del recolector de basura

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

Aquí vemos el mismo resultado que antes. Ahora las nuevas funciones son aproximadamente 35 veces más rápidas (R) y 40 veces más rápidas (Cpp) en comparación con sus primeras contrapartes. Curiosamente, la Rfunción optimizada sigue siendo muy ligeramente (0.3 ms o 4%) más rápida que mi Cppfunción optimizada . Mi mejor apuesta aquí es que hay algunos gastos generales delRcpp paquete, y si se eliminara, los dos serían idénticos o la R.

Del mismo modo, podemos verificar el rendimiento con 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

Una vez más, el resultado es el mismo.

Conclusión:

Como conclusión breve, vale la pena señalar que este es un ejemplo, donde la conversión de su código a Rcpp realmente no vale la pena. Este no es siempre el caso, pero a menudo vale la pena echar un segundo vistazo a su función, para ver si hay áreas de su código, donde se realizan cálculos innecesarios. Especialmente en situaciones en las que uno usa funciones vectorizadas buildin, a menudo no vale la pena convertir el código a Rcpp. Con mayor frecuencia, se pueden ver grandes mejoras si se usa un for-loopscódigo que no se puede vectorizar fácilmente para eliminar el bucle for.

Oliver
fuente
1
Puede tratar Obscomo IntegerVectoreliminar algunos moldes.
Ralf Stubner
Solo lo estaba incorporando antes de agradecerle por notar esto en su respuesta. Simplemente pasó por mí. Te he dado crédito por esto en mi respuesta @RalfStubner. :-)
Oliver
2
Como notó en este ejemplo de juguete (modelo mnl de solo intercepción), los predictores lineales ( beta) permanecen constantes durante las observaciones Obs. Si tuviéramos predictores variables en el tiempo, sería necesario un cálculo implícito de denomcada uno Obs, basado en el valor de una matriz de diseño X. Dicho esto, ya estoy implementando sus sugerencias en el resto de mi código con algunas ganancias realmente agradables :). ¡Gracias @RalfStubner, @Oliver y @thc por sus muy perspicaces respuestas! Ahora pasando a mi próximo cuello de botella!
smildiner
1
Me alegra que pudiéramos ayudar. En el caso más general, calcular el sustrato de resta en cada paso del segundo for-loopque le dará la mayor ganancia. También en el caso más general, sugeriría usar model.matrix(...)para crear su matriz para la entrada en sus funciones.
Oliver
9

Su función C ++ se puede hacer más rápido usando las siguientes observaciones. Al menos el primero también podría usarse con su función R:

  • La forma de calcular denom[i]es la misma para todos i. Por lo tanto, tiene sentido usar ay double denomhacer este cálculo solo una vez. También factorizo ​​restando este término común al final.

  • Sus observaciones son en realidad un vector entero en el lado R, y también las está usando como números enteros en C ++. Usando unIntegerVector para empezar hace que muchos lanzamientos sean innecesarios.

  • También puede indexar NumericVectorusando un IntegerVectoren C ++. No estoy seguro de si esto ayuda al rendimiento, pero hace que el código sea un poco más corto.

  • Algunos cambios más que están más relacionados con el estilo que el rendimiento.

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 mí, esta función es aproximadamente diez veces más rápida que su función R.

Ralf Stubner
fuente
Gracias por su respuesta Ralph, no detectó el tipo de entrada. También he incorporado esto a mi respuesta para darle el crédito. :-)
Oliver
7

Puedo pensar en cuatro optimizaciones potenciales sobre las respuestas de Ralf y Olivers.

(Debería aceptar sus respuestas, pero solo quería agregar mis 2 centavos).

1) uso // [[Rcpp::export(rng = false)]] como encabezado de comentario para la función en un archivo C ++ separado. Esto conduce a una velocidad de ~ 80% en mi máquina. (Esta es la sugerencia más importante de las 4).

2) Prefiero cmath cuando sea posible. (En este caso, no parece hacer la diferencia).

3) Evite la asignación siempre que sea posible, por ejemplo, no cambie beta a un nuevo vector.

4) Objetivo de estiramiento: use SEXPparámetros en lugar de vectores Rcpp. (Izquierda como ejercicio para el lector). Los vectores Rcpp son envoltorios muy delgados, pero siguen siendo envoltorios y hay una pequeña sobrecarga.

Estas sugerencias no serían importantes, si no fuera por el hecho de que está llamando a la función en un ciclo cerrado optim. Por lo tanto, cualquier sobrecarga es muy 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 es la respuesta de Oliver con rng=false. v4 incluye las sugerencias n. ° 2 y n. ° 3.

La función:

#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
fuente