Эффективно умножить разные матрицы в трехмерном массиве с разными векторами в R - PullRequest
Новая
       88

Эффективно умножить разные матрицы в трехмерном массиве с разными векторами в R

0 голосов
/ 18 апреля 2019

У меня есть Z, который представляет собой n=1000 x m=500,000 x p=3 массив Y, который представляет собой p=3 x n=1000 матрицу. Для каждого i в 1:n я хочу умножить матрицу: Z[i, ,] %*% Y[, i]

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

У меня есть конкретная форма для матрицы Y, как описано в mwe ниже. Меня интересуют решения, использующие эту структуру, а также более общие решения, которые этого не делают.

Память также вызывает беспокойство, так как кажется, что вам потребуется как минимум удвоенный размер объекта Z. Я устанавливаю m = 10000, чтобы облегчить жизнь, но на самом деле m = 500000. 

n=1000; m=10000; p=3
Z = array( rnorm(n*m*p) , dim = c(n, m, p))
Y = matrix(NA, nrow=p, ncol = n)
index = rbinom(n,1,0.5)==1
Y[,index] = c(0,1,2)
Y[,!index] = c(2,1,0)
pryr::object_size(Z)
pryr::object_size(Y)

# methods that do not rely on the structure of Y
method1 = do.call(cbind,matrix(mapply(plyr::alply(.data = Z, .margins = 1), plyr::alply(Y, .margins = 2), FUN = function(x,y) x %*% y , SIMPLIFY = FALSE )))
method2 = do.call(cbind,matrix(map2(.x = plyr::alply(.data = Z, .margins = 1), .y = plyr::alply(Y, .margins = 2 ), ~ .x %*% .y ) ))
method3 =  do.call(cbind,lapply(seq_len(dim(Z)[1]),function(i) Z[i,,]%*%Y[,i]))
identical(method1, method2)
identical(method2, method3)

# method that does rely on the structure of Y
fn4 = function(Z, index){
  m = dim(Z)[2]; n = dim(Z)[1]
  res = matrix(NA, m, n)
  method4.1 = apply( Z[index,,], 1, function(x) x %*% c(0,1,2) )
  method4.2 = apply( Z[!index,,], 1, function(x) x %*% c(2,1,0) )
  res[,index] = method4.1
  res[,!index] = method4.2
  res
}
method4 = fn4(Z, index)
identical(method4, method3)

bench::mark(
  method1 = do.call(cbind,matrix(mapply(plyr::alply(.data = Z, .margins = 1), plyr::alply(Y, .margins = 2), FUN = function(x,y) x %*% y , SIMPLIFY = FALSE ))),
  method2 = do.call(cbind,matrix(map2(.x = plyr::alply(.data = Z, .margins = 1), .y = plyr::alply(Y, .margins = 2 ), ~ .x %*% .y ) )),
  method3 = do.call(cbind,lapply(seq_len(dim(Z)[1]),function(i) Z[i,,]%*%Y[,i])),
  method4 = fn4(Z, index),
  iterations = 4
)
...