从 Rcpp (Armadillo) 调用 glmnet

Question

我想在 Rcpp Armadillo 中提取 glmnet 的系数估计值（交叉验证后），以便在 Armadillo 的另一个函数中使用它们。我搜索了类似的问题，但找不到解决方案。

附上我的尝试。（不工作）即使我能得到 cv.glmnet 的列表结果，我也无法使用 coef 函数来获取系数。

R码

library(glmnet)

set.seed(1)
X = matrix(rnorm(1e3 * 201), 1e3, 201)
beta = -100:100
y = X%*%beta + rnorm(1e3)
cvfit = cv.glmnet(X, y, alpha = 1)
coefs = coef(cvfit, s = "lambda.min")
coefs                                   # get these coefficients from Rcpp

cv.glmnet

的参数

args(cv.glmnet)
> function (x, y, weights, offset = NULL, lambda = NULL, type.measure = c("mse", "deviance", "class", "auc", "mae"), nfolds = 10, foldid, 
    alignment = c("lambda", "fraction"), grouped = TRUE, keep = FALSE, 
    parallel = FALSE, ...) 
NULL

C++代码

#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]

// [[Rcpp::export]]
Rcpp::List f_cpp(const arma::mat &x, const arma::vec &y, 
                 const arma::vec &weights, 
                 const arma::vec &lambda, double alpha, 
                 int nfolds = 10){

  Rcpp::Environment pkg = Rcpp::Environment::namespace_env("glmnet");

  Rcpp::Function f_R = pkg["cv.glmnet"];

  Rcpp::Nullable<arma::vec> offset = pkg["offset"];
  Rcpp::CharacterVector type_measure = pkg["type.measure"];
  arma::vec foldid = pkg["foldid"];
  Rcpp::CharacterVector alignment = pkg["alignment"];
  bool grouped = pkg["grouped"];
  bool keep = pkg["keep"];
  bool parallel = pkg["parallel"];

  return f_R(x, y, weights, offset, lambda, 
             type_measure, nfolds, foldid, 
             alignment, grouped, keep, parallel, alpha = alpha);
  // coef(f_R(...)) ???
}

Answer 1

从 C++ 调用像 cv.glmnet 这样的函数很复杂（甚至可能是不可能的），因为它使用了 R 提供的许多可能性，使函数签名非常灵活。但是，可以在 R 中定义一个使用实际使用的签名的包装函数。与其从（全局）环境中获取此函数，我更喜欢将其作为函数参数传递：

library(glmnet)
#> Loading required package: Matrix
#> Loading required package: foreach
#> Loaded glmnet 2.0-16

set.seed(1)
X = matrix(rnorm(1e3 * 201), 1e3, 201)
beta = -100:100
y = X%*%beta + rnorm(1e3)


# set seed since cv.glmnet uses random numbers
set.seed(1)
cvfit = cv.glmnet(X, y, alpha = 1)
coefs = coef(cvfit, s = "lambda.min")

# set seed since cv.glmnet uses random numbers
set.seed(1)
my.glmnet <- function(x, y, alpha) {
    cvfit <- cv.glmnet(x, y, alpha = alpha)
    coef(cvfit, s = "lambda.min")
}
Rcpp::cppFunction(depends = "RcppArmadillo", "
arma::sp_mat f_cpp(const arma::mat &x, const arma::vec &y, double alpha, Rcpp::Function f_R) {
    arma::sp_mat coef = Rcpp::as<arma::sp_mat>(f_R(x, y, alpha));
    return coef;
}")
coefs2 <- f_cpp(X, y, alpha = 1, my.glmnet)

all(coefs - coefs2 == 0)
#> [1] TRUE

^{由 reprex package (v0.3.0)}

创建于 2019-06-12

当然，您可以对计算出的系数做比将它们返回到 R 更有趣的事情。显式 Rcpp::as 是必需的，因为 C++ 无法知道 R 函数的参数类型 returns.在这种情况下它是一个稀疏矩阵，可以转换为 arma::sp_mat。顺便说一句，这会丢失矩阵的 Dimnames，这就是为什么不能使用 all.equal 进行比较的原因。

从 Rcpp (Armadillo) 调用 glmnet

Call glmnet from Rcpp (Armadillo)

r

rcpp

glmnet

method-call