在 Rcpp(Eigen) 中的 NumericVector/Matrix 和 VectorXd/MatrixXd 之间转换以执行 Cholesky 求解

Posted 2023-02-16

【中文标题】在 Rcpp(Eigen) 中的 NumericVector/Matrix 和 VectorXd/MatrixXd 之间转换以执行 Cholesky 求解

【英文标题】：Converting between NumericVector/Matrix and VectorXd/MatrixXd in Rcpp(Eigen) to perform Cholesky solve

【发布时间】：2013-12-06 23:42:01

【问题描述】：

编辑：从下面 Dirk 的回答中得到一些线索，我解决了这个问题，现在问题正文中的解决方案。

我确信这必须记录在某处，但我的谷歌技能让我失望。

我正在开发一个我认为不需要的 Rcpp 包 依赖于 Eigen，所以我非常广泛地使用了NumericVector/Matrix。 但是，我现在需要一个 Cholesky decomp/solve：我知道该怎么做 使用 RcppEigen 但我需要VectorXd/MatrixXd's。假设我有 P.S.D. NumericMatrix, A 和 NumericVector, b。我尝试了以下各种变化：

  Rcpp::NumericVector b(bb);
  Rcpp::NumericMatrix A(AA);
  Eigen::Map<Eigen::MatrixXd> A_eigen = as<Eigen::Map<Eigen::MatrixXd> >(A);
  Eigen::Map<Eigen::VectorXd> b_eigen = as<Eigen::Map<Eigen::VectorXd> >(b);
  Eigen::LLT<Eigen::MatrixXd> llt(A_eigen);
  Eigen::VectorXd x = llt.solve(b_eigen);
  Rcpp::NumericVector xx(wrap(x));
  return xx;

然后可以使用inline 包编译它（这里用codeString 表示），如下所示：

f = cxxfunction(signature(AA="matrix",bb="numeric"),codeString,plugin="RcppEigen")

我知道我可以直接从 SEXP 实例中获取 Eigen::MatrixXd/VectorXd（例如使用 Eigen::Map），但在我的使用中我需要从 NumericMatrix/Vector 实例中获取这些。

我的 A 会很小，所以我不会太小 担心是否会创建完整副本。如果有人能提供一个笨重的 解决方案，这将是一个很好的解决方案，一个优雅的解决方案将是 太好了！

在此先感谢您，也感谢您在 Rcpp(Eigen) 上所做的所有出色工作，

大卫

【问题讨论】：

您能否提供一个完整的可重现示例。并且可能显示你得到的编译器错误。

嗨罗曼 - 感谢您的回复。希望我添加的内容有所帮助。

【参考方案1】：

查看 RcppEigen 源和文件 src/fastLm.h 和 src/fastLm.cpp。

他们通过多种不同的分解方法建立了一个非常好的工厂。请注意，对于每个已实施的求解器，结果始终是受保护的 VectorXd m_coef（即 Eigen 类型），它仅在最后被转换：

        // Copy coefficients and install names, if any
NumericVector     coef(wrap(ans.coef()));
List          dimnames(NumericMatrix(Xs).attr("dimnames"));

完全可以想象，我们可以用不同的方式来做这件事——但是这个例子已经运行了很长一段时间了。我只会跟着它。

【讨论】：

谢谢！结合使用 `as' 去另一个方向做我需要的。

【参考方案2】：

Choleski 分解是在 base 中实现的。 Here 是文档的链接。如果我理解错了，那么也许你可以自己写一个？查看thisRosetta 代码链接，了解多种语言的 Choleski 分解，为您提供建议。