Search code examples
reigenrcpp

How do you convert object of class Eigen::MatrixXd to class Rcpp::NumericMatrix


I'm working on a package that requires some very fast matrix multiplication so looking to use RcppEigen. For a variety of reasons though having to do with the need for multidimensional arrays, I need to convert a created object of class Eigen::MatrixXd to class Rcpp::NumericMatrix.

I tried reversing the steps listed in RcppEigen::FastLm.cpp, but that doesn't seem to work

e.g. instead of using

const Map<MatrixXd> X(as<Map<MatrixXd> >(Xs));

I tried

Rcpp:NumericMatrix X(as<Rcpp::NumericMatrix>(Xs));

where Xs is a matrix of class Eigen::MatrixXd but that didn't seem to work:" error: no matching function for call to 'as' return Rcpp::asRcpp::NumericMatrix(z);"

If this isn't at all possible I can try another direction.

Basically what I need to do in R speak is

a = matrix(1, nrow = 10, ncol = 10)

b = array(0, c(10,10,10))

b[,,1] = a

To give a clearer starting example

How would I go about storing an object of class MatrixXd in an object of class NumericMatrix?

#include <Rcpp.h>
#include <RcppEigen.h>
using namespace Rcpp;
using namespace Eigen;

// [[Rcpp::export]]
NumericMatrix sample_problem() {

  Eigen::MatrixXd x(2, 2); x << 1,1,2,2;

  Eigen::MatrixXd z(2, 2);

  Eigen::MatrixXd y(2,2); y << 3,3,4,4;

  z =  x * y; // do some eigen matrix multiplication

  Rcpp::NumericMatrix w(2,2);

  // what I'd like to be able to do somehow: 
  // store the results of the eigen object z in
  // a NumericMatrix w
  // w = z;

  return w;
} 


Solution

  • Thanks for posting code! It makes everything easier. I just rearranged you code the tiniest bit.

    The key changes is to "explicitly" go back from the Eigen representation via an RcppEigen helper to a SEXP, and to then instantiate the matrix. Sometimes ... the compiler needs a little nudge.

    Code

    #include <Rcpp.h>
    #include <RcppEigen.h>
    
    // [[Rcpp::depends(RcppEigen)]]
    
    // [[Rcpp::export]]
    Rcpp::NumericMatrix sample_problem() {
    
      Eigen::MatrixXd x(2, 2), y(2, 2);
      x << 1,1,2,2;
      y << 3,3,4,4;
    
      // do some eigen matrix multiplication
      Eigen::MatrixXd z =  x * y;
    
      // what I'd like to be able to do somehow:
      // store the results of the eigen object z in
      // a NumericMatrix w
      // w = z;
      SEXP s = Rcpp::wrap(z);
      Rcpp::NumericMatrix w(s);
    
      return w;
    }
    
    /*** R
    sample_problem()
    */
    

    Demo

    R> sourceCpp("demo.cpp)
    
    R> sample_problem()
         [,1] [,2]
    [1,]    7    7
    [2,]   14   14
    R> 
    

    With g++-9 I need -Wno-ignored-attributes or I get screens and screens of warnings...