advi feedback from users

[Alp Kucukelbir]

hi folks,

i'm planning on making some updates to the advi codebase in the next few weeks and would love your feedback.

those of you who have used advi (successfully or not) in the past: would you mind sharing your experiences?

+ what kind of data / model did you use advi on?

+ how big was your dataset? (alternatively: what made you want to try advi?)

+ were the default settings sufficient? did you have to play with any of the settings?

+ did you try the fullrank version of advi? 

 

i'm looking forward to making advi easier and more reliable for everyone to use.

many thanks in advance.

cheers

alp

[Stephen Martin]

I tried using it on a fairly weird model of mine.
The data are dichotomous responses assumed to be generated from a 2PL model.
However, the purpose of the model is not to just estimate a 2PL IRT model, but to detect latent groups which some items exhibit DIF between.
So it's a fairly convoluted model --- The likelihood for each item and case is a mixture of bernoulli likelihoods, one which assumes DIF one which does not; when it assumes DIF, it is another mixture model aiming to assign individuals to latent clusters.
Moreover, one can 'predict' latent cluster membership using softmax or logistic regression from other respondent data.

Needless to say, ADVI did not work well. Optimizing doesn't either, presumably because this sort of problem needs an EM approach, not a basic ML approach.

I wanted to try ADVI because the model takes about 1200 seconds for 300 cases and 24 items.
I did try the fullrank, it didn't work well either.

By 'didnt work well', I mean it completely fails to estimate the correct latent group memberships (which I know, because I generate the data) and fails to discern the DIF items from non-DIF items.

[Janne Sinkkonen]

I have tried ADVI lately on three different kinds of models:

- Small negative binomial regressions, just out of curiosity. Full rank, didn't always estimate the theta parameter well. Not many replications to estimate it, so I'm not surprised ADVI has problems. 

- A simple non-negative regression model with fairly big count data (~ 100000, seven covariates), where I needed only point estimates for a recommender system. Full rank, works as expected. 

- A nonnegative component model for count data, with some hierarchy in the model. 240k counts, 10k parameters. Mean field, works well and is surprisingly fast with the new indexing notation (4 minutes). Our goal is not to make exact inference, but to find descriptive components. Goes to local minima etc., but this is not ADVI's fault. Sometimes throws an overflow, don't remember which function it comes from. 

A colleague has tested mean-field ADVI for a Poisson Factorisation model, with about 3M lines of data and at least the same number of parameters. (log-)Gaussian may not fit well with the posterior, or maybe the default gamma parametrisation was bad for mean field. 

All these are quick tests made in the context of commercial projects or theses, so we don't have definite reports or papers to show. 

[Gilles Boulianne]

It would be really nice to have an easy way of implementing mini-batches (even nicer if it could be done in pyStan). 

The method described in https://github.com/stan-dev/stan/blob/adsvi/how_to_ADSVI.md works for me, but I automated it with a script that must know data variable names, and the right place in the compiled model's .hpp file, etc. so it is ugly and difficult to maintain when trying several models.

[John Hall]

I'm sort of new to playing around with ADVI.

+ what kind of data / model did you use advi on?

Simulated data. The model was a one-dimensional PPCA with a constraint on the weight matrix to be a unit vector.

+ how big was your dataset? (alternatively: what made you want to try advi?)

240X200, but I also tried 240X2000 on the mean field and it was still pretty fast

+ were the default settings sufficient? did you have to play with any of the settings?

I found myself often decreasing the tol_rel_obj to try to increase the accuracy.

+ did you try the fullrank version of advi?

Yes (see below with other comments).

I was really impressed with the default mean-field option and how fast it is. However, when playing with the fullrank option, it sometimes will diverge for a while (delta_ELBO_mean will be decreasing steadily and then suddenly jump to like 1 or so for a while before declining) and then take way longer to get to a reasonable convergence, like 60,000+ iterations, depending on the model. Even then, I still feel like running it longer to get more accurate results.

There was a thread earlier about accessing the fullrank covariance matrix. I may not be interested in the actual content of the matrix, but I could imagine that I would be interested in the eigenvalues of the matrix. For instance, if the distribution of the eigenvalues is consistent with what is suggested by random matrix theory, then there might not be any advantage to using fullrank vs. mean field.

I sometimes get the error

Error in sampler$call_sampler(c(args, dotlist)) : 
  stan::variational::advi::adapt_eta: All proposed step-sizes failed. Your model may be either severely ill-conditioned or misspecified.

after eta adaptation, but then I just run it again and again until it stops.

Finally, the sampler says

We recommend genuine 'sampling' from the posterior distribution for final inferences!

at the end. But what if the reason you're using ADVI is because Stan has a hard time fitting the distribution? I could use this output to init a call to Stan, but are there any other options?

On Saturday, November 19, 2016 at 3:53:06 PM UTC-5, Alp Kucukelbir wrote:

[Bob Carpenter]

> On Feb 14, 2017, at 2:55 PM, John Hall <john.mic...@gmail.com> wrote:
>
> I'm sort of new to playing around with ADVI.
>

...

Thanks for the feedback.

> Finally, the sampler says
> We recommend genuine 'sampling' from the posterior distribution for final inferences!
> at the end. But what if the reason you're using ADVI is because Stan has a hard time fitting the distribution? I could use this output to init a call to Stan, but are there any other options?

Is the reason NUTS is having a hard time due to run time and
scaling?

The reason this recommendation is in there is that we don't know
the conditions under which ADVI fails. Thus we at least
suggest testing with simulated data to make sure you're recovering
the parameters you think you are in a simulated case.

- Bob

[John Hall]

Is the reason NUTS is having a hard time due to run time and
scaling?  

The model (below) is similar to the PPCA model in the ADVI paper. So some of the reasons why ADVI worked well for those models sort of applies here. Stan fits it with NUTS pretty quickly actually. The problem is that N_eff is low for some of the variables, like 3 for mu_F and 10 for lp__ out of 1000 in my last run (also rhat is above 2 for mu_F and above 1.2 for lp__). One of the best things I did when I was just playing around with NUTS was constraining B to be a unit_vector. I think this is because it helps lock down what B is in some sort of range and helps F get identified.

data {
    //sizes
    int<lower=1> T;
    int<lower=1> N;
   
    //data
    matrix[T, N] Y;
}
transformed data {
    real<lower=0> one_over_sqrt_N;
   
    one_over_sqrt_N = 1 / sqrt(N);
}
parameters {
    real mu_F;
    real<lower=0> sd_F;
    vector[T] F;
    real<lower=0> sd_B;
    unit_vector[N] B;

    real<lower=0> sd_e;
}
model {
    mu_F ~ normal(0, 0.25);
    sd_F ~ normal(0, 6);
    F ~ normal(mu_F, sd_F);
    sd_B ~ normal(0, one_over_sqrt_N);
    B ~ normal(0, sd_B);

    to_vector(Y) ~ normal(to_vector(F * B'), sd_e);

[Ben Goodrich]

On Tuesday, February 14, 2017 at 4:46:02 PM UTC-5, John Hall wrote:

The model (below) is similar to the PPCA model in the ADVI paper. So some of the reasons why ADVI worked well for those models sort of applies here. Stan fits it with NUTS pretty quickly actually. The problem is that N_eff is low for some of the variables, like 3 for mu_F and 10 for lp__

When you get an estimated N_eff of 3, that means either the convergence of chains did not occur or they did not even move appreciably. But neither of those will happen if the posterior distribution is even roughly multivariate normal, which is what full-rank ADVI is using to approximate the posterior. In the meanfield case, the approximation is multivariate normal but with dimensions that are independent.

[Michael Betancourt]

Having 3 / 1000 effective sample per transition is hinting that NUTS

is encountering really strong correlations (i.e. your model is only

weakly identified) and you may need to add additional regularization

or reparameterize your model.

But these are exactly the problematic features that will cause ADVI

to return inaccurate results.  Remember that we have never found

a situation where ADVI returns the right answer but NUTS does not,

so when NUTS indicates an issue ADVI is especially suspect.

--
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.
To post to this group, send email to stan-...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[John Hall]

If this is a problem for the model I posted above, then it is also a problem for the PPCA/ARD model in Figure 23 in the ADVI paper. I just ran NUTS on that exact model with my simulated data. With M=1, I got 2 / 1000 effective samples for the alpha and lp__ parameters (rhats around 10). With M=2, I got 2 / 1000 for both alphas (and rhats > 10), but a reasonable number for lp__.

I take your point about really strong correlations and weak identification seriously. I had realized that this was the problem with fitting the model with NUTS. The problem is that z and w in the Figure 23 model are correlated.

I had assumed ADVI would be better able to handle this type of problem, especially since I knew that my model was essentially the same as what was in the paper and the paper was using this as a potential use-case of ADVI.

To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.

To post to this group, send email to stan-...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

--
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/AZIlR7qfqLY/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

[Bob Carpenter]

> On Feb 15, 2017, at 11:44 AM, John Hall <john.mic...@gmail.com> wrote:
>
> If this is a problem for the model I posted above, then it is also a problem for the PPCA/ARD model in Figure 23 in the ADVI paper.

Correct.

I'm looking at the most recent incarnation of that paper, a
draft submitted (?) to JMLR (always hard to tell when authors
put preprints in the format of a journal):

http://www.stat.columbia.edu/~gelman/research/unpublished/advi_journal.pdf

I also wouldn't recommend the models in that paper --- they don't follow
our prior recommendations or our parameterization recommendations,
despite Andrew's being listed as a co-author. I'm also amused that it
uses the term "automatic relevance determination", because Andrew objects
every time it's used in his presence.

> I just ran NUTS on that exact model with my simulated data. With M=1, I got 2 / 1000 effective samples for the alpha and lp__ parameters (rhats around 10). With M=2, I got 2 / 1000 for both alphas (and rhats > 10), but a reasonable number for lp__.

Might this be a label-switching problem? I don't
know what PPCA/ARD is, but if it's some kind of factor
model, then that can be an issue.

> I take your point about really strong correlations and weak identification seriously. I had realized that this was the problem with fitting the model with NUTS. The problem is that z and w in the Figure 23 model are correlated.
>
> I had assumed ADVI would be better able to handle this type of problem, especially since I knew that my model was essentially the same as what was in the paper and the paper was using this as a potential use-case of ADVI.

Be careful what the paper evaluates. Like most machine learning
papers, it's focused on prediction, not parameter estimation and
posterior uncertainty calibration. Some of the earlier plots show
how these things go wrong in simple cases (like figure 5) and right
in others (at least w.r.t. the posterior mean---we know the posterior
variance is underestimated by construction of the KL divergence measure
being optimized).

They do the predictive accuracy properly, using full Bayes for sampling
and approximate Bayes through ADVI for VI.

- Bob

[John Hall]

I have finally got to something that fits my data with reasonable N_eff using NUTS (300 of 1000 for lp__, 1000 of 1000 for all hyperparameters, rhats between 1 and 1.01) and thought I would share (below). ADVI also fits it reasonably well. Of the hyperparameters, sd_z was higher than NUTS and the standard deviations of mu_z and sd_z were bigger than NUTS. But overall, it looks like it fits it well.

As many folks have noted, identification is very important for these PPCA-like models. What I really needed was a way to constrain w to be a positive unit normal vector. This isn't in Stan. However, I was able to make a close enough approximation by noting that var(w) should equal 1/N - mean(w)^2 for all unit vectors. The approach I took below uses this formula to constrain the mean. Since I didn't mind that w>0, it worked out pretty well. I also put a reasonably strong prior on the standard deviation of w, combined with constraints so mean(w) didn't have any imaginary variables in it.

Another thing I tried was to constrain z to have mean of 0 AND standard deviation of 1. The Stan documentation does a good job explaining how to force mean of 0, but it's harder than you would think to also force standard deviation of 1 while having the whole thing be identified properly (The only way I was able to do it resulted in lp__ near zero because the raw z values were not identified. The approach usually fit the data well, but it can be hard to diagnose problems when you're not printing all the variables). I tried holding back a second data point and did the math so that it works out, but it's very easy to get imaginary terms that screw it all up. So I gave up on that.

I would really love if Stan had better support for these kinds of constraints. Positive unit vectors seems do-able, but I don't know about forced unit normals.

data {
    int<lower=0> T;
    int<lower=0> N;
    vector[T] x[N];
}
transformed data {
    real<lower=0> one_over_N;
    real<lower=0> one_over_square_root_N;
    real<lower=0> sd_w_prior;
   
    one_over_N = 1.0 / N;
    one_over_square_root_N = sqrt(one_over_N);
    sd_w_prior = 1 / (sqrt(2 * N) * inv_Phi(0.95)); //selected so 95% of
                                                    //distribution is below
                                                    //(1 / sqrt(N))
}
parameters {
    real mu_z;
    real<lower=0> sd_z;
    vector[T] z;
    real<lower=0, upper=one_over_square_root_N> sd_w;
    vector<lower=0, upper=1>[N] w;
    real<lower=0> sd_e;
}
transformed parameters {
    real<lower=0, upper=one_over_square_root_N> mu_w;

    mu_w = sqrt(one_over_N - square(sd_w));
}
model {
    mu_z ~ normal(0, 0.5);
    sd_z ~ normal(1, 0.5);
    z ~ normal(mu_z, sd_z);
    sd_w ~ normal(0, sd_w_prior);
    w ~ normal(mu_w, sd_w);
    sd_e ~ normal(0, 0.05);
   
    for (n in 1:N)
        x[n] ~ normal(z * w[n], sd_e);

[Bob Carpenter]

Thanks for reporting back. That's great to know that it lines
up pretty closely in your application.

Did you try non-centered parameterizations of z and w? Not
sure if it'd be better or worse.

We have a unit vector type. If you declare

unit_vector[K] x;

it ensures x' * x = 1 and the default distribution is uniform
over all unit vectors.

The way it does it is as if we'd written it out by hand like this:

parameters {
vector[K] x_raw;

transformed parameters {
vector[K] x = x / dot_self(x);

model {
x_raw ~ normal(0, 1);

This ensures that the distribution on x is uniform over
the unit vectors. See:

Marsaglia, G. (1972). Choosing a point from the surface of a sphere.
The Annals of Mathematical Statistics, 43(2):645–646.

I'm not sure about a fully standardized type. You might be
able to construct it the way we construct the unit vector.
The transform is just:

vector[K] x = (x - mean(x)) / sd(x);

but it's not going to identify on the original y scale. We
really need a function standardize(x) to replace that RHS with
efficient analytic derivatives.

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

[John Hall]

+1 on the standardize function!

In the early iterations I had tried non-centered parameterization for these, but not in the most recent version of the model. I have sufficient data so that there isn't much pooling going on. I calculate a lambda score using a formula in Gelman's book and it is near zero after fitting the model (I noticed that non-centered parameterization helps most when lambda is above 0.5, indicating lots of pooling).

On unit_normal and standardizing: the issue with the way you present them is that the x_raw's don't get identified. As a result, the N_eff's will be near zero for those, which is not a big deal, but it is a big deal for lp__. I really like being able to use lp__ as a diagnostic. I don't typically print all the variables to the screen, just the hyperparamters and lp__. If lp__ is near zero, then at that point I would inspect what's going on with the other variables, If not, then I assume they are fine. When I was using those x_raws, I would get into a habit of ignoring the lp__, so other things could be going wrong in the model and I wouldn't realize it until much later. The only thing I could think of is if the calculation of N_eff for lp__ understands to ignore transformed variables that are not part of the final model.

Also, I was very familiar with unit_vector. My point is that you can't say
unit_vector<lower=0>[N] w;

or some sort of

positive_unit_vector[N] w;

What I did in the model was enforce an approximate unit vector that was always positive (I don't think I can get it exact with this approach because I get imaginary numbers). The reason I included it is that the sign of z and the sign of w weren't identified in any way. Sometimes the model would output w>0 and some z[i] > 0 and sometimes it would output w<0 with that z[i] <0. Forcing w > 0 helped keep it identified in a way that putting priors on the mean of z could not. I got much better results when I included this constraint vs. doing the analysis with just the unit_vector. Indeed, the problem of nailing down the sign seems similar to the issue discussed in the Stan manual under "Inference under Label Switching" in the "Finite Mixtures" section. Unlike finite mixtures, you may not have the information necessary to put a prior that enforces the constraint in a more general PPCA model.

A more general approach would probably to be something like below. This ensures that the largest square value is always positive. I haven't tried it yet and I don't know how well Stan would handle it, but I think it would be more useful than the positive_unit_vector above since it would allow sign identification in all these PPCA models, not just my specific case where I know all the w's are positive.

parameters {

     unit_vector[N] raw_w;
}

transformed parameters {
     unit_vector[N] w;

     w = raw_w;

     if (abs(max(w)) < abs(min(w)))

          w = -w;

}

> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.

> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

--

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/AZIlR7qfqLY/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

[Michael Betancourt]

Actually, positive_unit_vector is _much_ easier to implement 

than unit_vector as the half sphere is isomorphic to R^{n}.  Not

sure if there’s enough use cases to warrant the work to include

it in Stan, though.  Maybe a useful intro project for new devs.

To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.

[John Hall]

Anything that helps identify these models would be helpful.

One thing that I found very helpful was sampling from a simplified model a number of times and plotting pairs of statistics from them.

I'm attaching an R script that simulates from just the factor component of a 1-dimensional PPCA model. I included three charts that illustrate the three main ways that correlations can impact the model:
1) The deterministic relationship between the standard deviation and mean that is caused by taking the norm. This is the only thing I was correcting for in the code above.
2) There is an inverse relationship between the standard deviation of W and Z. This is likely because there is a limit on the variance of WZ and a higher value of the variance of W means a lower value for the variance of Z. The standard deviation of the standard deviation of Z also increases in the inverse of the standard deviation of W.
3) There is a positive relationship between the standard deviation of Z and the standard deviation of the mean of Z. This chart doesn't show it, but there is also a positive relationship between the standard deviation of Z and the mean of Z, but it only becomes apparent if the mean of Z is bigger than zero.

Also, it's amazing how different the charts looks when you change the variables. For instance, if you make the standard on the mu_W_sim equal to 0.1 instead of 1, then there is no correlation between sd_Z_norm and mu_Z_norm. The correlation between sd_W_norm and sd_Z_norm is still there, but much less pronounced. Thus, I suspect the fact that the reason I had better results with the unit vector constrained to be positive was because it kept the standard deviation of W in such a range that it kept these correlations from widening out too much.

> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.

> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

--
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/AZIlR7qfqLY/unsubscribe.

To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

To post to this group, send email to stan-...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Bob Carpenter]

Positive unit vectors can be coded using:

parameters {
simplex[K] y_raw;

transformed parameters {
vector<lower=0>[K] y = sqrt(y_raw);

You have to take the Jacobian into account if you want to put
a distribution on y, but that's trivial for sqrt(), being
just |sqrt'(x)|, so on the log scale in Stan dropping constants
and vectorizing, it's just

target += -log(sum(y));

though you should check my work on that.

- Bob

>> > To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.

>> > To post to this group, send email to stan-...@googlegroups.com.
>> > For more options, visit https://groups.google.com/d/optout.
>>

>> --
>> 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/AZIlR7qfqLY/unsubscribe.

>> To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

[John Hall]

Thanks, I was a little confused about what you were saying because I figured you could put the distribution on y_raw instead of y.

For the square root transformation, I got

  target += -0.5 * log(sum(y));

[Bob Carpenter]

You can do that. I think the 1/2 is a multiplicative
term in the derivative before you take log, so it's
just a constant and drops out.

> On Feb 19, 2017, at 11:53 AM, John Hall <john.mic...@gmail.com> wrote:
>
> Thanks, I was a little confused about what you were saying because I figured you could put the distribution on y_raw instead of y.

Yes, you can do put a distribution on y_raw instead of y.

I'm thinking we should start thinking more about priors
directly on the unconstrained scales.

>
> For the square root transformation, I got
>
> target += -0.5 * log(sum(y));

Ah, you're right---I mistakenly dropped the multiplicative
0.5 alng with the additive one.
the multiplicative one

log sqrt'(x) = log( 0.5 * sqrt(x)^(-0.5) )

= log(0.5) - 0.5 * log(x)

and you can drop the first log(0.5).

- Bob

[John Hall]

I was able to get a one-dimensional PPCA working pretty well in some other constrained circumstances that would be applicable to some of my problems, but never was able to get it work (in the sense of having large N_eff for all variables and lp__) in a more general case that would work with any PPCA. Often, I would get a version working, then multiply the distribution by -1, re-run it, and get gibberish results.

One of my takeaways is the connection of PPCA to finite mixtures: any {w, z} will be equivalent to {-w, -z}. This will cause the mean of z to have a bimodal symmetric distribution a lot of the time. Sometimes also with the mean of w, depending on the assumptions. It can get pretty tricky.

I tend to think of it like: 50% of the time, the mean of w will be greater than or equal to zero, and w ~ N(mu(w+), sd_w), then the mean of z might be positive or negative depending on w, and z ~ N(mu(z(w+)), sd_z). And exactly the opposite the other 50% of the time.

> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.

> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

--
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/AZIlR7qfqLY/unsubscribe.

To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

[Bob Carpenter]

Michael's about to roll out some case studies about
identifying mixtures. These might be tricky.

The draws may be decent, but the n_eff hobbled because
of non-identifiabilities. In some case, you can order
or constrain components.

- Bob