在 R 中使用嵌套循环优化（矢量化？） For 循环答案

【问题标题】：Optimizing (Vectorizing?) For loop with nested loop in R在 R 中使用嵌套循环优化（矢量化？） For 循环
【发布时间】：2014-03-11 16:15:24
【问题描述】：

我正在迭代地使用 rdist 来计算大型数据集的最近邻。目前我有一个相当小的矩阵，包含 634,000 个向量，有 6 列。

如前所述，我使用 rdist 来计算每个向量到每个其他向量的距离，每个距离计算都是一个步骤。此外，我在每一步都运行一个函数，计算 k=1,2,3,4 最近的邻居并取和（实际上是 k=所有邻居）。

###My function to compute k nearest neighbours from distance vector

    knn <- function (vec,k) {
      sum((sort(vec)[1:k+1]))
    }

###My function to compute nearest neighbours iteratively for every vector
myfunc <- function (tab) {

  rowsums <- numeric(nrow(tab)) ###Here I will save total sums
  knnsums_log <- matrix(nrow=nrow(tab),ncol=4) ###Matrix for storing each of my kNN sums

  for(i in 1:nrow(tab)) { ###For loop to compute distance and total sums
    q<-as.matrix(rdist(tab[i,],tab))
    rowsums[i] <- rowSums(q)

     for (k in c(1:4)) { ###Nested loop to run my knn function
     knnsums[i,k] <- knn(q,k) 
    }

  }

  return(cbind(rowsums,knnsums_log))
}

数据样例（634k 行）

    X1  X2  X3  X4  X5  X6
1   0.00    0.02    0   0   0.02    -0.263309267
2   0.00    0.02    0   0   0.02    -0.171764667
3   0.00    0.02    0   0   0.02    -0.128784869
4   0.00    0.02    0   0   0.02    -0.905651733

对于那些不熟悉函数 rdist 的人，可以得到争论之间的欧几里得距离。它的工作速度比自定义编写的函数快得多。它比 dist 更适用，因为 dist 仅在矩阵距离内计算。从技术上讲，我知道这就是我正在做的事情，但 dist 试图将其存储在内存中，而且它太大了，甚至无法考虑这样做。

我怎样才能使上述工作更好？我试过弄乱应用函数，但没有得到任何有用的东西。我希望我已经清楚地解释了一切。如果我的数学是正确的，最坏的情况估计需要我一个多星期才能运行该代码。我有非常强大的服务器来处理它。但是没有 GPU。我还没有尝试过多核（应该有 12 个可用），但是我不知道如何为每个核委派。

感谢您的帮助。

【问题讨论】：

试试parallels 包和函数mcapply，它将为您完成分发。回到根本问题：如果有任何方法可以让您“智能猜测”您需要哪些最近邻集合，而不是计算所有这些集合，您或许可以减少所需的计算。

标签： r optimization for-loop distance neighbours

【解决方案1】：

几个提示：

0) 使用 Rprof 和 line.profiling 选项分析您的代码

1) R 中的矩阵是按列排列的。因为您比较它们之间的向量，所以如果将它们存储为矩阵的列会更快

2) 我不知道 rdist 函数来自哪里，但是你应该避免使用 as.matrix(rdist(tab[i,],tab)) 来复制和创建一个新矩阵

3) 你可以优化你的 knn() 函数，它对同一个向量进行 4 次排序

4) 为什么不只是 rdist(tab) ？

【讨论】：

感谢您的提示。一些额外的事情：Knn 比 rdist 减慢了函数的速度（顺便说一下，它来自 fields 包）。 as.matrix 实际上加速了代码。不知道为什么？我已经测试了有无（最初我虽然它会减慢速度）但执行实际上受到了重大打击。再次感谢您的提示！
rdist(tab[i,],tab) 的类型是什么？
这是一个矩阵，同样令人困惑，但基准测试显示出明显的差异。也许这是一个奇怪的代码打嗝。
您能否再次检查 as.matrix 是否加快了您的代码速度？因为我看不到任何合理的解释。你是如何对你的代码进行基准测试的？！
嗨，卡尔，很抱歉回复晚了，一直很忙。我已经更新了我的代码，看看我的答案！似乎是一系列巧合（不知何故在每次尝试的间隔）导致了这个结论。其他几次测试表明使用 as.matrix 可以提高速度。为误导而道歉，尽管这让我更舒服，因为没有理由让它像现在这样工作。

【解决方案2】：

所以我已经为此工作了一段时间并进行了测试。对于遇到类似问题的其他人，这里有两个更优化的代码版本。我已经显着减少了计算时间，但是它仍然会因为太多的数据条目而崩溃。我的下一步是尝试使用 Rcpp 来实现这一点，如果可能的话，利用我可用的 12 个内核（最终目标是在合理的时间范围内计算 1-2 百万个条目）。不确定在任何一点上进行的最佳方式，但这是我的代码。感谢您的帮助！

##################################
##############Optimized code
t.m<-t(test_euclid_log)

knn_log <- function (vec,k) {
  sum(vec[1:k+1])
}
knn_log <- cmpfun(knn_log)

distf <- function(x,t.m) sqrt(colSums((x - t.m)^2))
distf <- cmpfun(distf)

myfunc <- function (tab) {
  rowsums<-numeric(nrow(tab))
  knnsums_log <- matrix(nrow=nrow(tab),ncol=4)
  for(i in 1:nrow(tab)) {
    q<-apply(tab[i,],1,distf,t.m=t.m)
    rowsums[i] <- colSums(q)
    q<-sort(q)
    for (kn in 1:4) {
      knnsums_log[i,kn] <- knn_log(q,kn)             
    }
  }
  return(cbind(rowsums,knnsums_log))
}
myfunc <- cmpfun(myfunc)
system.time(output <- myfunc(t))

我尝试使用 apply：

###############Vectorized
myfuncvec <- function (tab) {
  kn<-c(1:4)
  q<-apply(tab,1,distf,t.m=t.m)
  rowsums <- colSums(q)
  q<-sort(q)
  knnsums_log <- vapply(kn,knn_log,vec=q,FUN.VALUE=c(0))        
  return(c(rowsums,knnsums_log))
}
myfuncvec <- cmpfun(myfuncvec)

t1<-split(t,row(t))
system.time(out <- vapply(t1,myfuncvec,FUN.VALUE=c(0,0,0,0,0)))
out <- t(out)

作为参考，第一个代码似乎更快。

【讨论】：