Gaussian Process as default setting

[taikichi-pc]

Thank you for your development.
I have a question about below.

1. What is difference between fit_gpytorch_model and fit_gpytorch_mll

2. When I use fit_gpytorch_model, I can use argument, "optimizer".
   When I set "optimizer = fit_gpytorch_torch", optimization iteration was improved (decreased).
   What is difference between fit_gpytorch_model(mll, optimizer = fit_gpytorch_torch) and fit_gpytorch_model(mll)

3. I would like to implement the Cholesky decomposition, How should I write code to it?
   　Does Botorch use a less computable Gaussian Process Regression than that?

I'm sorry for my elementary questions.　
Please let me know about details and document link.

Thank you for your kindness.
Best regards.

[saitcakmak]

What is difference between fit_gpytorch_model and fit_gpytorch_mll

fit_gpytorch_model is the deprecated API for fitting models. It has been replaced with fit_gpytorch_mll but we're keeping it around for a while for backwards compatibility reasons. It calls fit_gpytorch_mll under the hood.

What is difference between fit_gpytorch_model(mll, optimizer = fit_gpytorch_torch) and fit_gpytorch_model(mll)

You can pass an optimizer to fit_gpytorch_mll as well. We generally use scipy optimizers for fitting the Exact GP models, and use torch optimizers for approximate inference.

I would like to implement the Cholesky decomposition, How should I write code to it?
　Does Botorch use a less computable Gaussian Process Regression than that?

What is the need for implementing your own Cholesky decomposition? BoTorch uses Cholesky decomposition during model fitting and for sampling from the posteriors. We use the implementation from PyTorch with some wrappers from GPyTorch to handle possible non-PSD matrices.

[taikichi-pc]

Thank you for your quick reply.

I'd like to know more details about following.

You can pass an optimizer to fit_gpytorch_mll as well. We generally use scipy optimizers for fitting the Exact GP models, and use torch optimizers for approximate inference.

When we use "fit_gpytorch_model(mll)" (optimizer = None),
what kind of optimizer would be used?

And dose "optimizer = fit_gpytorch_torch" more precisely optimize GP-model's parameter more than "optimizer = None" in fit_gpytorch_model?

In short, "optimizer = fit_gpytorch_torch" optimize more well?

[taikichi-pc]

I think it may depend on data,
For example, when using "optimizer = fit_gpytorch_torch",
is it possible that the parameters of the Gaussian process are better fitted and the variance is smaller?

[saitcakmak]

When we use "fit_gpytorch_model(mll)" (optimizer = None), what kind of optimizer would be used?

For most models it will use fit_gpytorch_scipy / fit_gpytorch_mll_scipy (first one being deprecated). If you're using an ApproximateGP with many training points, then it'll use fit_gpytorch_torch / fit_gpytorch_mll_torch.

In short, "optimizer = fit_gpytorch_torch" optimize more well?

In most cases, no. It has its uses in training ApproximateGP models or some other over-parameterized models such as HigherOrderGP, but for common GP models used in BO, we prefer the scipy optimizers.

is it possible that the parameters of the Gaussian process are better fitted and the variance is smaller?

To answer this question, we would have to define "better". By default, we use MAP estimation, and maximize the corresponding log likelihood (+ prior terms). Scipy optimizers reliably find the maximizer in most cases. It is possible that you'd consider a set of hyper-parameters that do no maximize the likelihood "better", in which case this is not a question of the optimizer choice but the choice of the priors & the likelihood or more generally the method used to pick the hyper parameters.

[taikichi-pc]

Thank you a lot for your kindness and detail explanation.
And I have some questions.

I am using SingleTaskGP at 2 patterns bellow.
1

def initialize_model(train_x, train_obj):
    gp_models = []
    for i in range(train_y.shape[-1]):
        gp_models.append(
            SingleTaskGP(train_x, train_y[..., i:i+1], outcome_transform=Standardize(m=1)))
    gp_model = ModelListGP(*gp_models)
    mll = SumMarginalLogLikelihood(gp_model.likelihood, gp_model)
    return mll, gp_model

2

def initialize_model_light_calc(
    train_x,
    train_obj,
    use_keops = False,
):
    gp_models = []
    kernel_kwargs = {"nu": 2.5, "ard_num_dims": train_x.shape[-1]}
    base_kernel = (
        KMaternKernel(**kernel_kwargs) if use_keops else MaternKernel(**kernel_kwargs)
    )
    covar_module = ScaleKernel(base_kernel)
    likelihood = GaussianLikelihood(noise_constraint=Interval(1e-8, 1e-3))

    for i in range(train_y.shape[-1]):
        gp_models.append(
            SingleTaskGP(train_x, train_y[..., i:i+1], outcome_transform=Standardize(m=1), likelihood=likelihood, covar_module=covar_module))
    gp_model = ModelListGP(*gp_models)
    mll = SumMarginalLogLikelihood(gp_model.likelihood, gp_model)
    return mll, gp_model

When first one used, mll is like below. (only remarkable part)
1

      (0): SingleTaskGP(
        (likelihood): GaussianLikelihood(
          (noise_covar): HomoskedasticNoise(
            (noise_prior): GammaPrior()
            (raw_noise_constraint): GreaterThan(1.000E-04)
          )
        )
        (mean_module): ConstantMean()
        (covar_module): ScaleKernel(
          (base_kernel): MaternKernel(
            (lengthscale_prior): GammaPrior()
            (raw_lengthscale_constraint): Positive()
          )
          (outputscale_prior): GammaPrior()
          (raw_outputscale_constraint): Positive()
        )
        (outcome_transform): Standardize()

When second one used, mll is like below. (only remarkable part)
2

      (0): SingleTaskGP(
        (likelihood): GaussianLikelihood(
          (noise_covar): HomoskedasticNoise(
            (raw_noise_constraint): Interval(1.000E-08, 1.000E-03)
          )
        )
        (mean_module): ConstantMean()
        (covar_module): ScaleKernel(
          (base_kernel): MaternKernel(
            (raw_lengthscale_constraint): Positive()
          )
          (raw_outputscale_constraint): Positive()
        )
        (outcome_transform): Standardize()
      )

According to above difference between 1 and 2,
These are different point.

• (raw_noise_constraint): GreaterThan(1.000E-04) or Interval(1.000E-08, 1.000E-03)
• (lengthscale_prior): GammaPrior() or nothing
• (outputscale_prior): GammaPrior() or nothing

What is different do you think?
Q1 If lengthscale_prior and outputscale_prior is nothing (like 2nd one), what kind of prior is set?
Q2 When mll fitted, What is (distance_module)?
Q3 Why did you choose gamma prior to lengthscale and outputscale? (briefly ok)
Q4 (raw_noise_constraint) is controling about GPs variance in object space?

Actually, when 2nd one and fit_gpytorch_model(mll, optimizer = fit_gpytorch_torch) are used, next point suggestion is better and more efficiently converged than 1st one.
(and 1st one used with fit_gpytorch_model(mll, optimizer = fit_gpytorch_torch) is worse suggestion.)

And 2nd one used with fit_gpytorch_model(mll) without fit_gpytorch_model(mll, optimizer = fit_gpytorch_torch)
one error happened like following.
RuntimeError: 60 elements of the 60 element gradient array gradf are NaN. This often indicates numerical issues. 

Q5 Why 2nd one(used with fit_gpytorch_model(mll, optimizer = fit_gpytorch_torch)) is better(good next point) do you think?
Q6 Why that only used 2nd one with fit_gpytorch_model(mll) happened error like that?

I'm sorry for my complex questions and my poor knowledge.
I appreciate your cooperations.

Thank you.

[taikichi-pc]

@saitcakmak
Could you answer this?
If you need more detail for answer, Please let me know.

[saitcakmak]

Q1 If lengthscale_prior and outputscale_prior is nothing (like 2nd one), what kind of prior is set?

No prior is set. You can think of this as being equivalent to using a prior that assigns equal value to all possible values.

Q3 Why did you choose gamma prior to lengthscale and outputscale? (briefly ok)

I wasn't involved in choosing the default priors. They were chosen based on prediction accuracy on a set of example datasets. The default priors are actually not great if you have many parameters. We are planning on updating them in near future. cc @dme65

Q4 (raw_noise_constraint) is controling about GPs variance in object space?

It constrains the inferred output noise to be in the given range (e.g., noise should be non-negative). I am not too familiar with priors / constraints in GPyTorch.

And 2nd one used with fit_gpytorch_model(mll) without fit_gpytorch_model(mll, optimizer = fit_gpytorch_torch)
one error happened like following.
RuntimeError: 60 elements of the 60 element gradient array gradf are NaN. This often indicates numerical issues.

This means your covariance matrix was most likely (numerically) singular. It can happen if you have repeated observations / observations that are too close in the search space, or if the covar_module was somehow misconfigured, leading to singular covariance matrices.

Q5 Why 2nd one(used with fit_gpytorch_model(mll, optimizer = fit_gpytorch_torch)) is better(good next point) do you think?

I don't think one method is better than another necessarily. It will be problem dependent, as explained above.

Q6 Why that only used 2nd one with fit_gpytorch_model(mll) happened error like that?

The optimizer likely lead to the hyper parameters to converge to some values that produced the singular covariance matrices. It is hard to say why it would happen with one optimizer and not the other.

[taikichi-pc]

@saitcakmak -san
Thank you for your reply and details in spite of my abstract question.

I'd like to ask more about
RuntimeError: 60 elements of the 60 element gradient array gradf are NaN. This often indicates numerical issues.

This means your covariance matrix was most likely (numerically) singular. It can happen if you have repeated observations / observations that are too close in the search space, or if the covar_module was somehow misconfigured, leading to singular covariance matrices.

Q1. In my code, it was probably due to repeated observations/observations too close together in the search space.
To avoid this error, what code should (could) be added or modified?
Maybe my current code is optimized to local solution.

Q2. You suggest that "the covar_module was somehow misconfigured"
Generally, what does misconfigured mean in this situation?
(surely, when I don't use Gamma prior for lengthscale&outputscale prior (use no_prior), this error often happen. And I avoided it to use fit_gpytorch_torch)
Maybe, Because fit_gpytorch_torch is supposed to use ApproximateGP, the gradient is unlikely to be NaN..?

[saitcakmak]

Q1. In my code, it was probably due to repeated observations/observations too close together in the search space. To avoid this error, what code should (could) be added or modified?

You could check each generated candidate against previously evaluated points and decide on how to proceed. If it's an identical / near identical point, you could manually perturb it a bit to ensure a bit more exploration. You could also replace it with a random point or use a different acquisition function / model to see if you get a different point. You can pretty much use any heuristic you come up with here.

Q2. You suggest that "the covar_module was somehow misconfigured". Generally, what does misconfigured mean in this situation?

It could be that some of the hyper parameters converged to some very large or very small values, which can lead to numerically unstable covariance matrices. This is where the prior & constraints on the length scales come into play. They bias the hyper parameters into certain ranges, which tends to lead to more stable results.

As far as I know, fit_gpytorch_torch should work with priors on the length scales (or ExactGP models in general). I don't know why you'd be getting errors with the Gamma prior there.