tensor product

[bob-carpenter]

Summary:

From @betanalpha C_{ij} = sum{m, n = 1}^{N} A_{imn} B_{jmn}

Description:

In Stan code:

matrix tensor_product(matrix[] A, matrix[] B) {
  matrix[size(A), size(B)] C;
  for (j in 1:J)
    for (i in 1:I)
      c[i, j] = sum(A[i] .* B[j]);
  return c;
}

with calling:

matrix[M, N] A[I];
matrix[M, N] B[J];
matrix[I, J] C = tensor_product(A, B);

We could write analytic derivatives for sum(matrix .* matrix) to cut down on memory usage. I don't know if we can piggyback on any of the Eigen tensor operations to make it more efficient. Lots of memory blocking issues here in how to do that sum and elementwise product efficiently.

Current Version:

v2.15.0

[spinkney]

Apparently the derivatives are really easy, dC/dA = B and dC/dB = A.

[t4c1]

Could we implement this by just as reshaping A into a row vector and B into a vector and matrix-multiplying them (calculating outer product)?

[spinkney]

@t4c1 yes you can do that. Just to be clear, you mean that each I, J matrix is converted to a vector. So these two functions are equal (in R code)

k = 2; n=3; m = 4
A <- array(rnorm(k * m * n), c(n,m,k))
A

j = 3
B  <- array(rnorm(j * m * n), c(n,m,j))
B

tensor_product <- function(A, B) {
  I <- dim(A)[3]
  J <- dim(B)[3]
  c <- matrix(NA, I, J)
  for (j in 1:J) {
    for (i in 1:I) {
      c[i, j] <- sum(A[,,i] * B[,,j])
    }
  }
    return(c)
}

tensor_product_t4c1 <- function(A, B){
  I <- dim(A)[3]
  J <- dim(B)[3]
  c <- matrix(NA, I, J)
  for (j in 1:J) {
    for (i in 1:I) {
      c[i, j] <- matrix(c(A[,,i]), nrow = 1) %*%  matrix(c(B[,,j]), ncol = 1)
    }
  }
  return(c)
}

all.equal(tensor_product(A, B), tensor_product_t4c1(A, B) )
[1] TRUE

[mathDR]

was this never implemented?

[bob-carpenter]

Nope, this was never implemented. It shouldn't be that hard to implement, but we don't have a good tensor data structure other than an array of matrices.