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?
简单回答:这项工作最昂贵的部分是任务.因此,解决方案是坚持矩阵,并使用Rcpp和C++修改矩阵到位.下面有两个很好的答案和例子.[对于那些应用于其他问题的人,请务必阅读解决方案中的免责声明!].滚动至问题底部,了解更多经验教训.
这是我的第一个堆栈溢出问题——我非常感谢您抽出时间来查看,如果我遗漏了什么,我向您道歉.我正在开发一个R软件包,在这个软件包中,我有一个性能瓶颈,需要对矩阵的某些部分进行子集设置和写入(对于统计学家来说,需要注意的是,在处理每个数据点之后,应用程序正在更新足够的统计数据).单个操作的速度非常快,但数量之多要求它尽可能快.这个 idea 的最简单版本是一个维数为K×V的矩阵,其中K通常在5到1000之间,V可以在1000到1000000之间.
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的"更改时复制"语义,它可能会非常缓慢(但请参阅下面的经验教训,在某些情况下,修改实际上可能会发生).
对于我的问题,对象不需要是矩阵(我对这里列出的差异很敏感).我总是想要完整的列,所以数据帧的列表 struct 很好.我的解决方案是使用Matthew Dowle的惊人数据.桌子Package.使用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在注释中指出,一个选项是转换为三重表示(行号、列号、值)并使用数据.表二进制搜索以加快检索速度.我手动测试了它,虽然速度很快,但它的内存需求大约是矩阵的三倍.如果可能的话,我希望避免这种情况.
下面的问题是: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]]如果没有方法分派,就我所知,无法一次获取多个列,尽管它看起来确实是可能的.与前一个问题一样,似乎因为我们可以快速覆盖这些值,所以我们应该能够快速读取它们.
有什么建议吗?另外,请告诉我是否有比数据更好的解决方案.表来解决这个问题.我意识到它在很多方面都不是预期的用例,但我试图避免将整个系列的操作移植到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 指出,作业(job)中最昂贵的部分是分配过程.这两种解决方案都基于C代码给出答案,C代码就地修改矩阵,打破了R惯例,即不修改函数的参数,但提供了更快的结果.
哈德利·威克姆(Hadley Wickham)在 comments 中提到,只要对象mat没有在其他地方引用,矩阵修改实际上是在原地完成的(见http://adv-r.had.co.nz/memory.html#modification-in-place).这指向了一个有趣而微妙的点.我在RStudio上做判断.正如哈德利在书中指出的那样,RStudio为函数之外的每个对象创建了一个额外的引用.因此,在函数的性能情况下,修改会在适当的位置发生,而在命令行上,它会产生一种"更改时复制"的效果.Hadley的软件包pryr有一些很好的功能来跟踪内存的引用和地址.