faster adaptation

[Bob Carpenter]

[retitled and moved to stan-dev]

I'd suggest starting with a Hessian, but then we're back
to O(N^2) time even if it's diagonal. But it's very easy
to parallelize. Can we look at the gradients and get something
out of that (I know you dismissed those approaches in the SoftAbs
paper)?

What I don't understand is why it's so slow to adapt
the dimensions, even if it's not as bad as 10^6. Is there
some theory saying why it'd be so slow or some insight? I'd
think we'd start diagonal, find a very small step size is
necessary, run NUTS and see ourselves hitting the tree depth
limit straight off in the small-scale dimensions. Would reducing
tree depth during warmup help for an initial attempt to determine
scales? I know there's a chicken-and-egg problem with not
being in the typical set to start, so we don't want to adapt
too agressively out in the tails (or at the mode).

- Bob

> On Aug 10, 2016, at 6:49 PM, Michael Betancourt <betan...@gmail.com> wrote:
>
> The fix was rescaling the parameters by 10^{6}, which
> implied that the original model had a condition number
> of at least 10^{6} which is going to be hard for _any_
> adaptive algorithm.
>
> On Aug 10, 2016, at 4:15 AM, Bob Carpenter <ca...@alias-i.com> wrote:
>
>> Thanks for reporting back. We really need to improve our
>> adaptation because we should be able to learn these marginal
>> scales during warmup.
>>
>> We generally recommend 50% of the iterations be devoted to
>> warmup. At most that'll cost 50% efficiency, but it has
>> a higher chance of adaptation converging, especially in
>> ill conditioned problems.
>>
>> - Bob
>>
>>> On Aug 9, 2016, at 4:18 PM, Jeroen CORSTJENS <jeroen.c...@uhasselt.be> wrote:
>>>
>>>
>>> Thanks for the suggestions. I tried rescaling by standardizing the numerical predictor variables (dividing by their standard deviation as they were already centred) and this made the sampling go a lot faster. It took almost 27 hours to run 4 chains in parallel with 1000 warm-up iterations and 3000 post-warmup iterations. However, the model did not converge, all group-level effects had Rhats larger than 1.1. Specifying tighter priors gave convergency problems as well (also when specifying a prior for all population-level effects).
>>>
>>> I did discover an alternative that runs faster and converges perfectly. Simply multiplying the 1/Y values with a large constant (1 000 000 in my case) did the trick. This time, it took nearly 12 hours to complete and gave no convergency problems.
>>>
>>>
>>> --
>>> 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.
>>
>> --
>> 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.
>
> --
> 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.
>
>

[Michael Betancourt]

> I'd suggest starting with a Hessian, but then we're back
> to O(N^2) time even if it's diagonal. But it's very easy
> to parallelize. Can we look at the gradients and get something
> out of that (I know you dismissed those approaches in the SoftAbs
> paper)?

No — absent of the full Hessian, or a finite difference approximation
thereof, there isn’t much that can be done to infer the covariance
structure accurately.

> What I don't understand is why it's so slow to adapt
> the dimensions, even if it's not as bad as 10^6. Is there
> some theory saying why it'd be so slow or some insight? I'd
> think we'd start diagonal, find a very small step size is
> necessary, run NUTS and see ourselves hitting the tree depth
> limit straight off in the small-scale dimensions. Would reducing
> tree depth during warmup help for an initial attempt to determine
> scales? I know there's a chicken-and-egg problem with not
> being in the typical set to start, so we don't want to adapt
> too agressively out in the tails (or at the mode).

The argument is as follows — the longest length scale determines
the integration time and the shortest length scale determines the
step size, so the number of leapfrog steps goes as longest length
scale / shortest length scale. A well-chosen mass matrix can
compensate for this, but initially we don’t know the ideal mass matrix
so we have to take that computational hit initially.

Because the variance quantifies the breadth of the posterior, we
need complete trajectories to get accurate estimates of the variance.
Prematurely terminating the trajectories gives stunted trajectories
and a chain with large autocorrelations, which will in general
underestimate variances and make adaptation even slower.

The only real possibility for improvement is using Hessian information,
although some care has to be taken there if the curvature varies
strongly with position, or sharing information across chains, or even
both at the same time.

[Bob Carpenter]

Thanks. This is super helpful for me.

Couldn't we look at the diagonal of the Hessian as a guess
of the mass matrix? That we can do in parallel, at least
in theory. And we can do it without much memory overhead,
even if we don't do it in parallel.

I thought early termination would be a problem. I knew
it wouldn't provide the proper covariance, but I thought it
might get moving in the right direction faster and help
the chicken-and-egg nature of the problem in less time.

Varying curvature's a problem no matter what we do with EHMC!
I'm not even sure what the "right" mass matrix is for the
funnel in EHMC. One's optimal for the neck and one for the
mouth, but I suppose we need the one closer to the neck and
just have to suck it up with lots of leapfrog steps for
the flat part.

- Bob

> --
> You received this message because you are subscribed to the Google Groups "stan development mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+u...@googlegroups.com.

[Michael Betancourt]

The problem with using the Hessian naively is that it might improve things locally while making it worse everywhere else. What we really want is an average of the Hessian across parameter space which brings us back to the global covariance which we're estimating empirically already. Remember, there's only ever one option which is an expectation with respect to the posterior!

[Krzysztof Sakrejda]

I know the implementation might be a headache but it seems like one place where using multiple chains to estimate one mass matrix might make sense. Maybe even just using the average from multiple chains would help. Krzysztof

[Bob Carpenter]

First, I just meant using the Hessian locally to help
us bootstrap our way to the typical set. It seems like
it could help with conditioning, even out in the tails.
I also understand that with parallelization we can explore
more parts of the space at initialization, but we still won't
be able to expect random inits to come anywhere near the
typical set due to the combinatorics.

> On Aug 11, 2016, at 12:38 AM, Michael Betancourt <betan...@gmail.com> wrote:
>
> The problem with using the Hessian naively is that it might improve things locally while making it worse everywhere else.

Exactly, but what we need to do is get from some random init
to the typical set. I'd think the local Hessian would make
sense for that use case.

> What we really want is an average of the Hessian across parameter space which brings us back to the global covariance which we're estimating empirically already.

If our goal is to reduce divergences to zero, why would
the posterior expectation of the Hessian be the best.
I'd have thought something that worked in the highly curved
parts of the space would be best for eliminating divergences,
even if it's not optimal for exploring the less constrained
volumes of the posterior.

> Remember, there's only ever one option which is an expectation with respect to the posterior!

Don't worry, that's etched into my brain permanently. I just
don't see why minimizing squared error over the posterior is
ideal for avoiding divergences.

- Bob

[Michael Betancourt]

You have to be careful — because HMC trajectories can be so
long the Hessian at one end can be a really, really bad approximation
of the Hessian at the other end! This is why full Riemannian HMC
is so important and other adaptive approaches tend to work very poorly.

The problem with divergences isn’t the Hessian so much but the minimal
length scale, which can be controlled with a threshold a la SoftAbs. The
Hessian also has the rotation which can also screw up trajectories and
global stability of the integrator.