RCPP：我的距离矩阵程序比包中的函数慢

Question

我想计算成对的欧式距离矩阵。我根据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，有太多的东西需要被检查（如method，diag）。

Answer 1

你快到了。但是你的内部循环体试图在一行中做太多事情。模板编程已经够难了，有时最好将指令分散一点，给编译器一个更好的机会。所以我只做了五个语句，并立即构建了它。

新代码：

#include <Rcpp.h>

using namespace Rcpp;

double dist1 (NumericVector x, NumericVector y){
  int n = y.length();
  double total = 0;
  for (int i = 0; i < n ; ++i) {
    total += pow(x(i)-y(i),2.0);
  }
  total = sqrt(total);
  return total;
}

// [[Rcpp::export]]
NumericMatrix calcPWD (NumericMatrix x){
  int outrows = x.nrow();
  int outcols = x.nrow();
  NumericMatrix out(outrows,outcols);

  for (int i = 0 ; i < outrows - 1; i++){
    for (int j = i + 1  ; j < outcols ; j ++) {
      NumericVector v1 = x.row(i);
      NumericVector v2 = x.row(j-1);
      double d = dist1(v1, v2);
      out(j-1,i) = d;
      out(i,j-1)= d;
    }
  }
  return (out) ;
}

/*** R
M <- matrix(log(1:9), 3, 3)
calcPWD(M)
*/

运行它：

R> sourceCpp("/tmp/mikebrown.cpp")

R> M <- matrix(log(1:9), 3, 3)

R> calcPWD(M)
         [,1]     [,2] [,3]
[1,] 0.000000 0.740322    0
[2,] 0.740322 0.000000    0
[3,] 0.000000 0.000000    0
R> 

不过，您可能想检查索引逻辑。看起来你错过了更多的比较。

编辑：对于踢球，这是距离函数的更紧凑版本：

// [[Rcpp::export]]
double dist2(NumericVector x, NumericVector y){
  double d = sqrt( sum( pow(x - y, 2) ) );
  return d;
}

Answer 2

Rcpp与内部R函数（C / Fortran）

首先，仅因为您使用Rcpp编写算法并不一定意味着它将击败R等效项，尤其是如果R函数调用C或Fortran例程来执行大量计算时。在其他情况下，仅用R编写函数，则很有可能在Rcpp中对其进行转换会产生所需的速度增益。

请记住，在重写内部函数时，将与绝对疯狂的C程序员的R Core团队抗衡。

基本实施 dist()

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

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

，这是dist()函数R源的最后一行。与C ++相比，它具有一点优势，因为它更精细而不是模板化。

此外，如果可用，C实现将使用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;
}