LOOCV with Rcpp (R file)

Rcpp::sourceCpp("RcppLOOCV.cpp")

N <- 600L
p <- 150L
X <- matrix(rnorm(N*p), N)
b <- rnorm(p, rnorm(2L), rgamma(2L,3,7))
y <- 1.5 + as.vector(X %*% b) + rnorm(N)

system.time({
  out <- sapply(1L:N, function(i){
    idx <- setdiff(1L:N, i)
    b <- coef(lm(y[idx] ~ X[idx , ]))
    return((y[i] - crossprod(c(1, X[i, ]), b))**2)
  })
})
#  user  system elapsed
#  8.83    0.28    9.20

# check results
for (i in 1L:5L) stopifnot(abs(mean(out) - arma_fastLOOCV1(y, X)) < 1e-8)
for (i in 1L:5L) stopifnot(abs(mean(out) - arma_fastLOOCV2(y, X)) < 1e-8)
for (i in 1L:5L) stopifnot(abs(mean(out) - eigen_fastLOOCV1(y, X)) < 1e-8)
for (i in 1L:5L) stopifnot(abs(mean(out) - eigen_fastLOOCV2(y, X)) < 1e-8)
for (i in 1L:5L) stopifnot(abs(mean(out) - eigen_fastLOOCV3(y, X)) < 1e-8)
for (i in 1L:5L) stopifnot(abs(mean(out) - eigen_fastLOOCV4(y, X)) < 1e-8)

library(microbenchmark)
microbenchmark(
  arma_fastLOOCV1(y, X),  # RcppArmadillo with RcppPrallel
  arma_fastLOOCV2(y, X),  # RcppArmadillo with openmp
  eigen_fastLOOCV1(y, X), # RcppEigen with RcppPrallel (Reduce)
  eigen_fastLOOCV2(y, X), # RcppEigen with RcppPrallel (For)
  eigen_fastLOOCV3(y, X), # RcppEigen with openmp
  eigen_fastLOOCV4(y, X), # RcppEigen without parallism
  times = 30L
)
# Unit: milliseconds
#                    expr       min        lq      mean    median        uq       max neval
#   arma_fastLOOCV1(y, X) 1009.9239 1037.4614 1052.9272 1045.9733 1068.5668 1101.5399    30
#   arma_fastLOOCV2(y, X) 1020.5068 1033.6180 1051.7171 1046.8751 1064.4512 1139.6371    30
#  eigen_fastLOOCV1(y, X)  623.7900  628.0837  639.1062  635.0657  644.8178  716.8571    30
#  eigen_fastLOOCV2(y, X)  620.1337  624.2227  634.4193  630.3217  641.1066  682.0993    30
#  eigen_fastLOOCV3(y, X)  620.0021  630.6396  643.2877  638.7085  655.9284  706.5759    30
#  eigen_fastLOOCV4(y, X) 2627.1159 2648.5870 2657.8932 2654.8032 2663.2681 2707.9168    30