※ 引述《jackhzt (巴克球)》之铭言：
: [问题类型]:
: 效能咨询(我想让R 跑更快)
: [软件熟悉度]:
: 使用者(已经有用R 做过不少作品)
: [问题叙述]:如何将以下的程式码跑快一点
: [程式范例]:R中的dist这function 因为想使用不同的计算方式，所以希望可以做到
: 和此function表现差不多的function
: 程式码可贴于以下网站：
: https://gist.github.com/anonymous/cf844933bb6858936e25
: 希望有提高效率的方法
1. 如果距离函数是常见的，通常建议用dist达成
2. 如果不常见，尽量考虑用矩阵运算求出来
EX: (以欧式距离来说)
library(magrittr)
x <- matrix(rnorm(50), 5, 10)
distMat <- sweep(-x %*% t(x) * 2, 2, rowSums(x^2), '+') %>%
sweep(1, rowSums(x^2), '+')
diag(distMat) <- 0
distMat %<>% sqrt
all.equal(distMat, as.matrix(dist(x)), check.attributes = FALSE) # TRUE
3. 用简单的平行，我可能会这样做：
library(magrittr)
library(foreach)
library(doSNOW)
library(plyr)
library(Matrix)
dis <- function(x, y) sum(abs(as.numeric(x)-as.numeric(y)))
x <- matrix(rnorm(50), 5, 10)
allCombinations <- combn(1:nrow(x), 2)
cl <- makeCluster(8, type = "SOCK")
registerDoSNOW(cl)
clusterExport(cl, list = c("dis", "x"))
res <- aaply(allCombinations, 2, function(v){
dis(x[v[1],], x[v[2],])
}, .parallel = TRUE)
stopCluster(cl)
distMat <- sparseMatrix(i = allCombinations[1,],
j = allCombinations[2,], x = res) %>%
rbind(0) %>% as.matrix
distMat[lower.tri(distMat)] <- res
4. 或是干脆用RcppArmadillo:
(下面是用之前kernel matrix估计的方法，之前有发过更快的方法...)
library(Rcpp)
library(RcppArmadillo)
## For windows user
# library(inline)
# settings <- getPlugin("Rcpp")
# settings$env$PKG_CXXFLAGS <- paste('-fopenmp', settings$env$PKG_CXXFLAGS)
# settings$env$PKG_LIBS <- paste('-fopenmp -lgomp', settings$env$PKG_LIBS)
# do.call(Sys.setenv, settings$env)
sourceCpp(code = '
// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>
#include <omp.h>
// [[Rcpp::plugins(openmp)]]
using namespace Rcpp;
using namespace arma;
// [[Rcpp::export]]
NumericMatrix kernelMatrix_cpp(NumericMatrix Xr, NumericMatrix Centerr,
double sigma) {
omp_set_num_threads(omp_get_max_threads());
uword n = Xr.nrow(), b = Centerr.nrow(), row_index, col_index;
mat X(Xr.begin(), n, Xr.ncol(), false);
mat Center(Centerr.begin(), b, Centerr.ncol(), false);
mat KerX(n, b);
#pragma omp parallel private(row_index, col_index)
for (row_index = 0; row_index < n; row_index++)
{
#pragma omp for nowait
for (col_index = 0; col_index < b; col_index++)
{
KerX(row_index, col_index) = exp(sum(square(X.row(row_index)
- Center.row(col_index))) / (-2.0 * sigma * sigma));
}
}
return wrap(KerX);
}')