NaN derivative in square(sqrt(x)) in model local variable

[syclik]

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

[syclik]

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)) to u.

I just checked that the boundary conditions in both sqrt() and square() 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

sqrt'(u) = 0.5 / sqrt(u)

square'(v) = 2 * v

yields

d/dx sqrt(square(x-x))

  = sqrt'(square(x-x)) * square'(x - x) * d/dx (x - x)

  = (+ infinity)  * 0 * 0

  = NaN.

I am going to add a discussion of this case to the manual, so I'm linking #

The problem is that (infinity * 0) produces NaN in floating-point arithmetic. There's no way given just the infinity and 0 that mulitplying them should produce 1.

[syclik]

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() for sqrt() to be:

void chain() {
  if (val_ > 1e-300)
    avi_->adj_ += adj_ / (2.0 * val_);
  else
    avi_->adj_ += adj_ / 2.0;
  }

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]

@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]

Great. This wasn't ever fixed. If you want to do as Marcus suggests, it's just going to be the sqrt.hpp in rev and fwd. I still don't get why Marcus wanted to do it that way but I trust he knows what he's doing.

[mbrubake]

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]

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]

I'm closing this because I don't think having square(sqrt(0)) produce NaN is a bug.