Kronecker matrix utilities

[bob-carpenter]

Need utilities for Kroneckers for Gaussian processes. I don't know much about it but want to get it down as an issue and point to Ben's R example:

https://groups.google.com/forum/#!msg/stan-users/iddShgJVwTQ/QgcDEP3h9JQJ

Talk to Ben to figure out how this should look in Stan.

[bgoodri]

I think the most immediate need is to have a user-facing function called something like quad_form_kronecker that takes 3 arguments (matrix X, matrix Y, and vector or matrix Z) and outputs the equivalent of Z' * kronecker_product(X, Y) * Z. Algebraiclly, that is a lot of subvector * scalar * matrix * subvector but with a smaller number of unique matrix * subvector expressions that could be performed once each internally.

It is often the case that either X or Y is diagonal or otherwise sparse, in which case there could be some more optimizations. I'm not sure we need any more than that because it is rare in statistics to see a Kronecker product that is not embedded in a quadratic form.

[flaxter]

Here's a summary of the main points from the Stan dev meeting on this topic (17 February 2015):

Eigen now has (experimental) support for Kronecker products: http://eigen.tuxfamily.org/dox-devel/unsupported/group__KroneckerProduct__Module.html and Eigen also has a fair bit of supported sparse types (some experimental) e.g.: http://eigen.tuxfamily.org/bz/show_bug.cgi?id=275 http://eigen.tuxfamily.org/dox/group__SparseQuickRefPage.html
http://eigen.tuxfamily.org/dox/group__Sparse__chapter.html
http://eigen.tuxfamily.org/dox-devel/group__SparseCholesky__Module.html

We discussed ways of implementing, a type kronecker_matrix (which would also be useful for new types tridiagonal_matrix, sparse_matrix and toeplitz_matrix and symmetric_matrix) which would be similar to cholesky_factor_cov.

Important design decision
If K = K1 ⊗ K2 ⊗ K3 and K1, K2, and K3 are different sizes (let's say n x n, m x m, and o x o) we need a new data structure to be able to handle this. There was a discussion of using view or map in C++/Eigen, with a vector of matrices all sized N x N where N = max(n,m,o). @bob-carpenter needs to weigh in on this, and others should make their cases.

Once we get that worked out, the main wish list for Kronecker matrix functions is simply everything in section 34.7 of the manual, Linear Algebra Functions and Solvers, e.g.:
row_vector operator/(row_vector b, kronecker_matrix A)
row_vector mdivide_right_kronecker(row_vector B, kronecker_matrix A)
vector eigenvalues_sym(kronecker_matrix A)

We also need matrix-vector multiplication (not sure the naming convention here):
row_vector matrix_multiply_vector(kronecker_matrix A, vector B)
row_vector vector_multiply_matrix(row_vector B, kronecker_matrix A)

I guess we need a function (although like inverse it shouldn't be the first choice for a user):
matrix kronecker_product(kronecker_matrix A)

Also very useful will be:
y ~ multi_normal_kronecker(vector mu, kronecker_matrix A)

Relatedly, the matrix normal distribution is already implemented, but not exposed. It should be exposed. (@mbrubake may have thoughts.)

[bgoodri]

As a start, we might think about exposing a Triplet http://eigen.tuxfamily.org/dox/classEigen_1_1Triplet.html That would make it possible for a C++ function to take a Triplet and do sparse stuff internally. In addition, that would be useful for (non-)missing data.

[syclik]

This issue was moved to stan-dev/math#60