Trouble reproducing sklearn Mauna Loa

[bcolloran]

Hi folks, I'm struggling to reproduce the sklearn Mauna Loa model, and I'd love any pointers folks here might have. In particular, I'm having trouble fitting the periodic kernel component.

You can see the complete notebook here. It's based on the Simple_GP_regression.ipynb tutorial notebook -- I've tried to stick close to that core to keep things simple and (hopefully) remove extraneous errors.

It seems probable that the place I'm going astray is my translation of the kernel from the sklearn api to Gpytorch. In sklearn, the kernel is:

# Kernel with optimized parameters
k1 = 50.0**2 * RBF(length_scale=50.0)  # long term smooth rising trend
k2 = 2.0**2 * RBF(length_scale=100.0) \
    * ExpSineSquared(length_scale=1.0, periodicity=1.0,
                     periodicity_bounds="fixed")  # seasonal component
# medium term irregularities
k3 = 0.5**2 * RationalQuadratic(length_scale=1.0, alpha=1.0)
k4 = 0.1**2 * RBF(length_scale=0.1) \
    + WhiteKernel(noise_level=0.1**2,
                  noise_level_bounds=(1e-5, np.inf))  # noise terms

kernel = k1 + k2 + k3 + k4

Which I've translated as:

#### long term
k_long_term_RBF = RBFKernel()
k_long_term_RBF.lengthscale = 50.0

k_long_term = ScaleKernel(k_long_term_RBF)
k_long_term.outputscale = 50.0**2

#### seasonal
k_seasonal_RBF = RBFKernel()
k_seasonal_RBF.lengthscale = 90

k_seasonal_periodic = PeriodicKernel()
k_seasonal_periodic.period_length = 1.0
# this must be squared, see https://github.com/cornellius-gp/gpytorch/issues/1020
k_seasonal_periodic.lengthscale = 1.0**2

k_seasonal = ScaleKernel(k_seasonal_RBF * k_seasonal_periodic)
k_seasonal.outputscale = 2.0**2

### medium term
k_medium_RQ = RQKernel()
k_medium_RQ.alpha = 1.0
k_medium_RQ.lengthscale = 1.0

k_medium = ScaleKernel(k_medium_RQ)
k_medium.outputscale = 0.5**2

#### noise term
k_noise_RBF = RBFKernel()
k_noise_RBF.lengthscale = 0.1

k_noise = ScaleKernel(k_noise_RBF)
k_noise.outputscale = 0.1**2

sklearn_reproduction_kernel = (
    k_long_term
    + k_seasonal
    + k_medium
    + k_noise
)

I'm also not sure how to reproduce the WhiteKernel(noise_level=0.1**2, noise_level_bounds=(1e-5, np.inf)) component from sklearn. Based on a pointer in this issue, I've tried likelihood.initialize(noise=0.1**2), but setting the noise level that low causes "not positive definite" crashes with this kernel.

Thanks for any pointers you can give me. I've tried to do as much googling and looking through examples as I can, but at this point I'm getting diminishing returns, so any ideas are much appreciated!

(PS: if I can get this working I'd be happy to submit a PR with a demo notebook on this data. I suspect it might help other folks since it's kind of a classic test case and provides a direct comparison with a scikit learn demo.)

[Balandat]

Have you tried running this in double rather than float dtype? The default float will often cause numerical issues on ill-conditioned systems that often appear in GPs.

[bcolloran]

Wow thanks! So simple, but I'm sure it would have taken me ages to check that on my own! Very grateful for the help :-)

For posterity and any who may need it later, here's an updated working notebook:
https://gist.github.com/bcolloran/8869c54f265606b66c09a56eb1241631
Your pointer along with a couple tweaks on the sklearn side of things, and I'm getting the expected alignment between gpytorch and sklearn output for pre-trained models. Thanks again!