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有一些很好的功能来跟踪内存的引用和地址.

推荐答案

Fun with 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级别重复调用函数,从而提高性能.

R相关问答推荐

使用lares::corr_var函数在for循环中分配变量的问题

在特定列上滞后n行,同时扩展框架的长度

更新合适的R mgcv::bam模型报告无效类型(关闭).'';错误

使用R的序列覆盖

如果列中存在相同的字符串,则对行值进行总和

更改编号列表的 colored颜色

在GGPLATE中将突出的点放在前面

在R中使用download. file().奇怪的URL?

在另存为PNG之前隐藏htmlwidget绘图元素

将饼图插入条形图

R spatstat Minkowski Sum()返回多个边界

函数可以跨多个列搜索多个字符串并创建二进制输出变量

根据另一列中的值和条件查找新列的值

警告消息";没有非缺失的参数到min;,正在返回数据中的inf";.表分组集

使用来自嵌套列和非嵌套列的输入的PURRR:MAP和dplyr::Mariate

按两个因素将观测值分组后计算单独的百分比

在r中整理图例和堆叠图的问题

长/纬点继续在堪萨斯-SF结束,整齐的人口普查

抽样变换-REXP与RWEIBUR

使用LAG和dplyr执行计算,以便按行和按组迭代