Жадная оптимизация в R

Question

Я пытаюсь повторить метод Каруана и др. Для Ensemble selection from libraries of models (pdf). В основе метода лежит жадный алгоритм добавления моделей в ансамбль (модели могут быть добавлены более одного раза). Я написал реализацию этого жадного алгоритма оптимизации, но это очень медленно:

library(compiler) 
set.seed(42) 
X <- matrix(runif(100000*10), ncol=10) 
Y <- rnorm(100000) 

greedOpt <- cmpfun(function(X, Y, iter=100){ 
    weights <- rep(0, ncol(X)) 

    while(sum(weights) < iter) { 

    errors <- sapply(1:ncol(X), function(y){ 
     newweights <- weights 
     newweights[y] <- newweights[y] + 1 
     pred <- X %*% (newweights)/sum(newweights) 
     error <- Y - pred 
     sqrt(mean(error^2)) 
    }) 

    update <- which.min(errors) 
    weights[update] <- weights[update]+1 
    } 
    return(weights/sum(weights)) 
}) 

system.time(a <- greedOpt(X,Y)) 

Я знаю, что R не делает петлю хорошо, но я не могу думать о каком-либо способе сделать этот вид ступенчатого поиск без цикла.

Любые предложения по улучшению этой функции?

Answer 1

Вот реализация R, которая на 30% быстрее, чем ваша. Не так быстро, как ваша версия Rcpp, но, возможно, это даст вам идеи, которые в сочетании с Rcpp ускорят работу. Два основных усовершенствования:

1. петля sapply была заменена на матричной композиции
2. умножение матриц была заменена рекурсии

---

greedOpt <- cmpfun(function(X, Y, iter = 100L){ 

    N   <- ncol(X) 
    weights  <- rep(0L, N) 
    pred  <- 0 * X 
    sum.weights <- 0L 

    while(sum.weights < iter) { 

     sum.weights <- sum.weights + 1L 
     pred   <- (pred + X) * (1L/sum.weights) 
     errors  <- sqrt(colSums((pred - Y)^2L)) 
     best   <- which.min(errors) 
     weights[best] <- weights[best] + 1L 
     pred   <- pred[, best] * sum.weights 
    } 
    return(weights/sum.weights) 
}) 

Кроме того, я поддерживать вы должны попробовать перейти на библиотеку атласа. Вы можете увидеть значительные улучшения.

Answer 2

Я выстрелил в написании версии Rcpp этой функции:

library(Rcpp) 
cppFunction(' 
    NumericVector greedOptC(NumericMatrix X, NumericVector Y, int iter) { 
    int nrow = X.nrow(), ncol = X.ncol(); 
    NumericVector weights(ncol); 
    NumericVector newweights(ncol); 
    NumericVector errors(nrow); 
    double RMSE; 
    double bestRMSE; 
    int bestCol; 

    for (int i = 0; i < iter; i++) { 
     bestRMSE = -1; 
     bestCol = 1; 
     for (int j = 0; j < ncol; j++) { 
     newweights = weights + 0; 
     newweights[j] = newweights[j] + 1; 
     newweights = newweights/sum(newweights); 

     NumericVector pred(nrow); 
     for (int k = 0; k < ncol; k++){ 
      pred = pred + newweights[k] * X(_, k); 
     } 

     errors = Y - pred; 
     RMSE = sqrt(mean(errors*errors)); 

     if (RMSE < bestRMSE || bestRMSE==-1){ 
      bestRMSE = RMSE; 
      bestCol = j; 
     } 
     } 

     weights[bestCol] = weights[bestCol] + 1; 
    } 

    weights = weights/sum(weights); 
    return weights; 
    } 
') 

Это более чем в два раза быстрее, чем в версии R:

set.seed(42) 
X <- matrix(runif(100000*10), ncol=10) 
Y <- rnorm(100000) 
> system.time(a <- greedOpt(X, Y, 1000)) 
    user system elapsed 
    36.19 6.10 42.40 
> system.time(b <- greedOptC(X, Y, 1000)) 
    user system elapsed 
    16.50 1.44 18.04 
> all.equal(a,b) 
[1] TRUE 

Не плохо, но я надеялся на большее ускорение при совершении прыжка с R на Rcpp. Это одна из первых функций Rcpp, которые я когда-либо писал, поэтому возможно дальнейшая оптимизация.