如何将用于排列的for循环向量化？

Question

我正在使用 R 进行数据分析，希望进行一次置换检验。为此，我使用了一个 for 循环，但是速度非常慢，我想尽可能地加快这段代码。 我认为矢量化是解决问题的关键。 然而，经过几天的尝试之后，我仍然没有找到一个合适的解决方案来重新编码它。 我非常感谢您的帮助！

我有一个对称矩阵，其中包含种群之间的生态距离（"dist.mat"）。 我想随机打乱这个距离矩阵的行和列，以生成一个置换后的距离矩阵（"dist.mat.mix"）。 然后，我想保存这个置换后的距离矩阵中的上三角值（大小为 "nr.pairs"）。 这个过程应该重复几次（"nr.runs"）。结果应该是一个矩阵("result")，包含多次重复运行的置换上三角值，其维度为nrow=nr.runs和ncol=nr.pairs。以下是一个使用 for 循环完成所需操作的示例 R 代码：

# example number of populations
nr.pops <- 20

# example distance matrix
dist.mat <- as.matrix(dist(matrix(rnorm(20), nr.pops, 5)))

# example number of runs
nr.runs <- 1000

# find number of unique pairwise distances in distance matrix
nr.pairs <- nr.pops*(nr.pops-1) / 2

# start loop
result <- matrix(NA, nr.runs, nr.pairs)
for (i in 1:nr.runs) {
  mix <- sample(nr.pops, replace=FALSE)
  dist.mat.mix <- dist.mat[mix, mix]
  result[i, ] <- dist.mat.mix[upper.tri(dist.mat.mix, diag=FALSE)]
}

# inspect result
result

我已经使用base::replicate函数做了一些笨拙的向量化尝试，但这并没有加速。实际上它更慢了：

# my for loop approach
my.for.loop <- function() {
  result <- matrix(NA, nr.runs, nr.pairs)
  for (i in 1:nr.runs){
    mix <- sample(nr.pops, replace=FALSE)
    dist.mat.mix <- dist.mat[mix ,mix]
    result[i, ] <- dist.mat.mix[upper.tri(dist.mat.mix, diag=FALSE)]
  }
}

# my replicate approach
my.replicate <- function() {
  results <- t(replicate(nr.runs, {
    mix <- sample(nr.pops, replace=FALSE)
    dist.mat.mix <- dist.mat[mix, mix]
    dist.mat.mix[upper.tri(dist.mat.mix, diag=FALSE)]
  }))
}

# compare speed
require(microbenchmark)
microbenchmark(my.for.loop(), my.replicate(), times=100L)

# Unit: milliseconds
# expr           min     lq      mean    median  uq      max       neval
# my.for.loop()  23.1792 24.4759 27.1274 25.5134 29.0666 61.5616   100
# my.replicate() 25.5293 27.4649 30.3495 30.2533 31.4267 68.6930   100    

如果您知道如何使用简洁的向量化解决方案加速我的for循环，我将不胜感激。这种方法是否可行？

- HaRa

---

你的真实数据有多大？你想运行多少个 nr.runs？ - minem

nr.pops 是 20，而 nr.runs 是 1000。我已经在问题中更新了这个。 - HaRa

这段程序哪里慢？你做了什么导致它变慢了？一次 my.for.loop() 的运行时间为 0.01 秒。你是在多次调用它吗？如果是，调用了多少次？调用该函数1000次需要大约20秒。对我来说这似乎不算慢。 - minem

这个for循环是一个重复数亿次的大型计算的一部分。这个大型计算运行了几周时间。即使只是减少这个for循环的计算时间几微秒，也可以为我节省几个小时甚至几天的运行时间。我发现这个for循环是我整个计算中最大的瓶颈，因此我想尽可能地让它快速运行。 - HaRa

“replicate”并不是向量化，它基本上与sapply/lapply相同。而这两者基本上都是“for”循环。https://dev59.com/QlgQ5IYBdhLWcg3wqF2l - minem

Answer 1

略微更快：

minem <- function() {
  result <- matrix(NA, nr.runs, nr.pairs)
  ut <- upper.tri(matrix(NA, 4, 4)) # create upper triangular index matrix outside loop
  for (i in 1:nr.runs) {
    mix <- sample.int(nr.pops) # slightly faster sampling function
    result[i, ] <- dist.mat[mix, mix][ut]
  }
  result
}
microbenchmark(my.for.loop(), my.replicate(), minem(), times = 100L)
# Unit: microseconds
# expr               min      lq      mean   median       uq      max neval cld
# my.for.loop()   75.062  78.222  96.25288  80.1975 104.6915  249.284   100   a
# my.replicate() 118.519 122.667 152.25681 126.0250 165.1355  495.407   100   a
# minem()         45.432  48.000 104.23702  49.5800  52.9380 4848.986   100   a

更新： 我们可以稍微不同的方式获得必要的矩阵索引，这样我们就可以一次性地对元素进行子集操作。

minem4 <- function() {
  n <- dim(dist.mat)[1]
  ut <- upper.tri(matrix(NA, n, n))
  im <- matrix(1:n, n, n)
  p1 <- im[ut]
  p2 <- t(im)[ut]
  dm <- unlist(dist.mat)

  si <- replicate(nr.runs, sample.int(nr.pops))
  p <- (si[p1, ] - 1L) * n + si[p2, ]
  result2 <- matrix(dm[p], nr.runs, nr.pairs, byrow = T)
  result2
}

microbenchmark(my.for.loop(), minem(), minem4(), times = 100L)
# Unit: milliseconds
# expr                min        lq     mean    median        uq       max neval cld
# my.for.loop() 13.797526 14.977970 19.14794 17.071401 23.161867  29.98952   100   b
# minem()        8.366614  9.080490 11.82558  9.701725 15.748537  24.44325   100  a 
# minem4()       7.716343  8.169477 11.91422  8.723947  9.997626 208.90895   100  a 

更新2： 我们可以使用dqrng样本函数获得额外的加速。

minem5 <- function() {
  n <- dim(dist.mat)[1]
  ut <- upper.tri(matrix(NA, n, n))
  im <- matrix(1:n, n, n)
  p1 <- im[ut]
  p2 <- t(im)[ut]
  dm <- unlist(dist.mat)

  require(dqrng)
  si <- replicate(nr.runs, dqsample.int(nr.pops))
  p <- (si[p1, ] - 1L) * n + si[p2, ]
  result2 <- matrix(dm[p], nr.runs, nr.pairs, byrow = T)
  result2
}

microbenchmark(my.for.loop(), minem(), minem4(), minem5(), times = 100L)
# Unit: milliseconds
# expr                min        lq      mean    median        uq      max neval  cld
# my.for.loop() 13.648983 14.672587 17.713467 15.265771 16.967894 36.18290   100    d
# minem()        8.282466  8.773725 10.679960  9.279602 10.335206 27.03683   100   c 
# minem4()       7.719503  8.208984  9.039870  8.493231  9.097873 25.32463   100  b  
# minem5()       6.134911  6.379850  7.226348  6.733035  7.195849 19.02458   100 a  

- minem

---

非常感谢，这正是我想要的！ - HaRa

请您能否详细说明一下矩阵索引背后的逻辑？p1和p2的作用是什么，您是如何得出p以及最终的矩阵results2的呢？ - HaRa