，使用R

Question

向量化双求和

我通过矢量化技术使用翻译此功能为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) 
} 

这是你的向量化版本，functionised：

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