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

[Giovanni Petris]

Hi Doug,

Thank you for your answer and for your patience reading my confusing code. I have read the jss paper on RcppEigen and the Eigen tutorial. I am very excited about the opportunities that RcppEigen (and Rcpp) gives to the R programmer.  Thank you so much for the very useful contribution!

To get back quickly to the point in my post, the reason for accessing and modifying an R object from within C++ code is that the final application will involve large matrices and many iterations of a Gibbs sampler. Therefore, I would like to have some workspace allocated once and for all at the beginning of the calculations. In C I would probably declare a static double *workmat, invisible from R, but permanent for the different calls to my C function. 

Best,
Giovanni
________________________________________
From: dmbates at gmail.com [dmbates at gmail.com] on behalf of Douglas Bates [bates at stat.wisc.edu]
Sent: Wednesday, October 24, 2012 10:34 AM
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 5:52 PM, Giovanni Petris <gpetris at uark.edu> wrote:
> 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?

This code is very confusing and breaks the fundamental rule of
functional semantics in R.  It is very dangerous to reach out into the
global environment, grab the value of a symbol and change it.  R does
copy on write but this code changes the wz1 object without R being
aware that it needs to copy it.

In other words, you are going about things in an awkward way.  You
have the ability to create the object to be returned within the C++
code so this is a roundabout and dangerous way of creating the result.

To stay within the Eigen framework, I would create the array to be
returned as a matrix of dimension (M, M*K) then change the dimension
in R to (M, M, K) upon return.

There are several idioms for dealing with matrices passed from R as
Eigen mapped matrices in C++ as Dirk and I explained in our vignette
on RcppEigen.  It would be worthwhile reading that vignette if you
want to use RcppEigen, and also read at least the tutorial in the
Eigen3 documentation.

I would start with something like the enclosed but that code doesn't
compile and gives the usual page after page of reports when compiling
Eigen-based C++ code with g++.  Regrettably I don't have time to debug
it.

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