Greedy optimization in R
I am trying to replicate Caruana et al's method of Selecting an Ensemble from Model Libraries (pdf) . The method is based on a greedy algorithm for adding models to an ensemble (models can be added more than once). I wrote 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 iterative search without a loop.
Any suggestions for improving this feature?
source to share
Here's an R implementation that's 30% faster than yours. Not as fast as your version of Rcpp, but maybe it will give you ideas that when combined with Rcpp will make things faster. The two main improvements are:
- the cycle has
sapply
been replaced by matrix formulation - matrix multiplication has been replaced by 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 submit that you should try switching to the atlas library. You can see significant improvements.
source to share
I made an attempt to write 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;
}
')
This is 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 more acceleration when jumping from R to Rcpp. This is one of the first Rcpp functions I have ever written, so further optimization is possible.
source to share