Laplace Bug when passing Eigen::Map in tuple of functor arguments

[SteveBronder]

From @avehtari on discourse

I'm trying to simplify the result here but have not been successful yet. i.e. this gives correct values. @WardBrian can you see anything here that is meaningfully different than what the compiler generates?

https://github.com/stan-dev/math/blob/e53e95b7fcdfad20e1d0e7df24333e4d5a8b5510/test/unit/math/laplace/aki_ex_test.cpp

[WardBrian]

@SteveBronder I think you misread Aki's post, the data you've hardcoded in here is what works.

The other data, like datat <- list(N=1, y=183, mu=0.5, sigmaz=2.5, run_laplace=0), causes problems

[SteveBronder]

Oh I did misread it my apologies

[WardBrian]

Changing the data to

  const std::vector<int>    y{183};
  const std::vector<double> mu{0.5};
  const double sigmaz = 2.5;

does indeed lead to the test failing with an exception thrown

C++ exception with description "Error in laplace_marginal_density: theta contains NaN values for all values." thrown in the test body.

[WardBrian]

By adding some prints to the likelihood it seems the issue here is that theta is shooting off to massive values in the optimization:

[0]
[0]
[100.265]
[100.265]
[100.265]
[4.78573e+30]
[4.78573e+30]
[4.78573e+30]

Exception: Error in laplace_marginal_density: theta contains NaN values for all values.

Inside integrand, the value of theta never exceeds ~7000 it seems

If you play around with theta_0, you can find values that work (for example, { "N": 1, "y": [183], "mu": [0.5], "sigmaz": 2.5, "run_laplace": 1, "theta_0": [3] }) seems fine! Theta never gets larger than 7.5 ), but most of them shoot off into either massively positive or massively negative numbers and then you get the exception

cc @avehtari

[avehtari]

I added print also for the log likelihood. Running with

datat <- list(N=1, y=183, mu=0.5, sigmaz=2.5, run_laplace=1)

interestingly theta does come back from e+31 range

Chain 1 [0] 
Chain 1 -684.009 
Chain 1 [0] 
Chain 1 -684.009 
Chain 1 [100.265] 
Chain 1 -5.77624e+43 
Chain 1 [100.265] 
Chain 1 -5.77624e+43 
Chain 1 [100.265] 
Chain 1 -5.77624e+43 
Chain 1 [-1.58456e+31] 
Chain 1 -2.89975e+33 
Chain 1 [-1.58456e+31] 
Chain 1 -2.89975e+33 
Chain 1 [-1.58456e+31] 
Chain 1 -2.89975e+33 
Chain 1 [1143.75] 
Chain 1 -inf 
Chain 1 [1143.75] 
Chain 1 -inf 
Chain 1 [1143.75] 
Chain 1 -inf 

One thing the optimization should do is that when ever the log likelihood is evaluated to be -inf or nan, the step should be halved. I'll try to have a look at the current code.

[avehtari]

I had a quick look at the code in laplace_marginal_density.hpp. It could be made safer by checking that every time when laplace_likelihood::log_likelihood() or laplace_likelihood::diff() are called, the return value is finite (not -inf or nan) and if a non-finite return value is encountered then the step size is halved. Now it is possible to end up with theta which has non-finite log likelihood and/or derivative.

[WardBrian]

the step should be halved.

Reading the current code, it looks like the stepsize is always 0.5 -- is that what you were referring to @avehtari?

There used to be a comment from Charles about changing this but it looks like it got deleted at some point

math/stan/math/mix/functor/laplace_marginal_density.hpp

Line 219 in a23525a

= (a + a_prev) * 0.5; // TODO(Charles) -- generalize for any factor

[SteveBronder]

@avehtari @WardBrian I have a branch here showing what is happening. From the debugging log below we can see that with a mean of 0.5 the gradient for theta shoots off to 180'ish which then causes a and b for updating theta to become huge. That's an issue because in the poisson_log function we take an exponential of that value.

__________
ITER: 0 / 100
diff theta gradient: 181.351
curr theta: 
	0
W:  
	Nonzero entries:
	(1.64872,0) 
	Outer pointers:
	0  $
	1.64872

W_r:  
	1.28403
theta_grad:  
	181.351
a: 
	16.0424
b: 
	181.351
--------------
theta: 
	100.265
--------------
objective_old: 
	-1.79769e+308
objective_new: 
	-5.77624e+43
__________
ITER: 1 / 100
diff theta gradient: 
	-5.77624e+43
curr theta: 
	100.265
W:  
Nonzero entries:
(5.77624e+43,0) 
Outer pointers:
0  $
5.77624e+43

W_r:  
	7.60016e+21
theta_grad:  
	-5.77624e+43
a: 
	-2.5353e+30
b: 
	5.73378e+45

This ran went of to infinity in the objective, but other runs would just go off to some huge but wrong number.

I changed laplace here so that we take a partial step in the direction of the new a instead of a full step. We try a new partial step and if that fails to give us a better theta we halve the step size and keep trying until we to a step size of 1e-8 or find a better theta. This seems to be better and the new version can solve this problem with a step size of 0.9. @avehtari do you know of a better way to fix this?

If we do this step approach do you both think we should expose step_size to the user in tol so they can control it?

[WardBrian]

When you say "halve the step size and keep trying", from where do you keep going? The original point, or the point immediately before failure?

[SteveBronder]

The point immediately before failure. So same a and b, but a smaller step size for the update.

[WardBrian]

I would expect it to perform better if it was trying again from the original point, because then you're not having to compensate for the initial too-large-but-not-yet-failing steps. Is that hard to try?