matrix multiplication, tensor product, block-diagonal and fast computation

Dear R-users,

I am searching to the "best" way to compute a series of n matrix
multiplications between each matrix (mXm) in an array (mXmXn), and each
column of a matrix (mXn).

Please find below an example with four possible solutions.
The first is a simple for-loop which one might avoid; the second
solution employs the tensor product but then manually selects the right
outcomes. The third solution uses too many (and too time-consuming)
functions in order to profit of mapply. The last solution creates a
block-diagonal from the array and multiply such matrix by the vectorized
version of the first matrix. Here, often, the block-diagonal matrix may
be too large and a specific list is needed (at least AFAIK).

Does anyone have a further and possibly more effective way of computing
such operation?

Thanks in advance for any suggestion,
Carlo Giovanni Camarda


library(tensor)
library(Matrix)
A <- array(seq(0,1,length=48), dim=c(4,4,3))
M <- matrix(1:12, nrow=4, ncol=3)

# first solution (avoid the for-loop)
M1 <- matrix(nrow=4, ncol=3)
for(i in 1:3){
   M1[,i] <- A[,,i] %*% M[,i]
}
# second solution (direct picking of the right cols)
A1 <- tensor(A, M, 2, 1)
M2 <- cbind(A1[,1,1],A1[,2,2],A1[,3,3])
# third solution (avoid as.data.frame and as.matrix)
Adf0 <- apply(A, 3, as.data.frame)
Adf1 <- lapply(X=Adf0, FUN=as.matrix, nrow=4, ncol=4)
M3 <- mapply(FUN="%*%", Adf1, as.data.frame(M))
# fourth solution (often too large block-diagonal matrix)
Alist <- NULL
for(i in 1:3){ # better way to create such list for bdiag
   Alist[[i]] <- A[,,i]
}
Abd <- bdiag(Alist)
M4 <- matrix(as.matrix(Abd %*% c(M)), nrow=4, ncol=3)

----------
This mail has been sent through the MPI for Demographic ...{{dropped:10}}