tensor product

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

The issue names no repository files, tests, or entry points. Start by locating the Stan tensor_product pattern and the relevant Eigen or automatic-differentiation operations; done would require an agreed implementation with analytic derivatives plus evidence that memory use and performance improve.