How to pick length-scale bounds for RBC kernels in Gaussian Process Regression? I am trying to fit GP regression models to several thousand $x, y$ pairs independantly. I am using Python's sklearn implementation with a constant kernel plus an RBF kernel plus a white noise kernel. Usually it goes well and I get good results (red is the function, blue the GPR predictions):

However sometimes it doesn't work. I get an error message saying the L-BFGS-B optimizer terminated in an abnormal state. The problem seems to be the bounds of the length scale in the RBF kernel. In the image above it was between $10^{-1}$ and $10^2$. If I change the upper bound to $10^3$ the optimizer fails. If I change it to $10^4$ it works but now the function looks like this:

How do I automatically select an appropriate set of bounds?
 A: When "fitting" a GP with a certain kernel to a dataset, what we're essentially trying to do is to infer the most suitable values of the hyperparameters (of which the lengthscale is an important one), given the data. This usually takes the form of an optimisation problem, where we try to maximise the probability of observing the data, given the model (or, more precisely, maximise the log marginal likelihood).
As GPs are very flexible models, this optimisation can be quite tricky if it is not set up properly, and this includes the viable search space. I believe the behaviour you're observing is due to the optimiser (L-BFGS-B) being initialised in (relatively) flat areas of the search space, and thus getting "stuck" there. By increasing the bounds of the search space, you are making it more likely that the lengthscale is initialised at a large value, which we can anticipate causes quite a flat log marginal likelihood function (a GP with an RBF kernel will "encode" asymptotically linear functions as $l \rightarrow \infty$). I'm not sure why the optimiser isn't failing in the case of the $10^3$ lengthscale upper bound, but you can see the predictions are flat, which suggests it has returned a nonsensically large lengthscale.
I would recommend reading section 3 of this article by Michael Betancourt to understand this problem in a little more detail, but a quick explanation is that suitable lengthscale bounds depend on the support of your function. Though it does depend on the parametrisation (I believe Sklearn uses a fairly standard one), the RBF kernel is unable to "distinguish properly between" (in the context of changes to the opimisation objective) lengthscales shorter than the minimum spacing in your data, and larger than the maximum spacing of your data.
As a starting point, I would recommend placing an upper bound on $l$ equal to the maximum spacing between any two points in your dataset, and a lower bound equal to the larger of either the minimum spacing between any two points or some small-ish value (say 1e-3; to prevent the minimum spacing being zero in the case of duplicate observations).
Note that I don't think the above holds precisely for kernels with anisotropic lengthscales (that is, "ARD" kernels with multiple different lengthscales), but is nevertheless a useful starting point.
