[Bug] Equality constraint violated in optimization of ACQF

[jduerholt]

🐛 Bug

In some cases the optimizer within optimize_acqf spits out candidates that violate the defined equality constraints provided to optimize_acqf. I build an MWE for it. You can find it below. If one reduces the number of points provided to the acqf, the chances that the constraint is fulfilled increase. Perhaps you have to run the MWE several times to observe the mentioned behavior, because the constraint is sometimes fulfilled and sometimes not. I attach the example data (totally made up data) as an csv file.

To reproduce

import torch
from torch import tensor
from botorch.models import SingleTaskGP
from botorch.fit import fit_gpytorch_model
from gpytorch.mlls import ExactMarginalLogLikelihood
from botorch.models.transforms.input import Normalize
from botorch.models.transforms.outcome import Standardize
from botorch.optim import optimize_acqf
from botorch.models import ModelListGP
from botorch.acquisition.multi_objective import qNoisyExpectedHypervolumeImprovement
import pandas as pd

torch.manual_seed(1)

tkwargs = {'dtype': torch.float64, 'device': 'cpu'}

df_data = pd.read_csv("equality_constraint.csv")
#df_data = pd.read_csv("equality_constraint.csv").loc[:10]
train_X = torch.from_numpy(df_data[["a", "b", "c", "d"]].values).to(**tkwargs)
train_Y = torch.from_numpy(df_data[["alpha", "beta"]].values).to(**tkwargs)

lower = tensor([0.1,0.3,0.1,30.])
upper = tensor([0.6,0.7,0.7,70.])

bounds = torch.stack((lower, upper)).to(**tkwargs)

models = []
for i, feat in enumerate(["alpha", "beta"]):
    gp = SingleTaskGP(train_X, train_Y[:,i].unsqueeze(-1), input_transform=Normalize(d=4), outcome_transform=Standardize(m=1))
    mll = ExactMarginalLogLikelihood(gp.likelihood, gp)
    fit_gpytorch_model(mll)
    models.append(gp)
model = ModelListGP(*models)

acqf = qNoisyExpectedHypervolumeImprovement(model=model, ref_point=df_data[["alpha", "beta"]].min().values.tolist(), X_baseline=train_X)

candidate, acq_value = optimize_acqf(
    acqf, bounds=bounds, q=1, num_restarts=8, raw_samples=1024, equality_constraints=[(tensor([1, 2, 0]), tensor([1., 1., 1.]).to(**tkwargs), 1.0)])

print(candidate)
print(candidate[0,:3].sum())

** Stack trace/error message **

// Paste the bad output here!

Expected Behavior

The expected behavior would be that candidate[0,:3].sum() equals 1, but often it returns numbers larger than 1.

Additional context

This is the csv file with the data.
equality_constraint.csv

[Balandat]

Thanks for the repro! We'll look into what's going on here. I vaguely remember we had some issues in the past with now always properly generating feasible initial conditions in the presence of inequality and equality constraints, not sure if that got fixed (I'll start there).

[jduerholt]

Thanks. Could you reproduce the behavior?

[Balandat]

Yes, I was able to repro the issue.

[jduerholt]

Heyho, any news on this issue? Best, Johannes

[Balandat]

sorry, not yet - we'll have to take a closer look at this, let me see if anyone has bandwidth to investigate

[esantorella]

I think this is happening because optimization is failing within scipy.optimize.minimize.

If I add options={"disp": True} to the call to optimize_acqf in @jduerholt's example, it prints

Positive directional derivative for linesearch    (Exit mode 8)
            Current function value: -147055.3436559537
            Iterations: 6
            Function evaluations: 2
            Gradient evaluations: 2

This is telling us that optimization terminated because the linesearch failed. I reproduced this issue a bunch of times with different seeds, and whenever optimization terminated successfully, the solution did obey the constraints (up to machine precision).

That this is poor error handling in two ways:

1. Scipy's solutions can violate the constraints, and we don't do anything when this happens. We could handle this the way we do with nonlinear equality constraints, by returning the starting point if there is a violation.
2. Scipy's optimization can fail silently, returning a bad solution. It could be worth investigating whether that is a more widespread issue. We could raise an error whenever optimization is not a success, but I'm reluctant to do that in case there are a lot of decent runs being marked as failures.

What do you think we should do about this @Balandat @saitcakmak ?

As for why the optimization failed and how to make it succeed, I don't know at this point -- that would probably be a much more involved task.

[saitcakmak]

Thanks for investigating this @esantorella! In fit_gpytorch_model, we check for failures in minimize and re-try the model fitting with a new set of initial conditions. We do this by parsing the result object in fit_gpytorch_scipy and raising a warning if the optimization failed, then catching this same warning in fit_gpytorch_model. Do you think something like this would work here? At the very least, I think we should raise an error or a warning if the minimize fails rather than letting it go quietly.

[Balandat]

Do we know whether the ICs we pass to the optimizer are feasible to start with? We use our polytope samplers in case non-bound constraints are specified. In the past this sampler has had some bugs that caused samples to be infeasible (not sure if this is still an issue). Either way, it would make sense to double check that whatever we're passing as initial conditions to the optimizer is indeed feasible (o/w I wouldn't blame the optimizer to fail during line search).

Somewhat related: If we used the sample_around_best option to generate starting points close to the best observed points, then it looks like this will cause potentially (in the case of equality constraints likely) infeasible points to be passed as ICs if any non-bound constraints are specified (cc @sdaulton).

[esantorella]

@saitcakmak :

In fit_gpytorch_model, we check for failures in minimize

Oh cool, thanks. For my own future reference this is here and here.

Do you think something like this would work here? At the very least, I think we should raise an error or a warning if the minimize fails rather than letting it go quietly.

It makes sense to me to raise a warning within gen_candidates_scipy so anything downstream can decide how to handle it. Right now functions calling gen_candidates_scipy can't tell if it failed.

@Balandat :

Do we know whether the ICs we pass to the optimizer are feasible to start with? We use our polytope samplers in case non-bound constraints are specified.

They are feasible in this case, and we have a test for that. minimize starts with a feasible solution and then jumps away. I think this can happen if SLSQP evaluates a non-feasible point as part of a line search or numerical differentiation and then fails, but it is hard to be sure what this Fortran code is doing 😬

If we used the sample_around_best option to generate starting points close to the best observed points, then it looks like this will cause potentially (in the case of equality constraints likely) infeasible points to be passed as ICs if any non-bound constraints are specified (cc @sdaulton).

In this case we're not using sample_around_best. But someone could use sample_upon_best in combination with equality constraints, so maybe we should enforce that this is an invalid option.

[Balandat]

They are feasible in this case, and we have a test for that. minimize starts with a feasible solution and then jumps away. I think this can happen if SLSQP evaluates a non-feasible point as part of a line search or numerical differentiation and then fails, but it is hard to be sure what this Fortran code is doing 😬

Hmm interesting so the issue appears to lie with SLSQP then rather than on our end. I think catching the optimization failure, issuing a warning, and retrying with a new set of ICs is probably the right approach at least in the short term.

In this case we're not using sample_around_best. But someone could use sample_upon_best in combination with equality constraints, so maybe we should enforce that this is an invalid option.

Yes, let's do that. Would you mind putting up a PR for this that raises an informative UnsupportedError (or NotImplementedError if we want to be ambitious and support sampling form the feasible subspace in the future)?

[esantorella]

Sounds like a plan!

[jduerholt]

Hi all,

thanks for investigating this issue.

Best,

Johannes