使用最终乘积而不是求和的快速矩阵乘积R方法

Question

使用最终乘积而不是求和的快速矩阵乘积R方法

rperformance

4

在R中，可以使用%*%在两个矩阵M1：n x p和M2：p x d之间执行交叉乘积，其中一个维度长度是相同的。

要执行交叉乘积，需要将M1中的每一行1..n与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])
        })
    }))

- Bakaburg

1

请注意，乘法是可结合的和可交换的。有了这个想法，您可以在一次计算中获取每个矩阵感兴趣的行/列的积，然后您将使用两个向量进行操作。 matrixStats具有rowProds和colProds函数，这些函数非常快。 - lmo

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- lmo · Accepted Answer

根据我的评论，以下是使用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