n_eff=1 in RStan 2.5.0 (but not in 2.4.0)

[Matthew Spencer]

We're working on a discrete-time multivariate Gaussian state space model, in which we have time series from a number of different sites, with a random intercept for sites. The process model for time t at site i is

y_{i,t} ~ MVN(mu + a_i + A*y_{i,t-1},Sigma)

where mu is a common intercept, a_i is a random intercept, A is a matrix, and Sigma is a covariance matrix. The observation model for the jth observation at the ith site is

z_j ~ MVN(y_j,Eta)

where Eta is a covariance matrix. There are some time points at which some sites aren't observed (I've skipped over the details of how we match up observations with corresponding true states).

In RStan 2.4.0 on 64-bit Ubuntu, everything seemed to work OK (although I'm sure it's not very efficient code). In RStan 2.5.0, if the warmup is long enough, we get n_eff=1 for all parameters (their values remain constant). Maybe we did something wrong that leads to an overflow or underflow, that in turn causes a gradient problem? And somehow we got away with it in 2.4.0?

R script and stan code attached, along with some example data files (ddatasim.txt, tmatch.txt, yrefs.txt, all needed to run the code) and output of sessionInfo() in R.

[Bob Carpenter]

I don't know much about this kind of model, but we have been
tinkering with the multivariate error checking, so that may
be causing problems. But that's probably not the root cause of
the problem.

For both statistical efficiency (high n_eff per iteration)
and computational efficiency (fast iterations), this model
needs to be rewritten with Cholesky factor parameters. The
problem is that we can't literally do that for your model because
we don't (yet) have (inverse) Wishart densities with Cholesky-factor
parameters.

We'd recommend using a correlation matrix parameterization and
an LKJ prior on correlation matrixes. There's a section in the regression
chapter of the 2.5 manual that explains how to code it.

As written, the y[k] where firstyear[k] == 1 will be problematic because
the model doesn't provide the same support as the parameter
declaration for y, which is unconstrained. This will kill any trajectory
explored by the sampler that moves into a region where the density
is 0 (outside min_y and max_y).

If you really need that uniform density with limited support, you
need to declare a different parameter vector with constraint
<lower=min_y,upper=max_y> for the parameters y[k] for which
firstyear[k] == 1. But it'll probably be easier just to put
a normal density such as

normal((max_y - min_y)/2, (max_y - min_y)/4);

For computational efficiency, the best thing to do would be to use
the vector form of all the densities. It'll probably speed up the
code by a factor of 2 to 10 depending on dimensionality. It'd
be a bit tricky with the switch on firstyear[k] == 1 unless you
define two local vectors (or even better a separate parameter vectors)
for the k such that firstyear[k] == 1; this is a separate motivation
to the constraint).

- 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.
> <ddatasim.txt><coralss7.stan><coralss_short.R><sessionInfo.txt><tmatch.txt><yrefs.txt>

[Matthew Spencer]

Hi Bob

Thanks for those suggestions. The only one I've implemented so far is changing the uniform density for y[k] (where firstyear[k]==1 is true) to a normal density (mean 0, standard deviation 10). The uniform density is often used for this kind of model, but a normal with a high standard deviation is also fine. That did sound as though it might be the issue, but turns out it didn't make much difference. The others will take a bit more effort. I'm sure they will make the code more efficient, but I also think that it may be changes to multivariate error checking that are the main reason why we get n_eff=1 for all parameters in STAN 2.5, but not 2.4. I've attached representative trace plots for four of the parameters, under both versions. These have the same pseudorandom number seeds, and it looks like everything is behaving the same for about the first 1800 iterations, and then suddenly the STAN 2.5 version stops changing.

Matt

[Daniel Lee]

Hi Matt,

About how long does it take you to run this model? I'm trying to run it now and just wondering about how much time it takes to run a single chain for you.



Daniel

[Andrew Gelman]

Hi, these traceplots look wrong.  Shouldn’t each chain be a different color on the plot?

<coralssSTAN2.4.pdf><coralssSTAN2.5.pdf>

[Matthew Spencer]

Hi Andrew

I think the trace plots are right. I only ran one chain in this case. The four traces are four different parameters (the four elements of a 2 x 2 matrix). The warmup period is only shown on the 2.5 version, but it's the same on the 2.4 version (I saved the 2.4 figure using dev.copy2eps, which doesn't support transparency. My postdoc made the 2.5 version on a different computer, and saved it in some other way).

Matt

[Matthew Spencer]

Hi Daniel

Sorry for the slow reply (time difference). For me, in version 2.4, it took about an hour to run 10000 iterations with one chain (on a 3.2 GHz Xeon, Ubuntu 12.04 64-bit). We've been looking for ways to reproduce the problem with a smaller data set, but haven't been able to do so yet.

Matt

[Matthew Spencer]

Traceplots made using

traceplot(stanres,pars="A")

in RStan

[Matthew Spencer]

Updated code attached (using Bob's suggestion of a normal distribution for initial observations in each time series: didn't fix the problem of n_eff=1 in RStan 2.5 on Ubuntu).

I also installed RStan 2.5 on a machine running 64-bit Windows 7, and this code ran fine (output of sessionInfo() on that computer attached). So either it's a problem that doesn't occur with the Windows version of RStan, or it's just a problem with the installation of RStan 2.5 on one particular Ubuntu machine.

[Bob Carpenter]

That's puzzling. It's always possible the other machine has
a memory or bus problem in its hardware. That used to be a lot
more common than it seems to be these days.

Has anyone reproduced the problem on a different machine, either
Linux or Mac?

Daniel found another problem with the model seg-faulting/bus-erroring
in diagnostic mode, which is something else we need to fix.

- Bob

> <coralss7.stan><coralss_short.R><sessionInfo_Windows.txt>

[Andrew Gelman]

Hi, Matt.  I’d recommend running 4 chains because I think this will give a better sense of mixing.  Not that this will solve your problem but it might give you some insight into the convergence problems.

A

[Matthew Spencer]

I think now it must be either a hardware problem or a problem with the installation on that machine. It's reproducible on that machine, but I installed RStan 2.5 on another Ubuntu 12.04 machine, and the same model seemed to run fine.

[Daniel Lee]

Matt,

Thanks for that info. I'm not really sure how to track this down. Can you try a different compiler or even different optimization level while compiling?

Daniel

[Matthew Spencer]

Reinstalled, switched to clang++, and it's working again.

[Daniel Lee]

Thanks for checking. I'm glad it's working. 

[Matthew Spencer]

Looks like I spoke too soon. It worked with that particular simulated data set, but with real data, the same problem recurred. There seem to be some other strange things happening on that computer, so it may well be something that isn't directly connected to STAN at all.

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

[Matthew Spencer]

Now reproduced the same problem with a slightly different model (coralss7mgmt.stan) on another computer. After many iterations, lp__ suddenly becomes -Inf (and so no longer shows on the trace plot). Then n_eff=1 if this happens before the end of the warmup period, and greater than 1 but small if it happens not long after the end). Attached traceplot of lp__: oldchain.pdf. This model contains only normal, multivariate normal and inverse wishart distributions, so the only obvious way this could happen is an overflow. Did it just take a catastrophically bad step? If so, what could I look at to debug it?

[Bob Carpenter]

Does it run into the same problem no matter what the initialization
is?

What happens to the step sizes at points where it fails and what
happens to n_divergent? The n_divergent variable reports on whether
there was an exception or -inf log density that caused the problem.

One thing you can try to do is reduce the initial step size and up
the target acceptance rate. That can help you avoid "catastrphically
bad steps".

We've also relaxed some error checking in some of our multivariate
stuff, which may be causing you problems.

If you sent us (a) the model in question, (b) the seed you're using
to get a bad chain, and (c) the data in a standalone file, we can
try running it in the 2.6 release candidate to see if too-strict error checking
is the problem.

There's also a lot you can do to improve the efficiency of the model as
written, but let's get it working first.

- Bob

> <oldchain.pdf><coralss7mgmt.stan>

[Matthew Spencer]

Plots of lp__, log stepsize, and n_divergent attached. Should it look like that?

The red vertical lines on the bottom two panels are the point where lp__ becomes -Inf. The step size doesn't change at that point. Although n_divergent goes to 1 at that point, there have been lots of 1s before. We get a couple of warning messages:

Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
cholesky_decompose(N4stan5agrad3varE): mm is not symmetric. m[1,2] is -nan:0, but m[2,1] element is -nan:0

Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
stan::prob::inv_wishart_log(N4stan5agrad3varE): LDLT_Factor of random variableunderlying matrix is not positive definite. LDLT_Factor of random variablelast conditional variance is -nan:0.

I was ignoring those, because they only appear at the start.

Will email data and code (and I hope will implement your earlier suggestions for improvements to efficiency when term ends). I'm checking some other initial conditions now.

Thanks

Matt

[Matthew Spencer]

Reproducible from different initial conditions and seed. Similar behaviour in lp__, step size, and n_divergent.

[Matthew Spencer]

Similar behaviour with target acceptance increased to 0.95. We get many fewer iterations for which n_divergent__ is 1, but we still get lp__ of _Inf after enough iterations (not exactly the same number of iterations as with the default target acceptance).

[Matthew Spencer]

Target acceptance increased to 0.99 and initial step size reduced to 1e-5: still essentially the same behaviour. Assuming I'm altering these in the right way:

stanres <- stan(file="../STAN/coralss7mgmt.stan",data=standata,iter=1e4,chains=1,seed=54321,control=list("adapt_delta",0.99,"stepsize",1e-5))

[Rob Trangucci]

Matthew,

Try dialing up the acceptance rate to something like 0.999999 and see if the number of n_divergent__ reduces. I had a hierarchical linear model with some nasty geometry that required the adapt_delta higher than 0.99. Also, try leaving the "stepsize" parameter unset, Stan should optimize that in warmup according to the target acceptance rate in adapt_delta.

Rob

[Bob Carpenter]

This should be an R list, as in:

control=list(adapt_delta=0.99,stepsize =1e-5)

- Bob

[Ben Goodrich]

On Friday, December 19, 2014 2:50:38 PM UTC-5, Matthew Spencer wrote:

Target acceptance increased to 0.99 and initial step size reduced to 1e-5: still essentially the same behaviour. Assuming I'm altering these in the right way:

stanres <- stan(file="../STAN/coralss7mgmt.stan",data=standata,iter=1e4,chains=1,seed=54321,control=list("adapt_delta",0.99,"stepsize",1e-5))

No, it should be

stanres <- stan(file="../STAN/coralss7mgmt.stan",data=standata,iter=1e4,chains=1,seed=54321,control=list(adapt_delta = 0.99, stepsize = 1e-5))

but possibly without the stepsize thing.

Ben

[Matthew Spencer]

Attached: target acceptance increased to 0.999999 (initial step size left at default, seed 54321, random initial conditions). We get many fewer n_divergent, but the same problem occurs. The sample mean of accept_stat__ was 0.9919052.

This plot was based on a simulated data set (coralsim.Rdump, attached: the real data with multivariate normal noise with mean zero added, because I don't own those data), but the real data give essentially the same outcome.

[Rob Trangucci]

Hi Matthew,

Can you send the R code you used to generate the coralsim.Rdump?

Rob

[Matthew Spencer]

The code used to generate coralsim.Rdump isn't likely to be very informative for debugging (all it does is take part of a real data set that I can't distribute, and add some multivariate normal noise to it). Attached anyway so you can see what it was doing (lpdebug5.R).

Probably more useful: a proper simulator that generates data sets with the same structure as a subset of the real data (bivariate time series from each of four coral reefs, each with some missing observations). File coralsim.R, gets the structure of the data set from coralsimexplanatory.Rdata, and uses parameter estimates based on a preliminary run with real data to simulate from a VAR(1) model with Gaussian process error and Gaussian observation error for each reef. Also reproduces the same problem.

It'll certainly be the case that mu (intercept vector), Zeta (covariance matrix for random effect of reef) and b (effect vector for the mgmt explanatory variable, which takes the values 0 or 1) will be very poorly estimated from data sets of this size (four reefs, two with mgmt=1 and two with mgmt=0). Likely to perform much better on the full data set, which contains many more reefs, but runs much more slowly. In any case, poorly identified parameters shouldn't (I think) lead to a proposal with lp_ = -Inf being accepted?

Matt

[Bob Carpenter]

Thanks for sticking with this.

We definitely should not be accepting any proposals
with lp__ = -Inf.

Is that what's happening? The sampler can get stuck for
other reasons having to do with detailed balance preservation --- there
are cases where Euclidean Hamiltonian Monte Carlo just won't be
able to move.

- Bob

> <coralsimexplanatory.Rdata><coralsim.R><lpdebug5.R>

[Matthew Spencer]

lp__ is being reported as jumping to and sticking at -Inf. It's not obvious to me how such a value could occur other than by overflow, but the parameters that give this value don't look dramatically different from those in earlier iterations, where lp__ was on the order of 4 x 10^{-2}, I think.

Calling log_prob from R seg faults, unfortunately.

I wonder whether it's a memory allocation problem? I noticed that the amount of memory being used by R when running this code (as reported by top in Ubuntu) was increasing rapidly, from about 7.5% around 10 minutes after starting the run, to about 77% after almost 2 hours (I didn't see how much memory was being used just before the end of the run, which took a bit less than 2.5h).

[Bob Carpenter]

Matthew --- could you send me:

* Stan program

* R dump file of the data

* seed that causes a problem

and I can try to debug in CmdStan. I tried gathering the
data from this thread and failed. The last rdump

coralsimexplanatory.Rdata

was missing

objects not found: z, min_y, max_y


NUTS shouldn't be able to jump to a -inf value because the Metropolis
accept step should reject a -inf value.

Is the jump always during warmup?

Memory for an iteration only depends on the size of the
expression graph and that memory gets reused, so it shouldn't
grow dramatically over the run.

I'm pretty sure RStan preallocates the space it needs for results,
so that also shouldn't be causing a jump.

I have no idea why calling log_prob() from R would seg fault --- any
seg fault caused by RStan is a bug. If you have a reproducible example,
you (or I) can open an issue on the RStan issue tracker. Unless there's
a reproducible example, there won't be any traction to try and fix it.

- Bob

[Matthew Spencer]

On Tuesday, December 30, 2014 7:02:50 PM UTC, Bob Carpenter wrote:

Matthew --- could you send me:

  * Stan program

  * R dump file of the data

  * seed that causes a problem

and I can try to debug in CmdStan. 

Attached: R dump file (coralpropersim.Rdump) and current version of stan program (coralss7mgmt.stan). I used seed 54321, but the value of the seed doesn't appear to be critical.
 

NUTS shouldn't be able to jump to a -inf value because the Metropolis
accept step should reject a -inf value.  

Yes. That's what's puzzling, and made me wonder whether it's a memory allocation problem.
 

Is the jump always during warmup?

No. Sometimes during warmup, sometimes after warmup.
 

Memory for an iteration only depends on the size of the
expression graph and that memory gets reused, so it shouldn't
grow dramatically over the run.

I'm pretty sure RStan preallocates the space it needs for results,
so that also shouldn't be causing a jump.  

I have no idea why calling log_prob() from R would seg fault --- any
seg fault caused by RStan is a bug.  If you have a reproducible example,
you (or I) can open an issue on the RStan issue tracker.  Unless there's
a reproducible example, there won't be any traction to try and fix it.

I'll first check whether there's something obviously wrong with my R code for calling log_prob(), and that it's reproducible on a different machine with a different OS.

Matt

[Matthew Spencer]

Attached: code that reproducibly seg faults R on call to log_prob(). Tested on 64-bit ubuntu and 64-bit Windows 7, both with R 3.1.2 and rstan 2.5.0. I was intending to create an issue on RStan issue tracker until I saw how complicated it is to attach files to issues in github. Is there an easy way?

[Daniel Lee]

Matthew,

Hopefully this issue was fixed with issue 1200: https://github.com/stan-dev/stan/issues/1200. It won't be available in RStan until the next release, but it has been fixed. If you're wanting to run it asap, you can grab CmdStan and update the Stan repository.

There isn't an easy way to attach files to issues in github, but you can create gists and share them: https://gist.github.com/

Daniel

[Matthew Spencer]

Thanks. I definitely want to try it right away, but what do I have to do? Do I rebuild Stan from stan-develop.zip (from the develop branch on github) and put it where my RStan can see it (/usr/local/lib/R/site-library/rstan/include/stansrc, on my Ubuntu installation)?

[Matthew Spencer]

I managed to build RStan from the develop branch, but I'm still getting the same seg fault problem with log_prob() for this particular model. I guess that I need to update the stan submodule from a different branch on github before building, but how do I do that exactly?

This is what I did to build:

git clone https://github.com/stan-dev/rstan.git rstan.repo
cd rstan.repo
git submodule update --init
git submodule update
cd rstan
sudo make install

Then in R,

library(rstan)

shows me that I now have the version I built:

rstan (Version 2.5.0, packaged: 2015-01-25 16:45:02 UTC, GitRev: 7d6bf44c5b45)

Also, git status shows that I have the develop branch of RStan.

[Matthew Spencer]

Pleased to report that both problems (lp__ getting stuck on -Inf, and log_prob seg faulting) seem to be fixed in 2.6.0. I still have plenty of work to do on making my code run faster, but at least it's running.

[Bob Carpenter]

Glad to hear it and thanks for reporting back.

At least now we won't have to debug these curious bugs.
I still don't know what would have caused them or why
2.6 would've fixed them. Perhaps the stability fixes Rob made
to the multivariate code.

- Bob