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

Question

对于长度大于1E8的两个逻辑向量x和y，计算2x2交叉表最快的方法是什么？

我猜想答案是用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. 上述方法之一，结合multicore包中的parallel（或新的parallel包）

我已经尝试了前三种选项（请参见我的答案），但我觉得一定有更好更快的方法。

我发现table的速度非常慢。对于一对逻辑向量来说，bigtabulate似乎有点过头了。最后，使用普通的逻辑操作似乎有点笨拙，并且它对每个向量进行了太多次的操作（3倍？7倍？），更不用说在处理过程中填充了大量的额外内存，这是一个巨大的时间浪费。

向量乘法通常不是一个好主意，但当向量是稀疏的时候，将其存储为稀疏向量，然后使用向量乘法可能会有优势。

---

为什么要使用这么多数据，当一个子样本可能足以用于创建交叉表呢？数据来自于两个变量中真实观测非常罕见的情况。其中一个是数据异常的结果，另一个可能是代码中的一个错误（可能是因为我们只看到了计算结果 - 将变量x视为“垃圾输入”，变量y视为“垃圾输出”）。因此，问题是输出中由代码引起的问题是否仅限于数据异常的情况，还是还有其他一些良好数据变坏的情况？（这就是为什么我问了一个关于遇到NaN、NA或Inf时停止的问题的原因。）

这也解释了为什么我的示例中TRUE值的概率很低；这些情况真正发生的时间少于0.1%。

这是否暗示了一条不同的解决方案路径？是的：它暗示我们可以使用两个索引（即每个集合中TRUE的位置）并计算集合的交集。我避免使用集合的交集，是因为我之前在Matlab中遇到过问题，它在执行交集操作之前会先对集合元素进行排序。（我模糊地记得复杂度甚至更令人尴尬：像O(n^2)而不是O(n log n)。）

那么我该如何在R中实现呢？

- Iterator

---

我很困惑为什么table在你看来速度很慢。在我使用它时，它总是很快的。（尽管在你的任务上花费了5分钟时间。） - IRTFM

是的，我也很惊讶。我的逻辑向量版本是sum(x>y)，sum(x<y)，sum(x&y)，然后从总数中减去这些。速度快多了。 - IRTFM

@DWin 没错 - 第四个操作可以通过从向量的长度中减去来简单地移除。尽管如此，考虑到我们知道我们将要进行制表操作，这仍然是多余的两个逻辑操作。这种简单性令人困扰。 - Iterator

@DWin 嗯，我们现在比“table()”有大约80倍的加速。 :) 我现在相信，至少还可以再提高10倍。 - Iterator

我知道这是一个关于大数据编程的问题...但是如果我建议做出一些假设并从向量中进行抽样，并接受一定程度的置信度而不是试图完美无缺，那么这是否会让人感到荒谬？例如，像table(sample(x,3e6),sample(y,3e6))这样的东西。 - Brandon Bertelsen

Answer 1

如果您需要对大量的逻辑向量进行操作，请查看bit包。它通过将布尔值存储为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

- Joshua Ulrich

---

谢谢！这也是应该实现的严肃改进。进行位向量转换需要消耗大约50%的时间（相对于使用已提供的位向量），但为了严重加速相对于“逻辑”版本，这是值得付出的代价。 - Iterator

Answer 2

这里是逻辑方法、table和bigtabulate的结果，当N等于3亿时：

         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提供已经是整数的值，即在bigtabulate之外进行类型转换1 * cbind(v1,v2)，那么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字节。

我从更新后的基准测试（N=3e8）中删除了表和bigtabulate变体。请注意，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秒的时间内完成此操作。毕竟，即使进行3个（或4个）不同的逻辑操作，它们也可以独立完成，尽管这仍然浪费了一些计算周期。

最相似的挑战者中，logical3B2比table快80倍左右。它比朴素的逻辑操作快10倍左右。而且它仍有很大的提升空间。

以下是生成上述内容的代码。请注意：我建议注释掉一些操作或向量，除非您有大量RAM——创建x、x1和xb以及相应的y对象将占用相当多的内存。

另外，请注意：对于bigtabulate，我应该使用1L作为整数乘数，而不仅仅是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)}
)

- Iterator

---

非常有趣。顺便说一下，tabulate()比table()快得多（实际上是由其调用）。 tabulate(1 + 1L*x + 2L*y)与您的func_logical()相比具有竞争力，但仍然比后者慢5-10％。 - Josh O'Brien

将整数映射是我还没有尝试的另一个技巧。在某个点之后等待太长变得痛苦，而且我对table的谱系感到紧张。 :) 整个制表过程在现代处理器上确实应该少于一秒钟。 此外，正如DWin指出的那样，func_logical()具有说明性 - 第4个求和可以被省略，我们只需从向量长度中减去其他3个值的总和即可。 - Iterator

@JoshO'Brien 我刚刚运行了你的建议。请看我的第二个答案，主要是关于追踪索引，而不是直接处理向量本身。速度很快，但现在最好的方法已经比之前快了10倍以上。 :) 我仍然认为这里有一个10倍的改进空间，但1秒几乎可以交互了。 - Iterator

Answer 3

这是一个使用Rcpp sugar的答案。

N <- 1e8
x <- sample(c(T,F),N,replace=T)
y <- sample(c(T,F),N,replace=T)

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


library(Rcpp)
library(inline)

doCrossTab1 <- cxxfunction(signature(x="integer", y = "integer"), body='
  Rcpp::LogicalVector Vx(x);
  Rcpp::LogicalVector Vy(y);
  Rcpp::IntegerVector V(4);

  V[0] = sum(Vx*Vy);
  V[1] = sum(Vx*!Vy);
  V[2] = sum(!Vx*Vy);
  V[3] = sum(!Vx*!Vy);
  return( wrap(V));
  '
, plugin="Rcpp")

system.time(doCrossTab1(x,y))

require(bit)
system.time(
{
xb <- as.bit(x)
yb <- as.bit(y)
func_logical(xb,yb)
})

这会导致：

> system.time(doCrossTab1(x,y))
   user  system elapsed 
  1.067   0.002   1.069 
> system.time(
+ {
+ xb <- as.bit(x)
+ yb <- as.bit(y)
+ func_logical(xb,yb)
+ })
   user  system elapsed 
  1.451   0.001   1.453 

因此，我们可以在比特包上稍微加速一下，尽管我对竞争时间的激烈程度感到惊讶。

更新：为了纪念Iterator，这里提供了一个Rcpp迭代器解决方案：

doCrossTab2 <- cxxfunction(signature(x="integer", y = "integer"), body='
  Rcpp::LogicalVector Vx(x);
  Rcpp::LogicalVector Vy(y);
  Rcpp::IntegerVector V(4);
    V[0]=V[1]=V[2]=V[3]=0;
  LogicalVector::iterator itx = Vx.begin();
  LogicalVector::iterator ity = Vy.begin();
  while(itx!=Vx.end()){
    V[0] += (*itx)*(*ity);
    V[1] += (*itx)*(!*ity);
    V[2] += (!*itx)*(*ity);
    V[3] += (!*itx)*(!*ity);    
    itx++;
    ity++;
  }
  return( wrap(V));
  '
, plugin="Rcpp")

system.time(doCrossTab2(x,y))
#   user  system elapsed 
#  0.780   0.001   0.782 

- Ian Fellows

---

这些数字非常有趣。Rcpp sugar看起来确实很有帮助；我之前曾尝试过使用逻辑运算来改进（我认为*操作会产生整数类型转换和整数乘法的成本），但今天没有太多时间去调整。我很好奇Dirk或Romain是否会提供一些进一步改进的见解。 - Iterator

你的迭代器解决方案非常吸引人和有教育性。此时我不禁想知道还有什么能够胜过一个迭代器呢？;-) 我期待着去尝试它。根据编译器的智能程度，可以通过避免一些计算（即利用样本非常稀疏的事实，因此在遇到 FALSE-FALSE 时只更新1个条目），并删除 V[3] 来改进平均时间。 因为可以根据其他值和长度确定其值。 - Iterator

Answer 4

另一种策略是仅考虑集合交集，使用TRUE值的索引，利用样本非常偏向（即大多数为FALSE）的优势。

为此，我介绍了func_find01和使用bit包的翻译（func_find01B）；所有不出现在上面答案中的代码都粘贴在下面。

我重新运行了完整的N=3e8评估，但忘记使用func_find01B；我对其进行了更快的方法重新运行，在第二次通过中。

          test replications elapsed   relative user.self sys.self
6   logical3B1            1   1.298   1.000000      1.13     0.17
4    logicalB1            1   1.805   1.390601      1.57     0.23
7   logical3B2            1   2.317   1.785054      2.12     0.20
5    logicalB2            1   2.820   2.172573      2.53     0.29
2       find01            1   6.125   4.718798      4.24     1.88
9 bigtabulate2            1  22.823  17.583205     21.00     1.81
3      logical            1  23.800  18.335901     15.51     8.28
8  bigtabulate            1  27.674  21.320493     24.27     3.40
1        table            1 183.467 141.345917    149.01    34.41

只需要使用“快速”方法：

        test replications elapsed relative user.self sys.self
3     find02            1   1.078 1.000000      1.03     0.04
6 logical3B1            1   1.312 1.217069      1.18     0.13
4  logicalB1            1   1.797 1.666976      1.58     0.22
2    find01B            1   2.104 1.951763      2.03     0.08
7 logical3B2            1   2.319 2.151206      2.13     0.19
5  logicalB2            1   2.817 2.613173      2.50     0.31
1     find01            1   6.143 5.698516      4.21     1.93

所以，在不使用预先转换的位向量的方法中，find01B 是最快的，仅仅比其它方法快了一点点（2.099秒对比 2.327秒）。那么 find02 是从哪里来的呢？之后我写了一个使用预先计算的位向量的版本，这现在是最快的。

通常，“指数法”方法的运行时间可能会受到边际和联合概率的影响。我怀疑当概率更低时，它将特别具有竞争力，但人们必须事先知道这一点，或通过子样本知道。

更新1。我还测试了Josh O'Brien的建议，使用tabulate()而不是table()。结果，经过12秒的时间，大约比find01快2倍，比bigtabulate2慢一半。现在最好的方法接近1秒，这也相对较慢。

 user  system elapsed 
7.670   5.140  12.815 

---

代码：

func_find01 <- function(v1, v2){
    ix1 <- which(v1 == TRUE)
    ix2 <- which(v2 == TRUE)

    len_ixJ <- sum(ix1 %in% ix2)
    len1    <- length(ix1)
    len2    <- length(ix2)
    return(c(len_ixJ, len1 - len_ixJ, len2 - len_ixJ,
             length(v1) - len1 - len2 + len_ixJ))
}

func_find01B <- function(v1, v2){
    v1b = as.bit(v1)
    v2b = as.bit(v2)

    len_ixJ <- sum(v1b & v2b)
    len1 <- sum(v1b)
    len2 <- sum(v2b)

    return(c(len_ixJ, len1 - len_ixJ, len2 - len_ixJ,
             length(v1) - len1 - len2 + len_ixJ))
}

func_find02 <- function(v1b, v2b){
    len_ixJ <- sum(v1b & v2b)
    len1 <- sum(v1b)
    len2 <- sum(v2b)

    return(c(len_ixJ, len1 - len_ixJ, len2 - len_ixJ,
             length(v1b) - len1 - len2 + len_ixJ))
}

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

func_tabulate01 <- function(v1,v2){
    return(tabulate(1L + 1L*x + 2L*y))
}

benchmark(replications = 1, order = "elapsed", 
    table = {res <- func_table(x,y)},
    find01  = {res <- func_find01(x,y)},
    find01B  = {res <- func_find01B(x,y)},
    find02  = {res <- func_find01B(xb,yb)},
    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)},

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

- Iterator

Answer 5

这里有一个使用Rcpp的答案，仅将不是0的条目进行制表。我猜应该有几种方法可以改进这个程序，因为它通常很慢；这也是我第一次尝试Rcpp，所以在移动数据时可能存在一些明显的低效率。我写了一个故意非常简单的示例，让其他人演示如何改进这个程序。

library(Rcpp)
library(inline)
doCrossTab <- cxxfunction(signature(x="integer", y = "integer"), body='
  Rcpp::IntegerVector Vx(x);
  Rcpp::IntegerVector Vy(y);
  Rcpp::IntegerVector V(3);
  for(int i = 0; i < Vx.length(); i++) {
    if( (Vx(i) == 1) & ( Vy(i) == 1) ){ V[0]++; } 
    else if( (Vx(i) == 1) & ( Vy(i) == 0) ){ V[1]++; } 
    else if( (Vx(i) == 0) & ( Vy(i) == 1) ){ V[2]++; } 
 }
  return( wrap(V));
  ', plugin="Rcpp")

N = 3E8 的计时结果：

   user  system elapsed 
 10.930   1.620  12.586 

这比我第二个答案中的func_find01B慢了6倍以上。

- Iterator