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.

Squalene

To speed up matrix multiplication within Armadillo, link with OpenBLAS instead of standard BLAS. This makes a massive difference. Also, your code is using the exp() function, which can be time consuming. You can speed it up by enabling OpenMP in Armadillo. – Intensive

Tanks I am not very beginner here. Please can you give a simple example ? – Squalene

Read through the README.txt file that comes with the Armadillo archive (download page). See also the Armadillo FAQ page. – Intensive

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

Distinguishing

Nice. Thank you And how can I convert Eigen to Armadillo – Squalene

Not a bad Q and a nice A. Addition to Gallery, maybe? – Flow

Aye. Probably this weekend with the 3x other posts I need to migrate... – Distinguishing

@DirkEddelbuettel: I'm here because a reference in R-SIG-Mac to coatless's blog posting describing how to make make an openMP-aware version of R. thecoatlessprofessor.com/programming/openmp-in-r-on-os-x – Frankfurter

@42 If I were you I'd use the rcpp-devel mailing list – Flow

OK. I'll subscribe. BTW, what's the Gallery? (Or perhaps this is already in the listinfo for rcpp-devel and I will see it in a couple of minutes?) – Frankfurter

It the widely-referenced contrubuted site over here with its source here at GH. – Flow

I'll admit that I probably should have googled it rather than question-commenting. First hit on the obvious strategy. So in penance, I'll now buy your book. – Frankfurter

@Distinguishing Thank you once for your answer. Your remark is very important If I intend to update the converted matrix. I found that if I do not want the update affects the first matrix, I can use arma(ptr_aux_mem, n_rows, n_cols, true, false) instead of arma(ptr_aux_mem, n_rows, n_cols, false, false). But that is for converting from Eigen to Armadillo. Do you know a way to convert from Armadillo to Eigen and allow update without changing Armadillo matrix? – Squalene

@Distinguishing I'm having an issue with this. The first element -i.e. at (0,0)- of the conversion from Eigen to Armadillo does not get converted correctly. I'm not sure whether I should post a new question though. – Dominate

@Dominate Just re-ran above and everything worked. Double check your compilation settings, package and R versions. – Distinguishing

pastebin.com/7AgTHPNf here's what I get. from R I also get it to work correctly – Dominate

Your code is creating a new copy on call instead of using a reference on call. Need to do a reference update. This is a separate question. – Distinguishing

That is, you need to use arma::mat matrixxd_to_armamat(Eigen::MatrixXd& eigen_A). – Distinguishing

@Distinguishing I'll post a question, I get the same result even if I use the reference – Dominate

@Distinguishing here: #58176522 – Dominate

Answer 2

I just spend a couple of hours trying to convert Eigen sparse matrix to Armadillo sparse matrix and I'll post the code here if someone else find a need to do the same.

I was doing this because I could not find an eigensolver for the sparse complex matrices, and Armadillo was the only one that had it, but the rest of my code was already done in Eigen so I had to do the conversion.

#include <Eigen/Sparse>
#include <armadillo>
using namespace std;
using namespace arma;
int main() {
    auto matrixA = new SparseMatrix<complex<double>>(numCols, numRows); //your Eigen matrix

    /*
    SOME CODE TO FILL THE Eeigen MATRIX

    */
    //  now create a separate vectors for row indeces, first non-zero column element indeces and non-zero values
    //  why long long unsigned int, because armadilo will expect that type when constructing sparse matrix
    vector<long long unsigned int> rowind_vect((*matrixA).innerIndexPtr(),
                                               (*matrixA).innerIndexPtr() + (*matrixA).nonZeros());
    vector<long long unsigned int> colptr_vect((*matrixA).outerIndexPtr(),
                                               (*matrixA).outerIndexPtr() + (*matrixA).outerSize() + 1);
    vector<complex<double>> values_vect((*matrixA).valuePtr(),
                                        (*matrixA).valuePtr() + (*matrixA).nonZeros());

    //  you can delete the original matrixA to free up space
    delete matrixA;

    //new Armadillo vectors from std::vector, we set the flag copy_aux_mem=false, so we don't copy the values again
    cx_dvec values(values_vect.data(), values_vect.size(), false);
    uvec rowind(rowind_vect.data(), rowind_vect.size(), false);
    uvec colptr(colptr_vect.data(), colptr_vect.size(), false);

    //  now create Armadillo matrix from these vectors
    sp_cx_dmat arma_hamiltonian(rowind, colptr, values, numCols, numRows);

    //  you can delete the vectors here if you like to free up the space

    return 0; 

}    

Consequence