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

Question

我正在处理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

Answer 1

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

代码看起来像这样(它仍然很慢)。

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万个元素需要花费1s。

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