NUTS step, when selected explicitly, requires the second derivative

The NUTS sampler, when auto-assigned, requires only the first derivative (i.e. calls grad() method only for the operation itself).

For example, this code for Interpolated distribution works:

import numpy as np
import pymc3 as pm

with pm.Model():
    uniform = pm.Interpolated('uniform', np.linspace(0, 1, 100), np.ones(100))
    pm.sample(1000)

However, when I try to specify NUTS or HamiltonianMC step explicitly, it fails because Interpolated distribution doesn't provide second-order derivatives (the variable returned by grad() doesn't have grad() method):

import numpy as np
import pymc3 as pm

with pm.Model():
    uniform = pm.Interpolated('uniform', np.linspace(0, 1, 100), np.ones(100))
    step = pm.NUTS()
    pm.sample(1000, step=step) # fails

According to the NUTS paper, it should require only first-order derivatives, but maybe I miss something.

So my question is: is it a problem or is it for some reason the expected behavior? In the second case I can add a second order gradient implementation to theInterpolated distribution that for example always returns zeros, but I don't understand why is it needed in the first place.

[ColCarroll]

It also fails with Metropolis -- looks like there is some scaling that is getting set as soon as you initialize the step function (the line above # fails).

@ColCarroll You are right, it actually fails at the initialization of the step function.

But Metropolis works for me:

import numpy as np
import pymc3 as pm

with pm.Model():
    uniform = pm.Interpolated('uniform', np.linspace(0, 1, 100), np.ones(100))
    step = pm.Metropolis()
    pm.sample(1000, step=step)

doesn't produce any errors. I run it with Python 3 from the latest version of Anaconda3 and master branch of PyMC3.

So I found that the error can be suppressed if I implement thegrad() method of the gradient operation that throws NotImplementedError.

In this case the NotImplementedError is caught in guess_scaling here, so that it falls back to fixed_hessian and samples.

[aseyboldt]

This is because nuts needs a mass matrix (the scaling argument). If this is not supplied by the user it estimates this by looking at the diagonal of the hessian at the initial position. You should usually specify it explicitly (and also set is_cov accordingly) if you assign the step method yourself.
Not sure I like this behaviour by the way. I don't think there really is a good default, but maybe the identity matrix would be better. The hessian makes sense for a model with an approximate gaussian posterior if the initial position is the map, but that's about it. It doesn't even have to be positive definite in general.
CC @twiecki @fonnesbeck

[twiecki]

@aseyboldt Yeah, the hessian seems to cause more problems than it's worth so I agree that we should default to an identity that we adapt during tuning using the trace.

[fonnesbeck]

Yes, that makes sense.

The referenced above PR #2211 solved the problem by throwing the correct type of exception for not implemented second derivatives and allowing guess_scaling to fall correctly back to fixed_hessian.

However, @aseyboldt pointed out that it is a good idea to replace guess_scaling by just a fixed identity matrix. I opened a PR #2232 for it.