differences in dual-averaging

[Daniel Lee]

Michael,

I'm seeing some massive differences in the implementation of the dual-averaging. This could be the cause of some of our problems. Let me know when you want to talk through it.

Daniel

[Michael Betancourt]

In what sense? The additional step size adaptation window should significantly improve performance - is it massively better?

My Internet situation is sketchy right now but I can correspond via email.  

--
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.
For more options, visit https://groups.google.com/groups/opt_out.
 
 

[Daniel Lee]

Before we get to dual-averaging, I'm seeing differences in:

v1.3.0: find_reasonable_parameters()

master: init_stepsize()

For the simple normal(0,1) model (src/models/command.stan), I'm seeing an epsilon for v1.3.0 at 2 and a nom_epsilon for master at 3.2.

Are these supposed to behave identically?

[Daniel Lee]

So, before we get into improving performance, I'm more concerned that we can't replicate things between mac/windows AND current behavior is different from the past.

[Michael Betancourt]

Yes, I added stepwise jitter in between.  So where before we just had stepsize now we have a nominal stepsize and the sampled stepsize.

[Daniel Lee]

In init_stepsize() in base_hmc, should:

        double delta_H = H0 - this->_hamiltonian.H(this->_z);
        int direction = delta_H > log(0.5) ? 1 : -1;

match:

        double H = logp - lastlogp;
        int direction = H > log(0.5) ? 1 : -1;

from find_reasonable_parameters() in hmc_base.hpp?

Right now, delta_H in the new code is:
H0: 1.257892961
H:  1.257924408
delta_H: -3.144732402e-05

where as in the old code:
lastlogp: -1.99975
logp: -0.171807
H: 1.82794

These are with the same seeds, same init values, etc.

[Michael Betancourt]

The last commit was just to fix a poor implementation in the refactored adaptation.  The performance should be comparable to before the refactor.  

The difference between Windows and OS X still baffles me.  It almost looks like a small floating point operator difference or rounding difference.

[Michael Betancourt]

No, they shouldn't match.  For some reason Matt used the log prob when determining the HMC acceptance rate here, which is not the best way to determine the optimal HMC stepsize.  The new implementation uses the full Hamiltonian, which is correct.

[Daniel Lee]

Ok, so the epsilons shouldn't match then. (Let me know if that's correct.)

If we're past this step and the epsilons do match, should we be expecting identical results?

[Michael Betancourt]

Right, in general they will be different.  If they are the same and the metrics are the same (unlikely since they are factored differently) then the results will still be slightly different due to the jitter, but you can suppress that with --epsilon_pm=0 on both versions.  If you want to focus on the stepsize adaptation then just use --equal_step_sizes

[Daniel Lee]

I found an issue:

the initial value isn't passed to sampler.init_stepsize() (base_hmc). So, the q value of _z is always 0 to start. The p is sampled properly. Michael, what's the best way of passing that into the sampler? seed()?

[Michael Betancourt]

init_stepsize is fine - the sampler itself just isn't initialized to the right position. A call to seed before init_stepsize in gm/command.hpp will fix that.  

[Daniel Lee]

do you want to add it in? I'm not sure I'm in tune with all the places it should / shouldn't be added.

[Bob Carpenter]

This is why it's better to construct objects in a
sound state and leave them that way! Otherwise, their
usage contracts become more complex (construct, then call
function f() before ever calling function g() or h(), etc.)

- Bob

On 5/6/13 3:13 PM, Michael Betancourt wrote:
> init_stepsize is fine - the sampler itself just isn't initialized to the right position. A call to seed before
> init_stepsize in gm/command.hpp will fix that.
>

> On May 6, 2013, at 8:02 PM, Daniel Lee <bea...@alum.mit.edu <mailto:bea...@alum.mit.edu>> wrote:
>
>> I found an issue:
>>
>> the initial value isn't passed to sampler.init_stepsize() (base_hmc). So, the q value of _z is always 0 to start. The
>> p is sampled properly. Michael, what's the best way of passing that into the sampler? seed()?
>>
>>
>>
>> On Mon, May 6, 2013 at 1:49 PM, Michael Betancourt <betan...@gmail.com <mailto:betan...@gmail.com>> wrote:
>>
>> Right, in general they will be different. If they are the same and the metrics are the same (unlikely since they
>> are factored differently) then the results will still be slightly different due to the jitter, but you can
>> suppress that with --epsilon_pm=0 on both versions. If you want to focus on the stepsize adaptation then just use
>> --equal_step_sizes
>>

>>>>> stan-dev+u...@googlegroups.com <mailto:stan-dev+u...@googlegroups.com>.

>>>>> For more options, visit https://groups.google.com/groups/opt_out.
>>>>>
>>>>>
>>>> --
>>>> 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 <mailto:stan-dev%2Bunsu...@googlegroups.com>.

>>>> For more options, visit https://groups.google.com/groups/opt_out.
>>>>
>>>>
>>>>
>>>>
>>>>
>>>> --
>>>> 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 <mailto:stan-dev+u...@googlegroups.com>.

>>>> For more options, visit https://groups.google.com/groups/opt_out.
>>>>
>>>>
>>> --
>>> 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 <mailto:stan-dev%2Bunsu...@googlegroups.com>.

>>> For more options, visit https://groups.google.com/groups/opt_out.
>>>
>>>
>>>
>>> --
>>> 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 <mailto:stan-dev+u...@googlegroups.com>.

>>> For more options, visit https://groups.google.com/groups/opt_out.
>>>
>>>
>> --
>> 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 <mailto:stan-dev%2Bunsu...@googlegroups.com>.

>> For more options, visit https://groups.google.com/groups/opt_out.
>>
>>
>>
>> --
>> 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 <mailto:stan-dev+u...@googlegroups.com>.

[Matt Hoffman]

Some questions:

Is jittering turned on by default? If so, why?

What exactly are you passing to the step size adaptation to tell it
whether the step size is too small or too large?

How are you deciding whether to use a sample from the new half-tree or
the old half-tree?

I'm concerned that the "refactor" has made some significant changes to
the actual algorithm, and that these changes haven't been adequately
tested.

Matt

[Michael Betancourt]

> Is jittering turned on by default? If so, why?

Defaults to zero jitter, just as before.

> What exactly are you passing to the step size adaptation to tell it
> whether the step size is too small or too large?

The average acceptance probability of the tree, just as before.

> How are you deciding whether to use a sample from the new half-tree or
> the old half-tree?

Uniformly based on total valid points in the subtree, just as before.

> I'm concerned that the "refactor" has made some significant changes to
> the actual algorithm, and that these changes haven't been adequately
> tested.

I may very well have made an error somewhere, but every attempt was made to reproduce the old algorithm exactly unless I found an explicit error (as in the stepsize initialization - I can't recall another instance). In particular I did a bunch of tests on the interaction between stepsize and tree structure and everything behaved as expected.

[Michael Betancourt]

Sure, but we don't normally use that state outside of actual sampling because it was decided that the sampler would only transform an input sample and not carry the state of the Markov chain internally (technically the state is stored, but it's always overwritten). Having adaptation really motivates the other approach, where the constructor you note would be much more natural.

[Michael Betancourt]

Let me see if I can find a strong enough wireless signal to get my laptop up (only chance is the bathroom...)

[Matt Hoffman]

Ok, I'll take a closer look at the code and double-check. I got a
little confused by the util.sum_prob, but now I see that it's not
being used to decide which branch to sample from but only as a measure
of integrator accuracy. One minor difference, though, is that in the
old version we only computed the average probability for the most
recent half of the tree, excluding the states in the half-tree nearest
the starting point. This is because the first state is always
acceptable, so including states near the start state biases our
estimate of the average acceptance probability upwards (though
probably only slightly). I think the original behavior could be
brought back by moving

util.n_tree = 0;
util.sum_prob = 0;

into the while loop.

BTW, can someone remind me how to set off code so that it looks nice
in these emails? Ben's code snippets always look so much nicer than
mine...

Matt

[Ben Goodrich]

On Monday, May 6, 2013 3:34:25 PM UTC-4, Matt Hoffman wrote:

BTW, can someone remind me how to set off code so that it looks nice
in these emails? Ben's code snippets always look so much nicer than
mine...

I just use the Google Groups interface, which has an icon that looks like {}. There used to be the same thing in Gmail itself, but it disappeared. I think the raw email may use the wiki syntax of triple braces. So, this

e = mc^2

should definitely look formatted to you, and this

{{{
e = mc^2
}}}

might be.

Ben

[Ben Goodrich]

On Monday, May 6, 2013 1:14:38 PM UTC-4, Daniel Lee wrote:

So, before we get into improving performance, I'm more concerned that we can't replicate things between mac/windows AND current behavior is different from the past.

What is a short, simple thing I can do to look for differences (I see the branch)? Is the Linux behavior known to be the same as Mac? Is in known that the differences stem from Windows the operating system vs. differences in how mingw64 compiles it?

Ben

[Daniel Lee]

Here is what I'm doing to generate the files:

make models/command
./models/command --data=src/test/models/command1.data.R --seed=1000 --iter=1000 > out_mac

If you want to do that and pipe to out_linux, we can see if there are differences there too.

[Ben Goodrich]

On Monday, May 6, 2013 5:00:03 PM UTC-4, Daniel Lee wrote:

Here is what I'm doing to generate the files:

make models/command
./models/command --data=src/test/models/command1.data.R --seed=1000 --iter=1000 > out_mac

If you want to do that and pipe to out_linux, we can see if there are differences there too.

Yes, there is a compiler difference on Linux. In pure g++4.6 with 64bit code, I get

Inference for Stan model:
1 chains: each with iter=(500); warmup=(0); thin=(1); 500 iterations saved.

            mean se_mean   sd 2.5%  25%  50%  75% 97.5% n_eff Rhat
y           -0.1     0.1  0.9 -2.1 -0.7 -0.0  0.5   1.6   260  1.0
lp__        -0.4     0.1  0.7 -2.4 -0.5 -0.2 -0.0  -0.0   167  1.0
stepsize__   2.0     0.0  0.0  2.0  2.0  2.0  2.0   2.0   500  1.0

after cross-compiling with 32bit i686-w64-mingw32-g++, I get

Inference for Stan model:
1 chains: each with iter=(500); warmup=(0); thin=(1); 500 iterations saved.

            mean se_mean   sd 2.5%  25%  50%  75% 97.5% n_eff Rhat
y            0.0     0.1  1.0 -1.9 -0.6  0.0  0.5   2.0   357  1.0
lp__        -0.5     0.0  0.8 -2.3 -0.6 -0.1 -0.0  -0.0   314  1.0
stepsize__   1.5     0.0  0.0  1.5  1.5  1.5  1.5   1.5   500  1.0

Clang++-3.3 agrees with pure g++-4.6, except that if I turn on -fsanitize=undefined I do get two (non-fatal) runtime errors:

src/stan/mcmc/hmc/nuts/base_nuts.hpp:119:70: runtime error: division by zero
src/stan/mcmc/hmc/nuts/base_nuts.hpp:212:56: runtime error: division by zero

Ben

[Daniel Lee]

Does the 64 bit version match out_mac and the 32 bit version match out_win exactly?

[Daniel Lee]

These are two files generated using the same seeds in mac and windows using v1.3.0.

If you take a look at the files, they don't differ until iteration 241. I was mistaken thinking they stayed in lock-step through the whole simulation, but they do stay in synch for much longer.

With the current master, we can't get the first samples to match. That's what I'm trying to track down.

[Ben Goodrich]

On Monday, May 6, 2013 5:46:01 PM UTC-4, Daniel Lee wrote:

Does the 64 bit version match out_mac and the 32 bit version match out_win exactly?

Aside from the elapsed times, my two files match each other, although I don't have actual results on Mac or Windows to compare them to.

Ben
 

[Daniel Lee]

In the diff_win_mac branch, there are two files: out_win and out_mac

[Ben Goodrich]

On Monday, May 6, 2013 5:28:31 PM UTC-4, Ben Goodrich wrote:

Clang++-3.3 agrees with pure g++-4.6, except that if I turn on -fsanitize=undefined I do get two (non-fatal) runtime errors:

src/stan/mcmc/hmc/nuts/base_nuts.hpp:119:70: runtime error: division by zero
src/stan/mcmc/hmc/nuts/base_nuts.hpp:212:56: runtime error: division by zero

 
I can be more specific about these if I don't redirect stdout:

Iteration:   1 / 1000 [  0%]  (Warmup)

src/stan/mcmc/hmc/nuts/base_nuts.hpp:119:70: runtime error: division by zero

Iteration:   2 / 1000 [  0%]  (Warmup)
Iteration:   3 / 1000 [  0%]  (Warmup)

src/stan/mcmc/hmc/nuts/base_nuts.hpp:212:56: runtime error: division by zero

Iteration:   4 / 1000 [  0%]  (Warmup)
...

Ben

[Daniel Lee]

What seed?

[Daniel Lee]

Nevermind. You're absolutely right,

n_valid_subtree + n_valid == 0

I get that with seed=1000

[Ben Goodrich]

On Monday, May 6, 2013 6:36:09 PM UTC-4, Daniel Lee wrote:

Nevermind. You're absolutely right,

n_valid_subtree + n_valid == 0

I get that with seed=1000

Yes, I was always using the seed you gave me.

Ben

[Ben Goodrich]

On Monday, May 6, 2013 6:27:07 PM UTC-4, Daniel Lee wrote:

In the diff_win_mac branch, there are two files: out_win and out_mac

Okay, I was just doing it on the master branch before. There are differences between out_mac and out_linux_g++.txt (attached), but at first only in the last (few) decimal places. They gradually accumulate up to more relevant decimal places and on iteration 165, you get a much different "transition's init_sample" and everything is different after that.

Ben

[Daniel Lee]

This line in base_nuts.hpp looks suspect:

 return (util.log_u + (h - util.H0) < 0);

That's the return for the base case in build_tree(...).

I think something like util.n_tree should be returned, but maybe not?

[Bob Carpenter]

On 5/6/13 10:54 PM, Daniel Lee wrote:
> This line in base_nuts.hpp looks suspect:
>
> return (util.log_u + (h - util.H0) < 0);
>
> That's the return for the base case in build_tree(...).
>
> I think something like util.n_tree should be returned, but maybe not?

I don't think so. This is the base case, so the
answer should be 0 or 1. The boolean returned by
operator<() will be cast to an int with value 0 or 1.
So the above is equivalent to:

return (util.log_u < util.H0 - h) ? 1 : 0;

- Bob

>
>
>
> On Mon, May 6, 2013 at 6:42 PM, Ben Goodrich <goodri...@gmail.com <mailto:goodri...@gmail.com>> wrote:
>
> On Monday, May 6, 2013 6:36:09 PM UTC-4, Daniel Lee wrote:
>
> Nevermind. You're absolutely right,
>
> n_valid_subtree + n_valid == 0
>
> I get that with seed=1000
>
>
> Yes, I was always using the seed you gave me.
>
> Ben
>
>
> --
> 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 <mailto:stan-dev%2Bunsu...@googlegroups.com>.

[Ben Goodrich]

Since I get different behavior on the same machine with different compilers, I should be able to bisect it. However, there are so many commits where the attached test is inconclusive due to the model failing to build or run, so there is nothing to compare. Thus, we can only narrow it down to a range of commits, that seem to all pertain to the refactoring:

There are only 'skip'ped commits left to test.
The first bad commit could be any of:
f7b426b10b0a13d76675057ec3c9aacea7d9e0c7
1dc4aed16e63c44562ae29dc6f3bc517c324794e
aaf821ba2387a629e911a6a0406185b70af1a27c
3b6353b932af74489570bf7a9886429159892f61
d0e8d45df36416290e167b81b92cd525505ba160
04b8f207985873eb8f2c9e0cf5fd9e63fa65feb7
c946070d58fd3646574c6a312d80f2322d97aff3
50f7b32642ebaa3e498360f932ad759c10de06f3
f4f0713cfc6c98751c902d1859989fef184df808
6daaea6b3bfbf566539aa1bac5ecb36136cc3437
fab3c0841d0755aaf81a48a2c729b0b3b2b959bf
86bd0e446cd34d258384201218c2f2ea3812e423
ca85e0e6a8bcffd90ca3047a0e540c674fc8e09d
99429292d821d14b1a1c83ef4885701e3d409c7f
e11396f1703a402ea248653ec8e1157d197166af
4385f7ff2af3a8651b3ed3152e8552ab1e27dc0e
c61335e3a8dda8a777c14e326ea86eda1916568e
0fe3f288f05d228a98c68c1cb60ce2805a73ec69
6bea89ee6ac669d3f7d9d85eaa0fa86411e8c601
0df04955917fe8ea69a2d4d80d1918c24be432b7
434b44743e1914093a682614ece759f35d8cda84
d22105bc4eb35d124de3a13a3b08e9775db6913e
We cannot bisect more!
bisect run cannot continue any more

goodrich@CYBERPOWERPC:/opt/stan$ git bisect log
git bisect start
# bad: [951b28226bc16a0c30f45355b1fb9960c958f4df] calculate average n_eff in rstan/example/demo.R
git bisect bad 951b28226bc16a0c30f45355b1fb9960c958f4df
# good: [1cde932d15a348df17e9dafce17a79299b9d8a2d] upgrade rstan verion to 1.2.0
git bisect good 1cde932d15a348df17e9dafce17a79299b9d8a2d
# good: [bac02df9a2d5d6f6be79b4a69cf5d60768961df4] Merge branch 'master' of https://github.com/stan-dev/stan
git bisect good bac02df9a2d5d6f6be79b4a69cf5d60768961df4
# good: [c4a2253bb7159992e5de557343fabf4f26f2d236] merged conflict in cloglog
git bisect good c4a2253bb7159992e5de557343fabf4f26f2d236
# bad: [3a78cab458537ab1bc611fd03e10707828b4e63b] added multiply_lower_tri_self_transpose to fvar/matrix
git bisect bad 3a78cab458537ab1bc611fd03e10707828b4e63b
# bad: [f9ea47fca6731b096bcf9f128869211eca598c41] Switched to Stan's Eigen header
git bisect bad f9ea47fca6731b096bcf9f128869211eca598c41
# skip: [6ed834b3b0a75db62ebb7cc880fc249e72f855ce] Update var/covar to use Welford algorithm, transform to scale-free output
git bisect skip 6ed834b3b0a75db62ebb7cc880fc249e72f855ce
# bad: [d3fddd6b0bab425508ce2f008689a4623bd8c772] Refactored default update into base_hamiltonian, catch exceptions thrown by grad_log_prob
git bisect bad d3fddd6b0bab425508ce2f008689a4623bd8c772
# skip: [6a87a6033a7fe56562f269f7ee0e55abefd3079f] Added diagonal and dense samplers, refactored static HMC algorithm into base class
git bisect skip 6a87a6033a7fe56562f269f7ee0e55abefd3079f
# bad: [e930ba046d55859bac4c8dccb85450287ebe234f] Hopefully NUTS is now working
git bisect bad e930ba046d55859bac4c8dccb85450287ebe234f
# skip: [e11396f1703a402ea248653ec8e1157d197166af] Interface/output tweaks
git bisect skip e11396f1703a402ea248653ec8e1157d197166af
# skip: [1dc4aed16e63c44562ae29dc6f3bc517c324794e] Removed deprecated code
git bisect skip 1dc4aed16e63c44562ae29dc6f3bc517c324794e
# skip: [434b44743e1914093a682614ece759f35d8cda84] Refactoring...
git bisect skip 434b44743e1914093a682614ece759f35d8cda84
# skip: [3b6353b932af74489570bf7a9886429159892f61] Added base MCMC sampler class, sample object for i/o
git bisect skip 3b6353b932af74489570bf7a9886429159892f61
# skip: [4385f7ff2af3a8651b3ed3152e8552ab1e27dc0e] Removed deprecated includes
git bisect skip 4385f7ff2af3a8651b3ed3152e8552ab1e27dc0e
# bad: [d22105bc4eb35d124de3a13a3b08e9775db6913e] Updated for mcmc refactor
git bisect bad d22105bc4eb35d124de3a13a3b08e9775db6913e
# skip: [f7b426b10b0a13d76675057ec3c9aacea7d9e0c7] Start testing the new MCMC code
git bisect skip f7b426b10b0a13d76675057ec3c9aacea7d9e0c7
# skip: [d0e8d45df36416290e167b81b92cd525505ba160] Added apapter class
git bisect skip d0e8d45df36416290e167b81b92cd525505ba160
# skip: [ca85e0e6a8bcffd90ca3047a0e540c674fc8e09d] Fixed bug in random number generator instatiation, moved to holding only a reference to the external RNG, various const and pass-by-reference fixes
git bisect skip ca85e0e6a8bcffd90ca3047a0e540c674fc8e09d
# skip: [aaf821ba2387a629e911a6a0406185b70af1a27c] Separated adapter into hiearchical classes
git bisect skip aaf821ba2387a629e911a6a0406185b70af1a27c
# skip: [86bd0e446cd34d258384201218c2f2ea3812e423] Merged master, preserving development code in mcmc and gm
git bisect skip 86bd0e446cd34d258384201218c2f2ea3812e423
# skip: [04b8f207985873eb8f2c9e0cf5fd9e63fa65feb7] sample() returns acceptance statistic (exact statistic will depend on algorithm).
git bisect skip 04b8f207985873eb8f2c9e0cf5fd9e63fa65feb7
# skip: [99429292d821d14b1a1c83ef4885701e3d409c7f] Development command file
git bisect skip 99429292d821d14b1a1c83ef4885701e3d409c7f
# skip: [c946070d58fd3646574c6a312d80f2322d97aff3] Debugging
git bisect skip c946070d58fd3646574c6a312d80f2322d97aff3
# skip: [fab3c0841d0755aaf81a48a2c729b0b3b2b959bf] Added debugging model
git bisect skip fab3c0841d0755aaf81a48a2c729b0b3b2b959bf
# skip: [50f7b32642ebaa3e498360f932ad759c10de06f3] Debugging
git bisect skip 50f7b32642ebaa3e498360f932ad759c10de06f3
# skip: [0fe3f288f05d228a98c68c1cb60ce2805a73ec69] Templated out phase space point so that it can also hold temporary variables for more sophisticated Hamiltonians
git bisect skip 0fe3f288f05d228a98c68c1cb60ce2805a73ec69
# skip: [f4f0713cfc6c98751c902d1859989fef184df808] Compute gradient only after position has been updated, unify variable naming to camel case convention.
git bisect skip f4f0713cfc6c98751c902d1859989fef184df808
# skip: [c61335e3a8dda8a777c14e326ea86eda1916568e] Restored bulk of previous content, ready for other samplers
git bisect skip c61335e3a8dda8a777c14e326ea86eda1916568e
# skip: [6daaea6b3bfbf566539aa1bac5ecb36136cc3437] Store a reference to the external model, not an explicity copy
git bisect skip 6daaea6b3bfbf566539aa1bac5ecb36136cc3437
# skip: [6bea89ee6ac669d3f7d9d85eaa0fa86411e8c601] Template fixes
git bisect skip 6bea89ee6ac669d3f7d9d85eaa0fa86411e8c601
# skip: [0df04955917fe8ea69a2d4d80d1918c24be432b7] Refactored, separating each class into its own file
git bisect skip 0df04955917fe8ea69a2d4d80d1918c24be432b7
goodrich@CYBERPOWERPC:/opt/stan$

Ben

[Michael Betancourt]

After staring at 15th decimal places for the last 3 hours I found the problem -- I'm going to go grab lunch before fixing it, so if anyone is up early don't bother looking into anything.

Ben

--
You received this message because you are subscribed to the Google Groups "stan development mailing list" group.

[Bob Carpenter]

Awesome!

- Bob

> stan-dev+u...@googlegroups.com <mailto:stan-dev%2Bunsu...@googlegroups.com>.

[Michael Betancourt]

So....I may have been premature.  I found two somewhat serious bugs in the NUTS code -- it's amazing they didn't manifest as more serious sampling problems (and frustrating, because I had extensively checked for those kinds of problems when first writing up the new code -- I have no idea what happened in between).  It seemed like the errors were introducing larger gradients which allowed very small differences in floating point calculations to leak into the sampler, accumulate, and eventually become nontrivial.  But even after fixing those problems and verifying that NUTS behaves as expected, Mac/Windows are different.  I'm still neck deep in code, hopefully something will turn up.

Incidentally, I can at least explain why the Mac/Windows code seems to "diverge".  The numerical differences eventually affect the NUTS criterion, and at a given iteration the tree on one machine reaches a deeper depth than on the other.  This misaligns the random number generators and the resulting sequence of samples.

    stan-dev+unsubscribe@googlegroups.com <mailto:stan-dev%2Bunsu...@googlegroups.com>.

    For more options, visit https://groups.google.com/groups/opt_out.



--
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+unsubscribe@googlegroups.com.

For more options, visit https://groups.google.com/groups/opt_out.

--
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+unsubscribe@googlegroups.com.

[Michael Betancourt]

Pushed the implementation fixes, I'll spend one more day looking for the source of the Mac/Windows discrepancy tomorrow.

> stan-dev+u...@googlegroups.com <mailto:stan-dev%2Bunsu...@googlegroups.com>.

>
> For more options, visit https://groups.google.com/groups/opt_out.
>
>
>
> --
> 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.

> For more options, visit https://groups.google.com/groups/opt_out.
>
>
>
> --
> 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.

[Marcus Brubaker]

I'm not sure it's worth looking too closely for the Mac/Windows discrepancy.

There are a number of reasons that different compilers/optimization
levels/CPUs/platforms/etc would give slightly different floating point
results. E.G., is one using x87 instructions while another is using
SSE instructions at some point? With x87 instructions, some floating
point ops are done at 80bit (long double) precision and the result
gets truncated down to 64bit when they're written to memory or
registers. This can result in slightly different results on different
platforms. Unfortunately, there's not too much we can do about that.

Cheers,
Marcus

[Ben Goodrich]

On Tuesday, May 7, 2013 2:54:14 PM UTC-4, Michael Betancourt wrote:

Pushed the implementation fixes, I'll spend one more day looking for the source of the Mac/Windows discrepancy tomorrow.

What exactly are we looking for at this point? As of d097e74, the test.sh file I was doing before is passing on master.

Ben

[Daniel Lee]

We're actually trying to understand the fundamental differences in code from v1.3.0 and current master. This is manifesting itself as differences in Mac/Windows right now. There's something else going on and if we knew what it was, we'd be ok with it.

We're not actually too concerned with floating point differences. Now, with that said, we can't draw the same first sample in Mac and Windows in master, which we were able to do in v1.3.0. We should either be able to reproduce that behavior or explain why we have differences.

[Peter Li]

On 5/6/13 3:34 PM, Matt Hoffman wrote:
> Ok, I'll take a closer look at the code and double-check. I got a
> little confused by the util.sum_prob, but now I see that it's not
> being used to decide which branch to sample from but only as a measure
> of integrator accuracy. One minor difference, though, is that in the
> old version we only computed the average probability for the most
> recent half of the tree, excluding the states in the half-tree nearest
> the starting point. This is because the first state is always
> acceptable, so including states near the start state biases our
> estimate of the average acceptance probability upwards (though
> probably only slightly). I think the original behavior could be
> brought back by moving
>
> util.n_tree = 0;
> util.sum_prob = 0;
>
> into the while loop.
>
> BTW, can someone remind me how to set off code so that it looks nice
> in these emails? Ben's code snippets always look so much nicer than
> mine...

I just found this out, but there is an add-on in thunderbird called
"Markdown Here" that let's you insert code blocks. It also allows you to
add in LaTeX (haven't figured out how to do it yet though, but there's
another add-on called LaTeXIt! that does this).

[Bob Carpenter]

That's almost exactly what I said during the meeting today :-).

What's a bit worrying is that we used to get identical results.

But the algebra's changed, so no reason to believe we still would.

Daniel's working on a test suite, so we can at least do functional
and speed tests cross-platform. If those work out, I wouldn't be
too worried either (as I also said in today's meeting).

- Bob

[Michael Betancourt]

Yeah, it's not that there is a difference but rather that the difference is new and so significant.

[Bob Carpenter]

This is why we need unit tests! Daniel's working on it.
I can join in when the dust clears from DARPA and
I get the traits cleared up so Peter can move forward.

- Bob

> f7b426b10b0a13d76675057ec3c9aa__cea7d9e0c7
> 1dc4aed16e63c44562ae29dc6f3bc5__17c324794e
> aaf821ba2387a629e911a6a0406185__b70af1a27c
> 3b6353b932af74489570bf7a988642__9159892f61
> d0e8d45df36416290e167b81b92cd5__25505ba160
> 04b8f207985873eb8f2c9e0cf5fd9e__63fa65feb7
> c946070d58fd3646574c6a312d80f2__322d97aff3
> 50f7b32642ebaa3e498360f932ad75__9c10de06f3
> f4f0713cfc6c98751c902d1859989f__ef184df808
> 6daaea6b3bfbf566539aa1bac5ecb3__6136cc3437
> fab3c0841d0755aaf81a48a2c729b0__b3b2b959bf
> 86bd0e446cd34d258384201218c2f2__ea3812e423
> ca85e0e6a8bcffd90ca3047a0e540c__674fc8e09d
> 99429292d821d14b1a1c83ef488570__1e3d409c7f
> e11396f1703a402ea248653ec8e115__7d197166af
> 4385f7ff2af3a8651b3ed3152e8552__ab1e27dc0e
> c61335e3a8dda8a777c14e326ea86e__da1916568e
> 0fe3f288f05d228a98c68c1cb60ce2__805a73ec69
> 6bea89ee6ac669d3f7d9d85eaa0fa8__6411e8c601
> 0df04955917fe8ea69a2d4d80d1918__c24be432b7
> 434b44743e1914093a682614ece759__f35d8cda84
> d22105bc4eb35d124de3a13a3b08e9__775db6913e

> We cannot bisect more!
> bisect run cannot continue any more
> |
>
> |

> goodrich@CYBERPOWERPC:/opt/__stan$ git bisect log
> git bisect start
> # bad: [__951b28226bc16a0c30f45355b1fb99__60c958f4df] calculate average n_eff in rstan/example/demo.R
> git bisect bad 951b28226bc16a0c30f45355b1fb99__60c958f4df
> # good: [__1cde932d15a348df17e9dafce17a79__299b9d8a2d] upgrade rstan verion to 1.2.0
> git bisect good 1cde932d15a348df17e9dafce17a79__299b9d8a2d
> # good: [__bac02df9a2d5d6f6be79b4a69cf5d6__0768961df4] Merge branch 'master' of
> https://github.com/stan-dev/__stan <https://github.com/stan-dev/stan>
> git bisect good bac02df9a2d5d6f6be79b4a69cf5d6__0768961df4
> # good: [__c4a2253bb7159992e5de557343fabf__4f26f2d236] merged conflict in cloglog
> git bisect good c4a2253bb7159992e5de557343fabf__4f26f2d236
> # bad: [__3a78cab458537ab1bc611fd03e1070__7828b4e63b] added multiply_lower_tri_self___transpose to fvar/matrix
> git bisect bad 3a78cab458537ab1bc611fd03e1070__7828b4e63b
> # bad: [__f9ea47fca6731b096bcf9f12886921__1eca598c41] Switched to Stan's Eigen header
> git bisect bad f9ea47fca6731b096bcf9f12886921__1eca598c41
> # skip: [__6ed834b3b0a75db62ebb7cc880fc24__9e72f855ce] Update var/covar to use Welford algorithm, transform to
> scale-free output
> git bisect skip 6ed834b3b0a75db62ebb7cc880fc24__9e72f855ce
> # bad: [__d3fddd6b0bab425508ce2f008689a4__623bd8c772] Refactored default update into base_hamiltonian,

> catch exceptions
> thrown by grad_log_prob

> git bisect bad d3fddd6b0bab425508ce2f008689a4__623bd8c772
> # skip: [__6a87a6033a7fe56562f269f7ee0e55__abefd3079f] Added diagonal and dense samplers, refactored static HMC
> algorithm into base class
> git bisect skip 6a87a6033a7fe56562f269f7ee0e55__abefd3079f
> # bad: [__e930ba046d55859bac4c8dccb85450__287ebe234f] Hopefully NUTS is now working
> git bisect bad e930ba046d55859bac4c8dccb85450__287ebe234f
> # skip: [__e11396f1703a402ea248653ec8e115__7d197166af] Interface/output tweaks
> git bisect skip e11396f1703a402ea248653ec8e115__7d197166af
> # skip: [__1dc4aed16e63c44562ae29dc6f3bc5__17c324794e] Removed deprecated code
> git bisect skip 1dc4aed16e63c44562ae29dc6f3bc5__17c324794e
> # skip: [__434b44743e1914093a682614ece759__f35d8cda84] Refactoring...
> git bisect skip 434b44743e1914093a682614ece759__f35d8cda84
> # skip: [__3b6353b932af74489570bf7a988642__9159892f61] Added base MCMC sampler class, sample object for i/o
> git bisect skip 3b6353b932af74489570bf7a988642__9159892f61
> # skip: [__4385f7ff2af3a8651b3ed3152e8552__ab1e27dc0e] Removed deprecated includes
> git bisect skip 4385f7ff2af3a8651b3ed3152e8552__ab1e27dc0e
> # bad: [__d22105bc4eb35d124de3a13a3b08e9__775db6913e] Updated for mcmc refactor
> git bisect bad d22105bc4eb35d124de3a13a3b08e9__775db6913e
> # skip: [__f7b426b10b0a13d76675057ec3c9aa__cea7d9e0c7] Start testing the new MCMC code
> git bisect skip f7b426b10b0a13d76675057ec3c9aa__cea7d9e0c7
> # skip: [__d0e8d45df36416290e167b81b92cd5__25505ba160] Added apapter class
> git bisect skip d0e8d45df36416290e167b81b92cd5__25505ba160
> # skip: [__ca85e0e6a8bcffd90ca3047a0e540c__674fc8e09d] Fixed bug in random number generator instatiation,

> moved to
> holding only a reference to the external RNG, various const and pass-by-reference fixes

> git bisect skip ca85e0e6a8bcffd90ca3047a0e540c__674fc8e09d
> # skip: [__aaf821ba2387a629e911a6a0406185__b70af1a27c] Separated adapter into hiearchical classes
> git bisect skip aaf821ba2387a629e911a6a0406185__b70af1a27c
> # skip: [__86bd0e446cd34d258384201218c2f2__ea3812e423] Merged master, preserving development code in mcmc
> and gm
> git bisect skip 86bd0e446cd34d258384201218c2f2__ea3812e423
> # skip: [__04b8f207985873eb8f2c9e0cf5fd9e__63fa65feb7] sample() returns acceptance statistic (exact

> statistic will
> depend on algorithm).

> git bisect skip 04b8f207985873eb8f2c9e0cf5fd9e__63fa65feb7
> # skip: [__99429292d821d14b1a1c83ef488570__1e3d409c7f] Development command file
> git bisect skip 99429292d821d14b1a1c83ef488570__1e3d409c7f
> # skip: [__c946070d58fd3646574c6a312d80f2__322d97aff3] Debugging
> git bisect skip c946070d58fd3646574c6a312d80f2__322d97aff3
> # skip: [__fab3c0841d0755aaf81a48a2c729b0__b3b2b959bf] Added debugging model
> git bisect skip fab3c0841d0755aaf81a48a2c729b0__b3b2b959bf
> # skip: [__50f7b32642ebaa3e498360f932ad75__9c10de06f3] Debugging
> git bisect skip 50f7b32642ebaa3e498360f932ad75__9c10de06f3
> # skip: [__0fe3f288f05d228a98c68c1cb60ce2__805a73ec69] Templated out phase space point so that it can also hold

> temporary variables for more sophisticated Hamiltonians

> git bisect skip 0fe3f288f05d228a98c68c1cb60ce2__805a73ec69
> # skip: [__f4f0713cfc6c98751c902d1859989f__ef184df808] Compute gradient only after position has been

> updated, unify
> variable naming to camel case convention.

> git bisect skip f4f0713cfc6c98751c902d1859989f__ef184df808
> # skip: [__c61335e3a8dda8a777c14e326ea86e__da1916568e] Restored bulk of previous content, ready for other
> samplers
> git bisect skip c61335e3a8dda8a777c14e326ea86e__da1916568e
> # skip: [__6daaea6b3bfbf566539aa1bac5ecb3__6136cc3437] Store a reference to the external model, not an
> explicity copy
> git bisect skip 6daaea6b3bfbf566539aa1bac5ecb3__6136cc3437
> # skip: [__6bea89ee6ac669d3f7d9d85eaa0fa8__6411e8c601] Template fixes
> git bisect skip 6bea89ee6ac669d3f7d9d85eaa0fa8__6411e8c601
> # skip: [__0df04955917fe8ea69a2d4d80d1918__c24be432b7] Refactored, separating each class into its own file
> git bisect skip 0df04955917fe8ea69a2d4d80d1918__c24be432b7
> goodrich@CYBERPOWERPC:/opt/__stan$

> |
>
> Ben
>
> --
> 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+unsubscribe@__googlegroups.com <mailto:stan-dev%2Bunsu...@googlegroups.com>
> <mailto:stan-dev%__2Buns...@googlegroups.com <mailto:stan-dev%252Buns...@googlegroups.com>__>.
>
> For more options, visit https://groups.google.com/__groups/opt_out <https://groups.google.com/groups/opt_out>.

>
>
>
> --
> 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+unsubscribe@__googlegroups.com <mailto:stan-dev%2Bunsu...@googlegroups.com>.
> For more options, visit https://groups.google.com/__groups/opt_out <https://groups.google.com/groups/opt_out>.

>
>
>
> --
> 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+unsubscribe@__googlegroups.com <mailto:stan-dev%2Bunsu...@googlegroups.com>.
> For more options, visit https://groups.google.com/__groups/opt_out <https://groups.google.com/groups/opt_out>.

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

[Ben Goodrich]

On Tuesday, May 7, 2013 3:00:51 PM UTC-4, Marcus Brubaker wrote:

I'm not sure it's worth looking too closely for the Mac/Windows discrepancy.

There are a number of reasons that different compilers/optimization
levels/CPUs/platforms/etc would give slightly different floating point
results.  

That is all true, but now I can get identical output with g++-4.6 and clang++-3.3, although there are floating point differences that starting in the first few iterations when cross-compiled with mingw64.

However, it was not entirely true that there were no differences as of v1.3.0. The attached files are identical through iteration 358. On 359, the first parameter differs in the string sense (-6.58535e-05 vs. -6.58535e-005) but that seems to be irrelevant. That happens a few more times. On iteration 748, there is a difference in the first parameter in the 8th decimal place. And this happens again in some subsequent iterations that accumulate and eventually cause genuine differences.

Ben

[Marcus Brubaker]

You need something like 16 digits of printed precision to actually
check for floating point equality of a double with a string. Perhaps
one compounding factor here is that I pushed a change recently which
upped the precision of the CSV output. I think that happened after
1.3.0, but I'm not 100% certain. How are people comparing CSVs?

Cheers,
Marcus

[Ben Goodrich]

On Tuesday, May 7, 2013 7:44:23 PM UTC-4, Marcus Brubaker wrote:

You need something like 16 digits of printed precision to actually
check for floating point equality of a double with a string.  Perhaps
one compounding factor here is that I pushed a change recently which
upped the precision of the CSV output.  I think that happened after
1.3.0, but I'm not 100% certain.  How are people comparing CSVs?

What commit was that? It might be possible to temporarily apply that at the v1.3.0 tag. But with current master the differences show up it <= 8 decimal places within a few dozen iterations whereas they were only visible after a few hundred iterations at v1.3.0.

Ben

[Marcus Brubaker]

This was in commits d88b257a and f06ed309

Cheers,
Marcus

[Ben Goodrich]

On Tuesday, May 7, 2013 8:25:33 PM UTC-4, Marcus Brubaker wrote:

This was in commits d88b257a  and f06ed309

If you do

git checkout v1.3.0
git cherry-pick d88b257a
git cherry-pick f06ed309

and run the test.sh script, there are floating point differences (when cross-compiling) in the 14th or 15th digit after just a few iterations (attached). It seems like the floating point differences are accumulating faster after the refactoring, but I don't know if it should block a release.

Ben

[Bob Carpenter]

In a way, it's reassuring we had differences before,
so bravo for tracking them down, Ben and Marcus!

It would be nice if the arithmetic were more stable,
but I agree with Ben that it shouldn't necessarily
block a release.

What would convince me we're OK is a set of functional tests that
show similar or better performance in terms of mean and variance of
effective sample sizes per second.

I also second Daniel's vote for some more unit testing
at the sampler level. I see code's already being checked
in along those lines! As soon as I crawl out from under
the DARPA proposal, I'll lend a hand.

I also agree with Daniel that it'd be very nice to make
the arithmetic more stable.

Daniel suggested blegging [begging on the blog] for help
from someone who's more familiar with floating point issues.
I'm not too proud to beg, but I'd like to wait until things have
stabilized a bit with more unit tests.

- Bob

[Michael Betancourt]

I went through all of the code again this morning, focusing on tracing differences between 1.3.0 NUTS and master NUTS. Everything looks to be equivalent except for a small difference in the NUTS criterion that can lead to significant differences (not necessarily errors, just differences) in the resulting trees, and hence samples, when we aggressively push the integrator step size.

If we had a perfect integrator then the two algorithms would be equivalent with

rho (new algorithm) = \delta x (old algorithm) \propto \int dt M^{-1} p(t)

When we appeal to the leapfrog integrator, however, these end up being slightly different,

rho = sum_{i} M^{-1} p^{i}

\delta x = sum_{i} M^{-1} p^{i - 1/2}

i.e. the old algorithm aggregates the momenta at the half step whereas the new algorithm uses the momenta after the full step.

When the accuracy of the integrator is limited by a large step size (which may be totally fine for HMC) this difference can be nontrivial, especially when the tree depth is big. Again, this just explains differences in the exact sequence of samples and possibly different sensitivities to floating point arithmetic; both should yield fine inferences.

-Mike