2nd & 3rd derivative auto-dif for diagonal RM-HMC

[Bob Carpenter]

We discussed using auto-dif at the meeting today
for second- and third-order derivatives. The
chain rule's a bit messier and also involves the
lower-order derivatives:

http://en.wikipedia.org/wiki/Chain_rule#Higher_derivatives

The derivative propagation can still be done in one
efficient backward pass one edge of the graph at a time.

So we'd need to extend the auto-dif vars to have
values, first derivatives, second derivatives, and
third derivatives. So it'll take about twice as much
space. One approach would be to write chain1(), chain2()
and chain3() derivative propagation rules, running each
in turn (each is required for the next). These can then
be consolidated into chain() which just calls them in order.

The nice part is that it won't mess with vectorization
in our probabilty functions.

The not-so-nice part is that it requires new auto-dif
variables (easy) and a whole lot of chainN() method writing
(a whole lot of slogging).

- Bob

[Matt Hoffman]

I think we can probably get by without the 3rd derivatives using an
appropriate finite differencing scheme (or something equivalent but
more numerically stable), same as you can do a Hessian-vector product
without computing the Hessian explicitly.

I'd actually say that being able to compute Hessian-vector products
using Pearlmutter's algorithm or whatever is more valuable than using
the higher-order chain rule to compute the full Hessian, since you can
get the Hessian easily via a series of Hv products.

Also, if we had efficient and stable Hv products, then we could do
some really cool Hessian-free optimization to get point estimates
efficiently.

Also, local curvature information isn't always a reliable guide—the
Laplacian distribution has no second derivative, for example. So as
much as I'd love to have the option of using higher-order derivatives
in Stan, I'm not sure it should be as high-priority as, say,
vectorization. (Although it might be more fun to implement.)

Something else that might be cool and relatively easy to implement
would be a user-defined "metric" block for RMHMC. If undefined, we'd
just do what we're doing. But if the block were defined, then we'd do
RM-NUTS with the metric the user specifies. An expert user (admittedly
there aren't many yet!) would then be able to play with different
Riemannian metrics that might dramatically improve their sampling
efficiency. (And Mike and I would be able to more easily test out
research ideas.)

Just some thoughts.

> --
>
>

[Andrew Gelman]

M3,
We definitely want to vectorize, but Mike was claiming that his RM-NUTS algorithm would be able to adapt step sizes and so better fit various models. Regarding 3rd derivs etc., I will defer to the rest of you. I just remember that _first_ derivs work better analytically than numerically.
A

> --
>

[Michael Betancourt]

Matt,

We were talking about a diagonal Hessian, or diagonal outer-product
approximation; something simple just to try and adjust to the local
step size.

I'm all for efficient Hessian-vector products and abstracting the code
to admit various metrics, as well. I imagine that we're going to have
to tackle that problem eventually.

> --
>
>

[Matt Hoffman]

Yup, I've been giving a lot of thought to such methods over the last
couple of months. I actually wrote most of a paper about one for NIPS,
before realizing that it didn't work at all for target distributions
that are Laplacian (zero curvature everywhere) or leptokurtotic (e.g.,
exp{-θ^4}, which has zero curvature around its mode) and might not
work for heavy-tailed distributions (since any continuously
differentiable heavy-tailed distribution must by definition include a
region with zero curvature). The issue is that the optimal choice for
a normal target distribution (step size proportional to 1/sqrt(largest
eigenvalue of Hessian)) will blow up to infinity at some point for
almost any target distribution that departs too far from normality.
One can introduce a maximum step size (as in Mike's "B-S" metric), but
then you've got to fit that somehow.

All of which is just to say that I'm becoming a little skeptical of
the existence of a one-size-fits-all solution to this problem. Not
that I've given up hope—I've got some ideas that I'm reasonably
excited about—but I'm not as convinced as I once was that local
curvature information offers a universal solution to the step size
problem.

Matt

> --
>
>

[Marcus Brubaker]

On Tue, Sep 18, 2012 at 5:36 PM, Matt Hoffman <mdho...@cs.princeton.edu> wrote:
> Yup, I've been giving a lot of thought to such methods over the last
> couple of months. I actually wrote most of a paper about one for NIPS,
> before realizing that it didn't work at all for target distributions
> that are Laplacian (zero curvature everywhere) or leptokurtotic (e.g.,
> exp{-θ^4}, which has zero curvature around its mode) and might not
> work for heavy-tailed distributions (since any continuously
> differentiable heavy-tailed distribution must by definition include a
> region with zero curvature). The issue is that the optimal choice for
> a normal target distribution (step size proportional to 1/sqrt(largest
> eigenvalue of Hessian)) will blow up to infinity at some point for
> almost any target distribution that departs too far from normality.
> One can introduce a maximum step size (as in Mike's "B-S" metric), but
> then you've got to fit that somehow.
>
> All of which is just to say that I'm becoming a little skeptical of
> the existence of a one-size-fits-all solution to this problem. Not
> that I've given up hope—I've got some ideas that I'm reasonably
> excited about—but I'm not as convinced as I once was that local
> curvature information offers a universal solution to the step size
> problem.

Isn't this precisely the reason that RMHMC proposes the Fisher
information matrix as the metric rather than the Hessian of the log
posterior? Even for problems where the Hessian is not positive
definite, the Fisher information matrix still is. It may be
practically impossible to compute (unlike the Hessian) but it should
be well behaved in theory.

Cheers,
Marcus

[Michael Betancourt]

The Fisher information is indeed guaranteed to be positive definite
but it's somewhat artificial to HMC. Firstly it's really a metric on
an entirely different manifold (where each point is a different
likelihood distribution) than the manifolds we use for HMC (where each
point is just a point in the support of a single distribution).
Secondly, it requires a likelihood-prior separation so it can't be
used on arbitrary distributions (perhaps there's a general
augmentation/integration scheme to be found?). Then, of course, is
the computational problems.

In most cases the Hessian looks like it should be the right metric,
but everything falls apart of the potential isn't convex everywhere
across its support. Which begs the question, what properties of a
distribution cause its potential to be concave? And what does that
inform about possible metrics.

The Laplacian distribution is an interesting case in that it's almost
--too-- easy to sample from. Instead of trying to fit a
regularization coefficient one could always just default to a constant
mass matrix when the Hessian vanishes identically. I think there's
some promise in metric that decompose into a sum of other metrics, but
like most of the ideas floating around there's nothing yet concrete.

On Tue, Sep 18, 2012 at 5:43 PM, Marcus Brubaker
<marcus.a...@gmail.com> wrote:
> On Tue, Sep 18, 2012 at 5:36 PM, Matt Hoffman <mdho...@cs.princeton.edu> wrote:
>> Yup, I've been giving a lot of thought to such methods over the last
>> couple of months. I actually wrote most of a paper about one for NIPS,
>> before realizing that it didn't work at all for target distributions
>> that are Laplacian (zero curvature everywhere) or leptokurtotic (e.g.,
>> exp{-θ^4}, which has zero curvature around its mode) and might not
>> work for heavy-tailed distributions (since any continuously
>> differentiable heavy-tailed distribution must by definition include a
>> region with zero curvature). The issue is that the optimal choice for
>> a normal target distribution (step size proportional to 1/sqrt(largest
>> eigenvalue of Hessian)) will blow up to infinity at some point for
>> almost any target distribution that departs too far from normality.
>> One can introduce a maximum step size (as in Mike's "B-S" metric), but
>> then you've got to fit that somehow.
>>
>> All of which is just to say that I'm becoming a little skeptical of
>> the existence of a one-size-fits-all solution to this problem. Not

>> that I've given up hope--I've got some ideas that I'm reasonably
>> excited about--but I'm not as convinced as I once was that local

>> curvature information offers a universal solution to the step size
>> problem.
>
> Isn't this precisely the reason that RMHMC proposes the Fisher
> information matrix as the metric rather than the Hessian of the log
> posterior? Even for problems where the Hessian is not positive
> definite, the Fisher information matrix still is. It may be
> practically impossible to compute (unlike the Hessian) but it should
> be well behaved in theory.
>
> Cheers,
> Marcus
>

> --
>
>

[Matt Hoffman]

Just to clarify Mike's point, the metric used by G&C in their RMHMC
paper is the sum of the Fisher matrix of p(y|θ) and the Hessian of
-log p(θ). It doesn't really make sense in this context to talk about
the Fisher matrix of the prior p(θ), which I think is what Mike meant
by the likelihood-prior separation implied by the Fisher matrix.

Matt

> --
>
>

[Michael Betancourt]

Well, there's the decomposition into two terms, but also the addition
of some variables as "data" which one can integrate the likelihood
against to yield the Fisher metric. How would one define the Fisher
metric for any distribution _not_ coming from a posterior?

> --
>
>

[Marcus Brubaker]

The definitional problems with the Fisher information matrix for a
generic distribution is a good point. My point was (though perhaps
poorly made) mostly that the Hessian is not what we want exactly even
if it is tempting.

I actually think something similar to the Fisher information matrix
could be defined for a generic posterior with a little bit of thought.
Plus there's been a fair amount of work in the NA community on
constructing positive definite approximations to local Hessians that
could be really fruitful to explore in this context. In some ways,
this is a big motivation for implementing a user specific metric block
in Stan: to allow people like us to fiddle with possible definitions.

Cheers,
Marcus

On Tue, Sep 18, 2012 at 6:05 PM, Michael Betancourt

> --
>
>

[Matt Hoffman]

> could be really fruitful to explore in this context. In some ways,
> this is a big motivation for implementing a user specific metric block
> in Stan: to allow people like us to fiddle with possible definitions.

That was certainly my main motivation for suggesting it!

Matt