Speeding up Stan for large data sets and hard problems - approximations etc

[Julian King]

Hi all,

Stan is very useful because of its flexibility as an MCMC langauge.  It is possible to get something working quickly and reasonably effectively because the user has to express intentions in a fairly high-level way, whilst all the difficult sampling is taken care of. The more useful Stan can be in a wider range of situations, the better.

One of the reasons why one might move away from Stan to custom code is because the problem is "too hard" in some sense - perhaps the number of data points is too large and requires special techniques.  Perhaps there are linear algebra operations which scale as some poor function of the number of data points.  Gaussian Process methods are the best example of this, where naive implementations scale as O(N^3).

Perhaps it is useful to start a wishlist of features which can help extend the usefulness of Stan in this way, by making it feasible to a wider range of problems.

In the context of Gaussian Process regression, there are a number of ways of getting around the O(N^3) method, but one way is the use of a conjugate gradient approximation to find the approximate inverse of the covariance matrix.  Perhaps a Stan function which finds the approximate matrix inverse of a matrix C would be useful, eg approxinverse(matrix C, tolerance e), where e is the fractional tolerance.   If you are sampling the hyperparameters as well, you will need a way of getting the determinant of the covariance matrix fast as well, or the call to multi_normal will spend O(N^3) time calculating the determinant. 

Any other thoughts?

[N Sanders]

Hi Julian,

I think this is a great discussion to have, and I would add to the list the ability to optimize for models with "blocky" conditional independence structures. 

I've encountered this in non-linear hierarchical regression models, where bottom level parameters in different groups are defined to be conditionally independent of each other, yet Stan cannot know about this structure and may do unnecessary gradient computations.  It's tempting to use a Gibbs sampler in such cases, but of course it's still valuable to do joint transitions between the hyperparameters and lower level parameters that do have dependencies.

I'm not sure if this has been discussed on the list already, perhaps we've touched on it here.  I know Michael B. has thought about this and may be able to weigh in with a better description of the problem or comments on whether or not it's worth further discussion. 

Best,

~Nathan

[Bob Carpenter]

We're in the middle of writing a proposal for funding
to extend Stan to include, among other things, a number of
more scalable approaches, including:

* data-parallel EP

* data-streaming VB

* maximum marginal likelihood (MML)

* structured matrices (particularly banded and blocked)
for GPs and other applications

* sparse and ragged data structures

* stiff ODE solving

We're very open to suggestions.

Keep in mind that for both our approach to sampling (HMC) and
optimization (L-BFGS), we need gradients of any operations we
employ.

Even without moving to Riemannian HMC, which introduces local
covariance matrix scaling issues, we wouldn't mind being able to
robustly esitmate a global covariance matrix for Euclidean HMC.
L-BFGS does this in its own way already. And we'd very much like
to do that in a way that isn't O(n^3) for n x n covariance
matrices.

We've been encouraging people to use Cholesky-factor parameterizations
of covariance estimators as it cuts down on some of the algebra and
is more numerically stable.

- Bob

> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Michael Betancourt]

Bitter “old” man rant: there’s really no utility of having a discussion about
extending Stan without having actual goals in mind. EP, VP, and MML
can be great data exploration tools but they will never admit robust
uncertainty modeling in complex systems, nor are they designed to.
Full Bayes will always be the most challenging goal but also provide
the biggest gains; adding new data structures (sparse or banded
matrices, etc) and functions (stiff ODEs) will help certain problems
but there are many such improvements and we’ll have to take
them one at a time as we see the motivation (developer help is
always welcome!).

Consequently I think making vague statements is not particularly
helpful. If there are real models motivated by real problems that
don’t fit in Stan then let’s talk specifics. Otherwise it’s just
un-actionable speculation.

For example, when does the N^{3} scaling of a GP really start mattering?
When N is large? But when is N large? Are you trying to fit a GP over
many dimensions? What happens to GPs that try to fit on many-dimensional
latent spaces? Turns out lots of problems that can dwarf the N^{3}
complexity.

Stan has much room to grow — let’s do it right.

[Andrew Gelman]

Hi, Michael. Perhaps it would help in our proposal if we were to list specific examples of where we want these sorts of algorithms. One example is pharmacokinetic/pharmacodynamic models with hundreds or thousands of patients. These models can run slowly and it would be useful to have faster approximations. Another example is hierarchical logistic regressions such as what we use to do multilevel regression and poststratification for survey data. Again, full Bayes can be slow and people are running lmer (which has lots and lots of problems, most notably an extreme inflexibility when it comes to specifying priors or otherwise regularizing the inferences) as alternatives. Another example are Gaussian process decomposition models such as we used in the birthday example on the cover of BDA3. Again, this is a model we can write down but which we simply can’t fit in reasomable time in Stan. It’s a little embarrassing for the example on the cover of my book to not be fittable!
I hope this helps (both for you and for the proposal)
See you
A

On Oct 17, 2014, at 1:27 PM, Michael Betancourt <betan...@gmail.com> wrote:

> Bitter "old" man rant: there's really no utility of having a discussion about
> extending Stan without having actual goals in mind. EP, VP, and MML
> can be great data exploration tools but they will never admit robust
> uncertainty modeling in complex systems, nor are they designed to.
> Full Bayes will always be the most challenging goal but also provide
> the biggest gains; adding new data structures (sparse or banded
> matrices, etc) and functions (stiff ODEs) will help certain problems
> but there are many such improvements and we'll have to take
> them one at a time as we see the motivation (developer help is
> always welcome!).
>
> Consequently I think making vague statements is not particularly
> helpful. If there are real models motivated by real problems that
> don't fit in Stan then let's talk specifics. Otherwise it's just
> un-actionable speculation.
>
> For example, when does the N^{3} scaling of a GP really start mattering?
> When N is large? But when is N large? Are you trying to fit a GP over
> many dimensions? What happens to GPs that try to fit on many-dimensional
> latent spaces? Turns out lots of problems that can dwarf the N^{3}
> complexity.
>

> Stan has much room to grow -- let's do it right.