Multiple constraints for single objective optimization

[jan-99]

Hello,

I was wondering how to implement multiple constraints for single objective optimization. Following the Closed-loop batch, constrained BO in BoTorch with qEI and qNEI tutorial:

def outcome_constraint(X):
    """L1 constraint; feasible if less than or equal to zero."""
    return X.sum(dim=-1) - 3

def weighted_obj(X):
    """Feasibility weighted objective; zero if not feasible."""
    return neg_hartmann6(X) * (outcome_constraint(X) <= 0).type_as(X)

Given that outcome_constraint returns the value it violates (or doesn't violate) the constraints and weighted_obj gives 0 if infeasible or the actual objective value if feasible, what is the proper way to extend this to use multiple constraints?

Assuming I add another synthetic constraint x1 + x2 ≤ 1, would it be appropiate to do the following?

def outcome_constraint(X):
    c1 = X.sum(dim=-1) - 3
    c2 = X[...,0] + X[...,1] - 1
    return torch.stack([c1, c2], dim=-1)

In this case, how should weighted_obj be defined? Should it return 0 as long as any constraint is violated and only return the actual objective value if all constraints are active?

Would greatly appreciate any help, thank you!

[Balandat]

Correct.

With your definition of outcome_constraint, you'd define weighted objective as something like

def weighted_obj(X):
    """Feasibility weighted objective; zero if not feasible."""
    return neg_hartmann6(X) * (outcome_constraint(X) <= 0).all(dim=-1).type_as(X)

[jan-99]

Hi @Balandat,

Thank you so much for your prompt reply!

Unfortunately I've met with difficulty in determining the correct shapes of the tensors for the optimization after using the above definitions of outcome_constraint and weighted_obj.

When executing this cell in the tutorial, I get the following error:

---------------------------------------------------------------------------
BotorchTensorDimensionError               Traceback (most recent call last)
~\AppData\Local\Temp/ipykernel_11532/2131489745.py in <module>
     29     # call helper functions to generate initial training data and initialize model
     30     train_x_ei, train_obj_ei, train_con_ei, best_observed_value_ei = generate_initial_data(n=10)
---> 31     mll_ei, model_ei = initialize_model(train_x_ei, train_obj_ei, train_con_ei)
     32 
     33     train_x_nei, train_obj_nei, train_con_nei = train_x_ei, train_obj_ei, train_con_ei

~\AppData\Local\Temp/ipykernel_11532/2027348400.py in initialize_model(train_x, train_obj, train_con, state_dict)
     20     # define models for objective and constraint
     21     model_obj = FixedNoiseGP(train_x, train_obj, train_yvar.expand_as(train_obj)).to(train_x)
---> 22     model_con = FixedNoiseGP(train_x, train_con, train_yvar.expand_as(train_con)).to(train_x)
     23     # combine into a multi-output GP model
     24     model = ModelListGP(model_obj, model_con)

~\Anaconda3\envs\concrete_opt\lib\site-packages\botorch\models\gp_regression.py in __init__(self, train_X, train_Y, train_Yvar, covar_module, outcome_transform, input_transform, **kwargs)
    215         if outcome_transform is not None:
    216             train_Y, train_Yvar = outcome_transform(train_Y, train_Yvar)
--> 217         self._validate_tensor_args(X=transformed_X, Y=train_Y, Yvar=train_Yvar)
    218         validate_input_scaling(
    219             train_X=transformed_X, train_Y=train_Y, train_Yvar=train_Yvar

~\Anaconda3\envs\concrete_opt\lib\site-packages\botorch\models\gpytorch.py in _validate_tensor_args(X, Y, Yvar, strict)
     81                         f" {Y.dim()}."
     82                     )
---> 83                 raise BotorchTensorDimensionError(message)
     84         else:
     85             warnings.warn(

BotorchTensorDimensionError: Expected X and Y to have the same number of dimensions (got X with dimension 2 and Y with dimension 3.

I've tried adjusting the shapes for generate_initial_data by adding dimensions to train_x and exact_obj:

def generate_initial_data(n=10):
    # generate training data
    train_x = torch.rand(10, 6, device=device, dtype=dtype)
    exact_obj = neg_hartmann6(train_x).unsqueeze(-1).unsqueeze(-1)  # add output dimensionx2
    exact_con = outcome_constraint(train_x).unsqueeze(-1)  # add output dimension
    train_obj = exact_obj + NOISE_SE * torch.randn_like(exact_obj)
    train_con = exact_con + NOISE_SE * torch.randn_like(exact_con)
    best_observed_value = weighted_obj(train_x).max().item()
    return train_x.unsqueeze(-1), train_obj, train_con, best_observed_value

but received this error:

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
~\AppData\Local\Temp/ipykernel_11532/2131489745.py in <module>
     45 
     46         # fit the models
---> 47         fit_gpytorch_model(mll_ei)
     48         fit_gpytorch_model(mll_nei)
     49 

~\Anaconda3\envs\concrete_opt\lib\site-packages\botorch\fit.py in fit_gpytorch_model(mll, optimizer, **kwargs)
     65         for mll_ in mll.mlls:
     66             fit_gpytorch_model(
---> 67                 mll=mll_, optimizer=optimizer, max_retries=max_retries, **kwargs
     68             )
     69         return mll

~\Anaconda3\envs\concrete_opt\lib\site-packages\botorch\fit.py in fit_gpytorch_model(mll, optimizer, **kwargs)
    123                 sample_all_priors(mll.model)
    124             try:
--> 125                 mll, _ = optimizer(mll, track_iterations=False, **kwargs)
    126             except NotPSDError:
    127                 retry += 1

~\Anaconda3\envs\concrete_opt\lib\site-packages\botorch\optim\fit.py in fit_gpytorch_scipy(mll, bounds, method, options, track_iterations, approx_mll, scipy_objective, module_to_array_func, module_from_array_func)
    245             jac=True,
    246             options=options,
--> 247             callback=cb,
    248         )
    249         iterations = []

~\AppData\Roaming\Python\Python37\site-packages\scipy\optimize\_minimize.py in minimize(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options)
    608     elif meth == 'l-bfgs-b':
    609         return _minimize_lbfgsb(fun, x0, args, jac, bounds,
--> 610                                 callback=callback, **options)
    611     elif meth == 'tnc':
    612         return _minimize_tnc(fun, x0, args, jac, bounds, callback=callback,

~\AppData\Roaming\Python\Python37\site-packages\scipy\optimize\lbfgsb.py in _minimize_lbfgsb(fun, x0, args, jac, bounds, disp, maxcor, ftol, gtol, eps, maxfun, maxiter, iprint, callback, maxls, **unknown_options)
    343             # until the completion of the current minimization iteration.
    344             # Overwrite f and g:
--> 345             f, g = func_and_grad(x)
    346         elif task_str.startswith(b'NEW_X'):
    347             # new iteration

~\AppData\Roaming\Python\Python37\site-packages\scipy\optimize\lbfgsb.py in func_and_grad(x)
    293     else:
    294         def func_and_grad(x):
--> 295             f = fun(x, *args)
    296             g = jac(x, *args)
    297             return f, g

~\AppData\Roaming\Python\Python37\site-packages\scipy\optimize\optimize.py in function_wrapper(*wrapper_args)
    325     def function_wrapper(*wrapper_args):
    326         ncalls[0] += 1
--> 327         return function(*(wrapper_args + args))
    328 
    329     return ncalls, function_wrapper

~\AppData\Roaming\Python\Python37\site-packages\scipy\optimize\optimize.py in __call__(self, x, *args)
     63     def __call__(self, x, *args):
     64         self.x = numpy.asarray(x).copy()
---> 65         fg = self.fun(x, *args)
     66         self.jac = fg[1]
     67         return fg[0]

~\Anaconda3\envs\concrete_opt\lib\site-packages\botorch\optim\utils.py in _scipy_objective_and_grad(x, mll, property_dict)
    215         loss = -mll(*args).sum()
    216     except RuntimeError as e:
--> 217         return _handle_numerical_errors(error=e, x=x)
    218     loss.backward()
    219     param_dict = OrderedDict(mll.named_parameters())

~\Anaconda3\envs\concrete_opt\lib\site-packages\botorch\optim\utils.py in _handle_numerical_errors(error, x)
    242     ):
    243         return float("nan"), np.full_like(x, "nan")
--> 244     raise error  # pragma: nocover
    245 
    246 

~\Anaconda3\envs\concrete_opt\lib\site-packages\botorch\optim\utils.py in _scipy_objective_and_grad(x, mll, property_dict)
    213         output = mll.model(*train_inputs)
    214         args = [output, train_targets] + _get_extra_mll_args(mll)
--> 215         loss = -mll(*args).sum()
    216     except RuntimeError as e:
    217         return _handle_numerical_errors(error=e, x=x)

~\Anaconda3\envs\concrete_opt\lib\site-packages\gpytorch\module.py in __call__(self, *inputs, **kwargs)
     28 
     29     def __call__(self, *inputs, **kwargs):
---> 30         outputs = self.forward(*inputs, **kwargs)
     31         if isinstance(outputs, list):
     32             return [_validate_module_outputs(output) for output in outputs]

~\Anaconda3\envs\concrete_opt\lib\site-packages\gpytorch\mlls\exact_marginal_log_likelihood.py in forward(self, function_dist, target, *params)
     60         # Get the log prob of the marginal distribution
     61         output = self.likelihood(function_dist, *params)
---> 62         res = output.log_prob(target)
     63         res = self._add_other_terms(res, params)
     64 

~\Anaconda3\envs\concrete_opt\lib\site-packages\gpytorch\distributions\multivariate_normal.py in log_prob(self, value)
    145     def log_prob(self, value):
    146         if settings.fast_computations.log_prob.off():
--> 147             return super().log_prob(value)
    148 
    149         if self._validate_args:

~\Anaconda3\envs\concrete_opt\lib\site-packages\torch\distributions\multivariate_normal.py in log_prob(self, value)
    207         if self._validate_args:
    208             self._validate_sample(value)
--> 209         diff = value - self.loc
    210         M = _batch_mahalanobis(self._unbroadcasted_scale_tril, diff)
    211         half_log_det = self._unbroadcasted_scale_tril.diagonal(dim1=-2, dim2=-1).log().sum(-1)

RuntimeError: The size of tensor a (2) must match the size of tensor b (6) at non-singleton dimension 1

I've also tried adjusting the shape of the constraints instead by omitting unsqueeze(-1) for exact_con:

def generate_initial_data(n=10):
    # generate training data
    train_x = torch.rand(10, 6, device=device, dtype=dtype)
    exact_obj = neg_hartmann6(train_x).unsqueeze(-1)  # add output dimension
    exact_con = outcome_constraint(train_x)#.unsqueeze(-1)  # DO NOT add output dimension
    train_obj = exact_obj + NOISE_SE * torch.randn_like(exact_obj)
    train_con = exact_con + NOISE_SE * torch.randn_like(exact_con)
    best_observed_value = weighted_obj(train_x).max().item()
    return train_x, train_obj, train_con, best_observed_value

but received this error:

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
~\AppData\Local\Temp/ipykernel_11532/2131489745.py in <module>
     67 
     68         # optimize and get new observation
---> 69         new_x_ei, new_obj_ei, new_con_ei = optimize_acqf_and_get_observation(qEI)
     70         new_x_nei, new_obj_nei, new_con_nei = optimize_acqf_and_get_observation(qNEI)
     71 

~\AppData\Local\Temp/ipykernel_11532/2918212314.py in optimize_acqf_and_get_observation(acq_func)
     18         num_restarts=NUM_RESTARTS,
     19         raw_samples=RAW_SAMPLES,  # used for intialization heuristic
---> 20         options={"batch_limit": 5, "maxiter": 200},
     21     )
     22     # observe new values

~\Anaconda3\envs\concrete_opt\lib\site-packages\botorch\optim\optimize.py in optimize_acqf(acq_function, bounds, q, num_restarts, raw_samples, options, inequality_constraints, equality_constraints, fixed_features, post_processing_func, batch_initial_conditions, return_best_only, sequential, **kwargs)
    191             options=options,
    192             inequality_constraints=inequality_constraints,
--> 193             equality_constraints=equality_constraints,
    194         )
    195 

~\Anaconda3\envs\concrete_opt\lib\site-packages\botorch\optim\initializers.py in gen_batch_initial_conditions(acq_function, bounds, q, num_restarts, raw_samples, fixed_features, options, inequality_constraints, equality_constraints)
    182                     end_idx = min(start_idx + batch_limit, X_rnd.shape[0])
    183                     Y_rnd_curr = acq_function(
--> 184                         X_rnd[start_idx:end_idx].to(device=device)
    185                     ).cpu()
    186                     Y_rnd_list.append(Y_rnd_curr)

~\Anaconda3\envs\concrete_opt\lib\site-packages\torch\nn\modules\module.py in _call_impl(self, *input, **kwargs)
   1100         if not (self._backward_hooks or self._forward_hooks or self._forward_pre_hooks or _global_backward_hooks
   1101                 or _global_forward_hooks or _global_forward_pre_hooks):
-> 1102             return forward_call(*input, **kwargs)
   1103         # Do not call functions when jit is used
   1104         full_backward_hooks, non_full_backward_hooks = [], []

~\Anaconda3\envs\concrete_opt\lib\site-packages\botorch\utils\transforms.py in decorated(cls, X, **kwargs)
    268         if cls.X_pending is not None:
    269             X = torch.cat([X, match_batch_shape(cls.X_pending, X)], dim=-2)
--> 270         return method(cls, X, **kwargs)
    271 
    272     return decorated

~\Anaconda3\envs\concrete_opt\lib\site-packages\botorch\utils\transforms.py in decorated(acqf, X, *args, **kwargs)
    227             # add t-batch dim
    228             X = X if X.dim() > 2 else X.unsqueeze(0)
--> 229             output = method(acqf, X, *args, **kwargs)
    230             if assert_output_shape and not _verify_output_shape(
    231                 acqf=acqf,

~\Anaconda3\envs\concrete_opt\lib\site-packages\botorch\acquisition\monte_carlo.py in forward(self, X)
    154             model and input `X`.
    155         """
--> 156         posterior = self.model.posterior(X)
    157         samples = self.sampler(posterior)
    158         obj = self.objective(samples, X=X)

~\Anaconda3\envs\concrete_opt\lib\site-packages\botorch\models\gpytorch.py in posterior(self, X, output_indices, observation_noise, **kwargs)
    564                 mvn_gen = zip(output_indices, mvns)
    565             else:
--> 566                 mvns = self(*transformed_X)
    567                 if observation_noise is not False:
    568                     mvnX = [(mvn, transformed_X[i]) for i, mvn in enumerate(mvns)]

~\Anaconda3\envs\concrete_opt\lib\site-packages\gpytorch\models\model_list.py in __call__(self, *args, **kwargs)
     81 
     82     def __call__(self, *args, **kwargs):
---> 83         return [model.__call__(*args_, **kwargs) for model, args_ in zip(self.models, _get_tensor_args(*args))]
     84 
     85     @property

~\Anaconda3\envs\concrete_opt\lib\site-packages\gpytorch\models\model_list.py in <listcomp>(.0)
     81 
     82     def __call__(self, *args, **kwargs):
---> 83         return [model.__call__(*args_, **kwargs) for model, args_ in zip(self.models, _get_tensor_args(*args))]
     84 
     85     @property

~\Anaconda3\envs\concrete_opt\lib\site-packages\gpytorch\models\exact_gp.py in __call__(self, *args, **kwargs)
    297                     train_input = train_input.expand(*batch_shape, *train_input.shape[-2:])
    298                 if batch_shape != input.shape[:-2]:
--> 299                     batch_shape = _mul_broadcast_shape(batch_shape, input.shape[:-2])
    300                     train_input = train_input.expand(*batch_shape, *train_input.shape[-2:])
    301                     input = input.expand(*batch_shape, *input.shape[-2:])

~\Anaconda3\envs\concrete_opt\lib\site-packages\gpytorch\utils\broadcasting.py in _mul_broadcast_shape(error_msg, *shapes)
     18             if any(size != non_singleton_sizes[0] for size in non_singleton_sizes):
     19                 if error_msg is None:
---> 20                     raise RuntimeError("Shapes are not broadcastable for mul operation")
     21                 else:
     22                     raise RuntimeError(error_msg)

RuntimeError: Shapes are not broadcastable for mul operation

Could you please assist in showing how the shapes of the tensors should be adjusted? Thank you very much!

[Balandat]

Looks like you're somehow getting an additional dimension, but I'm not sure where it's coming from without a full e2e runnable example.

What is the shape of X that you're using to evaluate the constraints? Is this a n x d tensor (usually what the external interface uses) or a (batch_shape) x n x d tensor (which is used internally in most of the computations in botorch)?

[jan-99]

This is the notebook I'm trying to run: closed_loop_botorch_only.txt

It is the same as the tutorial Closed-loop batch, constrained BO in BoTorch with qEI and qNEI, but with outcome_constraint and weighted_obj defined as such:

def outcome_constraint(X):
    c1 = X.sum(dim=-1) - 3
    c2 = X[...,0] + X[...,1] - 1
    return torch.stack([c1, c2], dim=-1)

def weighted_obj(X):
    """Feasibility weighted objective; zero if not feasible."""
    return neg_hartmann6(X) * (outcome_constraint(X) <= 0).all(dim=-1).type_as(X)

X is a n x d tensor.

[Balandat]

Sorry I didn't give you the right advice here the first time around.

Basically if you have multiple outcome constraints on black box functions of which you can observe the value (in your case those are X.sum(dim=-1) - 3 and X[...,0] + X[...,1] - 1), then what you'd do is model each of them individually in separate models.

So instead of one model_con GP, what you'd do is have two of them, one for each outcome constraint (using the same code as for one outcome constraint in the tutorial). You would then combine the three models into a ModelListGP via ModelListGP(model_obj, model_con_1, model_con_2)

For the weighted objective you then can do sth like

def weighted_obj(X):
    """Feasibility weighted objective; zero if not feasible."""
    feasible = (outcome_constraint_1(X) <= 0) & (outcome_constraint_2(X) <= 0)
    return neg_hartmann6(X) * feasible.type_as(X)

And for the ConstrainedMCObjective you'd have to pick out the two separate constraints as follows"

def obj_callable(Z):
    return Z[..., 0]

def constraint_callable_1(Z):
    return Z[..., 1]

def constraint_callable_2(Z):
    return Z[..., 2]


# define a feasibility-weighted objective for optimization
constrained_obj = ConstrainedMCObjective(
    objective=obj_callable,
    constraints=[constraint_callable_1, constraint_callable_2],
)

Hope that helps.

[jan-99]

It works now, thanks so much for your assistance!