1

我正在尝试尽快计算高斯核评估向量。我在 R^p 中有一个数据点 x，以及一个由 n 个向量 x_i 组成的矩阵 X。我想为每个 x_i 计算 exp( -||x-x_i||^2 / t) 并将结果作为向量返回。

我已经尝试通过以下代码在 R 和 RcppArmadillo 中实现它

代码：

kernel <- function(x, Data, sigma){
  if(sigma <= 0 ) stop('Gaussian kernel parameter <= 0.')

  DiffPart <- (t(t(Data) - x))^2  ## Computes the distance squared of the data and point x
  DiffPart <- rowSums(DiffPart)  # Sum of squares
  exp( - DiffPart / sigma)  #Divide by kernel parameter and evluate exponential function
}

Rcpp犰狳：

arma::Col<double> kernelCPP(arma::Row<double> x, arma::Mat<double> Data, double sigma){
  arma::Mat<double> Diff=Data.each_row()-x;
  int n = Data.n_rows;
  arma::Col<double> kern(n);
  for(int k = 0 ; k < n; k++){
    kern(k) = exp(-arma::accu(square(Diff.row(k)))/sigma);
  }
  return(kern);
}

不幸的是，我的 RcppArmadillo 代码并不比原始 R 代码快多少。我将在未来的代码/计算中计算数十万次内核向量，所以我希望这个过程尽可能快。

在进行微基准测试时，我得到以下结果：

> microbenchmark(
+   kernel(x= TrainX1[1,], Data = TrainX1, sigma = 100)
+ )
Unit: milliseconds
      min       lq   mean   median
  2.223359 2.274559 2.5199 2.308052
       uq     max neval
 2.575144 4.73301   100

和

> microbenchmark(
+   kernelCPP(x= TrainX1[1,], Data = TrainX1, sigma = 100)
+ )
Unit: milliseconds
    min       lq     mean
 1.697706 1.732053 1.826743
   median       uq      max neval
 1.775786 1.871786 2.493439   100

快一点，但不是很多。

4

0 回答 0