Greedy optimization in R

Question

I am trying to replicate Caruana et al.'s method for Ensemble selection from libraries of models (pdf). At the core of the method is a greedy algorithm for adding models to the ensemble (models can be added more than once). I've written an implementation for this greedy optimization algorithm, but it is very slow:

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))

I know R doesn't do loops well, but I can't think of any way to do this type of stepwise search without a loop.

Any suggestions for improving this function?

Answer 1

Here is an R implementation that is 30% faster than yours. Not as fast as your Rcpp version but maybe it will give you ideas that combined with Rcpp will speed things further. The two main improvements are:

1. the sapply loop has been replaced by a matrix formulation
2. the matrix multiplication has been replaced by a recursion

---

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

Also, I maintain you should try upgrading to the atlas library. You might see significant improvements.

Answer 2

I took a shot at writing an Rcpp version of this function:

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

It's more than twice as fast as the R version:

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

Not bad, but I was hoping for a bigger speedup when making the leap from R to Rcpp. This is one of the first Rcpp functions I've ever written, so perhaps further optimization is possible.