structured linear algebra and GP covariance functions

[bob-carpenter]

Moved from former wiki page of the same title (Wiki page was removed).

Introduction

Stan's linear algebra functions assume dense matrices with no special structure. We intend to add specialized matrix types in order to speed up matrix computations when matrices adhere to specific structures such as symmetric positive definite, Toeplitz, banded, and sparse matrices. As the new types are added, we will specialize the GP covariance functions as well to return the appropriate structured matrix type based on user specifications. We do not plan to do automated detection of specialized structure.

TODO

The current to do list is as follows:

• Implement symmetric, positive definite Toeplitz type (can be banded SPD Toeplitz)
• Specialize addition, matrix multiplication for SPD Toeplitz type
• Specialize Cholesky and inverse for SPD Toeplitz type
• Specialize log determinant for SPD Toeplitz
• Specialize cov_exp_quad to return SPD Toeplitz
• Add Kronecker and Hadamard products for arbitrary matrices
• Implement SPD type
• Specialize cov_exp_quad to return SPD
• Specialize Eigendecomposition for SPD matrices
• Implement sparse matrix type
• Implement SPD sparse matrix type
• Specialize matrix algebra functions for sparse matrix type
• Specialize matrix algebra functions for SPD sparse matrix type
• Specialize decompositions for sparse matrix type
• Specialize decompositions for SPD sparse matrix type

[spinkney]

I posted on the Stan blog about correlation matrices and the second part will arrive this week. I took the code in corr_constrain and modified it in Stan for my use case. Two extensions that would be really nice to have.

1. I wanted to go from the unconstrained space to the constrained and so needed to add to the lp but I wanted this in transformed parameters. Luckily, the first value of a cholesky factor correlation is 1 so I could store the log-abs-jacobian there, pass it to a real var, put 1 back in place and head to the model block. It would be nice to be able to calltarget += within a function in transformed parameters.

2. Creating a bunch of corr_constrain functions for known correlation types that are exposed to Stan. My next post shows the block structure I'm interested in (block compound symmetric). There's also Toeplitz (and an open issue for that https://github.com/stan-dev/math/issues/356) and first order AR. See https://blogs.sas.com/content/iml/2012/11/05/constructing-common-covariance-structures.html. I'm sure there are plenty of others. Building the template for a couple of these would be awesome and make it easier to add more in the future plus a benefit for time-series models. Lastly, the log-abs-jacobian is the same calculation which is nice.

3. I know I said two. This is related to having ragged-arrays. Slicing and indexing on the row scan for triangular matrices with block structure (where the blocks are possibly varying sizes) is a nightmare. Hopefully this comes soon!

This issue is possibly related to https://github.com/stan-dev/math/issues/2189#issue-741349189