Representing functions that are best presented programmatically

Modeling
1

This feels to me like a potential beginner question, but I have little experience of probabilistic programming to understand the right route though (though I have lots of experience of programming, as well as reasonable experience of Bayesian modelling and methods).

My model can be presented in matrix-vector form, but the matrices are huge but sparse, with the structure encoded such that the various matrices can be implemented through a fast compressed implementation. There are various different matrices that are in sparse in different bases (for example, I have a skyline matrix, and a matrix which is diagonal with sparse elements in the DT-CWT domain, but is necessarily transformed in the model).

All these transformations are implemented programmatically with the full matrices never being explicitly calculated and stored, only their action on an input.

I daresay it is possible to implement the full set of operations in Stan, but there are various subtleties that are tricky to get right. If one proceeds down that route, is it possible to verify the implementation on actual data?

Ideally, I’d be able plug in the function as a black box transform in whatever language. Is this possible in Stan?

The ultimate objective is to do variational inference, and it would save lots of pain if there is some automated technique that can bypass a lot of paper analysis.

2

Yes. You can write functions in Stan and test them. There are functions in RStan that make this easy, but it’s not too hard to do manually elsewhere. You can test values by applying in transformed data and you can test it works with autodiff by pulling out gradients (though you’d need to have a way to test they were correct). We have sparse matrix multiply implemented in CSR form already.

You can also write templated functions in C++ and plug them into Stan. You can generate some functions in Stan and see what the generated C++ looks like in the .hpp file that’s generated from compiling a Stan file. But it’s a bit of a pain to get external functions hooked up, especially if you want to supply custom gradients for efficiency.