Как улучшить время обработки для расчета евклидова расстояния - PullRequest
Новая
       13

Как улучшить время обработки для расчета евклидова расстояния

1 голос
/ 09 мая 2019

Я пытаюсь вычислить взвешенное евклидово расстояние (в квадрате) между двумя кадрами данных, которые имеют одинаковое количество столбцов (переменных) и разное количество строк (наблюдений).

Расчет производится по формуле: 

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

Мне специально нужно умножить каждую посылку соматика на определенный вес (лямбда).

Приведенный ниже код работает правильно, но если я использую его в сотнях итерацийэто занимает много времени обработки.Вчера мне потребовалось 18 часов, чтобы создать графику, используя несколько итераций функции, содержащей этот расчет.Используя библиотеку (profvis) profvis ({мой код}), я увидел, что эта конкретная часть кода занимает примерно 80% времени обработки.

Я много читал о том, как сократить время обработки, используяпараллельные и векторизованные операции, но я не знаю, как реализовать их в данном конкретном случае, из-за веса lamb #.

Может ли кто-нибудь помочь мне сократить время обработки с помощью этого кода?

Более подробную информацию о коде и структуре данных можно найти в приведенном ниже коде в виде комментариев. 

# 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 (для длины (лямбда> 1) и @F. Privé (для длины (лямбда == 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)
  }

Ответы [ 2 ]

4 голосов
/ 10 мая 2019

Перепишите проблему, используя линейную алгебру и векторизацию, которая намного быстрее, чем циклы.

Если у вас нет 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, '*'))
1 голос
/ 10 мая 2019

Здесь альтернативный способ использования 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
...