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!
fuente
ObscomoIntegerVectoreliminar algunos moldes.beta) permanecen constantes durante las observacionesObs. Si tuviéramos predictores variables en el tiempo, sería necesario un cálculo implícito dedenomcada unoObs, basado en el valor de una matriz de diseñoX. 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!for-loopque le dará la mayor ganancia. También en el caso más general, sugeriría usarmodel.matrix(...)para crear su matriz para la entrada en sus funciones.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 todosi. Por lo tanto, tiene sentido usar aydouble 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 un
IntegerVectorpara empezar hace que muchos lanzamientos sean innecesarios.También puede indexar
NumericVectorusando unIntegerVectoren 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:
Para mí, esta función es aproximadamente diez veces más rápida que su función R.
fuente
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
cmathcuando sea posible. (En este caso, no parece hacer la diferencia).3) Evite la asignación siempre que sea posible, por ejemplo, no cambie
betaa 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:
v3 es la respuesta de Oliver con
rng=false. v4 incluye las sugerencias n. ° 2 y n. ° 3.La función:
fuente