Efficiently computing sum of cross product for two 3D arrays in R

Question

For two 3D arrays in R, for example,

N <- 1000
x <- rnorm(N*3*3);   dim(x) <- c(N,3,3)
y <- rnorm(N*3*3);   dim(y) <- c(N,3,3)

I can do the following cross product by loop:

gg <- 0
for (n in 1:dim(x)[1]){
    gg <- gg + t(x[n,,]) %*% y[n,,]
}

My question is can we do it more efficiently (e.g., by vectorization or rcpp) for very large N instead of using loop?

there is matrix multiplication implemented in this post using Rcpp https://stackoverflow.com/questions/37191673/matrix-multiplication-in-rcpp — Ronak Shah, Mar 20 '19 at 05:24
Thank you for the comment! I've written rcpp and tested against loop, but it is 2 times slower than loop... I think F. Privé's solution below is most fast according to my tests. — user7283235, Mar 20 '19 at 06:49

score 5 · Accepted Answer · answered Mar 20 '19 at 06:15

5

If you rewrite your problem mathematically, you can show that it is equivalent to:

dim(x) <- c(3 * N, 3)
dim(y) <- c(3 * N, 3)
gg2 <- crossprod(x, y)

which should be very fast and should not make any copy.

answered Mar 20 '19 at 06:15

F. Privé

Thank you very much! It runs about 12 times faster than loop on my computer! – user7283235 Mar 20 '19 at 06:45

1 Answers1