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

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; 

}