big multinomial models

[Ross Boylan]

I added more interactions to my multinomial model and convergence got
slow and poor; I was able to overcome it by running a long time (~16 CPU
hours, although I didn't use all the results since they were too big).
I am looking for feedback on ways to do better, since the model I
estimated is not the most complex model I've contemplated. I have two
particular strategies under consideration; there may well be others.

I somehow got the impression that Bayesian methods would be more
forgiving of models with lots of interactions than frequentist ones,
though I can no longer find any specific references for that! P. 549 of
Gelman et al's Bayesian Data Analysis (2nd edition) says "We fit a
series of models, starting with the simplest, then adding interactions
until we pass the point where the existing data could estimate them
effectively." It's not clear to me what criteria for effectiveness is
intended, though perhaps my experience indicates I have passed that
point with my model and data.

I've been working on a multinomial model with 5 possible categories as
outcomes. Among the covariates are 3 category respondent "race" and 2
category HIV status (probably that will have 3 categories in the end).
I've been using the categorical_logit distribution in the beta code.

My initial model had 29 columns in the model matrix, for 5x29 = 145 beta
coefficients, as well as ~ 5 x 650 = 3250 cluster effects and 2 sigmas.
This behaved badly until I added the sigmas as estimable parameters.

The model that didn't work so great had 50 columns in the model matrix
-> 240 beta parameters.

I have 2 thoughts about how to do better.

First, although fixing the coefficients for one categorical outcome at 0
is not necessary with appropriate priors, would doing so improve
convergence?

Second, perhaps it would be helpful to incorporate more structure into
the different coefficient estimates. Right now the model doesn't make
any allowance for coefficients being related to the same covariate.

Notational digression to clarify what I mean: Two model terms are R*H*X2
and R*H*X2 for R=race, H=HIV status, and X1/2 = covariates. Let
r(i)*h(j)*X2 refer to the coefficient on X2 for level i of R and level j
of h. The model currently thinks there is no more relation between
r(1)*h(1)*X1 and r(1)*h(2)*X1 than between the first and r(3)*h(2)*X2.
In fact, although we want to allow the coefficients for X1 and X2 to
vary be R and H, we generally don't have strong reasons to think they
will (esp for variations by H).

So perhaps it would be better to consider all the coefficients on, e.g.,
X1 as a family with their own mean and sd, and similarly for X2. A
higher level process could model the means and sd's across covariates.

At first blush the second strategy would be a bit involved to implement.

Any advice about whether these or other approaches are worth pursuing?

Thanks.
Ross Boylan

P.S. Yes, I used thinning with the long run; I guess I should have
thinned a bit more.

[Andrew Gelman]

Interactions are good; sometimes it helps to assign them strong priors to pull them toward 0.

> --
> 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/groups/opt_out.

[Bob Carpenter]

On 9/26/13 10:03 AM, Andrew Gelman wrote:
> Interactions are good; sometimes it helps to assign them strong priors to pull them toward 0.
> On Sep 26, 2013, at 1:16 AM, Ross Boylan wrote:
>
>> I added more interactions to my multinomial model and convergence got
>> slow and poor; I was able to overcome it by running a long time (~16 CPU
>> hours, although I didn't use all the results since they were too big).
>> I am looking for feedback on ways to do better, since the model I
>> estimated is not the most complex model I've contemplated. I have two
>> particular strategies under consideration; there may well be others.
>>
>> I somehow got the impression that Bayesian methods would be more
>> forgiving of models with lots of interactions than frequentist ones,
>> though I can no longer find any specific references for that!

They should be in theory because:

-- they don't have the point estimation problem of controlling
hierarchical parameters to not go to variance 0 (though see some
of the work Andrew, Vince Dorie, and Sophia Rabe-Hesketh did on
this issue by using maximal marginal likelihood with boundary-avoiding
priors)

-- they take the uncertainty of their estimates into account for inference.

>> P. 549 of
>> Gelman et al's Bayesian Data Analysis (2nd edition) says "We fit a
>> series of models, starting with the simplest, then adding interactions
>> until we pass the point where the existing data could estimate them
>> effectively." It's not clear to me what criteria for effectiveness is
>> intended, though perhaps my experience indicates I have passed that
>> point with my model and data.

It's not precise because it's a practical issue. You need to be
able to fit the model and then you need to be able to see that the
more complicated model makes a difference to whatever inference you
really care about. You might care about some parameter but not be
able to estimate it accurately due to either software or data constraints.

>> I've been working on a multinomial model with 5 possible categories as
>> outcomes. Among the covariates are 3 category respondent "race" and 2
>> category HIV status (probably that will have 3 categories in the end).
>> I've been using the categorical_logit distribution in the beta code.
>>
>> My initial model had 29 columns in the model matrix, for 5x29 = 145 beta
>> coefficients, as well as ~ 5 x 650 = 3250 cluster effects and 2 sigmas.
>> This behaved badly until I added the sigmas as estimable parameters.
>>
>> The model that didn't work so great had 50 columns in the model matrix
>> -> 240 beta parameters.
>>
>> I have 2 thoughts about how to do better.
>>
>> First, although fixing the coefficients for one categorical outcome at 0
>> is not necessary with appropriate priors, would doing so improve
>> convergence?

Strong enough priors to identify the scales of the
coefficients should do the job.

>> Second, perhaps it would be helpful to incorporate more structure into
>> the different coefficient estimates. Right now the model doesn't make
>> any allowance for coefficients being related to the same covariate.

Pooling can sometimes help, too. Especially if there are lots of
related little groups without much data.

>> Notational digression to clarify what I mean: Two model terms are R*H*X2
>> and R*H*X2 for R=race, H=HIV status, and X1/2 = covariates. Let
>> r(i)*h(j)*X2 refer to the coefficient on X2 for level i of R and level j
>> of h.

You seem to wind up with a multiplicative non-identifiability problem here.
If you cut all the r(i) in half and double all the h(j), you get the
same coefficient. You can usually control this kind of thing by putting
a strong prior on either r or h and letting the other control the overall
scale.

>> The model currently thinks there is no more relation between
>> r(1)*h(1)*X1 and r(1)*h(2)*X1 than between the first and r(3)*h(2)*X2.
>> In fact, although we want to allow the coefficients for X1 and X2 to
>> vary be R and H, we generally don't have strong reasons to think they
>> will (esp for variations by H).
>>
>> So perhaps it would be better to consider all the coefficients on, e.g.,
>> X1 as a family with their own mean and sd, and similarly for X2. A
>> higher level process could model the means and sd's across covariates.
>>
>> At first blush the second strategy would be a bit involved to implement.
>>
>> Any advice about whether these or other approaches are worth pursuing?
>>
>> Thanks.
>> Ross Boylan
>>
>> P.S. Yes, I used thinning with the long run; I guess I should have
>> thinned a bit more.

The command-line version of Stan outputs incrementally rather than storing
everything in memory. We should add that kind of option to RStan, too.

- Bob

[Ross Boylan]

On Thu, Sep 26, 2013 at 01:32:18PM -0400, Bob Carpenter wrote:
>
>
> On 9/26/13 10:03 AM, Andrew Gelman wrote:
>> Interactions are good; sometimes it helps to assign them strong priors to pull them toward 0.
>> On Sep 26, 2013, at 1:16 AM, Ross Boylan wrote:
>>
>>> I added more interactions to my multinomial model and convergence got
>>> slow and poor; I was able to overcome it by running a long time (~16 CPU
>>> hours, although I didn't use all the results since they were too big).
>>> I am looking for feedback on ways to do better, since the model I
>>> estimated is not the most complex model I've contemplated. I have two
>>> particular strategies under consideration; there may well be others.

[snip]

>>> Second, perhaps it would be helpful to incorporate more structure into
>>> the different coefficient estimates. Right now the model doesn't make
>>> any allowance for coefficients being related to the same covariate.
>
> Pooling can sometimes help, too. Especially if there are lots of
> related little groups without much data.
>

By pooling do you mean something different than I was talking about?

>>> Notational digression to clarify what I mean: Two model terms are R*H*X2
>>> and R*H*X2 for R=race, H=HIV status, and X1/2 = covariates. Let
>>> r(i)*h(j)*X2 refer to the coefficient on X2 for level i of R and level j
>>> of h.
>
> You seem to wind up with a multiplicative non-identifiability problem here.
> If you cut all the r(i) in half and double all the h(j), you get the
> same coefficient. You can usually control this kind of thing by putting
> a strong prior on either r or h and letting the other control the overall
> scale.

The notation r(i)*h(j)*X2 refers to a single coefficient (maybe I
should have used a ":" instead of a "*" to fit the usual R model
notation better), which could be called beta[i, j, k] where k indexes
explanatory variables like X1 and X2. I'm currently estimating all
the beta, but not individual terms like r(i).

Now if I take some of the advice about pooling this could become an
issue, e.g., I would estimate one coefficient for the average effect
of X2 by race, r(i)*X2 in the notation I've been using, and then I
would estimate differences around that by group h, where I would
assume, or allow the model to estimate, that such variation was
smaller than the variation between, e.g., r(1)*X2 and r(3)*X1.

To be more concrete about what I was thinking, the model would have
beta[i, k] coefficient for race x covariate, with beta[i, k] ~
normal(0, sd1)
sd1
and
gamma[i, j, k] coefficients for the interaction, distributed normal(0, sd2),
so the total effect for i,j, k = beta[i, k] + gamma[i, j, k]

where I suspect that sd2 < sd1.

Ross

>
>>> The model currently thinks there is no more relation between
>>> r(1)*h(1)*X1 and r(1)*h(2)*X1 than between the first and r(3)*h(2)*X2.
>>> In fact, although we want to allow the coefficients for X1 and X2 to
>>> vary be R and H, we generally don't have strong reasons to think they
>>> will (esp for variations by H).
>>>
>>> So perhaps it would be better to consider all the coefficients on, e.g.,
>>> X1 as a family with their own mean and sd, and similarly for X2. A
>>> higher level process could model the means and sd's across covariates.
>>>
>>> At first blush the second strategy would be a bit involved to implement.
>>>
>>> Any advice about whether these or other approaches are worth pursuing?
>>>
>>> Thanks.
>>> Ross Boylan
>>>
>>> P.S. Yes, I used thinning with the long run; I guess I should have
>>> thinned a bit more.
>
> The command-line version of Stan outputs incrementally rather than storing
> everything in memory. We should add that kind of option to RStan, too.
>
> - Bob
>

[Ross Boylan]

On Thu, Sep 26, 2013 at 03:17:48PM -0400, Bob Carpenter wrote:

>
>
> On 9/26/13 2:12 PM, Ross Boylan wrote:
>> On Thu, Sep 26, 2013 at 01:32:18PM -0400, Bob Carpenter wrote:
>>>
>>>
>>> On 9/26/13 10:03 AM, Andrew Gelman wrote:
>>>> Interactions are good; sometimes it helps to assign them strong priors to pull them toward 0.
>>>> On Sep 26, 2013, at 1:16 AM, Ross Boylan wrote:

> ...

>>>>> Second, perhaps it would be helpful to incorporate more structure into
>>>>> the different coefficient estimates. Right now the model doesn't make
>>>>> any allowance for coefficients being related to the same covariate.
>>>
>>> Pooling can sometimes help, too. Especially if there are lots of
>>> related little groups without much data.
>>>
>> By pooling do you mean something different than I was talking about?
>

> Nope -- same thing.

Good. So I guess I should make the effort to try this.

>
>>>>> Notational digression to clarify what I mean: Two model terms are R*H*X2
>>>>> and R*H*X2 for R=race, H=HIV status, and X1/2 = covariates. Let
>>>>> r(i)*h(j)*X2 refer to the coefficient on X2 for level i of R and level j
>>>>> of h.
>>>
>>> You seem to wind up with a multiplicative non-identifiability problem here.
>>> If you cut all the r(i) in half and double all the h(j), you get the
>>> same coefficient. You can usually control this kind of thing by putting
>>> a strong prior on either r or h and letting the other control the overall
>>> scale.
>>
>> The notation r(i)*h(j)*X2 refers to a single coefficient (maybe I
>> should have used a ":" instead of a "*" to fit the usual R model
>> notation better), which could be called beta[i, j, k] where k indexes
>> explanatory variables like X1 and X2. I'm currently estimating all
>> the beta, but not individual terms like r(i).
>

> Ah --- so that's like what lmer uses for interaction terms in
> their regression notation? I can never figure out R notation,
> but when I saw "term" I figured "*" was multiplication.
>

R models typically look like
y~a+b
meaning outcome y is a function of a and b. If a is continuous that
means estimate one coefficient; if it is categorical, estimate several
(one can pick the exact contrast style). Estimating an interaction is
indicated with a ":" (using spacing to make the binding clearer)
y ~ a + b + a:b
where a:b denotes the interaction of a and b, which is computed
appropriately for the measurement level of each.

The previous formula can be shortened to
y ~ a*b

where "*" means the interaction term and all lower-order combinations
of interactions and main effects.


The model.matrix() function in R creates an appropriately expanded model
matrix so that each column corresponds to one coefficient to be
estimated. I use it to pre-proces the data before I hand it to rstan.

Before I do the model.matrix() I do model.frame() to drop rows with
missing data.

So to get the nesting we were discussing above I will need to figure
out which columns "go together" and pass that info into my stan model.

Ross

[Bob Carpenter]

On 9/26/13 2:12 PM, Ross Boylan wrote:

> On Thu, Sep 26, 2013 at 01:32:18PM -0400, Bob Carpenter wrote:
>>
>>
>> On 9/26/13 10:03 AM, Andrew Gelman wrote:
>>> Interactions are good; sometimes it helps to assign them strong priors to pull them toward 0.
>>> On Sep 26, 2013, at 1:16 AM, Ross Boylan wrote:

...

>>>> Second, perhaps it would be helpful to incorporate more structure into
>>>> the different coefficient estimates. Right now the model doesn't make
>>>> any allowance for coefficients being related to the same covariate.
>>
>> Pooling can sometimes help, too. Especially if there are lots of
>> related little groups without much data.
>>
> By pooling do you mean something different than I was talking about?

Nope -- same thing.

>>>> Notational digression to clarify what I mean: Two model terms are R*H*X2
>>>> and R*H*X2 for R=race, H=HIV status, and X1/2 = covariates. Let
>>>> r(i)*h(j)*X2 refer to the coefficient on X2 for level i of R and level j
>>>> of h.
>>
>> You seem to wind up with a multiplicative non-identifiability problem here.
>> If you cut all the r(i) in half and double all the h(j), you get the
>> same coefficient. You can usually control this kind of thing by putting
>> a strong prior on either r or h and letting the other control the overall
>> scale.
>
> The notation r(i)*h(j)*X2 refers to a single coefficient (maybe I
> should have used a ":" instead of a "*" to fit the usual R model
> notation better), which could be called beta[i, j, k] where k indexes
> explanatory variables like X1 and X2. I'm currently estimating all
> the beta, but not individual terms like r(i).

Ah --- so that's like what lmer uses for interaction terms in
their regression notation? I can never figure out R notation,
but when I saw "term" I figured "*" was multiplication.

- Bob

[Alexander Franks]

I've been having some trouble fitting multinomial logit in stan and have read a few things about issues with categorical(softmax()) here that suggested there may be a memory leak?  

I have a very large multinomial model that I am trying to fit: about 50 categories and a couple thousand parameters.  The model matrix is very sparse since I include "person-specific" random effect which each observation (with partial pooling on the params).  Each row of the model matrix has roughly 30 non-zero entries (out of a few thousand).  I have about 75,000 data points.  

When I tried to fit this model in stan everything compiled nicely but it hung once it began sampling.  Eventually it segfaulted.

I also have a few questions about implementing this efficiently in stan.  With a very sparse model matrix, am I better off specifying response[n] ~ categorical(softmax(beta*x[n]) as in the documentation, or rather should I specify eta[n] in the transformed parameters section using a series of for loops to index and sum up the relevant parameters?   E.g. in a loop something like

eta[n,k] <- eta[n,k]+beta[index[n,k]]

Finally, is it better at this stage to use softmax or to specify probs[n] <- exp(eta[n])/sum(eta[n]) in transformed parameters section?  categorical_logit will be available in stan 2.0?

Thanks,

Alex

[Bob Carpenter]

On 9/27/13 3:43 PM, Alexander Franks wrote:
> I've been having some trouble fitting multinomial logit in stan and have read a few things about issues with
> categorical(softmax()) here that suggested there may be a memory leak?

As far as we know, there's no memory leak for softmax
in 1.3 or in the coming 2.0 versions. Having said that,
it is memory intensive (more about that below).

> I have a very large multinomial model that I am trying to fit: about 50 categories and a couple thousand parameters.
> The model matrix is very sparse since I include "person-specific" random effect which each observation (with partial
> pooling on the params). Each row of the model matrix has roughly 30 non-zero entries (out of a few thousand). I have
> about 75,000 data points.
>
> When I tried to fit this model in stan everything compiled nicely but it hung once it began sampling. Eventually it
> segfaulted.

Segfaulted or just ran out of memory with a bad_alloc
exception? If it really was a segfault, could you send
us the output and some description of the platform you're
running on?

> I also have a few questions about implementing this efficiently in stan. With a very sparse model matrix, am I better
> off specifying response[n] ~ categorical(softmax(beta*x[n]) as in the documentation, or rather should I specify eta[n]
> in the transformed parameters section using a series of for loops to index and sum up the relevant parameters? E.g. in
> a loop something like
> eta[n,k] <- eta[n,k]+beta[index[n,k]]

I'm not exactly sure what the two alternatives are you're
considering or what the index[n,k] variable is doing. Could
you share the whole loop?

In general, you don't want to be multiplying by known
0 values if you can help it.

But rather than doing a loop like this:

real foo;
foo <- 0;
for (n in 1:N)
foo <- foo + bar[n];

it's better to do

foo <- sum(bar);

This is true even if everything's not in a simple array or vector
like bar to begin with (that is, if you have to allocate and
populate bar). The reason is because the sum operation uses
less memory than repeated application of + and performs fewer
(virtual) function applications for computing derivatives. The
allocation and loop for the array are relatively cheap.

> Finally, is it better at this stage to use softmax or to specify probs[n] <- exp(eta[n])/sum(eta[n]) in transformed
> parameters section? categorical_logit will be available in stan 2.0?

Use softmax. It should be have exactly the same way
in 1.3, which just defines softmax that way.

In Stan 2 (which we hope to launch a beta for next week),
has custom derivatives for softmax, which will be much faster
and much less memory intense. The issue with memory is that
Stan stores the whole expression tree for the log probability,
which can get big. For instance, softmax does a lot of work
for each n because it needs to be normalized.

The categorical_logit that's part of 2.0 gets around that
problem AND has a custom derivative. So it's much faster
and much less memory intensive.

We would've done this sooner if we'd known there were
so many multi-logit fitters out there. We're trying to
optimize according to perceived demand.

- Bob

[Ross Boylan]

On Fri, Sep 27, 2013 at 05:04:04PM -0400, Bob Carpenter wrote:
>
>
> On 9/27/13 3:43 PM, Alexander Franks wrote:
>> I've been having some trouble fitting multinomial logit in stan and have read a few things about issues with
>> categorical(softmax()) here that suggested there may be a memory leak?
>
> As far as we know, there's no memory leak for softmax
> in 1.3 or in the coming 2.0 versions. Having said that,
> it is memory intensive (more about that below).
>

I thought the problem I encountered was a memory leak (rstan development code).

Memory use grew steadily as the job ran, whereas I thought allocation all happened at the beginning.
And the memory use seemed out of whack with its expected value.
The memory use was not in R, as evidenced by gc().
Finally, the memory use did not go back down after deleting the return objects.

The problem went away when I stopped using categorical(softmax()).

Ross

[Alexander Franks]

Not much output to report:

TRANSLATING MODEL 'mlogit' FROM Stan CODE TO C++ CODE NOW.

 COMPILING THE C++ CODE FOR MODEL 'mlogit' NOW.

 SAMPLING FOR MODEL 'mlogit' NOW (CHAIN 1).

 Error: segfault from C stack overflow

 > 

The segfault comes after running for approx 15 min.  I'm running OSX 10.8.5.

Also, I can only get it to run with init=0.  When I use the default I get this error message:

Error : Error in function stan::prob::categorical_log(N4stan5agrad3varE): Probabilities parameter is not a valid simplex. The element at 0 is nan:0, but should be greater than or equal to 0 error occurred during calling the sampler; sampling not done

I don't see any reason why I would get any nans.

By the way, here is the loop where I specify eta[n]. It's pretty messy, but all I'm doing is specifying the logit probabilities as a sum of 6 terms.  O_id and D_id specify the index into alpha and beta (which are each of dimension 300x10)

transformed parameters {

  vector[K*5+1] eta[N];

  for(n in 1:N) {

        for(o in 1:5) {

              for(k in 1:K) {

                    eta[n,o*k] <- alpha[O_id[n,o],k];

                    for(d in 1:5)

                          eta[n,o*k] <- eta[n,o*k]+beta[D_id[n,d],k]*Frac[n,5*(d-1)+o];

                     eta[n,o*k] <- eta[n,o*k];

              }

        }

        eta[n,(K*5+1)] <- 0;

  }

}

[Ross Boylan]

I estimated a model in which the prior for the 3-way interactions were
allowed to have a different std dev than the other model terms (the
previous version treated all beta as being drawn from the same
distribution). I allowed stan to fit what this standard deviation was
using a hyperprior.

This worked much better. It took about 45 minutes with 3 processors and
produced a minimum effective sample size of 115; the first quartile was
over 1,200.

It may also be relevant that I changed the coding scheme to contr.sum so
that all the coefficients summed to zero (e.g., for the A*B interaction
the sum of the coefficients for each combination of levels was zero).

The prior for the lower order beta had a median sd of 0.08 and for the
3-way terms an sd of 0.002!

So, all the 3 way interactions are trivial (the posteriors are what
you'd expect from such a tight prior).

This is another example of estimation improving by allowing an sd to be
fit, in this case as an extremely small value. Conversely, another
example of convergence problems with excessively diffuse priors.

I worry a little that if one of the 3-way interactions really were
important it might be dragged down by their general smallness. I
suppose the model would not end up with such a tight prior if that were
the case.

Ross

[Bob Carpenter]

On 9/27/13 6:03 PM, Alexander Franks wrote:
> Not much output to report:
>
> TRANSLATING MODEL 'mlogit' FROM Stan CODE TO C++ CODE NOW.
> COMPILING THE C++ CODE FOR MODEL 'mlogit' NOW.
> SAMPLING FOR MODEL 'mlogit' NOW (CHAIN 1).
> Error: segfault from C stack overflow

That's really odd because there's nothing deeply
recursive in Stan.

A quick search on the web seems like it may be related to
R running out of memory.

What can happen if the model's complicated is that you'll
need more memory for certain iterations than others if
more leapfrog steps are necessary.

> The segfault comes after running for approx 15 min. I'm running OSX 10.8.5.

Can you run fewer iterations and get it to complete?

> Also, I can only get it to run with init=0. When I use the default I get this error message:
>
> Error : Error in function stan::prob::categorical_log(N4stan5agrad3varE): Probabilities parameter is not a valid
> simplex. The element at 0 is nan:0, but should be greater than or equal to 0 error occurred during calling the sampler;
> sampling not done

That may signal a problem with how you're initializing
values. It can also be a problem if you're doing a large
regression and a random init takes the likelihood into
NaN.

> I don't see any reason why I would get any nans.

Are you providing it with a simplex or with the result
of softmax?

Unexpected NaN usually occur due to overflows or underflows.

Stan 2 should be more stable around categorical distributions
because I added a log_sum_exp-like implementation.

> By the way, here is the loop where I specify eta[n]. It's pretty messy, but all I'm doing is specifying the logit
> probabilities as a sum of 6 terms. O_id and D_id specify the index into alpha and beta (which are each of dimension 300x10)
>
> transformed parameters {
> vector[K*5+1] eta[N];
> for(n in 1:N) {
> for(o in 1:5) {
> for(k in 1:K) {
> eta[n,o*k] <- alpha[O_id[n,o],k];
> for(d in 1:5)
> eta[n,o*k] <- eta[n,o*k]+beta[D_id[n,d],k]*Frac[n,5*(d-1)+o];
> eta[n,o*k] <- eta[n,o*k];
> }
> }
> eta[n,(K*5+1)] <- 0;
> }
> }

Is eta what gets fed into categorical as a parameter?
If so, how do you normalize?

Also, it looks like this line after the for (d ...) loop was
a mistake:

> eta[n,o*k] <- eta[n,o*k];

- Bob

>
>
> On Friday, September 27, 2013 5:45:41 PM UTC-4, Ross Boylan wrote:
>
> On Fri, Sep 27, 2013 at 05:04:04PM -0400, Bob Carpenter wrote:
> >
> >
> > On 9/27/13 3:43 PM, Alexander Franks wrote:
> >> I've been having some trouble fitting multinomial logit in stan and have read a few things about issues with
> >> categorical(softmax()) here that suggested there may be a memory leak?
> >
> > As far as we know, there's no memory leak for softmax
> > in 1.3 or in the coming 2.0 versions. Having said that,
> > it is memory intensive (more about that below).
> >
>
> I thought the problem I encountered was a memory leak (rstan development code).
>
> Memory use grew steadily as the job ran, whereas I thought allocation all happened at the beginning.
> And the memory use seemed out of whack with its expected value.
> The memory use was not in R, as evidenced by gc().
> Finally, the memory use did not go back down after deleting the return objects.
>
> The problem went away when I stopped using categorical(softmax()).
>
> Ross
>

[Bob Carpenter]

Is there a reproducible example? (If you already sent
it to me, sorry for asking again, but can you give me
a hint as to when or send again.)

The memory will jump every chain when running in R.

The memory can also jump if you need a lot of leapfrog
iterations, but that shouldn't jump enough to cause
an out-of-memory error.

- Bob

[Ross Boylan]

On Fri, Sep 27, 2013 at 08:25:11PM -0400, Bob Carpenter wrote:
> Is there a reproducible example? (If you already sent
> it to me, sorry for asking again, but can you give me
> a hint as to when or send again.)

It always failed for me :)
The dataset has confidential info; I can investigate behavior with fake data.

Ross

>
> The memory will jump every chain when running in R.
>
> The memory can also jump if you need a lot of leapfrog
> iterations, but that shouldn't jump enough to cause
> an out-of-memory error.
>
> - Bob
>
> On 9/27/13 5:45 PM, Ross Boylan wrote:
>> On Fri, Sep 27, 2013 at 05:04:04PM -0400, Bob Carpenter wrote:
>>>
>>>
>>> On 9/27/13 3:43 PM, Alexander Franks wrote:
>>>> I've been having some trouble fitting multinomial logit in stan and have read a few things about issues with
>>>> categorical(softmax()) here that suggested there may be a memory leak?
>>>
>>> As far as we know, there's no memory leak for softmax
>>> in 1.3 or in the coming 2.0 versions. Having said that,
>>> it is memory intensive (more about that below).
>>>
>>
>> I thought the problem I encountered was a memory leak (rstan development code).
>>
>> Memory use grew steadily as the job ran, whereas I thought allocation all happened at the beginning.
>> And the memory use seemed out of whack with its expected value.
>> The memory use was not in R, as evidenced by gc().
>> Finally, the memory use did not go back down after deleting the return objects.
>>
>> The problem went away when I stopped using categorical(softmax()).
>>
>> Ross
>>
>

[Bob Carpenter]

On 9/27/13 8:35 PM, Ross Boylan wrote:
> On Fri, Sep 27, 2013 at 08:25:11PM -0400, Bob Carpenter wrote:
>> Is there a reproducible example? (If you already sent
>> it to me, sorry for asking again, but can you give me
>> a hint as to when or send again.)
>
> It always failed for me :)
> The dataset has confidential info; I can investigate behavior with fake data.

That'd be very helpful. Pretty much all of
softmax and categorical is changed for Stan 2.

We should run it through a memory checker before
release, though, just to be sure. (The problem there
is we get swamped with false positives from the way
our auto-diff manages memory, which is where any
leak would be.)

- Bob