如何优化读取和写入R中的matrix的子节（可能使用data.table）
TL; DR
R中用于读取和写入非常大matrix列的子集的最快方法是什么？ 我试图用data.table解决scheme，但需要一个快速的方式来提取一系列的列？
简答：操作的昂贵部分是分配。 因此，解决scheme是坚持使用matrix，并使用Rcpp和C ++来修改matrix。 下面有两个很好的答案，其中有些例子适用于其他问题，请务必阅读解决scheme中的免责声明！ 滚动到问题的底部，了解更多的经验教训。
这是我的第一个堆栈溢出问题，我非常感谢你的时间在看，我很抱歉，如果我什么都没有留下。 我正在研究一个R包，其中我有一个性能瓶颈，从子集化和写入到matrix的某些部分（对统计学家来说，应用程序在处理每个数据点后更新足够的统计信息）。 单独的操作非常快，但它们的数量要求尽可能快。 这个想法的最简单版本是一个维度K乘以V的matrix，其中K一般在5到1000之间，V可以在1000到1000000之间。
set.seed(94253) K <- 100 V <- 100000 mat <- matrix(runif(K*V),nrow=K,ncol=V) 然后我们结束对列的一个子集进行计算并将其添加到完整的matrix中。 因此天真的看起来像
Vsub <- sample(1:V, 20) toinsert <- matrix(runif(K*length(Vsub)), nrow=K, ncol=length(Vsub)) mat[,Vsub] <- mat[,Vsub] + toinsert library(microbenchmark) microbenchmark(mat[,Vsub] <- mat[,Vsub] + toinsert)
因为这样做很多次，所以R的copy-on-change语义会很慢（但是看到下面的经验教训，在一些情况下实际上可以进行修改）。
对于我的问题，该对象不需要是matrix（我对这里所描述的区别非常敏感） 将一个matrix分配给data.table的一个子集 。 我总是想要整列，所以数据框的列表结构是好的。 我的解决scheme是使用Matthew Dowle的真棒data.table包。 写入可以使用set（）非常快速地完成。 不幸的是，获得价值有点复杂。 我们必须使用= FALSE来调用variables设置，这会显着降低速度。
library(data.table) DT <- as.data.table(mat) set(DT, i=NULL, j=Vsub,DT[,Vsub,with=FALSE] + as.numeric(toinsert))
在set（）函数中，使用i = NULL来引用所有行的速度非常快，但是（可能是因为事物被存储在引擎盖下），j没有可比的选项。 @Roland在评论中指出，一种select是转换为三重表示（行号，列号，值）并使用data.tables二进制search来加速检索。 我手动testing了这一点，虽然速度很快，但它的内存需求大约是matrix的三倍。 如果可能，我想避免这种情况。
在这里的问题之后： 从data.table和data.frame对象获取单个元素的时间 。 哈德利·韦翰（Hadley Wickham）为单一指数提供了一个非常快速的解决scheme
Vone <- Vsub[1] toinsert.one <- toinsert[,1] set(DT, i=NULL, j=Vone,(.subset2(DT, Vone) + toinsert.one))
但是因为.subset2（DT，i）只是DT [[i]]而没有方法调度，所以根据我的理解，没有办法一次抓取多个列，尽pipe它看起来应该是可能的。 就像上一个问题一样，我们可以快速覆盖这些值，所以我们应该能够快速读取这些值。
有什么build议么？ 也请让我知道，如果有比这个问题data.table更好的解决scheme。 我在很多方面意识到它并不是真正意义上的用例，但我试图避免将整个系列的操作移植到C中。
下面是讨论的元素的时序 – 前两个是所有列，后两个只是一列。
microbenchmark(mat[,Vsub] <- mat[,Vsub] + toinsert, set(DT, i=NULL, j=Vsub,DT[,Vsub,with=FALSE] + as.numeric(toinsert)), mat[,Vone] <- mat[,Vone] + toinsert.one, set(DT, i=NULL, j=Vone,(.subset2(DT, Vone) + toinsert.one)), times=1000L) Unit: microseconds expr min lq median uq max neval Matrix 51.970 53.895 61.754 77.313 135.698 1000 Data.Table 4751.982 4962.426 5087.376 5256.597 23710.826 1000 Matrix Single Col 8.021 9.304 10.427 19.570 55303.659 1000 Data.Table Single Col 6.737 7.700 9.304 11.549 89.824 1000
答案和经验教训：
评论将操作中最昂贵的部分确定为分配过程。 这两个解决scheme都给出了基于C代码的答案，这些代码修改了matrix，而不是修改函数的参数，而是提供了更快的结果。
哈德利·威克姆（Hadley Wickham）在评论中停下来注意到，只要对象垫没有在别处被引用，matrix修改实际上已经完成了（参见http://adv-r.had.co.nz/memory.html#modification-就地 ）。 这指出了一个有趣和微妙的一点。 我正在RStudio中执行我的评估。 作为哈德利在他的书中注意到的RStudio为不在函数内的每个对象创build一个额外的参考。 因此，在一个函数的执行情况下，修改会发生，在命令行它会产生一个“复制时间变化”的效果。 哈德利的软件包pryr有一些很好的function来跟踪内存的引用和地址。
乐趣与Rcpp：
你可以使用Eigen的Map类来修改一个R对象。
library(RcppEigen) library(inline) incl <- ' using Eigen::Map; using Eigen::MatrixXd; using Eigen::VectorXi; typedef Map<MatrixXd> MapMatd; typedef Map<VectorXi> MapVeci; ' body <- ' MapMatd A(as<MapMatd>(AA)); const MapMatd B(as<MapMatd>(BB)); const MapVeci ix(as<MapVeci>(ind)); const int mB(B.cols()); for (int i = 0; i < mB; ++i) { A.col(ix.coeff(i)-1) += B.col(i); } ' funRcpp <- cxxfunction(signature(AA = "matrix", BB ="matrix", ind = "integer"), body, "RcppEigen", incl) set.seed(94253) K <- 100 V <- 100000 mat2 <- mat <- matrix(runif(K*V),nrow=K,ncol=V) Vsub <- sample(1:V, 20) toinsert <- matrix(runif(K*length(Vsub)), nrow=K, ncol=length(Vsub)) mat[,Vsub] <- mat[,Vsub] + toinsert invisible(funRcpp(mat2, toinsert, Vsub)) all.equal(mat, mat2) #[1] TRUE library(microbenchmark) microbenchmark(mat[,Vsub] <- mat[,Vsub] + toinsert, funRcpp(mat2, toinsert, Vsub)) # Unit: microseconds # expr min lq median uq max neval # mat[, Vsub] <- mat[, Vsub] + toinsert 49.273 49.628 50.3250 50.8075 20020.400 100 # funRcpp(mat2, toinsert, Vsub) 6.450 6.805 7.6605 7.9215 25.914 100
我想这基本上就是@Joshua Ulrich提出的。 他关于打破R的function范式的警告适用。
我在C ++中进行了添加，但将函数更改为仅执行赋值操作是微不足道的。
显然，如果你可以在Rcpp中实现你的整个循环，你可以避免在R级重复函数调用，并且会获得性能。
这是我的想法。 这可能是Rcpp和朋友更性感，但我不熟悉这些工具。
#include <Rh> #include <Rinternals.h> #include <Rdefines.h> SEXP addCol(SEXP mat, SEXP loc, SEXP matAdd) { int i, nr = nrows(mat), nc = ncols(matAdd), ll = length(loc); if(ll != nc) error("length(loc) must equal ncol(matAdd)"); if(TYPEOF(mat) != TYPEOF(matAdd)) error("mat and matAdd must be the same type"); if(nr != nrows(matAdd)) error("mat and matAdd must have the same number of rows"); if(TYPEOF(loc) != INTSXP) error("loc must be integer"); int *iloc = INTEGER(loc); switch(TYPEOF(mat)) { case REALSXP: for(i=0; i < ll; i++) memcpy(&(REAL(mat)[(iloc[i]-1)*nr]), &(REAL(matAdd)[i*nr]), nr*sizeof(double)); break; case INTSXP: for(i=0; i < ll; i++) memcpy(&(INTEGER(mat)[(iloc[i]-1)*nr]), &(INTEGER(matAdd)[i*nr]), nr*sizeof(int)); break; default: error("unsupported type"); } return R_NilValue; }
把上面的函数放到addCol.c ，然后运行R CMD SHLIB addCol.c。 然后在R：
addColC <- dyn.load("addCol.so")$addCol .Call(addColC, mat, Vsub, mat[,Vsub]+toinsert)
这种方法对罗兰德的轻微优势是，这只是做这个任务。 他的function为你增加了更快，但也意味着你需要一个单独的C / C ++函数为你需要做的每一个操作。
