主页 / 问题 / 607323
Accepted
mr.T
Asked:2020-12-24 19:01:24 +0000 UTC2020-12-24 19:01:24 +0000 UTC

如何加速 dist() 函数

  • 772

有一个向量“x”和一个矩阵“y”

您需要快速计算“x”与矩阵“y”的每一行之间的欧氏距离

我通过编写自己的函数取代了标准函数“dist()”

штатная

system.time(for(i in 1:nrow(m)) {dist.ve[i] <- dist(rbind(x,m[i,]))})

   user  system elapsed 
   4.38    0.00    4.39

самописная

system.time(for(i in 1:nrow(m)) {dist.ve[i] <- euc.dist(x,m[i,])})
   user  system elapsed 
   0.65    0.00    0.67

但这还不够，我想加速到第二个零 0.0....

你还能想到什么？编码：

x <- rnorm(10)
m <- matrix(data = rnorm(1000000),ncol = 10)

euc.dist <- function(x1, x2) sqrt(sum((x1 - x2) ^ 2))

dist.ve <- rep(0,nrow(m)) # distance vector
system.time(for(i in 1:nrow(m)) {dist.ve[i] <- dist(rbind(x,m[i,]))})
system.time(for(i in 1:nrow(m)) {dist.ve[i] <- euc.dist(x,m[i,])})
r

1 个回答

  1. Best Answer

    您可以通过将函数编译成字节码来加速您的版本。

    您也可以用 C++ 重写代码。在这种情况下，糖功能就足够了Rcpp。

    #include <Rcpp.h>

using namespace Rcpp;

// [[Rcpp::export]]
NumericVector euc_dist3(const NumericVector& x, const NumericMatrix& y) {
    size_t n = y.nrow();
    if (x.size() != y.ncol())
        stop("Length 'x' and ncol 'y' must be equal.");
    NumericVector res = no_init(n);
    for (size_t i = 0; i < n; ++i)
        res[i] = sqrt(sum(pow(x - y.row(i), 2.0)));
    return res;
}

    不使用句法“糖”的等效代码：

    #include <Rcpp.h>

using namespace Rcpp;

// [[Rcpp::export]]
NumericVector euc_dist3(const NumericVector& x, const NumericMatrix& y) {
    size_t n = y.nrow(), m = y.ncol();
    if (x.size() != m)
        stop("Length 'x' and ncol 'y' must be equal.");
    NumericVector res = no_init(n);
    for (size_t i = 0; i < n; ++i) {
        double tmp = 0;
        for (size_t j = 0; j < m; ++j)
            tmp += std::pow(x[j] - y[i + n * j], 2.0);
        res[i] = std::sqrt(tmp);
    }
    return res;
}

    性能比较代码：

    # Данные для сравнения
x <- rnorm(100)
m <- matrix(data = rnorm(1000000), ncol = 100)

euc_dist <- function(x, m) {
    res <- numeric(nrow(m))
    for(i in 1:nrow(m))
        res[i] <- dist(rbind(x ,m[i,]))
    res
}

euc_dist2 <- function(x, m) {
    res <- numeric(nrow(m))
    for(i in seq_len(nrow(m)))
        res[i] <- sqrt(sum((x - m[i, ]) ^ 2))
    res
}

all.equal2 <- function(...) {
    l <- list(...)
    all(sapply(l[-1], all.equal, l[[1]]))
}

# Комплируем функции в байт-код
library(compiler)
euc_dist_c <- cmpfun(euc_dist)
euc_dist2_c <- cmpfun(euc_dist2)

# Убедимся, что функции возвращают одинаковый результат
all.equal2(euc_dist(x, m),
           euc_dist_c(x, m),
           euc_dist2_c(x, m),
           euc_dist2(x, m),
           euc_dist3(x, m))

# Сравниваем производительность функций
library(benchr)
benchmark(euc_dist(x, m),
          euc_dist_c(x, m),
          euc_dist2_c(x, m),
          euc_dist2(x, m),
          euc_dist3(x, m))

    比较结果：

    R> # Убедимся, что функции возвращают одинаковый результат
R> all.equal2(euc_dist(x, m),
..            euc_dist_c(x, m),
..            euc_dist2_c(x,  .... [TRUNCATED] 
[1] TRUE

R> # Сравниваем производительность функций
R> library(benchr)

R> benchmark(euc_dist(x, m),
..           euc_dist_c(x, m),
..           euc_dist2_c(x, m),
..           euc_dist2(x, m),
..           euc_dist3(x, m) .... [TRUNCATED] 
Benchmark summary:
Time units : milliseconds 
             expr n.eval   min  lw.qu median   mean  up.qu    max total relative
   euc_dist(x, m)    100 187.0 191.00 193.00 195.00 194.00 244.00 19500    99.70
 euc_dist_c(x, m)    100 180.0 183.00 184.00 186.00 185.00 237.00 18600    95.10
euc_dist2_c(x, m)    100  15.9  16.50  16.70  17.80  19.80  23.00  1780     8.65
  euc_dist2(x, m)    100  24.9  25.80  28.70  27.90  29.40  31.00  2790    14.80
  euc_dist3(x, m)    100   1.6   1.89   1.93   1.91   1.97   2.09   191     1.00

    如您所见，尽管仅使用了本机 R 代码，但编译版本euc_dist2( ) 本身非常有价值。euc_dist2_c

    如果您需要大大加快速度，可以使用RcppParallel.

    • 3

相关问题