I could not figure out, how to write a vectorized version of the following matrix and array multiplication:
v = rep(0.1, 2)
v2 = matrix(0, nrow = 2, ncol = 3)
A = matrix(0.25, nrow = 2, ncol = 3)
R = matrix(- 1, nrow = 2, ncol = 3)
P = array(1:12, dim = c(2, 2, 3))
for (i in 1:3) {
v2[, i] = A[, i] * (R[, i] + P[, , i] %*% v)
}
v = rowSums(v2)
Can someone help?
The problem that we are faced with is that we want to do matrix multiplication over P, a 3d array. As this is not efficiently possible, it is better to transform this array to a 2d column matrix.
P <- array(1:12, dim = c(2, 2, 3))
P <- matrix(P, nrow = 2)
The A and R matrices don't have to be adjusted for this to work
A <- matrix(0.25, nrow = 2, ncol = 3)
R <- matrix(-1, nrow = 2, ncol = 3)
Matrix multiplication with the v vector can then be accomplished by transforming the column vector into a sparse block diagonal matrix with the bdiag(...) function from the Matrix package.
library(Matrix)
v <- rep(0.1, 2)
v <- bdiag(v, v, v)
## [1,] 0.1 . .
## [2,] 0.1 . .
## [3,] . 0.1 .
## [4,] . 0.1 .
## [5,] . . 0.1
## [6,] . . 0.1
This gives the same result as before, however, this time in a vectorized manner.
v <- rowSums(A * (R + (P %*% v)))
## [1] 0.15 0.30
Edit: In the above equation I use two types of multiplication: *, the dot-product, which computes the element-wise products between two matrices of equal size, and %*%, the ordinary matrix product, which computes the matrix product of two conformable matrices, or in other words, the number of columns of the first, and the number of rows of the second matrix should be equal.
Edit2: As requested I have included how to transform the given P array to a matrix.
Edit3: This method can be generalized by introducing a variable n for the number of groups. The key adjustment is in the declaration of the block diagonal matrix, as next to individual matrices, this function also takes a list of matrices as input. This results in the following solution
library(Matrix)
n <- 3
P <- array(1:12, dim = c(2, 2, n))
P <- matrix(P, nrow = 2)
A <- matrix(0.25, nrow = 2, ncol = n)
R <- matrix(-1, nrow = 2, ncol = n)
v <- rep(0.1, 2)
v <- bdiag(rep(list(v),n))
v <- rowSums(A * (R + (P %*% v)))
## [1] 0.15 0.30