有效地计算 R 中大向量中成对差异的直方图？

Posted 2023-02-22

【中文标题】有效地计算 R 中大向量中成对差异的直方图？

【英文标题】：Efficiently compute histogram of pairwise differences in a large vector in R?

【发布时间】：2012-03-20 13:35:32

【问题描述】：

我正在使用 R 中的一个大整数向量（大约 1000 万个整数），我需要从这个向量中找到相差 500 或更少的每一对不同的整数，并制作它们差异的直方图（即每对，第二个减去第一个）。

这是完全非矢量化的代码，用于非常缓慢地完成我想做的事情：

# Generate some random example data
V <- round(rnorm(100) * 1000)

# Prepare the histogram
my.hist <- rep(0, 500)
names(my.hist) <- as.character(seq(1,500))
for (x1 in V) 
    for (x2 in V) 
        difference = x2 - x1
        if (difference > 0 && difference <= 500) 
            my.hist[difference] = my.hist[difference] + 1
        
    

（假设每个整数都是唯一的，所以difference > 0 位是可以的。这是允许的，因为我实际上并不关心差异为零的任何情况。）

下面是一些矢量化内循环的代码：

my.hist2 <- rep(0, 500)
names(my.hist2) <- as.character(seq(1,500))
for (x1 in V) 
    differences <- V[V > x1 & V <= x1+500] - x1
    difftable <- table(differences)
    my.hist2[names(difftable)] = my.hist2[names(difftable)] + difftable

这肯定比第一个版本快。然而，当V 仅包含 500000 个元素（半百万）时，即使这个变体也已经太慢了。

我可以在没有任何显式循环的情况下执行此操作，如下所示：

X <- combn(V, 2)
# X is a matrix with two rows where each column represents a pair
diffs <- abs(X[2,] - X[1,])
my.hist3 <- table(diffs[diffs <= 500])

但是矩阵 X 将包含 10e6 * (10e6 - 1) / 2，或大约 50,000,000,000,000 列，这可能是个问题。

那么有没有办法在不显式循环（太慢）或扩展所有对的矩阵（太大）的情况下做到这一点？

如果你想知道为什么我需要这样做，我正在实现这个：http://biowhat.ucsd.edu/homer/chipseq/qc.html#Sequencing_Fragment_Length_Estimation

【问题讨论】：

我是否正确理解您有 50,000,000,000,000 列，每列包含 10,000,000 个元素？我很想看到这个问题得到解决 :) 规模使它成为一个了不起的问题。

@user306848：我认为 Ryan 的意思是 combn 将返回一个 2 x 50E12 矩阵对（所有可能的组合）。

是的，@jbaums 是正确的。

感谢您的澄清 :)

你的向量的范围/分布是什么样的？有多少个唯一整数？

【参考方案1】：

一种可能的改进是对数据进行排序： 距离小于 500 的对 (i,j) 然后将接近对角线， 而且您不必探索所有的价值。

代码看起来像这样（它仍然很慢）。

n <- 1e5
k <- 500
V <- round(rnorm(n) * n * 10)
V <- as.integer(V)
V <- sort(V)
h <- rep(0,k)

for(i in 1:(n-1)) 
  for(j in (i+1):n) 
    d <- V[j] - V[i]
    if( d > k ) break
    if( d > 0 ) h[d] <- h[d]+1
  

编辑：如果你想要更快的东西，你可以用 C 编写循环。 1000 万个元素需要 1 秒。

n <- 10e6
k <- 500
V <- round(rnorm(n) * n * 10)
V <- as.integer(V)
V <- sort(V)
h <- rep(0,k)

library(inline)
sig <- signature(n="integer", v="integer", k="integer", h="integer")
code <- "
  for( int i = 0; i < (*n) - 1; i++ ) 
    for( int j = i + 1; j < *n; j++ ) 
      int d = v[j] - v[i];
      if( d > *k ) break;
      if( d > 0 ) h[d-1]++;
    
  
"
f <- cfunction( sig, code, convention=".C" )
h <- f(n,V,k,h)$h

【讨论】：

希望我能投票更多。我想到了相同的策略，但您的实施要好得多。

C 内联功能是否不支持函数参数和返回语句？

实际上，我对事情是如何从 R 传递到 C 有点困惑。在 C 代码中，变量 n、v、k 和 h 是指向 C 整数和数组的指针整数，还是 R 向量？还是 R 整数向量的内存表示与 C 整数数组兼容，所以实际上两者兼而有之？

在 R 和 C 之间交换数据有两种约定：.Call，它使用 R 对象（C 中的 SEXPR 类型），和.C，它使用指针（指向 C 数组）。如果您只有数字数组，.C 约定更易于使用，但您还必须传递数组的大小（这通常在包装函数中完成，因此最终用户永远不会看到它）。在大多数情况下，在 R 中创建包含结果的对象并在 C 中填充它们会更容易。Writing R extensions 手册中有更多信息。

好吧，看了那本手册，我想我明白了。当使用.C 约定传递给cfunction 时，R 向量将转换为C 数组。返回值呢？从我创建的一个玩具示例来看，它看起来像是返回签名中的变量列表，并转换回 R 向量。这是正确的吗？