Re: [R] speeding up a special product of three arrays

From: Vincent Goulet <vincent.goulet_at_act.ulaval.ca>
Date: Fri, 09 May 2008 01:58:45 -0400

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

______________________________________________
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.
```
Received on Fri 09 May 2008 - 06:05:23 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Fri 09 May 2008 - 06:30:34 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.