[R] speeding up a special product of three arrays

[Vincent Goulet]

Le ven. 09 mai à 03:44, Dimitris Rizopoulos a écrit :

> try this:
>
> A <- matrix(rnorm(10*4), 10, 4)
> B <- matrix(rnorm(3*4), 3, 4)
> C <- matrix(rnorm(5*4), 5, 4)
>
> nrA <- nrow(A); nrB <- nrow(B); nrC <- nrow(C)
> ind <- as.matrix(expand.grid(1:nrA, 1:nrB, 1:nrC))
> D <- rowSums(A[ind[, 1], ] * B[ind[, 2], ] * C[ind[, 3], ])
> dim(D) <- c(nrA, nrB, nrC)
> D

Dimitris,

We basically have the same solution. Actually, I first went exactly  
the same path as you, but I didn't think of expand.grid(), so I was  
creating the index matrix manually (easy, but ugly).

Now, just out of curiosity, I made a few tests with ~ 100-row matrices  
and my solution is slightly faster (to the order of 1 sec. vs 1.2  
sec). :-) This is mostly due to slice.index() being faster than  
expand.grid(). In fact, if one doesn't mind the ugliness, building the  
indexes manually brings your solutions on par with mine in terms of  
speed. Doing this is left as an exercise to the reader...

Does this matter? No. It's just that it's a sunny Friday...

Cheers

Vincent

>
>
>
> I hope it helps.
>
> Best,
> Dimitris
>
> ----
> Dimitris Rizopoulos
> Biostatistical Centre
> School of Public Health
> Catholic University of Leuven
>
> Address: Kapucijnenvoer 35, Leuven, Belgium
> Tel: +32/(0)16/336899
> Fax: +32/(0)16/337015
> Web: http://med.kuleuven.be/biostat/
> http://www.student.kuleuven.be/~m0390867/dimitris.htm
>
>
> ----- Original Message ----- From: "Giuseppe Paleologo" <paleologo at gmail.com 
> >
> To: <r-help at r-project.org>
> Sent: Thursday, May 08, 2008 11:20 PM
> Subject: [R] speeding up a special product of three arrays
>
>
>> I am struggling with R code optimization, a recurrent topic on this  
>> list.
>>
>> I have three arrays, say A, B and C, all having the same number of  
>> columns.
>> I need to compute an array D whose generic element is
>>
>> D[i, j, k] <- sum_n A[i, n]*B[j, n]*C[k, n]
>>
>> Cycling over the three indices and subsetting the columns won't do.  
>> Is there
>> any way to implement this efficiently in R or should I resign to do  
>> this in
>> C?
>>
>> Thanks,
>>
>> Giuseppe
>>
>> [[alternative HTML version deleted]]
>>
>> ______________________________________________
>> 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.
>
>
> Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm
>
> ______________________________________________
> 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.