RcppParallel calculation of lm coef

library(Rcpp)
library(RcppArmadillo)
library(RcppParallel)

sourceCpp(code = "
// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>
using namespace Rcpp;
using namespace arma;
// [[Rcpp::depends(RcppParallel)]]
#include <RcppParallel.h>
using namespace RcppParallel;

struct coef_compute: public Worker {
  const mat X;
  const vec& y;
  const umat& combnMat;
  mat& coefMat;
  coef_compute(const mat& X, const vec& y, const umat& combnMat,mat& coefMat):
    X(X), y(y), combnMat(combnMat), coefMat(coefMat) {}
  void operator()(std::size_t begin, std::size_t end)
  {
    for (uword i = begin; i < end; i++)
    coefMat.row(i) = solve(join_rows(ones<vec>(y.n_elem),
                           X.cols(combnMat.col(i))), y).t();
  }
};


// [[Rcpp::export]]
List fastCoef_allCombns(arma::vec y, arma::mat x, int k)
{
  Rcpp::Function combn_rcpp(\"combn\");
  umat combnMat = conv_to<umat>::from(as<mat>(combn_rcpp(x.n_cols, k))) - 1;
  mat coefMat = zeros<mat>(combnMat.n_cols, k + 1);
  coef_compute coefResults(x, y, combnMat, coefMat);
  parallelFor(0, combnMat.n_cols, coefResults);
  return List::create(combnMat.t() + 1, coefMat);
}")

n <- 30
X <- matrix(rnorm(120), n)
coefs <- c(2, 5, -7, 8, 0.5)
y <- as.vector(cbind(1, X) %*% coefs) + rnorm(n, 0, 5)

outRes <- lapply(1:ncol(X), function(i){
  fastCoef_allCombns(y, X, i)
})

outRes[[2]]
# # 哪些X的組合
# [[1]]
#      [,1] [,2]
# [1,]    1    2   # X1 跟 X2的組合
# [2,]    1    3
# [3,]    1    4
# [4,]    2    3
# [5,]    2    4
# [6,]    3    4
#
# # 每一個組合的X，其係數，BY列去看
# [[2]]
#           [,1]      [,2]       [,3]
# [1,] 0.4210033  6.445250 -7.5248928
# [2,] 0.7895878  7.378273  8.7397441
# [3,] 0.6180885  6.864235  0.9617974
# [4,] 2.6213473 -6.677914  6.3579295
# [5,] 2.3471788 -7.834296  1.4096750
# [6,] 2.9802756  8.212602  0.9764216

# benchmark
n <- 3000; p <- 10
X <- matrix(rnorm(n*p), n)
coefs <- sample(-10:10, p+1, TRUE)
y <- as.vector(cbind(1, X) %*% coefs) + rnorm(n, 0, 5)

st <- proc.time()
outRes <- lapply(1:ncol(X), function(i){
  fastCoef_allCombns(y, X, i)
})
proc.time() - st
# user  system elapsed
# 0.14    0.08    0.07

st <- proc.time()
outRes2 <- lapply(1:ncol(X), function(i){
  combnMat <- combn(1:ncol(X), i)
  return(list(t(combnMat), t(apply(combnMat, 2, function(v){
    coef(lm(y ~ X[,v]))
  }))))
})
proc.time() - st
# user  system elapsed
# 3.57    0.19    3.56