Почему поэлементное произведение матрицы намного быстрее с apply?

Question

Я хочу умножить все элементы матрицы вместе.Я могу сделать это с двумя для циклов или с применением.Моя интуиция заключалась в том, что циклы for будут быстрее.Применить необходимо создать временный вектор для хранения результатов произведения строк, а затем применить к нему произведение.

Он все еще должен выполняться для циклов, чтобы умножить все элементы, так что это просто дополнительная операция сохранения промежуточных результатов, которую не должен делать подход для циклов.И все же это примерно в 4 раза быстрее.Почему это?

cols <- 1000
rows <- 1000

a <- matrix(runif(cols * rows, 1, 2), nrow = rows)

system.time({
  result <- 1
  for(i in 1:nrow(a)) {
    for(j in 1:ncol(a)) {
      result <- result * a[i, j]
    }
  }
})
# 0.09s

system.time(result <- prod(apply(a, 1, prod)))
# 0.01s

Answer 1

Для векторизации вместо apply вы можете использовать rowProds из matrixStats:

library(matrixStats)

microbenchmark::microbenchmark({AA = prod(rowProds(a))}, times = 10)

занимает около 18 миллисекунд

Answer 2

Вот что я получил, пытаясь сравнить различные методы.У меня есть некоторые опасения по поводу того, что Inf был результатом многих вычислений, и мне интересно, может ли ограничение в диапазоне 0-1 быть другим.Как и @badmax, я был удивлен, что prod(a) был относительно медленным.Мне показалось, что это должно быть закодировано в C и быть более эффективным.Я также рассуждал, что подход, ориентированный на столбцы, может быть быстрее, чем подход, ориентированный на строки, поскольку именно так хранятся матрицы в R и он был правильным:

library(microbenchmark)
cols <- 1000
rows <- 1000
a <- matrix(runif(cols * rows, 1, 2), nrow = rows)

microbenchmark(loop1 ={
                   result <- 1
                   for(i in 1:nrow(a)) {
                      for(j in 1:ncol(a)) {
                        result <- result * a[i, j]
                         }               } }, 
          loop2 ={result <- 1
                    for(j in 1:ncol(a)) {
                         result <- result * prod(a[ , j])
                         }               },
          loop3 = {
                     result <- 1
                     for(i in 1:nrow(a)) {
                        result <- result * prod( a[i, ])
                         }                }, 
         apply_test = {result <- prod(apply(a, 1, prod))},
         prod_test = {result <- prod(a) },
         Reduce_test = {result <- Reduce("*", a)},
         log_sum = { result<- exp( sum(log(a)))}) #since sum of logs == log of prod
#====================
Unit: milliseconds
        expr        min         lq       mean     median         uq       max neval   cld
       loop1  58.872740  59.782277  60.665321  60.246169  61.156176  67.33558   100   c  
       loop2   5.314437   5.843748   7.316167   6.024948   6.626402  57.36532   100 a    
       loop3   9.614727  10.248335  11.521343  10.541872  10.947829  45.08280   100 ab   
  apply_test   8.336721   8.924148   9.960122   9.166424   9.429118  17.35621   100 ab   
   prod_test  94.314333  95.438939  95.956394  95.911858  96.286444  98.54229   100    d 
 Reduce_test 292.907175 312.754959 389.959756 354.369616 511.151578 545.80829   100     e
     log_sum  19.258281  19.916965  22.333617  20.134510  20.551704 180.18492   100  b   

Я думаю, что результат apply_test по сути делает то же самоекак loop2, возможно, с небольшим дополнительным штрафом за apply_test.Вот результаты для тестового случая, когда диапазон случайных значений был ограничен [0-1] (вместо [1-2]), и они подтверждают мое подозрение, что некоторая разница заключается в обработке значений Inf:

Unit: milliseconds
        expr        min         lq       mean     median         uq        max neval  cld
       loop1  56.693831  58.846847  59.688896  59.448108  60.208619  63.005431   100   c 
       loop2   5.667955   5.907125  10.090634   6.109151   6.532552 183.985617   100 ab  
       loop3   9.533779  10.080330  12.760057  10.431867  10.734991 183.262217   100 ab  
  apply_test   8.144015   8.622861   9.940263   8.904425   9.962390  17.470028   100 ab  
   prod_test   1.327710   1.371449   1.411990   1.394160   1.432646   1.677596   100 a   
 Reduce_test 293.697339 312.384739 385.437743 356.997439 500.446356 557.253762   100    d
     log_sum  22.444015  23.224879  24.064932  23.539085  24.210656  29.774315   100  b  

Функция prod теперь спасена из своего нижнего положения относительно циклов и применяет методы.

Answer 3

Похоже, apply быстрее, потому что продолжается некоторая векторизация.Рассмотрим этот for цикл:

system.time({
  result <- 1
  for(i in 1:nrow(a)) {
      result <- result * prod(a[i,])
    }
  })

, который, для меня, идет так же быстро, как apply.