15

对于长度 > 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. 以上一些，带有parallel来自multicore包（或新parallel包）

我已经尝试了前三个选项（请参阅我的答案），但我觉得必须有更好更快的东西。

我发现它的table工作非常缓慢。 bigtabulate对于一对逻辑向量来说似乎有点矫枉过正。最后，执行普通的逻辑运算似乎是一个杂项，它查看每个向量的次数太多（3X？7X？），更不用说它在处理过程中会占用大量额外的内存，这是一个巨大的时间浪费。

向量乘法通常是一个坏主意，但是当向量稀疏时，将其存储起来，然后使用向量乘法可能会获得优势。

随意改变N，p如果这将展示制表函数的任何有趣行为。:)

更新1。我的第一个答案给出了三种幼稚方法的时间，这是相信table缓慢的基础。然而，要认识到的关键是“逻辑”方法效率极低。看看它在做什么：

  • 4个逻辑向量运算
  • 4 种类型转换（逻辑到整数或 FP - for sum）
  • 4个向量求和
  • 8个赋值（1个用于逻辑运算，1个用于求和）

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

更新 2. 以免有人指出 R 中的逻辑向量支持NA，并且这将成为这些交叉表系统中的扳手（在大多数情况下都是如此），我应该指出我的向量来自is.na()or is.finite()。:) 我一直在调试NA和其他非有限值——它们最近让我很头疼。如果您不知道您的所有条目是否都是，您可以在采用本问答中出现的一些想法之前进行NA测试- 这将是明智的。any(is.na(yourVector))

更新 3. Brandon Bertelsen 在评论中提出了一个非常合理的问题：为什么在子样本（初始集毕竟是样本 ;-)）可能足以创建交叉样本时使用这么多数据制表？不要在统计数据中走得太远，但数据来自TRUE两个变量的观察非常罕见的情况。一个是数据异常的结果，另一个是由于代码中可能的错误（可能是错误，因为我们只看到计算结果 - 将变量x视为“Garbage In”和y“Garbage Out”。结果，问题是代码导致的输出问题是否仅仅是数据异常的情况，还是有其他一些好的数据变坏的情况？（这就是为什么我问一个关于遇到NaN,NA或时停止Inf。）

这也解释了为什么我的示例的值概率很低TRUE；这些确实发生的几率远低于 0.1%。

这是否暗示了不同的解决方案？是的：这表明我们可以使用两个索引（即TRUE每个集合中的位置）并计算集合交点。我避免设置交叉点，因为我被 Matlab 烧了一段时间（是的，这是 R，但请耐心等待），它会在交叉点之前先对集合的元素进行排序。（我隐约记得复杂性更令人尴尬：likeO(n^2)而不是O(n log n).）

4

5 回答 5

11

如果您对巨大的逻辑向量进行大量操作，请查看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
于 2012-02-07T11:16:51.393 回答
9

以下是N = 3E8的逻辑方法table和的结果：bigtabulate

         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 字节。

我已经放弃了更新基准（N = 3e8）table的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 秒的时间内完成此操作应该是可行的，更改包括以下一项或多项：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)}
)
于 2012-02-07T06:05:40.090 回答
3

这是使用 Rcpp 糖的答案。

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

因此，我们可以在 bit 包上稍微加快速度，尽管我对时代的竞争程度感到惊讶。

更新：为了纪念迭代器，这是一个 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
于 2012-02-14T03:12:45.590 回答
2

一种不同的策略是考虑只设置交点，使用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 秒后的结果是大约 2 倍find01和大约一半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)}
)
于 2012-02-08T15:42:27.087 回答
0

这是一个答案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

这比我的第二个答案要多 6 倍func_find01B。

于 2012-02-13T13:18:53.380 回答