将矩阵列与R中的向量元素相乘的最快方法

Question

我有一个矩阵m和矢量v。我想将第一列矩阵m乘以向量v的第一个元素，并将第二列矩阵m乘以向量v的第二个元素，依此类推。我可以使用以下代码完成它，但我正在寻找一种不需要两个转置调用的方法。我怎样才能在R中更快地做到这一点？

m <- matrix(rnorm(120000), ncol=6)
v <- c(1.5, 3.5, 4.5, 5.5, 6.5, 7.5)

system.time(t(t(m) * v))

#   user  system elapsed 
#   0.02    0.00    0.02 

Answer 1

使用一些线性代数并执行矩阵乘法，这在R中非常快。

例如

m %*% diag(v)

一些基准测试

m = matrix(rnorm(1200000), ncol=6)

 v=c(1.5, 3.5, 4.5, 5.5, 6.5, 7.5)
library(microbenchmark)
microbenchmark(m %*% diag(v), t(t(m) * v))
##    Unit: milliseconds
##            expr      min       lq   median       uq      max neval
##   m %*% diag(v) 16.57174 16.78104 16.86427 23.13121 109.9006   100
##     t(t(m) * v) 26.21470 26.59049 32.40829 35.38097 122.9351   100

Answer 2

如果列数较多，则t（t（m）* v）解决方案的性能优于矩阵乘法解决方案。但是，有一个更快的解决方案，但它在内存使用方面的成本很高。使用rep（）和乘法元素创建一个大到m的矩阵。这是比较，修改mnel的例子：

m = matrix(rnorm(1200000), ncol=600)
v = rep(c(1.5, 3.5, 4.5, 5.5, 6.5, 7.5), length = ncol(m))
library(microbenchmark)

microbenchmark(t(t(m) * v), 
  m %*% diag(v),
  m * rep(v, rep.int(nrow(m),length(v))), 
  m * rep(v, rep(nrow(m),length(v))), 
  m * rep(v, each = nrow(m)))

# Unit: milliseconds
#                                   expr        min         lq       mean     median         uq       max neval
#                            t(t(m) * v)  17.682257  18.807218  20.574513  19.239350  19.818331  62.63947   100
#                          m %*% diag(v) 415.573110 417.835574 421.226179 419.061019 420.601778 465.43276   100
#  m * rep(v, rep.int(nrow(m), ncol(m)))   2.597411   2.794915   5.947318   3.276216   3.873842  48.95579   100
#      m * rep(v, rep(nrow(m), ncol(m)))   2.601701   2.785839   3.707153   2.918994   3.855361  47.48697   100
#             m * rep(v, each = nrow(m))  21.766636  21.901935  23.791504  22.351227  23.049006  66.68491   100

正如您所看到的，在rep（）中使用“each”可以为了清晰起见而牺牲速度。 rep.int和rep之间的区别似乎可以忽略不计，两种实现都在重复运行microbenchmark（）时交换位置。请记住，ncol（m）== length（v）。

Answer 3

正如@Arun指出的那样，我不知道你会在时间效率方面打败你的解决方案。在代码可理解性方面，还有其他选择：

一种选择：

> mapply("*",as.data.frame(m),v)
      V1  V2  V3
[1,] 0.0 0.0 0.0
[2,] 1.5 0.0 0.0
[3,] 1.5 3.5 0.0
[4,] 1.5 3.5 4.5

而另一个：

sapply(1:ncol(m),function(x) m[,x] * v[x] )

Answer 4

为了完整起见，我将sweep添加到基准测试中。尽管它有一些误导性的属性名称，但我认为它可能比其他替代品更具可读性，而且速度非常快：

n = 1000
M = matrix(rnorm(2 * n * n), nrow = n)
v = rnorm(2 * n)

microbenchmark::microbenchmark(
  M * rep(v, rep.int(nrow(M), length(v))),
  sweep(M, MARGIN = 2, STATS = v, FUN = `*`),
  t(t(M) * v),
  M * rep(v, each = nrow(M)),
  M %*% diag(v)
)

Unit: milliseconds
                                       expr         min          lq        mean
    M * rep(v, rep.int(nrow(M), length(v)))    5.259957    5.535376    9.994405
 sweep(M, MARGIN = 2, STATS = v, FUN = `*`)   16.083039   17.260790   22.724433
                                t(t(M) * v)   19.547392   20.748929   29.868819
                 M * rep(v, each = nrow(M))   34.803229   37.088510   41.518962
                              M %*% diag(v) 1827.301864 1876.806506 2004.140725
      median          uq        max neval
    6.158703    7.606777   66.21271   100
   20.479928   23.830074   85.24550   100
   24.722213   29.222172   92.25538   100
   39.920664   42.659752  106.70252   100
 1986.152972 2096.172601 2432.88704   100

Answer 5

正如bluegrue所做的那样，一个简单的rep也足以执行逐元素乘法。

乘法和求和的数量大幅减少，就像执行使用diag()的简单矩阵乘法一样，在这种情况下可以避免许多零乘法。

m = matrix(rnorm(1200000), ncol=6)
v=c(1.5, 3.5, 4.5, 5.5, 6.5, 7.5)
v2 <- rep(v,each=dim(m)[1])
library(microbenchmark)
microbenchmark(m %*% diag(v), t(t(m) * v), m*v2)

Unit: milliseconds
          expr       min        lq      mean    median        uq      max neval cld
 m %*% diag(v) 11.269890 13.073995 16.424366 16.470435 17.700803 95.78635   100   b
   t(t(m) * v)  9.794000 11.226271 14.018568 12.995839 15.010730 88.90111   100   b
        m * v2  2.322188  2.559024  3.777874  3.011185  3.410848 67.26368   100  a