Is it possible to pass a list of Matrices to vector<Eigen::Map<Eigen::MatrixXd>> inline?

Question

I am attempting to conduct operations on rows of matrices in a list format to calculate portions of a Hessian matrix (mixed partial derivatives matrix) for some software I am writing. I found that I could only do this so fast in R (even with parallelization), and so have switched over to Rcpp for quicker speed and RcppEigen for the high level matrix operations provided. When I have relied on the List type for representing lists of matrices/vectors passed from R, my Cpp code slows down tremendously as the length of the lists (each element being a matrix or vector) increases. I am not sure exactly why, but it may be because of dynamically sized objects? My questions is: Can I pass a list from R into a vector container from the Standard Template Library (STL) using RcppEigen via something like the following?

vector<Eigen::Map<Eigen::MatrixXd>> A(as<vector<Eigen::Map<Eigen::MatrixXd>> >(AA)) 

The reason I want to do so is because I have read that accessing vectors is much faster than accessing lists. However, I may have misinterpreted this and I apologize if so.

The idea is to pass in a list of vectors (B2) and a list of matrices (A2). Within each index of these lists, I iterate over the rows of the matrix (A) in the current index of A2 and the vector (b) in the current index of B2, calculating:

b[j] * t(A[j,]) %*% A[j,]

for j from 0 to rows-1. I would end up with a list with size equal to the number of rows of the matrix in that index, then move on to the next index of the outer limit, etc.

Here is a reproducible example of what I have been able to do using List:

library(inline)
library(RcppEigen)
library(microbenchmark)    

## Create function which takes list into Rcpp and does all manipulations internally (no lapply outside)
A2 <- lapply(1:2, function(t) matrix(rnorm(10 * t), nrow = t, ncol = 10))
B2 <- lapply(1:2, function(t) rnorm(t))

## This becomes slower relative to R as the size increases.
## Something is not right in how I am programming this.
retLLMat <- "using Eigen::VectorXd;
         typedef Eigen::Map<Eigen::MatrixXd> MapMatd;
         typedef Eigen::Map<Eigen::VectorXd> MapVecd;
         List A(AA), B(BB);
         int listSize = A.size(), ncol, sublistSize;
         List outList;
         double sub;
         for (int i = 0; i < listSize; i++)
         {
           List subList;
           MapMatd subMat(as<MapMatd >(A[i]));
           MapVecd subVec(as<MapVecd >(B[i]));
           ncol = subMat.cols();
           VectorXd currRow(ncol);
           sublistSize = subMat.rows();
           for (int j = 0; j < sublistSize; j++)
           {
             currRow = subMat.row(j);
             sub = subVec[j];
             subList[String(j)] = (sub * currRow) * currRow.transpose();
           }
           outList[String(i)] = subList; 
         }
         return wrap(outList);"

## Compile Cpp code
retLLMatC <- cxxfunction(signature(AA = "List", BB = "List"), retLLMat, plugin = "RcppEigen")
## R version
retLLMat <- function(A, B) mapply(function(a, b) mapply(function(a, b) b * a, lapply(apply(a, 1, function(t) list(tcrossprod(t))), "[[", 1), b, SIMPLIFY = FALSE), A, B, SIMPLIFY = FALSE)

## Test R vs Rcpp version
microbenchmark(retLLMat(A2, B2), retLLMatC(A2, B2))

The above works, but as I increase the length of A2 and B2 to be closer to what I have in my actual real application, which is more than 1000, the Cpp version slows down relative to the R implementation and eventually is slower. To overcome this, I thought of trying to use the standard template library vector format. I didn't know how to do this, so I figured I'd start simple, and just pass a List, try to convert to

vector<Eigen::Map<Eigen::MatrixXd>> 

and then send back to R. This is what I tried:

## Testing using a vector of Map<MatrixXd>
## Simplified by trying to return the list after reading it in
VectorMat <- "using Eigen::MatrixXd;
              using std::vector;
              typedef Eigen::Map<Eigen::MatrixXd> MapMatd;
              vector<MapMatd> A(as<vector<MapMatd> >(AA);
              return wrap(A);"

## This produces errors 
test <- cxxfunction(signature(AA = "List"), VectorMat, plugin = "RcppEigen")             

I thank all in advance for insights, and sincerely apologize if this question is not designed well for StackOverflow. I looked around at previous StackOverflow questions, read the posting guide, and googled beforehand for quite a while to try to find an answer to my problem, but it seemed I was just outside of the scope of what had been already asked. I am more than willing to make any changes necessary to make this example more reproducible and also to make what I want to do more clear. I know that you all are very busy, and I do not want to waste your time.

At Dirk's suggestion, I tried passing a List, then defining a ListOf, which should be faster given the explicit nature of the objects within (please correct me if this is wrong!)

Here's what the snippet then looked like:

ListMat1 <- "using Eigen::MatrixXd;
typedef Eigen::Map<Eigen::MatrixXd> MapMatd;
ListOf<MapMatd> A(as<ListOf<MapMatd> >(AA));
return wrap(A);"

This passed on my machine by calling:

ListMat <- cxxfunction(signature(AA = "List"), ListMat1, plugin = "RcppEigen")
res <- ListMat(A2)

However, it seems to still be problematic to access elements within this ListOf. I know this normally works, because I tested it with a numeric vector, i.e.

b2 <- lapply(1:5, function(t) numeric(t))

vecTest <- "ListOf<NumericVector> b(as<ListOf<NumericVector> >(bb)); 
NumericVector res = b[0];
return wrap(res);"   

vecTestfn <- cxxfunction(signature(bb = "List"), vecTest, plugin = "RcppEigen")
vecTestfn(b2)

If I try to do the same thing with ListOf where each element is a MatrixXd, I seem to have an issue:

ListMatInd <- "using Eigen::MatrixXd;
typedef Eigen::Map<Eigen::MatrixXd> MapMatd;
ListOf<MapMatd> A(as<ListOf<MapMatd> >(AA));
MatrixXd res = A[0];
return wrap(res);"

This creates an error upon trying:

ListMatIndfn <- cxxfunction(signature(AA = "List"), ListMatInd, plugin = "RcppEigen")

I will continue trying things. I just want to update regarding where I am right now. At this moment, before having read the Rcpp book in full (I am about to order it with my personal development fund), this is the only way I knew to pass each element of the list as a Eigen Map MatrixXd. Thank you for your time!

Answer 1

You may be able to do it the other way around: use a standard List (which we know passes, obviously) where each element (which has to be SEXP anyway) passes an Eigen Map<MatrixXd> (which we also know pass individually).

So I would start simple and make it more complicated step by step til it breaks.