Outcome constraints using constraint parameter in CEI

[GStechschulte]

In a lot of applications there is an outcome constraint (lower and or upper bound) on the black box function that you are wanting to optimize lower < f(x) < upper. For example, with industrial machine processes this may be minimum and maximum belt speeds, heating elements, watts, etc. I am trying to implement a constrained expected improvement policy with an upper bound of 0.0 on a simulated example to demonstrate the use of outcome constraints.

Using the forrester function as the objective function, I am attempting to find a set of points that minimize this function.

import torch

from botorch.acquisition import ConstrainedExpectedImprovement
from botorch.fit import fit_gpytorch_model
from botorch.models import SingleTaskGP
from botorch.optim import optimize_acqf
from gpytorch.mlls import ExactMarginalLogLikelihood

def objective_fn(x):
    return -((x + 1) ** 2) * torch.sin(2 * x + 2) / 5 + 1 + x / 3

lb, ub = -5, 5
bounds = torch.tensor([[lb], [ub]], dtype=torch.float)
xs = torch.linspace(lb, ub, 100).unsqueeze(1)
ys = -objective_fn(xs) # minimize --> flip the sign

n = 10
train_x = bounds[0] + (bounds[1] - bounds[0]) * torch.rand(n, 1)
train_y = -objective_fn(train_x) # minimize --> flip the sign

plt.figure(figsize=(4, 2))
plt.plot(xs, ys, label="Objective")
plt.scatter(train_x, train_y, color="black", label="Observations")
plt.axhline(y=0, color="k", linestyle="--", label="Upper Bound")
plt.legend();

Since BoTorch assumes maximization I flip the objective and set an upper bound outcome constraint of 0.0. Thus, I use ConstrainedExpectedImprovement with the constraint constraints = {0: (None, 0.0)} to model the non-positivity constraint.

I am getting this idea from the docstring of ConstrainedExpectedImprovement .

N_ITERATIONS = 20
for iteration in range(N_ITERATIONS):
    model = SingleTaskGP(train_x, train_y)
    mll = ExactMarginalLogLikelihood(model.likelihood, model)
    fit_gpytorch_model(mll)

    constraints = {0: (None, 0.0)}  # 0th output should be non-negative

    cEI = ConstrainedExpectedImprovement(
        model, 
        best_f=train_y[train_y < 0].max(), 
        objective_index=1, 
        constraints=constraints
    )

    new_x, _ = optimize_acqf(
        acq_function=cEI,
        bounds=bounds,
        q=1,  
        num_restarts=5,
        raw_samples=20,
    )

    new_y = -objective_fn(new_x)

    train_x = torch.cat([train_x, new_x])
    train_y = torch.cat([train_y, new_y])

However, when running this example, I get the following error:

---------------------------------------------------------------------------
{
	"name": "RuntimeError",
	"message": "zero-dimensional tensor (at position 0) cannot be concatenated",
	"stack": "---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)

File ~/miniforge3/envs/bayesopt/lib/python3.11/site-packages/botorch/optim/optimize.py:547, in optimize_acqf(acq_function, bounds, q, num_restarts, raw_samples, options, inequality_constraints, equality_constraints, nonlinear_inequality_constraints, fixed_features, post_processing_func, batch_initial_conditions, return_best_only, gen_candidates, sequential, ic_generator, timeout_sec, return_full_tree, retry_on_optimization_warning, **ic_gen_kwargs)
    524     gen_candidates = gen_candidates_scipy
    525 opt_acqf_inputs = OptimizeAcqfInputs(
    526     acq_function=acq_function,
    527     bounds=bounds,
   (...)
    545     ic_gen_kwargs=ic_gen_kwargs,
    546 )
--> 547 return _optimize_acqf(opt_acqf_inputs)

File ~/miniforge3/envs/bayesopt/lib/python3.11/site-packages/botorch/optim/optimize.py:569, in _optimize_acqf(opt_inputs)
    566     return _optimize_acqf_sequential_q(opt_inputs=opt_inputs)
    568 # Batch optimization (including the case q=1)
--> 569 return _optimize_acqf_batch(opt_inputs=opt_inputs)

File ~/miniforge3/envs/bayesopt/lib/python3.11/site-packages/botorch/optim/optimize.py:266, in _optimize_acqf_batch(opt_inputs)
    263     batch_initial_conditions = opt_inputs.batch_initial_conditions
    264 else:
    265     # pyre-ignore[28]: Unexpected keyword argument `acq_function` to anonymous call.
--> 266     batch_initial_conditions = opt_inputs.get_ic_generator()(
    267         acq_function=opt_inputs.acq_function,
    268         bounds=opt_inputs.bounds,
    269         q=opt_inputs.q,
    270         num_restarts=opt_inputs.num_restarts,
    271         raw_samples=opt_inputs.raw_samples,
    272         fixed_features=opt_inputs.fixed_features,
    273         options=options,
    274         inequality_constraints=opt_inputs.inequality_constraints,
    275         equality_constraints=opt_inputs.equality_constraints,
    276         **opt_inputs.ic_gen_kwargs,
    277     )
    279 batch_limit: int = options.get(
    280     \"batch_limit\",
    281     opt_inputs.num_restarts
    282     if not opt_inputs.nonlinear_inequality_constraints
    283     else 1,
    284 )
    286 def _optimize_batch_candidates() -> Tuple[Tensor, Tensor, List[Warning]]:

File ~/miniforge3/envs/bayesopt/lib/python3.11/site-packages/botorch/optim/initializers.py:393, in gen_batch_initial_conditions(acq_function, bounds, q, num_restarts, raw_samples, fixed_features, options, inequality_constraints, equality_constraints, generator)
    391         Y_rnd_list.append(Y_rnd_curr)
    392         start_idx += batch_limit
--> 393     Y_rnd = torch.cat(Y_rnd_list)
    394 batch_initial_conditions = init_func(
    395     X=X_rnd, Y=Y_rnd, n=num_restarts, **init_kwargs
    396 ).to(device=device)
    397 if not any(issubclass(w.category, BadInitialCandidatesWarning) for w in ws):

RuntimeError: zero-dimensional tensor (at position 0) cannot be concatenated"
}

I think the issue is because this class is expecting the model to be multi-outcome. Indeed, the Closed-loop batch, constrained BO in BoTorch uses two GPs: one for the objective function and one for the constraint function. However, the docstring of CEI makes it seem as though you can use a SingleTaskGP and simply define an outcome constraint with a lower and upper bound.

Most of the tutorials I see for outcome constraints assume we know the cost function and thus train a GP for both the objective and cost function. However, I am interested in cases where all I know about the outcome constraint is the lower and upper bound of the domain according to technical specifications. Is it possible to achieve this without having to train two GPs: one for the objective function and one for the cost function?

Thanks for your time and help!

[Balandat]

So to distill your setup into a single sentence: You want to optimize some black box function f but limit that optimization to some range lower < f(x) < upper(x). Is this correct?

If so, then this is something that we haven't really considered (typically you want to constrain on something other than the objective), which is likely why you're seeing this low-level failure. If you confirm that this is indeed what you want to do we can look into how to make that work.

[GStechschulte]

Thanks for the information 😄 However, when I use the constraint you provided above, I still get proposed data points $x$ that result in objective values that violate the constraint $f(x) < 0$.

Reading the docstring of SampleReducingMCAcquisitionFunction it appears that the callable (a lambda function in this case) should evaluate the posterior samples and return a negative value if the sample satisfies the constraint and positive otherwise. So I should be able to do

    EI = botorch.acquisition.qLogExpectedImprovement(
        model, 
        best_f=train_y[train_y < 0].max(), 
        constraints=[lambda Z: -1 * (Z.squeeze(-1) < 0)]
    )

where the lambda function returns -1 if an element in the posterior sample is < 0 and 0 otherwise.

Even using this lambda function results in proposed objective values that violate the constraint.

Thanks for your patience and help!

[sdaulton]

You'd want to keep the constraint differentiable w.r.t to the input Z. So the constraints here would just be constraints=[lambda Z: Z.squeeze(-1)]. Note that Z here would be passing in all outcomes, but typically one of your outcomes would be the objective and another outcome would be the constraint, so you would want to index into Z appropriately and make separate callables for each constraint outcome. e.g. constraints=[lambda Z: Z[..., constraint_outcome_idx]].

[GStechschulte]

Thanks @sdaulton. So constraint is used to index the outcome with its corresponding constraint. Since my problem only has one outcome, .squeeze(-1) selects the correct (and only) outcome dimension.

But now, I am back to to my problem of how to impose an outcome constraint of the form lower < f(x) < upper where f is some black-box function and lower and upper are bounds that I want to enforce such that the optimization routine only proposes values that result in lower < f(x) < upper.

On one hand, I think this should be solvable using an outcome constraint. On another, it is partly happening because the GP is "allowed" to make such predictions. Perhaps using a Truncated Normal likelihood with lower and or upper bounds is also feasible.

[Balandat]

But now, I am back to to my problem of how to impose an outcome constraint of the form lower < f(x) < upper where f is some black-box function and lower and upper are bounds that I want to enforce such that the optimization routine only proposes values that result in lower < f(x) < upper.

constraints=[
    lambda Z: Z.squeeze(-1) - upper,
    lambda Z: lower - Z.squeeze(-1),
]

should do the trick, no?

[GStechschulte]

Yes, this worked! Thanks a lot @Balandat and @sdaulton for the help and the work on the great package 👍🏼