Цикл for более эффективен, чем вы думаете
Ваша задача умножения пар n (A, B) не эквивалентна тензорному умножению в обычном смысле, хотяwhuber представил очень аккуратный способ превращения его в матричное умножение, складывая Bs как блоки в разреженной матрице.
Вы сказали, что хотите избежать цикла for, но подход for-loopна самом деле очень конкурентоспособен, когда запрограммирован эффективно, и я бы посоветовал вам пересмотреть его.
Я буду использовать те же обозначения, что и whuber, с A измерения nxk и B измерения kxkxn, например:
n <- 1e4
k <- 3
A <- array(rnorm(k*n),c(n,k))
B <- array(rnorm(k*k*n),c(k,k,n))
Простое и эффективное решение для цикла будет выглядеть следующим образом:
justAForLoop <- function(A,B) {
n <- nrow(A)
for (i in 1:n) A[i,] <- A[i,] %*% B[,,i]
A
}
, создающее матрицу результатов nxk.
Я изменил функцию f3 whuber для загрузки матрицыпакет, в противном случае функция sparseMatrix недоступна.Моя версия f3 немного быстрее, чем оригинал, потому что я удалил последнюю транспонирование матрицы перед возвратом результата.С этой модификацией он возвращает идентичные числовые результаты в justAForLoop.
f3 <- function(a, b) {
require(Matrix)
n <- dim(b)[3]
k <- dim(b)[2]
i0 <- (1:n-1)*k+1
i <- rep(i0, each=k)
j <- 1:(k*n)
aa <- sparseMatrix(i, j, x=c(t(a)), dims=c(n*k, n*k))
bb <- matrix(aperm(b, c(1,3,2)), nrow=n*k)
(aa %*% bb)[i0, ]
}
Теперь я перезапускаю симуляцию Вубера в новом сеансе R:
> k <- 3
> n <- 1e6
> a <- matrix(runif(k*n), ncol=k)
> b <- array(runif(k^2*n), dim=c(k,k,n))
>
> system.time(c1 <- f1(a,b))
user system elapsed
3.40 0.09 3.50
> system.time(c3 <- f3(a,b))
Loading required package: Matrix
user system elapsed
1.06 0.24 1.30
> system.time(c4 <- justAForLoop(a,b))
user system elapsed
1.27 0.00 1.26
Подход for-loop на самом делесамый быстрый с узким краем.Это намного быстрее, чем f1, который опирается на sapply.(Моя машина - ПК с Windows 10 с 32 ГБ ОЗУ под управлением R 3.6.0.)
Если я запускаю все три метода во второй раз, то f3 становится самым быстрым, потому что на этот раз пакет Matrix уже находится впуть поиска и не должен быть перезагружен:
> system.time(c1 <- f1(a,b))
user system elapsed
3.23 0.04 3.26
> system.time(c3 <- f3(a,b))
user system elapsed
0.33 0.20 0.53
> system.time(c4 <- justAForLoop(a,b))
user system elapsed
1.28 0.01 1.30
Однако f3 использует больше оперативной памяти, чем цикл for.На моем ПК я могу успешно запустить justAForLoop с n=1e8, тогда как f1 и f3 не хватает памяти и произойдет сбой.
Сводка
Подход с прямым циклом гораздо более эффективен, чем sapply.
Для вашей задачи с n = 10 000 умножений матриц выполнение цикла for является простым и эффективным, занимая <0,02 сек.В отличие от этого, простая загрузка пакета с разреженными матричными функциями требует около 2 / 3сек. </p>
Для n между 1-10 миллионами решение для разреженных матриц whuber начинает выигрывать, особенно если пакет Matrix уже загружен.
Цикл for использует наименьшую оперативную память из трех методов.Для n на 100 миллионов на моем ПК с 32 ГБ ОЗУ работает только подход for-loop.