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

Question

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;
} 

Manos

Can you focus and sharpen the question and actually show code? Ex ante a multidimensional array simply won't fit into a 2d NumericMatrix so maybe you can make it more concrete about 2-d slices of a multi-d array? And maybe once you have 2-d matrix slices or view in Eigen it will become more obvious how those become a SEXP which can of course become a NumericMatrix (and all can happen zero-copy). – Volant

Thanks, see hopefully sharpened question with code added to original post – Manos

Perfect, see below. – Volant

Answer 1

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...

Volant

That's exactly what I was looking for, thanks! Now to tackle assignment to a multidimensional array, but this should get me going in the right direction – Manos

Slight issue with this I think: running into an error "Not compatible with requested type: [type=any; target=double]." seemingly at random; running the exact same function again after encountering that error works just fine. There are no stochastic processes in the model. Seems potentially related to this issue? #32376037. Apologies that this is a vague question but I am stumped; any general thoughts/tips are appreciated! – Manos

No idea, doing Rcpp::sourceCpp("~/git/stackoverflow/62586950/answer.cpp") still works fine here (modulo the lots of Eigen compilation noise). – Volant

Sorry, should have clarified, the error is in the code I built based on this suggestion (github.com/DanOvando/marlin/blob/…), but guessing if it's not an obvious and known problem to you it's something very weird I need to do more digging on. thanks – Manos