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

[Camarda, Carlo Giovanni]

Dear Chuck,

	thanks a lot for your reply.
Nevertheless, and though not so elegant, the for-loop solution seems faster (see below).

Best,
Carlo Giovanni Camarda

m <- 50
n <- 40
nl <- 1000
A <- array(1:(m*m*n), dim=c(m,m,n))
M <- matrix(1:(m*n), nrow=m, ncol=n)

# for-loop solution
system.time(
for(j in 1:nl){
    M1 <- matrix(nrow=m, ncol=n)
    for(i in 1:n){
        M1[,i] <- A[,,i] %*% M[,i]
    }
})
#   user  system elapsed 
#   5.65    0.03    5.72

# rowSums solution
system.time(
for(j in 1:nl){
    M5 <- rowSums( aperm(A * as.vector( M[ rep(1:ncol(A),each=nrow(A)),]),c(1,3,2)), dims = 2 )
})
#   user  system elapsed 
#  12.00    1.08   13.13



-----------------------------------------------------
Dr. Carlo Giovanni Camarda
Research Scientist
Max Planck Institute for Demographic Research
Laboratory of Statistical Demography
Konrad-Zuse-Straße 1
18057 Rostock - Germany
Phone: +49 (0)381 2081 172
Fax: +49 (0)381 2081 472
camarda at demogr.mpg.de
http://www.demogr.mpg.de/en/staff/camarda/default.htm
-----------------------------------------------------





-----Original Message-----
From: Charles C. Berry [mailto:cberry at tajo.ucsd.edu] 
Sent: Thursday, March 12, 2009 4:25 PM
To: Camarda, Carlo Giovanni
Cc: r-help at stat.math.ethz.ch
Subject: Re: [R] matrix multiplication, tensor product, block-diagonal and fast computation

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



----------
This mail has been sent through the MPI for Demographic Research.  Should you receive a mail that is apparently from a MPI user without this text displayed, then the address has most likely been faked. If you are uncertain about the validity of this message, please check the mail header or ask your system administrator for assistance.