如何优化对R中矩阵子节的读写(可能使用data.table)

Question

TL; DR

R中用于从非常大的矩阵读取和写入列子集的最快方法是什么.我尝试使用data.table的解决方案,但需要一种快速的方法来提取一系列列？

简答:操作的昂贵部分是分配.因此,解决方案是坚持使用矩阵并使用Rcpp和C++来修改矩阵.下面有两个很好的答案和例子.[适用于其他问题的人一定要阅读解决方案中的免责声明!].滚动到问题的底部,了解更多经验教训.

---

这是我的第一个Stack Overflow问题 - 我非常感谢您抽出时间看看,如果我遗漏了任何内容,我会深表歉意.我正在研究一个R软件包,其中我有一个性能瓶颈,从子集化和写入矩阵的部分(统计学家的NB,应用程序在处理每个数据点后更新足够的统计数据).单个操作速度非常快,但是它们的数量绝对要求它尽可能快.该概念的最简单版本是尺寸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的复制变更语义,它可能会很慢(但是参见下面的经验教训,修改实际上可以在某些情况下发生).

对于我的问题,对象不必是一个矩阵(我对这里概述的差异很敏感.将一个矩阵分配给data.table的子集).我总是想要完整的列,因此数据框的列表结构很好.我的解决方案是使用Matthew Dowle的令人敬畏的data.table包.使用set()可以非常快速地完成写入.不幸的是,获得价值有点复杂.我们必须使用= 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二进制搜索来加速检索.我手动测试了它,虽然速度很快,但它大约是矩阵内存需求的三倍.如果可能的话,我想避免这种情况.

在这里提出问题:从data.table和data.frame对象获取单个元素的时间.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]],所以没有办法(据我所知)一次抓取几个列,尽管看起来它应该是可能的.与前一个问题一样,似乎因为我们可以快速覆盖这些值,所以我们应该能够快速读取它们.

有什么建议？如果这个问题有一个比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

---

答案和经验教训:

评论确定了作为分配过程的最昂贵的操作部分.两种解决方案都基于C代码给出答案,C代码修改了矩阵,破坏了R约定,即不修改函数的参数,但提供更快的结果.

Hadley Wickham在评论中停下来指出,只要对象垫没有在别处引用,矩阵修改实际上已经完成(参见http://adv-r.had.co.nz/memory.html#modification-就地).这指出了一个有趣而微妙的观点.我在RStudio进行评估.Hadley在他的书中指出的RStudio为不在函数内的每个对象创建了一个额外的引用.因此,在函数的性能情况下,修改将在适当的位置发生,在命令行,它产生了复制变换效果.Hadley的包pryr有一些很好的功能来跟踪引用和内存地址.

Answer 1

有趣的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的功能范式的警告适用.

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

显然,如果你可以在Rcpp中实现你的整个循环,你可以避免在R级重复的函数调用,并获得性能.

Answer 2

这就是我的想法.对于Rcpp和朋友来说,这可能会更加性感,但我对这些工具并不熟悉.

#include <R.h>
#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)

这种方法相对于Roland的优势在于,这只能完成任务.他的功能为您添加了更快的功能,但也意味着您需要为您需要执行的每项操作提供单独的C/C++功能.