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

问题描述

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

P = matrix(c(0.3,0.2,0.3,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 >10000T = 1000太慢了,有什么更好的主意来改善它吗?

解决方法

我认为我们可以做得更好,方法是只处理所需长度的向量,然后将其放入矩阵的最后。我尚未对此进行仔细检查,并且您的代码中没有供我比较意图的注释,因此请在信任它之前确保它是正确的。

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)
    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