r - 用最终产品代替总和的矩阵乘积的快速 R 方法

Question

在 R 中，可以使用 %*% 执行叉积。两个矩阵之间 M1: n x p和 M2: p x d ，即具有共同的一维长度。

对每一行做一个乘积的交叉乘积1..n在 M1 和 M2 中的第 1..d 列中，相对 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