在R中交叉制表两个大型逻辑向量的最快方法

Question

对于两个逻辑矢量,x并且y,长度> 1E8的,什么是计算2×2交叉表格最快的方法？

我怀疑答案是用C/C++编写它,但我想知道R中是否有一些关于这个问题已经非常聪明,因为它并不罕见.

示例代码,对于300M条目(如果3E8太大,可以让N = 1E8;我选择的总大小不到2.5GB(2.4GB).我的目标密度为0.02,只是为了让它更有趣(可以使用稀疏向量,如果这有帮助,但类型转换可能需要时间).

set.seed(0)
N = 3E8
p = 0.02
x = sample(c(TRUE, FALSE), N, prob = c(p, 1-p), replace = TRUE)
y = sample(c(TRUE, FALSE), N, prob = c(p, 1-p), replace = TRUE)

一些明显的方法:

1. table
2. bigtabulate
3. 简单的逻辑运算(例如sum(x & y))
4. 矢量乘法(嘘)
5. data.table
6. 上面的一些,parallel从multicore包(或新parallel包)

我已经尝试了前三个选项(请参阅我的回答),但我觉得必须有更好更好的东西.

我觉得这table很慢. bigtabulate对于一对逻辑向量来说似乎有些过分.最后,进行vanilla逻辑运算看起来像一个kludge,它看了每个向量太多次(3X？7X？),更不用说它在处理期间填充了大量额外的内存,这是一个巨大的时间浪费.

向量乘法通常是一个坏主意,但是当向量稀疏时,可能会因为存储它而获得优势,然后使用向量乘法.

随意改变N和p,如果将展示的制表功能,任何有趣的行为.:)

---

更新1.我的第一个答案给出了三种天真方法的时间,这是相信table速度缓慢的基础.然而,要意识到的关键是"逻辑"方法效率极低.看看它在做什么:

• 4个逻辑向量运算
• 4种类型转换(逻辑到整数或FP - 用于sum)
• 4矢量求和
• 8个赋值(1表示逻辑运算,1表示求和)

不仅如此,它甚至没有编译或并行化.然而,它仍然打败了裤子table.请注意bigtabulate,使用额外的类型转换(1 * cbind...)仍然是节拍table.

更新2.为了避免任何人指出R支持逻辑向量NA,并且这将成为系统中这些交叉表的扳手(在大多数情况下都是如此),我应该指出我的向量来自is.na()或is.finite().:)我一直在调试NA和其他非限定值 - 最近他们一直很头疼.如果你不知道你的所有参赛作品是否都是NA,你可以考试一下any(is.na(yourVector))- 在你采纳本问答中出现的一些想法之前,这是明智之举.

---

更新3. Brandon Bertelsen在评论中提出了一个非常合理的问题:为什么在子样本(最初的集合,毕竟是样本;-))中使用如此多的数据可能足以创建交叉制表？不要过于偏向统计数据,但数据来自TRUE两个变量的观察非常罕见的情况.一个是数据异常的结果,另一个是由于代码中可能存在的错误(可能的错误,因为我们只看到计算结果 - 将变量x视为"垃圾输入",以及y"垃圾输出".结果,问题是由代码输出的问题是否仅仅在数据异常的情况下,还是有哪里好数据变坏一些其他情况？(这就是为什么我问一个问题关于停止时NaN,NA或者Inf是遇到了.)

这也解释了为什么我的例子的TRUE值概率很低; 这些确实发生在0.1%以下的时间.

这是否表明了不同的解决方案路径？是的:它表明我们可以使用两个指数(即TRUE每组中的位置)和计数集交叉点.我避免设置交叉点,因为我被Matlab烧了一段时间(是的,这是R,但是请跟我一起),这会先对一组的元素进行排序.(我含糊地回忆起复杂性更令人尴尬的是:喜欢O(n^2)而不是O(n log n).)

Answer 1

如果您在巨大的逻辑向量上进行大量操作,请查看位包.通过将布尔值存储为真正的1位布尔值,可以节省大量内存.

这无济于事table; 它实际上使它变得更糟,因为由于它的构造方式,位向量中有更多的唯一值.但它确实有助于逻辑比较.

# N <- 3e7
require(bit)
xb <- as.bit(x)
yb <- as.bit(y)
benchmark(replications = 1, order = "elapsed", 
    bit = {res <- func_logical(xb,yb)},
    logical = {res <- func_logical(x,y)}
)
#      test replications elapsed relative user.self sys.self user.child sys.child
# 1     bit            1   0.129  1.00000     0.132    0.000          0         0
# 2 logical            1   3.677 28.50388     2.684    0.928          0         0

Answer 2

下面是用于逻辑方法,结果table,以及bigtabulate,对于N = 3E8:

         test replications elapsed relative user.self sys.self
2     logical            1  23.861 1.000000     15.36     8.50
3 bigtabulate            1  36.477 1.528729     28.04     8.43
1       table            1 184.652 7.738653    150.61    33.99

在这种情况下,table是一场灾难.

为了比较,这里是N = 3E6:

         test replications elapsed relative user.self sys.self
2     logical            1   0.220 1.000000      0.14     0.08
3 bigtabulate            1   0.534 2.427273      0.45     0.08
1       table            1   1.956 8.890909      1.87     0.09

在这一点上,似乎编写自己的逻辑函数是最好的,即使这种情况滥用sum,并且多次检查每个逻辑向量.我还没有尝试编译函数,但这应该会产生更好的结果.

更新1如果我们给出bigtabulate已经是整数的值,即如果我们在1 * cbind(v1,v2)bigtabulate之外进行类型转换,那么N = 3E6倍数是1.80,而不是2.4.相对于"逻辑"方法的N = 3E8倍数仅为1.21,而不是1.53.

更新2

正如Joshua Ulrich指出的那样,转换为位向量是一项重大改进 - 我们正在分配和移动少量数据:R的逻辑向量每个条目消耗4个字节("为什么？",你可能会问...... 嗯,我不知道,但答案可能会出现在这里.),而一个位向量消耗,一个位,每个条目 - 即1/32的数据.因此,x消耗1.2e9字节,而xb(下面代码中的位版本)仅消耗3.75e7字节.

我已经删除了更新基准测试table的bigtabulate变化(N = 3e8).请注意,logicalB1假设数据已经是位向量,而logicalB2对于类型转换的惩罚是相同的操作.由于我的逻辑向量是对其他数据的操作结果,因此我没有从位向量开始的好处.尽管如此,要支付的罚款相对较小.["logical3"系列仅执行3次逻辑运算,然后进行减法运算.由于它是交叉制表,我们知道总数,正如DWin所说的那样.

        test replications elapsed  relative user.self sys.self
4 logical3B1            1   1.276  1.000000      1.11     0.17
2  logicalB1            1   1.768  1.385580      1.56     0.21
5 logical3B2            1   2.297  1.800157      2.15     0.14
3  logicalB2            1   2.782  2.180251      2.53     0.26
1    logical            1  22.953 17.988245     15.14     7.82

我们现在已经把这个加速到1.8-2.8秒,即使有很多严重的低效率.有毫无疑问,应该这样做公不到1秒钟可行的,与变化,包括一种或多种的:C代码,编译和多核处理.在所有3(或4)个不同的逻辑操作可以独立完成之后,即使这仍然是计算周期的浪费.

最相似的最佳挑战者,logical3B2比大约快80倍table.它比天真的逻辑操作快约10倍.它仍然有很大的改进空间.

这是产生上述代码的代码. 注:我建议注释掉一些操作或载体的,除非你有很多的RAM -创建x,x1以及xb与相应的沿y物体时,会占用的内存公平一点.

另外,请注意:我应该使用1L整数乘数bigtabulate,而不仅仅是1.在某些时候,我将重新运行此更改,并建议更改为使用该bigtabulate方法的任何人.

library(rbenchmark)
library(bigtabulate)
library(bit)

set.seed(0)
N <- 3E8
p <- 0.02

x <- sample(c(TRUE, FALSE), N, prob = c(p, 1-p), replace = TRUE)
y <- sample(c(TRUE, FALSE), N, prob = c(p, 1-p), replace = TRUE)

x1 <- 1*x
y1 <- 1*y

xb <- as.bit(x)
yb <- as.bit(y)

func_table  <- function(v1,v2){
    return(table(v1,v2))
}

func_logical  <- function(v1,v2){
    return(c(sum(v1 & v2), sum(v1 & !v2), sum(!v1 & v2), sum(!v1 & !v2)))
}

func_logicalB  <- function(v1,v2){
    v1B <- as.bit(v1)
    v2B <- as.bit(v2)
    return(c(sum(v1B & v2B), sum(v1B & !v2B), sum(!v1B & v2B), sum(!v1B & !v2B)))
}

func_bigtabulate    <- function(v1,v2){
    return(bigtabulate(1*cbind(v1,v2), ccols = c(1,2)))
}

func_bigtabulate2    <- function(v1,v2){
    return(bigtabulate(cbind(v1,v2), ccols = c(1,2)))
}

func_logical3   <- function(v1,v2){
    r1  <- sum(v1 & v2)
    r2  <- sum(v1 & !v2)
    r3  <- sum(!v1 & v2)
    r4  <- length(v1) - sum(c(r1, r2, r3))
    return(c(r1, r2, r3, r4))
}

func_logical3B   <- function(v1,v2){
    v1B <- as.bit(v1)
    v2B <- as.bit(v2)
    r1  <- sum(v1B & v2B)
    r2  <- sum(v1B & !v2B)
    r3  <- sum(!v1B & v2B)
    r4  <- length(v1) - sum(c(r1, r2, r3))
    return(c(r1, r2, r3, r4))
}

benchmark(replications = 1, order = "elapsed", 
    #table = {res <- func_table(x,y)},
    logical = {res <- func_logical(x,y)},
    logicalB1 = {res <- func_logical(xb,yb)},
    logicalB2 = {res <- func_logicalB(x,y)},

    logical3B1 = {res <- func_logical3(xb,yb)},
    logical3B2 = {res <- func_logical3B(x,y)}

    #bigtabulate = {res <- func_bigtabulate(x,y)},
    #bigtabulate2 = {res <- func_bigtabulate2(x1,y1)}
)