Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Kronecker matrix utilities #60

Open
syclik opened this issue Jul 6, 2015 · 6 comments
Open

Kronecker matrix utilities #60

syclik opened this issue Jul 6, 2015 · 6 comments
Assignees
Milestone

Comments

@syclik
Copy link
Member

syclik commented Jul 6, 2015

From @bob-carpenter on January 3, 2015 2:22

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.

Copied from original issue: stan-dev/stan#1196

@syclik
Copy link
Member Author

syclik commented Jul 6, 2015

From @bgoodri on January 3, 2015 2:47

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.

@syclik
Copy link
Member Author

syclik commented Jul 6, 2015

From @flaxter on February 24, 2015 0:28

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.)

@syclik
Copy link
Member Author

syclik commented Jul 6, 2015

From @bgoodri on February 24, 2015 5:9

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.

@rtrangucci
Copy link
Contributor

@syclik @flaxter is there a branch open for this yet?

@syclik
Copy link
Member Author

syclik commented Sep 13, 2016

I don't have one.

@rtrangucci
Copy link
Contributor

I'm going open a new issue to add
vector kronecker_matrix_vector(matrix K1, matrix K2, vector V) to the library that references this issue rather than using this issue. Better to keep issues manageable.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

No branches or pull requests

2 participants