首页 文章

通过向量乘以矩阵行?

提问于
浏览
61

我正在优化一个函数,我想摆脱缓慢的循环 . 我正在寻找一种更快的方法来将矩阵的每一行乘以一个向量 .

有任何想法吗?

编辑:

我不是在寻找'经典'乘法 .

例如 . 我有一个有23列和25行的矩阵和一个长度为23的向量 . 结果我想要矩阵25x23,每行乘以向量 .

4 回答

  • 66

    我想你正在寻找 sweep() .

    > (mat <- matrix(rep(1:3,each=5),nrow=3,ncol=5,byrow=TRUE))
         [,1] [,2] [,3] [,4] [,5]
    [1,]    1    1    1    1    1
    [2,]    2    2    2    2    2
    [3,]    3    3    3    3    3
    > vec <- 1:5
    > sweep(mat,MARGIN=2,vec,`*`)
         [,1] [,2] [,3] [,4] [,5]
    [1,]    1    2    3    4    5
    [2,]    2    4    6    8   10
    [3,]    3    6    9   12   15
    

    它是R的核心功能之一,虽然多年来已经对其进行了改进 .

  • 35
    > MyMatrix <- matrix(c(1,2,3, 11,12,13), nrow = 2, ncol=3, byrow=TRUE)
    > MyMatrix
         [,1] [,2] [,3]
    [1,]    1    2    3
    [2,]   11   12   13
    > MyVector <- c(1:3)
    > MyVector
    [1] 1 2 3
    

    你可以使用:

    > t(t(MyMatrix) * MyVector)
         [,1] [,2] [,3]
    [1,]    1    4    9
    [2,]   11   24   39
    

    要么:

    > MyMatrix %*% diag(MyVector)
         [,1] [,2] [,3]
    [1,]    1    4    9
    [2,]   11   24   39
    
  • 23

    实际上, sweep 不是我电脑上最快的选择:

    MyMatrix <- matrix(c(1:1e6), ncol=1e4, byrow=TRUE)
    MyVector <- c(1:1e4)
    
    Rprof(tmp <- tempfile(),interval = 0.001)
    t(t(MyMatrix) * MyVector) # first option
    Rprof()
    MyTimerTranspose=summaryRprof(tmp)$sampling.time
    unlink(tmp)
    
    Rprof(tmp <- tempfile(),interval = 0.001)
    MyMatrix %*% diag(MyVector) # second option
    Rprof()
    MyTimerDiag=summaryRprof(tmp)$sampling.time
    unlink(tmp)
    
    Rprof(tmp <- tempfile(),interval = 0.001)
    sweep(MyMatrix ,MARGIN=2,MyVector,`*`)  # third option
    Rprof()
    MyTimerSweep=summaryRprof(tmp)$sampling.time
    unlink(tmp)
    
    Rprof(tmp <- tempfile(),interval = 0.001)
    t(t(MyMatrix) * MyVector) # first option again, to check order 
    Rprof()
    MyTimerTransposeAgain=summaryRprof(tmp)$sampling.time
    unlink(tmp)
    
    MyTimerTranspose
    MyTimerDiag
    MyTimerSweep
    MyTimerTransposeAgain
    

    这会产生:

    > MyTimerTranspose
    [1] 0.04
    > MyTimerDiag
    [1] 40.722
    > MyTimerSweep
    [1] 33.774
    > MyTimerTransposeAgain
    [1] 0.043
    

    除了是最慢的选项之外,第二个选项达到内存限制(2046 MB) . 但是,考虑到剩下的选项,在我看来,双转置似乎比_966774好多了 .


    Edit

    只是重复尝试较小的数据:

    MyMatrix <- matrix(c(1:1e3), ncol=1e1, byrow=TRUE)
    MyVector <- c(1:1e1)
    n=100000
    
    [...]
    
    for(i in 1:n){
    # your option
    }
    
    [...]
    
    > MyTimerTranspose
    [1] 5.383
    > MyTimerDiag
    [1] 6.404
    > MyTimerSweep
    [1] 12.843
    > MyTimerTransposeAgain
    [1] 5.428
    
  • 1

    对于速度,可以在乘法之前从矢量创建矩阵

    mat <-  matrix(rnorm(1e6), ncol=1e4)
    vec <- c(1:1e4)
    mat * matrix(vec, dim(mat)[1], length(vec))
    
    library(microbenchmark)
    microbenchmark(
      transpose = t(t(mat) * vec), 
      make_matrix = mat * matrix(vec, dim(mat)[1], length(vec), byrow = TRUE),
      sweep = sweep(mat,MARGIN=2,vec,`*`))
    #Unit: milliseconds
    #       expr      min        lq     mean    median       uq      max neval cld
    #  transpose 9.940555 10.480306 14.39822 11.210735 16.19555 77.67995   100   b
    #make_matrix 5.556848  6.053933  9.48699  6.662592 10.74121 74.14429   100   a 
    #      sweep 8.033019  8.500464 13.45724 12.331015 14.14869 77.00371   100   b
    

相关问题