MultiObjectBO using product of some PIs

[taikichi-pc]

Hi. Thank you for your development. I love botorch.
Today I would like to ask you how should I write a code about one method of MOBO.

I'd like to use product of each object PI. (sometimes sum of logPI)
Because I often use EHVI based acquisition function mainly.
But calculating hypervolume is sometimes heavy for me.

For example, I have two objects to optimize. And I want to maximize both of them.
I'd like to get next candidate x based on product of probability of improvement.

When I fulfill above condition, how should I code with your Botorch.

Thank you in advance for your consideration.

[Balandat]

Basically what I would do is subclass [MultiObjectiveAnalyticAcquisitionFunction](https://github.com/pytorch/botorch/blob/main/botorch/acquisition/multi_objective/analytic.py#L42 and implement a forward method that computes PI as done in the single-output case here and then multiplies the probability values together.

Another option would be to use the posterior_transform argument to ProbabilityOfImprovement to extract the individual outputs of he multi-output model and just compute two ProbabilityOfImprovement values that you multiply. But that has the downside that you'd call posterior on the multi-output model twice which will cause unnecessary compute.

[taikichi-pc]

Thank you a lot for quick reply.
It make sense a little to me.
I am understanding above [MultiObjectiveAnalyticAcquisitionFunction]

I'm sorry but, do you mind if you make a sample code?

Thank you for your kindness.

[taikichi-pc]

@Balandat
Is it, sample code, difficult?

[taikichi-pc]

@saitcakmak
Thank you for your reply!
I'll read and test it.

I have a additional question.
qPI's code(excerption of forward() part) is below, as you know,

    @concatenate_pending_points
    @t_batch_mode_transform()
    def forward(self, X: Tensor) -> Tensor:
        r"""Evaluate qProbabilityOfImprovement on the candidate set `X`.
        Args:
            X: A `batch_shape x q x d`-dim Tensor of t-batches with `q` `d`-dim design
                points each.
        Returns:
            A `batch_shape'`-dim Tensor of Probability of Improvement values at the
            given design points `X`, where `batch_shape'` is the broadcasted batch shape
            of model and input `X`.
        """
        posterior = self.model.posterior(
            X=X, posterior_transform=self.posterior_transform
        )
        samples = self.get_posterior_samples(posterior)
        obj = self.objective(samples, X=X)  # `sample_shape x batch_shape x q`-dim
        max_obj = obj.max(dim=-1)[0]  # `sample_shape x batch_shape`-dim
        impr = max_obj - self.best_f.to(max_obj)
        val = torch.sigmoid(impr / self.tau).mean(dim=0)
        return val

When we use qPI, we use reparameterizable acquisition functions(sigmoid func. should appear).
Your sample code above seems to be a little different from single-qPI.

Why is there a difference with single-qPI?
Using prod_val = torch.prod(val, 1) is wrong?

[saitcakmak]

I did not look at the implementation of qPI while writing the above example. Looks like there is an important piece that was missing in my implementation. The use of sigmoid is necessary to ensure that the acquisition function is differentiable, so that you can optimize it using optimize_acqf. So, something like this would work better:

class ProductOfProbabilityOfImprovement(MultiObjectiveMCAcquisitionFunction):
    def __init__(self, model: Model, best_f: Tensor, tau: float, sampler: MCSampler, objective: MCMultiOutputObjective) -> None:
        super().__init__(model=model, sampler=sampler, objective=objective)  # sampler & objective are optional
        self.tau = tau  # temperature parameter for sigmoid
        self.best_f = best_f

    def forward(self, X: Tensor) -> Tensor:
        posterior = self.model.posterior(X=X)  # Assuming a multi-output posterior
        samples = self.get_posterior_samples(posterior)  # Samples here are num_samples x batch_X x q x m (m being number of outputs)
        obj = self.objective(samples, X=X)  # If you use an objective, feel free to skip.
        # Assuming self.best_f to be an `m` dim tensor. Computing the sigmoid approximation of improvement
        improvement = obj - self.best_f.to(obj)
        max_improvement  = improvement.max(dim=-2).values  # Taking the maximum over the q dimension, assuming any one of them being an improvement is sufficient if q > 1 (e.g. in a parallel evaluation setup)
        poi = torch.sigmoid(max_improvement / self.tau).mean(dim=0)
        return poi.prod(dim=-1)  # return the product of PoI for each outcome

[taikichi-pc]

@saitcakmak
Thank you for your quick reply.
It's more like qPI.

But I guess below code in forward() is better(means more like single-qPI), right?
Is it the same as yours?

samples = self.get_posterior_samples(posterior)
obj = self.objective(samples, X=X)  # `sample_shape x batch_shape x q x models`-dim
max_obj = obj.max(dim=-2)[0]  # `sample_shape x batch_shape x models`-dim        
impr = max_obj - self.best_f.to(max_obj)
vals = torch.sigmoid(impr / self.eta).mean(dim=0) # mean of samples
prod_val = torch.prod(vals, 1) # batch_shape
return prod_val

[saitcakmak]

Yeah, looks like the two should produce the same results.

[taikichi-pc]

@saitcakmak
The sample code was properly executed.

I appreciate your help so much.

[An-xu]

Hi, may I ask how to use the EHVI acquisition function? I am really confused by the documents and the code I wrote has some errors.
After I fit the model, the code I wrote is:
ref_point = torch.tensor([10, 10])
EHVI = ExpectedHypervolumeImprovement(model=model, ref_point=ref_point, partitioning=NondominatedPartitioning(ref_point=ref_point) )

The error is: pareto_Y has not been initialized
ExpectedHypervolumeImprovement.init(self, model, ref_point, partitioning, objective)
125 if len(ref_point) != partitioning.num_outcomes:
126 raise ValueError(
127 "The length of the reference point must match the number of outcomes. "
128 f"Got ref_point with {len(ref_point)} elements, but expected "
129 f"{partitioning.num_outcomes}."
130 )
131 ref_point = torch.tensor(
132 ref_point,
--> 133 dtype=partitioning.pareto_Y.dtype,
134 device=partitioning.pareto_Y.device,
135 )
136 better_than_ref = (partitioning.pareto_Y > ref_point).all(dim=1)
137 if not better_than_ref.any() and partitioning.pareto_Y.shape[0] > 0:

Is there has any example that can show me how to use this acquisition function and how to fix this error?

Thanks a lot for helping me!!

[saitcakmak]

Hi @An-xu. You should pass in the Y observed so far into NondominatedPartitioning. The qEHVI parts of this tutorial are probably the closest reference we have for EHVI.