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

[Douglas Bates]

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.