c – Rcpp：我的距离矩阵程序比包中的函数慢

我想计算成对的欧几里德距离矩阵.我根据Dirk Eddelbuettel的建议写了Rcpp程序,如下所示

NumericMatrix calcPWD1 (NumericMatrix x){
  int outrows = x.nrow();
  double d;
  NumericMatrix out(outrows,outrows);

  for (int i = 0 ; i < outrows - 1; i++){
    for (int j = i + 1  ; j < outrows ; j ++){
      NumericVector v1= x.row(i);
      NumericVector v2= x.row(j);
      NumericVector v3=v1-v2;
      d = sqrt(sum(pow(v3,2)));
      out(j,i)=d;
      out(i,j)=d;
    }
  }
  return (out) ;
}

但我发现我的程序比dist函数慢.

> benchmark(as.matrix(dist(b)),calcPWD1(b))
                test replications elapsed relative user.self sys.self user.child sys.child
1 as.matrix(dist(b))          100  24.831    1.000    24.679    0.010          0         0
2        calcPWD1(b)          100  27.362    1.102    27.346    0.007          0         0

你们有什么建议吗？我的矩阵非常简单.没有列名或行名,只有普通矩阵(例如b = matrix(c(rnorm(1000 * 10)),1000,10)).
这是dist的程序

> dist
function (x, method = "euclidean", diag = FALSE, upper = FALSE, 
    p = 2) 
{
    if (!is.na(pmatch(method, "euclidian"))) 
        method <- "euclidean"
    METHODS <- c("euclidean", "maximum", "manhattan", "canberra", 
        "binary", "minkowski")
    method <- pmatch(method, METHODS)
    if (is.na(method)) 
        stop("invalid distance method")
    if (method == -1) 
        stop("ambiguous distance method")
    x <- as.matrix(x)
    N <- nrow(x)
    attrs <- if (method == 6L) 
        list(Size = N, Labels = dimnames(x)[[1L]], Diag = diag, 
            Upper = upper, method = METHODS[method], p = p, call = match.call(), 
            class = "dist")
    else list(Size = N, Labels = dimnames(x)[[1L]], Diag = diag, 
        Upper = upper, method = METHODS[method], call = match.call(), 
        class = "dist")
    .Call(C_Cdist, x, method, attrs, p)
}
<bytecode: 0x56b0d40>
<environment: namespace:stats>

我希望我的程序比dist快,因为在dist中,有太多东西需要检查(比如方法,diag).

解决方法:

Rcpp与内部R函数(C / Fortran)

首先,仅仅因为你使用Rcpp编写算法并不一定意味着它会击败R等价物,特别是如果R函数调用C或Fortran例程来执行大部分计算.在函数纯粹用R编写的其他情况下,很有可能在Rcpp中转换它将产生所需的速度增益.

请记住,当重写内部函数时,一个人正在面对R Core团队的绝对疯狂的C程序员很可能会胜出.

dist()的基本实现

其次,距离计算R使用在C中完成,如下所示：

.Call(C_Cdist, x, method, attrs, p)

,这是dist()函数的R源的最后一行.这使得它与C相比略有优势,因为它更精细而不是模板化.

此外,C implementation uses OpenMP可用于并行化计算.

建议的修改

第三,通过稍微改变子集顺序并避免创建附加变量,版本之间的时间减少.

#include <Rcpp.h>

// [[Rcpp::export]]
Rcpp::NumericMatrix calcPWD1 (const Rcpp::NumericMatrix & x){
  unsigned int outrows = x.nrow(), i = 0, j = 0;
  double d;
  Rcpp::NumericMatrix out(outrows,outrows);

  for (i = 0; i < outrows - 1; i++){
    Rcpp::NumericVector v1 = x.row(i);
    for (j = i + 1; j < outrows ; j ++){
      d = sqrt(sum(pow(v1-x.row(j), 2.0)));
      out(j,i)=d;
      out(i,j)=d;
    }
  }

  return out;
}

标签：c,r,rcpp
来源： https://codeday.me/bug/20191002/1843704.html