Here's what I want to do: I have two matrices A and B of dimensions N x k1 and N x k2. I now want to pointwise multiply each column of the matrix A with B.
Implementation one does this in a for loop.
For speed optimization purposes, I considered to vectorize the entire operation - but it turns out vectorization (as I have implemented it, via kronecker products) did not improve my runtime for larger problems.
Does anyone have a suggestion how to differently implement this operation, having runtime in mind?
The code below starts with a small example, then implements a loop-based and vectorized solution, then benchmarks on a larger problem.
# 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
