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.