Как увеличить цикл - PullRequest
Новая
       11

Как увеличить цикл

2 голосов
/ 13 июня 2019

Я страдаю от медленного выполнения цикла for в R. Здесь я предоставляю часть моего кода, которая вызывает задержку. 

## subsitutes for original data
DC <- matrix(rnorm(10), ncol=101, nrow=6400)
C <- matrix(rnorm(20), ncol=101, nrow=6400)


N <- 80
Vcut <- ncol(DC) 
V <- seq(-2.9,2.5,length=Vcut)
fNC <- matrix(NA, nrow=(N*N), ncol=Vcut)
fNDC <- matrix(NA, nrow=(N*N), ncol=Vcut)


Arbfunc <- function(dV){

b <- matrix(NA, nrow=1, ncol=Vcut)

  for(i in 1:(N*N)) {
    for (n in 1:Vcut) {
      for (k in 1:Vcut) {
        b[k] = (V[2]-V[1])*(exp((-1)*abs(V[k])))*exp(abs(V[n]-V[k])/dV)*(C[i,k]/V[k])
      }
      fNC[i,n] = exp(1*abs(V[n]))*(1/(2*dV))*(sum(b[]))
      fNDC[i,n] = DC[i,n]/fNC[i,n]
    }
  }   
}

Arbfunc(0.5)

Так как мне нужно сравнить результаты среди различных значений dV's , этот код должен быть запущен как минимум в течение нескольких секунд.Но результат - 

user   system  elapsed
40.15   0.03   40.24

, что слишком медленно для достаточного сравнения.Я попробовал несколько методов распараллеливания, но результат не был удовлетворительным (40 -> 25 секунд, хотя я использовал 11 потоков на моем компьютере).

Поэтому я предполагаю, что узким местом является сам цикл for, а не непараллельный код.Не могли бы вы дать мне несколько советов по улучшению этого цикла for или подсказку для распараллеливания?Был бы благодарен за короткий комментарий.

1 Ответ

6 голосов
/ 13 июня 2019

Большое спасибо @Mikko Marttila за исправление функций 3 и 4 и предоставление идеи для функции 5.

R лучше всего подходит с векторизованными параметрами вместо явных циклов. Например, внутренний цикл с k: 

for (k in 1:Vcut) {
  b[k] = (V[2]-V[1])*(exp((-1)*abs(V[k])))*exp(abs(V[n]-V[k])/dV)*(C[i,k]/V[k])
}

Это то же самое, что сказать 

(V[2]-V[1])*(exp((-1)*abs(V)))*exp(abs(V[n]-V)/dV)*(C[i,]/V)

Это небольшое изменение дает нам увеличение производительности в 500 раз для этой части функции: 

Unit: microseconds
         expr     min      lq      mean  median      uq     max neval
       k_loop 13186.7 13603.2 14605.471 13832.9 14517.8 41935.1   100
 k_vectorized    16.4    17.6    25.559    28.8    32.0    52.7   100

Теперь, если мы посмотрим на внешний цикл с i, мы увидим, что действительно нет необходимости циклически повторять каждую строку. Вместо этого мы могли бы создать матрицу для оператора sum(b[k]), превращающего это: 

(V[2]-V[1])*(exp((-1)*abs(V)))*exp(abs(V[n]-V)/dV)*(C[i,]/V)

В это: 

(V[2]-V[1])*(exp((-1)*abs(V)))*exp(abs(V[n]-V)/dV)*(t(C)/V)

Это только что спасло нас N*N*k петли. В вашем случае это 646400 петель.

В целом, у нас будет: 

Arbfunc3 <- function(dV){
    for (n in 1:Vcut) {
      sum_b = colSums((V[2]-V[1])*(exp((-1)*abs(V)))*exp(abs(V[n]-V)/dV)*(t(C)/V))
      fNC[, n] = exp(1*abs(V[n]))*(1/(2*dV))*(sum_b)
      fNDC[, n] = DC[,n]/fNC[,n]
    }
}

Мое среднее время микробенчмарка для этой альтернативы составляет 750 миллисекунд.

Для дальнейшего повышения производительности нам нужно обратиться к V[n] - V. К счастью, у R есть функция - outer(V, V, '-'), и это даст матрицу со всеми необходимыми комбинациями. 

Arbfunc4 <- function(dV) {
  sum_b = apply((V[2]-V[1])*(exp((-1)*abs(V)))*exp(abs(outer(V, V, '-')) / dV) / V, 2, function(x) colSums(x * t(C)))

  fNC = exp(1*abs(V))*(1/(2*dV))*t(sum_b)
  fNDC= DC/t(fNC)
  fNDC
}

Спасибо @Mikko Marttila за предложение избавиться от применения с точечным продуктом. 

Arbfunc5 <- function(dV) {
  a = (V[2] - V[1]) * exp(-abs(V)) * t(C) / V
  b = exp(abs(outer(V, V, "-")) / dV) %*% a

  fNC = exp(1*abs(V))*(1/(2*dV))*(b)
  fNDC= DC/t(fNC)
  fNDC
}

Вот system.time для каждого решения (Arbfunc2 - это исключение k_loop). Оптимизированное решение в 2600 раз быстрее оригинального. 

> system.time(Arbfunc(0.5))
   user  system elapsed 
  78.03    0.39   79.72 
> system.time(Arbfunc2(0.5))
   user  system elapsed 
  10.41    0.03   10.46 
> system.time(Arbfunc3(0.5))
   user  system elapsed 
   0.69    0.13    0.81 
> system.time(Arbfunc4(0.5))
   user  system elapsed 
   0.43    0.05    0.47 
> system.time(Arbfunc5(0.5))
   user  system elapsed 
   0.03    0.00    0.03

Окончательное редактирование: Вот полный код, который я запустил после перезапуска R и очистки моего окружения. Нет ошибок: 

## subsitutes for original data
DC <- matrix(rnorm(10), ncol=101, nrow=6400)
C <- matrix(rnorm(20), ncol=101, nrow=6400)

N <- 80
Vcut <- ncol(DC) 
V <- seq(-2.9,2.5,length=Vcut)

# Unneeded for Arbfunc4 adn Arbfunc5
# Corrected from NA to NA_real_ to prevent coercion from logical to numeric
# h/t to @HenrikB
fNC <- matrix(NA_real_, nrow=(N*N), ncol=Vcut)
fNDC <- matrix(NA_real_, nrow=(N*N), ncol=Vcut)

Arbfunc <- function(dV){
  b <- matrix(NA, nrow=1, ncol=Vcut)

  for(i in 1:(N*N)) {
    for (n in 1:Vcut) {
      for (k in 1:Vcut) {
        b[k] = (V[2]-V[1])*(exp((-1)*abs(V[k])))*exp(abs(V[n]-V[k])/dV)*(C[i,k]/V[k])
      }
      fNC[i,n] = exp(1*abs(V[n]))*(1/(2*dV))*(sum(b[]))
      fNDC[i,n] = DC[i,n]/fNC[i,n]
    }
  }
  fNDC
}

Arbfunc2 <- function(dV){
  b <- matrix(NA, nrow=1, ncol=Vcut)

  for(i in 1:(N*N)) {
    for (n in 1:Vcut) {
      sum_b = sum((V[2]-V[1])*(exp((-1)*abs(V)))*exp(abs(V[n]-V)/dV)*(C[i,]/V))
      fNC[i,n] = exp(1*abs(V[n]))*(1/(2*dV))*(sum_b)
      fNDC[i,n] = DC[i,n]/fNC[i,n]
    }
  }
  fNDC
}

Arbfunc3 <- function(dV){
  for (n in 1:Vcut) {
    sum_b = colSums((V[2]-V[1])*(exp((-1)*abs(V)))*exp(abs(V[n]-V)/dV)*(t(C)/V))
    fNC[, n] = exp(1*abs(V[n]))*(1/(2*dV))*(sum_b)
    fNDC[, n] = DC[,n]/fNC[,n]
  }
  fNDC
}

Arbfunc4 <- function(dV) {
  sum_b = apply((V[2]-V[1])*(exp((-1)*abs(V)))*exp(abs(outer(V, V, '-')) / dV) / V, 2, function(x) colSums(x * t(C)))

  fNC = exp(1*abs(V))*(1/(2*dV))*t(sum_b)
  DC/t(fNC)
}

Arbfunc5 <- function(dV) {
#h/t to Mikko Marttila for dot product
  a = (V[2] - V[1]) * exp(-abs(V)) * t(C) / V
  b = exp(abs(outer(V, V, "-")) / dV) %*% a

  fNC = exp(1*abs(V))*(1/(2*dV))*(b)
  DC/t(fNC)
}

#system.time(res <- Arbfunc(0.5))
system.time(res2 <- Arbfunc2(0.5))
system.time(res3 <- Arbfunc3(0.5))
system.time(res4 <- Arbfunc4(0.5))
system.time(res5 <- Arbfunc5(0.5))

all.equal(res2,res3,res4,res5)

Как упоминает @HenrikB, fNC и fNDC инициализируются как логические матрицы. Это означает, что мы получаем снижение производительности при приведении их к real матрицам. Неправильное выполнение - единовременное попадание в 1 мс для этого набора данных, но если бы это приведение было в цикле, оно могло бы действительно сложиться. 

mat_NA_real_ <- function() {
  mat = matrix(NA_real_, nrow = 6400, ncol = 101)
  mat[1,1] = 1
}

mat_NA <- function() {
  mat = matrix(NA, nrow = 6400, ncol = 101)
  mat[1,1] = 1
}
microbenchmark(mat_NA_real_(), mat_NA())

Unit: microseconds
           expr    min      lq     mean  median     uq     max neval
 mat_NA_real_()  979.5  992.25 1490.081  998.65 1021.1  7612.5   100
       mat_NA() 1865.8 1883.30 3793.119 1911.30 5335.4 53635.2   100
Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.

Нет похожих вопросов

...