r - hadamard 将所有矩阵列与另一个矩阵相乘的最快方法

Question

这是我想要做的:我有两个矩阵 A 和 B，维度分别为 N x k1 和 N x k2。我现在想将矩阵 A 的每一列与 B 逐点相乘。

实现一在 for 循环中执行此操作。

出于速度优化的目的，我考虑过对整个操作进行矢量化 - 但事实证明，矢量化(正如我通过 kronecker 产品实现的那样)并没有改进我的运行时间以解决更大的问题。

有没有人建议如何以不同的方式实现此操作，同时考虑到运行时？

下面的代码从一个小示例开始，然后实现一个基于循环的矢量化解决方案，然后对一个更大的问题进行基准测试。

# toy example: 

N <- 5
k1 <- 2
k2 <- 3


A <- matrix(rnorm(N*k1), N, k1)
B <- matrix(rnorm(N*k2), N, k2)

colmat_prod <- function(x, y){
  
  k2 <- ncol(y)
  k1 <- ncol(x)
  res <- array(NA, c(N, k2 , k1))
  for(i in 1:k1){
    res[, ,i] <- x[,i] * y
  }
  res
}

colmat_prod_vec <- function(x, y){
  k1 <- ncol(x)
  res_vec <- c(x) * (rep(1, k1) %x% y)
  res_vec
}

colmat_prod(A, B)
colmat_prod_vec(A, B)

# > colmat_prod(A, B)
# , , 1
# 
# [,1]        [,2]        [,3]
# [1,]  1.95468879  0.55206339  0.24713400
# [2,] -0.02678564 -0.03762645 -0.03144102
# [3,]  0.30964437  0.26912771 -0.49451656
# [4,] -1.40719543  0.77245522 -0.47236888
# [5,] -1.71485558  0.98348809  0.16569915
# 
# , , 2
# 
# [,1]        [,2]        [,3]
# [1,]  1.60358991  0.45290242  0.20274409
# [2,] -0.21009808 -0.29513001 -0.24661348
# [3,]  0.04069121  0.03536681 -0.06498577
# [4,] -2.89562745  1.58950383 -0.97200734
# [5,] -1.59504293  0.91477425  0.15412217
# 
# > colmat_prod_vec(A, B)
# [,1]        [,2]        [,3]
# [1,]  1.95468879  0.55206339  0.24713400
# [2,] -0.02678564 -0.03762645 -0.03144102
# [3,]  0.30964437  0.26912771 -0.49451656
# [4,] -1.40719543  0.77245522 -0.47236888
# [5,] -1.71485558  0.98348809  0.16569915
# [6,]  1.60358991  0.45290242  0.20274409
# [7,] -0.21009808 -0.29513001 -0.24661348
# [8,]  0.04069121  0.03536681 -0.06498577
# [9,] -2.89562745  1.58950383 -0.97200734
# [10,] -1.59504293  0.91477425  0.15412217


# speed: 
N <- 10000
k1 <- 1000
k2 <- 9
A1 <- matrix(rnorm(N*k1), N, k1)
B1 <- matrix(rnorm(N*k2), N, k2)

library(microbenchmark)
microbenchmark(colmat_prod(A1, B1), 
               colmat_prod_vec(A1, B1), 
               times = 10)
#Unit: seconds
#expr      min       lq      mean   median uq      max neval
#colmat_prod(A1, B1) 1.981737 2.179122  2.769812  2.32343 2.680407  4.96276    10
#colmat_prod_vec(A1, B1) 9.779629 9.955576 10.291264 10.21356 10.380702 11.70494    10

Answer 1

您可以尝试 apply(A, 2, '*', B) 并像 colmat_prod 一样使用 array(apply(A, 2、'*', B), c(dim(B), ncol(A))):

identical(array(apply(A, 2, '*', B), c(dim(B), ncol(A))), colmat_prod(A, B))
#[1] TRUE

另一种选择是对 A 的列使用 rep:

array(A[,rep(seq_len(ncol(A)), each=ncol(B))] * as.vector(B), c(dim(B), ncol(A)))

时间:

library(microbenchmark)
microbenchmark(colmat_prod(A1, B1), 
               colmat_prod_vec(A1, B1),
               array(apply(A1, 2, '*', B1), c(dim(B1), ncol(A1))),
               array(A1[,rep(seq_len(ncol(A1)), each=ncol(B1))] * as.vector(B1), c(dim(B1), ncol(A1))),
               times = 10)
#Unit: milliseconds
#                                                                                            expr      min        lq      mean    median        uq       max neval  cld
#                                                                             colmat_prod(A1, B1) 831.5437  857.0305  910.5694  878.6842  999.5354 1025.0915    10   c 
#                                                                         colmat_prod_vec(A1, B1) 981.9241 1010.9482 1174.1700 1162.7004 1319.3478 1444.6158    10    d
#                                              array(apply(A1, 2, "*", B1), c(dim(B1), ncol(A1))) 716.1469  725.7862  765.4987  732.2520  789.3843  907.4417    10  b  
# array(A1[, rep(seq_len(ncol(A1)), each = ncol(B1))] * as.vector(B1),      c(dim(B1), ncol(A1))) 404.8460  406.2848  430.4043  428.2685  458.9400  462.0634    10 a