speeding up a special product of three arrays

[Giuseppe Paleologo]

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@... 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.

[Vincent Goulet]

Le jeu. 8 mai à 17:20, Giuseppe Paleologo a écrit :

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

Hum, interesting one. Here's one solution relying on indexing. Will  
this be fast enough for your purposes?

 > set.seed(1)
 > (A <- matrix(sample(1:10, 4), 2, 2))
      [,1] [,2]
[1,]    3    5
[2,]    4    7
 > (B <- matrix(sample(1:10, 6), 3, 2))
      [,1] [,2]
[1,]    3    5
[2,]    9    4
[3,]    8    1
 > (C <- matrix(sample(1:10, 8), 4, 2))
      [,1] [,2]
[1,]    3    5
[2,]    2    7
[3,]    6    8
[4,]   10    4
 > nrA <- nrow(A)
 > nrB <- nrow(B)
 > nrC <- nrow(C)
 > res <- structure(numeric(prod(nrA, nrB, nrC)), dim = c(nrA, nrB,  
nrC))
 > res[] <- rowSums(A[slice.index(res, 1), ] * B[slice.index(res,  
2), ] * C[slice.index(res, 3), ])
 > res
, , 1

      [,1] [,2] [,3]
[1,]  152  181   97
[2,]  211  248  131

, , 2

      [,1] [,2] [,3]
[1,]  193  194   83
[2,]  269  268  113

, , 3

      [,1] [,2] [,3]
[1,]  254  322  184
[2,]  352  440  248

, , 4

      [,1] [,2] [,3]
[1,]  190  350  260
[2,]  260  472  348

HTH

---
   Vincent Goulet, Associate Professor
   École d'actuariat
   Université Laval, Québec
   Vincent.Goulet@...   http://vgoulet.act.ulaval.ca

______________________________________________
R-help@... 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.

[Dimitris Rizopoulos]

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


I hope it helps.

Best,
Dimitris

----
Dimitris Rizopoulos
Biostatistical Centre
School of Public Health
Catholic University of Leuven

Address: Kapucijnenvoer 35, Leuven, Belgium
Tel: [phone]
Fax: [phone]
Web: http://med.kuleuven.be/biostat/
     http://www.student.kuleuven.be/~m0390867/dimitris.htm


----- Original Message -----
From: "Giuseppe Paleologo" <paleologo@...>
To: <r-help@...>
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@... 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@... 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.

[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: [phone]
> Fax: [phone]
> Web: http://med.kuleuven.be/biostat/
> http://www.student.kuleuven.be/~m0390867/dimitris.htm
>
>
> ----- Original Message ----- From: "Giuseppe Paleologo" <paleologo@...
> >
> To: <r-help@...>
> 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@... 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@... 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.

______________________________________________
R-help@... 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.