"divergent gradient" when log prob and its gradient are both very finite

[Kevin Van Horn]

 I have some modeling code that I've been testing out, and in one of my test cases with synthetic data I get the error "Rejecting user-specified initialization because of divergent gradient."

This is the line of code that produces the error:

    fit <- rstan::sampling(sm, data=data, chains=1, warmup=1, iter=1, init=list(params), test_grad=TRUE)

When I subsequently call

    val <- rstan::log_prob(fit, rstan::unconstrain_pars(fit, params))

the value I get is reasonable (around -20000). When I repeatedly call rstan::log_prob to compute a finite-differences approximation to the gradient, using a delta of +/- 1e-4, all components of the gradient has absolute values less than 300.

This has me utterly confused. Can you get the "divergent gradient" error for any reason other than a very large gradient value?

[Kevin Van Horn]

I have some more information after digging into this some more. The problem occurs only if the following user-defined Stan function is called with a relatively large value for c:

    real ll_censored_using_psevere(
      int c, int has_severe, real alpha, real nb_rate, real logit_psevere)
    {
      real nbrate1;
      real psevere;
      real pnotsevere;
      real llnone;
      psevere <- inv_logit(logit_psevere);
      pnotsevere <- inv_logit(-logit_psevere);
      nbrate1 <- (nb_rate + psevere) / pnotsevere;
      llnone <- neg_binomial_ccdf_log(c - 1, alpha, nbrate1) -
                alpha * log1p(psevere / nb_rate);
      if (has_severe) {
        return log_diff_exp(neg_binomial_ccdf_log(c - 1, alpha, nb_rate),
                                  llnone);
      }
      else {
        return llnone;
      }
    }

In the original data set that exhibited the problem, c ranged between 29 and 119, while alpha, nb_rate, and logit_psevere were all moderate values resulting from 0-initialization of model parameters. When I modified the data set so that c=29 always when this function is called, with no other change, I no longer got the "divergent gradient" error.

[Michael Betancourt]

This is a problem with the negative binomial CDF implementation.

Essentially the gradient requires some crazy hypergeometric

function derivatives that are all implemented with power series.

The initial implementation was not robust to large values and

would never converge.  We added a max iteration threshold,

but then the series just terminate prematurely and the resulting

gradients are inaccurate and cause the integrator to diverge.

Coincidentally enough, I encountered the same problem in

a model on which I’m working and spent way too much work

to fix it.  It’s currently under a pull request at 

https://github.com/stan-dev/stan/pull/1433.  It may still be

unstable for super pathological arguments, but it should be

much more robust.

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

[Kevin Van Horn]

Thanks, Michael. Do you think it's worth me filing an issue, or do you think your existing issue 1426 covers it?

[Michael Betancourt]

I would wait for this pull request to be merged and see if it fixes your problem.  If it doesn't then create an issue.

[Kevin Van Horn]

When you say "not robust to large values," do you mean large values of the count, large values of alpha and/or beta, or both? Over what range is the current implementation stable (prior to your fix)?

On Wednesday, April 29, 2015 at 10:11:23 PM UTC-6, Kevin Van Horn wrote:

[Michael Betancourt]

> When you say "not robust to large values," do you mean large values of the count, large values of alpha and/or beta, or both?

All three. Formally problems arise in the current implementation
of the incomplete beta function when alpha is big and/or
when z is close to 1. Translating to the neg_binomial_2
parameterization this means large alpha, large beta, and/or
large counts.

> Over what range is the current implementation stable (prior to your fix)?

It’s not trivial to specify a stable region in the current implementations
as the convergence of the underlying power series depends on the
interactions of all of the parameters. Problems can arise in alpha
values as low as 100 if z is sufficiently close to 1 (say z > 0.99) and
it might be even worse.

The fix was checked on a bunch of values pushing into the millions
for z close to zero and one and seemed to be sufficiently accurate
everywhere.