Применить функцию к массиву, возвращающему исходное количество измерений

#r #arrays #apply

#r #массивы #применить

Вопрос:

Возьмем этот пример массива:

 set.seed(1)
rows <- 5
cols <- 4
dept <- 3
a <- array(sample(1:100, rows*cols*dept), dim = c(rows, cols, dept))
  

возврат

 > a
, , 1

     [,1] [,2] [,3] [,4]
[1,]   68   43   85   73
[2,]   39   14   21   79
[3,]    1   82   54   37
[4,]   34   59   74   83
[5,]   87   51    7   97

, , 2

     [,1] [,2] [,3] [,4]
[1,]   44   96   72   99
[2,]   84   42   80   91
[3,]   33   38   40   75
[4,]   35   20   69    6
[5,]   70   28   25   24

, , 3

     [,1] [,2] [,3] [,4]
[1,]   32   22  100   50
[2,]   94   92   62   65
[3,]    2   90   23   11
[4,]   45   98   67   17
[5,]   18   64   49   36
  

Для каждого измерения «dept» я хочу получить сумму по строкам, сохраняя при этом исходные три измерения массива. Я пытался

 b <- apply(a, c(2,3), sum)
> b
     [,1] [,2] [,3]
[1,]  229  266  191
[2,]  249  224  366
[3,]  241  286  301
[4,]  369  295  179
  

что дает правильный результат, но уменьшает его до матрицы 4 на 3, поскольку размер строки свернут до 1 и больше не является строго необходимым. Однако для моих вычислений неудобно, когда интерпретации измерений меняются каждый раз, когда я выполняю операцию, поэтому я хочу вместо этого получить массив размером 1x4x3:

 c <- array(b, dim = c(1, 4, 3))
> c
, , 1

     [,1] [,2] [,3] [,4]
[1,]  229  249  241  369

, , 2

     [,1] [,2] [,3] [,4]
[1,]  266  224  286  295

, , 3

     [,1] [,2] [,3] [,4]
[1,]  191  366  301  179
  

Это выполняет то, что я хочу, но я думаю, что это немного громоздко, и я не уверен, как обобщить его на различные операции с любым количеством измерений. Должен быть более компактный способ выполнения этих операций. Я нашел пакет `rray», но он несовместим с R 4.0.2. Обратите внимание, что мои фактические массивы намного больше, чем в этом примере, и мне придется много раз применять эти типы операций в задаче численной оптимизации, поэтому важна эффективность вычислений.

Ответ №1:

Чтобы обобщить и сохранить вычисления в одной строке, вы могли бы сделать:

 array(apply(a, 2:3, sum), c(1, dim(a)[-1]))
# , , 1
# 
# [,1] [,2] [,3] [,4]
# [1,]  229  249  241  369
# 
# , , 2
# 
# [,1] [,2] [,3] [,4]
# [1,]  266  224  286  295
# 
# , , 3
# 
# [,1] [,2] [,3] [,4]
# [1,]  191  366  301  179
  

Или, поскольку он векторизован и, следовательно, намного быстрее, используя colSums

 array(colSums(a, dims=1), c(1, dim(a)[-1]))
# , , 1
# 
# [,1] [,2] [,3] [,4]
# [1,]  229  249  241  369
# 
# , , 2
# 
# [,1] [,2] [,3] [,4]
# [1,]  266  224  286  295
# 
# , , 3
# 
# [,1] [,2] [,3] [,4]
# [1,]  191  366  301  179
  

Тест:

 set.seed(42)
A <- array(rnorm(5e4*100*10), dim=c(5e4, 100, 10))

library(rray)
microbenchmark::microbenchmark(apply=array(apply(A, 2:3, sum), c(1, dim(A)[-1])),
                               colSums=array(colSums(A, dims=1), c(1, dim(A)[-1])),
                               rray_sum=rray_sum(A, 1))  ## rray: see other answer
# Unit: milliseconds
#     expr        min         lq       mean     median         uq        max neval cld
#    apply 1273.51152 1381.72037 1416.33429 1395.84693 1433.72407 1848.88436   100   b
#  colSums   72.07086   73.02890   73.85052   73.63013   74.38916   79.70227   100  a 
# rray_sum   71.46261   72.50294   73.27564   73.00747   73.70348   80.36409   100  a 
  

Ответ №2:

Я смог остановить версию пакета, совместимую с R4.0 rray , используя

 remotes::install_github("r-lib/rray")
  

Затем желаемый результат достигается (намного быстрее) с

 # Increasing the array size for more realistic benchmarking
rows <- 500
cols <- 100
dept <- 10

draws <- rnorm(rows*cols*dept) # Standard normal draws instead of sampling from integers

a <- rray(draws, dim = c(rows, cols, dept))
b <- rray_sum(a, 1)
  

автозапуск микробенчмарка

Тестовый код:

 bm <- microbenchmark(
  base = {
    a <- array(draws, dim = c(rows, cols, dept))
    b <- apply(a, c(2,3), sum)
    c <- array(b, dim = c(1, 4, 3))
    c
  },
  rray = {
    a <- rray(draws, dim = c(rows, cols, dept))
    b <- rray_sum(a, 1)
    b
  }, times = 100)

> bm
Unit: microseconds
 expr    min     lq     mean  median      uq     max neval
 base 8619.9 8763.9 9245.898 8832.05 8984.25 20968.5   100
 rray  838.6  939.6 1186.008 1103.50 1134.40 13580.8   100