Быстрый метод R для матричного продукта с конечным продуктом вместо суммы

Question

В R можно выполнить перекрестное произведение, используя %*% между двумя матрицами M1: n x p и M2: p x d, которые имеют общую длину одного измерения.

Чтобы получить перекрестное произведение, нужно умножить для каждой строки 1..n в M1 и столбце 1..d в M2 относительное p_M1 x p_M2, а затем суммировать полученный вектор.

Но вместо суммы я бы хотел получить произведение prod(p_M1 x p_M2).

Я могу сделать это с помощью вложенных циклов в R, но это очень медленно и мои матрицы очень большие. Есть ли альтернатива так же быстро, как %*%?

Пример:

    set.seed(1)
    a <- matrix(sample((1:100) / 100, 15), ncol = 3)
    b <- matrix(sample((1:100) / 100, 15), ncol = 5)

    # This produces the usual cross-product...
    a %*% b

    # ...which can be done also using loops
    do.call('cbind', lapply(1:5, function(i) {
        sapply(1:5, function(j) {
            sum(a[i,] * b[,j])
        })
    }))

    # But I need to do the product of the paired vectors instead of the sum. I could use a nested loop but it takes hours.
    do.call('cbind', lapply(1:5, function(i) {
        sapply(1:5, function(j) {
            prod(a[i,] * b[,j])
        })
    }))

Answer 1

После моего комментария, вот метод с пакетом matrixStats и outer для выполнения расчета.

# nested loop
mat1 <- 
    do.call('cbind', lapply(1:5, function(i) {
        sapply(1:5, function(j) {
            prod(a[i,] * b[,j])
        })
    }))

# vectorized-ish
library(matrixStats)

mat2 <- outer(colProds(b), rowProds(a))

Теперь проверьте, что они численно эквивалентны.

all.equal(mat1, mat2)
[1] TRUE

Если вы хотите внешний вид %*%, вы можете изменить это на

mat2 <- colProds(b) %o% rowProds(a)

---

Вы можете придерживаться базы R, если хотите избежать посылок. Вот один из методов.

mat3 <- outer(
               vapply(seq_len(ncol(b)), function(x) prod(b[, x]), numeric(1L)),
               vapply(seq_len(nrow(a)), function(x) prod(a[x, ]), numeric(1L))
              ))

---

проверяя скорость этих двух, я получаю следующее

library(microbenchmark)

microbenchmark(nest=
                do.call('cbind', lapply(1:5, function(i) {
                        sapply(1:5, function(j) {
                                prod(a[i,] * b[,j])
                               })
                        })),
               vect=outer(colProds(b), rowProds(a)),
               baseVect=outer(
                  vapply(seq_len(ncol(b)), function(x) prod(b[, x]), numeric(1L)),
                  vapply(seq_len(nrow(a)), function(x) prod(a[x, ]), numeric(1L))
               ))

Unit: microseconds
  expr     min       lq      mean  median       uq      max neval
  nest    129.228 133.2225 172.43874 136.833 142.9640 3531.144   100
  vect     23.831  25.8690  28.38306  27.705  29.1815   94.546   100
 baseVect  27.223  29.8970  57.85946  31.471  32.8400 2647.373   100