[Rcpp-devel] Slices of 3d arrays as matrices

[Giovanni Petris]

Thank you Doug and Dirk.

I was afraid I had to move to Eigen... The problem now is that in the transition I am working with matrices both of type NumericMatrix and Map<MatrixXd>.

The following snippet is a loop that is giving me troubles...

   wz1 <- lapply(1:K, function(k) matrix(0.0, N, M)) # this is defined in R, but used by the C++ function

    Rcpp::NumericMatrix cpp_u(u);
    Rcpp::NumericMatrix cpp_mu(mu);
    Rcpp::NumericVector cpp_Sigma(Sigma);
    double *mm = cpp_Sigma.begin();
    int K = cpp_mu.ncol();
    int M = cpp_mu.nrow();
    int N = cpp_u.nrow();
    Mmat var1(NULL, M, M); 
    SelfAdjointEigenSolver<MatrixXd> var1_eigen(M);
    Environment glob = Environment::global_env(); 
    List cpp_wz1(glob.get("wz1"));
    NumericMatrix tmpMatrix1;

     for (int k=0; k < K; k++) {
        tmpMatrix1 = as<SEXP>(cpp_wz1[k]);
        for (int i=0; i < M; i++) 
            tmpMatrix1(_, i) = cpp_u(_, i) - cpp_mu(i, k);
        new (&var1)Mmat(mm + k * M * M, M, M);
        var1_eigen.compute(var1);
        tmpMatrix1 = tmpMatrix1 * var1_eigen.operatorInverseSqrt(); // <== here is the problem
       // more stuff here...
    }

Do you have any suggestion for an analog of the line 

  tmpMatrix1 = as<SEXP>(cpp_wz1[k]);

using an Eigen::Map<Eige::MatrixXd> object instead of the Rcpp::NumericMatrix object tmpMatrix1?

Thank you in advance for any insight.

Best,
Giovanni

________________________________________
From: dmbates at gmail.com [dmbates at gmail.com] on behalf of Douglas Bates [bates at stat.wisc.edu]
Sent: Tuesday, October 23, 2012 1:20 PM
To: Giovanni Petris
Cc: rcpp-devel at lists.r-forge.r-project.org
Subject: Re: [Rcpp-devel] Slices of 3d arrays as matrices

On Tue, Oct 23, 2012 at 12:56 PM, Giovanni Petris <gpetris at uark.edu> wrote:
> Hello,
>
> I have a 3d array and I want to compute the eigendecomposition of each slice [,, k]. I am trying to make Rcpp see the relevant slice of the 3d array as a matrix. My experiments, along the lines illustrated below, have so far been unsuccessful. Does anybody have a suggestion? Thank you in advance!

>    NumericVector cpp_Sigma(Sigma); // 3d array, with dim(Sigma) = c(M,M,K), as a vector
>    NumericMatrix tmpMatrix;
>    List var_eigen;
>    Function eigen("eigen");
>    int k=2; // 3rd slice, for example
>
>    tmpMatrix = cpp_Sigma[k * M * M];
>    tmpMatrix.attr("dim") = Dimension(M, M);
>    var_eigen = eigen(tmpMatrix);
>    return var_eigen;

To begin with, it is not a good idea to call R's eigen function from
within C++ code.  It would be much better to call code from the
RcppEigen or RcppArmadillo packages.  The eigen function in R is just
going to marshal its arguments and call a Lapack routine, which is
done more cleanly in RcppArmadillo.

Are the slices symmetric matrices?  If they are general matrices you
may need to allow for complex eigenvalues and eigenvectors.

Assuming that the slices are symmetric you could use something like

typedef Eigen::Map<Eigen::MatrixXd> Mmat;
const NumericVector cpp_Sigma(Sigma);
double *mm = cpp_Sigma.begin();
int k = 2;
const Mmat em(mm + k * M * M, M, M);
const Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> VLV(em);
return wrap(VLV.eigenvalues());

if you just want the eigenvalues.