Search code examples
rrcppcross-productrcpparmadillo

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

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?

Solution

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