为什么我的 Rcpp 代码比 glmnet 慢得多？

Question

我从this 站点编辑了套索代码，以将其用于多个 lambda 值。 我将 lassoshooting 包用于一个 lambda 值（此包适用于一个 lambda 值），并使用 glmnet 用于多个 lambda 值进行比较。

系数估计值不同，这是预期的，因为标准化和缩减到原始规模。这超出了范围，在这里并不重要。

对于一种参数情况，lassoshooting 快 1.5 倍。

这两种方法都在我的代码中使用了所有 100 个 lambda 值来处理多个 lambda 情况。但是 glmnet 比我的 cpp 代码快 7.5 倍。当然，我预计 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

套索循环和 cpp 的结果：

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

Answer 1

包 {glmnet} 使用热启动和特殊规则来丢弃大量预测变量，这使得整个“正则化路径”的拟合速度非常快。

见their paper。