no-U-turn condition on diagonal mass matrix?

[Bob Carpenter]

If I understand Michael B correctly, there's currently
a (serious?) issue in the way the no-U-turn condition
is defined for our default mode of operation (diagonally
adapted NUTS).

Quick question: does this destroy detailed balance or just
negatively impact efficiency?

I don't have an opinion on how this one gets resolved,
but I was hoping that Michael and Matt could put their
heads together and come to some kind of agreement on what
fix is needed. I'm happy to pitch in on the coding and
testing side.

Given this and the fact that smoothing hasn't been
finalized (or has it?) for the full mass matrix adaptation,
I haven't tried to merge Yuanjun's branch into the trunk.
But I'd really really like to get it in sooner rather than
later because of how well it seemed to work for Jiqiang's
multilevel GLM.

Given that we're working on more mass-matrix related approaches,
all else being equal it's simpler conceptually (both for me to
understand and for the doc) if we could use the mass matrix
approach everywhere.

- Bob

[Matt Hoffman]

Detailed balance is not destroyed. NUTS's detailed balance relies on
only two aspects of the leapfrog integrator:
1. reversibility
2. volume-preservation
both of which are preserved by nuts_diag (once it stops adapting, of course).

The criterion for nuts_diag is the classic U-turn criterion, which is
probably suboptimal—it should probably do the transformation implied
by using different step sizes when computing its dot products.

Matt

> --
>
>

[Marcus Brubaker]

I suspect that for nuts_diag (and nuts_nondiag) the correct U-turn criteria is to check the dot product between the gradient and the velocity vector (i.e., inverse mass matrix times momentum) which I think is what Matt is proposing.

For the identity mass matrix case, these two dot products would have the same sign (if not the same magnitude) so it didn't matter then. It should be a relatively easy fix.

--

[Michael Betancourt]

On Jan 18, 2013, at 4:44 PM, Matt Hoffman wrote:

> The criterion for nuts_diag is the classic U-turn criterion, which is
> probably suboptimal—it should probably do the transformation implied
> by using different step sizes when computing its dot products.

Correct, this is about what Yuanjun was concerned. Should be a simply fix for the current code (add the current step-size vector as a parameter to the criterion function and do an entry wise product before projecting onto p).

For dense mass matrix computations the multiple step-size trick is no longer valid and we'll have to appeal to the full P = (L^{T})^{-1} p, Q = L^{T} q transform.

[Bob Carpenter]

Thanks everyone -- that's clear enough that I can patch
the code for the diagonal case. I'll run the changes by
the list again when I'm done.

- Bob

On 1/18/13 5:05 PM, Michael Betancourt wrote:
>
> On Jan 18, 2013, at 4:44 PM, Matt Hoffman wrote:
>
>> The criterion for nuts_diag is the classic U-turn criterion, which is

>> probably suboptimal�it should probably do the transformation implied

[Bob Carpenter]

That took longer to e-mail about than to code
and test :-)

Here's what I did. I haven't committed yet:

inline static bool compute_criterion(std::vector<double>& xplus,
std::vector<double>& xminus,
std::vector<double>& mplus,
std::vector<double>& mminus,
std::vector<double>& step_sizes) {
std::vector<double> total_direction;
stan::math::sub(xplus, xminus, total_direction);
for (size_t i = 0; i < total_direction.size(); ++i)
total_direction[i] /= step_sizes[i];
return stan::math::dot(total_direction, mminus) > 0
&& stan::math::dot(total_direction, mplus) > 0;
}

...

criterion = compute_criterion(xplus, xminus, mplus, mminus,_step_sizes);
...
criterion &= compute_criterion(xplus, xminus, mplus, mminus, _step_sizes);

The alternative to passing step_sizes, which is a member variable
in the nuts_diag class, would be to have compute_criterion be
a member function (method) rather than a static member function and
using step_sizes_[i] in the definition. Any reason to prefer one approach
the other? I know I should favor local variables, but this involves
sending in a reference, whereas the member variable is just a plain
old std::vector<double>.

- Bob

[Matt Hoffman]

Looks right to me. I'd like to have some kind of empirical evidence
that it doesn't fail spectacularly for some unanticipated reason,
though.

I like the functional rather than object-oriented implementation.
Seems more reusable somehow.

Also, why haven't there been more intense holy wars between functional
and object-oriented programming partisans? Shouldn't the ML people be
at the Java people's throats and vice versa? Or is it an "enemy of my
enemy [c++] is my friend" thing?

Matt

>>>> probably suboptimal—it should probably do the transformation implied

>>>> by using different step sizes when computing its dot products.
>>>
>>>
>>> Correct, this is about what Yuanjun was concerned. Should be a simply
>>> fix for the current code (add the current
>>> step-size vector as a parameter to the criterion function and do an entry
>>> wise product before projecting onto p).
>>>
>>> For dense mass matrix computations the multiple step-size trick is no
>>> longer valid and we'll have to appeal to the
>>> full P = (L^{T})^{-1} p, Q = L^{T} q transform.
>>>
>>
>

> --
>
>

[Michael Betancourt]

> I like the functional rather than object-oriented implementation.
> Seems more reusable somehow.
>
> Also, why haven't there been more intense holy wars between functional
> and object-oriented programming partisans? Shouldn't the ML people be
> at the Java people's throats and vice versa? Or is it an "enemy of my
> enemy [c++] is my friend" thing?

I've only ever heard of functional implementations being really slow compared to the object-oriented stuff. Now granted much of that comes from trying to shove functional methods into places they don't belong (of course functional Python code is going to be slow) but that kind of anecdotal evidence seems to have soured people, perhaps prematurely.

[Yuanjun Gao]

Dear Bob,

I think the new code is the right one.

I am not quite familiar with the static member function. I am just a little worried because for nondiag_nuts and maybe Geonuts, we will need to use another version of stopping criterion, which uses a mass matrix rather than a step-size vector(although I didn't fully understand what kind of stopping criterion is Geonuts using...). Will these functions conflict with each other when both are claimed to be static and are using the same name? I am worried because when I try to use the new stopping criterion function which is claimed to be static, the compiler complains.

I do think that the Neal's trick can be generalized to densed mass matrix(the detail of which is attached). We just need substitute the step-sizes to the Cholesky decomposition(and sometimes the transposition of the decomposition) of the mass matrix. Originally I find this trick quite useful because this trick doesn't require any matrix inversion at all, but now it turns out that in order to compute the stopping criterion, a matrix inversion is inevitable, so maybe that trick won't bring efficiency to the algorithm.

Right now I changed my criterion in nondiag_nuts to be

inline bool compute_criterion(std::vector<double>& xplus,

std::vector<double>& xminus,
std::vector<double>& mplus,

std::vector<double>& mminus) {

std::vector<double> total_direction;
stan::math::sub(xplus, xminus, total_direction);

Eigen::Map<Eigen::VectorXd> total_direction_mat(&total_direction[0],total_direction.size());
Eigen::Map<Eigen::VectorXd> mplus_mat(&mplus[0],mplus.size());
Eigen::Map<Eigen::VectorXd> mminus_mat(&mminus[0],mminus.size());
_cov_L.triangularView<Eigen::Lower>().solveInPlace(total_direction_mat);
return total_direction_mat.dot(mplus_mat) > 0
&& total_direction_mat.dot(mminus_mat) > 0;
}

which I am quite sure is the right thing.

[Yuanjun Gao]

As for estimating the full mass matrix, I am really not sure how to "finalize" it because what I can do is just trying to tune the several parameters in many ways and try to fit several models... I begin to feel that maybe the "instability" problem can be partly attributed to the wrong stopping criterion I was using. For a diagonal mass matrix, the error is not that serious because no rotation is involved, but for nondiag_nuts, using the wrong stopping criterion seems to affect a lot.

Hopefully I can try to come up with a better version next week.

[Bob Carpenter]

You can't have two functions with the same fully
qualified name. The fully qualified name includes
the namespace, and if it's a member function (non-static
function defined within a class or struct), the name
of the class, plus the function name + the argument
type sequence. The return type doesn't matter.

So you just need to name your function something else.
Or make it a non-static member function.

You could probably make this a bit more efficient by
passing in a reference to an Eigen::VectorXd to
a slightly differently defined stan::math::sub() that populates
it directly.

void stan::math::sub(..., Eigen::VectorXd& total_direction) {
... write into total_direction ...

Then you can call it presized (which you should also do with total_direction
std::vecotr if you leave it the way you have it):

Eigen::VectorXd total_direction(xplus.size());

stan::math::sub(xplus, xminus, total_direction);

I'd also suggest renaming variables that you call _mat
to _vec, because although it's technically correct that
everything's a matrix in Eigen, I found it confusing reading
the code.

In fact, compute_criterion could itself be changed to
use Eigen vectors rather than std::vector. It would make
much of the new code cleaner (wouldn't require the
map wrappers and all that reference to first array business,
which is the only way to do things now).

We could actually rewrite everything to use Eigen vectors,
which would probably be more efficient. I wrote it the way
it is now thinking we could write almost everything without
a dependency on Eigen (for instance, if you look at the model
specifications, they only depend on std::vector, not on anything
in Eigen). But since then, for both efficiency
and for expressiveness, we've used more and more of Eigen
throughout the mcmc package. So I don't think it's worth
worrying about keeping all these things separate. The same
thing happened to intertwine the probability classes and
the auto-dif classes.

- Bob

[Bob Carpenter]

To lay some ugly rumors to rest,

* functional languages are not slow compared to
object-oriented languages,

* each language has its strengths and weaknesses ranging from
speed to library availability to portability to
size of developer community to steepness of learning curve,

* as of the later 1.6 JVMs, Java is pretty close to C++
in basic arithmetic, array operations, function dispatching,
etc.,

* Fortran and some functional languages (like ML) and some styles
of C++ programming with templates can be very fast because
the compiler can do so much static resolution,

* folded code can be surprisingly fast in Java due to the
just-in-time compiler's ability to dynamically unfold loops
based on run-time analysis, and

* interpreted languages (Python, R, MATLAB) are very slow to loop
and make function calls, but if one function call does a lot
of work (typically calling vectorized operations in Fortran or C(++)),
then the added efficiency of a compiled language may not even
be noticeable.

What is slow is

* unpredictable branch points,

* including virtual function calls (including most Java non-static
method calls, C++ virtual function calls, C function
pointer calls), and

* accessing memory out of order.

- Bob

P.S. I can almost always be baited into a programming language
discussion. So please, help me stop.

[Bob Carpenter]

On 1/19/13 11:02 PM, Matt Hoffman wrote:
> Looks right to me. I'd like to have some kind of empirical evidence
> that it doesn't fail spectacularly for some unanticipated reason,
> though.

Absolutely. We'll run through all the model tests for
the BUGS models.

What I'm hoping is that by adjusting for step sizes, we'll get
more balanced effective sample sizes across parameters. The danger
I've always worried about is that the larger scale parameters would
dominate the U-turn condition. I brought it up in the last meeting
and Michael explained why it wouldn't be an issue in RM-HMC.

> I like the functional rather than object-oriented implementation.
> Seems more reusable somehow.

Non-class functions seems to be the standard idiom in C++
for helper functions. The reasoning I've seen is
that it makes compilation more modular. I still have habits
hanging over from Java where there aren't any non-class
functions.

Anyway, I coded it that way. When in Rome and all.

- Bob

[Matt Hoffman]

I'm torn on this last question. I think there are definitely lots of
places where we use std::vector<double> where we should be using
Eigen's vectors*. On the other hand, I do like isolating dependencies
as much as possible. But maybe the benefits outweigh the aesthetic
costs, at least within the MCMC package. It's kind of silly that we
have our own addition, subtraction, dot product, etc. operations in
the leapfrog integrator, NUTS criterion, etc.

Of course, it'd be a decent amount of refactoring work.

Matt

* Normally I wouldn't possessivise like that, but it eliminates the
homophonic ambiguity of "Eigen vector".

> --
>
>

[Marcus Brubaker]

I would strongly support using Eigen::Matrix in place of std::vector<double> everywhere.  Eigen is a light, header only dependence so I'm hard pressed to see a reason to keep it separated.

Cheers,

Marcus

--

[Bob Carpenter]

On 1/22/13 2:07 PM, Marcus Brubaker wrote:
> I would strongly support using Eigen::Matrix in place of std::vector<double> everywhere. Eigen is a light, header only
> dependence so I'm hard pressed to see a reason to keep it separated.

I agree. My original motivation was that someone
might want to use different matrix libs with our auto-dif
or just not have to learn yet another set of interfaces
in addition to the STL.

But now that we're doing actual matrix operations with
the vectors and we've tied Eigen vectors into our multivariate
distributions, I support moving to Eigen data structures.

- Bob

[Ben Goodrich]

On Tuesday, January 22, 2013 2:07:00 PM UTC-5, Marcus Brubaker wrote:

I would strongly support using Eigen::Matrix in place of std::vector<double> everywhere.  Eigen is a light, header only dependence so I'm hard pressed to see a reason to keep it separated.

I agree in spirit. For a while toward the beginning of Stan, I was trying to maximize the use of Eigen, while Bob was trying to minimize it. The only thing I would add is that I think an Eigen::Matrix<> or Eigen::Vector<> should be used instead of std::vector<> iff linear algebra is involved. If it is just holding values or being utilized elementwise, then an Eigen::Array<> should be used.

Ben

[Bob Carpenter]

On 1/22/13 7:51 PM, Ben Goodrich wrote:
> On Tuesday, January 22, 2013 2:07:00 PM UTC-5, Marcus Brubaker wrote:
>
> I would strongly support using Eigen::Matrix in place of std::vector<double> everywhere. Eigen is a light, header
> only dependence so I'm hard pressed to see a reason to keep it separated.
>
>
> I agree in spirit. For a while toward the beginning of Stan, I was trying to maximize the use of Eigen, while Bob was
> trying to minimize it. The only thing I would add is that I think an Eigen::Matrix<> or Eigen::Vector<> should be used
> instead of std::vector<> iff linear algebra is involved.

There's always been a dot product involved. It'll let us
remove some functions if we use Eigen. And the Eigen ones
may be a bit more efficient because they're pretty careful
on the implementation details.

> If it is just holding values or being utilized elementwise,
> then an Eigen::Array<> should be used.

What's the advantage of Array? I thought Array and
Matrix were lightweight wrapper classes holding a pointer
to the data and dimensions. I know that Array supports
more elementwise ops and fewer matrix ops, but I'd think
they'd be of the same efficiency.

I really like how Eigen lets you map its data, so it's
pretty cheap to convert back and forth (just a constructor
call).

- Bob