在 Rcpp(Eigen) 中的 NumericVector/Matrix 和 VectorXd/MatrixXd 之间转换以执行 Cholesky 求解 [英] Converting between NumericVector/Matrix and VectorXd/MatrixXd in Rcpp(Eigen) to perform Cholesky solve

问题描述

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

Edit: with some clues from Dirk's answer below I worked this out, solution in the body of the question now.

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

I'm sure this must be documented somewhere but my google skills are failing me.

我正在开发一个我认为不需要的 Rcpp 包对 Eigen 的依赖，所以我非常广泛地使用 NumericVector/Matrix.但是，我现在需要一个 Cholesky decomp/solve:我知道如何做到这一点使用 RcppEigen 但我需要 VectorXd/MatrixXd.说我有一个 P.S.D.NumericMatrix, A 和一个 NumericVector, b.我尝试了以下各种变体:

I'm developing an Rcpp package where I didn't think I would need dependency on Eigen, so I used NumericVector/Matrix quite extensively. However, I now need a Cholesky decomp/solve: I know how to do this using RcppEigen but I need VectorXd/MatrixXd's. Say I have a P.S.D. NumericMatrix, A, and a NumericVector, b. I have tried various variations on the following:

  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 表示)如下:

One can then compile this with (denoted here by codeString) using the inline package as follows:

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

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

I am aware that I could directly get the Eigen::MatrixXd/VectorXd from the SEXP instances (e.g. using Eigen::Map) but in my use I need to get these from NumericMatrix/Vector instances.

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

My A will be pretty small so I'm not too worried if a full copy gets created. If anyone can offer a clunky solution, that would be a good, and an elegant solution would be great!

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

Thanks in advance, and also for all the great work on Rcpp(Eigen),

大卫

推荐答案

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

Look at the RcppEigen sources and the files src/fastLm.h and src/fastLm.cpp.

他们通过多种不同的分解方法建立了一个非常好的工厂.请注意，对于每个实现的求解器，结果始终是受保护的 VectorXd m_coef(即特征类型)，它仅在最后被转换:

They set up a really nice factory over a number of different decomposition methods. Note how for each of the implemented solvers, the results is always the protected VectorXd m_coef (ie an Eigen type) which gets converted only at the very end:

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

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

It is quite conceivable that we could have done this differently -- but this example has been working fine for quite some time. I would just follow it.

这篇关于在 Rcpp(Eigen) 中的 NumericVector/Matrix 和 VectorXd/MatrixXd 之间转换以执行 Cholesky 求解的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！