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

[Douglas Bates]

On Wed, Oct 24, 2012 at 11:14 AM, Giovanni Petris <gpetris at uark.edu> wrote:
> 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.

I know the problem.  A lot of the effort in the RcppEigen-based C++
code for the lme4 package is designed to allow for in-place updating
of large structures.

One approach that is legitimate in R, but only because of special
semantics for environments, is to make a new environment and install
the objects that you want to modify in there.  Environments behave
differently from other R objects in that there is only ever one copy
of an environment so changes in the objects in an environment made in
one place are reflected in all other views of the environment.

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