More about ic_generator and batch_initial_conditions

[madhavkrishnan94]

This may be considered as a continuation of the discussion from #1783 .

I have a few questions about the use of ic_generator and batch_initial_conditions that have to be used in addition to the nonlinear inequality constraint.

I took a look at the code in Scipy non-linear inequality constraints.ipynb.txt. Here, the get_batch_initial_conditions function invokes another function that generates a Sobol sequence.

def get_feasible_sobol_points(n):
    """Sobol sequence where we randomly set three of the parameters to zero to satisfy the constraint"""
    X = SobolEngine(dimension=6, scramble=True).draw(n).to(torch.double)
    inds = torch.argsort(torch.rand(n, 6), dim=-1)[:, :3]
    X[torch.arange(X.shape[0]).unsqueeze(-1), inds] = 0
    return X


def get_batch_initial_conditions(num_restarts, raw_samples, acqf):
    X = get_feasible_sobol_points(raw_samples).unsqueeze(1)
    acq_vals = acqf(X)
    return X[acq_vals.topk(num_restarts).indices]

I see that the dataset is initialized with 10 points generated using the get_feasible_sobol_points function (y-values calculated using the hartmann function). This data is then fed to a for loop where, for every iteration, the model is retrained and the acquisition function is optimized with the nonlinear inequality constraint and the initial conditions passed as arguments in the optimize_acqf function. And the generation of new points from the BO happens in an iterative process.

for i in range(40):
    gp = SingleTaskGP(X, Y, outcome_transform=Standardize(m=1))
    mll = ExactMarginalLogLikelihood(gp.likelihood, gp)
    fit_gpytorch_model(mll)

    EI = ExpectedImprovement(model=gp, best_f=Y.min(), maximize=False)

    batch_initial_conditions = get_batch_initial_conditions(num_restarts=20, raw_samples=512, acqf=EI)

    candidate, acq_value = optimize_acqf(
        EI,
        bounds=torch.cat((torch.zeros(1, 6), torch.ones(1, 6))),
        q=1,
        nonlinear_inequality_constraints=[ineq_constraint],
        batch_initial_conditions=batch_initial_conditions,
        num_restarts=20,
    )

Firstly, is this iterative method the only way to generate new points? Is generating Sobol points necessary always? Is there an equivalent method to impose the nonlinear inequality constraints in one shot as there is for linear inequality constraints?

Secondly, I would like to know how the ic_generator works. Can you let me know how they are different from the batch_initial_conditions. I see that the docs don't mention much beyond the fact that it

Defaults to gen_one_shot_kg_initial_conditions for qKnowledgeGradient acquisition functions and gen_batch_initial_conditions otherwise. Must be specified for nonlinear inequality constraints.

Just providing additional context, in case it helps: the data I work with is fully from experimental measurements. How do I generate initial conditions in my case?

Thanks in advance!

[Balandat]

Firstly, is this iterative method the only way to generate new points? Is generating Sobol points necessary always? Is there an equivalent method to impose the nonlinear inequality constraints in one shot as there is for linear inequality constraints?

This depends on the nature of the constraints. Maybe in some cases, but maybe not in others. This is going to be specific to the problem and not really something that we can handle in all generality on the botorch end. We just need the ic_generator to provide these feasible initial conditions.

Secondly, I would like to know how the ic_generator works. Can you let me know how they are different from the batch_initial_conditions. I see that the docs don't mention much beyond the fact that it

ic_generator is just a callable with a particular signature that generates the initial conditions. You can either pass the conditions explicitly (e.g. as a batch_initial_conditions) or provide the ic_generator callable. The latter will be useful if you want to generate batches of points in a sequential greedy fashion (the default way that optimize_acqf works) - in that case ic_generator will be called under the hood in a loop to generate the initial conditions for each new candidate point.

Just providing additional context, in case it helps: the data I work with is fully from experimental measurements. How do I generate initial conditions in my case?

What are the constraints you have on your parameters? That will dictate how to generate the initial conditions for optimization.

[madhavkrishnan94]

The only constraints I impose are the fact that I want the newly generated candidates to always lie in the feasible region.

I'm not sure what you mean by parameters but these are the only conditions.

Can you provide an example of an ic_generator implementation?

Thanks!

[madhavkrishnan94]

Hello @Balandat,

Please consider this query along the same lines as those discussed in issues #1783 and #1572 .

I'm trying to define a nonlinear inequality constraint and use it to generate new candidates inside the optimize_acqf function.

This is how I define the callable. x is a 3-dimensional tensor but I am interested in only the first two dimensions to evaluate the feasibility.

def nl_con(x):
        if x[0] == 64 and x[1] == 131:
            return 0 #binding
        elif x[0] >= 64 and x[1] > 131:
            return -1 #negative : infeasible
        else:
            return 1 #positive : feasible

Now, I follow the discussion in #1572 to define the ic_generator using

from functools import partial

_gen_batch_initial_conditions_nonlinear = partial(
    gen_batch_initial_conditions_nonlinear,
    nonlinear_constraint=nl_con,)

and

from botorch.optim.initializers import gen_batch_initial_conditions

def get_batch_initial_conditions_nonlinear(
    acq_function,
    bounds,
    q,
    num_restarts,
    raw_samples,
    **kwargs
):
    # your code for generating a `num_restarts x q x d` tensor of
    # feasible initial conditions (s.t. bounds and nonlinear constraints) goes here.
    need = num_restarts * q
    all_points = torch.tensor([])
    ndims = bounds.shape[1]
    while all_points.shape[0] <= need:
        # Using this as the initial guess seemed reasonable...
        res = gen_batch_initial_conditions(
            acq_function=acq_function,
            bounds=bounds,
            q=q,
            num_restarts=num_restarts,
            raw_samples=raw_samples,
            **kwargs
        )
        res = res.reshape(-1, ndims)
        where = torch.where(nl_con(X_rnd) >= 0)[0]
        all_points = torch.cat([all_points, res[where, :]], axis=0)
    
    return all_points[:need, :].reshape(num_restarts, q, ndims)

I run the optimize_acqf function in the following manner.

candidates, _ = optimize_acqf(
        acq_function=acq_func,
        bounds=constrained_bounds,
        q=BATCH_SIZE,
        num_restarts=NUM_RESTARTS,
        raw_samples=RAW_SAMPLES,  # used for intialization heuristic
        #options={"batch_limit": 5, "maxiter": 200}, #additional hyperparameters
        nonlinear_inequality_constraints=[nl_con],
        ic_generator=_gen_batch_initial_conditions_nonlinear,
        sequential=True)

However, this returns the following error: TypeError: gen_batch_initial_conditions() got an unexpected keyword argument 'nonlinear_constraint'.

Another doubt I have is regarding the use of X_rnd. How do we define this in the first place?

Any leads in helping me resolve these issues would be highly appreciated. Any other suggestions/error corrections are welcome too.

Thanks!

[Balandat]

You shouldn't make the constraints discontinuous, this will make it hard for the optimizer. See my comment on the issue you referenced: #1783 (comment)

As to the error - my suspicion is that in this block in your code you're passing nonlinear_constraint via **kwargs also to gen_batch_initial_conditions, which would throw this error:

        res = gen_batch_initial_conditions(
            acq_function=acq_function,
            bounds=bounds,
            q=q,
            num_restarts=num_restarts,
            raw_samples=raw_samples,
            **kwargs
        )

Generally, it's hard to tell what is going on without a full stack trace - can you please share those in future issues?

[madhavkrishnan94]

Thanks a lot, @Balandat . Allow me to divide my response into three parts.

---

1. I took at the comment referenced. The feasibility criterion is arrived upon empirically and I don't have an explicit function to compute the 0,+/-1 that I want to return. In this case, I wonder how I can ensure that the callable's returning values are not discontinuous.

2. Here's the stack trace as requested.

TypeError                                 Traceback (most recent call last)
/tmp/ipykernel_1426258/105757521.py in <module>
     14 
     15 #for i in range(40):
---> 16 candidates, _ = optimize_acqf(
     17             acq_function=acq_func,
     18             bounds=constrained_bounds,

~/.local/lib/python3.10/site-packages/botorch/optim/optimize.py 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)
    542         ic_gen_kwargs=ic_gen_kwargs,
    543     )
--> 544     return _optimize_acqf(opt_acqf_inputs)
    545 
    546 

~/.local/lib/python3.10/site-packages/botorch/optim/optimize.py in _optimize_acqf(opt_inputs)
    571 
    572     # Batch optimization (including the case q=1)
--> 573     return _optimize_acqf_batch(
    574         opt_inputs=opt_inputs, start_time=start_time, timeout_sec=timeout_sec
    575     )

~/.local/lib/python3.10/site-packages/botorch/optim/optimize.py in _optimize_acqf_batch(opt_inputs, start_time, timeout_sec)
    267     else:
    268         # pyre-ignore[28]: Unexpected keyword argument `acq_function` to anonymous call.
--> 269         batch_initial_conditions = opt_inputs.get_ic_generator()(
    270             acq_function=opt_inputs.acq_function,
    271             bounds=opt_inputs.bounds,

/tmp/ipykernel_1426258/3612278866.py in get_batch_initial_conditions_nonlinear(acq_function, bounds, q, num_restarts, raw_samples, **kwargs)
     16     while all_points.shape[0] <= need:
     17         # Using this as the initial guess seemed reasonable...
---> 18         res = gen_batch_initial_conditions(
     19             acq_function=acq_function,
     20             bounds=bounds,

TypeError: gen_batch_initial_conditions() got an unexpected keyword argument 'nonlinear_constraint'

3. Regarding the **kwargs : I removed the **kwargs from the function argument and ran it.

    while all_points.shape[0] <= need:
        # Using this as the initial guess seemed reasonable...
        res = gen_batch_initial_conditions(
            acq_function=acq_function,
            bounds=bounds,
            q=q,
            num_restarts=num_restarts,
            raw_samples=raw_samples,
            #**kwargs
        )
        res = res.reshape(-1, ndims)
        where = torch.where(nl_con(X_rnd) >= 0)[0]
        all_points = torch.cat([all_points, res[where, :]], axis=0)

It then throws up an error about X_rnd not being defined. Here's the corresponding stack trace.

NameError Traceback (most recent call last)
/tmp/ipykernel_1426258/105757521.py in
14
15 #for i in range(40):
---> 16 candidates, _ = optimize_acqf(
17 acq_function=acq_func,
18 bounds=constrained_bounds,

~/.local/lib/python3.10/site-packages/botorch/optim/optimize.py 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)
542 ic_gen_kwargs=ic_gen_kwargs,
543 )
--> 544 return _optimize_acqf(opt_acqf_inputs)
545
546

~/.local/lib/python3.10/site-packages/botorch/optim/optimize.py in _optimize_acqf(opt_inputs)
571
572 # Batch optimization (including the case q=1)
--> 573 return _optimize_acqf_batch(
574 opt_inputs=opt_inputs, start_time=start_time, timeout_sec=timeout_sec
575 )

~/.local/lib/python3.10/site-packages/botorch/optim/optimize.py in _optimize_acqf_batch(opt_inputs, start_time, timeout_sec)
267 else:
268 # pyre-ignore[28]: Unexpected keyword argument acq_function to anonymous call.
--> 269 batch_initial_conditions = opt_inputs.get_ic_generator()(
270 acq_function=opt_inputs.acq_function,
271 bounds=opt_inputs.bounds,

/tmp/ipykernel_1426258/3372450131.py in get_batch_initial_conditions_nonlinear(acq_function, bounds, q, num_restarts, raw_samples, **kwargs)
25 )
26 res = res.reshape(-1, ndims)
---> 27 where = torch.where(nl_con(X_rnd) >= 0)[0]
28 all_points = torch.cat([all_points, res[where, :]], axis=0)
29

NameError: name 'X_rnd' is not defined

This takes me back to the query in my previous question.

Another doubt I have is regarding the use of X_rnd. How do we define this in the first place?

I gather that X_rnd is the random starting point(s) that you generate.

---

Thanks in advance!

[Balandat]

In this case, I wonder how I can ensure that the callable's returning values are not discontinuous.

How about something like this?

def nl_con(x):
    d0 = x[0] - 64
    d1 = x[1] - 131
    if d0 > 0 and d1 > 0:
        return - (d0 * d1)
    elif d0 > 0 and d1 <= 0:
        return d1 - d0
    elif d0 <= 0 and d1 > 0:
        return d0 - d1
    else:
        return d0 * d1

This is positive in the feasible region, zero on the boundary, and negative in the infeasible region. It's also differentiable almost everywhere.

Another doubt I have is regarding the use of X_rnd. How do we define this in the first place?

I have no idea what X_rnd is in your code. Did you mean to set X_rnd = gen_batch_initial_conditions? Looks like those should be the unconstrained random ICs that you then filter to be the ones satisfying your constraints.

[madhavkrishnan94]

This is positive in the feasible region, zero on the boundary, and negative in the infeasible region. It's also differentiable almost everywhere.

Thanks, this seems to be a good fix.

Also, X_rnd is from your discussion in this thread. I guess, X_rnd generates candidates which are then subjected to feasibility test (under the nonlinear inequality constraint).

where = torch.where(nonlinear_constraint(X_rnd) >= 0)[0]