将矩阵乘以向量的最快方法

Question

我有一个矩阵mat 和一个向量v。我想将矩阵mat 的第一列乘以向量v 的第一个元素，并将矩阵mat 的第二列乘以向量v 的第二个元素。如图所示，我可以做到。既然我们得到了一个大矩阵，我怎样才能在 R 中更快地做到这一点？

    mat = matrix(rnorm(1500000), ncol= 100)
    v= rnorm(100)
    > system.time( mat %*% diag(v))
      user  system elapsed 
      0.02    0.00    0.02 

Answer 1

回收可以使其更快，但您在列内回收，而不是跨列回收，因此只需转置和转回即可。

t( t(mat) * v )

这应该比sweep 或%*% 更快。

microbenchmark(mat %*% diag(v),sweep(mat, 2, v, FUN = "*"), t(t(mat)*v))
Unit: milliseconds
            expr       min        lq    median        uq      max neval
             %*% 150.47301 152.16306 153.17379 161.75416 281.3315   100
           sweep  35.94029  42.67210  45.53666  48.07468 168.3728   100
   t(t(mat) * v)  16.50813  23.41549  26.31602  29.44008 160.1651   100

Answer 2

游戏有点晚了，但有人说最快吗？！这可能是Rcpp 的另一个好用处。默认情况下，此函数（称为mmult）将矩阵的每一列乘以向量的每个连续元素，但可以通过设置byrow = FALSE 选择按列执行此操作。它还检查m 和v 是否在给定byrow 选项的情况下具有适当的大小。无论如何，它快（比最好的原生 R 答案快大约 10-12 倍）...

编辑

@chris 向我提出的另一个问题提供了this great answer，试图让它与RcppArmadillo 一起工作。然而，我在这里发布的 Rcpp-only 函数似乎仍然比这快 8 倍左右，比 OP 方法快 70 倍左右。单击@chris 函数的代码链接 - 非常简单。

我会将基准测试放在首位..

require( microbenchmark )
m <- microbenchmark( mat %*% diag(v) , mmult( mat , v ) , sweep(mat, 2, v, FUN = "*") , chris( mat , v ) , t( t(mat) * v ) , times = 100L )
print( m , "relative" , order = "median" , digits = 3 )
Unit: relative
                        expr   min    lq median    uq   max neval
               mmult(mat, v)  1.00  1.00   1.00  1.00  1.00   100
               chris(mat, v) 10.74  9.31   8.15  7.27 10.44   100
               t(t(mat) * v)  9.65  8.75   8.30 15.33  9.52   100
 sweep(mat, 2, v, FUN = "*") 20.51 18.35  22.18 21.39 16.94   100
             mat %*% diag(v) 80.44 70.11  73.12 70.68 54.96   100

继续浏览以了解 mmult 的工作原理并返回与 OP 相同的结果...

require( Rcpp )

#  Source code for our function
func <- 'NumericMatrix mmult( NumericMatrix m , NumericVector v , bool byrow = true ){
  if( byrow );
    if( ! m.nrow() == v.size() ) stop("Non-conformable arrays") ;
  if( ! byrow );
    if( ! m.ncol() == v.size() ) stop("Non-conformable arrays") ;

  NumericMatrix out(m) ;

  if( byrow ){
    for (int j = 0; j < m.ncol(); j++) {
      for (int i = 0; i < m.nrow(); i++) {
        out(i,j) = m(i,j) * v[j];
      }
    }
  }
  if( ! byrow ){
    for (int i = 0; i < m.nrow(); i++) {
      for (int j = 0; j < m.ncol(); j++) {
        out(i,j) = m(i,j) * v[i];
      }
    }
  }
  return out ;
}'

#  Make it available
cppFunction( func )

#  Use it
res1 <- mmult( m , v )

#  OP function
res2 <- mat %*% diag(v)

#  Same result?
identical( res1 , res2 ) # Yes!!
[1] TRUE

Answer 3

sweep 在我的机器上运行起来似乎快了一点

sweep(mat, 2, v, FUN = "*")

一些基准测试：

> microbenchmark(mat %*% diag(v),sweep(mat, 2, v, FUN = "*"))

Unit: milliseconds
  expr       min        lq   median        uq      max neval
   %*% 214.66700 226.95551 231.2366 255.78493 349.1911   100
 sweep  42.42987  44.72254  62.9990  70.87403 127.2869   100