rcpp - Converting an Armadillo Matrix to an Eigen MatriXd and vice versa

Question

How can I convert from an Armadillo Matrix to an Eigen MatrixXd and vice versa?

I have nu as an arma::vec of size N, z as arma::mat of dimension N x 3. I want to compute a matrix P such as the entry P_ij is

Pij=exp(nu(i) + nu(j) + z.row(j)*z.row(j)))

Thus I used this code

int N=z.n_rows;
mat P= exp(nu*ones(1,N) + one(N,1)*(nu.t()) + z*(z.t()));

But the computation takes too long. In particular, for N = 50,000 the run time is far to high.

It seems that using Eigen can be faster. But my matrix are Armadillo. How can I use Eigen operations ? Or how can I do this operation faster.

Answer 1

Using armadillo's .memptr() class member function, we are able to extract the memory pointer. From here, we can use Eigen's Map<T>() constructor to create an Eigen matrix.

Now, we can go from the Eigen matrix using the .data() member function to extract a point to Eigen's memory structure. Then, using the advanced constructor options of arma::mat we can create an armadillo matrix.

For example:

#include <RcppArmadillo.h>
#include <RcppEigen.h>

// [[Rcpp::depends(RcppEigen)]]

// [[Rcpp::depends(RcppArmadillo)]]

// [[Rcpp::export]]
Eigen::MatrixXd example_cast_eigen(arma::mat arma_A) {

  Eigen::MatrixXd eigen_B = Eigen::Map<Eigen::MatrixXd>(arma_A.memptr(),
                                                        arma_A.n_rows,
                                                        arma_A.n_cols);

  return eigen_B;
}

// [[Rcpp::export]]
arma::mat example_cast_arma(Eigen::MatrixXd eigen_A) {

  arma::mat arma_B = arma::mat(eigen_A.data(), eigen_A.rows(), eigen_A.cols(),
                               false, false);

  return arma_B;
}

/***R
(x = matrix(1:4, ncol = 2))
example_cast_eigen(x)
example_cast_arma(x)
*/

Results:

(x = matrix(1:4, ncol = 2))
#      [,1] [,2]
# [1,]    1    3
# [2,]    2    4

example_cast_eigen(x)
#      [,1] [,2]
# [1,]    1    3
# [2,]    2    4

example_cast_arma(x)
#      [,1] [,2]
# [1,]    1    3
# [2,]    2    4

One quick remark: If you are using Eigen's Mapping function, then you should automatically have the change in the Armadillo matrix (and vice versa), e.g.

#include <RcppArmadillo.h>
#include <RcppEigen.h>

// [[Rcpp::depends(RcppEigen)]]

// [[Rcpp::depends(RcppArmadillo)]]

// [[Rcpp::export]]
void map_update(Eigen::MatrixXd eigen_A) {

  Rcpp::Rcout << "Eigen Matrix on Entry: " << std::endl << eigen_A << std::endl;

  arma::mat arma_B = arma::mat(eigen_A.data(), eigen_A.rows(), eigen_A.cols(),
                               false, false);

  arma_B(0, 0) = 10;
  arma_B(1, 1) = 20;

  Rcpp::Rcout << "Armadill Matrix after modification: " << std::endl << arma_B << std::endl;

  Rcpp::Rcout << "Eigen Matrix after modification: " << std::endl << eigen_A << std::endl;
}

Run:

map_update(x)

Output:

Eigen Matrix on Entry: 
1 3
2 4

Armadill Matrix after modification: 
   10.0000    3.0000
    2.0000   20.0000

Eigen Matrix after modification: 
10  3
 2 20