如何在R中高效地计算大向量中成对差异的直方图？

Question

如何在R中高效地计算大向量中成对差异的直方图？

10

我在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

- Ryan C. Thompson

我理解的对吗，你有 50,000,000,000,000 列，每列有 10,000,000 个元素？我很想看到这个问题的解决 :) 这个规模使它成为一个令人敬畏的问题。 - ma cılay

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

是的，@jbaums 是正确的。 - Ryan C. Thompson

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

假设没有两个整数相等，它们的范围从0到100,000,000。在这个范围内的分布预计会高度不均匀，有许多类似值的聚集体。（这些聚集体代表人类基因组中特定蛋白质附着的坐标。） - Ryan C. Thompson

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Vincent Zoonekynd · Accepted Answer

一种可能的改进是对数据进行排序：距离低于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