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

ssm49v7z  于 2023-11-14  发布在  其他
关注(0)|答案(5)|浏览(99)

对于两个长度> 1 E8的逻辑向量x和y，计算2x2交叉表的最快方法是什么？
我怀疑答案是用C/C++编写它，但我想知道R中是否有一些东西已经非常聪明地解决了这个问题，因为它并不罕见。
示例代码，用于300 M条目（如果3E 8太大，可以让N = 1 E8;我选择的总大小略低于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
    1.简单的逻辑运算（例如sum(x & y)）
    1.向量乘法（BOO）
  3. data.table
    1.以上部分，使用multicore包中的parallel（或新的parallel包）
    我已经尝试了前三个选项（见我的答案），但我觉得必须有更好更快的东西。
    我发现table的工作速度非常慢。bigtabulate对于一对逻辑向量来说似乎有点大材小用。最后，执行普通逻辑运算似乎是一种组装，它会多次查看每个向量（3X？7 X？），更不用说它在处理过程中会占用大量额外的内存，这是一个巨大的时间浪费。
    向量乘法通常是一个坏主意，但当向量稀疏时，可以将其存储为这样，然后使用向量乘法。
    既然一个次级样本可能足以制作一个交叉表格，为什么要使用这么多数据？这些数据来自两个变量的TRUE观测值非常罕见的情况。一个是数据异常的结果，另一个是由于代码中可能存在错误（可能的bug，因为我们只看到计算结果-将变量x视为“垃圾输入”，而y为“垃圾输出”。因此，问题是，由代码引起的输出中的问题是否仅仅是数据异常的情况，还是存在其他一些好数据变坏的情况？（这就是为什么我问了一个关于stopping when a NaN , NA , or Inf is encountered的问题。）
    这也解释了为什么我的例子中TRUE值的概率很低;这些值实际上发生的概率远远低于0.1%。
    这是否意味着一种不同的解决途径？是的：它表明我们可以使用两种指数（即TRUE在每个集合中的位置）和计数集合交集。我避免了集合交集，因为我被Matlab烧了一段时间，它会先对集合中的元素进行排序，然后再进行交集。（我隐约记得复杂性更令人尴尬：比如O(n^2)而不是O(n log n)。
    那么，在R中如何实现呢？
r

5条答案

按热度按时间
hjqgdpho

hjqgdpho1#

如果你要对巨大的逻辑向量进行大量的操作，看看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

字符串

赞(0）分享  回复(0）举报  2023-11-14
px9o7tmv

px9o7tmv2#

这个答案给出了三个简单方法的时间，这是相信table很慢的基础。然而，要认识到的关键是“逻辑”方法效率非常低。看看它在做什么：

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

不仅如此，它甚至没有被编译或并行化。然而，它仍然击败了table。注意，bigtabulate，加上 * 一个额外的类型转换 *（1 * cbind...）仍然击败了table。
下面是逻辑方法table和bigtabulate的结果，N = 3E 8：

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 = 3E 6：

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= 3E 6的倍数是1.80，而不是2.4。相对于“逻辑”方法的N= 3E 8的倍数只有1.21，而不是1.53。
更新2

正如约书亚乌尔里希所指出的，转换为位向量是一个重大的改进-我们分配和移动的数据少了很多：R的逻辑向量每个条目占用4个字节（“为什么？“，你可能会问... Well, I don't know, but an answer may turn up here.），而一个位向量消耗，嗯，一个位，每个条目-即1/32的数据。所以，x消耗1.2e9字节，而xb（下面代码中的位版本）仅消耗3.75e7字节。
我已经从更新的基准测试中删除了table和bigtabulate变体（N= 3e 8）。请注意，logicalB1假设数据已经是位向量，而logicalB2是相同的操作，但类型转换会受到惩罚。由于我的逻辑向量是对其他数据进行操作的结果，我没有从位向量开始的好处。尽管如此，付出的代价相对较小。[“logical 3”系列只执行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，比table快80倍左右，比朴素逻辑运算快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)}
)
赞(0）分享  回复(0）举报  2023-11-14
mi7gmzs6

mi7gmzs63#

这是一个使用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

所以，我们可以得到一个小的速度超过位包，虽然我很惊讶，如何竞争激烈的时代。
更新：为了荣誉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
赞(0）分享  回复(0）举报  2023-11-14
aurhwmvo

aurhwmvo4#

一种不同的策略是只考虑集合的交集，使用TRUE值的索引，利用样本非常有偏见（即主要是FALSE）。
为此，我介绍了func_find01和一个使用bit包的翻译（func_find01B）;所有没有出现在上面答案中的代码都粘贴在下面。
我重新运行了完整的N= 3e 8评估，除了忘记使用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奥布莱恩的建议进行了计时，使用tabulate()而不是table()。结果，在12秒过去后，大约是2X 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)}
)
赞(0）分享  回复(0）举报  2023-11-14
9w11ddsr

9w11ddsr5#

这里有一个关于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倍多的时间。

赞(0）分享  回复(0）举报  2023-11-14
我来回答

相关问题

    查看更多

    热门标签

    更多
    JavaquerypythonNode开发语言requestUtil数据库Table后端算法LoggerMessageElementParser

    最新问答

    更多
    • 蜀ICP备13028337号-1 大数据知识库 https://www.saoniuhuo.com © All rights reserved
    • 本站内容来源互联网,如果侵犯您的权益请联系我们删除, 联系方式：448109455@qq.com