Наиболее эффективный способ вычисления функции с большим количеством комбинаций параметров

Question

Минималистский пример того, что я пытаюсь сделать:

dX_i <- rnorm(100, 0, 0.0002540362)

p_vec <- seq(0, 1, 0.25)  
gamma_vec <- seq(1, 2, 0.25)     
a_vec <- seq(2, 6, 1)
sigma_hat_vec <- c(0.03201636, 0.05771143, 0.07932116, 0.12262327, 0.15074560)
delta_j_vec <- c(0.0000005850109, 0.0000011700217, 0.0000017550326, 0.0000035100651, 0.0000052650977)

parameters <- expand.grid("p" = p_vec, "gamma" = gamma_vec, "a" = a_vec, "sigma_hat" = sigma_hat_vec, "delta_j" = delta_j_vec)


result <- sapply(1:nrow(parameters), function(x) {
  tmp <- parameters[x,]
  p <- tmp$p
  a <- tmp$a
  gamma <- tmp$gamma
  sigma_hat <- tmp$sigma_hat
  delta_j <- tmp$delta_j

  B <- sum( (abs(dX_i)^p) * ( abs(dX_i) < gamma * a * sigma_hat * delta_j^(1/2) ))

  return(B)
})

Цель: мне нужно вычислить B для вектора dX, учитывая все комбинации p, a, gamma, sigma_hat, delta_j.

Однако в действительности сетка parameters имеет ~ 600 тыс. Строк, а dX_i имеет длину ~ 80 тыс. Более того, у меня есть список с ~ 1000 dX_i. Поэтому я хочу сделать этот расчет максимально эффективным. Другие подходы, например, преобразование parameters в data.table и запуск sapply в этом data.table, похоже, не дают ускорения.

Я попытался распараллелить функцию (я ограничен запуском сценария на виртуальной Windows машине):

cl <- makePSOCKcluster(numCores)
num.iter <- 1:nrow(parameters)
parSapply(cl, num.iter, function(x, parameters, dX_i) {
  tmp <- parameters[x,]
  p <- tmp$p
  a <- tmp$a
  gamma <- tmp$gamma
  sigma_hat <- tmp$sigma_hat
  delta_j <- tmp$delta_j
  sum( (abs(dX_i)^p) * ( abs(dX_i) < gamma * a * sigma_hat * delta_j^(1/2) ))
}, parameters, dX_i)
stopCluster(cl)

Хотя это дало мне ускорение, я все еще чувствую, что я ' на самом деле я не решаю эту проблему наиболее эффективным способом и буду признателен за любые предложения.

Answer 1

@ Josliber ответ очень хороший. Тем не менее, это выглядит так, будто R - это плохо ... и вам нужно переключиться на C ++ для повышения производительности.

В ответе реализованы три трюка:

• предварительное вычисление вектор порога
• предварительно вычисляет абсолютное значение dX_i
• сортирует эти значения, чтобы остановить сумму рано

Первые два трюка - просто трюк R называется "векторизация" -> в основном выполнять ваши операции (например, gamma * a * sigma_hat * delta_j^(1/2) или abs()) над целыми векторами, а не над одним элементом внутри al oop.

Это именно то, что вы делаете при использовании sum( dX_i^p * vec_boolean ); он векторизован (* и sum), так что он должен быть очень быстрым.

Если мы реализуем только эти два трюка (мы не можем действительно сделать третий трюк одинаково, потому что он ломается векторизация), это дает:

abs_dX_i <- abs(dX_i)
thresh <- with(parameters, gamma * a * sigma_hat * sqrt(delta_j))
p <- parameters$p
result3 <- sapply(1:nrow(parameters), function(i) {
  in_sum <- (abs_dX_i < thresh[i])
  sum(abs_dX_i[in_sum]^p[i])
})
all.equal(result, result3) # TRUE

Если мы сравним все три решения:

microbenchmark::microbenchmark(
  OP = {
    result <- sapply(1:nrow(parameters), function(x) {
      tmp <- parameters[x,]
      p <- tmp$p
      a <- tmp$a
      gamma <- tmp$gamma
      sigma_hat <- tmp$sigma_hat
      delta_j <- tmp$delta_j

      B <- sum( (abs(dX_i)^p) * ( abs(dX_i) < gamma * a * sigma_hat * delta_j^(1/2) ))

      return(B)
    })
  },
  RCPP = {
    result2 <- proc(sort(abs(dX_i)), parameters$gamma * parameters$a *
                      parameters$sigma_hat * parameters$delta_j^(1/2), parameters$p)
  },
  R_VEC = {
    abs_dX_i <- abs(dX_i)
    thresh <- with(parameters, gamma * a * sigma_hat * sqrt(delta_j))
    p <- parameters$p
    result3 <- sapply(1:nrow(parameters), function(i) {
      in_sum <- (abs_dX_i < thresh[i])
      sum(abs_dX_i[in_sum]^p[i])
    })
  },
  times = 10
)

Мы получим:

Unit: milliseconds
  expr      min       lq      mean   median       uq      max neval
    OP 224.8414 235.4075 289.90096 270.2767 347.1727 399.3262    10
  RCPP  14.8172  15.4691  18.83703  16.3979  20.3829  29.6624    10
 R_VEC  28.3136  29.5964  32.82456  31.4124  33.2542  45.8199    10

Это дает огромное ускорение всего за незначительное изменение исходного кода в R. Это менее чем в два раза медленнее, чем код R cpp, и его можно легко распараллелить, как вы делали ранее с parSapply().

Answer 2

Когда я хочу ускорить трудно векторизованный код, я часто обращаюсь к R cpp. В конце дня вы пытаетесь суммировать abs(dX_i)^p, ограничиваясь значениями abs(dX_i), меньшими порога gamma * a * sigma_hat * delta_j^(1/2). Вы хотите сделать это для пары пар p и порога. Вы можете выполнить sh это с помощью:

library(Rcpp)
cppFunction(
"NumericVector proc(NumericVector dX_i, NumericVector thresh, NumericVector p) {
  const int n = thresh.size();
  const int m = dX_i.size();
  NumericVector B(n);
  for (int i=0; i < n; ++i) {
    B[i] = 0;
    for (int j=0; j < m; ++j) {
      if (dX_i[j] < thresh[i]) {
        B[i] += pow(dX_i[j], p[i]);
      } else {
        break;
      }
    }
  }
  return B;
}"
)
result2 <- proc(sort(abs(dX_i)), parameters$gamma * parameters$a * parameters$sigma_hat * parameters$delta_j^(1/2), parameters$p)
all.equal(result, result2)
# [1] TRUE

Обратите внимание, что мой код сортирует абсолютное значение dX_i, поэтому он может остановить вычисление, как только встретит первое значение, превышающее порог.

На моей машине я вижу 20-кратное ускорение с 0,158 секунды для вашего кода до 0,007 секунды для кода R cpp (измерено с использованием system.time).

Answer 3

Одно наблюдение заключается в том, что у вас есть огромное количество повторов каждого значения p в вашем наборе параметров. Вы можете отдельно обработать каждое значение p; таким образом, вам нужно всего лишь суммировать dX_i, повышенное до определенного значения p один раз.

result4 <- rep(NA, nrow(parameters))
sa_dX_i <- sort(abs(dX_i))
thresh <- parameters$gamma * parameters$a * parameters$sigma_hat * parameters$delta_j^(1/2)
loc <- findInterval(thresh, sa_dX_i)
loc[loc == 0] <- NA  # Handle threshold smaller than everything in dX_i
for (pval in unique(parameters$p)) {
  this.p <- parameters$p == pval
  cs_dX_i_p <- cumsum(sa_dX_i^pval)
  result4[this.p] <- cs_dX_i_p[loc[this.p]]
}
result4[is.na(result4)] <- 0  # Handle threshold smaller than everything in dX_i
all.equal(result, result4)
# [1] TRUE

Чтобы увидеть это в действии, давайте увеличим исходный набор данных до того, что описано в вопросе. (~ 600 тыс. Строк параметров и ~ 80 тыс. Значений в dX_i):

set.seed(144)
dX_i <- rnorm(80000, 0, 0.0002540362)
p_vec <- seq(0, 1, 0.025)  
gamma_vec <- seq(1, 2, 0.025)     
a_vec <- seq(2, 6, 0.3)
sigma_hat_vec <- c(0.03201636, 0.05771143, 0.07932116, 0.12262327, 0.15074560)
delta_j_vec <- c(0.0000005850109, 0.0000011700217, 0.0000017550326, 0.0000035100651, 0.0000052650977)
parameters <- expand.grid("p" = p_vec, "gamma" = gamma_vec, "a" = a_vec, "sigma_hat" = sigma_hat_vec, "delta_j" = delta_j_vec)
dim(parameters)
# [1] 588350      5
length(unique(parameters$p))
# [1] 41

Ускорение довольно драматично c - на моем компьютере этот код занимает 0,27 секунды, а R cpp Код, размещенный в моем другом ответе на этот вопрос, занимает 655 секунд (ускорение 2400x, используя чистый R!). Очевидно, что это ускорение работает только в том случае, если в кадре данных parameters относительно мало p значений (каждое повторяется много раз). Если каждое значение p уникально, это, вероятно, будет намного медленнее, чем другие предложенные подходы.