Optimized logpt computation

[junpenglao]

I change the logpt computation in #2949 to fix #2948, however, it slows down the speed as some graph optimization is turned off (those optimization is originally cause the error in #2948). I am trying with a differen approach here.
@ColCarroll

[ColCarroll]

I don't quite understand the problem: it looks like you are trying a more explicit loop to give the compiler a better chance of optimizing? I guess I also don't understand what caused the edge case in the first place, and whether it is ok to let that fail (that is the failing test)

[junpenglao]

Thanks for having a look @ColCarroll.
So the error is a bit difficult to trace down as it happens (i suspect) in the graph optimization. In that particular edge case (edge because it only doest work with sd being a vector in logNormal distribution, other distribution works fine), if you set the ValueGradientFunction to compile in python (mode='FAST_COMPILE') it works fine. My initial fix is to set the logpt computation to the older way, but seems in effect it skip the optimization in the step of the logpt computation step (that's why no more error in the edge case). So now what I am trying to do is be more explicit at these similar edge case and compile a function that turns off the C optimization.
What do you think?

[twiecki]

Seems like a theano bug, no?

[junpenglao]

@twiecki not sure, I cannot really pin down where exactly is the error yet... compiling a theano function using pm.Lognormal.dist().logp()+... etc works without problem

[ColCarroll]

Ran the asv benchmarks locally, and this gets the hierarchical model running quicker again (160 ES/S vs 138 ES/S). Agree that it would be nice to understand more about why this happened, but I wanted 14% of my time back!

[junpenglao]

So, this PR makes me quite paranoid so I did a bit more digging.
Turns out doing logp = tt.add(*map(tt.sum, factors)) is actually slightly faster in terms of raw speed, as the graph is smaller, but surprisingly in practice (I tested a few other models that I am working on) doing logp = tt.sum([tt.sum(factor) for factor in factors]) seems to be numerically more stable (hence the higher effective sample size / sampling time ratio in the benchmark).

All the output is in the folder below if you would like to test on your own model https://github.com/junpenglao/Planet_Sakaar_Data_Science/tree/master/Miscellaneous/theano_graph_optim