[R] speeding up a special product of three arrays

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