Repository navigation
NaN derivative in square(sqrt(x)) in model local variable #21
Description
From @bob-carpenter on September 4, 2013 7:7
Simon Barthelmé reported via the stan-users list. I simplified the reported model to this minimal instance:
parameters {
real x;
}
model {
real y;
y <- sqrt(square(x-x));
x ~ normal(0,1);
}
As Simon correctly diagnosed, this model causes a problem with autodiff, as can be seen when running from RStan 1.3:
TESTING GRADIENT FOR MODEL 'barthelme2' NOW (CHAIN 4).
TEST GRADIENT MODE
Log probability=-0.0571533
param idx value model finite diff error
0 -0.338093 nan 0.338093 nan
First step is to move this into a unit test that displays the behavior.
Copied from original issue: stan-dev/stan#203
Metadata
Metadata
Assignees
Labels
No labels
Activity
syclik commentedon Jul 6, 2015
From @bob-carpenter on September 23, 2013 18:39
We can't really fix this problem without a huge amount of work moving to expression templates; with expression templates, we could follow Eigen in carrying out symbolic math at compile time, such as reducing
sqrt(square(u))tou.I just checked that the boundary conditions in both
sqrt()andsquare()are correct, as is the application of the chain rule.The problem is that the chain rule is applied literally with floating-point arithmetic, which given
yields
I am going to add a discussion of this case to the manual, so I'm linking #
The problem is that
(infinity * 0)producesNaNin floating-point arithmetic. There's no way given just the infinity and 0 that mulitplying them should produce 1.syclik commentedon Jul 6, 2015
From @bob-carpenter on November 20, 2013 5:20
Implement as Marcus suggested in e-mail:
Putting those two things together, maybe we can just rewrite
chain()forsqrt()to be:I asked: Will this do the right thing for HMC and BFGS? And Marcus responded:
This is effectively returning a subgradient for points where the gradient isn't defined. HMC won't care so long as same subgradient is always returned (which it is in this case). The regions where this will happen are extremely small, so, except for a few pathological cases, I wouldn't expect it to impact either method significantly.
sakrejda commentedon Mar 4, 2016
@bob-carpenter Am I right that this change did not make it into the appropriate file? If it's not on your short list I can go through it for the sake of cleaning up issues. If you have a guess at which additional files the change needs to touch that would help too.
bob-carpenter commentedon Mar 4, 2016
Great. This wasn't ever fixed. If you want to do as Marcus suggests, it's just going to be the
sqrt.hppinrevandfwd. I still don't get why Marcus wanted to do it that way but I trust he knows what he's doing.mbrubake commentedon Mar 4, 2016
What I suggested is something of a kludge. I'm not sure if it's really the best idea but it might resolve the issue in practice. That said, I really don't remember too much of the discussion around this. To put it another way, I think I might mistrust my past self.
sakrejda commentedon Mar 5, 2016
Ok, I'll figure out what to do here, thanks for the pointers. It sounds like it needs a test and some demonstration that whatever we do doesn't screw up results as far as sampling is concerned. Seems like a gentle enough introduction to the autodiff code base...
bob-carpenter commentedon Aug 17, 2016
I'm closing this because I don't think having
square(sqrt(0))produceNaNis a bug.