r - 计算Hawk过程梯度的有效方法

Question

我有兴趣计算以下数量

B(i) = \sum_{j < i}(x_i-x_j)exp^{-\beta(x_i - x_j)}

这是计算 Hawk 过程可能性参数之一的梯度的一部分(更多信息可在此处找到:http://www.ism.ac.jp/editsec/aism/pdf/031_1_0145.pdf)。

Beta 只是问题抖动的一个常数，x_i 是我的第 i 个数据点。

我正在尝试使用以下代码块在 RCPP 中计算上述数量:

for( int i = 1; i< x.size();i++) {
    double temp=0;
    for(int j=0; j<=i-1;j++){
      temp+=(x[i]-x[j])*exp(-beta*(x[i]-x[j]));

    }

但它非常低效且缓慢。关于如何加速这个公式的任何建议？

Answer 1

标准操作在 C++ 中非常快(+、- 等)。然而，exp计算起来更复杂，所以更慢。

因此，如果我们想要提高性能，则更有可能预先计算 exp计算。

在这里，B(i) = \sum_{j < i}(x_i-x_j)exp^{-\beta(x_i - x_j)}相当于B(i) = \sum_{j < i}(x_i-x_j) / exp^{\beta x_i} * exp^{\beta x_j}这样你就可以预先计算 exp仅针对每个索引(并将取决于 i 的索引排除在循环之外)。通过重构它，您可以进行其他预计算。因此，我将之前的两个解决方案放在这里，然后是我的增量解决方案:

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
Rcpp::NumericVector hawk_process_org(Rcpp::NumericVector x, double beta = 3) {

  int n = x.size();
  Rcpp::NumericVector B = Rcpp::no_init( n - 1);

  for (int i = 1; i < n; i++) {

    double temp = 0;

    for (int j = 0; j <= i - 1; j++) {
      temp += (x[i] - x[j]) * exp(-beta * (x[i] - x[j]));
    }

    B(i - 1) = temp;
  }

  return B;
}

// [[Rcpp::export]]
Rcpp::NumericVector hawk_process_cache(Rcpp::NumericVector x, double beta = 3) {

  int n = x.size();
  Rcpp::NumericVector B = Rcpp::no_init( n - 1);

  double x_i;
  for (int i = 1; i < n; ++i) {

    double temp = 0;
    x_i = x[i];

    for (int j = 0; j <= i - 1; ++j) {
      temp += (x_i - x[j]) * 1 / exp(beta * (x_i - x[j]));
    }

    B(i - 1) = temp;
  }

  return B;
}

// [[Rcpp::export]]
Rcpp::NumericVector hawk_process_cache_2(Rcpp::NumericVector x, 
                                         double beta = 3) {

  int i, j, n = x.size();
  Rcpp::NumericVector B(n);
  Rcpp::NumericVector x_exp = exp(beta * x);

  double temp;
  for (i = 1; i < n; i++) {

    temp = 0;
    for (j = 0; j < i; j++) {
      temp += (x[i] - x[j]) * x_exp[j] / x_exp[i];
    }

    B[i] = temp;
  }

  return B;
}

// [[Rcpp::export]]
Rcpp::NumericVector hawk_process_cache_3(Rcpp::NumericVector x, 
                                         double beta = 3) {

  int i, j, n = x.size();
  Rcpp::NumericVector B(n);
  Rcpp::NumericVector x_exp = exp(beta * x);

  double temp;
  for (i = 1; i < n; i++) {

    temp = 0;
    for (j = 0; j < i; j++) {
      temp += (x[i] - x[j]) * x_exp[j];
    }

    B[i] = temp / x_exp[i];
  }

  return B;
}

// [[Rcpp::export]]
Rcpp::NumericVector hawk_process_cache_4(Rcpp::NumericVector x, 
                                         double beta = 3) {

  Rcpp::NumericVector exp_pre = exp(beta * x);
  Rcpp::NumericVector exp_pre_cumsum = cumsum(exp_pre);
  Rcpp::NumericVector x_exp_pre_cumsum = cumsum(x * exp_pre);
  return (x * exp_pre_cumsum - x_exp_pre_cumsum) / exp_pre;
}

// [[Rcpp::export]]
Rcpp::NumericVector hawk_process_cache_5(Rcpp::NumericVector x, 
                                         double beta = 3) {

  int n = x.size();
  NumericVector B(n);
  double exp_pre, exp_pre_cumsum = 0, x_exp_pre_cumsum = 0;

  for (int i = 0; i < n; i++) {
    exp_pre = exp(beta * x[i]);
    exp_pre_cumsum += exp_pre;
    x_exp_pre_cumsum += x[i] * exp_pre;
    B[i] = (x[i] * exp_pre_cumsum - x_exp_pre_cumsum) / exp_pre;
  }

  return B;
}



/*** R
set.seed(111)

x = rnorm(1e3)

all.equal(
  hawk_process_org(x),
  hawk_process_cache(x)
)
all.equal(
  hawk_process_org(x),
  hawk_process_cache_2(x)[-1]
)
all.equal(
  hawk_process_org(x),
  hawk_process_cache_3(x)[-1]
)

all.equal(
  hawk_process_org(x),
  hawk_process_cache_4(x)[-1]
)

all.equal(
  hawk_process_org(x),
  hawk_process_cache_5(x)[-1]
)

microbenchmark::microbenchmark(
  hawk_process_org(x),
  hawk_process_cache(x),
  hawk_process_cache_2(x),
  hawk_process_cache_3(x),
  hawk_process_cache_4(x),
  hawk_process_cache_5(x)
) 
*/

x = rnorm(1e3) 的基准:

Unit: microseconds
                    expr       min         lq        mean     median         uq       max neval   cld
     hawk_process_org(x) 19801.686 20610.0365 21017.89339 20816.1385 21157.4900 25548.042   100    d 
   hawk_process_cache(x) 20506.903 21062.1370 21534.47944 21297.8710 21775.2995 26030.106   100     e
 hawk_process_cache_2(x)  1895.809  2038.0105  2087.20696  2065.8220  2103.0695  3212.874   100   c  
 hawk_process_cache_3(x)   430.084   458.3915   494.09627   474.2840   503.0885  1580.282   100  b   
 hawk_process_cache_4(x)    50.657    55.2930    71.60536    57.6105    63.5700  1190.260   100 a    
 hawk_process_cache_5(x)    43.373    47.0155    60.43775    49.6640    55.6235   842.288   100 a    

这比尝试从小的优化中获得纳秒要有效得多，小的优化可能会使您的代码更难阅读。

---

但是，让我们在我最后的解决方案中尝试@coatless 提出的优化:

// [[Rcpp::export]]
Rcpp::NumericVector hawk_process_cache_6(Rcpp::NumericVector x, 
                                         double beta = 3) {

  int n = x.size();
  NumericVector B = Rcpp::no_init(n);
  double x_i, exp_pre, exp_pre_cumsum = 0, x_exp_pre_cumsum = 0;

  for (int i = 0; i < n; ++i) {
    x_i = x[i];
    exp_pre = exp(beta * x_i);
    exp_pre_cumsum += exp_pre;
    x_exp_pre_cumsum += x_i * exp_pre;
    B[i] = (x_i * exp_pre_cumsum - x_exp_pre_cumsum) / exp_pre;
  }

  return B;
}

x = rnorm(1e6) 的基准:

Unit: milliseconds
                    expr      min       lq     mean   median       uq       max neval cld
 hawk_process_cache_5(x) 42.52886 43.53653 45.28427 44.46688 46.74129  57.38046   100   a
 hawk_process_cache_6(x) 42.14778 43.19054 45.93252 44.28445 46.51052 153.30447   100   a

还是不太有说服力..