R help - Mar 2009 - 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}}

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 at 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:cberry at tajo.ucsd.edu	            UC San Diego
http://famprevmed.ucsd.edu/faculty/cberry/  La Jolla, San Diego 92093-0901