c++ - 是否可以将矩阵列表传递给 vector> 内联？

Question

我正在尝试对列表格式的矩阵行进行运算，以计算我正在编写的某些软件的 Hessian 矩阵(混合偏导数矩阵)的部分。我发现我只能在 R 中如此快速地执行此操作(即使使用并行化)，因此切换到 Rcpp 以获得更快的速度，并切换到 RcppEigen 以提供高级矩阵运算。当我依赖 List 类型来表示从 R 传递的矩阵/vector 列表时，我的 Cpp 代码会随着列表长度(每个元素都是矩阵或 vector )的增加而大大减慢。我不确定到底是为什么，但这可能是因为动态大小的对象？我的问题是:我可以使用 RcppEigen 通过类似下面的方式将列表从 R 传递到标准模板库 (STL) 的 vector 容器中吗？

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

我想这样做的原因是因为我读到访问 vector 比访问列表快得多。但是，我可能误解了这一点，如果是这样，我深表歉意。

我们的想法是传入一个 vector 列表 (B2) 和一个矩阵列表 (A2)。在这些列表的每个索引中，我遍历 A2 的当前索引中的矩阵 (A) 的行和 B2 的当前索引中的 vector (b)，计算:

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

j 从 0 到 rows-1。我最终会得到一个列表，其大小等于该索引中矩阵的行数，然后移动到外部限制的下一个索引，等等。

这是我使用 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))

上面的工作，但随着我增加 A2 和 B2 的长度以接近我在实际应用程序中的长度，即超过 1000，Cpp 版本相对于 R 实现变慢并最终变慢.为了克服这个问题，我想到了尝试使用标准模板库 vector 格式。我不知道该怎么做，所以我想我会从简单开始，只传递一个列表，尝试转换为

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

然后发送回 R。这是我尝试过的:

## 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")             

我提前感谢大家的见解，如果这个问题没有为 StackOverflow 设计好，我深表歉意。我查看了以前的 StackOverflow 问题，阅读了发帖指南，并事先在谷歌上搜索了很长时间，试图找到我的问题的答案，但似乎我只是超出了已经提出的问题范围。我非常愿意进行任何必要的更改，以使该示例更具可重现性，并使我想做的事情更加清楚。我知道你们都很忙，我不想浪费你们的时间。

在 Dirk 的建议下，我尝试传递一个 List，然后定义一个 ListOf，鉴于其中对象的明确性质，这应该更快(如果这是错误的，请纠正我!)

下面是代码片段的样子:

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

这通过调用在我的机器上传递:

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

但是，访问这个 ListOf 中的元素似乎仍然存在问题。我知道这通常有效，因为我用数字 vector 测试了它，即

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)

如果我尝试对每个元素都是 MatrixXd 的 ListOf 做同样的事情，我似乎遇到了问题:

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

这会在尝试时产生错误:

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

我会继续尝试。我只想更新我现在的位置。此时此刻，在完整阅读 Rcpp 书(我即将用我的个人开发基金订购它)之前，这是我知道的将列表的每个元素作为 Eigen Map MatrixXd 传递的唯一方法。感谢您的宝贵时间!

Answer 1

你可以反过来做:使用标准 List (我们知道显然通过了)其中每个元素(无论如何都必须是 SEXP)通过本征 Map<MatrixXd> (我们也知道单独通过)。

所以我会从简单的开始，逐步变得复杂，直到它崩溃。