提问者:小点点

更快版本的梳理


有没有办法加快compn命令以获取从向量中获取的2个元素的所有唯一组合?

通常会这样设置:

# Get latest version of data.table
library(devtools)
install_github("Rdatatable/data.table",  build_vignettes = FALSE)  
library(data.table)

# Toy data
d <- data.table(id=as.character(paste0("A", 10001:15000))) 

# Transform data 
system.time({
d.1 <- as.data.table(t(combn(d$id, 2)))
})

但是,comn比使用data. table计算所有可能的组合慢10倍(在我的计算机上为23秒而不是3秒)。

system.time({
d.2 <- d[, list(neighbor=d$id[-which(d$id==id)]), by=c("id")]
})

处理非常大的向量,我正在寻找一种方法来节省内存,只计算唯一的组合(如comn),但以data. table的速度(参见第二个代码片段)。

我很感激任何帮助。


共3个答案

匿名用户

这是一种使用data. table函数foveraps()的方法,结果也很快!

require(data.table) ## 1.9.4+
d[, `:=`(id1 = 1L, id2 = .I)] ## add interval columns for overlaps
setkey(d, id1, id2)

system.time(olaps <- foverlaps(d, d, type="within", which=TRUE)[xid != yid])
#  0.603   0.062   0.717

请注意,foveraps()并不计算所有的排列。需要子集xid!=yid来删除self重叠。通过实现忽略自己参数,可以更有效地在内部处理子集-类似于IRanges::findOverlaps

现在只需使用获得的id执行子集:

system.time(ans <- setDT(list(d$id[olaps$xid], d$id[olaps$yid])))
#   0.576   0.047   0.662 

所以完全,~1.4秒。

优点是,即使您的data. tabled有超过1列您必须在其中获取组合,并且使用相同数量的内存(因为我们返回索引),您也可以以相同的方式执行此操作。在这种情况下,您只需执行以下操作:

cbind(d[olaps$xid, ..your_cols], d[olaps$yid, ..your_cols])

但它仅限于替换组合(.,2L)。不超过2L。

匿名用户

您可以使用gRbase中的compnPrim

source("http://bioconductor.org/biocLite.R")
biocLite("gRbase") # will install dependent packages automatically.
system.time({
 d.1 <- as.data.table(t(combn(d$id, 2)))
 })
#   user  system elapsed 
# 27.322   0.585  27.674 

system.time({
d.2 <- as.data.table(t(combnPrim(d$id,2)))
 })
#   user  system elapsed 
#  2.317   0.110   2.425 

identical(d.1[order(V1, V2),], d.2[order(V1,V2),])
#[1] TRUE

匿名用户

标题中带有快速一词的任何变体的帖子在没有基准的情况下都是不完整的。在我们发布任何基准之前,我只想提一下,自从这个问题发布以来,已经为R发布了两个高度优化的包,排列RcppAlgos(我是作者),用于生成组合。请注意,从RcppAlgos的版本2.3.0开始,我们可以利用多个线程来提高效率。

为了让您了解它们在ombnPrim和gRbase::c速度,这里有一个基本的基准:

## We test generating just over 3 million combinations
choose(25, 10)
[1] 3268760

microbenchmark(arrngmnt = arrangements::combinations(25, 10),
               combn = combn(25, 10),
               gRBase = gRbase::combnPrim(25, 10),
               serAlgos = RcppAlgos::comboGeneral(25, 10),
               parAlgos = RcppAlgos::comboGeneral(25, 10, nThreads = 4),
               unit = "relative", times = 20)
Unit: relative
    expr        min         lq       mean     median         uq        max neval
arrngmnt   2.979378   3.072319   1.898390   3.756307   2.139258  0.4842967    20
   combn 226.470755 230.410716 118.157110 232.905393 125.718512 17.7778585    20
  gRBase  34.219914  34.209820  18.789954  34.218320  19.934485  3.6455493    20
serAlgos   2.836651   3.078791   2.458645   3.703929   2.231475  1.1652445    20
parAlgos   1.000000   1.000000   1.000000   1.000000   1.000000  1.0000000    20

现在,我们针对生成组合选择2和生成data. table对象的非常具体的情况对发布的其他函数进行基准测试。

职能如下:

funAkraf <- function(d) {
    a <- comb2.int(length(d$id))      ## comb2.int from the answer given by @akraf
    setDT(list(V1 = d$id[a[,1]], V2 = d$id[a[,2]]))
}

funAnirban <- function(d) {
    indices <- combi2inds(d$id)
    ans2 <- setDT(list(d$id[indices$xid], d$id[indices$yid]))
    ans2
}

funArun <- function(d) {
    d[, `:=`(id1 = 1L, id2 = .I)] ## add interval columns for overlaps
    setkey(d, id1, id2)
    olaps <- foverlaps(d, d, type="within", which=TRUE)[xid != yid]
    ans <- setDT(list(d$id[olaps$xid], d$id[olaps$yid]))
    ans
}

funArrangements <- function(d) {
  a <- arrangements::combinations(x = d$id, k = 2)
  setDT(list(a[, 1], a[, 2]))
}

funGRbase <- function(d) {
  a <- gRbase::combnPrim(d$id,2)
  setDT(list(a[1, ], a[2, ]))
}

funOPCombn <- function(d) {
  a <- combn(d$id, 2)
  setDT(list(a[1, ], a[2, ]))
}

funRcppAlgos <- function(d) {
  a <- RcppAlgos::comboGeneral(d$id, 2, nThreads = 4)
  setDT(list(a[, 1], a[, 2]))
}

以下是OP给出的示例的基准:

d <- data.table(id=as.character(paste0("A", 10001:15000))) 

microbenchmark(funAkraf(d),
               funAnirban(d),
               funArrangements(d),
               funArun(d),
               funGRbase(d),
               funOPCombn(d),
               funRcppAlgos(d),
               times = 10, unit = "relative")
    Unit: relative
              expr       min        lq      mean    median        uq       max neval
       funAkraf(d)  3.220550  2.971264  2.815023  2.665616  2.344018  3.383673    10
     funAnirban(d)  1.000000  1.000000  1.000000  1.000000  1.000000  1.000000    10
funArrangements(d)  1.464730  1.689231  1.834650  1.960233  1.932361  1.693305    10
        funArun(d)  3.256889  2.908075  2.634831  2.729180  2.432277  2.193849    10
      funGRbase(d)  3.513847  3.340637  3.327845  3.196399  3.291480  3.129362    10
     funOPCombn(d) 30.310469 26.255374 21.656376 22.386270 18.527904 15.626261    10
   funRcppAlgos(d)  1.676808  1.956696  1.943773  2.085968  1.949133  1.804180    10

我们看到@AnirbanMukherjee提供的函数是这个任务最快的,其次是RcppAlgos/排列。对于这个任务,nThread没有任何影响,因为传递的向量是一个字符,这不是线程安全的。如果我们改为将id转换为因子呢?

dFac <- d
dFac$id <- as.factor(dFac$id)

library(microbenchmark)
microbenchmark(funAkraf(dFac),
               funAnirban(dFac),
               funArrangements(dFac),
               funArun(dFac),
               funGRbase(dFac),
               funOPCombn(dFac),
               funRcppAlgos(dFac),
               times = 10, unit = "relative")
Unit: relative
                 expr        min         lq      mean   median        uq       max   neval
       funAkraf(dFac)  10.898202  10.949896  7.589814 10.01369  8.050005  5.557014      10
     funAnirban(dFac)   3.104212   3.337344  2.317024  3.00254  2.471887  1.530978      10
funArrangements(dFac)   2.054116   2.058768  1.858268  1.94507  2.797956  1.691875      10
        funArun(dFac)  10.646680  12.905119  7.703085 11.50311  8.410893  3.802155      10
      funGRbase(dFac)  16.523356  21.609917 12.991400 19.73776 13.599870  6.498135      10
     funOPCombn(dFac) 108.301876 108.753085 64.338478 95.56197 65.494335 28.183104      10
   funRcppAlgos(dFac)   1.000000   1.000000  1.000000  1.00000  1.000000  1.000000      10 

现在,我们看到RcppAlgos比任何其他解决方案快大约2倍。特别是,RcppAlgos解决方案比Anirban给出的以前最快的解决方案快大约3倍。应该注意的是,这种效率的提高是可能的,因为因子变量实际上是底层的整数,带有一些额外的属性

它们也都会给出相同的结果。唯一需要注意的是gRbase解决方案不支持因子。也就是说,如果您传递因子,它将被转换为字符。因此,如果您传递dFac,除了gRbase解决方案之外,所有解决方案都会给出相同的结果:

identical(funAkraf(d), funOPCombn(d))
#[1] TRUE
identical(funAkraf(d), funArrangements(d))
#[1] TRUE
identical(funRcppAlgos(d), funArrangements(d))
#[1] TRUE
identical(funRcppAlgos(d), funAnirban(d))
#[1] TRUE
identical(funRcppAlgos(d), funArun(d))
#[1] TRUE

## different order... we must sort
identical(funRcppAlgos(d), funGRbase(d))
[1] FALSE
d1 <- funGRbase(d)
d2 <- funRcppAlgos(d)

## now it's the same
identical(d1[order(V1, V2),], d2[order(V1,V2),])
#[1] TRUE

感谢@Frank指出如何比较两个data. table,而无需经历创建新的data.table然后排列它们的痛苦:

fsetequal(funRcppAlgos(d), funGRbase(d))
[1] TRUE