parameter correlations in posterior?

[Bob Carpenter]

My understanding is that HMC with unit mass matrix
is invariant to rotations of the posterior, but not scaling.
Gibbs is the opposite (invariant to scale, but not rotation).

But when we do adaptation of a diagonal mass matrix, we run
into the problem that we break rotational symmetry in the algorithm.
Is that why correlated parameters can be problematic for HMC?
If we adapt a dense mass matrix, do we regain rotation invariance
(up to the formidable issue of estimating that mass matrix).

I'm not worried (yet) in this discussion about varying curvature.

- Bob

[Andrew Gelman]

I know from bitter experience that linear transformations matter for HMC. If you run a simple linear regression, y = a + bx + error, and x has a mean that is far from 0, it can run really really really slow. But if you center x, it runs fast.
A

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

There is a lot of confusion around these points, so let me
start from the beginning.

1) In theory Hamiltonian Monte Carlo is fully geometric,
which means it is independent to whatever choice of
parameters we use. Gibbs sampling does depend on
the parametrization, as the exact parametrization defines
the conditionals.

Operationally, if we change the parameters (potential
energy) then the momenta (kinetic energy) have to change
in the opposite way so that the Hamiltonian, and hence
the trajectories, remain invariant.

Provided that we transform the parameters and momenta
together then all versions of Hamiltonian Monte Carlo are
invariant to any well-behaved transformation, including
rotations and scalings.

Interestingly enough, Random walk Metropolis behaves
the same way! When we transform the parameters we
also have to transform the proposal covariance or _we’re
specifying a different transition_!

2) This is confusing because in practice people often
transform the parameters but leave the momenta (kinetic
energy) the same. Because we’re not transforming them
together, this implies that we’re specifying a _different_
kinetic energy.

For example, a “unit mass matrix” is defined only relative
to a given parameterization. If we transform the system
then both the parameters and momenta change, including
the mass matrix to something non-unit. Forcing the
mass matrix to stay unit while changing the parameters
is implicitly defining a different kinetic energy.

The same thing happens with Random walk Metropolis.
If we keep the covariance constant but change the
parameters then we’re effectively defining a completely
different transition and can expect different performance.

3) All of this is complicated by adaptation. If we define
the mass matrix to be the global covariance of the target
distribution then we have no freedom to transform the
parameters and momenta separately. Change the
parameters and the form of the global covariance
changes, hence the mass matrix changes. But this
leaves the total system invariant as in (1).

In summary, all of this is confusing because people often
talk about the invariance of the posterior _relative to a fixed
kinetic energy and mass matrix_. But that’s an odd construction
because the parameters and momenta/kinetic energy/mass matrix
are supposed to transform together, so fixing the kinetic energy
and changing the parameterization is just a sneaky way to specify
an mass matrix.

The performance of HMC depends on how the target distribution
interacts with the kinetic energy. If the target distribution is
correlated in a parameterization but the kinetic energy is defined
with a unit mass matrix then the trajectories will be highly
oscillatory and expensive to simulate. If the parameterization
is fixed then rotating and scaling the mass matrix can define
a more optimal kinetic energy that reduces the oscillatory nature
of these trajectories.

Ultimately one can show that the optimal Euclidean mass matrix
is given by the global covariance of the target distribution, which
is invariant to any particular parameterization as discussed above.

Okay, that probably didn’t make any sense so we’ll have to iterate.

[Daniel Lee]

Whoa.

That clarified a lot for me. Thanks!


Daniel

[Bob Carpenter]

> On Apr 5, 2016, at 4:36 PM, Michael Betancourt <betan...@gmail.com> wrote:
>

> ...

> Ultimately one can show that the optimal Euclidean mass matrix
> is given by the global covariance of the target distribution, which
> is invariant to any particular parameterization as discussed above.

Cool. That's exactly what I suspected, but hadn't put all the
coupled invariances together into the Hamiltonian.

> Okay, that probably didn’t make any sense so we’ll have to iterate.

Thanks!

I actually think I understood all that and it indirectly
answered my question in the affirmative. I was specifically
thinking about unit mass matrices in situations where posterior covariance
is not position dependent, and changing the parameterization
only by rotation. (I think that's kosher, but I wouldn't bet on
my knowledge of Hamiltonians, exponential families, or linear algebra).

- Bob

[Michael Betancourt]

Given a parameterization, the only unconstrained distribution

specified by a global covariance is the Gaussian,

so we can reduce this discussion to Gaussian targets.

The Hamiltonian in this case is given by

H = 0.5 * p^{T} M^{-1} p + 0.5 * (q - mu)^{T} Sigma^{-1} (q - mu)

Now we can apply an orthogonal transformation to rotate the q/mu,

but to maintain the same system we have to do the opposite

to the momenta,

q -> X q

p -> X^{T} p

which gives

H = 0.5 * p^{T} X M^{-1} X^{T} p + 0.5 * (q - mu)^{T} X^{T} Sigma^{-1} X (q - mu)

The optimal transform completely decouples this into N one-dimensional

(1 q + 1 p) systems.  We start by decoupling the potential energy — if 

Sigma = L L^{T} then we want X = L so that 

H = 0.5 * p^{T} L M^{-1} L^{T} p + 0.5 * (q - mu)^{T} (q - mu)

Then the kinetic energy decouples if M = Sigma.  In other words,

the only way to decouple a correlated Gaussian is to set the

mass matrix to the global covariance of the target distribution.

[Bob Carpenter]

Cool. Thanks. I had no idea that the Gaussian was
the only distribution with global covariance.

I was thinking unit Gaussian (hence the unit mass matrix).
Clearly that is invariant to rotation because the density
doesn't change in Euclidean coordinates under rotations.
I hadn't quite put all the pieces together to realize how
limited the case was that I was considering.

- Bob

[Michael Betancourt]

Well, no, it does change — the covariance rotates from being
diagonal to being dense!

The Gaussian is the only unconstrained distribution is that
_completely_ specified by the first and second moments,
which means that it’s the only distribution where the global
covariance will be equal to one of the natural parameters
in the density (all of this conditioned on a choice of parameters
or dominating measure or specific R^{n}, or specific Lebesgue
measure, etc, etc).

[Bob Carpenter]

> On Apr 6, 2016, at 4:33 PM, Michael Betancourt <betan...@gmail.com> wrote:
>
> Well, no, it does change — the covariance rotates from being
> diagonal to being dense!

I thought multi_normal(0, I) was rotationally symmetric.
Meaning the density doesn't change under rotation.
And isn't the set of rotations of the identity matrix just the
set of correlation matrices?

So if the density doesn't change, I'd think sampling wouldn't change.

But then when I think about it, if I have data y, then it
would seem that the correlation matrix should impact the
density. And clearly this Stan program doesn't always return
the same value for lpdf:

data {
vector[4] y[10];
}
parameters {
cholesky_factor_corr[4] Omega;
}
model {
}
generated quantities {
real lpdf;
lpdf <- multi_normal_log(y, rep_vector(0, 4), Omega * Omega');
}

So I clearly have the wrong end of some stick.

- Bob

[Michael Betancourt]

Covariance is the variance between _parameters_, which
depends on your chosen parameters. A rotation is a change
of parameters which changes the covariance.

What’s staying the same in all of this is the latent probability
distribution — we’re just changing the sample space that we
use to represent that distribution. This is why I had to go so
high-level in the intro draft.