Есть ли эффективный способ извлечь фрагмент трехмерного массива?

Подумайте об этом:

a = array(1:60, c(3,4,5))

# Extract (*,2,2)^th elements
a[cbind(1:3,2,2)]
                                
# ^^ returns a vector of length 3, c(16, 17, 18)
# Note this is asymptotically inefficient -- cbind(1:k,2,2) is a kx3 matrix
# So it's not making use of slicing
 
# Extract (1, y_i, z_i)^th elements
N = 100000
y = sample(1:4, N, replace = TRUE)
z = sample(1:5, N, replace = TRUE)
a[cbind(1, y, z)]

# ^^ returns a vector of length N

Как мне эффективно извлечь элементы (*, y_i, z_i), получив результат в виде матрицы Nx3 (т. е. двумерного массива)?

Обратите внимание этот вопрос аналогичен, но единственный ответ выполняется поэлементно + поэтому будет медленно.

r arrays performance matrix indexing

19.08.2024 11:29

Стоит ли изучать PHP в 2023-2024 годах?

Привет всем, сегодня я хочу высказать свои соображения по поводу вопроса, который я уже много раз получал в своем сообществе: "Стоит ли изучать PHP в...

Поведение ключевого слова "this" в стрелочной функции в сравнении с нормальной функцией

В JavaScript одним из самых запутанных понятий является поведение ключевого слова "this" в стрелочной и обычной функциях.

Приемы CSS-макетирования - floats и Flexbox

Здравствуйте, друзья-студенты! Готовы совершенствовать свои навыки веб-дизайна? Сегодня в нашем путешествии мы рассмотрим приемы CSS-верстки - в...

Тестирование функциональных ngrx-эффектов в Angular 16 с помощью Jest

В системе управления состояниями ngrx, совместимой с Angular 16, появились функциональные эффекты. Это здорово и делает код определенно легче для...

Концепция локализации и ее применение в приложениях React ⚡️

Локализация - это процесс адаптации приложения к различным языкам и культурным требованиям. Это позволяет пользователям получить опыт, соответствующий...

Пользовательский скаляр GraphQL

Листовые узлы системы типов GraphQL называются скалярами. Достигнув скалярного типа, невозможно спуститься дальше по иерархии типов. Скалярный тип...

Перейти к ответу Данный вопрос помечен как решенный

Ответы 4

Вы можете повторить каждый элемент 1-го вектора индексации N раз и связать его со 2-м и 3-м векторами индексации. Затем измените форму извлеченного фрагмента до нужного размера.

`dim<-`(a[cbind(rep(1:3, each = N), y, z)], c(N, 3))

#           [,1] [,2] [,3]
#      [1,]   40   41   42
#      [2,]   49   50   51
#      [3,]   34   35   36
#      [4,]    7    8    9
#      [5,]   46   47   48
#      [6,]   16   17   18
#      [7,]   22   23   24
#      [8,]   25   26   27
#      [9,]   49   50   51
#     [10,]   16   17   18
#  [ reached getOption("max.print") -- omitted 99990 rows ]

Вы также можете использовать sapply

sapply(1:3, \(i) a[cbind(i, y, z)])

Контрольный показатель

Я увеличиваю размер a до 30x40x50 и N=1000. Для меня удивительно, что sapply быстрее, чем матричное индексирование.

a = array(1:60000, c(30, 40, 50))
N = 1000
y = sample(1:40, N, replace = TRUE)
z = sample(1:50, N, replace = TRUE)

microbenchmark::microbenchmark(
  matrix_indexing = `dim<-`(a[cbind(rep(seq_len(nrow(a)), each = N), y, z)], c(N, nrow(a))),
  sapply_indexing = sapply(seq_len(nrow(a)), \(i) a[cbind(i, y, z)]),
  times = 100L,
  check = "identical",
  unit = "relative"
)

# Unit: relative
#             expr      min       lq     mean   median       uq       max neval cld
#  matrix_indexing 1.937193 1.930684 1.621237 1.908087 1.927041 0.1649859   100  a 
#  sapply_indexing 1.000000 1.000000 1.000000 1.000000 1.000000 1.0000000   100   b

Если я уменьшу N до 100, результат изменится.

N <- 100

# Unit: relative
#             expr      min       lq    mean   median       uq      max neval cld
#  matrix_indexing 1.000000 1.000000 1.00000 1.000000 1.000000 1.000000   100  a 
#  sapply_indexing 1.515337 1.519891 1.54679 1.516215 1.517535 1.756281   100   b

из-за очень маленького i?

— 19.08.2024 13:01

удивительно, что sapply здесь так великолепно выступает, сегодняшняя звезда!

— 19.08.2024 13:19

кажется скорость зависит от распределения габаритов, интересно!

— 19.08.2024 13:21

The values in ‘data’ are taken to be those in the array with the leftmost subscript moving fastest <?array> может быть здесь играет роль.

— 19.08.2024 14:36

19.08.2024 12:41

Для удобного кодирования, возможно, вас заинтересует asplit

do.call(rbind, asplit(a, -1)[cbind(y, z)])

но это определенно медленнее, чем `dim<-` от @Darren Tsai

Контрольный показатель

microbenchmark(
  darren = `dim<-`(a[cbind(rep(1:3, each = N), y, z)], c(N, 3)),
  tic = do.call(rbind, asplit(a, -1)[cbind(y, z)]),
  check = "identical",
  unit = "relative"
)

шоу

Unit: relative
   expr      min       lq     mean   median       uq      max neval
 darren 1.000000  1.00000  1.00000  1.00000  1.00000  1.00000   100
    tic 9.000844 11.42021 16.63028 13.42301 19.16755 25.68256   100

19.08.2024 12:58

Ответ принят как подходящий

Можно было бы вычислить линейный индекс массива, используя его размер.

outer((y-1)*dim(a)[1] + (z-1)*prod(dim(a)[1:2]), seq_len(nrow(a)), \(x, y) a[x+y])

do.call(cbind, lapply(seq_len(nrow(a)), \(x, y) a[x+y], (y-1)*dim(a)[1] + (z-1)*prod(dim(a)[1:2])))

Или уменьшите яркость и выберите целые столбцы или целые строки.

t(`dim<-`(a, c(dim(a)[1], prod(dim(a)[-1])))[,y + (z - 1)*dim(a)[2]])

t(`dim<-`(a, c(dim(a)[1], prod(dim(a)[-1]))))[y + (z - 1)*dim(a)[2],]

Тест:

m <- alist(
  matrix_indexing = `dim<-`(a[cbind(rep(seq_len(nrow(a)), each = N), y, z)], c(N, nrow(a))),
  sapply_indexing = sapply(seq_len(nrow(a)), \(i) a[cbind(i, y, z)]),
  linearIndex = outer((y-1)*dim(a)[1] + (z-1)*prod(dim(a)[1:2]), seq_len(nrow(a)), \(x, y) a[x+y]),
  linearIndex2 = do.call(cbind, lapply(seq_len(nrow(a)), \(x, y) a[x+y], (y-1)*dim(a)[1] + (z-1)*prod(dim(a)[1:2]))),
  asplit = do.call(rbind, asplit(a, -1)[cbind(y, z)]),
  asplit1 = do.call(rbind, asplit(a, -1)[cbind(y, z)]),
  aperm = `dim<-`(aperm(a, c(2, 3, 1)), c(prod(dim(a)[-1]), dim(a)[1]))[y + (z - 1)*dim(a)[2],],
  dimT = t(`dim<-`(a, c(dim(a)[1], prod(dim(a)[-1])))[,y + (z - 1)*dim(a)[2]]),
  tDim = t(`dim<-`(a, c(dim(a)[1], prod(dim(a)[-1]))))[y + (z - 1)*dim(a)[2],]
)

a = array(1:60, c(3,4,5))
N = 100000
y = sample(1:4, N, replace = TRUE)
z = sample(1:5, N, replace = TRUE)
bench::mark(exprs = m)
#  expression           min   median `itr/sec` mem_alloc `gc/sec` n_itr  n_gc
#  <bch:expr>      <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl> <int> <dbl>
#1 matrix_indexing   9.63ms  10.62ms     90.4     6.87MB     21.6    46    11
#2 sapply_indexing   3.73ms   4.05ms    216.      8.02MB     61.9   108    31
#3 linearIndex       4.75ms   5.36ms    154.      9.54MB     47.4    78    24
#4 linearIndex2      3.08ms   3.24ms    269.      7.25MB     63.7   135    32
#5 asplit           38.42ms  99.12ms     10.5     3.05MB     15.7     6     9
#6 asplit1         125.77ms 127.74ms      7.82    3.05MB     15.6     4     8
#7 aperm              1.3ms   1.33ms    622.      2.29MB     34.0   311    17
#8 dimT              2.03ms   2.09ms    431.      3.43MB     41.9   216    21
#9 tDim              1.29ms   1.33ms    683.      2.29MB     45.9   342    23

a = array(1:60000, c(30, 40, 50))
N = 1000
y = sample(1:40, N, replace = TRUE)
z = sample(1:50, N, replace = TRUE)
bench::mark(exprs = m)
#  expression           min   median `itr/sec` mem_alloc `gc/sec` n_itr  n_gc
#  <bch:expr>      <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl> <int> <dbl>
#1 matrix_indexing 951.03µs   1.01ms      960.     703KB    16.8    457     8
#2 sapply_indexing 521.71µs 554.22µs     1746.     825KB    35.3    792    16
#3 linearIndex     424.02µs 429.51µs     2255.     836KB    49.7    999    22
#4 linearIndex2    326.72µs 333.26µs     2894.     606KB    42.4   1296    19
#5 asplit             6.6ms    6.8ms      142.     637KB     9.44    60     4
#6 asplit1           6.64ms   7.56ms      128.     637KB     8.52    60     4
#7 aperm           318.52µs 338.26µs     2676.     363KB    18.9   1274     9
#8 dimT            154.78µs 162.86µs     5927.     481KB    56.0   2648    25
#9 tDim            198.04µs 210.11µs     4603.     598KB    51.8   2043    23

19.08.2024 14:27

Добавление еще пары опций в тест @GKi:

yz <- cbind(y, z)
sapply(1:dim(a)[1], \(i) a[i,,][yz])

`dim<-`(aperm(a, c(2, 3, 1)), c(prod(dim(a)[-1]), dim(a)[1]))[y + (z - 1)*dim(a)[2],]

Контрольный показатель

m <- alist(
  matrix_indexing = `dim<-`(a[cbind(rep(seq_len(nrow(a)), each = N), y, z)], c(N, nrow(a))),
  sapply_indexing = sapply(seq_len(nrow(a)), \(i) a[cbind(i, y, z)]),
  sapply_indexing2 = {yz <- cbind(y, z); sapply(1:dim(a)[1], \(i) a[i,,][yz])},
  linearIndex = outer((y-1)*dim(a)[1] + (z-1)*prod(dim(a)[1:2]), seq_len(nrow(a)), \(x, y) a[x+y]),
  linearIndex2 = do.call(cbind, lapply(seq_len(nrow(a)), \(x, y) a[x+y], (y-1)*dim(a)[1] + (z-1)*prod(dim(a)[1:2]))),
  asplit = do.call(rbind, asplit(a, -1)[cbind(y, z)]),
  aperm = `dim<-`(aperm(a, c(2, 3, 1)), c(prod(dim(a)[-1]), dim(a)[1]))[y + (z - 1)*dim(a)[2],]
)

a = array(1:60, c(3,4,5))
N = 100000
y = sample(1:4, N, replace = TRUE)
z = sample(1:5, N, replace = TRUE)
bench::mark(exprs = m)
#> # A tibble: 7 × 6
#>   expression            min   median `itr/sec` mem_alloc `gc/sec`
#>   <bch:expr>       <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl>
#> 1 matrix_indexing    3.83ms    4.7ms     210.     6.87MB     71.1
#> 2 sapply_indexing    2.32ms   3.73ms     271.    11.96MB     91.5
#> 3 sapply_indexing2    1.9ms    2.8ms     360.     5.41MB     78.6
#> 4 linearIndex        4.24ms   5.21ms     192.     9.57MB    123. 
#> 5 linearIndex2       1.89ms   3.16ms     322.     7.28MB    103. 
#> 6 asplit            26.15ms  34.32ms      31.3    3.08MB    135. 
#> 7 aperm             699.7µs   1.07ms     976.     2.29MB     62.1

Больший массив, меньше образцов:

a = array(1:60000, c(30, 40, 50))
N = 1000
y = sample(1:40, N, replace = TRUE)
z = sample(1:50, N, replace = TRUE)
bench::mark(exprs = m)
#> # A tibble: 7 × 6
#>   expression            min   median `itr/sec` mem_alloc `gc/sec`
#>   <bch:expr>       <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl>
#> 1 matrix_indexing   312.3µs  443.2µs     2223.     703KB     44.3
#> 2 sapply_indexing     305µs 495.25µs     2048.     825KB     42.6
#> 3 sapply_indexing2  493.5µs  651.8µs     1623.     729KB     30.4
#> 4 linearIndex       325.1µs  466.3µs     2246.     836KB     50.6
#> 5 linearIndex2        194µs  341.5µs     3060.     606KB     48.8
#> 6 asplit             4.57ms   4.73ms      204.     637KB     12.6
#> 7 aperm             199.7µs  277.9µs     3642.     363KB     24.8

В общем случае, не должна ли ваша версия aperm принимать трехэлементный вектор в качестве входных данных самой функции aperm()? Не то чтобы я ожидал какой-либо разницы в скорости

— 19.08.2024 17:34

Да, я считаю, что aperm() будет подходящим вариантом в общем случае (*, y_i, z_i), (y_i, *, z_i) или (y_i, z_i, *). dimT или tDim @GKi кажутся лучшими для конкретного случая (*, y_i, z_i) ОП (в зависимости от N и размеров a).