# [R] matrix multiplication, tensor product, block-diagonal and fast computation

Charles C. Berry cberry at tajo.ucsd.edu
Thu Mar 12 16:24:59 CET 2009

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)
> 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