如何提高欧几里得距离计算的处理时间答案

【问题标题】：How to improve processing time for euclidean distance calculation如何提高欧几里得距离计算的处理时间
【发布时间】：2019-09-27 16:16:17
【问题描述】：

我正在尝试计算具有相同列数（变量）和不同行数（观察值）的两个数据帧之间的加权欧几里得距离（平方）。

计算遵循公式：

DIST[m,i] <- sum(((DATA1[m,] - DATA2[i,]) ^ 2) * lambda[1,])

我特别需要将体细胞的每个部分乘以特定的重量 (lambda)。

下面提供的代码运行正确，但如果我在数百次迭代中使用它，则需要大量处理时间。昨天，我花了 18 个小时，使用包含此计算的函数的多次迭代来创建图形。使用 library(profvis) profvis({ my code }) 我看到代码的这个特定部分占用了大约 80% 的处理时间。

我阅读了很多关于如何使用并行和矢量化操作来减少处理时间的文章，但我不知道如何在这种特殊情况下实现它们，因为它很重#。

有人可以帮我缩短这段代码的处理时间吗？

有关代码和数据结构的更多信息可以在下面作为 cmets 提供的代码中找到。

# Data frames used to calculate the euclidean distances between each observation 
#   from DATA1 and each observation from DATA2.
# The euclidean distance is between a [600x50] and a [8X50] dataframes, resulting 
#   in a [600X8] dataframe.
DATA1 <- matrix(rexp(30000, rate=.1), ncol=50) #[600x50]
DATA2 <- matrix(rexp(400, rate=.1), ncol=50) #[8X50]

# Weights used for each of the 50 variables to calculate the weighted 
#   euclidean distance.
# Can be a vector of different weights or a scalar of the same weight 
#   for all variables.
lambda <- runif(n=50, min=0, max=10)   ## length(lambda) > 1
# lambda=1   ## length(lambda) == 1

if (length(lambda) > 1) {
  as.numeric(unlist(lambda))
  lambda <- as.matrix(lambda)
  lambda <- t(lambda)
}

nrows1 <- nrow(DATA1)
nrows2 <- nrow(DATA2)

# Euclidean Distance calculation
DIST <- matrix(NA, nrow=nrows1, ncol=nrows2 )  
for (m in 1:nrows1) {
  for (i in 1:nrows2) {
    if (length(lambda) == 1) { 
      DIST[m, i] <- sum((DATA1[m, ] - DATA2[i, ])^2) 
    }
    if (length(lambda) > 1){ 
      DIST[m, i] <- sum(((DATA1[m, ] - DATA2[i, ])^2) * lambda[1, ])
    }
    next
  }
  next
}

在所有建议之后，结合@MDWITT（对于长度（lambda > 1）和@F.Privé（对于长度（lambda == 1））的答案，最终解决方案只需要一分钟即可运行，而原始解决方案我花了一个半小时来运行一个具有该计算的更大代码。对于那些感兴趣的人来说，这个问题的最终代码是：

#Data frames used to calculate the euclidean distances between each observation from DATA1 and each observation from DATA2.
#The euclidean distance is between a [600x50] and a [8X50] dataframes, resulting in a [600X8] dataframe.
DATA1 <- matrix(rexp(30000, rate=.1), ncol=50) #[600x50]
DATA2 <- matrix(rexp(400, rate=.1), ncol=50) #[8X50]

#Weights used for each of the 50 variables to calculate the weighted euclidean distance.
#Can be a vector of different weights or a scalar of the same weight for all variables.
#lambda <- runif(n = 50, min = 0, max = 10)   ##length(lambda) > 1
lambda = 1   ##length(lambda) == 1

nrows1 <- nrow(DATA1)
nrows2 <- nrow(DATA2) 

#Euclidean Distance calculation
DIST <- matrix(NA, nrow = nrows1, ncol = nrows2)  

if (length(lambda) > 1){
  as.numeric(unlist(lambda))
  lambda <- as.matrix(lambda)
  lambda <- t(lambda)

  library(Rcpp)
  cppFunction('NumericMatrix weighted_distance (NumericMatrix x, NumericMatrix y, NumericVector lambda){

              int n_x = x.nrow();
              int n_y = y.nrow();


              NumericMatrix DIST(n_x, n_y);

              //begin the loop

              for (int i = 0 ; i < n_x; i++){
              for (int j = 0  ; j < n_y ; j ++) {
              double d = sum(pow(x.row(i) - y.row(j), 2)*lambda);
              DIST(i,j) = d;
              }
              }
              return (DIST) ;
  }')

    DIST <- weighted_distance(DATA1, DATA2, lambda = lambda)}


  if (length(lambda) == 1) { 
    DIST <- outer(rowSums(DATA1^2), rowSums(DATA2^2), '+') - tcrossprod(DATA1, 2 * DATA2)
  }

【问题讨论】：

可能不会更快，但为什么将 lambda 存储在单行矩阵中。只需存储在向量中并避免子集lambda[1,]。另外，你代码中的next 语句是没用的
可能也值得一试stackoverflow.com/questions/36829700/…>这将是向 Rcpp 的转移，但您需要添加 Lamba
那么，我会在length(lambda)==1时设置lambda = 1，完全去掉这两个if语句，这样就可以在第二个循环中直接使用一行公式
您的代码中的索引和大小有误。请修复它。
@MDEWITT 我会调查的！谢谢！

标签： r performance euclidean-distance

【解决方案1】：

这里使用Rcpp 的另一种方式只是为了拥有这个概念文档。在一个名为 euclidean.cpp 的文件中，我有

#include <Rcpp.h>
#include <cmath>

using namespace Rcpp;

// [[Rcpp::export]]

NumericMatrix weighted_distance (NumericMatrix x, NumericMatrix y, NumericVector lambda){

  int n_x = x.nrow();
  int n_y = y.nrow();


  NumericMatrix out(n_x, n_y);

  //begin the loop

  for (int i = 0 ; i < n_x; i++){
    for (int j = 0  ; j < n_y ; j ++) {
      double d = sum(pow(x.row(i) - y.row(j), 2)*lambda);
      out(i,j) = d;
    }
  }
  return (out) ;
}

在 R 中，我有

library(Rcpp)
sourceCpp("libs/euclidean.cpp")

# Generate Data
DATA1 <- matrix(rexp(30000, rate=.1), ncol=50) #[600x50]
DATA2 <- matrix(rexp(400, rate=.1), ncol=50) #[8X50]
lambda <- runif(n=50, min=0, max=10)

# Run the program

out <- weighted_distance(DATA1, DATA2, lambda = lambda)

当我使用以下方法测试速度时：

microbenchmark(
  Rcpp_way = weighted_distance(DATA1, DATA2, lambda = lambda),
other = {DIST <- matrix(NA, nrow=nrows1, ncol=ncols)  
for (m in 1:nrows1) {
  for (i in 1:nrows2) {
    if (length(lambda) == 1) { 
      DIST[m, i] <- sum((DATA1[m, ] - DATA2[i, ])^2) 
    }
    if (length(lambda) > 1){ 
      DIST[m, i] <- sum(((DATA1[m, ] - DATA2[i, ])^2) * lambda[1, ])
    }
    next
  }
  next
}}, times = 100)

你可以看到它是一个更快的剪辑：

Unit: microseconds
     expr       min        lq       mean    median         uq        max neval
 Rcpp_way   446.769   492.308   656.9849   562.667   846.9745   1169.231   100
    other 24688.821 30681.641 44153.5264 37511.385 50878.3585 200843.898   100

【讨论】：

非常感谢您的回复。我在 R 中运行它： library(Rcpp) cppFunction('NumericMatrix weighted_distance ... return (out) ;}') out
您的功能似乎没问题，但结果却不是，正如我上面所说。不知道怎么解决。
@AnaF.updated。我可以运行您的原始函数和更新后的 Rcpp 函数，并在 600 x 8 矩阵中得到完全相同的答案
现在完美运行！！问题解决了！！！我只花了一分钟来运行一个需要一个半小时才能运行的代码。
@AnaF。很高兴听见！当你有很多循环时，Rcpp 是一个很大的提升

【解决方案2】：

重写问题以使用线性代数和向量化，这比循环快得多。

如果你没有lambda，这只是

outer(rowSums(DATA1^2), rowSums(DATA2^2), '+') - tcrossprod(DATA1, 2 * DATA2)

有了lambda，就变成了

outer(drop(DATA1^2 %*% lambda), drop(DATA2^2 %*% lambda), '+') -
    tcrossprod(DATA1, sweep(DATA2, 2, 2 * lambda, '*'))

【讨论】：

lambda 解决方案导致Error in DATA1^2 %*% lambda : non-conformable arguments。非 lambda 解决方案比 OP 循环快 50 倍。
@Cole 是对的。 lambda 解决方案导致 Error in DATA1^2 %*% lambda : non-conformable arguments，但另一行正常工作。
使用lambda <- runif(n=50, min=0, max=10)，我没有收到任何错误。