0

I'm trying to find a much more efficient way to code in R the following matrix: Let A and C be two 3D array of dimension (n, n, m) and B a matrix of dimension (m, m), then M is an (n, n) matrix such that:

M_ij = SUM_kl A_ijk * B_kl * C_ijl

for (i in seq(n)) {
  for (j in seq(n)) {
    M[i, j] <- A[i,j,] %*% B %*% C[i,j,]
    }
  }

It is possible to write this with the TensorA package using i and j as parallel dimension, but I'd rather stay with base R object.

einstein.tensor(A %e% log(B), C, by = c("i", "j"))

Thanks!

Theo
  • 57,719
  • 8
  • 24
  • 41
saintclair
  • 1
  • isn't A[I,j,] a vector? can't you write it as a single loop by omitting I or j and then doing matrix multiplication? and if you want to "drop" the loop there should be a way with apply – MichaelChirico Aug 04 '19 at 10:33

1 Answers1

0

I don't know if this would be faster, but it would avoid one level of looping:

for (i in seq(n))
  M[i,] <- diag(A[i,,] %*% B %*% t(C[i,,]))

It gives the same answer as yours in this example:

n <- 2
m <- 3
A <- array(1:(n^2*m), c(n, n, m))
C <- A + 1
B <- matrix(1:(m^2), m, m)

M <- matrix(NA, n, n)
for (i in seq(n))
  M[i,] <- diag(A[i,,] %*% B %*% t(C[i,,]))
M
#      [,1] [,2]
# [1,] 1854 3216
# [2,] 2490 4032

Edited to add: Based on https://stackoverflow.com/a/42569902/2554330, here's a slightly faster version:

for (i in seq(n))
  M[i,] <- rowSums((A[i,,] %*% B) * C[i,,])

I did some timing with n <- 200 and m <- 300, and this was the fastest at 3.1 sec, versus my original solution at 4.7 sec, and the one in the question at 17.4 sec.

user2554330
  • 37,248
  • 4
  • 43
  • 90