Почему мой код Rcpp намного медленнее, чем у glmnet?

Я редактировал код Лассо из этот сайт, чтобы использовать его для нескольких значений лямбда.
Я использовал пакет lassoshooting для одного лямбда-значения (этот пакет работает для одного лямбда-значения) и glmnet для нескольких лямбда-значений для сравнения.

Оценки коэффициентов отличаются, и это ожидается из-за стандартизации и масштабирования до первоначального масштаба. Это выходит за рамки и не важно здесь.

Для случая с одним параметром lassoshooting в 1,5 раза быстрее.

Оба метода использовали все 100 лямбда-значений в моем коде для множественного лямбда-случая. Но glmnet в 7,5 раз быстрее, чем мой код cpp. Конечно, я ожидал, что glmnet будет быстрее, но эта сумма кажется слишком большой. Это нормально или мой код неверен?

РЕДАКТИРОВАТЬ

Я также приложил lshoot функция, которая вычисляет коэффициент пути в петле R. Это превосходит мой код cpp тоже.

Могу ли я улучшить свой код cpp?

Код C ++:

// [[Rcpp::depends(RcppArmadillo)]]

#include <RcppArmadillo.h>
using namespace Rcpp;
using namespace arma;// [[Rcpp::export]]
vec softmax_cpp(const vec & x, const vec & y) {
return sign(x) % max(abs(x) - y, zeros(x.n_elem));
}

// [[Rcpp::export]]
mat lasso(const mat & X, const vec & y, const vec & lambda,
const double tol = 1e-7, const int max_iter = 10000){
int p = X.n_cols; int lam = lambda.n_elem;
mat XX = X.t() * X;
vec Xy = X.t() * y;
vec Xy2 = 2 * Xy;
mat XX2 = 2 * XX;
mat betas = zeros(p, lam); // to store the betas

vec beta = zeros(p); // initial beta for each lambda

bool converged = false;
int iteration = 0;
vec beta_prev, aj, cj;

for(int l = 0; l < lam; l++){
while (!converged && (iteration < max_iter)){

beta_prev = beta;

for (int j = 0; j < p; j++){
aj = XX2(j,j);
cj = Xy2(j) - dot(XX2.row(j), beta) + beta(j) * XX2(j,j);
beta(j) = as_scalar(softmax_cpp(cj / aj, as_scalar(lambda(l)) / aj));
}
iteration = iteration + 1;
converged =  norm(beta_prev - beta, 1) < tol;
}
betas.col(l) = beta;
iteration = 0;
converged = false;
}
return betas;
}

Код R:

library(Rcpp)
library(rbenchmark)
library(glmnet)
library(lassoshooting)

sourceCpp("LASSO.cpp")

library(ElemStatLearn)
X <- as.matrix(prostate[,-c(9,10)])
y <- as.matrix(prostate[,9])
lambda_one <- 0.1
benchmark(cpp=lasso(X,y,lambda_one),
lassoshooting=lassoshooting(X,y,lambda_one)$coefficients,
order="relative", replications=100)[,1:4]

################################################
lambda <- seq(0,10,len=100)

benchmark(cpp=lasso(X,y,lambda),
glmn=coef(glmnet(X,y,lambda=lambda)),
order="relative", replications=100)[,1:4]

####################################################

РЕДАКТИРОВАТЬ

lambda <- seq(0,10,len=100)

lshoot <- function(lambda){
betas <- matrix(NA,8,100)
for(l in 1:100){
betas[, l] <- lassoshooting(X,y,lambda[l])$coefficients
}
return(betas)
}

benchmark(cpp=lasso(X,y,lambda),
lassoshooting_loop=lshoot(lambda),
order="relative", replications=300)[,1:4]

Результаты для одного параметра:

           test replications elapsed relative
2 lassoshooting          300    0.06      1.0
1           cpp          300    0.09      1.5

Результаты для случая с несколькими параметрами:

  test replications elapsed relative
2 glmn          300    0.70    1.000
1  cpp          300    5.24    7.486

Результаты для цикла lassoshooting и cpp:

                test replications elapsed relative
2 lassoshooting_loop          300    4.06    1.000
1                cpp          300    6.38    1.571