在Rcpp（Eigen）中，将NumericVector / Matrix和VectorXd / MatrixXd相互转换以执行Cholesky求解。

Question

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

我相信这一定有记录在案的，但是我的谷歌技能正在失败。

我正在开发一个 Rcpp 包，在那里我认为我不需要依赖 Eigen，因此我相当广泛地使用了 NumericVector/Matrix。然而，现在我需要进行 Cholesky 分解/求解：我知道如何使用 RcppEigen 进行操作，但我需要 VectorXd/MatrixXd。假设我有一个 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 :: Map），但在我的应用中，我需要从NumericMatrix / Vector实例中获取这些内容。

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

提前感谢您的回复，并且感谢您在Rcpp（Eigen）上做出的所有伟大工作，

戴维

- daknowles

---

你能提供一个完整的可重现的例子吗？也许可以展示一下你得到的编译器错误。 - Romain Francois

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

Answer 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"));

可以想象我们本可以采用不同的方法——但是这个例子已经运作良好很长一段时间了。我建议您跟随这个例子。

- Dirk Eddelbuettel

---

谢谢！这个结合使用 as 来实现双向转换正是我所需要的。 - daknowles

Answer 2

Cholesky分解已经在base中实现。这里是文档链接。如果我有误解，那么也许你可以自己写一个？请查看这个Rosetta Code链接以获取关于Cholesky分解的不同语言的示例。

- pstew