r - 如何优化 R 中矩阵子部分的读写(可能使用 data.table)

Question

TL;博士

What is the fastest method in R for reading and writing a subset of columns from a very large matrix. I attempt a solution with data.table but need a fast way to extract a sequence of columns?

Short Answer: The expensive part of the operation is assignment. Thus the solution is to stick with a matrix and use Rcpp and C++ to modify the matrix in place. There are two excellent answers below with examples.[for those applying to other problems be sure to read the disclaimers in the solutions!]. Scroll to the bottom of the question for some more lessons learned.

这是我的第一个 Stack Overflow 问题 - 我非常感谢您抽出时间来查看，如果我遗漏了任何内容，我深表歉意。我正在研究一个 R 包，我在子集和写入矩阵的部分时遇到了性能瓶颈(对于统计学家来说，应用程序在处理每个数据点后更新足够的统计数据)。单个操作的速度非常快，但它们的绝对数量要求它尽可能快。该想法的最简单版本是维度为 K 乘以 V 的矩阵，其中 K 通常介于 5 和 1000 之间，而 V 可以介于 1000 和 1,000,000 之间。

set.seed(94253)
K <- 100
V <- 100000
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
library(microbenchmark)
microbenchmark(mat[,Vsub] <- mat[,Vsub] + toinsert)

因为这样做了很多次，由于 R 的更改时复制语义，它可能会非常慢(但请参阅下面的经验教训，在某些情况下实际上可以进行修改)。

对于我的问题，对象不必是矩阵(并且我对此处概述的差异很敏感 Assign a matrix to a subset of a data.table )。我总是想要完整的列，因此数据框的列表结构很好。我的解决方案是使用 Matthew Dowle 很棒的 data.table 包。使用 set() 可以非常快速地完成写入。不幸的是，获得值(value)有点复杂。我们必须使用=FALSE 调用变量设置，这会大大减慢速度。

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 在评论中指出，一种选择是转换为三重表示(行号、列号、值)并使用 data.tables 二进制搜索来加速检索。我手动测试了它，虽然它很快，但它确实使矩阵的内存需求大约增加了三倍。如果可能的话，我想避免这种情况。

跟随这里的问题: Time in getting single elemets from data.table and data.frame objects . Hadley Wickham 为单个索引提供了令人难以置信的快速解决方案

Vone <- Vsub[1]
toinsert.one <- toinsert[,1]
set(DT, i=NULL, j=Vone,(.subset2(DT, Vone) + toinsert.one))

但是，由于 .subset2(DT,i) 只是 DT[[i]] 没有方法分派(dispatch)，因此(据我所知)无法一次抓取几列，尽管它看起来确实应该是可能的。与上一个问题一样，看起来既然我们可以快速覆盖这些值，我们应该能够快速读取它们。

有什么建议么？另外请让我知道是否有比 data.table 更好的解决方案来解决这个问题。我意识到它在许多方面并不是真正的预期用例，但我试图避免将整个系列的操作移植到 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

答案和经验教训:

Comments identified the most expensive part of the operation as the assignment process. Both solutions give answers based on C code which modify the matrix in place breaking R convention of not modifying the argument to a function but providing a much faster result.

Hadley Wickham stopped by in the comments to note that the matrix modification is actually done in place as long as the object mat is not referenced elsewhere (see http://adv-r.had.co.nz/memory.html#modification-in-place). This points to an interesting and subtle point. I was performing my evaluations in RStudio. RStudio as Hadley notes in his book creates an additional reference for every object not within a function. Thus while in the performance case of a function the modification would happen in place, at the command line it was producing a copy-on-change effect. Hadley's package pryr has some nice functions for tracking references and addresses of memory.

Answer 1

Rcpp 的乐趣:

您可以使用Eigen's Map class就地修改 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 的功能范式的警告适用。

我在 C++ 中进行了加法，但是将函数更改为仅进行赋值是微不足道的。

显然，如果您可以在 Rcpp 中实现整个循环，您就可以避免在 R 级别重复调用函数并获得性能。