从 Rcpp (Armadillo) 调用 glmnet答案

【问题标题】：Call glmnet from Rcpp (Armadillo)从 Rcpp (Armadillo) 调用 glmnet
【发布时间】：2019-10-25 02:05:37
【问题描述】：

我想在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(...)) ???
}

【问题讨论】：

从 C++ 调用 R 例程不会加速该过程。您是否考虑过只使用 R？如果不是，请考虑使用 coefs 结果调用 C++ 例程。
您对程序的速度是正确的。但我需要多次进行此计算。手动移动coefs 结果非常困难。
为什么手动移动coefs很难？你试过什么？
@Ralf Stubner 也许我让一个简单的问题变得困难。我不知道。但是我添加了我尝试过的代码。

标签： r rcpp glmnet method-call

【解决方案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 年 6 月 12 日创建}

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

【讨论】：