Эффективный способ расчета градиента процесса Хоука

Question

Я заинтересован в вычислении следующего количества

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

, которое является частью вычисления градиента по одному из параметров вероятности процесса Хоука (более подробную информацию можно найти здесь: http://www.ism.ac.jp/editsec/aism/pdf/031_1_0145.pdf).

Бета - это всего лишь константа для решения проблемы, а 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

Все еще не очень убедительно ..

Answer 2

Интересный вопрос.В моих тестах объединение двух ответов действительно дает дополнительное повышение производительности (тесты еще ниже):

#include <Rcpp.h>
using namespace Rcpp;

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

  int n = x.size();
  NumericVector B = Rcpp::no_init(n-1);
  double exp_pre(exp(beta * x[0]));
  double exp_pre_cumsum(exp_pre);
  double x_exp_pre_cumsum(x[0] * exp_pre);
  double x_i;

  for (int i = 1; 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-1] = (x_i * exp_pre_cumsum - x_exp_pre_cumsum) / exp_pre;
  }

  return B;
}

all.equal(
  hawk_process_org(x),
  hawk_process_cache_combined(x)
)
#> [1] TRUE

Теперь, хотя исходная формулировка «смущающе параллельна», это больше не относится к этому выражению.Однако алгоритмы сканирования префиксов, такие как cumsum, также могут быть распараллелены.А библиотеки типа ArrayFire предоставляют интерфейсы для таких алгоритмов с использованием графического процессора.Используя RcppArrayFire можно писать на основе Ф. Приве hawk_process_cached_4:

// [[Rcpp::depends(RcppArrayFire)]]
#include <RcppArrayFire.h>
// [[Rcpp::export]]
af::array hawk_process_af(RcppArrayFire::typed_array<f32> x, 
                          double beta = 3) {

  af::array exp_pre = exp(beta * x);
  af::array exp_pre_cumsum = af::accum(exp_pre);
  af::array x_exp_pre_cumsum = af::accum(x * exp_pre);
  af::array result = (x * exp_pre_cumsum - x_exp_pre_cumsum) / exp_pre;
  return result(af::seq(1, af::end));
}

Здесь результаты не совсем равны, поскольку мой драйвер / карта поддерживает только операции с плавающей запятой одинарной точности:

all.equal(
  hawk_process_org(x),
  hawk_process_af(x)
)
#> [1] "Mean relative difference: 3.437819e-07"

С двойной точностью написали бы f64 выше и получили бы идентичные результаты.Теперь для тестов:

set.seed(42)
x <- rnorm(1e3)
microbenchmark::microbenchmark(
  hawk_process_af(x),
  hawk_process_cache_combined(x),
  hawk_process_cache_5(x)[-1]
) 
#> Unit: microseconds
#>                            expr     min       lq      mean   median      uq      max neval
#>              hawk_process_af(x) 245.281 277.4625 338.92232 298.5410 346.576 1030.045   100
#>  hawk_process_cache_combined(x)  35.343  39.0120  43.69496  40.7770  45.264    84.242   100
#>     hawk_process_cache_5(x)[-1]  52.408  57.8580  65.55799  60.5265  67.965  125.864   100
x <- rnorm(1e6)
microbenchmark::microbenchmark(
  hawk_process_af(x),
  hawk_process_cache_combined(x),
  hawk_process_cache_5(x)[-1]
)
#> Unit: milliseconds
#>                            expr      min       lq     mean   median       uq       max neval
#>              hawk_process_af(x) 27.54936 28.42794 30.93452 29.20025 32.40667  49.41888   100
#>  hawk_process_cache_combined(x) 34.00380 36.84497 40.74862 39.03649 41.85902 111.51628   100
#>     hawk_process_cache_5(x)[-1] 47.02501 53.24702 57.94747 55.35018 58.42097 130.89737   100  

Таким образом, для небольших векторов комбинированный подход быстрее, а при более длительной разгрузке на GPU окупается.Все это не с мощным графическим процессором, а с простой встроенной графикой:

RcppArrayFire::arrayfire_info()
#> ArrayFire v3.5.1 (OpenCL, 64-bit Linux, build 0a675e8)
#> [0] BEIGNET: Intel(R) HD Graphics Skylake ULT GT2, 4096 MB

Answer 3

Это операция O (N ^ 2) без учета стоимости опыта.Любые изменения могут привести к минимальным улучшениям.

Несколько быстрых предложений:

• кеш значение x[i] во внешнем цикле, так как вы неоднократно поднастраиваете это во внутреннем цикле.
• переключитесь с использования exp(-beta * ..) на 1/exp(beta*(x ... ))
• используйте ++i вместо i++, чтобы избежать небольшого сбоя производительности , поскольку вы избегаете копирования iчто последний делает.

---

Оригинальный код:

#include<Rcpp.h>

// [[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;
}

Модифицированный код:

#include<Rcpp.h>

// [[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;
}

Тест

set.seed(111)

x = rnorm(1e4)

all.equal(
  hawk_process_org(x),
  hawk_process_cache(x)
)
#> [1] TRUE

bench_func = microbenchmark::microbenchmark(
  hawk_process_org(x),
  hawk_process_cache(x)
)

bench_func
#> Unit:milliseconds
#>                 expr      min       lq     mean   median       uq      max   neval
#> hawk_process_org(x)   436.5349 465.9674 505.9606 481.4703 500.6652 894.7477   100
#> hawk_process_cache(x) 446.0499 454.9098 485.3830 468.6580 494.9457 799.0940   100

Таким образом, вы получите несколько лучшие результаты по рекомендациям.