使用 R 向量化双求和

Question

我正在努力使用 via vectorization 技术将此函数转换为 R：

到目前为止我所能做的是：

c <- matrix(1:9, 3)
z <- 1:3

sum(abs(outer(z, z,"-")) * c)/sum(c)

但我认为它不一定正确。我尝试了一个 for-loop 版本，但它太长了，而且我的回答很可能是错误的。有人热衷于此吗？我错过了什么（或做错了什么）？任何帮助，将不胜感激。

Answer 1

这是一个双循环版本：

q = 
function(z,c){
num = 0
 for(i in 1:length(z)){
  for(j in 1:length(z)){
    num = num + abs(z[i]-z[j]) * c[i,j]
  }
 }
num/sum(c)
}

这是你的矢量化版本，功能化：

q2 = 
function(z,c){sum(c*abs(outer(z,z,'-')) /sum(c))}

对于一个小矩阵来说，它们之间的时间差异不大：

> microbenchmark::microbenchmark(q(z,c), q2(z,c))
Unit: microseconds
     expr    min      lq     mean  median      uq    max neval cld
  q(z, c) 15.368 15.7505 16.59644 16.0225 16.6290 30.346   100   b
 q2(z, c) 12.232 12.8885 13.79178 13.2225 13.6585 44.085   100  a

但对于更大的测试来说，这是一个巨大的胜利：

> c2 = matrix(runif(100*100),100,100)
> z2 = runif(100)
> microbenchmark::microbenchmark(q(z2,c2), q2(z2,c2))
Unit: microseconds
       expr      min       lq       mean   median       uq       max neval cld
  q(z2, c2) 7437.031 7588.131 8046.92272 7794.927 8332.104 10729.799   100   b
 q2(z2, c2)   74.742   78.647   94.20153   86.113  100.125   188.428   100  a 
>

数值差异在浮点公差范围内：

> q(z2,c2) - q2(z2,c2)
[1] 6.661338e-16

所以除非有人有更快的代码，否则我会坚持使用你的代码。

Answer 2

正如@Spacedman 完美解释的那样，您的方法非常有效，但如果您仍想更快，可以尝试 Rcpp :

library(Rcpp)

sourceCpp(code='
#include <Rcpp.h>

// [[Rcpp::export]]
double qRcpp(const Rcpp::NumericVector z, const Rcpp::NumericMatrix cm){
  int zlen = z.length();
  if(!(zlen == cm.nrow() && cm.nrow() == cm.ncol()))
    Rcpp::stop("Invalid sizes");

  double num = 0;
  for(int i = 0 ; i < zlen ; i++){
    for(int j = 0 ; j < zlen ; j++){
      num = num + std::abs(z[i]-z[j]) * cm(i,j);
    }
  } 
  return num / Rcpp::sum(cm);
}

')

基准：

c2 = matrix(runif(100*100),100,100)
z2 = runif(100)
microbenchmark::microbenchmark(q(z2,c2), q2(z2,c2),qRcpp(z2,c2))
# Unit: microseconds
#           expr       min         lq        mean    median         uq       max neval
#      q(z2, c2) 10273.035 10976.3050 11680.85554 11348.763 11765.2010 44115.632   100
#     q2(z2, c2)    64.292    67.9455    80.56427    75.543    86.3565   244.019   100
#  qRcpp(z2, c2)    21.042    21.9180    25.30515    24.256    26.8860    56.403   100

使用 R 向量化双求和

Vectorizing double summations using R

math

r

addition