On Wed, 11 Mar 2009, Camarda, Carlo Giovanni wrote:
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?
This is fairly quick:
rowSums( aperm(A * as.vector( M[ rep(1:ncol(A),each=nrow(A)),]),c(1,3,2)), dims
= 2 )
But my advice is to code this in C along the lines of your first solution
(using the BLAS routines to carry it out in the inner products). Your code
will be easier to read and debug and will probably be faster and easier on
memory, too.
Years ago I wrote a lot of stuff like this in native S. I 'optimized' the
heck out of it using tricks like the one above as I was debugging the
code. I had to rewrite the bulk of it in C anyway. The S code was so hard
to read that I could not migrate it to C bit by bit. I had to start over
and the effort spent debugging it in S was lost.
As an alternative you might try the 'jit' package on your first solution.
HTH,
Chuck
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}}
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Charles C. Berry (858) 534-2098
Dept of Family/Preventive Medicine
E mailto:cbe...@tajo.ucsd.edu UC San Diego
http://famprevmed.ucsd.edu/faculty/cberry/ La Jolla, San Diego 92093-0901
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
