R：是否可以矢量化/加速这个双循环？答案

【问题标题】：R: Is it possible to vectorise / speed-up this double loop?R：是否可以矢量化/加速这个双循环？
【发布时间】：2012-05-23 21:29:26
【问题描述】：

这是一个高层次的一般性问题。周围有一些类似的例子，但有不同的、更简洁的例子。或许无法回答。 conn 是一个矩阵。

     for (i in 2:dim(conn)[1]) {
        for (j in 2:dim(conn)[1]) {
          if ((conn[i, 1] == conn[1, j]) & conn[i, 1] != 0) {
              conn[i, j] <- 1
              conn[j, i] <- 1
              }
              else {
                conn[i, j] <- 0
                conn[j, i] <- 0
                }
           }
      }

这直接来自 clusterCons 包中的cluscomp。

我的问题很简单：是否可以加快循环或将其矢量化？作为 R 初学者，我看不到它，也不想以挫败感告终，因为这可能是不可能的。我会接受任何可以说是或否的答案，并暗示可能涉及的工作量。

【问题讨论】：

请发布示例数据和预期结果。另外，试着用文字描述代码应该做什么。
是的，在这种情况下很可能可以找到一个非常快速的解决方案。我看看能不能快点煮熟。
矩阵conn是否总是对称的？
您的意思是“cluscomp”吗？ clusterCons 中没有 conscomp。
查看 cluscomp 代码，似乎 (a) 这指的是 cluscomp，是的，并且 (b) 数据始终是对称的。嗯，我有个主意……

标签： r loops vectorization

【解决方案1】：

非矩阵解决方案 - 应该非常快，假设 conn 是非负的且对称的......

connmake = function(conn){
  ordering = order(conn[,1])
  breakpoints = which(diff(conn[ordering,1]) != 0)
  if (conn[ordering[1], 1] != 0){
    breakpoints = c(1, breakpoints + 1, nrow(conn) + 1)
  } else {
    breakpoints = c(breakpoints + 1, nrow(conn) +1)
  }
  output = matrix(0, nrow(conn), nrow(conn))

  for (i in 1:(length(breakpoints) - 1)){
    output[ ordering[breakpoints[i]:(breakpoints[i+1] -1)],
        ordering[breakpoints[i]:(breakpoints[i+1] -1)]] =  1
  }
  output[,1] = conn[,1]
  output[1,] = conn[,1]
  output
}

一些使用早期基准测试的测试代码。（原代码实现为orig()，f2()是较早的建议。）

size = 2000
conn  = matrix(0, size, size)
conn[1,] = sample( 1:20, size, replace = T)
conn[,1] = conn[1,]

system.time(orig(conn) -> out1)
#user  system elapsed 
#20.54    0.00   20.54 
system.time(f2(conn) -> out2)
#user  system elapsed
#0.39    0.02    0.41 
system.time(connmake(conn) -> out3)
#user  system elapsed 
#0.02    0.00    0.01 
identical(out1, out2)
#[1] TRUE
identical(out1, out3)
#[1] TRUE

请注意，对于包含 0 的 conn，f2 实际上会失败，但不是我的问题，是吗？带有负值的 conn 可以简单地处理，例如通过安全偏移量增加相关值。非对称连接需要更多思考，但应该是可行的......

一般的教训是排序比成对比较快。成对比较是 O(N^2)，而 R 中最慢的排序算法是 O(N^4/3)。数据排序后，比较就变得微不足道了。

【讨论】：

另外，如果它仍然不够快，可能有一种方法可以摆脱那里的循环，但我怀疑它是否值得。
+1 表示不错的算法。感谢您指出我忘记了0 的情况，我已经修复了我的代码。
我向您的 R 专业知识低头。事实上 f2 看起来很有前途，速度很快，但没有按预期工作，因此，还没有“答案”。明天早上将与 connmake 一起去。再次感谢！
我不确定哪些对象进出，但我必须添加rownames(output) <- row.names(x)和colnames(output) <- row.names(x)才能使其工作。

【解决方案2】：

我会这样写，用outer 代替双循环。请注意，它仍然进行了比需要更多的计算，但肯定更快。我假设conn 是一个方阵。

原代码：

f1 <- function(conn) {
   for (i in 2:dim(conn)[1]) {
      for (j in 2:dim(conn)[1]) {
         if ((conn[i, 1] == conn[1, j]) & conn[i, 1] != 0) {
            conn[i, j] <- 1
            conn[j, i] <- 1
         } else {
            conn[i, j] <- 0
            conn[j, i] <- 0
         }
      }
   }
   return(conn)
}

我的建议：

f2 <- function(conn) {
   matches <- 1*outer(conn[-1,1], conn[1,-1], `==`)
   matches[conn[-1,1] == 0, ] <- 0
   ind <- upper.tri(matches)
   matches[ind] <- t(matches)[ind]
   conn[-1,-1] <- matches
   return(conn)
}

一些示例数据：

set.seed(12345678)
conn <- matrix(sample(1:2, 5*5, replace=TRUE), 5, 5)
conn
#      [,1] [,2] [,3] [,4] [,5]
# [1,]    2    2    1    2    1
# [2,]    1    1    2    2    1
# [3,]    2    2    1    2    1
# [4,]    2    2    2    2    1
# [5,]    1    1    2    2    1

结果：

f1(conn)
#      [,1] [,2] [,3] [,4] [,5]
# [1,]    2    2    1    2    1
# [2,]    1    0    1    1    0
# [3,]    2    1    0    0    1
# [4,]    2    1    0    1    0
# [5,]    1    0    1    0    1

identical(f1(conn), f2(conn))
# [1] TRUE

一个更大的例子，时间比较：

set.seed(12345678)
conn <- matrix(sample(1:2, 1000*1000, replace=TRUE), 1000, 1000)

system.time(a1 <- f1(conn))
# user  system elapsed 
# 59.840   0.000  57.094 

system.time(a2 <- f2(conn))
# user  system elapsed 
# 0.844   0.000   0.950 

identical(a1, a2)
# [1] TRUE

也许不是你能得到的最快的方法（我相信这里的其他人可以使用编译器或 Rcpp 更快地找到），但我希望对你来说足够短和快。

编辑：由于已经指出（从提取此代码的上下文中）conn 是一个对称矩阵，我的解决方案可以缩短一点：

f2 <- function(conn) {
   matches <- outer(conn[-1,1], conn[1,-1],
                    function(i,j)ifelse(i==0, FALSE, i==j)) 
   conn[-1,-1] <- as.numeric(matches)
   return(conn)
}

【讨论】：

最后一个（哦，太优雅了）版本让我的机器崩溃了，但另一个运行正常。

【解决方案3】：

我想到了几件事。

首先，您可以通过仅遍历对角线下方或对角线上方的条目将时间缩短一半。如果矩阵是方形的，则任何一个都可以。如果是dim(conn)[1] > dim(conn)[2]，那么你会想要使用类似

的东西来遍历左下角的三角形

for (j in 2:dim(conn)[2]) {
  for (i in j:dim(conn)[1]) {
    ...
  }
}

其次，人们可能会尝试使用apply，它是同类，因为它们通常会显着减少时间。但是，在这种情况下，每个 [i,j] 单元格同时引用列头 [1,j] 和行头 [i,1]，这意味着我们不能只将单元格、行或列发送到 *pply。为了代码清晰，我可能会保留for 循环。任何有效的基于 *pply 的技巧都会非常聪明，以至于一年后我会忘记它是如何工作的。

最后，这似乎是一个经典示例，说明使用从 R 调用的 C 会快得多。这似乎需要做很多工作，但它比您想象的要容易得多，甚至（对于这个特定的例如）如果你不知道 C。第一个从 R 调用 C 对我来说有意义的简短示例是 here，但它没有利用 Rcpp，所以我不会停在那里。或者，如果您从 any 工作 Rcpp 代码的简单示例开始，那么您可以在此处对其进行修改以执行您想要的操作。如果你只是想修改别人的代码，请以this StackOverflow thread开头。

【讨论】：