Inquiry about Kronecker product

Modeling
1

Dear community,

I am trying to implement a spatial model in stan that needs the Kronecker product. My attempt is below, but I do not think is optimized. I tried with the kronecker_prod fuction given here State Space Models in Stan, but unfortunately, It looks like Stan does not allow to index a matrix, like, for example C[row_start:row_end, col_start:col_end]. Can you help me?

Thanks in advance.

matrix kronecker_product(matrix A, matrix B) {
    int m = rows(A);
    int n = cols(A);
    int p = rows(B);
    int q = cols(B);

    // Initialize the resulting matrix
    matrix[m * p, n * q] result;

    // Compute the Kronecker product
    for (i in 1:m) {
      for (j in 1:n) {
        for (k in 1:p) {
          for (l in 1:q) {
            result[(i - 1) * p + k, (j - 1) * q + l] = A[i, j] * B[k, l];
          }
        }
      }
    }
    return result;
  }
2

@mike-lawrence has a few implementations here: Kronecker Products · Issue #1454 · stan-dev/stanc3 · GitHub

2 Likes
3

Thanks!

4

You can’t really use those implementations at scale. And they traverse the matrices in non-local order.

The trick to working with Kronecker products is to never expand them out. You use the fact that things like inverses and Cholesky factors distribute, e.g.,

\text{cholesky}(A \otimes B) = \text{cholesky}(A) \otimes \text{cholesky}(B).

This is why we didn’t put them into Stan, because Stan always creates the full matrices it returns, so if we had a Kronecker product, it’d just be brute force like the implementations linked above.

Stan does allow that—the result is a submatrix of the original matrix. But it copies everything, so it’s inefficient.