生成多个邻接矩阵的更快方法答案

【问题标题】：Faster way to generate multiple adjacency matrix生成多个邻接矩阵的更快方法
【发布时间】：2020-10-27 03:01:19
【问题描述】：

假设我有一个任意概率矩阵P，如下所示，

P = matrix(c(0.3,0.2,0.2,0.2,0.3,0.2,0.2,0.2,0.3),3,3)
P 
      [,1] [,2] [,3]
[1,]  0.3  0.2  0.2
[2,]  0.2  0.3  0.2
[3,]  0.2  0.2  0.3

对于单个邻接矩阵，它的生成类似于（未加权，无自放样）

tem = matrix(runif(3^2), nrow = 3)
tmpG = 1 * (tmpmat < P)
tmpG[lower.tri(tmpG)] <- 0
tmpG <- t(tmpG) + tmpG - diag(diag(tmpG))

但是，如果我需要生成100个邻接矩阵怎么办，所以我写下以下代码

G = list()
for (i in 1:rep) {
  tmpmat = matrix(runif(n^2), nrow = n)
  tmpG = 1 * (tmpmat < P)
  tmpG[lower.tri(tmpG)] <- 0
  tmpG <- t(tmpG) + tmpG - diag(diag(tmpG))
  if (noloop) {
    diag(tmpG) = 0
  }
  G[[i]] = tmpG
}

在我的情况下，n >10000 和 T = 1000，所以它非常慢，有什么更好的改进方法吗？

【问题讨论】：

编写的代码（加上 noloop=FALSE、rep replicate 将成为您避免Second Circle of the R Inferno 的朋友
在我的例子中，n > 10000, T = 1000。谢谢
你的例子中的小错字，你有tem，我认为你的意思是tmpmat？另外，您能否澄清一下，当您说n > 10000 时，您当前使用的维度是3 吗？当您说T = 1000 时，是您当前使用的复制次数（未定义）rep？
Dubukay 是对的，replicate(rep, {<<code inside your for loop>>}) 会更快，但是当你的矩阵变大时，它可能已经无关紧要了。如果不改进算法，我认为您不会获得太多速度。也许你可以只操作上面的三角形？

标签： r igraph rcpp rcpparmadillo

【解决方案1】：

我认为我们可以做得更好，只使用所需长度的向量，并在最后将其放入矩阵中。我没有仔细检查过，你的代码没有可供我比较意图的 cmets，所以在信任它之前请确保它是正确的。

p_vec = P[upper.tri(P, diag = !noloop)]
nn = length(p_vec)

tmpG_vec = runif(nn) < p_vec
tmpG = matrix(0, n, n)
tmpG[upper.tri(tmpG, diag = !noloop)] = tmpG_vec
tmpG[lower.tri(tmpG, diag = !noloop)] = tmpG_vec
tmpG

然后我们可以将其包装在 replicate 中以进行迭代。

在更多维度/更高次数上进行基准测试，我们得到了大约 25% 的加速，但仍然很慢（我放弃了 n = 5000 的基准测试，因为我厌倦了等待）。通过并行运行，您可能会获得相当大的速度 - 如果您有 8 个内核，则可以说几乎是 8 倍的加速。参见，例如，this question，尽管可能有更现代的方法来做到这一点。

rep = 5L
n = 2000
noloop = TRUE

P = matrix(runif(n^2), n)
P = P %*% t(P)
P = P / colSums(P)

p_vec = P[upper.tri(P, diag = !noloop)]
nn = length(p_vec)


microbenchmark::microbenchmark(
  loop = {
    G = list()
    for (i in 1:rep) {
      tmpmat = matrix(runif(n^2), nrow = n)
      tmpG = 1 * (tmpmat < P)
      tmpG[lower.tri(tmpG)] <- 0
      tmpG <- t(tmpG) + tmpG - diag(diag(tmpG))
      if (noloop) {
        diag(tmpG) = 0
      }
      G[[i]] = tmpG
    }
  },
  diagonal = replicate(rep, {
    tmpG_vec = runif(nn) < p_vec
    tmpG = matrix(0, n, n)
    tmpG[upper.tri(tmpG, diag = !noloop)] = tmpG_vec
    tmpG[lower.tri(tmpG, diag = !noloop)] = tmpG_vec
    tmpG
  }),
  times = 5L
)

# Unit: seconds
#      expr      min       lq     mean   median       uq      max neval
#      loop 1.525028 1.614544 2.136637 2.148771 2.387423 3.007417     5
#  diagonal 1.312022 1.360457 1.592914 1.444902 1.602536 2.244652     5

【讨论】：