Converting between NumericVector/Matrix and VectorXd/MatrixXd in Rcpp(Eigen) to perform Cholesky solve

r, rcpp

Solution

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

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.

Problem

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. 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; ``` 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") ``` 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. 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! Thanks in advance, and also for all the great work on Rcpp(Eigen), David

Original source