Reproducing sklearn Mauna Loa part 2 (not getting good fits with gpytorch, looking at initialization and optimization)

[bcolloran]

This is a follow up to #1646, where @Balandat gave me some very helpful advice that allowed me to fit models successfully. Now that I'm able to fit models successfully, the next step is to get really great fits :-)

As before, I'm using the sklearn Mauna Loa example as a baseline to help me sanity check things as I go. What I'm finding now is that, though I'm able to reproduce sklearn outputs almost identically with gpytorch when I specify equivalent kernels and do no additional training (see this gist), when I initialize identical kernels in sklearn and gpytorch and try then try to optimize them, gpytorch I get substantially worse fits than I get with the .fit method of sklearn.

Since this is from the same initialization and with the same noise level, my next step was to check the optimizer -- I'd been following the standard gpytorch examples that use Adam, but sklearn uses LBFGS by default, so I tried dropping in gpytorch's LBFGS optimzer based on the example here.

Also, the sklearn example offers a pre-optimzed kernel specification as a point of departure, so I tried optimizing starting both from that kernel (which I think should represent a pretty good initialization) and from an unoptimized default kernel (which I specified to have the same initial parameters for both gpytorch and sklearn. Starting from those kernels, I found:

preoptimized kernel

• sklearn gives best fit
  *gpytorch with Adam and LBFGS both improve on the preoptimized fit (LFBGS slightly more so), but not to the same extent as sklearn

default kernel

• sklearn is able to produce a fit starting from the default kernel that is basically as good as when starting from the pre-optimized kernel
• Gpytorch with Adam produces a pretty poor fit
• Gpytorch with LBFGS crashes

I have a notebook at this gist that shows what I've tried. (Note that to run it you'll need to adjust the sys.path.append('../../gpytorch/examples') call in the imports cell to point to the gyptorch/examples folder on your machine.)

Many thanks for any advice anyone may have of what to try next, or pointers about things I'm doing wrong or have overlooked! :-)

[bcolloran]

Any ideas from anyone? I bet I'm missing something that would be obvious to someone who has worked with Gpytorch a lot, but that has not occurred to me to try yet. Suggestions would be most appreciated! Thanks!

[bcolloran]

Does anyone have any suggestions of other things I might try? I have the feeling that this should be trivially easy (I can't think of any reason why sklearn should give me better fits than Gpytorch) but I remain stuck on this and I'm not sure what I'm missing. Thanks in advance!

[wjmaddox]

Pretty sure that the issue is one of optimization: I substituted in @Balandat 's botorch.fit_gpytorch_model utility (which uses multiple random restarts of scipy's L-BFGS-B internally) and it gives reasonable fits:

That's basically a two line change to your optimize_LBFGS function to be just

    from botorch import fit_gpytorch_model
    fit_gpytorch_model(mll)
   return None

This makes me suspect that scipy's L-BFGS-B outperforms the pytorch version you're using (which isn't that surprising).

[Balandat]

which uses multiple random restarts of scipy's L-BFGS-B internally

Just a quick clarification that we are by default not actually using multiple random restarts (we fall back if there are issues with the optimization), so this is probably just the difference between scipy's L-BFGS-B and the implementation you're using (not sure where that one is coming from, btw).

[wjmaddox]

Looks like they're using the Pytorch-LBFGS repo (https://github.com/hjmshi/PyTorch-LBFGS), but could be incorrect.

[Balandat]

Hmm I see. That repo implements a stochastic quasi-second order method, but the MLL here is deterministic. That's likely why the scipy variant gets much better fits.

[bcolloran]

Thanks @Balandat! I was not aware of botorch, and I suspect that even if I'd known of it, never in a million years would I have thought to try swapping out LBFGS implementations, so thank you for your expertise and insight! :-) Using botorch.fit_gpytorch_model was indeed the missing piece.

---

In case it may be useful to anyone finding this thread in the future here's a complete notebook that uses gpytorch to more or less reproduce the sklearn Mauna Loa example.

[Balandat]

FWIW We're currently working on a custom SQP / L-BFGS-B optimizer that is able to exploit parallelism/batched evaluation (scipy's cannot do that). That solver could also be of interest for fitting models using a bunch of random restarts (currently the focus is on acquisition function optimization). I think what we'd need for this is modify the MLL computation so that it can happen in a batched fashion. - @vishwakftw, @dme65

[AustinT]

Just having a glance at your code, a key difference might be that botorch turns off fast computations and uses exact cholesky, whereas it looks like you were using the default CG method. If you wrapped your optimization in a with gpytorch.settings.fast_computation(False, False, False) I wouldn't be surprised if that made the difference (L-BFGS apparently doesn't work well with noisy estimates of MLL).

[wjmaddox]

I don't think that this is the case as the number of data points (~500) is less than the 800 data point threshold that's used for CG. Quickly re-running the Adam fit with that flag also doesn't improve the fit.

Yes, the noisy estimates of the MLL could also be playing a role in the fit - here the log determinant calculations? rather than the CG calls.

[jacobrgardner]

Since the data isn't normalized, I think you're just running into numerical issues with gigantic lengthscales (100+) and extremely small likelihood noises. For me, either just training for 75 iterations of LBFGS by passing training_iter=75 works to avoid the noise getting so small that you run into fp32 issues, which might explain why the fit_gpytorch_model function works since it presumably stops based on a stopping criterion instead of continuing to run.

Alternatively, you can add torch.set_default_dtype(torch.float64) to the top of your notebook, and then just running everything as in the original notebook works fine for me (note that sklearn / scipy will use fp64 by default, whereas torch uses fp32 by default).

I get the below: