R help - May 2008 - 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]]

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 at act.ulaval.ca   http://vgoulet.act.ulaval.ca

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: +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