A few ideas for future versions...

[Arthur Breitman]

I thought I'd just share a few ideas for Stan. I'm still fairly new to the software, so I may misunderstand some of the features. Some of those ideas are likely to have been floated before, and I appreciate that the development team doesn't have unlimited resources.

That said, perhaps this can be the start of a useful conversation.

General
------------

- Support for sparse covariance matrices / band diagonal Cholesky matrices would be swell for implementing some Gaussian processes, such as AR(p). Alternatively, supporting vector and array slicing could confer the same benefit.

- Support for Daubechies wavelets and FFT is also useful to model large arrays with few coefficients.

- Some functions such as sqrt aren't yet vectorized.

Sampling
---------------

- l-bfgs-hmc looks quite promising and seems much cheaper than Riemannian HMC.

- parallel tempering is big in molecular dynamics where finding where the mass of the distribution lies is half the battle. Any thoughts on it?

Modeling language
-------------------------------

- Why not use name mangling to limit accidental collisions with C++ keywords?

- Parameters and generated quantities could be annotated to specify the sampler's behavior
(@all: keep track of variable / @m[k]: keep track of the first k moments only / @category: keep track of the marginal categorical distribution of the parameter / @ignore: ignore)

- If I'm not mistaken, upper and lower bounds only seem to apply during initialization and aren't enforced later on, unless T[] is used. The semantic is a little bit confusing as the bound is applied to the parameter's type, suggesting it's an intrinsic property. If it is only used for initialization, it should appear on the right hand side.

Python specific
-------------------------

- It would be very useful to be able to build a python representation of the model, before compilation, for debugging purposes. Beats "print".

- Is there a way to pass a logger to the sampling method instead of having it write to stdout by default?

Best,
Arthur

[Nigel Goodwin]

I've looked at L-BFGS-HMC. It is very good when there is not much curvature - i.e. the Hessian is fairly constant - but is not at all useful otherwise. It is quite easy to see why - the chains are in different parts of the space, spread out, so any approximate hessian is rather global, and doesn't reflect local behaviour. Imagine a strong banana, and the chains at different points in the banana.....

Maybe an adaptive approach which detects how much curvature there is, and uses L-BFGS-HMC when appropriate.....

My experience is that a bad approximate hessian can be worse than just a unit matrix or diagonal  Hessian.

[Bob Carpenter]

On Sep 9, 2014, at 1:32 PM, Nigel Goodwin <goo...@essence-property.com> wrote:

> I've looked at L-BFGS-HMC. It is very good when there is not much curvature - i.e. the Hessian is fairly constant - but is not at all useful otherwise. It is quite easy to see why - the chains are in different parts of the space, spread out, so any approximate hessian is rather global, and doesn't reflect local behaviour. Imagine a strong banana, and the chains at different points in the banana.....
>
> Maybe an adaptive approach which detects how much curvature there is, and uses L-BFGS-HMC when appropriate.....
>
> My experience is that a bad approximate hessian can be worse than just a unit matrix or diagonal Hessian.

We've found the same thing in estimating dense mass matrices,
even when the curvature isn't quite so varying as in the banana.

If the curvature is bad and varying over the posterior, you
really need something like Riemannian manifold HMC.

- Bob

[Bob Carpenter]

Thanks for posting to the group. I'll respond again for everyone.

On Sep 9, 2014, at 12:28 AM, Arthur Breitman <arthur....@gmail.com> wrote:

> I thought I'd just share a few ideas for Stan. I'm still fairly new to the software, so I may misunderstand some of the features. Some of those ideas are likely to have been floated before, and I appreciate that the development team doesn't have unlimited resources.
>
> That said, perhaps this can be the start of a useful conversation.
>
> General
> ------------
>
> - Support for sparse covariance matrices / band diagonal Cholesky matrices would be swell for implementing some Gaussian processes, such as AR(p). Alternatively, supporting vector and array slicing could confer the same benefit.

Yup. That one's on our radar. We do vector and matrix slicing and blocking, but
not general N-dimensional array operators. Kronecker products also come up in
some of these settings.

> - Support for Daubechies wavelets and FFT is also useful to model large arrays with few coefficients.

No idea what this is, but it probably explains why someone
was asking us if we had an FFT function in Stan itself.

Discrete FFT's are pretty easy to compute and Eigen's library is
solid. We're using it to compute autocorrelations in order to
estimate effective sample sizes.

I don't know how well auto-diffing Eigen's FFT algorithm would work.
Are there analytic gradients?

> - Some functions such as sqrt aren't yet vectorized.

Most functions aren't generalized to apply to vectors. That's a
convenience issue as it won't be any faster in Stan to vectorize.

> Sampling
> ---------------
>
> - l-bfgs-hmc looks quite promising and seems much cheaper than Riemannian HMC.

See other response.

> - parallel tempering is big in molecular dynamics where finding where the mass of the distribution lies is half the battle. Any thoughts on it?

See Michael Betancourt's arXiv paper on thermodynamic sampling.
We plan to implement that sooner rather than later. As soon as
I can generalize the model to split out the priors and likelihoods
(or bases and deltas as I believe Michael called them in the paper).

> Modeling language
> -------------------------------
>
> - Why not use name mangling to limit accidental collisions with C++ keywords?

This is a must-do that's on me to code up. I've been meaning to
do it for ages. It'll also prevent collisions with our own
built-in functions.

> - Parameters and generated quantities could be annotated to specify the sampler's behavior
> (@all: keep track of variable / @m[k]: keep track of the first k moments only / @category: keep track of the marginal categorical distribution of the parameter / @ignore: ignore)

We don't keep track of any moments directly (if by that you mean storing
X^2 as well as X --- they're easy to compute after the fact).

RStan lets you specify which variables to keep on the outside, but
we'll probably move to something more efficient with annotations like
you suggest in the model. We're torn as to whether this is really part
of the model, so annotations seem reasonable.

>
> - If I'm not mistaken, upper and lower bounds only seem to apply during initialization and aren't enforced later on, unless T[] is used. The semantic is a little bit confusing as the bound is applied to the parameter's type, suggesting it's an intrinsic property. If it is only used for initialization, it should appear on the right hand side.

I'm not sure what you mean here. If there are upper and lower bounds
declared on a parameter:

parameters {
real<lower=0,upper=1> theta;
}

then what happens is that we sample an underlying real variable that
is the logit-transform of theta (i.e., log(theta/(1-theta))). Then
we inverse transform with inverse logit and apply the Jacobian to get
probabilities. It happens every time the log probability function is
called.

Bounds on data and transformed data are validated once when the data's
read in for the data block and after the definitions in the transformed
data block. Transformed parameters and generated quantities also both test
after their blocks have finished executing.

There are no constraints allowed on local variables because there's
no good way to keep them always satisfying their constraints.

>
> Python specific
> -------------------------
>
> - It would be very useful to be able to build a python representation of the model, before compilation, for debugging purposes. Beats "print".

You mean of the Stan program? Stan isn't restricted
to graphical models like BUGS or JAGS.

There's an underlying abstract syntax tree that's
generated by the parser for the Stan language and then
used as the basis for C++ code generation.

> - Is there a way to pass a logger to the sampling method instead of having it write to stdout by default?

No idea. Allen should know, though, so if he doesn't catch this note,
take it up again in an appropriately subjected message and I'm sure he'll
respond.

- Bob

[Allen Riddell]

On 09/09, Arthur Breitman wrote:
> Python specific
> -------------------------
>
> - It would be very useful to be able to build a python representation of
> the model, before compilation, for debugging purposes. Beats "print".

Can you explain what you'd like in a bit more detail? I'm intrigued.

> - Is there a way to pass a logger to the sampling method instead of having
> it write to stdout by default?

I did some research on this awhile back and I didn't find an easy solution. It
would certainly be nice. I'll create an issue to track this.

Thanks for the suggestions!

Best wishes,

-a

[Marcus Brubaker]

> - Support for Daubechies wavelets and FFT is also useful to model large arrays with few coefficients.

No idea what this is, but it probably explains why someone
was asking us if we had an FFT function in Stan itself.

Discrete FFT's are pretty easy to compute and Eigen's library is
solid.  We're using it to compute autocorrelations in order to
estimate effective sample sizes.

I don't know how well auto-diffing Eigen's FFT algorithm would work.
Are there analytic gradients?

FFT is easy as it is simply a linear operator.  Maybe we can add implementing FFT to my plate come October.

 

> - Some functions such as sqrt aren't yet vectorized.

Most functions aren't generalized to apply to vectors.  That's a
convenience issue as it won't be any faster in Stan to vectorize.

That said, if someone wanted to get their feet wet contributing to Stan, this would be a really natural way to do it.  Add elementwise vectorized versions of some of these functions, write unit tests, update docs and submit a pull request.  All of the pieces of contributing, but the technical component requires very little knowledge of the code.

 

> Sampling
> ---------------
>
> - l-bfgs-hmc looks quite promising and seems much cheaper than Riemannian HMC.

See other response.

Just to add to that, my memory of the l-bfgs-hmc paper is that it reads a lot more promising than it actually is.  As someone else noted, the overdispersed nature of the samples means that the estimated Hessian will be terrible on anything with non-constant curvature.

Cheers,

Marcus

[Arthur B.]

Thanks for your replies, my own replies are inlined

On Tuesday, September 9, 2014 1:49:10 PM UTC-4, Bob Carpenter wrote:

Thanks for posting to the group.  I'll respond again for everyone.

On Sep 9, 2014, at 12:28 AM, Arthur Breitman <arthur....@gmail.com> wrote:

> - Support for Daubechies wavelets and FFT is also useful to model large arrays with few coefficients.

No idea what this is, but it probably explains why someone
was asking us if we had an FFT function in Stan itself.

Well the simplest Daubechies wavelet would be the Haar transform. You can get creative with those,

for example you could impose a sparse prior on the Haar decomposition of the coefficients of a linear

regression. 

Discrete FFT's are pretty easy to compute and Eigen's library is
solid.  We're using it to compute autocorrelations in order to
estimate effective sample sizes.

I don't know how well auto-diffing Eigen's FFT algorithm would work.
Are there analytic gradients?

You bet, dX_k/dx_n = exp(-i 2 pi k n /N)... remember, it's just a linear operator!

  
> - Some functions such as sqrt aren't yet vectorized.

Most functions aren't generalized to apply to vectors.  That's a
convenience issue as it won't be any faster in Stan to vectorize.

Indeed, but suppose for example you're doing a bayesian neural net. It's very natural

to represent it using weight matrices and activation vectors... however at some point

you're going to need to do

for (i in 1:a.size())

   a[i] <- tanh( a[i] )

which isn't as a readable as tanh(a).

I looked into the code, and it looks that it's not super duper easy to get it for free unfortunately. It seems that

the exp function when called on a vector comes from the Eigen library and has a different implementation from the exp

function called on a scalar. I do think it should be possible to make all scalar function accept vector and matrix types

directly with enough templates.

 

> - Parameters and generated quantities could be annotated to specify the sampler's behavior
> (@all: keep track of variable / @m[k]: keep track of the first k moments only / @category: keep track of the marginal categorical distribution of the parameter / @ignore: ignore)

We don't keep track of any moments directly (if by that you mean storing
X^2 as well as X --- they're easy to compute after the fact).

No, I meant storing a running sum of X^2, because X is a vector with one million elements and you can't afford to record

every single one of them. 

 

I'm not sure what you mean here.  If there are upper and lower bounds
declared on a parameter:

I think I need to reread the documentation on bounds, I probably don't understand them correctly. 

 
> Python specific
> -------------------------
>
> - It would be very useful to be able to build a python representation of the model, before compilation, for debugging purposes. Beats "print".

You mean of the Stan program?  Stan isn't restricted
to graphical models like BUGS or JAGS.

No, I mean have a model interpreter written in python that replicates the behavior of the one compiled in C++. For instance

m = pystan.model('blah.model') 

# let's see if it's reasonable

print m.log_likelihood( data, parameters = m.sample_params() )

print m.sample( data, n=10 )

# ok it's all good, now compile to C++

m.compile() 

[Bob Carpenter]

On Sep 9, 2014, at 8:45 PM, Arthur B. <arthur....@gmail.com> wrote:

> Thanks for your replies, my own replies are inlined
>
> On Tuesday, September 9, 2014 1:49:10 PM UTC-4, Bob Carpenter wrote:
> Thanks for posting to the group. I'll respond again for everyone.
>
> On Sep 9, 2014, at 12:28 AM, Arthur Breitman <arthur....@gmail.com> wrote:

...

> Discrete FFT's are pretty easy to compute and Eigen's library is
> solid. We're using it to compute autocorrelations in order to
> estimate effective sample sizes.
>
> I don't know how well auto-diffing Eigen's FFT algorithm would work.
> Are there analytic gradients?
>
> You bet, dX_k/dx_n = exp(-i 2 pi k n /N)... remember, it's just a linear operator!

Cool. I think Marcus already said he could add it.
He understands all this much better than I do.

> > - Some functions such as sqrt aren't yet vectorized.
>
> Most functions aren't generalized to apply to vectors. That's a
> convenience issue as it won't be any faster in Stan to vectorize.
>
> Indeed, but suppose for example you're doing a bayesian neural net. It's very natural
> to represent it using weight matrices and activation vectors... however at some point
> you're going to need to do
>
> for (i in 1:a.size())
> a[i] <- tanh( a[i] )
>
> which isn't as a readable as tanh(a).

Agreed. Vectorizing everything's on our to-do list.

> I looked into the code, and it looks that it's not super duper easy to get it for free unfortunately. It seems that
> the exp function when called on a vector comes from the Eigen library and has a different implementation from the exp
> function called on a scalar. I do think it should be possible to make all scalar function accept vector and matrix types
> directly with enough templates.

With what we now know about template programming, it'd be
easy to vectorize just about all of our functions pretty
easily.

Eigen using a different exp() implementation shouldn't matter;
it's just an implementation detail and we compute the gradients
the same way either way.

> > - Parameters and generated quantities could be annotated to specify the sampler's behavior
> > (@all: keep track of variable / @m[k]: keep track of the first k moments only / @category: keep track of the marginal categorical distribution of the parameter / @ignore: ignore)
>
> We don't keep track of any moments directly (if by that you mean storing
> X^2 as well as X --- they're easy to compute after the fact).
>
> No, I meant storing a running sum of X^2, because X is a vector with one million elements and you can't afford to record
> every single one of them.

We've thought about this but in the end decided we weren't
likely to run into this problem. We may have to rethink as
Stan gets more scalable.

We don't anticipate people storing a million iterations!
The bigger deal is the width --- how many parameters do they
store.

We know how to do online mean and (co)variance calculations.
Quantiles are much harder (e.g., 10% quantiles require 10%
of the draws to be stored at a time). Boost has some nice
streaming accumulators.

> > Python specific
> > -------------------------
> >
> > - It would be very useful to be able to build a python representation of the model, before compilation, for debugging purposes. Beats "print".
>
> You mean of the Stan program? Stan isn't restricted
> to graphical models like BUGS or JAGS.
>
> No, I mean have a model interpreter written in python that replicates the behavior of the one compiled in C++. For instance

It'd be nice to have alternative implementations of the
Stan language, but an interpreted version in Python would
be a prohibitively large undertaking.

- Bob

[Arthur B.]

On Wednesday, September 10, 2014 1:28:47 AM UTC-4, Bob Carpenter wrote:

Agreed.  Vectorizing everything's on our to-do list.

Wait, that  .size() function, though mentioned in the documentation on page 213 - isn't actually there... 

I can't do  inv_logit( v )

I can't do 1 /. (1 + exp( -v ))  ( can't divide scalar by vector )

I can't do size( v ) because size only applies to arrays.

I can do exp(-log(1+exp(-v))) but that's such a waste of computation !

I can pass a length parameter known at compile time, but it's clunky

Or I can convert to an array first, but I don't know what the overhead is on that...z

Am I missing another way?

[Bob Carpenter]

On Sep 11, 2014, at 11:43 AM, Arthur B. <arthur....@gmail.com> wrote:

> On Wednesday, September 10, 2014 1:28:47 AM UTC-4, Bob Carpenter wrote:
>
> Agreed. Vectorizing everything's on our to-do list.
>
>
> Wait, that .size() function, though mentioned in the documentation on page 213 - isn't actually there...

No idea how that creeped in. I'm supposed to be making
sure all these functions parse! There's a size() function
for arrays and then rows() and columns() for matrices and
vectors. But nothing with syntax like .size()!

> I can't do inv_logit( v )

Not yet.

> I can't do 1 /. (1 + exp( -v )) ( can't divide scalar by vector )

You should be able to do this.

> I can't do size( v ) because size only applies to arrays.

cols(v) is what you want if v is a column vector.

We should probably extend size() to vectors (and even to
matrices, where size is (M x N) --- Eigen lets you index
matrices in column-major order with a single index, which
might also be convenient in places in Stan).

> I can do exp(-log(1+exp(-v))) but that's such a waste of computation !
> I can pass a length parameter known at compile time, but it's clunky
> Or I can convert to an array first, but I don't know what the overhead is on that...z

Not much overhead converting to an array. The bigger problem is
the clutter and lack of readability in the model. Most of Stan's time,
either for optimization or sampling, is spent in autodiff calculating
gradients.

> Am I missing another way?

But you can use the logit-parameterized log
probability functions, like bernoulli_logit and
binomial_logit, and equivalently categorical_logit. And
you can append _log to that to get the log probability
function if you need it in a mixture model or something.

- Bob

[Arthur B.]

Here's another potentially interesting feature.

Computing the likelihood, let alone the gradient, can be quite expensive in large datasets. Linear regression models on large datasets have benefited greatly from stochastic gradient descent. It turns out something similar can be achieved in MCMC. Instead of computing the likelihood of an entire i.i.d sample, this likelihood can be approximated by using a small, random, subsample.

Using the generalized poisson estimator introduced in Fearnhead, Papaspiliopoulos and Roberts [1], a subset of the data can be used as an unbiased estimate of the log likelihood and provide and unbiased acceptance / rejection step. 

Quiroz, Villani, and Kohn [2] suggest doing a simple convexity adjustment by assuming normality of the log likelihood distribution across the sample, and Barnet, Doucet, and Holmes [3] seem to take another approach altogether, but with the same goal of using a subsample.

In theory, this should translate to an HMC sampler, though it might decrease the acceptance rate. 

[1] http://arxiv.org/pdf/0710.4245.pdf

[2] http://arxiv.org/pdf/1404.4178v1.pdf
[3] http://jmlr.org/proceedings/papers/v32/bardenet14.pdf

[Bob Carpenter]

Thanks for the links.

For scalability, we're currently pursuing black-box VB and
EP algorithms (as well as maximum marginal likelihood).

One issue we have to deal with for the data subsampling
algorithms is linking parameters to subsets of the data,
as you would need to scale something like an IRT model or
really any model with random effects that grow with the data.

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

[Arthur Breitman]

While we're at it, there's also the crossentropy method.

If you have a good analytical approximation of your posterior for which maximum likelihood is tractable, you can alternate between importance sampling, and fitting it to the drawn samples.

It has the advantage that it's a global approach, it has the disadvantage that, well, the posterior is probably not well approximated by anything tractable.

You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/Dfa_mV3MyXA/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

[Andrew Gelman]

Yes, that’s what we’re doing with our (in-progress) EP implementation. We’re focusing on hierarchical models where the data subsetting corresponds to a parameter subsetting, thus reducing the number of parameters as well as the number of data points in each subset.

[Bob Carpenter]

Ditto for stochastic VB. EP is nice because it'll be data
parallel, whereas stochastic VB is nice because it is streaming.
Streaming can potentially converge before even looking at all the
data, as in Matt Hoffman's original Wikipedia LDA example.

- Bob

[Arthur Breitman]

If the model is a simple expression, the irrelevant parameters for a subset could be detected by making a fast dummy pass though the likelihood function and tainting the parameters as we go.

With if statements, recursion, etc, a real pass can be made without AD, followed by an AD pass that knows the support of the gradient.

After the parameters are updated (for a given batch), the support is unlikely to have changed, and it can be used as a guess for the next evaluation, refining it with more (or less) parameters if inadequate.

It gets murkier when working with matrices though.

You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/Dfa_mV3MyXA/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

[Daniel Lee]

Hi Arthur,

Can you provide some context to what you've said? I'd like to understand what you've said, but I think there are just too many terms that I don't know. For example:

- irrelevant parameters

- irrelevant parameters for a subset (perhaps if I understand what irrelevant meant, this would be clear)

- detected

- fast dummy pass

- tainting the parameters
...

Perhaps if I knew where you were coming from, that would put enough context around your point where I can follow along.

Daniel

[Arthur Breitman]

I'm replying to on old thread, which may provide more context.

The idea is that some techniques allow one to sample by computing the likelihood function on random subsets of the data. For this to be efficient however, one needs to know how parameters relate to these subsets.

One approach is to explicitly index some parameters (typically some hidden variables) along with the observations. However, I was pointing out that such a relationship could be detected automatically.

[Bob Carpenter]

This is the same data slicing idea we need for
EP and VB.

We've been thinking along similar lines for
detecting Markov blankets of discrete parameters.

As long as we (1) pin down the conditionals in the language
to not touch parameters, and (2) don't do conditionals inside
our autodiff function implementations, we get a static expression
tree.

It's very tempting to start thinking about code-generating
autodiff at that point.

But what we'll probably do is somehow link the parameters to
the data in the language. Haven't figured out how we'll annotate
that in the log density function though.

- Bob